# Physico-chimie et modélisation du frittage des combustibles d'oxydes d'actinides

Jacques Léchelle

► To cite this version:

Jacques Léchelle. Physico-chimie et modélisation du frittage des combustibles d'oxydes d'actinides. Matériaux. CEA, CEA/DEN/CAD/DEC/SPUA/LMPC, 2014. cea-01090740

HAL Id: cea-01090740

https://cea.hal.science/cea-01090740

Submitted on 4 Dec 2014

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.

Public Domain
CEA-R-6379

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

ISSN 0429 - 3460

COMMISSARIAT À L'ÉNERGIE ATOMIQUE
ET AUX ÉNERGIES ALTERNATIVES

DIRECTION DE L'ÉNERGIE NUCLÉAIRE

PHYSICO-CHIMIE ET MODÉLISATION
DU FRITTAGE DES COMBUSTIBLES D'OXYDES D'ACTINIDES

par
Jacques LÉCHELLE

CEA CADARACHE
DIRECTION DE L'ÉNERGIE NUCLÉAIRE
DÉPARTEMENT D'ÉTUDES DES COMBUSTIBLES
SERVICE PLUTONIUM URANIUM ET
ACTINIDES MINEURS

2014

DIRECTION DES SYSTÈMES
D'INFORMATION

CEA / SACLAY 91191 GIF-SUR-YVETTE CEDEX FRANCE

RAPPORT

RAPPORT

CEA-R-6379
<figure><img src="image_2.png" /></figure>
# RAPPORT CEA-R-6379 – Jacques LÉCHELLE

## «Physico-chimie et modélisation du frittage des combustibles d'oxydes d'actinides»

**Résumé** - Ce document synthétise les travaux que j'ai menés ou auxquels j'ai contribué en tant que support numérique entre 1996 et 2012 au Service Plutonium Uranium et Actinides mineurs du CEA Cadarache. Ces travaux sont essentiellement orientés vers l'étude et la modélisation du frittage du combustible nucléaire qui servent de fil conducteur au mémoire. La modélisation du frittage a été abordée à une échelle sub-granulaire afin d'une part de prendre en compte les spécificités physico-chimiques des oxydes d'actinides qui constituent le matériau et d'autre part de pouvoir combiner les divers phénomènes de transport à l'origine du frittage. Le modèle se limite actuellement au cas d'une composition chimique constante à l'échelle sub-granulaire, sans échanges avec la phase gaz. Son extension à une composition variable est prévue. Des travaux sont en cours vers une simulation à l'échelle d'un volume représentatif de la pastille. Des moyens de validation à l'échelle de la pastille (Diffusion de Neutrons aux Petits angles), à l'échelle des granulés de poudre (microtomographie X Haute Résolution) ont été testés sur l'oxyde (U,Pu)O₂; d'autres sont en cours d'étude à l'Institut de Chimie Séparative de Marcoule à l'échelle du grain monocristallin sur CeO₂ (Microscopie Electronique à Balayage Environnemental). L'homogénéité requise pour la microstructure du combustible a amené à étudier expérimentalement le mode d'action de certains adjuvants de frittage tels que Cr₂O₃. L'ensemble des travaux permet d'améliorer la compréhension que l'on peut avoir du frittage et, de là, à pouvoir en assurer une meilleure maîtrise.

**Mots clés** : frittage, actinides, modèle, échelle sub-granulaire, élasticité 3D, transport de matière, diffusion aux joints de grains, diffusion en surface, diffusion en volume, méthode des éléments finis, méthode des volumes finis, évolution du maillage.

## 2014 – Commissariat à l'Énergie Atomique et aux Énergies Alternatives – France

# RAPPORT CEA-R-6379 – Jacques LÉCHELLE

## « Physical chemistry and modelling of the sintering of actinide oxides »

**Abstract** - This report gives a synthesis of the work I have carried out or to which I have numerically contributed from 1996 up to 2012 in the Department of Plutonium Uranium and minor Actinides in Cadarache CEA Center. Their main goal is the study and the modeling of the sintering process of nuclear fuels which is the unifying thread of this document. Both in order to take into account the physical and chemical features of the actinide bearing oxide material and in order to combine the different transport phenomena leading to sintering, a sub-granular scale model is under development. Extension to a varying chemical composition as well as exchanges with the gaseous phase are foreseen. A simulation on a larger scale (pellet scale) is ongoing in the framework of a PhD thesis. Validation means have been tested with (U,Pu) O₂ material on the scale of the pellet (Small Angle Neutron Diffusion), on the scale of powder granules (X-Ray High Resolution Micro-Tomography) and with CeO₂ at the 'Institut de Chimie Séparative' in Marcoule on a single crystal scale (Environmental Scanning Electron Microscope). The required microstructure homogeneity for nuclear fuels has led to a campaign of experimental studies about the role of Cr₂O₃ as a sintering aid. Whole of these studies improve our understanding of fuel sintering and hence leads to an improved mastering of this process.

**Keywords:** sintering, actinides, model, sub-granular scale, 3D-Elasticity, mass transport, grain boundary diffusion, surface diffusion, bulk diffusion, finite element method, finite-volume method, mesh evolution.

## 2014 – Commissariat à l'Énergie Atomique et aux Énergies Alternatives – France
<figure><img src="image_3.png" /></figure>
CEA Cadarache
Direction de L'Énergie Nucléaire
Département d'Études des Combustibles
Service Plutonium Uranium et Actinides Mineurs

PHYSICO-CHIMIE ET MODÉLISATION
DU FRITTAGE DES COMBUSTIBLES D'OXYDES D'ACTINIDES

par

Jacques LÉCHELLE
<figure><img src="image_4.png" /></figure>
# Physico-chimie et modélisation du frittage des combustibles d'oxydes d'actinides

## Habilitation à Diriger des Recherches

présentée et soutenue publiquement le 20 décembre 2013

Aix-Marseille Université

(spécialité : mécanique, génie mécanique)

par

Jacques LÉCHELLE

## Composition du jury

**Président :** Jean-Marc HEINTZ Professeur Université de Bordeaux I

**Rapporteurs :**
Sylvie BORDÈRE Chargée de Recherche - Dr - HDR ICMCB
Didier BOUVARD Professeur Grenoble INP
Anne LERICHE Professeur Université de Valenciennes et du Hainaut-Cambrésis

**Examinateurs :**
Jean-Marie GATT Maître de Conférences - HDR INSTN
Pierre HALDENWANG Professeur Aix-Marseille Université
Jacques ROGEZ Directeur de Recherche - CNRS - Aix-Marseille Université
Mis en page avec la classe thloria.
Remerciements

Je remercie le Département d'Étude des combustibles, le Service Plutonium Uranium et Actinides Mineurs et le Laboratoire d'Étude des Microstructures et des Propriétés du Combustible et leurs responsables de m'avoir encouragé et aidé à préparer cette habilitation à diriger des recherches.

Ces démarches n'auraient pas pu être envisagées si je n'avais pas été aidé et soutenu au quotidien par les collègues du laboratoire, du service et aussi du département. A chaque occasion, j'ai pu bénéficier d'idées et d'aide de la part de ces collègues pour réaliser des expériences pas gagnées d'avance, et ce, dans les meilleures conditions.

Je remercie aussi Robert Boyer du Centre de Mathématiques et d'Informatique de Marseille de
son aide ininterrompue depuis 1997 dans la modélisation du frittage du point de vue mécanique,
mathématique et numérique. Je remercie Kacem Saikouk de son aide sur les aspects numériques
et ceux de la mise en œuvre informatique du modèle.

Je remercie Marie-Hélène Mathon du Laboratoire Léon Brillouin pour son aide quasiment inconditionnelle pour les essais de diffusion de neutrons aux petits angles et pour le temps qu'elle a pu m'accorder.

Je remercie également Christophe Valot, Méryl Brothier et Philippe Sornay pour leurs conseils et leur aide précieuse pour la préparation de la présentation orale.

Je remercie mon "tuteur" Pierre Haldenwang qui m'a soutenu dans la démarche auprès de l'Université de Marseille, avec l'aide de Géraldine Callier qui en maîtrise les rouages administratifs.

Je remercie Sylvie Bordère, Anne Leriche, et Didier Bouvard d'avoir accepté d'être rapporteurs de ce manuscrit. Je remercie Jean-Marc Heintz, Jacques Rogez et Jean-Marie Gatt, qui ont accepté de participer au jury.

Je remercie Sylvie Thabet pour l'organisation de ma soutenance.

Je remercie ma famille et mes amis de m'avoir soutenu et supporté, dans tous les sens du terme, et d'avoir été patients tout au long de ma préparation de HDR.
À mes Parents,
à Florence,
à François.
# Table des matières

Introduction générale

1

## Partie I Description classique du frittage en phase solide

Introduction - le frittage en phase solide

5

### Chapitre 1

Phénomènes de transport lors du frittage en phase solide

<table><tr><td>1.1</td><td>Diffusion en surface</td><td>9</td></tr><tr><td>1.1.1</td><td>Définition</td><td>9</td></tr><tr><td>1.1.2</td><td>Ordre de grandeur</td><td>10</td></tr><tr><td>1.2</td><td>Diffusion aux joints de grains</td><td>11</td></tr><tr><td>1.2.1</td><td>Définition - le joint de grain</td><td>11</td></tr><tr><td>1.2.2</td><td>Ordre de grandeur</td><td>12</td></tr><tr><td>1.3</td><td>Diffusion en volume</td><td>13</td></tr><tr><td>1.4</td><td>Évaporation/condensation</td><td>14</td></tr><tr><td>1.5</td><td>Conclusion</td><td>15</td></tr></table>

### Chapitre 2

Conséquences microstructurales des phénomènes de transport

<table><tr><td>2.1</td><td>Les différents stades du frittage</td><td>17</td></tr><tr><td>2.2</td><td>Lois du premier stade du frittage</td><td>17</td></tr><tr><td>2.2.1</td><td>Loi d'évolution des cous en l'absence de retrait</td><td>18</td></tr><tr><td>2.2.2</td><td>Loi d'évolution des cous avec retrait</td><td>19</td></tr><tr><td>2.3</td><td>Lois du second stade du frittage</td><td>19</td></tr><tr><td>2.3.1</td><td>Diffusion en volume</td><td>19</td></tr><tr><td>2.3.2</td><td>Diffusion aux joints de grains</td><td>20</td></tr><tr><td>2.3.3</td><td>Lois de densification</td><td>20</td></tr></table>
<table><tr><td>2.4</td><td>Lois du troisième stade du frittage</td><td>20</td></tr><tr><td>2.5</td><td>Limites des modèles classiques</td><td>21</td></tr><tr><td>2.6</td><td>Conclusion</td><td>21</td></tr></table>

## Chapitre 3

### Maîtrise de la microstructure - Cartes de microstructure 23

<table><tr><td>3.1</td><td>Définition</td><td>23</td></tr><tr><td>3.2</td><td>Exemple d'utilisation</td><td>24</td></tr><tr><td>3.3</td><td>Conclusion</td><td>25</td></tr></table>

## Chapitre 4

### Le frittage à l'échelle des milieux continus

<table><tr><td>4.1</td><td>Notion de contrainte de frittage</td><td>27</td></tr><tr><td>4.2</td><td>Cas d'un matériau biphasé - contrainte différentielle</td><td>27</td></tr><tr><td>4.3</td><td>Conclusion</td><td>30</td></tr></table>

## Chapitre 5

### Conlusion

# Partie II Les spécificités des oxydes d'actinides

## Chapitre 1

### Le caractère radioactif

<table><tr><td>1.1</td><td>La décroissance des isotopes des éléments majoritaires U et Pu</td><td>37</td></tr><tr><td>1.1.1</td><td>Vecteur isotopique typique du combustible vierge</td><td>37</td></tr><tr><td>1.1.2</td><td>Isotopes de U dans les combustibles</td><td>38</td></tr><tr><td>1.1.3</td><td>Isotopes de Pu dans les combustibles</td><td>38</td></tr><tr><td>1.2</td><td>La décroissance des isotopes des actinides mineurs Am, Np et Cm</td><td>39</td></tr><tr><td>1.3</td><td>La thermodynamique pour les isotopes radioactifs</td><td>39</td></tr><tr><td>1.3.1</td><td>Les fonctions d'état</td><td>39</td></tr><tr><td>1.3.2</td><td>1<sup>er</sup> principe de la thermodynamique</td><td>40</td></tr><tr><td>1.3.3</td><td>2<sup>nd</sup> principe de la thermodynamique</td><td>40</td></tr><tr><td>1.3.4</td><td>3<sup>eme</sup> principe de la thermodynamique</td><td>41</td></tr><tr><td>1.3.5</td><td>Conclusion</td><td>42</td></tr></table>

## Chapitre 2

### Propriétés chimiques des oxydes d'actinide
<table><tr><td>2.1</td><td>Variations de comportement des électrons de valence dans la série des actinides</td><td>43</td></tr><tr><td>2.2</td><td>Écarts à la stœchiométrie</td><td>44</td></tr><tr><td>2.3</td><td>Voisinage de la stoechiométrie - défauts ponctuels</td><td>45</td></tr><tr><td>2.3.1</td><td>Équations de diffusion chimique pour les oxydes</td><td>46</td></tr><tr><td>2.3.2</td><td>Simulation de l'interdiffusion lors du frittage</td><td>47</td></tr></table>

## Chapitre 3

### Spécificités liées à la cristallographie du matériau

<table><tr><td>3.1</td><td>Structure cristalline</td><td>49</td></tr><tr><td>3.2</td><td>Symétrie des propriétés élastiques</td><td>49</td></tr><tr><td>3.3</td><td>Symétrie des coefficients de diffusion</td><td>50</td></tr><tr><td>3.4</td><td>Symétrie de la tension de surface</td><td>50</td></tr><tr><td>3.5</td><td>Le facettage des grains</td><td>50</td></tr><tr><td>3.6</td><td>Conséquences sur les équations de conservation - spécificité des cristaux ion-niques et de la symétrie cubique</td><td>52</td></tr></table>

## Chapitre 4

### Conclusion

# Partie III Modèle développé à l'échelle sub-granulaire

## Chapitre 1

### Cas d'une composition chimique constante

<table><tr><td>1.1</td><td>Choix du modèle</td><td>59</td></tr><tr><td>1.2</td><td>Hypothèses du modèle</td><td>61</td></tr><tr><td>1.3</td><td>Effets de la courbure de surface</td><td>62</td></tr><tr><td>1.3.1</td><td>Calcul numérique de la courbure</td><td>62</td></tr><tr><td>1.3.2</td><td>Effet mécanique de la tension de surface</td><td>62</td></tr><tr><td>1.3.3</td><td>Conditions aux limites de type Neumann</td><td>63</td></tr><tr><td>1.3.4</td><td>Conditions algébriques</td><td>63</td></tr><tr><td>1.3.5</td><td>Résolution numérique du problème</td><td>65</td></tr><tr><td>1.3.6</td><td>Comportement élastique</td><td>65</td></tr><tr><td>1.4</td><td>Potentiel chimique généralisé</td><td>69</td></tr><tr><td>1.5</td><td>Évolution des interfaces entre sous-domaines</td><td>70</td></tr><tr><td>1.5.1</td><td>Modélisation des flux</td><td>70</td></tr><tr><td>1.5.2</td><td>Évolution des interfaces grain/gaz</td><td>70</td></tr><tr><td>1.5.3</td><td>Évolution des joints de grains</td><td>71</td></tr></table>
<table><tr><td>1.5.4</td><td>Calcul de la vitesse d'interface</td><td>72</td></tr><tr><td>1.5.5</td><td>Problème de la conservation de la matière</td><td>73</td></tr><tr><td>1.5.6</td><td>Problème de l'évolution de l'angle dièdre</td><td>74</td></tr><tr><td>1.5.7</td><td>Qualité du maillage de surface</td><td>76</td></tr><tr><td>1.6</td><td>Choix du pas de temps</td><td>76</td></tr><tr><td>1.7</td><td>Conclusion</td><td>77</td></tr></table>

## Chapitre 2

### Résultats du modèle - composition constante

<table><tr><td>2.1</td><td>Rayons de courbure principaux</td><td>81</td></tr><tr><td>2.2</td><td>Champ de déplacement</td><td>82</td></tr><tr><td>2.3</td><td>Évolution de la taille du cou en fonction du temps</td><td>82</td></tr><tr><td>2.3.1</td><td>Diffusion en surface et aux joints de grains</td><td>83</td></tr><tr><td>2.3.2</td><td>Diffusion en volume</td><td>84</td></tr><tr><td>2.4</td><td>Conclusion</td><td>86</td></tr></table>

## Chapitre 3

### Composition chimique variable - interdiffusion

<table><tr><td>3.1</td><td>Spécificités du combustible (U, Pu)O<sub>2</sub></td><td>87</td></tr><tr><td>3.2</td><td>Interdiffusion et modèle à l'échelle sub-granulaire</td><td>87</td></tr><tr><td>3.3</td><td>Études visant à rendre le (U, Pu)O<sub>2</sub> plus homogène</td><td>88</td></tr></table>

## Chapitre 4

### Conclusion

# Partie IV Acquisition des paramètres d'entrée du modèle à l'échelle sub-granulaire

## Chapitre 1

### Tension de surface

## Chapitre 2

### Constantes élastiques isothermes

## Chapitre 3

### Coefficients de diffusion

<table><tr><td>3.1</td><td>Diffusion en volume</td><td>101</td></tr></table>
# Introduction générale

Le combustible nucléaire est fabriqué, en France, par métallurgie des poudres, avec, parmi les étapes clés du procédé, une étape dite de frittage. Elle consiste en la consolidation d'une pastille de poudre sous l'effet de la chaleur en dessous de la température de fusion de ses constituants. La présente étude vise à apporter des outils de compréhension du rôle des phénomènes physico-chimiques qui ont lieu durant le frittage du combustible à base de $UO_2$ et de $PuO_2$ en vue de pouvoir agir sur cette étape de fabrication que ce soit par un choix approprié d'additifs, de mode d'introduction, d'atmosphère, de cycle thermique voire de morphologie ou de distribution de taille d'agglomérats de poudre. Elle correspond au travail effectué de 1996 à 2012 au CEA de Cadarache à la Direction de l'Énergie Nucléaire au Département d'Étude des Combustibles, Service Plutonium Uranium et Actinides Mineurs, Laboratoire des Microstructures et Propriétés des Combustibles.

Le document est découpé en cinq parties et plusieurs annexes :

- Une description classique du frittage est donnée à l'échelle sub-granulaire, avec les phénomènes qui interviennent à l'échelle atomique (diffusion en surface, aux joints de grains, en volume, évaporation/condensation) et leurs conséquences à l'échelle micrométrique : grossissement des grains, élimination des macro et méso porosités, ainsi qu'à l'échelle de la pastille : retrait, élimination de la porosité ouverte, élimination d'une partie de la porosité fermée.

- Une analyse est faite des spécificités des matériaux mis en jeu. Les oxydes d'actinides présentent un caractère radioactif avec des conséquences sur la mise en œuvre de la thermodynamique. Leurs électrons de valence 6d/5f présentent une variation du caractère ionique le long de la série des actinides. Ces oxydes sont ainsi sensibles à la pression partielle d'oxygène et sont non-stœchiométriques. La structure cristalline dans laquelle ils s'organisent a des conséquences sur la symétrie des propriétés (élastiques, diffusionnelles, de tension de surface) qui se reflète dans le facetage des grains d'oxyde lors de leur croissance, ainsi que sur la prise en compte mathematique de la conservation de la matière lors des phénomènes de diffusion (conservation des sites par unité de masse ou de volume de matériau suivant que le matériau soit ou pas considéré comme compressible).

- Les deux aspects précédents ayant été pris en considération, un modèle numérique, à l'échelle sub-granulaire, de l'évolution de quelques grains au contact a été développé. Les variations de rayon de courbure le long de la surface des grains induit un champ de contraintes dans le grain qui va engendrer des variations locales de potentiel chimique, puis des flux de matière. Ces flux vont déformer les surfaces des grains (surfaces libres et joints de grains) et modifier leur volume. Ce modèle est illustré dans le cas de deux grains et a été testé pour ce cas de figure avec des grains de même taille ou de taille différente et une composition chimique constante.

- Les moyens d'acquisition des données utilisées par le modèle, essentiellement les coefficients de diffusion, ont été étudiés et mis en œuvre par d'autres équipes du laboratoire pour un nombre restreint de conditions thermodynamiques. Une des spécificités de ces matériaux étant leurs faibles coefficients de diffusion, une méthode numérique spécifique a été développée et appliquée
pour le dépouillement de profils d'interdiffusion. Afin de pouvoir ajuster leur valeur par le biais d'additifs ou en contrôlant l'atmosphère de frittage, un dispositif de mesure de la conductivité électrique nucléarisable a été développé, et des simulations numériques ont permis d'adapter l'interprétation de la technique de mesure à la géométrie de la pastille et à celle des électrodes. Cet outil permettra de déterminer la nature des porteurs de charge, ainsi que leur concentration et, éventuellement, à terme, de relier les coefficients d'interdiffusion aux concentrations en défauts ponctuels.

- Enfin, la dernière partie de ce document vise à définir les moyens nécessaires à la validation du modèle à l'échelle sub-granulaire. C'est, notamment, la tomographie haute résolution, qui, grâce à l'énergie des X des sources synchrotrons, peuvent traverser jusqu'à 200µm de UO₂ prélevés sur un matériau fritté à différents stades du frittage pour cartographier en 3D la porosité en mode transmission ou le réseau de joints de grains en contraste de phase. C'est aussi la microscopie électronique à balayage environnementale en 2D qui permet de suivre l'évolution de la surface extérieure de deux grains au contact (observation dans le plan de focalisation du MEB). Pour des échelles macroscopiques de modèles de frittage, la diffusion de neutrons aux petits angles permet de traverser une pastille centimétrique de ces matériaux en renseignant de manière non destructive sur l'aire totale de l'interface grain/gaz ainsi que sur les rayons de courbure des pores lorsque ces derniers sont nanométriques.

Les annexes présentent la liste des publications et des brevets, les thèses et stages encadrés, les activités d'enseignement et d'animation scientifique, trois articles représentatifs de l'étude, le rapport de soutenance de thèse et les diagrammes de décroissance radioactive des isotopes de U, Pu et Am.
Première partie

Description classique du frittage en
phase solide
<figure><img src="image_5.png" /></figure>
Introduction - le frittage en phase solide
<figure><img src="image_6.png" /></figure>
Le frittage en phase solide est la densification d'un matériau initialement pulvérulent sous l'effet de la chaleur à une température inférieure à celle de son point de fusion. On se limite ici, plus particulièrement, au cas où le matériau est cristallisé, ce qui est le cas des oxydes d'uranium et de plutonium.

## La vision des chimistes

Pour les matériaux considérés, dans les conditions de température et de pression d'usage, l'ensemble des ions les constituant, à leur position d'équilibre dans le cristal, sont plus stables que les atomes correspondants à l'infini les uns des autres (l'énergie libre de l'ensemble est plus faible). Les ions de surface ne profitent pas de ce gain de stabilité car il leur manque des ions voisins. La différence d'énergie libre entre le solide de dimensions finies isolé et le même solide considéré à l'intérieur d'un milieu infini constitué de ce même matériau, rapportée à la surface du solide est l'énergie de surface. Les ions en surface, selon la courbure de la surface (en marches d'escalier à l'échelle des ions vont être dans des environnements différents et avoir un potentiel chimique variable le long de cette surface. Ils vont alors diffuser pour uniformiser le potentiel chimique.

## La vision des mécaniciens

Tout se passe comme si le solide, de petites dimensions, considéré comme un milieu continu était recouvert d'une membrane. Cette dernière met en compression le solide là où la surface est convexe et en tension les zones concaves. Si l'on coupe sur une petite longueur cette membrane, les bords de la coupure vont s'écarter du fait de la force exercée par le reste de la membrane sur la coupure. Cette force sera tangente à la surface et perpendiculaire à la coupure. La tension de surface est alors la norme de cette force rapportée à l'unité de longueur de coupure. L'énergie libre emmagasinée dans le solide du fait de ces contraintes locales en surface va induire des variations de potentiel chimique généralisé que la diffusion va tendre à faire disparaître.

## La vision pour la modélisation : thermodynamique linéaire des processus irréversibles

Autant la vision de la chimie du matériau permet la compréhension des phénomènes, autant elle est inadaptée à la modélisation de l'évolution de la taille de grain et de l'élimination de la porosité du fait de son échelle trop petite. La vision mécanicienne du matériau permet de traiter un grain comme un milieu continu et est plus adaptée. Elle induit, cependant, des limitations ou des difficultés. La première limitation est que le tout début d'un contact entre deux grains qui ne l'étaient pas préalablement au frittage ne peut pas être traité simplement : le contact ponctuel sera difficile à gérer numériquement. Un autre aspect est que la notion de joint de grains est liée à l'organisation cristalline qui n'existe pas explicitement en mécanique des milieux continus. On donne alors aux grains, à l'instant initial, chacun une orientation qu'ils gardent tant qu'ils existent (i.e. tant qu'ils n'ont pas été consommés par la diffusion de leur matière constitutive vers les grains adjacents). La matière qui se dépose sur un grain du fait des flux de diffusion va se déposer sur le grain en s'alignant sur les positions cristallines des atomes voisins déjà présents (épitaxie). Le système étudié étant hors équilibre, les processus diffusionnels sont pris en compte et décrits par la thermodynamique linéaire des processus irréversibles.
1

Phénomènes de transport lors du
frittage en phase solide

En fonction de la température (et éventuellement pour les oxydes, dans une moindre me-
sure, de la pression partielle d'oxygène) la diffusion peut ne pas être le seul moyen de transport
permettant d'uniformiser le potentiel chimique, en effet, le transport de matière peut aussi se
produire via les phénomènes d'évaporation/condensation.

Cette diffusion est thermiquement activée et très lente dans les dioxydes d'actinides. Sur des
durées aussi limitées que celles mises en jeu dans le frittage (seulement quelques heures pour un
frittage classique de (U, Pu)O₂), elle ne peut intervenir de façon notable qu'à haute température
(au delà d'environ 1200°C). En dessous de cette température l'évaporation/condensation (que
l'on devrait plutôt appeler sublimation/condensation) permet, via la phase gazeuse, de trans-
porter les ions sous forme moléculaire des zones en compression pour les déposer dans les zones
en traction, en constituant de nouvelles mailles cristallines en lissant la surface. Ce phénomène
nécessite que la pression partielle des espèces en phase gaz soit non négligeable et intervient,
typiquement, entre 800°C et 1200°C. En dessous de 800°C les pressions partielles d'espèces
contenant les éléments constitutifs du solide est négligeable, au dessus de 1200°C l'aire de l'in-
terface solide/gaz devient plus petite ce qui diminue les échanges entre phases solide et gazeuse,
et la fermeture de la porosité limite la distance sur laquelle le transport en phase gaz peut avoir
lieu. La phase du cycle thermique de frittage dans laquelle il intervient est très limitée dans le
temps pour le frittage standard du combustible au laboratoire, comme le montre la figure 1.1
(2h40' pour une durée totale, au dessus de la température de 800°C, de 16h), mais pourrait ne
pas l'être pour d'autres cycles thermiques.

On traite, par la suite, des seuls phénomènes de diffusion à l'état solide.

1.1 Diffusion en surface

1.1.1 Définition

En microscopie tunnel à balayage, on peut observer des terrasses d'autant plus larges que la
surface se rapproche d'un plan cristallographique de haute symétrie, séparées par des marches de
hauteur atomique. Elles ne suivent elles-mêmes pas une direction cristallographique et présentent
des crans. On observe, sur ces terrasses, des lacunes qui permettent de voir la terrasse sous-jacente
appelées ad-lacunes ainsi que des atomes isolés positionnés sur ces terrasses, appelés ad-atomes.
<figure><img src="image_7.png" /><figcaption>FIGURE 1.1 – cycle thermique typique d’un frittage de (U, Pu)O₂ au laboratoire</figcaption></figure>

On peut définir le coefficient de diffusion $D_S$ d'un ad-atome à partir du carré de la distance moyenne qu'il parcourt sur une durée $t$, la moyenne étant obtenue en choisissant des instants initiaux différents $t'$. Soit $\vec{r}(t)$ sa position à l'instant $t$, $d$ la dimension de l'espace dans laquelle la diffusion a lieu (ici $d = 2$) :

$$ \left\langle \left[ \vec{r}(t' + t) - \vec{r}(t') \right]^2 \right\rangle = 2dD_S t \qquad (1.1) $$

$D_S$ peut être mesuré en déposant une mono-couche de traceur sur une moitié de surface et en mesurant son étalement après diffusion. De même, on peut définir un coefficient de diffusion en surface pour une ad-lacune. Celui-ci est noté $\Lambda_S$. On montre, [JV1995], que le coefficient de diffusion de la loi de Fick (flux causés par des gradients de potentiel chimique et non pas des gradients de concentration) en surface $\tilde{D}$ est relié aux coefficients précédents par :

$$ \tilde{D}_S = \Lambda_S \sigma + D_S \rho \qquad (1.2) $$

où $\sigma$ et $\rho$ désignent respectivement les concentrations en ad-lacunes et ad-atomes. On voit que ce coefficient va dépendre, a priori, de l'orientation de la surface.

### 1.1.2 Ordre de grandeur

Henney [JH1968] a étudié la diffusion en surface de U sur $UO_2$ à haute température de manière indirecte à partir de la vitesse de déplacement d'un joint de grains mesurée par inter-férométrie optique ainsi que d'un modèle de type Mullins [WM1957] de déplacement du joint de grain. L'étude de Zhou et Olander de 1983 [SYZ1983] est la plus récente trouvée dans la littérature pour déterminer un coefficient de diffusion en surface pour $UO_2$, elle a été menée en
utilisant l'isotope $^{234}$U. La figure 1.2 montre que le coefficient de diffusion en surface de U vaut
$0,1 \text{ cm}^2\text{s}^{-1}$ à $1700^{\circ}\text{C}$ tandis qu'il est proche de $1 \text{ cm}^2\text{s}^{-1}$ à $2100^{\circ}\text{C}$.

<figure><img src="image_8.png" /><figcaption>FIGURE 1.2 – Coefficient de diffusion en surface de U sur UO₂</figcaption></figure>

Les valeurs du coefficient de diffusion de $O^{2-}$ qui diffuse bien plus rapidement que $U^{n+}$ et qui n'est pas limitant dans la cinétique du frittage ne sont pas discutées.

1.2 Diffusion aux joints de grains

1.2.1 Définition - le joint de grain

Le joint de grain est une zone d'accommodation des deux réseaux cristallographiques des grains qu'il sépare. La description de la structure cristallographique des joints de grains dans le cas des métaux, leur nomenclature et l'étude de leurs propriétés date du début de la microscopie électronique [MB1975]. Les simulations récentes ont aussi fait progresser la connaissance des joints de grains des céramiques (la figure 1.3 montre l'agencement des ions de part et d'autre d'un joint de grain de UO₂ pour lequel des coefficients de diffusion ont été calculés).

On donne, à titre d'exemple, en figure 1.4, l'énergie de joint de grain calculée pour l'aluminium lorsqu'on fait varier l'angle de rotation autour d'un axe [011]. Sur cet exemple, parmi les joints de forte désorientation, les macles (111) et (113) présentent des minima d'énergie très marqués, même à haute température, alors que la macle (112), de faible indice elle aussi, se comporte comme un joint quelconque. On voit que l'énergie et les propriétés des joints de grains varient
<figure><img src="image_9.png" /><figcaption>FIGURE 1.3 – Projection selon la direction de tilt [001] d'un joint de grain symétrique (310)/[001]Σ₅ pour UO₂. L'interface est en pointillés. Les grosses sphères représentent les ions uranium, les petites les ions oxydes [EVA2009].</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>FIGURE 1.4 – Comparaison des énergies calculées et mesurées, dans le cas des joints de flexion d'axe < 110 >, dans l'aluminium d'après [MB1975]</figcaption></figure>

fortement avec la désorientation entre les deux sous-réseaux.

## 1.2.2 Ordre de grandeur

Les ions qui diffusent aux joints de grains sont freinés par les réseaux cristallins des deux grains de part et d'autre du joint. Cette diffusion, plus lente que la diffusion en surface est souvent étudiée à haute température. Yajima [SY1966] a étudié la diffusion en volume et aux joints de grains des ions uranium dans UO₂ polycristallin dans l'intervalle de température 1900–2150°C en utilisant comme traceur l'isotope <sup>237</sup>U et la méthode de sectionnement. La courbe de pénétration était constituée de deux parties. Dans la première partie, proche de la surface, la diffusion en volume et la diffusion aux joints de grains contribuent toutes deux à la pénétration des ions uranium diffusants. Dans la seconde partie, qui correspond à une profondeur de pénétration plus
grande, la diffusion intergranulaire devient prépondérante et la diffusion en volume est négligeable comparée à la première. Les données expérimentales sont analysées à l'aide d'un modèle de source instantanée. D'autre part, en appelant a la demi-largeur du joint de grain, le produit $D' \cdot 2a$, où $D'$ est le coefficient de diffusion aux joints de grains, suit une loi de variation avec la température $D' \cdot 2a = 5,19 \cdot 10^{-9} e^{-47200/RT} cm^3 sec^{-1}$. Pour un joint de grain d'une épaisseur de 10 Å on obtient $D' = 0,474 cm^2 s^{-1}$ à 2100°C.

Là encore, la simulation à partir de la dynamique moléculaire permet de calculer avec une bonne précision les valeurs des coefficients de diffusion de U aussi bien que de O aux joints de grains de $UO_2$ en fonction de la température (de 300 K à 2100 K dans l'étude citée). Vincent-Aublant et al. [EVA2009] ont ainsi étudié deux joints de grains particuliers : $\Sigma_5$ et $\Sigma_{41}$ Ils ont montré que l'autodiffusion en volume de U et de O est augmentée de plusieurs ordres de grandeurs par rapport au volume du grain et ce sur quelques nanomètres de distance du joint. Les mécanismes de saut ont été identifiés pour le joint $\Sigma_{41}$ dans toutes les directions et il a été montré que pour un joint $\Sigma_5$ le mécanisme de saut n'a lieu que le long du joint de grain. Les valeurs des coefficients de diffusion sont données sur la figure 1.5.

<figure><img src="image_11.png" /><figcaption>FIGURE 1.5 – Coefficients d'autodiffusion de l'uranium et de l'oxygène en fonction de $\frac{1}{k_B T}$ calculés pour des atomes de coeur de joint de grains $\Sigma_5$ et $\Sigma_{41}$</figcaption></figure>

Mendez [SM1995-2],[SM1995], puis plus récemment Noyau [SN2012] ont mesuré sur $(U, Pu) O_{2-x}$ des coefficients de diffusion aux joints de grains entre 1500°C et 1700°C pour un potentiel d'oxygène fixé. Mendez avait utilisé comme technique expérimentale la seule microsonde électronique, tandis que Noyau avait utilisé la spectrométrie de masse à ionisation secondaire (SIMS) avec un dépôt préalable d'un traceur de Pu ainsi que la microsonde électronique sur des profils obliques.

## 1.3 Diffusion en volume

L'expérience de Yajima [SY1966] donne, pour le coefficient de diffusion en volume $D_V$ de l'uranium, $D_V = 5,82 \cdot 10^{-5} e^{-72700/RT} cm^2 sec^{-1}$, soit à 1750°C $7,72 10^{-7} cm^2 s^{-1} = 7,72 10^{-11} m^2 s^{-1}$. Mendez [SM1995-2],[SM1995] a mesuré des coefficients d'interdiffusion globaux (volume et joints de grains) à partir de l'observation à la microsonde électronique de l'élément Pu dans des couples d'interdiffusion $UO_2/PuO_2$, pour en déduire, par une méthode de dépouillement adaptée, le coefficient de diffusion en volume. Elle obtient :

$$ \tilde{D}_V = 6.2 10^{-14} e^{-\frac{130000}{R}(\frac{1}{T} - \frac{1}{2023})} cm^2 \cdot s^{-1} \quad (1.3) $$
Noyau [SN2012] a, de plus, mesuré des coefficients d'autodiffusion en volume à l'aide du SIMS en déconvolutant dans ses profils la part de diffusion en volume de celle aux joints de grains. Ce coefficient, très faible, est très sensible à la concentration en défauts ponctuels intrinsèques du matériau et donc à la pression partielle d'oxygène.

## 1.4 Évaporation/condensation

La température de frittage du combustible est choisie à partir de l'oxyde majoritaire; celui de $UO_2$ dont la température de frittage est prise proche des deux tiers de son point de fusion soit $T_{frittage} = \frac{2}{3}T_F = \frac{2}{3}3120 = 2080K = 1807^\circ C$. Dans la pratique la température choisie est de $1750^\circ C$. Mais à cette température les autres oxydes du combustibles ont une pression de vapeur saturante non négligeable car leur point de fusion est inférieur à celui de $UO_2$ comme on le voit dans la tableau 1.1.

<table><caption>TABLE 1.1 – Température de fusion des oxydes d'actinides composant le combustible</caption><thead><tr><th>Oxyde</th><th>ΔH<sup>0</sup><sub>formation</sub>kJ·mol<sup>-1</sup></th><th>T<sub>Fusion</sub>(K)</th></tr></thead><tbody><tr><td>UO<sub>2</sub></td><td>-1085</td><td>3120</td></tr><tr><td>NpO<sub>2</sub></td><td>-1074</td><td>2820</td></tr><tr><td>PuO<sub>2</sub></td><td>-1055,8</td><td>2691</td></tr><tr><td>AmO<sub>2</sub></td><td>-932,3</td><td>2450</td></tr></tbody></table>

La loi recommandée pour les études de réacteurs pour la pression de vapeur de $UO_2(g)$ en équilibre avec $UO_2$ solide repose sur les mesures d'Ackermann, Rauh, et Rand dans la gamme 1800 à 2600 K [IAEA2006]. L'équation de pression d'équilibre de $UO_2$ qu'ils donnent est la suivante :

$$log_{10}(P) = 66,53672 + 4,382 10^{-3}T - 4,411 10^{-7}T^2 - \frac{37090}{T} - 19,070 log_{10}(T) \quad (1.4)$$

où $T$ est la température en Kelvin. Mais, du fait de l'existence d'autres espèces contenant de l'uranium en phase gaz, la pression totale des espèces uranifères est plus élevée que celle donnée par 1.4. Plus généralement, pour $(U, Pu)O_2$ la sublimation est non-congruente et la phase vapeur a un $\frac{Pu}{U+Pu}$ différent de celui de la phase solide. Un modèle classiquement utilisé pour décrire cette vaporisation est celui de Rand et Markin [HR1968]. Les espèces prises en compte sont : $O, O_2, U, UO, UO_2, UO_3, Pu, PuO, PuO_2$ :

$$ \left\{ \begin{array}{ll} \frac{1}{2} O_2(g) & \rightleftharpoons O(g) \\ U(g) + \frac{1}{2} O_2(g) & \rightleftharpoons UO(g) \\ UO(g) + \frac{1}{2} O_2(g) & \rightleftharpoons UO_2(g) \\ UO_2(g) + \frac{1}{2} O_2(g) & \rightleftharpoons UO_3(g) \\ Pu(g) + \frac{1}{2} O_2(g) & \rightleftharpoons PuO(g) \\ PuO(g) + \frac{1}{2} O_2(g) & \rightleftharpoons PuO_2(g) \end{array} \right. \quad (1.5) $$

La figure 1.6 montre les variations des pressions partielles de ces espèces, à 1970 K, en équilibre avec la phase solide $PuO_{2-x}$, sur la base d'ajustements plus récents de données thermodynamiques (C. Guéneau [CG2008]).
<figure><img src="image_12.png" /><figcaption>FIGURE 1.6 – Variations des pressions totale et partielles en fonction de O/Pu à 1970 K pour un système fermé - Les données expérimentales de Battles et Ackermann sont reportées dans le domaine bi-phasé [CG2008]</figcaption></figure>

Le modèle de frittage développé à l'échelle sub-granulaire incluait initialement l'inter-diffusion chimique en phase solide ainsi que les échanges avec la phase gazeuse [JL2001]. En 2003, dans le cadre de la thèse de Mounia Ajdour [MA2006], il a été décidé de progresser, dans un premier temps, plus spécifiquement sur le frittage lui-même avec une composition chimique constante. Les échanges avec la phase gaz ont ainsi, eux aussi, été provisoirement supprimés de ce modèle simplifié. Les développements liés à une composition chimique variable et à des échanges avec la phase gazeuse pourront être repris une fois le modèle simplifié validé afin de prendre en compte les phénomènes d'évaporation/condensation.

## 1.5 Conclusion

Les phénomènes élémentaires intervenant lors du frittage des oxydes d'actinides constituant les combustibles nucléaires ont été étudiés expérimentalement et quantifiés dans la littérature ou dans notre laboratoire. Cela permet de les mettre en œuvre, au moins dans un certain domaine de température et de pression partielle d'oxygène, dans le modèle de frittage à l'échelle sub-granulaire que nous avons développé.
# Conséquences microstructurales des phénomènes de transport

## 2.1 Les différents stades du frittage

On distingue classiquement [DBA2005] trois stades pour décrire le frittage sur la base de l'évolution de la microstructure du matériau :

- le premier stade pendant lequel les "particules" initiales (cristallites monocristallines de 0,1 à 0,3µm dans le cas de $UO_2$ et $PuO_2$ qui laisseront place au terme du frittage à des grains au sens métallurgique) en contact se soudent les unes aux autres. On observe une formation de "ponts" ou de "cous" entre particules. Un des rayons de courbures principaux des cous étant négatif la matière va diffuser vers ces cous avec l'apparition d'une phase poreuse interconnectée. Le joint de grain s'étend au niveau du cou, tandis qu'en surface le rayon des cous augmente en valeur absolue.

- le second stade du frittage montre des grains polyédriques et constitue l'étape dans laquelle la porosité va diminuer. Les grains grossissent pour atteindre 5 à 10 fois leur taille initiale. Ce stade continue jusqu'à atteindre une densité théorique de 95% alors que typiquement les matériaux pour le combustible nucléaire lorsqu'ils sont "crus", i.e. non-frittés, présentent une densité initiale (à l'issue du pressage) de 60% de la densité théorique.

- le troisième stade permet l'élimination de la majeure partie des pores isolés intra et intergrahnulaires ainsi qu'une augmentation de la taille des grains.

Dans les modèles classiques de frittage, les cristallites (on n'utilise pas le terme "grains" car à ce stade la taille de ces particules monocristallines est petite) sont représentées par des objets géométriques simples monodisperses de composition chimique constante et dont la forme varie selon le stade considéré de frittage. Les entités qui frittent sont bien supposées monocristallines puisque ces modèles distinguent pour les coefficients de diffusion mis en jeu ceux aux joints de grains, en volume et en surface.

## 2.2 Lois du premier stade du frittage

La figure 2.1 montre la géométrie sphérique utilisée pour représenter les grains (qui sont plutôt des cristallites à ce stade) lors du premier stade de frittage. Ce sont des sphères monodisperses de rayon $r$ qui sont initialement tangentes les unes aux autres et qui évoluent soit :

- en restant tangentes si les cous (de rayon $x$) croissent sans retrait des cristallites,
- en étant sécantes si un retrait ($\frac{\Delta L}{L_0}$) apparaît.

<figure><img src="image_13.png" /><figcaption>"Les différents déplacements de matière au cours du frittage de deux grains (les chiffres en bleu indiquent les mécanismes de transport de matière)"</figcaption></figure>

FIGURE 2.1 - Premier stade du frittage - Géométrie des cristallites (qualifiées de "grains" dans [DBA2005]) et mécanismes de diffusion [DBA2005]

Les mécanismes mis en jeu sont numérotés sur la figure 2.1, avec, à gauche, ceux ne donnant pas lieu à un retrait et à droite ceux conduisant au retrait :

- l'évaporation/condensation (1) qui ne sera pas considéré dans le modèle développé,

- la diffusion en surface (2),

- la diffusion en volume (de la surface des grains vers le pont : 3 et du joint de grains vers le pont : 4),

- la diffusion intergranulaire (5).

- la déformation plastique qui n'est pas considérée dans notre modélisation du frittage libre et qui amène la matière du s volume vers le pont (6).

### 2.2.1 Loi d'évolution des cous en l'absence de retrait

Des lacunes en sur-concentration au niveau des cous vont migrer en surface vers les zones en compression, et inversement des atomes vont venir faire croître les cous. On suppose que les cristallites restent sphériques. En notant $\delta_s$ l'épaisseur, de quelques nanomètres, de la couche de surface dans laquelle la diffusion a lieu, on obtient alors la loi suivante de croissance des cous [LON2009] :

$$ \frac{x^7}{r^3} = \left( 56\gamma\Theta \frac{D_s \delta_s}{kT} \right) t \qquad (2.1) $$

avec $\gamma$ la tension de surface, $\Theta$ le volume molaire (pour $UO_2 \Theta = 4.09 \ 10^{-23} cm^3$) et $t$ le temps. Si le transport, au lieu d'avoir intégralement lieu en surface, passe par le volume du grain pour aller de la surface au cou, la loi de croissance s'écrit :

$$ \frac{x^5}{r^2} = \left( \frac{5\pi}{2} \gamma \Theta \frac{D_v}{kT} \right) t \qquad (2.2) $$
### 2.2.2 Loi d'évolution des cous avec retrait

Si le transport, a lieu par le volume du grain :

$$ \frac{x^5}{r^2} = \left( 20\pi\gamma\Theta \frac{D_v}{kT} \right) t \qquad (2.3) $$

Si la diffusion intergranulaire est mise en jeu :

$$ \frac{x^6}{r^2} = \left( 96\gamma\Theta \frac{D_{gb}\delta_{gb}}{kT} \right) t \qquad (2.4) $$

La géométrie qu'on se donne en hypothèse (figure 2.1) permet de déduire des lois précédentes un retrait linéaire en remarquant que :

$$ \frac{\Delta L}{L} \approx \frac{x^2}{4r^2} \qquad (2.5) $$

Les lois de retrait s'écrivent :

$$ \left\{ \begin{array}{l} \frac{\Delta L}{L} = \left( 2 \frac{\gamma \Theta D_v}{kT r^3} \right)^{2/5} t^{2/5} \quad (\text{diffusion en volume}) \\ \\ \frac{\Delta L}{L} = \left( 3 \frac{\gamma \Theta D_{gb} \delta_{gb}}{2kT r^4} \right)^{1/3} t^{1/3} \quad (\text{diffusion intergranulaire}) \end{array} \right. \qquad (2.6) $$

## 2.3 Lois du second stade du frittage

Dans le second stade du frittage seules les diffusion en volume et aux joints de grains interviennent. Les grains étant cette fois polyédriques (cuboctaèdres tronqués) comme représenté sur la figure 2.2, la géométrie idéale des modèles de frittage repose plutôt sur le rhombododécaèdre dont les arêtes sont remplacées par de fins cylindres qui représentent la porosité ouverte [RLC1958] tandis que les pores sphériques représentent la porosité fermée.

<figure><img src="image_14.png" /><figcaption>FIGURE 2.2 – Microstructure du deuxième stade de frittage [DBA2005]</figcaption></figure>

### 2.3.1 Diffusion en volume

Chacune des faces du polyèdre est un joint entouré d'un pore cylindrique de rayon $r_p$. On fait l'approximation que ce joint est un disque de rayon $x$. Les atomes quittent le joint et viennent se répartir uniformément sur la surface latérale du cou :

$$ r^n - r_0^n = r_{p0}^n \alpha t \qquad (2.7) $$

avec $r_0$ et $r_{p0}$ les rayons initiaux des grains et des pores, $\alpha$ une fonction de la température et $2 \le n \le 3$.
### 2.3.2 Diffusion aux joints de grains

Les déplacements des atomes s'effectuent dans un disque cylindrique de rayon $x$ et de hauteur $\delta_{gb}$ qui est l'épaisseur du joint de grain :

$$p^{1/2} \frac{dp}{dt} = -28 \frac{\gamma \Theta D_{gb} \delta_{gb}}{kTl^4} \qquad (2.8)$$

avec $p$ la porosité et $l$ la longueur de l'arête des joints de grains qui n'a pas d'expression analytique mais peut être estimée moyennant des hypothèses.

### 2.3.3 Lois de densification

Les lois précédentes peuvent s'exprimer en terme de densité relative $\rho$ :

$$\begin{cases} \rho^{3/2} - \rho_0^{3/2} = k_{gb} \frac{D_{gb}}{r^4} t & (\text{diffusion aux joints de grains}) \\ \rho - \rho_0 = k_v \frac{D_v}{r^3} t & (\text{diffusion en volume}) \end{cases} \qquad (2.9)$$

## 2.4 Lois du troisième stade du frittage

Dans le troisième stade du frittage il ne reste que des pores isolés (porosité fermée) situés soit dans le volume des grains (porosité intragranulaire) soit dans les joints de grains (porosité intergranulaire). Ces pores peuvent disparaître du fait de la diffusion intergranulaire ou volumique ou bien coaslescer. Les grains évoluent en taille, ce qui peut aussi se traduire en terme de mouvement des joints de grains.

### Disparition d'un pore intragranulaire

Un pore intragranulaire ne peut disparaître que du fait de la diffusion en volume avec une vitesse lorsqu'il est dans un grain de taille infinie $\xi_{inf}$ donnée par :

$$\xi_{inf} = \frac{8\pi D_v \gamma_s}{kT} \qquad (2.10)$$

Cette vitesse dépend peu de la distance aux joints de grains et du rapport de la taille du pore à celle du grain. Moyennant certaines hypothèses la cinétique de disparition d'un pore isolé peut s'écrire :

$$r_p^3 - r_{p0}^3 = - \frac{6\gamma_s \Theta D_v}{kT} t \qquad (2.11)$$

### Détachement/épinglement d'un pore

Un pore peut ou pas rester intergranulaire. S'il se détache du joint de grain où il était situé initialement, un pore devient intragranulaire. C'est un point positif si l'on souhaite garder de la porosité en fin de frittage, négatif si l'on souhaite atteindre de fortes densités.

Dans le cas du combustible $(U, Pu)O_2$, par exemple, on veut pouvoir piéger les produits de fission gazeux sans que ce gaz produit n'aille se localiser entre la pastille et la gaine où il créerait une barrière thermique néfaste au refroidissement du combustible. Dans ce cas un pore intragranulaire est plus intéressant pour piéger le gaz car ce dernier devra diffuser en volume (donc très lentement) jusqu'au joint de grain d'où il pourra ensuite gagner rapidement l'espace entre la pastille et la gaine par diffusion aux joints des grains. La séparation des pores des joints de grains augmentent aussi la vitesse des joints de grain et rallonge ainsi le chemin que devra parcourir le gaz de frittage
pour rejoindre le joint de grain.

Si l'on souhaite atteindre de fortes densités un pore intragranulaire sera plus difficile à éliminer
car les mécanismes liés à sa disparition reposeront sur la diffusion en volume bien plus lente que
la diffusion aux joints de grains. De plus, les joints de grains n'étant plus ancrés par les pores
intergranulaires, leur mouvement peut conduire à de la croissance anormale.

## 2.5 Limites des modèles classiques

Les modèles classiques de frittage font l'hypothèse qu'un seul phénomène a lieu à un instant donné alors que dans le second stade du frittage la densification se poursuit en même temps que les grains croissent et que la porosité devient sphérique. Les phénomènes de diffusion intergranulaire et de diffusion en volume interviennent aussi bien dans le premier que dans le second stade du frittage. Sur certains matériaux ces deux stades sont même confondus. De même, dans le second stade, certaines études [DLJ1970] ont nécessité de prendre en compte simultanément la diffusion en volume et celle aux joints de grains.

## 2.6 Conclusion

Du fait des limites des modèles classiques de frittage nous avons proposé un modèle qui ne fait
pas d'hypothèse sur un mécanisme prépondérant à un instant donné et ce modèle est à l'échelle
sub-granulaire.

Un second objectif initial était que les paramètres mis en jeu par le modèle puissent être
mesurés expérimentalement par des expériences spécifiques dédiées ne faisant pas intervenir le
frittage.
# Maîtrise de la microstructure - Cartes de microstructure

## 3.1 Définition

On a vu que les phénomènes élémentaires de diffusion à l'état solide pouvaient servir de base à la compréhension de l'évolution de la microstructure (taille de grains et de pores, type inter ou intra-granulaire pour les pores). Il est commode, en fonction des buts à atteindre (présence de pores intragranulaires par détachement des pores intergranulaires...), de choisir, de représenter les phénomènes microstructuraux ayant lieu lors du frittage en fonction des deux variables taille des grains et taille des pores. A cette fin on peut dresser une carte de microstructure (ou carte de frittage) qui délimite les zones de prépondérance des phénomènes qui contrôlent l'évolution de la microstructure (croissance des grains ...). On peut alors, représenter l'effet d'un adjuvant de frittage sur une telle carte comme le montre la figure 3.1 :

<figure><img src="image_15.png" /><figcaption>FIGURE 3.1 – Influence sur une solution solide d'un ajout de 1% lorsque les pores se déplacent par diffusion superficielle [DBA1993]</figcaption></figure>
Ashby [MFA1974] avait réalisé de telles cartes surtout pour le premier stade du frittage en mentionnant les processus dominant la croissance des cous. Les lignes dans ces diagrammes délimitent les conditions expérimentales pour lesquelles deux mécanismes contribuent à part égale à la vitesse de croissance des cous. Ces diagrammes, dont un exemple est illustré en figure 3.2, montrent aussi les contours iso-durée de frittage avec les valeurs des cous à différentes températures.

<figure><img src="image_16.png" /><figcaption>FIGURE 3.2 – Diagramme de frittage d'un agrégat de sphères de 38 µm d'Ag pur [SJK2005]</figcaption></figure>

Pour $UO_2$ Mohan et al. [AM1979] ont développé un logiciel pour construire automatiquement ces cartes de frittage sur la base des équations de grossissement de cou de chaque mécanisme pris isolément.

## 3.2 Exemple d'utilisation

Dans le cadre de la thèse d'Erick Girard (cf Annexe C page 133), que j'avais co-encadrée de 2001 à 2004 avec le Laboratoire de Thermodynamique et Physicochimie Métallurgiques (J.M. Chaix) et l'École des Mines de Saint-Étienne (P. Goeuriot, F. Valdivieso), le frittage de $UO_2$ avait été étudié avec des hétérogénéité de densité (inclusions de $UO_2$ cru tamisé-forcé - dits "amas" - notés $am_i$, agrégats pré-frittés de $UO_2$ - notés $ag_i$, avec $i = 1$ pour des échantillons prélevés à 1600°C lors de la montée en température, $i = 2$ pour ceux prélevés au bout de 15 minutes sur le palier à 1700°C et $i = 3$ pour ceux prélevés au bout de 2h sur le palier à 1700°C et pores générés avec des billes calibrées d'une dizaine de micromètres initialement en latex). Après attaque thermique de révélation des grains, des cartographies étendues autour de ces défauts sur une distance de 40µm avaient été acquises en "recollant" des images prises au Microscope Électronique à Balayage à effet de champ (MEB FEG du Laboratoire Bernard François de Cadarache, puis MEB FEG de l'École des Mines de Saint-Étienne). Ces cartographies de grandes dimensions avaient été traitées à l'aide de macro-commandes Aphelion au SIMAP par E. Girard et J.-M.Chaix afin de tracer les distributions surfaciques des pores intra et inter-granulaires ainsi que la taille des grains en fonction de l'abscisse (négative dans le défaut et positive à l'extérieur).
rieur du défaut) l'origine étant prise sur la frontière externe du défaut dans le plan d'observation.

De là, on a pu tracer des cartes 3D de frittage avec pour axes la taille des grains, la taille des pores intergranulaires et le taux de porosité. Des tendances ont pu se dégager avec une répartition des amas sur une nappe tandis que les agrégats se répartissent sur une autre nappe (figure 3.3). De plus, les défauts et leur voisinage ont été observés à différents stades du frittage, de sorte que l'évolution des groupes de points dans ces axes a pu être suivie en fonction du temps (figure 3.4).

On retrouve, par exemple, que les amas sont localisés sur une même nappe avec moins de porosité résiduelle car ces amas ont été tamisés-forcés et ne contenaient, ainsi, vraisemblablement, initialement, que très peu de porosité, tandis que les agrégats contenaient de la porosité et peu de surface de contact entre leurs particules ce qui a conduit à une porosité résiduelle non négligeable.

## 3.3 Conclusion

On peut retenir que la constitution de cartes de frittage, qu'elles soient 2D ou 3D, est potentiellement utile comme une aide à l'interprétation.

<figure><img src="image_17.png" /><figcaption>FIGURE 3.3 – Regroupement des échantillons sur des nappes de frittage - thèse E. Girard</figcaption></figure>
<figure><img src="image_18.png" /><figcaption>FIGURE 3.4 – Trajectoires de frittage - thèse E. Girard</figcaption></figure>
# Le frittage à l'échelle des milieux continus

Reid et al. [CRR1990], [CRR1994] traduisent le comportement d'un objet macroscopique lors du frittage en loi de comportement mécanique :

$$ \bar{\bar{T}} = \left\{ -p_0 - \sigma \left[ \bar{H} - S_\nu + \nu \frac{\partial S_\nu}{\partial \nu} \right] + \lambda \operatorname{tr}(\bar{D}) \right\} \bar{I} + 2\mu \bar{\bar{D}} \quad (4.1) $$

où $\bar{\bar{T}}$ est le tenseur des contraintes, $\bar{\bar{D}}$ celui des déformations, $p_0$ la pression appliquée lors du frittage (pour un frittage dit « libre » $p_0 = 1 \text{ atm}$), $\sigma$ la tension de surface de l'interface grain/gaz (en général de l'ordre $1N/m$), $\bar{H}$ la courbure moyenne de l'interface, $\nu$ le rapport de la densité à la densité « théorique », $S_\nu$ l'aire de l'interface solide/gaz pour une densité relative donnée $\nu$, $\lambda$ et $\mu$ les coefficients de Lamé du matériau à la fin de la densification.

Cette approche est spécifique aux matériaux monophasés, son avantage est la simplicité du calcul des contraintes ou des déformations. L'inconvénient est de savoir, à un instant donné de la transformation, quelles sont les distributions de taille de porosité et de rayon de courbure. Cependant, lorsque ces deux derniers paramètres varient peu l'approche de Reid peut être utile.

## 4.1 Notion de contrainte de frittage

La contrainte de frittage est définie comme la contrainte nécessaire à empêcher le frittage spontané du matériau qui fritte, i.e., la valeur de $\bar{T}$ telle que $\bar{\bar{D}} = 0$. C'est une grandeur qui permet de situer la propension d'un matériau à évoluer spontanément vers un solide plus dense lors d'un traitement thermique. Elle peut permettre de classer des matériaux entre eux en terme de frittabilité.

## 4.2 Cas d'un matériau biphasé - contrainte différentielle

Comme $UO_2$ et $(U, Pu)O_2$ présentent des vitesses de retrait au frittage différentes et que ces deux matériaux forment une solution solide continue, leur traitement thermique va être tel que la jonction entre les deux matériaux sera très cohésive (même solution solide) mais un champ de contrainte va se développer dans les pastilles lors du frittage du fait de la différence de cinétique de retrait des deux matériaux (frittage différentiel).

Dans le cadre de la thèse de Grégoire Toury de 1999 à 2001, (cf Annexe C page 133), j'ai participé
à la mise en œuvre d'un modèle de retrait de matériaux biphasés relevant de la mécanique des milieux continus [AT1995]. Il a ainsi été mis en œuvre dans le cas d'un composé céramique-céramique (CERCER) constitué d'une matrice $UO_2$ avec des amas préfrittés de $(U, Pu)O_2$ afin d'estimer en contraintes planes l'ordre de grandeur des contraintes lors du frittage [AT1995], [GT2001]. En toute rigueur les deux compositions $UO_2$ et $(U, Pu)O_2$ correspondent à une seule phase du point de vue de la thermodynamique, toutefois leur paramètre de maille (et éventuellement leur coefficient de dilatation thermique) est légèrement différent et le coefficient d'interdiffusion $U/Pu$ est très faible. Typiquement la distance d'interdiffusion est de l'ordre de $5\mu m$ à $1700°C$ pour une durée de 4h. On peut considérer que les deux constituants du CERCER resteront inchangés si les inclusions sont de dimension assez grandes devant cette distance d'interdiffusion. En général c'est le cas pour les fabrications liées aux REP : les inclusions de $(U, Pu)O_2$ sont de l'ordre de $20$ à $100 \mu m$ de diamètre.

L'étude expérimentale a révélé que les retraits différentiels mis en jeu pouvaient, dans le cas d'un frittage suffisamment retardé des inclusions plutonifères, provoquer une mise en compression des amas par la matrice pour aller jusqu'à la fissuration des amas. La simulation 2D en contraintes planes a eu pour but d'essayer d'évaluer le champ de contraintes dans ce matériau. La démarche a consisté à calculer les déformations élasto-viscoplastiques relatives à chacun des constituants du composite, puis à appliquer les lois de l'élasticité à la composante élastique pour estimer les contraintes. La déformation relative globale s'écrivant alors :

$$ \epsilon_{rs} = \epsilon_{rs}^{\text{frittage}} + \epsilon_{rs}^{\text{dil}} + \epsilon_{rs}^{\text{el}} + \epsilon_{rs}^{\text{viscoplast}} \quad (4.2) $$

avec *dil* pour dilatation thermique, *el* pour élastique et *viscoplast* pour viscoplastique, *r* et *s* désignent les directions d'espace. La déformation $\epsilon_{rs}^{\text{frittage}}$ due au frittage du matériau est assimilée, localement, au retrait d'un des constituants du composite pris isolément, et, de plus, supposée localement isotrope :

$$ \epsilon_{rs}^{\text{frittage}} = \frac{\Delta l}{l_0} \delta_{rs} \quad (4.3) $$

avec $\delta_{rs}$ le symbole de Kronecker. La composante viscoplastique de déformation a été calculée à l'aide d'un paramètre de fluidité $\gamma$ et d'un potentiel viscoplastique de dissipation $\Omega$ estimés :

$$ \dot{\epsilon}^{\text{viscoplast}} = \gamma \langle \Phi(F) \rangle \frac{\partial \Omega}{\partial \sigma} \quad (4.4) $$

Le matériau n'étant pas ductile, le critère de Tresca a été choisi, pour lequel :

$$ \langle \Phi(F) \rangle = \frac{2\bar{\sigma}\cos(\theta) - C_Y}{C_Y} \quad (4.5) $$

On a choisi un modèle de plasticité associée, i.e. $\Omega = F$, tandis que $\gamma$ a été estimé à partir de la connaissance de la loi de fluage. L'équilibre des forces amène alors au problème suivant :

$$ \left\{ \begin{array}{l} \sigma_{rs} = \lambda u_{q,q} \delta_{rs} - (3\lambda + 2\mu) \left(\alpha \Delta T + \frac{\Delta l}{l_0}\right) \delta_{rs} - \lambda \epsilon_{qq}^{\text{fluage}} \delta_{rs} - 2\mu \epsilon_{rs}^{\text{fluage}} + \mu (u_{r,s} + u_{s,r}) \\ \sigma_{rs,s} = -\rho b_r \end{array} \right. \quad (4.6) $$

Le système 4.6 est résolu sur les deux sous-domaines $\Omega_1$ (matrice $UO_2$) et $\Omega_2$ (amas plutonifères). En notant $\partial\Omega_{1/2}$ la frontière entre ces deux sous-domaines, $\partial\Omega_1$ la frontière entre le sous-domaine $\Omega_1$ et l'extérieur et $\partial\Omega_2$ la frontière entre le sous-domaine $\Omega_2$ et l'extérieur. Les conditions aux limites sont :

$$ \text{pour } \partial\Omega_{1/2} \left\{ \begin{array}{l} \vec{u}_{\Omega_1} = \vec{u}_{\Omega_2} \quad (\text{continuite des deplacements}) \\ \sigma_{\Omega_1}(\vec{n}) = \sigma_{\Omega_2}(\vec{n}) \quad (\text{equilibre des forces}) \end{array} \right. \quad (4.7) $$
$$
\text{pour } \partial\Omega_1 \cup \partial\Omega_2 \left\{
\begin{array}{ll}
\text{a la base de la pastille qui repose sur son support} & : \quad \begin{cases} \sigma(\vec{n}) \cdot \vec{n} & \le 0 \\ \vec{u} \cdot \vec{n} & \le 0 \\ (\vec{u} \cdot \vec{n}) \sigma(\vec{n}) & = 0 \end{cases} \\
\text{sur le reste du contour de la pastille} & : \quad \sigma_{rs}(\vec{n}) = -p\delta_{rs}\vec{n}
\end{array}
\right.
\tag{4.8}
$$

La microstructure a été simulée comme étant une distribution d'ellipsoïdes (représentant les amas plutonifères) en 3D dans un cylindre (la matrice de $UO_2$). Les amas sont distribués aléatoirement avec une densité (nombre d'amas par unité de volume) qui suit une loi de Poisson. La microstructure biphasée amas/matrice est alors générée suivant un schéma booléen : à chaque centre implanté de façon aléatoire est associé un amas de forme ellipsoïdale de manière à prendre en compte la déformation observée due au pastillage ou à l'étape de granulation. La valeur du rapport du petit axe au grand axe de l'ellipsoïde suit une distribution lognormale observée expérimentalement. On ne garde alors qu'une coupe 2D de la microstructure simulée (figure 4.1), pour ensuite reconnaître les contours de ses amas, générer un fichier de directives de maillage pour le mailleur netgen [JS1997] et finalement mailler avec des éléments P1 en 2D (figure 4.2).

<figure><img src="image_19.png" /><figcaption>FIGURE 4.1 – Représentation du plan de coupe transversal d'une pastille composite crue modélisée</figcaption></figure>

Pour chacun des composites simulés, on a calculé en chaque point de l'interface amas/matrice
les contraintes normale et tangentielle qui sont ensuite moyennées le long de ces mêmes interfaces
(cf tableau 4.1) :

$$
\langle \sigma \rangle = \frac{\sum_{\text{tous les amas}} \left( \int_{C_i} \sigma^* ds \right)}{\sum_{\text{tous les amas}} L_i} \quad (4.9)
$$

Bien que les résultats soient très qualitatifs eu égard aux nombreuses hypothèses faites, on
peut trouver des tendances, par exemple le fait que les amas sont d'autant plus mis en compression
par la matrice qu'ils ont été pré-frittés (tableau 4.1) et que le l'intensité des contraintes de von
Mises dans la matrice, cartographiée dans la figure 4.3, augmente lorsque la distance moyenne
entre amas diminue.
<figure><img src="image_20.png" /><figcaption>FIGURE 4.2 – Distribution en taille du pas *h* de maillage et zoom sur le maillage, en vert : la matrice, en rouge : les amas de $(U, Pu)O_2$</figcaption></figure>

<table><caption>TABLE 4.1 – Contraintes aux interfaces amas-matrice à 1700°C</caption><thead><tr><th>Température de pré-frittage</th><th>20°C</th><th>1000°C</th><th>1200°C</th><th>1400°C</th></tr></thead><tbody><tr><td>retrait linéaire matrice %</td><td>-14.6</td><td>-14.6</td><td>-14.6</td><td>-14.6</td></tr><tr><td>retrait linéaire amas %</td><td>-13.1</td><td>-11.4</td><td>-8.7</td><td>-3.0</td></tr><tr><td>retrait linéaire global expérimental %</td><td>-14.4</td><td>-14.2</td><td>-13.5</td><td>-10.2</td></tr><tr><td>retrait linéaire global simulé %</td><td>-14.0</td><td>-13.7</td><td>-13.2</td><td>-12.2</td></tr><tr><td>σ<sub>N</sub></td><td>-5</td><td>-9</td><td>-32</td><td>-60</td></tr><tr><td>σ<sub>t</sub></td><td>12</td><td>18</td><td>28</td><td>50</td></tr></tbody></table>

## 4.3 Conclusion

Les modèles mécaniques macroscopiques peuvent aider à dégager des tendances à l'échelle de la pastille, tendances qui pourraient ne pas être dégagées à l'aide d'une simple loi des mélanges lorsque, par exemple, les amas de mélange-mère de la pastille sont proches les uns des autres et interfèrent.
<figure><img src="image_21.png" /><figcaption>FIGURE 4.3 – Simulation de l’effet de la taille des amas sur l’intensité des contraintes de von Mises au sein d’un composite CERCER UO₂ – (U, Pu)O₂ fritté à 1700°C sous atmosphère réductrice (fraction volumique constante de 20%)</figcaption></figure>
5

Conlusion

Dans cette première partie les phénomènes élémentaires qui sont à l'origine du frittage ont été rappelés avec leurs conséquences sur la microstructure à l'échelle sub-granulaire (évolution du diamètre des ponts entre particules rapprochement du centre de gravité des grains). Ces phénomènes élémentaires peuvent être concomitants selon les rampes de température ou, pour les oxydes, aussi selon la pression partielle d'oxygène.

Les domaines de prépondérance de ces phénomènes peuvent être délimités dans des représentations bi ou tridimensionnelles dont les axes sont des caractéristiques microstructurales telles la taille de grains, la fraction volumique de porosité et la densité.

Dans la pratique il est rare qu'une plage de conditions expérimentales ne présente qu'un unique phénomène prépondérant. Aussi, afin de pouvoir restituer l'évolution de cette micro-structure lors du frittage, on a développé un modèle à l'échelle sub-granulaire dans lequel on a cherché à rendre compte de chacun de ces phénomènes. Certains, tels l'évaporation-condensation et l'interdiffusion chimique, pris en compte à l'origine du modèle ont été temporairement laissés de côté afin de valider, dans un premier temps, le modèle par comparaison aux modèles classiques analytiques.

La seconde partie qui va suivre souligne les spécificités physico-chimiques des actinides à prendre en compte pour rendre compte numériquement de leur frittage.
# Deuxième partie

## Les spécificités des oxydes d'actinides
<figure><img src="image_22.png" /></figure>
# 1

# Le caractère radioactif

## Sommaire

<table><tr><td>1.1</td><td>La décroissance des isotopes des éléments majoritaires U et Pu</td><td>37</td></tr><tr><td>1.1.1</td><td>Vecteur isotopique typique du combustible vierge</td><td>37</td></tr><tr><td>1.1.2</td><td>Isotopes de U dans les combustibles</td><td>38</td></tr><tr><td>1.1.3</td><td>Isotopes de Pu dans les combustibles</td><td>38</td></tr><tr><td>1.2</td><td>La décroissance des isotopes des actinides mineurs Am, Np et Cm</td><td>39</td></tr><tr><td>1.3</td><td>La thermodynamique pour les isotopes radioactifs</td><td>39</td></tr><tr><td>1.3.1</td><td>Les fonctions d'état</td><td>39</td></tr><tr><td>1.3.2</td><td>1<sup>er</sup> principe de la thermodynamique</td><td>40</td></tr><tr><td>1.3.3</td><td>2<sup>nd</sup> principe de la thermodynamique</td><td>40</td></tr><tr><td>1.3.4</td><td>3<sup>eme</sup> principe de la thermodynamique</td><td>41</td></tr><tr><td>1.3.5</td><td>Conclusion</td><td>42</td></tr></table>

Une partie majeure de l'interprétation des phénomènes irréversibles se produisant lors du frittage repose sur la thermodynamique linéaire des processus irréversibles. Aussi, le caractère radioactif des matériaux est à considérer car sa comptabilité avec la thermodynamique classique n'est pas évidente. Cette compatibilité dépend de la période des isotopes et donc du vecteur isotopique des éléments chimiques mis en jeu. C'est à ce titre que les périodes et les modes de décroissance des isotopes des éléments chimiques majoritaires du combustible vierge sont rappelés. Cet aspect serait à aborder d'une manière bien plus cruciale si l'on s'intéressait au re-frittage du combustible en réacteur.

## 1.1 La décroissance des isotopes des éléments majoritaires U et Pu

### 1.1.1 Vecteur isotopique typique du combustible vierge

L'élément plutonium, produit à partir d'un combustible vierge (i.e. non irradié) $UO_2$, présente un rapport entre ses isotopes (leur nombre de masse variant de 238 à 242) particulier en fonction de l'irradiation que le combustible dont il est issu a subi en réacteur (taux de combustion) et de la durée de son refroidissement après déchargement du réacteur. Lorsque ce plutonium sera recyclé pour être à son tour incorporé à du combustible vierge (combustible $(U, Pu)O_2$ - mélange d'oxydes d'Uranium et de Plutonium), il aura en sortie de réacteur et après refroidissement un
nouveau rapport entre ses isotopes; un nouveau "vecteur isotopique". Le taux de combustion dans un réacteur est ajusté en fonction des bilans de plutonium et d'actinides mineurs avec des contraintes de sûreté et de non-prolifération. Un vecteur isotopique typique pour les combustibles vierges de première et seconde génération REP (respectivement à 3,3% de Pu irradié à $33GWjt^{-1}$ et irradié à $45GWjt^{-1}$) refroidis 7 ans est donné dans le tableau 1.1 [TA2007] :

<table><caption>TABLE 1.1 – Un vecteur isotopique typique du $(U, Pu)O_2$ vierge [TA2007]</caption><thead><tr><td>Isotope</td><td>% 33GWjt<sup>-1</sup></td><td>% 45GWjt<sup>-1</sup></td></tr></thead><tbody><tr><td>238</td><td>1,4</td><td>2,6</td></tr><tr><td>239</td><td>60,0</td><td>54,5</td></tr><tr><td>240</td><td>21,7</td><td>28,1</td></tr><tr><td>241</td><td>11,8</td><td>8,0</td></tr><tr><td>242</td><td>5,2</td><td>6,8</td></tr></tbody></table>

### 1.1.2 Isotopes de U dans les combustibles

Les isotopes de l'uranium utilisés lors de la fabrication du combustible vierge sont $^{235}$U pour le combustible $UO_2$ et $^{238}$U pour le $(U, Pu)O_2$. Certains pays, à l'issue le retraitement de l'uranium, utilisent aussi $^{236}$U [TA2007].

$^{235}$U donne par décroissance $\alpha$ du $^{231}$Th. L'énergie de la particule $\alpha$ émise est de 4,6787 MeV et sa période de 710 millions d'années. Sa faible radioactivité, additionnée à celle de ses descendants dans sa chaîne de désintégration, génère une puissance de 0,1 Watt par tonne. Tout se passe comme si c'était un isotope stable. Son schéma de décroissance est donné en annexe G figure G.1.

$^{238}$U conduit, par décroissance $\alpha$ à l'isotope $^{234}$Th en émettant une particule $\alpha$ de 4,270 MeV avec une période de 4,5 milliards d'années (figure 1.1 ou G.2 en annexe G).

### 1.1.3 Isotopes de Pu dans les combustibles

Les isotopes de Pu mis en œuvre dans les combustibles sont les isotopes 238, 239, 240, 241 et 242. Contrairement aux isotopes de l'uranium, ils ont une période plus courte qui fait que la puissance thermique de certains [DWO1976] est non négligeable (tableau 1.2). Leurs modes de décroissance ainsi que leur période sont donnés figures G.3, G.4, G.5, G.6 et G.7.

<table><caption>TABLE 1.2 – Puissance thermique des isotopes de Pu</caption><thead><tr><td>isotope</td><td>puissance thermique (µWg<sup>-1</sup>)</td><td></td></tr></thead><tbody><tr><td><sup>238</sup>Pu</td><td>0,5672</td><td>10<sup>+6</sup></td></tr><tr><td><sup>239</sup>Pu</td><td>1,9293</td><td>10<sup>+3</sup></td></tr><tr><td><sup>240</sup>Pu</td><td>7,097</td><td>10<sup>+3</sup></td></tr><tr><td><sup>241</sup>Pu</td><td>3,390</td><td>10<sup>+3</sup></td></tr><tr><td><sup>242</sup>Pu</td><td>115,94</td><td>10<sup>+3</sup></td></tr></tbody></table>

On voit que cet élément possède des isotopes dont les 238 et 241 de demi-vie respective 86a et 14,35a sont bien moins stables que ceux de l'uranium.
<figure><img src="image_23.png" /><figcaption>FIGURE 1.1 – Période et modes de décroissance de $^{238}$U d’après [ENSDF2006]</figcaption></figure>

## 1.2 La décroissance des isotopes des actinides mineurs Am, Np et Cm

Pour les combustibles du futur certains isotopes de Am, Np et Cm seront incorporés au combustible : $^{237}$Np ($T = 2,14 \cdot 10^6$ a), $^{241}$Am ($T = 432$ a), $^{242m}$Am ($T = 141$ a), $^{243}$Am ($T = 7370$ a), $^{242}$Cm ($T = 162,8$ j, $P = 120$ W g$^{-1}$) et $^{243}$Cm ($T = 29,1$ a). Leur décroissance et/ou celle de leurs descendants donne lieu aussi à une émission $\gamma$. Les figures G.8, G.9 et G.10 donnent les modes de décroissance des isotopes de l’américium.

La puissance thermique de $^{241}$Am d’une valeur de $114,3 \cdot 10^{+3} \mu$ W g$^{-1}$ est comparable à celle de $^{242}$Pu.

## 1.3 La thermodynamique pour les isotopes radioactifs

### 1.3.1 Les fonctions d’état

Il n’est pas possible de définir une fonction d’état pour un système constitué d’atomes présentant une décroissance radioactive. Un tel système ne peut jamais revenir à son état initial. Pour cela, il faudrait que la transformation qu’il subisse ait lieu instantanément, ce qui exclut, en particulier, les changements de phases.

Si l’on voulait étendre la définition de telles fonctions, il faudrait en définir qui, lorsque la
période de l'isotope de l'élément chimique tend vers l'infini, tendent vers les fonctions thermodynamiques usuelles, afin de disposer d'un cadre cohérent pour faire des calculs. Et dans ce cadre, on pourrait ensuite faire des calculs d'incertitude liées à la durée des expériences par rapport aux périodes des éléments pour évaluer l'impact sur les résultats des calculs. Les logiciels de calcul d'équilibre seraient aussi à adapter.

Comme les durées de frittage sont généralement courtes vis-à-vis des périodes des isotopes considérés que l'on met en œuvre dans le combustible vierge et que l'on ne dispose, de toutes façons, pas d'autre cadre théorique, on utilisera dans la suite la thermodynamique classique malgré ses limites pour des matériaux radioactifs.

### 1.3.2 1er principe de la thermodynamique

Pour un système présentant une décroissance radioactive, si la masse était conservée il n'y aurait pas d'énergie libérée. Il faudrait définir le système fermé comme contenant les descendants, les neutrons et particules alpha émis et les photons γ émis. Il faudrait décrire le système en descendant aux particules élémentaires, le premier principe serait alors encore valable, mais pas en ne considérant que les atomes pères et fils, c'est-à-dire qu'il ne s'appliquerait pas à la thermochimie.

On illustre cette conséquence sur deux exemples : $^{241}$PuO$_2$ et $^{244}$CmO$_2$.

La décroissance d'une mole de $^{241}$PuO$_2$, selon le schéma représenté en figure 1.2 suit une

<figure><img src="image_24.png" /><figcaption>FIGURE 1.2 – Schéma de décroissance de $^{241}$Pu</figcaption></figure>

loi $N = N_0e^{-\lambda t}$ avec $\lambda = \frac{\ln(2)}{\tau_{1/2}}$ avec $\tau_{1/2} = 14,35$ ans. Si l'on veut que cette énergie soit inférieure à 5% de $\Delta H$ de fusion de PuO$_2$, il faut que $\Delta E \cdot (N_0 - N) \le 0,05 \cdot 70000$ avec $\Delta E = 0,0208$ MeV = $2 \cdot 10^9$ J·mol$^{-1}$. Il faut donc que la durée des phénomènes mis en jeu et dont on veut mesurer l'enthalpie durent moins de 794s. Ceci ne représente pourtant qu'une faible variation en nombre de $^{241}$Pu : $\frac{\Delta N}{N_0} = -1,21 \cdot 10^{-6}$.

Pour $^{244}$CmO$_2$ avec $\tau_{1/2} = 18,10$ ans (figure 1.3) et $\Delta H^0_{fusion} = 45,6$ kJ·mol$^{-1}$ on obtient une durée $\Delta t$ inférieure à 66s.

### 1.3.3 2nd principe de la thermodynamique

Toute transformation qui se produit en un laps de temps non nul, est irréversible (puisqu'il y a eu pendant ce temps une décroissance radioactive).
Donc la relation :

$$dS = \frac{\delta Q_{reversible}}{T} \qquad (1.1)$$
<figure><img src="image_25.png" /><figcaption>FIGURE 1.3 – Schéma de décroissance de $^{244}$Cm</figcaption></figure>

où $\delta Q_{reversible}$ serait la quantité de chaleur reçue par le système (constitué du matériau radioac- tif) lors d'une transformation réversible, ne peut pas servir à définir une entropie.

### 1.3.4 3ème principe de la thermodynamique

À la température de 0K du fait de la décroissance radioactive une partie des atomes de dé- part sont remplacés par leur(s) descendant(s), ainsi le désordre augmente donc l'entropie usuelle, même si elle était nulle à un instant donné ne pourrait pas l'être dans les instants ultérieurs.

On l'illustre sur la cas de $^{241}$Pu, pour lequel $^{241}$Pu est substitué, du fait de sa décroissance, par $^{241}$Am. Il se forme ainsi une solution solide de substitution dont l'entropie de configuration peut s'exprimer en fonction de la fraction $x$ de $^{241}$Pu par :

$$ \Delta S_c = -R(x\ln(x) + (1-x)\ln(1-x)) \quad (1.2) $$

qui est maximale à la demi-vie $t = \tau_{1/2} = 14,35$ ans pour laquelle elle atteint la valeur $\Delta S_{cmax} = R\ln(2)$.

L'impact de ce phénomène pour la demi-vie, i.e. quand le combustible a été extrait du réac- teur depuis 14,35 ans, est une contribution de $11,5kJmol^{-1}$ à l'incertitude de l'enthalpie libre à 1750°C, à comparer, par exemple, à l'enthalpie de fusion de $PuO_2$ : 66,5 à 70,3kJmol$^{-1}$.

On peut s'attendre à des comportements analogues pour $^{244}$Cm : $\tau_{1/2} = 18,10$ ans, $^{243}$Cm : $\tau_{1/2} = 28$ ans et dans une moindre mesure pour $^{238}$Pu : $\tau_{1/2} = 86$ ans.

Il en résulte :

- un décalage de l'origine de l'entropie ($S$ à 0K) pour un élément chimique donné qui est variable dans le temps et qui est fonction de l'âge du lot mis en oeuvre (typiquement de l'ordre de la dizaine d'années pour le lot de plutonium actuellement manipulé au laboratoire),
- un décalage de la valeur de l'entropie à 0K variable d'un élément à l'autre, puisque leurs périodes diffèrent,

- l'isotopie serait donc à donner avec les résultats thermodynamiques expérimentaux, et les modèles de frittage, qui s'appuient sur la thermodynamique auront une précision limitée.

## 1.3.5 Conclusion

Il serait possible de traiter le système présentant une décroissance radioactive comme un système ouvert mais il faudrait, en plus, que les états de référence des éléments soient reliés par la différence d'énergie liée à la transformation d'un élément en l'autre (exemple : $^{241}$Pu → $^{241}$Am, hν = 20,78keV). On ne pourrait plus garder la convention usuelle $\Delta H_{formation}^0 = 0$ à 298.15K pour tous les éléments, mais leurs $\Delta H_f^0$ seraient reliés. Par exemple, $\Delta H_{f_{241}Pu}^0$ (298,15K) – $\Delta H_{f_{241}Am}^0$ (298,15K) = $N \cdot 20,78 \cdot 10^3 (1,60217657 \cdot 10^{-19}) = (6,02214129 \cdot 10^{23}) \cdot 20,78 \cdot 10^3 (1,60217657 \cdot 10^{-19}) \approx 2 \cdot 10^9 J \cdot mol^{-1}$. Cette approche conviendrait aussi lors de l'émission d'un noyau α, mais pas d'une réaction de fission spontanée (FS) si les neutrons émis ne sont pas réabsorbés par les autres noyaux de l'échantillon (constituant le système).

Une manière de contourner ces difficultés est de dire qu'on applique la thermodynamique aux isotopes stables des éléments constituant le matériau à fritter.
est variable et, sous leur forme oxyde, ils vont présenter des écarts à la stœchiométrie.

## 2.2 Écarts à la stœchiométrie

L'écart à la stœchiométrie est relié, par la thermodynamique, à la pression partielle d'oxygène, $p_{O_2}$ et à la température $T$. Cette relation repose sur le modèle utilisé pour décrire le potentiel de $O$ en phase solide et, ainsi, sur l'enthalpie libre de la phase solide. Je décris ici un cas particulier sur lequel j'ai travaillé dans le cadre de la thèse de Damien Prieur (cf Annexe C, page 133). La solution solide est décrite avec comme pôles : $UO_2$, $U_{1/3}^1$, $U_4O_9$ ou $U_3O_7$, $AmO_2$, $Am_{5/4}O_2$ et $U_{4/3}Am_{4/3}O_{16/3}$ caractérisés thermodynamiquement par une enthalpie libre de formation $\Delta G_f^0(T)$ correspondant à l'enthalpie libre de la réaction qui amène $O_2(g)$ et le métal $M$ à l'infini l'un de l'autre vers la position qu'occupent les ions correspondants ($O^{-2}$, $M^{3+,4+,5+}$, avec $M = U$ ou $Am$) à leur position d'équilibre dans le cristal à la température considérée. Ces enthalpies libres sont des fonctions affines de la température ($A + BT$). Une première approximation consiste à considérer le mélange mécanique simple qui conduit à une contribution à l'enthalpie libre de la forme :

$$ \Delta G_f^0(T) = \sum_{i=1}^{6} x_i \Delta G_{f_{pole} n^\circ i}^0(T) \qquad (2.1) $$

avec $x_i$ la fraction molaire du pôle $n^\circ i$ (avec $\sum_{i=1}^6 x_i = 1$). Vient ensuite se greffer le fait que la présence de plusieurs constituants amène plusieurs possibilités d'arrangements et ainsi une entropie de configuration, qui va aussi contribuer :

$$ \Delta G_f^0(T) = \sum_{i=1}^{6} x_i \Delta G_{f_{pole} n^\circ i}^0(T) + RT \sum_{i=1}^{6} x_i L n(x_i) \qquad (2.2) $$

<footnote n="1">1. Ce pôle a été choisi par analogie aux travaux de Lindemer et Besmann sur $UO_{2\pm x}$ [TLI1985] et $(U, Pu)O_{2-x}$ [TBE1985]. Ces derniers ont calculé l'expression du potentiel d'oxygène en fonction de l'écart à la stœchiométrie $x$ dans l'hypothèse d'une solution régulière de paramètre d'interaction $E$ avec les couples de pôles suivants pour le domaine sous-stœchiométrique :</footnote>

<footnote>$UO_2 - U_2O_3$ équilibre de formation à partir du pôle $UO_2$ s'écrit : $4[UO_2] = 2 [U_2O_3] + O_2$ réaction (1)</footnote>

<footnote>$UO_2 - U_{4/3}O_2$ équilibre : $4[UO_2] = 3[U_{4/3}O_2] + O_2$ réaction (2)</footnote>

<footnote>$UO_2 - UO_{3/2}$ équilibre : $4[UO_2] = 4[UO_{3/2}] + O_2$ réaction (3)</footnote>

<footnote>Et d'un manière analogue avec un pôle au degré d'oxydation +II :</footnote>

<footnote>$UO_2 - UO$</footnote>

<footnote>$UO_2 - U_{2/3}O_{2/3}$</footnote>

<footnote>$UO_2 - U_{1/2}O_{1/2}$</footnote>

<footnote>puis au degré d'oxydation 0 :</footnote>

<footnote>$UO_2 - U_{1/4}$</footnote>

<footnote>$UO_2 - U_{1/3}$</footnote>

<footnote>$UO_2 - U_{1/2}$</footnote>

<footnote>Les courbes potentiel d'oxygène en fonction de $x$ sont différentes dans les cas (1), (2) et (3) par exemple. La meilleure adéquation a été obtenue pour $U_{1/3} - UO_2$. Comme ils interprètent aussi la solution solide comme une espèce $U$ au degré d'oxydation 0 diluée dans $UO_2$ de structure fluorine, l'ajustement des constantes du modèle leur permet d'obtenir l'enthalpie libre de formation de ce qu'ils notent $[U_{1/3}]CFC$ qui représente l'enthalpie libre de formation de la solution solide $UO_2$ extrapolée à la composition $U_{1/3}$ avec une structure CFC (le sens cristallographique de cette description est limité, mais il signifie seulement la présence d'atomes $U$ dans la structure pour atteindre la sous-stœchiométrie. Aucune étude n'a été faite pour le démontrer, peut-être des expériences de mesures de conductivité électrique pourraient y aider.</footnote>
Cela correspond au cas où les constituants n'interagissent pas entre eux. Quand ce n'est plus le cas, que la solution n'est plus idéale mais régulière ou sub-régulière, viennent s'ajouter des termes d'interaction. Dans une première approximation ces interactions n'interviennent que sur des paires $(i, j)$ de pôles (interactions de paires). Des termes vont alors s'ajouter à l'expression de l'enthalpie libre, et ce sont ces termes que l'on a ajustés à partir des données expérimentales. On détaille l'expression de ces termes qui sont les polynômes standard en thermodynamique dits de Redlich-Kister-Muggianu (RKM) :

$$x_i x_j E_{ij} = x_i x_j (E_{ij}^0 (x_i - x_j)^0 + E_{ij}^1 (x_i - x_j)^1 + \dots + E_{ij}^n (x_i - x_j)^n) \quad (2.3)$$

Chacun des termes $E_{ij}^k$ dépend lui-même de la température selon une loi en $A + BT$. L'adéquation aux données expérimentales a nécessité d'ajuster les termes d'ordre 0 et 1 pour les interactions entre les pôles $Am_{5/4}O_2$ et $U_{4/3}Am_{4/3}O_{16/3}$ ainsi qu'entre $UO_2$ et $U_{4/3}Am_{4/3}O_{16/3}$. De plus l'enthalpie libre de formation du pôle $U_{4/3}Am_{4/3}O_{16/3}$ qui est cité comme pôle dans la littérature mais n'existe pas comme composé défini dans la nature a du être estimée.

La relation entre les valeurs du rapport $Am/(U + Am)$, celle de l'écart à la stœchiométrie $x$ ($x = O/(Am + U)$) et le potentiel d'oxygène (dérivée de l'enthalpie libre précédente par rapport au nombre de moles de $O$ en phase solide) sont complexes et obtenues numériquement à l'aide d'un logiciel de thermochimie minimisant l'enthalpie libre du système avec contrainte de conservation de la matière pour chaque élément. Pour ajuster les constantes, on a écrit un logiciel avec lequel, partant d'un jeu de données pour ces constantes, on déplace le point représentatif de ces données dans leur espace de sorte à minimiser le $\chi^2$ qui somme les carrés des écarts entre le potentiel calculé et le potentiel expérimental sur le jeu de données expérimentales disponible. L'algorithme utilisé est de type Hooke et Jeeves [RH1961], de sorte qu'aucune dérivée partielle de l'enthalpie libre n'a besoin d'être connue.

Le modèle est en bonne adéquation avec les données expérimentales, il s'est même trouvé être correct de manière prédictive sur deux points expérimentaux supplémentaires, mais les spécia listes du domaine ont estimé que la valeur ajustée pour l'enthalpie libre de formation du pôle $U_{4/3}Am_{4/3}O_{16/3}$ était anormalement élevée. On a, au moins, abouti à une piste d'investigation pour la suite de ces études.

## 2.3 Voisinage de la stoechiométrie - défauts ponctuels

Les écarts à la stœchiométrie induisent la présence, dans le matériau, de défauts ponctuels [RC1971]. En sur-stœchiométrie, on trouve des ions oxydes interstitiels (notés $O_i''$ en notation de Kröger-Vink) avec, dans leur voisinage, des cations plus oxydés que ceux du réseau idéal (par exemple un ion $U^{5+}$ sur un site cationique - noté $U_M'$). En sous-stœchiométrie on rencontre des lacunes anioniques (□O) accompagnées de cations réduits ($Pu'_M$). Ces défauts peuvent se regrouper en défauts plus complexes [MK2009], et vont avoir un impact sur de nombreuses propriétés du matériau y-compris sur les coefficients de diffusion.

La majorité de ces défauts étant chargée, la mesure de la conductivité électrique du matériau peut apporter des informations sur leur concentration. Des études de ce type ont été menées dans le département sur le matériau $UO_2$ (thèse de P. Ruello) [PR2001] puis $UO_2$ dopé chrome (thèse de M. Fraczkiewicz [MF2010]).

En vue de mieux connaître les défauts intervenant dans le $(U, Pu)O_2$, le stage de Mickael Grimoin, que j'ai encadré (cf Annexe C, p. 133), avait pour but de modifier un porte-échantillons pour un four en boîte-à-gants du laboratoire afin de pouvoir y faire des mesures de conductivité électrique en température par une méthode 4 points. L'échantillon ayant une géométrie correspondant à
celle que l'on peut aisément fabriquer au laboratoire, de type pastille cylindrique, des simulations numériques ont été nécessaires pour calculer la densité de courant dans la pastille et, de là, sa résistance équivalente. La résolution numérique de l'équation de Poisson qui donne le potentiel électrique a été faite en éléments finis (de type P2) sous la forme d'une application parallèle C++ utilisant la librairie DiffpackTM. La parallélisation des calculs, nécessaire pour ces maillages d'environ 8 millions de nœuds, a consisté à découper les maillages à l'aide du logiciel METIS [GK1995] implémenté dans la "parallel Toolbox" de Diffpack. On a ainsi pu tester l'influence de paramètres tels que l'écartement entre les fils, la largeur de la zone de contact entre les fils et la pastille, ... sur la valeur de la résistance équivalente.

<figure><img src="image_26.png" /><figcaption>FIGURE 2.1 – À gauche : système étudié représenté avec les conditions aux limites. Au dessus : schéma électrique équivalent du système. À droite : situation expérimentale -Stage M. Grimald</figcaption></figure>

Le dispositif avait été testé en inactif (hors boîte-à-gants). Il a été repensé par un post-doctorant (S. Quetel-Weben) de manière à ce que les électrodes servant à faire passer un courant dans l'échantillon présentent une aire de contact avec l'échantillon supérieure à celle des deux autres électrodes.

Pour le modèle à l'échelle sub-granulaire que l'on a développé, le seul impact est que les coefficients de diffusion, les tensions de surface et les constantes élastiques auront une valeur qui sera celle fixée par l'écart à la stœchiométrie sans que cette grandeur n'apparaisse explicitement.

### 2.3.1 Équations de diffusion chimique pour les oxydes

En reprenant les termes de Tannhauser [DST1984], i.e. en appelant diffusion chimique la variation de l'écart à la stœchiométrie en oxygène seule et interdiffusion l'échange de cations entre deux oxydes, la diffusion chimique ne peut pas être traitée à l'aide des équations de Darken :

$$ \left\{ \begin{aligned} D^{Dar} &= N_A D_B^I + N_B D_A^I = (N_A D_B^I + N_B D_A^I) R \frac{d \ln(a_A)}{d \ln(N_A)} \\ R &= 1 + \frac{(1-f) N_A N_B (D_A^{Tr} - D_B^{Tr})^2}{f(N_A D_B^{Tr} + N_B D_A^{Tr}) (N_A D_A^{Tr} + N_B D_B^{Tr})} \end{aligned} \right. \quad (2.4) $$

où $D^I$ est le coefficient de diffusion intrinsèque tandis que $D^{Tr}$ est celui du traceur. En effet, Tannhauser montre que leur utilisation conduit à des aberrations du fait de l'hypothèse de Darken qui s'applique seulement à un système pour lequel le nombre total d'atomes par unité de volume est constant.
Dans le cas de $(U_{1-y}Pu_y)O_{2-x}$, le paramètre de maille, donc le nombre de sites par unité de volume, varie avec x et y, et la relation de Darken ne sera pas appliquable pour une interdiffusion entre $(U_{1-y_A}Pu_{y_A})O_{2-x_A}$ et $(U_{1-y_B}Pu_{y_B})O_{2-x_B}$. Tannhauser propose un traitement correct des profils de diffusion par de tels couples.

### 2.3.2 Simulation de l'interdiffusion lors du frittage

Bien qu'actuellement le modèle à l'échelle subgranulaire ne simule pas la diffusion des éléments U, Pu et O, ce serait un aspect à prendre en compte pour le cas des combustibles contenant plusieurs actinides dont les précurseurs ne seraient pas de composition homogène (cas le plus fréquent aujourd'hui hormis les co-précipités de U et de Pu).
# Spécificités liées à la cristallographie du matériau

## Sommaire

<table><tr><td>3.1</td><td>Structure cristalline</td><td>49</td></tr><tr><td>3.2</td><td>Symétrie des propriétés élastiques</td><td>49</td></tr><tr><td>3.3</td><td>Symétrie des coefficients de diffusion</td><td>50</td></tr><tr><td>3.4</td><td>Symétrie de la tension de surface</td><td>50</td></tr><tr><td>3.5</td><td>Le facettage des grains</td><td>50</td></tr><tr><td>3.6</td><td>Conséquences sur les équations de conservation - spécificité des cristaux ioniques et de la symétrie cubique</td><td>52</td></tr></table>

## 3.1 Structure cristalline

$(U, Pu)O_2$, à haute température ($\ge 800^\circ C$), forme une solution solide continue qui cristallise dans une structure de type fluorine $(CaF_2)$ de groupe d'espace $Fm\bar{3}m$. Bien que le système réticulaire soit cubique et corresponde en général à des propriétés proches de l'isotropie, il s'en écarte pour $UO_2$ en ce qui concerne les propriétés élastiques et de tension de surface. Il est probable que ces spécificités influent sur le premier stade du frittage.

## 3.2 Symétrie des propriétés élastiques

Le tenseur $c_{ijkl}$ des constantes élastiques à 81 composantes se réduit à 21 pour tous les matériaux du fait de la symétrie du tenseur des déformations.

En symétrie cubique, l'invariance du réseau cristallin vis-à-vis des symétries du groupe ponctuel $O_h$ réduit ce nombre à 3 qui sont notées $c_{11}$, $c_{12}$ et $c_{44}$.

Cependant, comme mentionné au paragraphe 1.3.6 p. 65, pour le matériau $(U, Pu)O_2$, la symétrie n'est pas isotrope.

Dans le modèle à l'échelle sub-granulaire de frittage, on a choisi de garder 3 constantes élastiques pour pouvoir représenter le caractère non isotrope.
## 3.3 Symétrie des coefficients de diffusion

Les coefficients de diffusion, sont liés aux fréquences de saut des atomes vers les sites libres voisins de leur position courante [JP1990]. Leur déplacement est thermiquement activé du fait de la barrière d'énergie à franchir pour qu'un atome ou ion passe de sa position courante à une nouvelle position. Ils sont ainsi de la forme :

$$D = D_0 e^{-\frac{E_a}{RT}} \qquad (3.1)$$

Pour le coefficient de diffusion en volume la fréquence de saut pourrait varier selon la direction malgré la symétrie cubique, comme c'est le cas pour certains métaux cubiques centrés (exemple d'un mécanisme à une lacune où les modes de vibration qui rapprochent les plans (111) rap-prochent les atomes des lacunes vers lesquelles ils se dirigent). Mais pour répondre à cette inter-rogation il faudrait d'une part connaître le spectre de phonons du matériau (comme par exemple pour $UO_2$ [ZAK1993]), et d'autre part, connaître le mécanisme de diffusion, ce qui est bien moins le cas, malgré quelques pistes données dans la littérature [SN2012]. Cependant, l'étude du fac-teur de Debye-Waller en diffraction de rayons X en fonction de la température [MC2012] montre que l'amplitude des vibrations est isotrope dans $UO_2$, et ce, de la température ambiante jusqu'à 2000K, on peut donc s'attendre, pour $UO_2$ au moins, à une isotropie des coefficients de diffusion.

<table><caption>TABLE 3.1 – Amplitude des vibrations $U_{ii}$ dans la direction $i$ pour $UO_2$ en fonction de la température</caption><thead><tr><th>T(K)</th><th>481</th><th>2003</th></tr></thead><tbody><tr><td>$U_{11}$</td><td>0,02724</td><td>0,0654</td></tr><tr><td>$U_{22}$</td><td>0,02724</td><td>0,0654</td></tr><tr><td>$U_{33}$</td><td>0,02724</td><td>0,0654</td></tr></tbody></table>

## 3.4 Symétrie de la tension de surface

La littérature [ASB2012] donne les valeurs expérimentales et calculées de la tension de surface pour les familles de plans {100} et {111} (notée $\sigma$ dans ce paragraphe pour reprendre les notations de [ASB2012]). Elles sont proches de l'unité mais peuvent varier dans un rapport compris entre 1 et 2.

<table><caption>TABLE 3.2 – Tension de surface en fonction de l'orientation et de la température</caption><thead><tr><th>T(K)</th><th>$\sigma_{100}$ Jm<sup>-2</sup></th><th>$\sigma_{111}$ Jm<sup>-2</sup></th><th>$\frac{\sigma_{100}}{\sigma_{111}}$</th></tr></thead><tbody><tr><td>2200 (calcul)</td><td>1,602 ± 0,016</td><td>1,137 ± 0,032</td><td>1,408 ± 0,042</td></tr><tr><td>1770 (exp.)</td><td></td><td>0,9 ± 0,2</td><td></td></tr><tr><td>300 (exp.)</td><td></td><td>1,8 ± 0,3</td><td></td></tr></tbody></table>

## 3.5 Le facettage des grains

Les observations au MEB à effet de champ de $UO_2$ cru (avant frittage) montre des ponts préexistants entre particules avec des grains dont la surface est délimitée par des facettes qui
<figure><img src="image_27.png" /><figcaption>FIGURE 3.1 – $UO_2$ observé au MEB FEG (x200000</figcaption></figure>

reflètent la cristallographie du matériau. Cette morphologie particulière peut éventuellement être simplement due au précurseur de l'oxyde. Tant que la porosité ouverte reste importante, on peut penser que la morphologie des grains va avoir une influence sur leur vitesse de croissance. Au terme du frittage (lorsque la porosité occupe une faible fraction volumique et que les surfaces grain/grain sont plus importantes que les surfaces grain/gaz), il n'est pas évident que la forme de cuboctaèdre de Kelvin figure 3.2 observée, pour les grains de la dizaine à quelques dizaines de micromètres dans le matériau dense, sur les coupes céramographiques 3.3, soit liée aux propriétés cristallographiques plutôt qu'au simple nombre de premiers voisins d'un grain pour un matériau dans lequel les grains grossissent à la même vitesse dans toutes les zones du matériau (avec des modulations en fonction des variations de rayons des grains mais sans croissance anormale).

<figure><img src="image_28.png" /><figcaption>FIGURE 3.2 – Cuboctaèdre de Kelvin [MC2003]</figcaption></figure>
<figure><img src="image_29.png" /><figcaption>FIGURE 3.3 – Observation des grains de UO₂ au voisinage d'un macropore - thèse E.Girard</figcaption></figure>

## 3.6 Conséquences sur les équations de conservation - spécificité des cristaux ioniques et de la symétrie cubique

Le solide est ionique et cristallisé, les cations sont localisés sur un sous-réseau cationique et les anions sur un, voire deux, sous-réseaux si l'on prend en compte les positions des ions O²⁻ données par Willis en diffraction de neutrons [BTW1964]. De même, les sous-réseaux de lacunes octaédriques et tétraédriques peuvent recevoir des ions interstitiels.

La composition chimique de la maille d'ordre 4 peut s'écrire sous une forme du type de celle utilisée dans [CG2008] :

$$ \{U'_M, U_M^x, U_M, U_{\ddot{M}}, Pu_M^x, Pu'_M, \Box_M^{4'}\}_4 \{O''_{i_1}, \Box_{i_1}^x\}_{64} \{U_{i_2}^4, Pu_{i_2}^4, \Box_{i_2}^x\}_4 \{O_O^x, \Box_{\ddot{O}}\}_8 \quad (3.2) $$

avec :

$U'_M$ : uranium 3⁺ sur un site cationique,

$U_M^x$ : uranium 4⁺ sur un site cationique,

$U_M$ : uranium 5⁺ sur un site cationique,

$U_{\ddot{M}}$ : uranium 6⁺ sur un site cationique,

$Pu_M^x$ : plutonium 4⁺ sur un site cationique,

$Pu_M$ : plutonium 3⁺ sur un site cationique,

$\Box_M^{4'}$ : lacune cationique,

$O_{i_1}$ : ion oxyde sur un site interstitiel anionique (on confond les deux types de sites interstitiels anioniques décrits par Willis),

$\Box_{i_1}^x$ : site interstitiel anionique vide,

$U_{i_2}^4$ : uranium 4⁺ sur un site interstitiel cationique (la littérature ne fait, cependant, état que de Pu en site interstitiel en conditions réductrices et à haute température),

$Pu_{i_2}^4$ : plutonium 4⁺ sur un site interstitiel cationique,

$\Box_{i_2}^x$ : site interstitiel cationique vide,
$O_O^x$ : ion oxyde sur un site anionique, $\square_O^\ddot{O}$ : site anionique vide, $H_{\ddot{O}}^\ddot{O}$ : hydrogène atomique sur un site d'oxygène. Cette espèce, signalée par Sherman [DFS1989], est prise en compte : on suppose que l'hydrogène diffuse sous forme atomique dans les sites anioniques vacants dans le volume des grains de $UO_2$.

Les équilibres solide-gaz, liés par exemple aux échanges d'oxygène, doivent alors être écrits en terme d'unité de structure. Un jeu d'unités de structure de base doit être choisi pour écrire ces bilans de réactions. Elles doivent être indépendantes et pouvoir représenter l'intégralité de la maille.

Par exemple, on peut choisir d'utiliser la base suivante d'unités de structure :

$$
\begin{align*}
U_1 &\equiv \{U_M + U'_M - 2U_M^x\} \\
U_2 &\equiv \{U_{\dot{M}} + 2U'_M - 3U_M^x\} \\
U_3 &\equiv \{Pu_M^x - U_M^x\} \\
U_4 &\equiv \{-U'_M + Pu'_M\} \\
U_5 &\equiv \{-4U'_M + 3U_M^x + V_M^{4'}\} \\
U_6 &\equiv \{-2U'_M + 2U_M^x + O_{i1}'' - V_{i1}^x\} \tag{3.3} \\
U_7 &\equiv \{4U'_M - 4U_M^x + U_{i2}^4 - V_{i2}^x\} \\
U_8 &\equiv \{4U'_M - 4U_M^x + Pu_{i2}^4 - V_{i2}^x\} \\
U_9 &\equiv \{2U'_M - 2U_M^x - O_O^x + V_{\ddot{O}}^\ddot{O}\} \\
U_{10} &\equiv \{2U'_M - 2U_M^x - O_O^x + H_{\ddot{O}}^\ddot{O}\}
\end{align*}
$$

Cette description permet d'éviter des erreurs d'exposant dans la dépendance de la concen-
tration des défauts ponctuels à la pression partielle d'oxygène.

Par ailleurs, une seconde conséquence du caractère cristallisé du matériau intervient quand on souhaite résoudre les équations de diffusion chimique ou d'interdiffusion. En effet, lorsqu'un joint de grain se déplace, le flux des espèces doit être nul au travers de cette interface dans le référentiel du joint de grain. Mais le matériau étant cristallisé, il faut que ce soit le flux des unités de constructions [JSK1987] qui soit nul. Une autre façon de le dire est que non seulement la conservation de la matière et de l'électroneutralité doivent être respectées mais aussi la conservation du nombre de sites par unité de volume. Cependant, jusqu'à présent, l'interdiffusion a été laissée de côté dans le modèle à l'échelle sub-granulaire.
4

Conclusion

Les spécificités des actinides qui peuvent influer sur leur frittage sont soulignées. Ce sont leur caractère radioactif avec ses conséquences sur les bilans thermochimiques, leurs propriétés chimiques liées à leur nombre de charge élevé avec ses conséquences sur l'écart à la stœchiométrie en oxygène et, de là, sur les valeurs des coefficients d'interdiffusion, mais aussi la symétrie cristalline avec ses conséquences sur l'anisotropie mécanique (tout au moins à température ambiante sur UO$_2$). Une dernière propriété importante pour les équations de conservation mises en jeu pour traduire le bilan des flux de matière en mouvements d'interfaces est la cristallochimie avec la conservation du nombre de sites par unité de volume.
Troisième partie

Modèle développé à l'échelle
sub-granulaire
<figure><img src="image_30.png" /></figure>
# 1

## Cas d'une composition chimique constante

### Sommaire

<table><tr><td>1.1</td><td>Choix du modèle</td><td>59</td></tr><tr><td>1.2</td><td>Hypothèses du modèle</td><td>61</td></tr><tr><td>1.3</td><td>Effets de la courbure de surface</td><td>62</td></tr><tr><td>1.3.1</td><td>Calcul numérique de la courbure</td><td>62</td></tr><tr><td>1.3.2</td><td>Effet mécanique de la tension de surface</td><td>62</td></tr><tr><td>1.3.3</td><td>Conditions aux limites de type Neumann</td><td>63</td></tr><tr><td>1.3.4</td><td>Conditions algèbriques</td><td>63</td></tr><tr><td>1.3.5</td><td>Résolution numérique du problème</td><td>65</td></tr><tr><td>1.3.6</td><td>Comportement élastique</td><td>65</td></tr><tr><td>1.4</td><td>Potentiel chimique généralisé</td><td>69</td></tr><tr><td>1.5</td><td>Évolution des interfaces entre sous-domaines</td><td>70</td></tr><tr><td>1.5.1</td><td>Modélisation des flux</td><td>70</td></tr><tr><td>1.5.2</td><td>Évolution des interfaces grain/gaz</td><td>70</td></tr><tr><td>1.5.3</td><td>Évolution des joints de grains</td><td>71</td></tr><tr><td>1.5.4</td><td>Calcul de la vitesse d'interface</td><td>72</td></tr><tr><td>1.5.5</td><td>Problème de la conservation de la matière</td><td>73</td></tr><tr><td>1.5.6</td><td>Problème de l'évolution de l'angle dièdre</td><td>74</td></tr><tr><td>1.5.7</td><td>Qualité du maillage de surface</td><td>76</td></tr><tr><td>1.6</td><td>Choix du pas de temps</td><td>76</td></tr><tr><td>1.7</td><td>Conclusion</td><td>77</td></tr></table>

### 1.1 Choix du modèle

Sous irradiation, les produits de fission gazeux peuvent migrer lentement en volume du cœur des grains vers les joints de grains, puis, plus rapidement, vers l'espace entre la pastille et la gaine. Dans ces conditions, le gaz joue alors le rôle d'une barrière thermique et freine le refroidissement du crayon combustible ou de l'aiguille par le caloporteur. La taille des grains du combustible à l'issue de la fabrication est ainsi un paramètre important de la microstructure du combustible.
Par ailleurs, ces produits de fission peuvent être piégés dans des porosités intragranulaires et leur répartition entre les pores et la matrice cristalline résulte d'un équilibre chimique. Cette porosité provient en partie de la porosité résiduelle à l'issue de la fabrication (~ 5 % en volume) et la maîtrise de sa distribution de taille est nécessaire.

Aussi a-t-on choisi de développer un modèle à l'échelle sub-granulaire du frittage. Par ailleurs, le moteur du frittage, qui est la variation du rayon de courbure à la surface des grains (surface libre et joints de grains) ne peut être appréhendée qu'à l'échelle sub-granulaire. Un des objectifs initiaux a aussi été d'utiliser des paramètres d'entrée pour le modèle qui soient des grandeurs mesurables par ailleurs sur la base d'expériences dédiées ne mettant pas elles-mêmes en jeu le frittage. C'est aussi une façon de s'assurer que l'on prend en compte les spécificités liées aux matériaux étudiés.

Les modèles classiques de frittage reposent sur l'hypothèse d'un seul phénomène (diffusion en volume, en surface, aux joints de grains, évaporation-condensation) à un instant donné. Tou- tefois, pendant le second stade du frittage, la densification a lieu alors que les grains grossissent et que la forme des pores s'arrondit. Pour certains matériaux les premier et second stades sont confondus. Aussi, le modèle proposé ici ne repose pas sur l'hypothèse d'un unique phénomène prépondérant à un instant donné.

À cette échelle, plusieurs types d'approches sont possibles, notamment des approches de type Monte-Carlo [VT2003], [MB2005], [SB2002], [SB2005] dans lesquelles le matériau est découpé en pixels en 2D, ou en voxels en 3D. Dans ces modèles, une énergie est attribuée à chacune des configuration de voxels servant à représenter la microstructure. Les voxels peuvent être déplacés selon certaines règles pour modifier l'interface grain/gaz. Une configuration est gardée lorsqu'elle conduit à une diminution de l'énergie du système. La courbure, par exemple, est prise en compte en comptant le nombre de voxels de type gaz contenus dans un cube de $b^3$ voxels où b est la taille de la zone servant à lisser la courbure. La relation entre les itérations de Monte-Carlo et le temps peut-être déterminée lorsqu'un seul phénomène diffusionnel intervient, elle est moins évidente lorsque ce n'est pas le cas.

Bruchon [JB2010], [JB2011], [JB2012], [JB2013] a développé un modèle pour un système mono-constituant de grains qui frittent, mais en laissant de côté la notion de joints de grains (i.e. les joints de grains existent à l'instant initial, mais ne sont plus pris en compte dans les instants ultérieurs).

La littérature contient d'autres approches, par exemple de type champ de phase, dans lesquelles on définit, pour chacun des grains une fonction qui vaut 1 à l'intérieur du grain et 0 à l'extérieur et qui passe continûment de 1 à 0 sur l'épaisseur $\delta_{GB}$ du joint de grain. Leur évolution reposent sur les équations de Ginzburg-Landau dépendant du temps [LQC1994]. Cette approche, qui utilise les dérivées partielles par rapport aux paramètres d'ordre de l'énergie libre totale, a été étendue récemment par Nestler [BN2005] au cas de la diffusion multicomposants avec des énergies de surface et des coefficients de diffusion anisotropes en prenant en compte l'apparition éventuelle de nouvelles phases.

Wakai a récemment proposé un modèle à l'échelle sub-granulaire [FW2011] qui considère à la fois la diffusion aux joints de grains et en surface pour simuler l'évolution de la forme des grains lors du frittage. Cette évolution repose sur le "surface evolver" de Brakke [KB1992]. Wakai a
proposé une approche théorique de l'évolution d'un agrégat de grains en tant que solution d'une équation tensorielle du viriel [FW2013].

Le modèle que nous avons développé [JL1999], [JL2001], [MA2006] est présenté dans sa formule forte dans ce qui suit. Les algorithmes et les choix numériques reposent sur l'aide de Robert Boyer et de Kacem Saikouk du Centre de Mathématiques et d'Informatique de Marseille dont j'ai bénéficié au cours de ces années. Dans ce cadre, un contrat de collaboration scientifique avait été établi entre le Centre de Mathématiques et d'Informatique de Marseille et le CEA. Ces algorithmes ont été mis en place en vue d'une parallélisation ultérieure possible, et l'algorithme de décomposition de domaines de type Dirichlet-Neumann proposé a été programmé, utilisé dans une thèse encadrée (M. Ajdour [MA2006]), mais actuellement laissé provisoirement de côté au profit d'une résolution globale qui servira de référence aux autres méthodes de résolution.

1.2 Hypothèses du modèle

La simulation ne prend pas en compte le réarrangement des cristallites qui correspondrait à
un glissement ou à un roulement des grains aux niveau des cous qui les connectent à leurs voisins.
Ceci signifie qu'une fois le contact entre deux grains établi, celui-ci demeure tout au long de la
simulation, ou d'un frittage comportent un grand nombre de grains avec de fortes inhomogénéités
initiales de répartition des grains. Ce phénomène serait essentiellement à prendre en compte dans
le cas d'un frittage sous charge.

Les grains sont supposés être monocristallins même s'ils contiennent des défauts ponctuels
éventuellement organisés en dislocations coins ou vis, en boucles de dislocations ou en joints de
flexion. Lors de la montée en température lors du cycle de frittage ces défauts sont partiellement
recuits sous l'effet de la diffusion en volume. Comme ce dernier phénomène est le plus lent des
phénomènes de diffusion en phase solide, il est probable qu'une partie de ces défauts subsiste,
cependant, lors du frittage. Les défauts ponctuels peuvent provenir par exemple d'un écart à la
stœchiométrie dû à la pression partielle d'oxygène, leur association en défauts plus complexes
à un dopage (par exemple $Cr^{3+}$ dans la structure fluorine $UO_2$). Les dislocations peuvent être
induites par les étapes de broyage et de pressage (pour cette dernière, typiquement, la pression
appliquée est de l'ordre de 450MPa pour le cas du combustible nucléaire). Ainsi, aux grains
monocristallins peuvent être associés 3 angles $(\phi, \theta, \psi)$ pour relier l'orientation de leur réseau
cristallin $(\vec{a}, \vec{b}, \vec{c})$ aux axes du laboratoire $(\vec{e}_1, \vec{e}_2, \vec{e}_3)$. On suppose que le frittage est libre, i.e.
que la contrainte appliquée est suffisamment faible pour que les déformations des grains restent
élastiques, même si la diffusion conduit à une déformation irréversible de l'ensemble. Le gradient
d'énergie mécanique emmagasinée localement dans les grains est à l'origine de flux irréversibles de
matière. Bien que les céramiques telles que $UO_2$ ou $(U_{1-y}Pu_y)O_{2-x}$ soient denses $(\sim 11g.cm^{-3})$
le poids n'a pas été pris en compte dans la simulation.

Pour illustrer le modèle et pour les tests actuels, le cas de deux grains au contact est considéré.
La surface libre du grain 1 est notée $\Gamma_1$, celle du grain 2 $\Gamma_2$, le joint de grain $\Gamma_{1/2}$ et la ligne
triple $\Delta_{1/2}$. Le volume des grains est noté $\Omega_i$ ($i = 1$ ou 2).

Chaque grain est considéré comme étant un milieu continu mais caractérisé par trois angles
qui font le lien avec l'échelle atomique.
## 1.3 Effets de la courbure de surface

La courbure aux frontières d'un grain, que ce soient des surfaces libres ou des joints de grains, est à l'origine des variations de contraintes le long des contours des grains. La courbure agit sur la surface en exerçant une force $d\vec{f}$ de type Laplace sur un élément de surface $dS$. Sa résultante est normale à la surface.

$$d\vec{f} = \gamma \left( \frac{1}{R_1} + \frac{1}{R_2} \right) \vec{n} dS \qquad (1.1)$$

où $R_1$ et $R_2$ sont les rayons de courbure principaux et $\vec{n}$ la normale à la surface.

Dans le matériau constituant le combustible nucléaire cru, à l'issue du pressage, les agrégats sont constitués de 1 à 3 cristallites et la largeur des raies de diffraction montre que leur taille est de 0,1 à 0,3 micromètres. Au terme du frittage la taille des grains est de l'ordre de 10 micromètres. La tension de surface est de l'ordre de $1J \cdot m^{-1}$, ainsi, dans le matériau cru la contrainte en surface des grains est de l'ordre de $10MPa$. On voit que le moteur de la diffusion va être très sensible à la valeur de la courbure moyenne locale $H$ définie par

$$H = \frac{1}{2} \left( \frac{1}{R_1} + \frac{1}{R_2} \right) \qquad (1.2)$$

### 1.3.1 Calcul numérique de la courbure

Dans le cas d'une surface discrétisée l'évaluation de la courbure locale moyenne $H$ est critique pour la bonne simulation de l'évolution de la surface. Diverses techniques ont été mises en œuvre et testées [MA2006] puis [MW2007]. Ce problème est très délicat et a fait l'objet de nombreux travaux relatés dans la littérature [TS2003] [DCS2003]; c'est la pierre d'angle de l'évolution de la surface. Au final, nous nous sommes tournés vers une bibliothèque numérique spécialisée, la bibliothèque Computational Geometry Algorithms Library (CGAL [AF1996]) en utilisant plus particulièrement l'application *Jet\_fitting\_3D* [FC2008]. Elle ajuste, au sens des moindres carrés, pour chacun des nœuds pour lequel on souhaite calculer $H$ (ceux appartenant à la frontière des grains), à partir des coordonnées d'un ensemble de plus proches voisins, une forme de Monge. Cette dernière s'écrit dans un repère où l'axe des $z$ est dirigé par la normale sortante du grain et dont l'origine est le nœud considéré.

$$z(x, y) = \frac{1}{2}(k_1 x^2 + k_2 y^2) + \frac{1}{6}(b_0 x^3 + 3b_1 x^2 y + 3b_2 xy^2 + b_3 y^3) + \sqrt{x^2 + y^2} \epsilon(x, y) \qquad (1.3)$$

Le nombre de couches de premiers voisins à prendre en compte pour un nœud donné est un paramètre critique. Si le nombre de couches de voisins est trop petit les informations de courbure liées au calcul à l'issue d'un pas de temps seront bruitées, s'il est trop grand ces informations de courbure seront trop moyennées. Dans les deux cas les flux de matière liés à ces variations de courbure ne seront pas corrects.

### 1.3.2 Effet mécanique de la tension de surface

L'équilibre quasi-statique des forces s'écrit dans l'hypothèse des petites déformations à l'aide de l'équation de Navier-Lamé :

$$\sum_{i=1}^{3} \frac{\partial \sigma_{ij}}{\partial x_j} + \rho f_i = 0, i \in \{1, ..., 3\} \qquad (1.4)$$
où $\sigma_{ij}$ désigne la composante $ij$ du tenseur des contraintes et $\rho$ la masse volumique du solide. Cette approche, de type petites déformations, est adaptée car il s'agit seulement de rendre compte de la densité d'énergie élastique emmagasinée dans le solide à chaque instant, la déformation globale, qui elle est importante ($\frac{\Delta L}{L_0} \simeq -15\%$) entre le début et la fin du frittage est due non pas à l'élasticité mais à la déformation liée aux flux irréversibles de matière.

### 1.3.3 Conditions aux limites de type Neumann

Une partie des conditions de fermeture du problème sont des conditions aux limites de type von Neumann.

Pour une interface grain/gaz elles traduisent une force de Laplace exercée localement à la surface du fait des contraintes liées à la courbure ainsi que la pression du gaz-procédé :

$$ \overline{\overline{\sigma}}(\vec{n}) \cdot \vec{n} = -p_{gaz} - 2\gamma_{SV}H_{SV} \quad (1.5) $$

L'indice $SV$ désignant l'interface grain/gaz (solid/vapor en anglais).

Pour une interface grain/grain elles traduisent une force de Laplace exercée localement à la surface du fait des contraintes liées à la courbure ainsi que la contrainte exercée par le grain voisin au contact :

$$ \overline{\overline{\sigma_1}}(\vec{n}_1) \cdot \vec{n}_1 = \overline{\overline{\sigma_2}}(\vec{n}_1) \cdot \vec{n}_1 - 2\gamma_{GB}H_{GB} \quad (1.6) $$

L'indice $GB$ désignant l'interface grain/grain (grain boundary en anglais), et l'indice de $\sigma$ désignant le domaine (au sens mathématique du terme, ce qui correspond au grain) pour lequel $\sigma$ a été calculé.

### 1.3.4 Conditions algébriques

Les autres conditions de fermeture sont de type algébrique et permettent d'éliminer les mouvements de corps rigide de l'ensemble des solutions du problème.

Les forces extérieures au système, constitué d'un ensemble de grains (ici deux) au contact en cours de frittage, sont le poids, la réaction du support et les forces de pression en surface, les deux premières ne sont pas prises en compte ici.

#### Effet des forces de pression

On suppose que la zone de contact entre l'assemblée de grains et le support est de mesure nulle de sorte que le gaz soit au contact sur toute la surface de l'ensemble des deux solides. Cela peut être vrai au début du frittage, par contre, lorsque celui-ci est suffisamment avancé, la zone de contact avec le support est probablement étendue.

On montre, sur l'exemple 2D, que la somme des forces de pression appliquées à la surface est de résultante nulle.

On considère un solide (une assemblée de grains au contact) de contour continu fermé (les joints de grains sont, pour cette partie du raisonnement, des frontières internes au solide) avec une origine O sur le contour et une abscisse curviligne s. La force de pression exercée sur un élément de contour de longueur $ds$ s'écrit :

$$ d\vec{f} = - \left( \gamma \frac{1}{\rho(s)} + p_0 \right) ds \overrightarrow{n(s)} = -\gamma \overrightarrow{t(s)} - p_0 \overrightarrow{n(s)} \quad (1.7) $$
soit, en sommant sur tout le contour :

$$ \vec{f} = -\gamma \int_{0}^{L} \overrightarrow{t(s)} ds - p_{0} \int_{0}^{L} \overrightarrow{n(s)} ds \quad (1.8) $$

Le contour étant fermé le premier terme est nul, tandis que pour le second :

$$ -p_0 \int_0^L \overrightarrow{n(s)} ds = \iint_S \overrightarrow{\nabla}(p) dS \quad (1.9) $$

Si le solide est incompressible, la pression est constante $p = p_0$ dans le solide et le dernier terme s'annule.

La somme des forces appliquées à l'ensemble des grains en contact est nulle. On démontre aussi que le moment, par rapport à un point, des forces de pression et de tension de surface s'annule si ce point est le centre de gravité du solide.

## Conséquences sur le champ de déplacements

Si l'on considère une transformation élastique du solide dans laquelle le champ de déplacements passe de $\vec{0}$ à $\vec{u}$, la somme des forces étant nulle à chaque instant, le centre de gravité garde sa vitesse initiale qui est nulle.

Par ailleurs, le moment des déplacements par rapport au centre de gravité reste lui aussi nul.

## Propriétés du champ de déplacements

L'absence de translation du centre de gravité des deux grains donne en notant :
$G$ le centre de gravité après déformation devient $G'$

$M$ un point du solide qui devient après déformation $M'$ tel que $\overrightarrow{MM'} = \vec{u}$

$$ \left\{ \begin{array}{l} \iiint_{\Omega_1 \cup \Omega_2} \rho \overrightarrow{GM} d\tau = \vec{0} \quad (\text{definition de } G) \\ \iiint_{\Omega_1 \cup \Omega_2} \rho \overrightarrow{G'M'} d\tau = \vec{0} \quad (\text{definition de } G') \\ \overrightarrow{GG'} = \vec{0} \quad (\text{centre de gravité immobile}) \end{array} \right. \quad (1.10) $$

où $\rho$ est la masse volumique du matériau que l'on peut supposer pratiquement constante (solide quasi-incompressible). Ainsi :

$$ \iiint_{\Omega_1 \cup \Omega_2} \vec{u} \, d\tau = \vec{0} \quad (1.11) $$

L'absence de rotation des deux grains autour du centre de gravité donne :

$$ \iiint_{\Omega_1 \cup \Omega_2} \overrightarrow{GM} \wedge \vec{u} \, d\tau = \vec{0} \quad (1.12) $$
Autre justification des conditions algébriques de fermeture

Lorsque l'on s'intéresse à une assemblée de grains au contact qui frittent, on peut imaginer
que pour observer l'état de ce système à un instant donné on le tourne et/ou on l'approche ou
on l'éloigne de l'observateur. Il n'empêche que les distances en angles relatifs entre deux points
du systèmes resteront inchangés. Pour traduire cette invariance on peut définir une relation
d'équivalence R entre deux solides S et S' disant qu'ils sont identiques s'il existe un torseur de
vecteur de translation I et de vecteur de rotation ω tels que :

$$
\overrightarrow{G'M'} = \vec{I} + \overrightarrow{GM} \wedge \vec{\omega} \qquad (1.13)
$$

$\vec{I}$ et $\vec{\omega}$ sont alors donnés par :

$$
\left\{
\begin{array}{l}
\vec{I} = \\
\vec{\omega} \in \left\{ -\frac{\overrightarrow{GM} \wedge (\overrightarrow{MM'} - \vec{I})}{GM^2} + k \overrightarrow{GM}, k \in \mathbb{R} \right\}
\end{array}
\right.
\qquad (1.14)
$$

En fait, on s'intéresse seulement à l'ensemble quotient de solutions modulo la relation d'équivalance $\mathcal{R}$ disant qu'il existe un torseur $(\vec{I}, \vec{\omega})$ amenant le solide solution du problème d'élasticité sur un autre solide lui aussi solution.

1.3.5 Résolution numérique du problème

Le champ de déplacements $\vec{u}$ est calculé par la méthode des éléments finis, à l'aide d'éléments de type $P_2$, afin que les dérivées secondes de $\vec{u}$, qui sont à l'origine des flux, soient non nulles.

La formulation faible pour un grain isolé avec des conditions aux limites de type Neumann s'obtient en intégrant par parties 1.4 :

$$
\forall w \in \{\text{champs admissibles}\} \iint\int_{\Omega_1 \cup \Omega_2} (\mathrm{div}(\sigma) + \rho \vec{f}) \cdot \vec{w} \, d\tau = 0 \quad (1.15)
$$

soit :

$$
\forall w \in \{\text{champs admissibles}\} \iint\int_{\Omega_1 \cup \Omega_2} (\bar{\sigma} : \bar{e}[\bar{w}]) \, dV = \iint_{\partial(\Omega_1 \cup \Omega_2)} (\bar{\sigma} \cdot \bar{n}) \cdot [\bar{w}] \, dV \quad (1.16)
$$

1.3.6 Comportement élastique

Le comportement élastique du matériau cubique (groupe d'espace $Fm\bar{3}m$) se traduit par un
tenseur élastique à 3 coefficients indépendants $c_{11}$, $c_{12}$ et $c_{44}$ dans le repère des axes du cristal,
i.e. d'un grain donné (en notations de Voïgt) :

$$
\overline{\overline{c}} = \begin{pmatrix} c_{11} & c_{12} & c_{12} & 0 & 0 & 0 \\ c_{12} & c_{11} & c_{12} & 0 & 0 & 0 \\ c_{12} & c_{12} & c_{11} & 0 & 0 & 0 \\ 0 & 0 & 0 & c_{44} & 0 & 0 \\ 0 & 0 & 0 & 0 & c_{44} & 0 \\ 0 & 0 & 0 & 0 & 0 & c_{44} \end{pmatrix} \tag{1.17}
$$
Pour $UO_2$ [SRP2009], à température ambiante, le cas cubique est très éloigné du cas isotrope, pour lequel on aurait $c_{44} = \frac{1}{2} (c_{11} - c_{12})$.

En effet,

$$ \left\{ \begin{array}{lcl} c_{11} & = & 389 - 396 \, \text{GPa} \\ c_{12} & = & 119 - 121 \, \text{GPa} \\ c_{44} & = & 59,7 - 64,1 \, \text{GPa} \end{array} \right. \qquad (1.18) $$

Le cas cubique a été implémenté dans le logiciel que l'on a développé ("Salammbo" - écrit en C++), bien que la constante $c_{44}$ ait été ajustée à partir des valeurs de $c_{11}$ et de $c_{12}$ pour correspondre au cas isotrope, ceci à des fins de test.

La loi de comportement élastique s'écrit :

$$ \bar{\bar{\sigma}} = \bar{\bar{c}} : \bar{\bar{\epsilon}} \qquad (1.19) $$

$$ \text{avec : } \epsilon_{ij} = \frac{1}{2} \left( \frac{\partial u_i}{\partial x_j} + \frac{\partial u_j}{\partial x_i} \right) $$

Chaque grain $g$ a une orientation $(\phi^g, \theta^g, \psi^g)$ qui lui est propre, ce qui va impacter, dans la cas d'une symétrie cubique non isotrope, les équations de l'équilibre mécanique (quasi-statique). Ces dernières s'écrivent, en fonction des déplacements seuls, dans le référentiel des axes du cristal :

$$ \left\{ \begin{array}{l} c_{11} \frac{\partial^2 u_1^g}{\partial^2 x_1^g} + c_{12} \frac{\partial^2 u_2^g}{\partial x_1^g \partial x_2^g} + c_{12} \frac{\partial^2 u_3^g}{\partial x_1^g \partial x_3^g} + \frac{1}{2} c_{44} \left( \frac{\partial^2 u_1^g}{\partial x_1^g \partial x_2^g} + \frac{\partial^2 u_2^g}{\partial x_2^g \partial x_3^g} \right) + \frac{1}{2} c_{44} \left( \frac{\partial^2 u_3^g}{\partial x_3^g \partial x_1^g} + \frac{\partial^2 u_3^g}{\partial x_3^g \partial x_2^g} \right) = 0 \\ \frac{1}{2} c_{44} \left( \frac{\partial^2 u_2^g}{\partial x_2^g \partial x_1^g} + \frac{\partial^2 u_3^g}{\partial x_3^g \partial x_1^g} \right) + c_{12} \frac{\partial^2 u_1^g}{\partial x_1^g \partial x_3^g} + c_{11} \frac{\partial^2 u_2^g}{\partial x_2^g \partial x_3^g} + \frac{1}{2} c_{44} \left( \frac{\partial^2 u_2^g}{\partial x_2^g \partial x_3^g} + \frac{\partial^2 u_3^g}{\partial x_3^g \partial x_1^g} \right) = 0 \\ \frac{1}{2} c_{44} \left( \frac{\partial^2 u_3^g}{\partial x_3^g \partial x_1^g} + \frac{\partial^2 u_3^g}{\partial x_3^g \partial x_2^g} \right) + \frac{1}{2} c_{44} \left( \frac{\partial^2 u_3^g}{\partial x_3^g \partial x_2^g} + \frac{\partial^2 u_3^g}{\partial x_3^g \partial x_3^g} \right) + c_{12} \frac{\partial^2 u_1^g}{\partial x_1^g \partial x_3^g} + c_{12} \frac{\partial^2 u_2^g}{\partial x_2^g \partial x_1^g} + c_{11} \frac{\partial^2 u_3^g}{\partial x_3^g \partial x_2^g} + c_{11} \frac{\partial^2 u_3^g}{\partial x_3^g \partial x_3^g} = 0 \end{array} \right. \qquad (1.20) $$

L'expression du système d'équations utile pour le code développé est celle dans le repère du laboratoire.

<figure><img src="image_31.png" /><figcaption>FIGURE 1.1 – Axes du laboratoire et du cristal</figcaption></figure>

On note $c'_{\alpha\beta\gamma\delta}$ les composantes du tenseur élastique du grain $g$ dans le repère du laboratoire et
$^gP$ la matrice de passage du repère du laboratoire à celui du grain $g$ (pour les changements de référentiel les notations de Voïgt ont été évitées).

$$
c'_{\alpha\beta\gamma\delta} = c_{ijkl} {}^gP_{\alpha}{}^i {}^gP_{\beta}{}^j {}^gP_{\gamma}{}^k {}^gP_{\delta}{}^l \quad (1.21)
$$

où les *i*, *j*, *k*, *l* sont les indices des lignes de la matrice de passage. On précise la matrice de passage
[LA2002] :

$$
\begin{pmatrix}
\cos(\psi^g)\cos(\phi^g) - \sin(\phi^g)\sin(\psi^g)\cos(\theta^g) & -\sin(\psi^g)\cos(\phi^g) - \sin(\psi^g)\cos(\phi^g)\cos(\theta^g) & \sin(\theta^g)\sin(\psi^g) \\
\cos(\phi^g)\sin(\psi^g) + \cos(\psi^g)\sin(\phi^g)\cos(\theta^g) & -\sin(\phi^g)\sin(\psi^g) + \cos(\phi^g)\cos(\theta^g)\cos(\psi^g) & -\sin(\theta^g)\cos(\psi^g) \\
\sin(\phi^g)\sin(\theta^g) & \cos(\psi^g)\sin(\theta^g) & \cos(\theta^g)
\end{pmatrix}
\tag{1.22}
$$

L'angle $\psi^g = (\vec{x}_0, \vec{u})$, compté autour de $\vec{z}$, est la précession, l'angle $\theta^g = (\vec{z}, \vec{c})$, compté autour de $\vec{u}$, la mutation, l'angle $\phi^g = (\vec{u}, \vec{a})$, compté autour de $\vec{c}$, la rotation propre.
Le système d'équations différentielles s'écrit, lors de ce changement de repère :

$$
\forall M \in \Omega^g, \quad \left\{
\begin{array}{@{}l@{\quad}l@{}}
\displaystyle \sum_{i=1}^{3} \sum_{j=1}^{3} \sum_{k=1}^{3} C_{ijk1} \frac{\partial^2 u_i}{\partial x_j \partial x_k} \Big|_M \, d\tau & = 0 \\[2ex]
\displaystyle \sum_{i=1}^{3} \sum_{j=1}^{3} \sum_{k=1}^{3} C_{ijk2} \frac{\partial^2 u_i}{\partial x_j \partial x_k} \Big|_M \, d\tau & = 0 \\[2ex]
\displaystyle \sum_{i=1}^{3} \sum_{j=1}^{3} \sum_{k=1}^{3} C_{ijk3} \frac{\partial^2 u_i}{\partial x_j \partial x_k} \Big|_M \, d\tau & = 0
\end{array}
\right.
\qquad (1.23)
$$

avec $C_{ijkl} = \frac{1}{2} c'_{ljik} + \frac{1}{2} c'_{ljki}$

Formulation faible du problème mécanique

La composante $i \in \{1, ..., 3\}$ du champ $u$ est développée sur une base de polynômes du second degré $N^c$ avec comme coordonnées les valeurs $U_i^c$ de la fonction aux nœuds du maillage :

$$
u_i = \sum_{c=1}^{10} U_i^c N^c (x_1, x_2, x_3) \quad (1.24)
$$

L'intégration sur un tétraèdre K d'un maillage J du sous-domaine Ω<sup>g</sup>, (g ∈ {1, 2} représentant
le numéro du grain), conduit au système :

$$
\forall K \in \mathcal{T}_{\Omega^g}, \forall M \in K \quad \left\{
\begin{array}{@{}l@{\quad}l@{}}
\displaystyle \sum_{i=1}^{3} \sum_{j=1}^{3} \sum_{k=1}^{3} \iiint_{\mathcal{T}} C_{ijk1} N^b \sum_c U_i^c \frac{\partial^2 N^c}{\partial x_j \partial x_k} \Big|_M \, d\tau & = 0 \\[2ex]
\displaystyle \sum_{i=1}^{3} \sum_{j=1}^{3} \sum_{k=1}^{3} \iiint_{\mathcal{T}} C_{ijk2} N^b \sum_c U_i^c \frac{\partial^2 N^c}{\partial x_j \partial x_k} \Big|_M \, d\tau & = 0 \\[2ex]
\displaystyle \sum_{i=1}^{3} \sum_{j=1}^{3} \sum_{k=1}^{3} \iiint_{\mathcal{T}} C_{ijk3} N^b \sum_c U_i^c \frac{\partial^2 N^c}{\partial x_j \partial x_k} \Big|_M \, d\tau & = 0
\end{array}
\right.
\quad (1.25)
$$

en intégrant par parties :

$$
\forall b \in \{1, \dots, 10\} \quad \left\{
\begin{aligned}
& \sum_{i=1}^{3} \sum_{j=1}^{3} \sum_{k=1}^{3} \sum_{c} C_{ijk1} N^b U_i^c \begin{bmatrix} \iiint_K \frac{\partial N^b}{\partial x_j} \frac{N^c}{\partial x_k} - \iint_K N^b \frac{\partial N^c}{\partial n} \\ \iiint_K \frac{\partial N^b}{\partial x_j} \frac{N^c}{\partial x_k} - \iint_K N^b \frac{\partial N^c}{\partial n} \\ \iiint_K \frac{\partial N^b}{\partial x_j} \frac{N^c}{\partial x_k} - \iint_K N^b \frac{\partial N^c}{\partial n} \end{bmatrix} \, d\tau = 0 \\
& \sum_{i=1}^{3} \sum_{j=1}^{3} \sum_{k=1}^{3} \sum_{c} C_{ijk2} N^b U_i^c \begin{bmatrix} \iiint_K \frac{\partial N^b}{\partial x_j} \frac{N^c}{\partial x_k} - \iint_K N^b \frac{\partial N^c}{\partial n} \\ \iiint_K \frac{\partial N^b}{\partial x_j} \frac{N^c}{\partial x_k} - \iint_K N^b \frac{\partial N^c}{\partial n} \\ \iiint_K \frac{\partial N^b}{\partial x_j} \frac{N^c}{\partial x_k} - \iint_K N^b \frac{\partial N^c}{\partial n} \end{bmatrix} \, d\tau = 0 \\
& \sum_{i=1}^{3} \sum_{j=1}^{3} \sum_{k=1}^{3} \sum_{c} C_{ijk3} N^b U_i^c \begin{bmatrix} \iiint_K \frac{\partial N^b}{\partial x_j} \frac{N^c}{\partial x_k} - \iint_K N^b \frac{\partial N^c}{\partial n} \\ \iiint_K \frac{\partial N^b}{\partial x_j} \frac{N^c}{\partial x_k} - \iint_K N^b \frac{\partial N^c}{\partial n} \\ \iiint_K \frac{\partial N^b}{\partial x_j} \frac{N^c}{\partial x_k} - \iint_K N^b \frac{\partial N^c}{\partial n} \end{bmatrix} \, d\tau = 0
\end{aligned}
\right. \quad (1.26)
$$
Les éléments de matrice élémentaire $^eA_{rs}^{ij}$ sont calculés dans le logiciel que l'on a développé :

$$ ^eA_{ij}^{rs} = \frac{1}{2} \left\{ \sum_{s=1}^{d} \sum_{l,m=1}^{d} \sum_{j=1}^{n} c_{rslm} u_j^l N_{j,m} N_{i,s} + \sum_{s=1}^{d} \sum_{l,m=1}^{d} \sum_{k=1}^{n} c_{rslm} u_k^m N_{k,l} N_{i,s} \right\} \quad (1.27) $$

ainsi que les termes du membre de droite associés :

$$ ^eF_i = \sum_{s=1}^{d} \iint_{\Gamma} \sum_{l,m=1}^{d} \frac{1}{2} c_{rslm} \left( \sum_{j=1}^{n} u_j^l N_{j,m} + \sum_{k=1}^{n} u_k^m N_{k,l} \right) N_{i,s} d\Omega \quad (1.28) $$

Après assemblage le système est du type :

$$ Au = F \quad (1.29) $$

ou, si l'on sépare les termes ne mettant en jeu que les valeurs de $\vec{u}$ sur les noeuds de $\Omega_i$ qui n'appartiennent pas au joint de grain $\Gamma$ :

$$ \begin{pmatrix} A_{11} & A_{1\Gamma} \\ A_{1\Gamma}^t & A_{\Gamma\Gamma} \end{pmatrix} \begin{pmatrix} U_1 \\ U_\Gamma \end{pmatrix} = \begin{pmatrix} F_1 \\ F_\Gamma \end{pmatrix} \quad (1.30) $$

On calcule à tour de rôle et séparément, pour chacun des deux grains, les éléments de matrice A et F dans notre application, en utilisant les fonctions de la bibliothèque d'éléments finis diffpack, en s'arrêtant, à chaque fois, au stade de l'assemblage. Les conditions de non-translation 1.11 et de non-rotation 1.12 s'expriment aussi sous forme matricielle avec $\vec{U}$ le vecteur des inconnues aux noeuds du maillage :

$$ (B) (U) = (0) \quad (1.31) $$

Cette dernière relation correspond à 6 équations scalaires qui sont des contraintes imposées au champ $\vec{u}$.

On a choisi ici de les imposer au sens des moindres carrés, en se donnant un vecteur d'inconnues supplémentaires (des paramètres de Lagrange) $L \in \mathbb{R}^6$, pour obtenir le système d'équations :

$$ \begin{pmatrix} A_{11} & A_{1\Gamma} & 0 & B_{11}^t \\ A_{1\Gamma}^t & A_{\Gamma\Gamma} & A_{2\Gamma}^t & B_{1\Gamma}^t + B_{2\Gamma}^t \\ 0 & A_{2\Gamma} & A_{22} & B_{22}^t \\ B_{11} & B_{1\Gamma} + B_{2\Gamma} & B_{22} & 0 \end{pmatrix} \begin{pmatrix} U_1 \\ U_\Gamma \\ U_2 \\ L \end{pmatrix} = \begin{pmatrix} F_1 \\ F_\Gamma \\ F_2 \\ 0 \end{pmatrix} \quad (1.32) $$

Afin d'améliorer le conditionnement du système, une méthode de Lagrangien augmenté pourrait être mise en œuvre (méthode de Fortin [MF1982]) où $r$ est un facteur de pondération. Le système devient alors :

$$ \left\{ \begin{pmatrix} A_{11} & A_{1\Gamma} & 0 \\ A_{1\Gamma}^t & A_{\Gamma\Gamma} & A_{2\Gamma}^t \\ 0 & A_{2\Gamma} & A_{22} \end{pmatrix} + r \begin{pmatrix} B_{11}^t B_{11} & B_{11}^t (B_{1\Gamma} + B_{2\Gamma}) & B_{11}^t B_{22} \\ (B_{1\Gamma}^t + B_{2\Gamma}^t) B_{11} & (B_{1\Gamma}^t + B_{2\Gamma}^t) (B_{1\Gamma} + B_{2\Gamma}) & (B_{1\Gamma}^t + B_{2\Gamma}^t) B_{22} \\ B_{22}^t B_{11} & B_{22}^t (B_{1\Gamma} + B_{2\Gamma}) & B_{22}^t B_{22} \end{pmatrix} \right\} \begin{pmatrix} U_1 \\ U_\Gamma \\ U_2 \end{pmatrix} + \begin{pmatrix} B_{11}^t \\ B_{1\Gamma} + B_{2\Gamma} \\ B_{22}^t \end{pmatrix} (L)^n = \begin{pmatrix} F_1 \\ F_\Gamma \\ F_2 \end{pmatrix} \quad (1.33) $$
$L^0$ est obtenu en résolvant le système sans Lagrangien augmenté de solution $U^0$. Ensuite $L^{n+1}$ est donné par :

$$L^{n+1} = L^n + \rho_n B U^n \text{ with } 0 < \rho_n < 2r \quad (1.34)$$

Le Lagrangien non augmenté a été utilisé jusqu'à présent.

La traction exercée sur un élément $\sigma$ ($\sigma = K \cap (\partial\Omega_1 \cup \Omega_2)$) de surface s'écrit :

$$\vec{F}_{\sigma} = \iint_{\sigma} \left( -p - \gamma_{\Sigma} \left( \frac{1}{R_1} + \frac{1}{R_2} \right) \right) \vec{n}_{\Sigma} dS \quad (1.35)$$

## 1.4 Potentiel chimique généralisé

La partie élastique du potentiel chimique généralisé peut être exprimée au moyen du tenseur des contraintes [MA2006], c'est ce qui est présenté dans ce paragraphe. Le travail des forces élastiques par unité de volume a pour expression :

$$dW = \sigma_{ij} du_{ik} \quad (1.36)$$

La quantité de chaleur reçue par ce même élément de volume :

$$\delta Q = T dS \quad (1.37)$$

L'énergie interne par unité de volume s'écrit :

$$dU = \delta Q + dW = T dS + \sigma_{ij} du_{ik} \quad (1.38)$$

et comme $F=U-TS$ il vient :

$$dF = dU - TdS - SdT = -SdT + \sigma_{ij} du_{ik} \quad (1.39)$$

De

$$G = F - \sigma_{ij} u_{ik} \quad (1.40)$$

on déduit :

$$dG = -\sigma_{ik} du_{ik} - u_{ik} d\sigma_{ik} - SdT + \sigma_{ik} du_{ik} = -u_{ik} d\sigma_{ik} - SdT \quad (1.41)$$

Comme on travaille à température constante à chaque pas de temps (équilibre quasi-isotherme) :

$$dG = -u_{ik} d\sigma_{ik} \quad (1.42)$$

$$\left. \frac{\partial G}{\partial n} \right|_T = \mu_{\text{elastique}} = -u_{ik} \frac{\partial \sigma_{ik}}{\partial n} \quad (1.43)$$

où n représente le nombre de moles de l'espèce chimique dont on calcule le potentiel chimique. Le potentiel chimique généralisé d'origine élastique s'écrit alors :

$$\mu_{\text{elastique}} = -\frac{1}{2} \Theta(\epsilon : \sigma) \quad (1.44)$$

soit, en terme de contraintes seulement à l'aide de la matrice de complaisance $s$ :

$$\mu_{\text{elastique}} = -\frac{1}{2} \Theta(\sigma : s : \sigma) \quad (1.45)$$
Le potentiel chimique généralisé total prend non seulement en compte l'énergie élastique mais aussi le travail des forces de pression soit :

$$ \mu_{total} = -\frac{tr(\sigma)}{3}\Theta - \frac{1}{2}\Theta(\sigma : s : \sigma) \qquad (1.46) $$

Les relations d'Onsager ([LO1931-1], [LO1931-2]) donnent une proportionnalité entre les flux et leur cause :

$$ \vec{J} = -L\vec{\nabla} \left(\frac{\mu}{T}\right) \qquad (1.47) $$

En identifiant (1.47) à la première loi de Fick exprimée en termes de concentration, qui fait intervenir le coefficient de diffusion D :

$$ \vec{J} = -D\vec{\nabla}c \qquad (1.48) $$

et en faisant l'hypothèse d'une loi de Raoult pour le comportement de la solution solide ($\mu = \mu^0 + RT\ln(c)$), on obtient une relation entre L et D [MA2006] :

$$ L = \frac{D}{R\Theta} \qquad (1.49) $$

où R est la constante des gaz parfaits et Θ le volume molaire. On peut ainsi, par la suite, utiliser dans 1.47 cette expression de L avec, pour $\mu$, le potentiel chimique généralisé total.

# 1.5 Évolution des interfaces entre sous-domaines

## 1.5.1 Modélisation des flux

On définit un flux en volume dans chacun des grains, noté $\vec{j}_V^g$. Il est défini dans le cœur du grain sauf la demi-épaisseur $\delta_{GB}/2$ du joint de grain et l'épaisseur $\delta_S$ de la surface libre, comme représenté sur la figure 1.2. De plus, on a défini des flux au joint de grains $\vec{j}_{V,GB}^g$, pour chacun des grains en vis-à-vis, liés chacun au gradient de potentiel chimique dans la demi-épaisseur du joint de grains. On décompose ce flux au joint de grains pour le grain $g$, en deux termes; l'un tangential au joint de grains $\vec{j}_{GB}^g$ et l'autre perpendiculaire : $\vec{j}_{GB}^{\perp g}$. Sa partie tangentielle vérifie :

$$ \vec{j}_{GB}^g = -\frac{\delta_{GB}}{2} D_{GB} \vec{\nabla}_S c^g \qquad (1.50) $$

où $c^g$ est la concentration dans la couche de surface et définit ainsi le coefficient de diffusion au joint de grains.

## 1.5.2 Évolution des interfaces grain/gaz

L'équation de conservation de la matière (seconde loi de Fick) consiste à dire que le flux est nul au travers de l'interface dans le repère de l'interface. Les flux $\vec{j}_S$ et $\vec{j}_S^\perp$ sont définis de manière analogue à $\vec{j}_{GB}$ et $\vec{j}_{GB}^\perp$. Soit $\vec{n}$ le vecteur normal sortant à la surface lissée en un point M donné, pour un solide incompressible la seconde loi de Fick conduit à [MA2006] :

$$ \iint_S \vec{V}_{int} \cdot \vec{n} dS = \Theta \iint_S \left( \vec{j}_S^\perp \cdot \vec{n} + \text{div}_S (\vec{j}_S) \right) dS \qquad (1.51) $$
<figure><img src="image_32.png" /><figcaption>FIGURE 1.2 – Flux utilisés dans le modèle</figcaption></figure>

avec Θ le volume molaire.

Dans le cas où le solide serait compressible la conservation de la masse s'écrirait :

$$ \iint_S \rho \overrightarrow{V_{int}} \cdot \vec{n} dS = M \iint_S \left( \overrightarrow{j_S} \cdot \vec{n} + \text{div}_S \left( \overrightarrow{j_S} \right) \right) dS \quad (1.52) $$

avec M la masse molaire.

L'opérateur div_S met en jeu des dérivées partielles par rapport à une surface pour laquelle on ne dispose pas d'équation (sauf éventuellement à t = 0) et qui n'est pas paramétrée ; la relation entre la métrique de surface et celle d'espace n'est pas évidente à prendre en compte. Les flux de surface et aux joints de grains font intervenir eux-mêmes un second opérateur de surface ; $\overrightarrow{\nabla}_S$ :

$$ \overrightarrow{J_S} = -\delta_S D_S \overrightarrow{\nabla}_S (c_S) \quad (1.53) $$

Ces opérateurs sont utilisés en particulier en mécanique des coques [JG1999].

### 1.5.3 Évolution des joints de grains

Aux joints de grains, l'équation de conservation de la matière devient, pour un solide incompressible, [MA2006] :

$$ \iint_S \rho \overrightarrow{V_{int}} \cdot \vec{n} dS = \Theta \iint_S \left( \left( \overrightarrow{j_{GB}^{\perp 1}} + \overrightarrow{j_{GB}^{\perp 2}} \right) \cdot \vec{n} + \text{div}_S \left( \overrightarrow{j_{GB}^{\perp 1}} - \overrightarrow{j_{GB}^{\perp 2}} \right) \right) dS \quad (1.54) $$

avec Θ le volume molaire.

Soit, pour un solide compressible :

$$ \iint_S \rho \overrightarrow{V_{int}} \cdot \vec{n} dS = M \iint_S \left( \left( \overrightarrow{j_{GB}^{\perp 1}} + \overrightarrow{j_{GB}^{\perp 2}} \right) \cdot \vec{n} + \text{div}_S \left( \overrightarrow{j_{GB}^{\perp 1}} - \overrightarrow{j_{GB}^{\perp 2}} \right) \right) dS \quad (1.55) $$
### 1.5.4 Calcul de la vitesse d'interface

Les équations de la vitesse des interfaces montrent qu'il est nécessaire de savoir calculer le flux volumique $\vec{J}_V$ ainsi que la divergence surfacique du flux de surface. Il n'est donc pas nécessaire de savoir exprimer les flux de surface eux-mêmes.

L'expression du flux en volume est simple [MA2006] :

$$ \forall k \in \{1, 2, 3\}, j_k = \frac{D_V}{3RT} \sum_{i=1}^{3} \frac{\partial \sigma_i}{\partial x_k} + \frac{D_V}{2RT} \sum_{i=1}^{6} \sum_{j=1}^{6} s_{ij} \left( \sigma_i \frac{\partial \sigma_j}{\partial x_k} + \sigma_j \frac{\partial \sigma_i}{\partial x_k} \right) \quad (1.56) $$

où s est la matrice de complaisance (inverse de celle d'élasticité).

Le flux de la divergence surfacique est donné par la circulation du vecteur de flux en surface. Si l'on s'intéresse à un des termes de cette circulation, celui correspondant à l'intégration le long de l'arête [AB] du triangle de la figure 1.3, et pour seulement les termes d'origine élastique, alors, [MA2006] :

<figure><img src="image_33.png" /><figcaption>FIGURE 1.3 – Représentation d'un tétraèdre de l'interface</figcaption></figure>

$$ \int_{[AB]} \vec{j}_{\text{interface}} = l_{AB} \frac{D_{\text{interface}}}{2RT} \sum_{i=1}^{6} \sum_{j=1}^{6} s_{ij} \underbrace{\overrightarrow{\nabla_S (\sigma_i \sigma_j)} \cdot \overrightarrow{n_{AB}}}_{C_{AB}} \quad (1.57) $$

Une méthode qui avait été utilisée dans le cadre des travaux de M. Ajdour pour estimer $C_{AB}$ était d'utiliser la valeur de $\sigma_i\sigma_j$ au centre du cercle circonscrit au triangle, au milieu de l'arête [AB] et au centre du cercle circonscrit au triangle voisin d'arête commune [AB].

Un des théorèmes de Duduchava [RD2006] facilite le calcul de ce gradient puisqu'il peut être simplement calculé en prenant le produit scalaire du gradient usuel avec un vecteur unitaire dans la direction de dérivation. Cette dernière est donnée, aux nœuds du maillage, par le produit vectoriel de la normale à la surface lissée donnée par la bibliothèque CGAL par le vecteur tangent à l'arête du triangle de surface. Comme la circulation du vecteur est à calculer le long de l'arête, le vecteur normal est interpolé entre les deux nœuds, il suffit de savoir estimer le gradient en chaque point le long de l'arête et on peut ainsi ne travailler que sur un élément de surface à la fois.
### 1.5.5 Problème de la conservation de la matière

Les deux relations 1.51 et 1.54 sont des relations intégrales qui donnent l'évolution d'un élément de surface en conservant la matière. Avec la discrétisation liée au maillage, on sait calculer une vitesse pour un triangle $\sigma$ de surface, mais on ne sait pas la calculer en un nœud. Lorsque plusieurs triangles partageant un sommet se déplacent, il est nécessaire de pouvoir donner les nouvelles coordonnées des nœuds. C'est lors de cette dernière étape que l'on risque de perdre ou de gagner de la matière. Plusieurs approches avaient été testées et comparées par M. Wolff [MW2007]. S. Martin a récemment proposé un algorithme pour répondre à ce problème [SM2013]. L'objectif est de remailler le solide déformé à l'issue du pas de temps, à nombre de nœuds de surface constant et en conservant la connectivité.

La solution adoptée ici consiste à déplacer la position d'un sommet le long de la normale à la surface lissée.

En première approximation, si la surface lissée était plane, le déplacement $h \cdot \vec{n}$ du triangle, serait donné par :

$$h = \frac{\|\vec{V}_{int}\| \times \Delta t}{S} \qquad (1.58)$$

Si l'on veut pouvoir utiliser des triangles de surface de dimensions non négligeables devant le rayon de courbure locale, comme au voisinage du joint de grain, les nœuds sont déplacés le long de la normale à la surface de l'instant précédent modifiée par la déformation élastique du pas de temps en cours. Les triangles sont de taille non négligeable devant le rayon local de courbure parce que le raffinement initial du maillage à l'aide de netgen [JS1997] n'a pas été fait en fonction de la courbure moyenne. Il existe une possibilité de le faire dans netgen, pour le maillage initial, en utilisant des fichiers de raffinement .msz qui donnent la valeur du rapport du raffinement local au raffinement moyen au voisinage de points de l'espace (aussi d'arêtes ou de surfaces d'un premier maillage non adapté) dont on donne les coordonnées.

On note $\vec{n}_i$ le vecteur normal au nœud $i$ du triangle, $\theta_i$ l'angle entre $\vec{n}$ et $\vec{n}_i$ comme indiqué sur la figure 1.4.

La nouvelle position des nœuds est calculée comme suit :

$$c'_{i} = c_{i} + \frac{h}{\cos(\theta_{i})} \times \vec{n}_{i} \qquad (1.59)$$

Ensuite, la nouvelle surface du triangle, $S'$ est déduite de la nouvelle position des nœuds $c'_{i}$. $h$ peut alors être corrigé :

$$h' = \frac{\|\vec{V}_{int}\| \times \Delta t}{\frac{S + S'}{2}} \qquad (1.60)$$

Après le calcul du déplacement du triangle $h'$, le déplacement de chacun des nœud est calculé comme la moyenne des déplacements des triangles dont il est le sommet. Par exemple, si un nœud $i$ appartient à $k$ triangles de surface, sa nouvelle position à la fin du pas de temps sera donnée par :

$$c_{i}^{new} = c_{i} + \frac{1}{k} \sum_{l=1}^{k} \frac{h'_{l}}{\cos(\theta_{il})} \vec{n}_{i} \qquad (1.61)$$
<figure><img src="image_34.png" /><figcaption>FIGURE 1.4 – Vitesse des nœuds à la surface des grains</figcaption></figure>

avec $\theta_{il}$ l'angle entre le vecteur normal à la surface du triangle $k$ et le vecteur normal au nœud $i$ et $h'_l$ le déplacement du triangle de surface $l$ calculé par 1.60.

Au final, le déplacement des nœuds est calculé comme suit :

- Tout d'abord, 1.58 est appliquée pour obtenir une valeur initiale de $q$.

- Puis, 1.59 and 1.60 sont utilisées de manière itérative jusqu'à obtenir une conservation acceptable de la matière.

- Finalement, la nouvelle position des nœuds est calculée en utilisant 1.61.

### 1.5.6 Problème de l'évolution de l'angle dièdre

Lors de la résolution des équations de l'élasticité pendant un pas de temps la méthode des éléments finis interdit d'imposer des forces de liaison : les efforts de liaison (force aux points triples) sont obtenus à partir des forces nodales qui suivent le principe fondamental de la statique. Ces forces ne peuvent résulter que des conditions aux limites de type Dirichlet et Neumann imposées à la solution du problème [LG2003].

Ainsi, dans notre modèle, cet angle dièdre n'est pas fixé et est libre d'évoluer. On explique ici quelles sont les condition de Neumann appliquées à la résolution du déplacement qui influent sur la valeur de l'angle dièdre et on montre les résultats de simulation pour cette grandeur.

#### Conditions de Neumann pour les triangles dont une arête appartient à la ligne triple

L'équation 1.35 montre comment est calculée la traction sur une facette qui appartient à l'intérieur d'une surface que ce soit celle du joint de grain ou une surface libre, i.e. dont aucun nœud n'appartient à la ligne triple ou un seul nœud y appartient. Il s'agit simplement d'appliquer la loi de Laplace.
Lorsque la facette σ pour laquelle on veut calculer la traction présente une arête appartenant
à la ligne triple, il est important de calculer sa traction de manière aussi précise que possible. On
note T₁ et T₂ ses nœuds P1 aux extrémités. On considère alors la force de Laplace qui s'exercerait
dessus si cette facette était intérieure à la surface à laquelle elle appartient puis on y soustrait
la contribution le long de l'arête appartenant à la ligne triple (puisqu'on est en bord de surface
et que la partie manquante de la surface ne peut pas tirer sur la facette via cette arête) et on
y ajoute les forces exercées via cette arête triple par les deux autres surfaces partageant cette
arête. L'expression de la traction devient :

$$
\vec{F}_{\sigma} = \iint_{\sigma} \left( -p - \gamma_{\Sigma} \left( \frac{1}{R_1} + \frac{1}{R_2} \right) \right) \overrightarrow{n_{\Sigma}} dS - \int_{T_1 \to T_2} \gamma_{\Sigma} \overrightarrow{n_{T_1 \to T_2}} dl - \int_{T_1 \to T_2} \gamma_{\Sigma'} \overrightarrow{n'_{T_1 \to T_2}} dl - \int_{T_1 \to T_2} \gamma_{\Sigma''} \overrightarrow{n''_{T_1 \to T_2}} dl \quad (1.62)
$$

avec $n_{T_1 \to T_2}$ la normale à l'arête $[T_1 T_2]$ dans le plan tangent à la surface $\Sigma$ et orientée par la
normale sortante $\vec{n}_{\Sigma}$. Les exposants '$\prime$ et '$''$ servent à désigner les deux autres surfaces mises en
jeu.

Cette traction résultante agit sur la globalité de la facette en question et est une condition de
Neumann. Aucune force nodale n'est imposée à la solution numérique.

Géométrie des modèles classiques de frittage et angle dièdre

Les modèles classiques de frittage présentent souvent la géométrie du cou avec un rayon de courbure négatif comme représenté sur la figure 2.1, or, cette géométrie redonnée ici en figure 1.5 est incompatible avec l'existence d'un joint de grain à l'équilibre.

<figure><img src="image_35.png" /><figcaption>FIGURE 1.5 – Représentation de l’angle dièdre des modèles classiques de frittage [DBA2005]</figcaption></figure>

En effet, l'angle dièdre $\Theta$ serait alors de 180° et l'équilibre des forces au point triple : $\gamma_{GB} = 2\gamma_S \cos(\Theta/2)$ conduirait à une tension de joint de grain nulle i.e. à l'absence de joint de grain.

Résultats de simulation : angle dièdre

Bien que les résultats du modèle soient présentés au chapitre suivant, comme on ne peut
pas donner d'expression analytique de l'évolution de l'angle dièdre, on donne ici l'allure de son
évolution via des résultats numériques. La figure 1.6 montre les résultats d'une simulation de
cette grandeur sur deux grains sphériques de cuivre en conditions isothermes (à 1100°C) avec
notre modèle. On voit que la valeur d'équilibre est atteinte au bout de quelques pas de temps.
On remarque, sur cette figure, la dispersion des points de simulation. Elle est due à la taille du
maillage et diminue lorsque la finesse de se dernier augmente.

On voit que, partant d'une géométrie qui n'est pas à l'équilibre thermodynamique, les surfaces
vont évoluer jusqu'à annuler les flux de matière, et à chaque étape l'équilibre mécanique élastique
est atteint (équilibre quasi-statique). Ici $\gamma_{GB} = 0,854103J \cdot m^{-2}$ et $\gamma_{SV} = 1,65J \cdot m^{-2}$ pour le
<figure><img src="image_36.png" /><figcaption>FIGURE 1.6 – Évolution de l'angle dièdre dans notre modèle, cas du cuivre à 1100°C</figcaption></figure>

cuivre à 1100°C [WZ1995]. On atteint bien l'angle théorique d'équilibre des tensions de surface soit $\Theta = 150°$.

### 1.5.7 Qualité du maillage de surface

Une fois la position des nœuds mise à jour, elle est passée, dans le format "STereoLithography" (STL), au mailleur GMSH [CG2009] pour mailler en volume. Ce remaillage ne modifie ni la position ni la connectivité des nœuds de surface.

Aussi, l'étape de remaillage devient critique lorsque les triangles deviennent irréguliers. Ceci se produit au voisinage des joints de grains là où les déformations sont les plus importantes. Une solution provisoire a été développée et mise en œuvre par Sylvain Martin (thèse en cours cf Annexe C p.133) pour modifier la triangulation de surface lorsque les triangles deviennent aplatis.

Les nœuds de surface sont classés en couches correspondant à leur distance à la ligne triple.

Si l'un des angles d'un triangle testé est en dessous d'une valeur limite, le sommet opposé, le plus loin du joint de grain est déplacé. Pour conserver la matière et la forme, le déplacement du nœud est calculé dans le plan du triangle, puis projeté sur le plan défini par la normale au nœud considéré.

Une autre approche consisterait à remailler en surface à partir des nœuds seuls, sans leur connectivité, et en imposant une continuité de la courbure d'un nœud à l'autre. Un outil potentiel serait le "surface evolver" de Brakke [KB1992].

D'autres auteurs ont opté pour une méthode de type level-set [JB2011].

## 1.6 Choix du pas de temps

La résolution se faisant en explicite, la stabilité du schéma en temps a été étudiée par Marc Wolff dans le cadre de son stage de Master I sur la modélisation du frittage au sein du LMPC ([MW2007] - annexe C, p.133) :

$$ \Delta t \le \frac{1}{2} \frac{h^2}{D} \qquad (1.63) $$
Ces conditions dites de Courant, Friedrichs, Lewy (CFL) [CFL1928], dites CFL devraient assurer
la convergence pour la diffusion. Comme plusieurs coefficients de diffusion interviennent D<sub>V</sub>, D<sub>S</sub>
et D<sub>GB</sub> le plus grand parmi D<sub>V</sub>, δD<sub>S</sub> et δD<sub>GB</sub> , où δ est l'épaisseur de la couche de surface ou
la demi-épaisseur du joint de grain, sera pris en compte pour limiter le pas de temps.

Une seconde limitation du pas de temps intervient de sorte que les maillages de deux grains
au contact ne puissent pas se rencontrer en dehors du joint de grain.

L'enchaînement global des calculs est illustré par la figure 1.7 avec le détail du calcul des
déplacements à l'aide de diffpack donnés par la figure 1.7.

1.7 Conclusion

Un certain nombre de points sont à étudier pour que les résultats du modèle puissent refléter
au mieux les phénomènes que l'on essaie de simuler.

Une question est de savoir comment extrapoler correctement les rayons de courbure aux
points de la ligne triple et aussi, de manière générale, ailleurs sur la surface, définir la taille du
voisinage d'un point à utiliser pour avoir une précision visée sur la courbure moyenne locale
calculée.

L'algorithme de décomposition de domaines (de type Marini-Quarteroni), qui a été programmé et utilisé dans la thèse d'Ajdour [MA2006] pourra être remis en œuvre pour s'assurer qu'il donne les mêmes résultats que sur deux grains, en vue de passer à quelques grains, et de pouvoir comparer aux résultats obtenus par éléments discrets.
La formulation faible du problème global en espace et en temps serait nécessaire pour formaliser correctement le problème.

Bien que la chimie du problème n'ait pas été exposée ici, la progression par étapes voulant
que l'on commence par des modèles plus simplifiés à composition constante, un point n'a pas
été appréhendé : il s'agit de la gestion de la composition locale de la phase gaz permettant
de simuler l'évaporation/condensation. La difficulté dans cette perspective serait d'une part de
savoir déterminer automatiquement le caractère ouvert ou fermé d'un pore et d'autre part de
gérer les pas de temps, la diffusion s'effectuant bien plus rapidement en phase gazeuse qu'en
phase solide. Pour le premier point, on peut penser qu'en décrivant l'ensemble des grains en
termes d'intersection et d'union de demi-espaces le mailleur netgen puisse déterminer le nombre
de domaines et ainsi le nombre de pores fermés.
<figure><img src="image_37.png" /><figcaption>FIGURE 1.7 – Enchaînement des étapes de calcul</figcaption></figure>
<figure><img src="image_38.png" /><figcaption>FIGURE 1.8 – Utilisation des classes de Diffpack pour le calcul des déplacements et de leurs dérivées lissées</figcaption></figure>
# Résultats du modèle - composition constante

## Sommaire

<table><tr><td>2.1</td><td>Rayons de courbure principaux</td><td>81</td></tr><tr><td>2.2</td><td>Champ de déplacement</td><td>82</td></tr><tr><td>2.3</td><td>Évolution de la taille du cou en fonction du temps</td><td>82</td></tr><tr><td>2.3.1</td><td>Diffusion en surface et aux joints de grains</td><td>83</td></tr><tr><td>2.3.2</td><td>Diffusion en volume</td><td>84</td></tr><tr><td>2.4</td><td>Conclusion</td><td>86</td></tr></table>

On a donné un nombre limité de résultats de simulation qui montrent que les résultats ne semblent pas aberrants.

## 2.1 Rayons de courbure principaux

Les rayons de courbure principaux $R_1$ et $R_2$ calculés ont été visualisés à l'aide du logiciel paraview. Ils sont représentés sur la figure 2.1. L'écart aux valeurs théoriques est calculable à

<figure><img src="image_39.png" /><figcaption>FIGURE 2.1 – Rayons de courbure principaux $R_1$ et $R_2$</figcaption></figure>
l'instant $t = 0$ sur cette géométrie très simple. Il peut être réduit en augmentant la finesse du maillage.

## 2.2 Champ de déplacement

Ce champ est calculé en résolvant le système 1.32 p. 68. La méthode itérative utilisée est la "minimum residual method for symmetric, indefinite systems" [HPL2003], bien que dans notre cas le système soit défini<fnref n="2" />. Le moniteur de convergence utilisé est $\rho_k = \frac{|B(b-Ax_k)|_B}{|B(b-Ax_0)|_B}$ en se fixant un seuil de $10^{-3}$. La figure 2.2 donne l'évolution du moniteur de convergence avec le nombre d'itérations.

<figure><img src="image_40.png" /><figcaption>FIGURE 2.2 – Vitesse de convergence de l’algorithme de calcul des déplacements – cas d’un maillage avec 4700 noeuds par grain</figcaption></figure>

On a représenté, figure 2.3, l'opposé du champ de déplacements $\vec{u}$ afin de mieux le visualiser. On voit que les sphères sont en compression tandis que le cou n'a pratiquement pas de contrainte car ses rayons de courbure principaux se compensent et conduisent à une courbure presque nulle. Mais il ne s'agit que du champ de déplacements élastiques pas de celui des déformations irréversibles à l'issue d'un pas de temps. Une simulation faite en forçant $c_{44} = \frac{1}{2}(c_{11} - c_{12})$, afin de se placer un cas isotrope, présente (figure 2.3) une symétrie de révolution qui semble correcte.

## 2.3 Évolution de la taille du cou en fonction du temps

Des simulations ont été menées sur deux particules de même taille en contact. Le maillage initial, pour une géométrie de deux sphères tronquées dont le rapport du rayon $r$ du cou à celui

<footnote n="2">2. Elle consiste à résoudre un système du type $BAx = Bb$, où $A$ est symétrique non-inversible et $B$ est symétrique définie positive. Le système linéaire est symétrique par rapport au produit scalaire $((\cdot, \cdot)_B = (B^{-1}\cdot, \cdot))$. Dans notre cas, $A$ est inversible mais la méthode travaille par blocs ce qui est intéressant. L'itérée $x_k$ est déterminée en minimisant la norme du résidu $|B(b - Ay)|_B$ sur l'espace vectoriel de Krylov $K_k(AB, Bb)$ engendré par $Bb$ et ses $k-1$ images par les puissances successives de $BA : Vect\{Bb, BABb, ..., (BA)^{k-1}Bb\}$. Ici, la norme est définie par le produit scalaire $((\cdot, \cdot)_B)$. Le moniteur de convergence par défaut est $\rho_k = |B(b - Ax_k)|_B = (B(b - Ax_k), b - Ax_k)$. Le résidu $(b - Ax_k)$ n'est pas calculé durant l'itération, aussi, un calcul direct de cette quantité nécessite un produit matrice-vecteur supplémentaire. Dans l'algorithme, cette quantité est calculée de manière récursive. Malheureusement, ces calculs accumulent de l'erreur et il est nécessaire, à intervalles réguliers, de calculer le résidu exact.</footnote>
<figure><img src="image_41.png" /><figcaption>FIGURE 2.3 – Opposé du champ de déplacements élastiques $\vec{u}$ à gauche : cas isotrope, à droite : cas cubique - rayon initial des grains : 1,1µm</figcaption></figure>

$R$ de la sphère est proche de 0,245, a été généré par netgen. Le joint de grains initial est le plan de symétrie des deux grains perpendiculaire à l'axe de révolution. Les résultats qui suivent sont issus des travaux de thèse de S. Martin [SM2013] dont je suis un co-encadrant.

### 2.3.1 Diffusion en surface et aux joints de grains

La modélisation du frittage, pour deux particules sphériques de même taille, a été étudiée par différents auteurs en considérant la diffusion en surface couplée à la diffusion aux joints de grains [JP1998], [JS1995], [WZ1998], [DB1996]. Ces travaux ont débouché sur une loi analytique de la forme $\frac{r}{R} = (\alpha \frac{t}{\tau_{gb}})^{\frac{1}{6}}$ [DB1996], [FP1998], [FP1999], [PC2008], donnant une bonne description des phénomènes, utilisable pour la validation de notre modèle. Les principaux paramètres en sont le rayon du cou $r$, le rayon du grain $R$ et le temps adimensionnel $\frac{t}{\tau_{gb}}$ avec $\tau_{gb} = \frac{R^4 k T}{\Theta \delta_{GB} D_{GB} \gamma_S}$. Pour un rapport $\frac{\delta_{GB} D_{GB}}{\delta_S D_S} = 1$ qui correspond aux conditions de notre simulation, Parhami [FP1999] a donné $\alpha = 192$.

Comme on peut le voir sur la figure 2.4, nos simulations donnent des résultats très proches de ceux de Parhami. De petits écarts sont, cependant, observables et peuvent provenir de la valeur trop élevée du rayon initial du cou ($r_0 = 0,049\mu m$) vis-à-vis du rayon initial des grains ($R_0 = 0,200\mu m$) ou du manque de finesse de notre maillage. Pour rendre compte correctement de l'évolution de la surface, il est nécessaire de faire un compromis entre la finesse du maillage, qui impacte directement le temps de calcul, et la préservation de la forme des grains. De plus, la valeur du rayon de courbure est très sensible à la finesse du maillage. Au tout début du frittage la surface au voisinage du joint de grains subit une modification de courbure qui peut seulement être obtenue au moyen d'un maillage fin, le raffinement local du maillage serait à ajuster en fonction de la courbure. Au cours de la simulation, les nœuds de surface se déplacent conduisant à des éléments de surface (triangles $P_2$) irréguliers qui ne sont pas convenables à une méthode des éléments finis. L'algorithme actuel de contrôle du remaillage, implémenté dans le logiciel, modifie les nœuds de surface en les déplaçant à connectivité constante et retarde l'apparition de ce phénomène.

Le retrait $h$ peut être évalué d'un point de vue géométrique. Avec les hypothèses de Coble [RLC1958], l'évolution peut être décrite à l'aide de sphères tronquées. Avec l'hypothèse d'un
<figure><img src="image_42.png" /><figcaption>FIGURE 2.4 – Comparaison de l'évolution du rayon du cou (travaux thèse S. Martin) à une solution analytique [WC1980] ($\alpha$ = 192)</figcaption></figure>

rapprochement entre grains petit devant leur rayon initial, une relation entre le rapprochement des centres, noté ici $\dot{h}$ (la distance entre centres étant $2R - h$) et le rayon $r$ du cou est donnée par [FP1999]-équation (46) :

$$ \frac{r - r_0}{R} = 2\sqrt{\frac{h - h_0}{R}} \qquad (2.1) $$

L'évolution du rayon du cou en fonction de la distance d'approche centre à centre est donnée en Figure 2.5 et on retrouve le cou prédit par le modèle géométrique. On note, cependant, de petits écarts :

- d'une part $h$ est légèrement négatif sur les premiers pas de temps,

- d'autre part, le rayon du cou est légèrement sous-estimé vis-à-vis du retrait.

On peut attribuer ces deux effets à un manque de finesse du maillage.

### 2.3.2 Diffusion en volume

La diffusion en volume a généralement peu d'impact sur le comportement des matériaux homogènes, car elle est bien plus lente que la diffusion en surface ou aux joints de grains. Néanmoins elle est à considérer dans le cas des matériaux de composition variable. Son effet a été étudié dans le cas d'une composition homogène par Bruchon [JB2010] et Pino-Muñoz [JB2013], [DHP2012], mais ces auteurs ont considéré un milieu continu sans joint de grain.
<figure><img src="image_43.png" /><figcaption>FIGURE 2.5 – Evolution du rayon du cou avec la distance centre à centre. Comparaison à Parhami [FP1999] (travaux thèse S. Martin)</figcaption></figure>

Dans nos simulations, l'évolution du rayon du cou est en bon accord avec la loi cinétique donnée par Coble [RLC1958] pour la diffusion en volume quand la source de matière est le joint de grains :

$$ \left(\frac{r}{R}\right)^n = \frac{FD_v\gamma\Theta}{R^m k T} \times t \qquad (2.2a) $$

$$ \left(\frac{h}{2R}\right)^p = \frac{BD_v\gamma\Theta}{R^m k T} \times t \qquad (2.2b) $$

avec $F = 32, m = 3, n = 4, p = 2$ et $B = 2$.

L'évolution du rayon du cou est représentée dans la figure 2.6 et l'approche centre à centre (h) sur la figure 2.7. L'évolution du rayon du cou suggère que le phénomène prédominant est la diffusion en volume (du joint de grains vers le cou), cependant, la courbe de retrait ne correspond pas bien à ce qui est attendu. Il faut noter que h représente le retrait des grains, qui correspond au déplacement des centres de sphères tronquées dans le cas de Coble [RLC1958] alors que, dans notre modèle, c'est le rapprochement relatif des centres de gravité des grains. Ceci explique en partie que le retrait soit plus petit dans nos résultats.
<figure><img src="image_44.png" /><figcaption>FIGURE 2.6 – Évolution du rayon du cou. Comparaison au modèle de Coble [RLC1958]</figcaption></figure>

<figure><img src="image_45.png" /><figcaption>FIGURE 2.7 – Évolution de la distance centre à centre en fonction du temps. Comparaison au modèle de Coble [RLC1958]</figcaption></figure>

## 2.4 Conclusion

Les premiers résultats donnant le rayon relatif du cou en fonction du temps sont cohérents avec les lois analytiques de la littérature (travaux de Coble) ainsi qu'avec les résultats des modèles numériques de frittage comme ceux des travaux de Bouvard et de Parhami. Néanmoins, l'évolution du cou des grains est sensible au raffinement local du maillage et à la forme initiale du grain. Lorsque les nœuds de surface se déplacent lors de la simulation, la conservation de leur nombre et de leur connectivité induit des éléments de surface irréguliers. Pour retarder l'apparition de ce problème un algorithme de modification du maillage de surface a été mis en place [SM2013]. La surface des grains nécessiterait d'être remaillée à partir des nouvelles positions des nœuds avec un raffinement adapté à la courbure locale. La façon dont la surface est remaillée à la fin d'un pas de temps nécessite des améliorations pour rendre compte des déformations importantes qui surviennent dans les premiers instants du frittage.
dans les fours de frittage et la littérature [DFS1989] indique que c'est sous forme atomique que H diffuserait dans $UO_2$. Sa présence n'est pas anodine car il semblerait qu'utiliser $H_2/H_2O$ plutôt que $CO/CO_2$ donnerait, sous certaines conditions, des microstructures à plus gros grains.
A l'intérieur des grains $g$ et $g'$ au contact on a :

$$ \forall i \in \{1, \dots, 5\} \quad \left\{ \begin{array}{l} \mathrm{div} \left( \vec{j}_{V_i}^g \right) + \frac{\partial c_i^g}{\partial t} = 0 \\ \mathrm{div} \left( \vec{j}_{V_i}^{g'} \right) + \frac{\partial c_i^{g'}}{\partial t} = 0 \end{array} \right. \qquad (3.1) $$

Il y a par ailleurs égalité des potentiels de part et d'autre de l'interface :

$$ \forall i \in \{1, \dots, 5\} \quad \mu_i^{\text{solide}} (c_1^g, \dots, c_5^g, T_1^g, \dots, T_5^g) = \mu_i^{\text{solide}} (c_1^{g'}, \dots, c_5^{g'}, T_1^{g'}, \dots, T_5^{g'}) \quad (3.2) $$

où $T_i$ désigne la i-ième composante du tenseur des contraintes en notations de Voïgt. A cela viennent s'ajouter les équations de flux nuls au travers des interfaces. Ces dernières équations aux dérivées partielles sont résolues par un schéma de type volumes finis 2D sur les variétés que sont les surfaces des grains avec des conditions aux limites de type Dirichlet données par la résolution de l'égalité des potentiels aux interfaces.

Claire Négrel (Annexe C, p.133) a proposé un schéma de discrétisation de l'ensemble mécanique et chimie dans son rapport de stage de maîtrise de maths appliquées effectué au laboratoire.

Le schéma effectivement programmé est celui proposé par Robert Boyer en 2000, mais il n'a jamais été testé ni débogué jusqu'à présent. La conservation de la matière devrait, de plus, tenir compte de la variation des volumes atomiques avec la composition [YCC2004].

## 3.3 Études visant à rendre le $(U, Pu)O_2$ plus homogène

Des études expérimentales ont été menées et sont en cours dans des cas de figure beaucoup plus complexes que ceux qui ont pu être décrits dans le modèle précédent. On étudie les mécanismes de certains adjuvants de frittage tel $Cr_2O_3$.

J'ai encadré de 2004 à 2007 la thèse d'Adrien Pieragnoli sur l'étude de l'action de $Cr_2O_3$ sur la microstructure des pastilles de combustible, notamment, au niveau de leur répartition en plutonium.

Son étude a pu mettre en évidence la nature des mécanismes d'action de $Cr_2O_3$ sur l'évolution, lors du frittage, de la microstructure d'une pastille $(U, Pu)O_2$. Dans un premier temps, notre observation des phénomènes d'interaction intervenant au sein d'une pastille $(U, Pu)O_2$ do-pée $Cr_2O_3$, que l'on peut se représenter en cru comme un matériau composite constitué d'une matrice de $UO_2$ contenant des amas de 50 à 100 µm de mélange-mère de $(U, Pu)O_2$ avec, dans un des deux constituants ou dans les deux, de très faibles quantités (de quelques dixièmes de pour-cent) de $Cr_2O_3$, a été simplifiée en réalisant un système modèle constitué des mêmes matériaux que ceux du composite, traités thermiquement dans des conditions similaires de frittage, mais de dimensions millimétiques au lieu de micrométriques, afin de mieux pouvoir en étudier les interfaces (pastilles dites gigognes).

On a ainsi pu mettre en évidence la formation du composé de type pérovskite $(Pu_{1-\epsilon}U_{\epsilon})CrO_3$ aux joints de grains de la solution solide $(U, Pu)O_2$. L'étude en Spectroscopie d'Absorption X de ce composé à l'European Synchrotron Radiation Facility (ESRF) a montré que les deux cations
de la pérovskite étaient au degré d'oxydation +III.

On a pu montrer que le coefficient d'interdiffusion en volume entre $UO_2$ et $(U, Pu)O_2$, bien qu'estimé d'une façon peu précise, était augmenté d'un ordre de grandeur en présence de chrome.

La solubilité du chrome a été mesurée dans les conditions des essais dans la matrice $(U, Pu)O_2$ de structure fluorine. Ces conditions représentaient une gamme restreinte de conditions thermodynamiques $(\mu_{O_2}, T)$. Une telle gamme a été élargie dans le cadre de la thèse de Régis Thomas.

Dans un second temps, l'ensemble des informations acquises a permis de proposer une interprétation de l'évolution, au cours du frittage, de la microstructure d'une pastille $(U, Pu)O_2$ en présence de chrome. L'interprétation que l'on en a faite est la suivante : il serait possible que le chrome forme un film liquide de composition $CrO$ au niveau des joints de grains de $(U, Pu)O_2$ les ions $Pu^{4+}$ au contact de ce film sont réduits en $Pu^{3+}$, tandis que $Cr^{2+}$ dans la phase liquide est oxydé en $Cr^{3+}$ en formant $PuCrO_3$ qui précipite. La surface du grain de $(U, Pu)O_2$ étant appauvrie en Pu, un gradient de teneur Pu se créerait qui va entraîner une interdiffusion de U et Pu et une homogénéisation de la concentration en Pu. Le film liquide en question n'a, cependant, pas été mise en évidence. Cependant, en vue de répondre à cette interrogation on a préparé deux échantillons fortement dopés en $Cr_2O_3$ puis on a réalisé deux fractures : l'une intergranulaire (réalisée à basse température), l'autre intergranulaire (réalisée à haute température), afin de pouvoir observer une éventuelle différence entre ces surfaces de fracture (volume du grain et joints de grains) en spectroscopie XPS-UPS. Le caractère intra ou intergranulaire a été vérifié après révélation des joints de grains sur l'autre demi-pastille. La spectroscopie XPS-UPS devrait est réalisable au CEA de Valduc. Le $PuCrO_3$ formé peut, lors de la poursuite du traitement thermique, se redissoudre par diffusion du chrome de cette phase vers le reste de la matrice $UO_2$ qui n'en contient pas. Le Pu sous forme 3+ diffuse plus vite aux joints de grains de $UO_2$ que ne le ferait $Pu^{4+}$, mais diffuse ensuite moins vite en volume dans les grains de $UO_2$.

Ces travaux ont été poursuivis dans le cadre d'une seconde thèse, celle de Régis Thomas [RT2013], dont le but a été de déterminer plus précisément les propriétés thermodynamiques de $PuCrO_3$ en déterminant expérimentalement les limites de son domaine d'existence. Cette thèse a aussi permis à une collègue du Département - Chantal Riglet-Martial - d'établir un modèle de solubilité du chrome dans une matrice $(U, Pu)O_2$ avec, suivant les conditions thermodynamiques, une précipitation de $Cr$ métallique ou $Cr_2O_3$ ou $PuCrO_3$ lorsque la limite de solubilité est atteinte. J'ai participé à la rédaction d'un protocole de transport entre l'Installation Nucléaire de Base qu'est notre laboratoire et celle qui contient la ligne INE sur le synchrotron ANKA à Karlsruhe en Allemagne. Le but des essais à ANKA est de déterminer par Spectroscopie d'Absorption X les degrés d'oxydation des cations $U, Pu$ et $Cr$ dans nos échantillons, ainsi que la nature des ions dans les couches de premiers voisins ainsi que leur distance à l'ion central considéré. La localisation du chrome dans les sous-réseaux ioniques de la matrice $(U, Pu)O_2$ est d'un intérêt capital pour la modélisation de sa solubilité. L'intérêt de cette ligne de rayons X est le fait qu'elle accepte des quantités assez importantes de matière nucléaire (des fragments de pastilles de $(U, Pu)O_2$). Les résultats de SAX ont montré que le chrome dans la matrice était sous forme de défauts complexes. Bien que la transformée de Fourier n'ait pu être menée que sur une étendue limitée, on a pu voir que les ions chromes avaient trois autres chromes à égale distance. Ils se disposent donc aux sommets d'un tétraèdre et sont entourés de plus proches voisins qui sont des ions oxydes. Si cette structure était plus étendue, il s'agirait probablement de feuillets d'oxyde de chrome insérés dans la matrice fluorine et probablement en épitaxie avec celle-ci. Dans
le cadre de cette thèse, j’ai pu contribuer à une quantification de l’homogénéité du matériau à partir de cartographies semi-quantitatives des éléments U et Pu et de profils quantitatifs pris dans celles-ci en déterminant des distributions de niveaux de teneur en élément Pu.

La mise en œuvre du chrome dans le frittage est trop complexe pour être modélisée dans le cadre du modèle à l'échelle sub-granulaire de frittage développé du fait des cinétiques très différentes entre l'interdiffusion entre U et Pu et la diffusion de l'élément Cr dans la matrice (U, Pu)O₂ et de la nécessité de pouvoir simuler l'apparition de phases secondaires telles que PuCrO₃, Cr, Cr₂O₃ voire CrO.
Conclusion

Un modèle de type Lagrangien a été développé à une échelle sub-granulaire pour modéliser
le frittage du combustible. Initialement (en 1999) il avait été programmé par Dassault Data
Services (DDS) sur la base du modèle qui avait été établi et incluait d'une part des algorithmes
de décomposition de domaines de type Dirichlet-Neumann préconisés par le Centre de Mathé-
matiques et d'Informatique de Marseille (CMI)/Laboratoire d'Analyse, Topologie, Probabilités
(LATP)/Équipe équations aux Dérivées Partielles (EDP) avec une composition chimique va-
riable. Le calcul de déplacement fonctionnait sur plusieurs grains avec la mécanique seule, et la
partie diffusion chimique en phase solide avait été programmée mais pas mise en œuvre.

Dans le cadre de la thèse de Mounia Ajour, de 2003 à 2006, en collaboration avec l'École
des Mines de Saint-Étienne (ENSME-SE)/Département Science et Génie des Matériaux (SGM),
nous sommes revenus à deux grains seulement et une composition fixe afin de pouvoir valider
les flux de matière et la déformation irréversible des contours des grains en comparant avec les
résultats des modèles classiques.

Des travaux effectués dans le cadre d'une collaboration avec le laboratoire Science des Pro-
cédés Céramiques et de Traitement de Surface (SCPTS) de Limoges et l'Université de Sofia en
Bulgarie, thèse de Radostin Blazhev, nous ont permis d'étendre le modèle au cas du frittage
Spark Plasma Sintering (SPS<fnref n="3" />) en calculant le champ de potentiel électrique, le champ de tem-
pérature en résultant du fait de l'effet Joule, puis la déformation due au frittage dans ce champ
de température. Cette simulation est faite sur deux grains de composition constante.

Dans le cadre de la thèse de Sylvain Martin (UTC-Université Jules Verne de Picardie - 2011-

<footnote n="3">3. La technique SPS n'est classiquement pas utilisée pour le combustible nucléaire. Toutefois, récemment, trois publications relatives à l'utilisation de la méthode SPS pour le frittage de UO₂, d'un cermet de UO₂ et d'un mélange UO₂ - SiC est apparue dans la littérature :</footnote>

<footnote>- Densification of uranium dioxide fuel pellets prepared by spark plasma sintering (SPS), Lihao Ge, Ghatu Subhash, Ronald H. Baney, James S. Tulenko, Edward McKenna, Journal of Nuclear Materials, Volume 435, Issues 1-3, April 2013, Pages 1-9, http://dx.doi.org/10.1016/j.jnucmat.2012.12.010</footnote>

<footnote>- Spark Plasma Sintering of W - UO₂ cermets, O'Brien, R. C. ; Jerred, N. D., Journal of Nuclear Materials, Volume 433, Issue 1-3, p. 50-54, DOI : 10.1016/j.jnucmat.2012.08.044</footnote>

<footnote>- Fabrication Strategies and Thermal Conductivity Assessment of High Density UO2 Pellet Incorporated with SiC, 2012 MRS Spring Meeting, Sunghwan Yeo, Edward Mckenna, Ronald Baney, Ghatu Subhash and James Tulenko, in Actinides and Nuclear Energy Materials Volume 1444, Éditeur Corwin H. Booth, David Andersson, Peter C. Burns, Materials Research Society, 2012</footnote>
2014), les bogues ou erreurs dans le modèle à l'échelle sub-granulaire sur deux grains ont été corrigés afin d'en déduire une force de frittage qu'il a intégrée à un modèle de contacts discrets (MULTICOR) en vue de passer à une échelle supérieure. Dans le cadre de cette thèse, des progrès substantiels ont été réalisés concernant la prise compte correcte de la conservation de la matière liée aux échanges entre grains et aux raffinements dynamiques du maillage.

L'ensemble de ces démarches est schématisé en figure 4.1.

<figure><img src="image_46.png" /><figcaption>FIGURE 4.1 – Évolution du développement du code de frittage SALAMMBO</figcaption></figure>

Les résultats à composition chimique constante ont été présentés.

La chimie des phénomènes liés au frittage du combustible a été étudiée en expérimentalement. Ces études expérimentales ont été mesurées en collaboration avec l'École de Céramiques de Limoges (ENSCI), l'Institut de Chimie de la matière condensée de Bordeaux (ICMCB), l'Ecole Nationale Supérieure de Chimie, de Biologie et de Physique (ENSCBP). Elle se révèle suffisamment complexe pour ne pas pouvoir être représentée par une simple interdiffusion.

La partie suivante présente les moyens d'acquisition des grandeurs d'entrée du modèle à l'échelle sub-granulaire.
# Quatrième partie

Acquisition des paramètres d'entrée du modèle à l'échelle sub-granulaire
<figure><img src="image_47.png" /></figure>
Une étude bibliographique avait été menée en 1998 afin de recenser les moyens d'acquisition au laboratoire des grandeurs d'entrée du modèle [JL1998]. On liste ici ces grandeurs, leurs principales méthodes d'acquisition et, pour les coefficients de diffusion en volume et aux joints de grains ce qui a été réalisé au laboratoire depuis.
1

Tension de surface

Dans sa revue bibliographique Fink [JKF2000] recommande pour $UO_{2\pm x}$ entre 0°C et 2850°C l'expression suivante pour la tension de surface solide/gaz :

$$ \gamma_{SV} = 0,85 - 1,40 \ 10^{-4} (T - 273) \qquad (1.1) $$

avec T la température en Kelvin et, pour le cas non-stoechiométrique :

$$ \gamma_{SV}|_{x} - \gamma_{SV}|_{0} = 6,8x \qquad (1.2) $$

avec $0 \le x \le 0,05$ et $T \le 2170K$.

Plus récemment, Castell [MC2003] a étudié expérimentalement le diagramme de Wulff de $UO_2$, en fabriquant pour cette étude, des monocristaux centimétriques. Certains clivages de ces monocristaux le long des plans (111) montrent des porosités de quelques micromètres. L'utilisation du MEB à basse tension montre qu'ils sont délimités par des facettes {111} et {001} et certains ont atteint leur forme liée à l'équilibre thermodynamique (forme de Wulff). Une analyse de la surface des facettes montre que les pores les plus petits sont proches de la forme de Wulff, et que le rapport de leurs énergies de surface vaut $\frac{\gamma_{001}}{\gamma_{111}} = 1.42 \pm 0.05$. Les pores les plus larges montrent une accumulation de marches à l'intersection des facettes avec des structures ordonnées de dizaines de nanomètres de période.

Boyarchenkov [ASB2012] a récemment calculé, à l'aide de la dynamique moléculaire, ces grandeurs à 1927°C et 2027°C pour $UO_2$ et retrouve les valeurs de Castell, à savoir : $\gamma_{111} = 0.875 \pm 0.008 Jm^{-2}$.

Aucune étude expérimentale n'a, à ma connaissance, été réalisées sur le $(U, Pu)O_2$. Une telle étude pourrait être réalisée par indentation de Hertz [TI1981] au microscope optique ou avec un suivi acoustique des fissures générées par l'indenteur. Les études sur $UO_2$ ont montré que la température avait peu d'impact sur ce paramètre, mais que l'écart à la stœchiométrie en avait. Il semble probable que la substitution de U par Pu ait un impact sur la valeur de la tension de surface.

Il est à noter que l'étude de la tension de surface n'intéresse pas seulement le frittage, mais aussi le comportement des bulles de produits de fission gazeux, donc le comportement du combustible lors de sa fabrication et sous irradiation.
Constantes élastiques isothermes

Les constantes élastiques isothermes ont été mesurées à température ambiante sur $UO_2$ ainsi que, plus récemment, sur du combustible irradié [DL2002]. A haute température les études sont plus anciennes [MOM1969]. Les calculs récents de dynamique moléculaire jusqu’à 3500K donnent d’excellents accords avec l’expérience sur des données dérivées telles le module d’Young, le coefficient de Poisson et le module de Coulomb. Les modules élastiques obtenus par calcul seraient donc suffisants, même à haute température pour prévoir les propriétés du matériau [SIP2011].
côtés gauche et droite respectivement. On a préféré utiliser $ln(λ_5)$ et $ln(λ_3)$ plutôt que $λ_5$ et $λ_3$ dans le logiciel développé car les coefficients de diffusion sont très petits et la perte de précision de la représentation machine amenait sinon à des aberrations. En effet, plus le coefficient de diffusion est petit et plus $λ_3$ et $λ_5$ sont grands. On ne détaille pas ici les propriétés de cette

<figure><img src="image_48.png" /><figcaption>FIGURE 3.1 – Concentrations brute et expérimentale dans une expérience d'interdiffusion entre $UO_2$ et $(U, Pu)O_2$</figcaption></figure>

fonction qui ont été utilisées, l'ensemble est donné dans une publication [JL2012]. On précise seulement qu'en utilisant cette fonction ajustée sur nos échantillons, on a pu mettre en œuvre trois méthodes d'estimation des coefficients d'interdiffusion : celle de Boltzman-Matano, celle de Den Broeder et celle de Hall spécifique aux extrémités du profil.
Conclusion

Toutes les grandeurs d'entrée du modèle sont mesurables expérimentalement et certaines donnent un accord suffisamment bon avec l'expérience pour pouvoir être calculées pour d'autres conditions expérimentales.

La partie suivante présente les techniques qui ont été testées comme outils de validation du
modèle de frittage à différentes échelles.
Cinquième partie

Étude de moyens de validation
expérimentale
<figure><img src="image_49.png" /></figure>
# 1

# Moyens de validation à l'échelle sub-granulaire

## Sommaire

<table><tr><td>1.1</td><td>μ-XANES et μ-tomographie</td><td>107</td></tr><tr><td>1.2</td><td>MEB environnemental</td><td>111</td></tr><tr><td>1.3</td><td>Autres pistes</td><td>112</td></tr></table>

Ce chapitre décrit les voies qui ont été explorées afin de déterminer des techniques qui puissent permettre de valider le modèle à l'échelle d'un agglomérat de poudre de quelques micromètres cubes ainsi qu'à l'échelle de la pastille.

## 1.1 μ-XANES et μ-tomographie

Des expériences de μ-XANES et de μ-tomographie X haute résolution en transmission ont été réalisées à l'ESRF (European Synchrotron Radiation Facility) sur la ligne ID22 afin de tester ces techniques comme outils d'étude de l'évolution d'agglomérats isolés pris à différents stades du frittage du point de vue de :

- leur densité (validation du frittage lui-même),

- la nature des cations (U ou Pu) (validation du frittage avec chimie),

- le degré d'oxydation des cations.

Deux groupes d'échantillons ont été utilisés :

- des agglomérats de 50 µm constitués de cristallites de 0.3µm (agglomérats de poudre crue oxydés par l'oxygène et la vapeur d'eau à température ambiante puis frittés à 1600°C pendant 1 heure, et à 1700°C pendant 2 heures dans une atmosphère 95%Ar - 5%H₂,

- des agglomérats de la même poudre tamisés-forcés à 50 µm et frittés dans les mêmes conditions.

## μ-XANES

Le μ-XANES a permis de vérifier les degrés d'oxydation par comparaison des spectres d'absorption des échantillons aux étalons UO₂, U₃O₇ et U₃O₈. C'est aussi un outil pour déterminer si le contraste est suffisamment élevé pour mener à bien une μ-tomographie X haute résolution.
<figure><img src="image_50.png" /><figcaption>FIGURE 1.1 – Description de la ligne ID22</figcaption></figure>

<figure><img src="image_51.png" /><figcaption>FIGURE 1.2 – échantillon 1</figcaption></figure>

Les spectres de µ-XANES ont été enregistrés entre 17.100 et 17.350 keV avec une résolution de 2eV. Le seuil d'absorption $L_{III}$ de l'échantillon cru de "UO$_2$" (UO$_2$ oxydé à l'air à température ambiante) est constitué de deux sous-spectres comme le laisse voir la dissymétrie du pic de la raie blanche de la figure 1.3 (deux maxima sont observés : $E_{A1}$ = 17.179keV et $E_{A2}$ = 17.182keV). Cet écart de 3eV correspond au décalage entre l'abosrption de $U(+IV)$ et celle de $U(+V)$.

Ainsi, une µ-tomographie subséquente menée à différentes énergies adaptées, après reconstruc-
tion et soustraction entre-elles devrait permettre de visualiser en 3D les degrés d'oxydation de U.

De plus, l'acquisition des spectres µ-XANES en transmission permet de vérifier que l'in-
tensité transmise présente suffisamment de contraste pour mener à bien des expériences de µ-
tomographie.

Il faut, cependant, tenir compte du fait que plusieurs phénomènes peuvent fausser les résultats en entraînant une distorsion des spectres XAS :

- un mouvement de faisceau,

- un changement de taille de faisceau.
<figure><img src="image_52.png" /><figcaption>FIGURE 1.3 – XANES en transmission d’agglomérats cru et fritté de UO₂</figcaption></figure>

µ-tomographie X haute résolution

La µ-tomographie a été testée pour étudier la possibilité de mesurer la densité d’échantillons de UO₂ et d’en suivre l’évolution lors du frittage.

Dans ces expériences l'énergie utilisée (22keV) est supérieure à celle du seuil $L_{III}$ de U occultant les divers degrés d'oxydation de U. Les images sont obtenues en transmission en faisant tourner l'échantillon sur 180° par pas de 0.288° et enregistrées à l'aide d'une caméra FRELON2000. L'image, de 2048*2048 pixels avait une résolution limitée à 0.7µm. Ces images sont corrigées du bruit de fond et du "flat-field" pour conduire à des reconstructions 3D.

Les reconstructions 3D sont elles-même découpées en tranches 2D parallèles d'épaisseur 0.7µm afin de déterminer, à l'intérieur de l'agglomérat, dans chacune de ces tranches, le niveau de gris (densité). Ce dernier est ensuite tracé en fonction de la distance à l'interface. Le passage du
<figure><img src="image_53.png" /><figcaption>FIGURE 1.4 – Reconstruction 3D des agglomérats cru et fritté (diamètres respectifs de 70µm et 50µm)</figcaption></figure>

niveau de gris à la densité pour un agglomérat a été fait par des mesures de porosité mercure sur les agglomérats au laboratoire Bernard François du CEA Cadarache. Les résultats sont donnés figure 1.5 : on voit que pour les agglomérats naturels la densité est plus élevée sur la couche la plus extérieure du grain du fait du procédé de fabrication (atomisation-séchage). Cette densité élevée est observée tout au long du frittage. Pour les échantillons tamisés-forçés la densité élevée à la surface des agglomérats est probablement liée au fait que le tamisage-forçage tasse plus les agglomérats en surface qu'au cœur. Là encore cet effet ne disparaît pas au cours du frittage. Il faut toutefois noter qu'il n'y a pas eu de corrections de faites sur les niveaux de gris et qu'on les a supposés reliés à l'absorption seule de $UO_2$, donc à la densité.

<figure><img src="image_54.png" /><figcaption>FIGURE 1.5 – Densité intra-agglomérat en fonction de la distance à la surface de l'agglomérat</figcaption></figure>

Un échantillon de $(U, Pu)O_2$ avait pu être transporté à l'ESRF (ce qui avait fait l'objet préalable d'une Commission de Sûreté des Installations) mais sa double enveloppe (ayant pour fonction une double barrière contre la contamination) ne permettait pas de bien le localiser sur le porte-échantillons. Étant mal centré (sur le côté de la tige d'acier au lieu d'être au dessus) la reconstruction 3D n'a pas pu se faire, on a pu seulement obtenir quelques images en 2D en
transmission entre le seuil de U et de Pu et au dessus du seuil de Pu. Elle laisse, cependant,
distinguer les agrégats de UO₂ d'un agrégat de PuO₂.

<figure><img src="image_55.png" /><figcaption>FIGURE 1.6 – Image en transmission d'un échantillon entre le seuil de U et celui de Pu (E=17,5keV - image du haut) et au dessus du seuil de Pu (E=18,5keV - image du bas)</figcaption></figure>

**1.2 MEB environnemental**

Dans le cadre d'une collaboration entre le Département d'Études du Combustible (DEC)
et les Laboratoires de Microscopie Électronique et Environnementale (LM2E) et des Interfaces
de Matériaux en Évolution (LIME) UMR 5257 - R. Podor, N. Clavier) de l'Institut de Chimie
Séparative de Marcoule, G.I. Nkou Bouala a étudié dans son stage de master I puis de Master II
de Chimie Séparative, Matériaux et Procédés (CSMP), intitulé "Etude du frittage de systèmes
simplifiés par microscopie électronique à balayage à haute température " l'évolution d'agglomé-
rats denses de cristallites de CeO₂ et de Nd₂O₃ en température au microscope environnemental
[GINB2012].

La morphologie initiale des agglomérats denses de CeO₂, i.e. avant frittage, consiste en des
sphères de 200 ou de 600 nm de diamètre en fonction du protocole d'obtention, tandis que celle
des agglomérats denses de Nd₂O₃ était soit cylindrique (3 µm de longueur, 1.5 µm de diamètre),
soit en prismes hexagonaux (6 µm d'arête), soit cubique polydisperse ( arêtes de 2 à 9 µm).
Un agglomérat dense de CeO₂ isolée a été suivie à différentes températures (900, 1100 et 1200°C)
et l'analyse d'image (électrons secondaires) a permis de suivre l'évolution du nombre de grains
dans cet agglomérat au cours du temps à une température fixée. À 1100°C et 1200°C le nombre
de grains dans l'agglomérat passe de 250 à 1 en 210 (resp. 25) minutes. Les grains finaux gardent
une forme sphérique.

De tels grains finaux de même diamètre ou de diamètre différent pourraient servir de matière première pour une nouvelle expérience de frittage en vue d'une validation du modèle à l'échelle sub-granulaire.

Dans ces expériences la taille de la zone de contact entre "particules" (i.e. agglomérats denses) a été suivie, ainsi que l'angle dièdre apparent jusqu'à l'obtention d'un grain unique pour chacun de ses agglomérats denses, configuration dans laquelle les grandeurs suivies correspondent effectivement respectivement au cou et à l'angle dièdre entre deux grains.

1.3 Autres pistes

L'étude bibliographique montre que récemment le frittage de systèmes hétérogènes de poudres
a été étudié par Olmos-Navarette [LRO2010] pour une poudre de cuivre atomisée avec une large
distribution de taille de particules comme matériau modèle, avec des défauts modèles (porosités
artificielles comme défauts d'empilement et particules d'alumine comme inclusions inertes). Le
frittage a pu être réalisé in-situ sous le faisceau de rayons X à l'ESRF et être suivi en 3D. Cependant, certaines caractéristiques microstructurales telle la taille de cous n'a pas pu être mesurée
faute d'une résolution suffisante. La distribution de taille de pores a toutefois pu être mesurée.

Parmi les autres pistes à explorer, une autre possibilité, toujours en μ-tomographie, au lieu d'observer en transmission, i.e. d'obtenir une information sur la porosité, serait de travailler en contraste de phase afin d'obtenir une reconstruction 3D des joints de grains.
# Moyens de validation à l'échelle des milieux continus

## Sommaire

<table><tr><td>2.1</td><td>Interface grain/gaz - aire et courbure moyenne</td><td>113</td></tr><tr><td>2.1.1</td><td>Les microscopies optique et électronique à balayage</td><td>113</td></tr><tr><td>2.1.2</td><td>La diffusion de neutrons aux petits angles</td><td>114</td></tr><tr><td>2.2</td><td>Texture 3D</td><td>118</td></tr><tr><td>2.3</td><td>Autres techniques</td><td>120</td></tr></table>

## 2.1 Interface grain/gaz - aire et courbure moyenne

Les modèles de frittage, qui considèrent la pastille comme un milieu continu, reposent sur la notion de contrainte de frittage qui est fonction de l'aire de l'interface solide/gaz par unité de volume ainsi que de la courbure moyenne de ces mêmes interfaces.

Dans ce cadre, trois techniques ont été utilisées pour étudier l'évolution de la taille des pores à différents stades du frittage de pastilles de combustible $(U, Pu)O_2$ dopé au chrome :

- la microscopie optique,

- la microscopie électronique à balayage,

- la diffusion de neutrons aux petits angles.

### 2.1.1 Les microscopies optique et électronique à balayage

La microscopie optique est utilisée au grossissement "x2" enregistrée sur 12 images dont les pixels font $2.8\mu m$ d'arête, "x10" sur 300 images. Cette technique couvre la gamme d'observation de $30\ \mu m$ à $1\ cm$.

La microscopie électronique à balayage fournit, au grossissement "x500", une série de 10 champs de $256x256$ pixels couvrant des dimensions linéaires de 1 à 20 ou $30\ \mu m$ et, au grossissement "x2000", 25 champs de $128x128$ pixels correspondant à une gamme d'échelle d'observation allant de 0.2 à $1\ \mu m$.

On ne descend pas avec cette dernière technique en dessous de 200 $nm$. Par ailleurs, ces distributions de taille 2D nécessitent une correction afin de remonter à une distribution 3D. La littérature [CCT1983] fournit des démarches numériques pour y arriver.

La figure 2.1 montre, sur un exemple pour lequel on a effectué cette démarche $((U, Pu)O_2)$ frittée
à 1700°C pendant 4 heures), la densité de population obtenue en fonction du diamètre des pores. Chacun des estimateurs de densité de population est calculé à partir des échantillons par la méthode des noyaux (ici des noyaux gaussiens ont été choisis). On a représenté en gras la courbe résultante qui utilise, dans chaque gamme de diamètre, les données issues de la technique la plus adaptée.

<figure><img src="image_56.png" /><figcaption>FIGURE 2.1 – Densité de population de pores en fonction de leur diamètre</figcaption></figure>

## 2.1.2 La diffusion de neutrons aux petits angles

Les matériaux à base d'actinides étant peu perméables aux rayons X du fait de leur nombre de charge Z élevé; les X traversent tout au plus 200 µm de matériau, seuls les neutrons peuvent traverser un objet centimétrique du type d'une pastille dans le cadre d'une analyse non destructive.

La diffusion aux petits angles permet de déterminer la forme et l'organisation de particules, dans notre cas des pores, dans un milieu homogène, dans notre cas la matrice (U, Pu)O₂ [JPC1999]. Pour mesurer des tailles de 0,5 à 50nm avec des neutrons de longueur d'onde comprise entre 0,1 et 0,2nm il faut que l'angle $θ$ du vecteur de diffusion $\vec{q}$ avec $q = \frac{4\pi}{\lambda} \sin(\theta/2)$ soit suffisamment petit. Un faisceau de neutrons est envoyé sur l'échantillon et on étudie les variations de l'intensité diffusée en fonction de $θ$.

Par ailleurs, on va voir que par diffusion de neutrons aux petits angles, on peut, lorsque leur longueur d'onde est négligeable par rapport aux rayons de courbure des pores ($ρ \gg 100nm$), obtenir l'aire totale de l'interface solide/gaz, tandis que lorsque ces pores sont nanométriques, cette technique donne aussi la distribution de taille des pores.

Cette technique a déjà été utilisée dans les études de matériaux nucléaires ([ADH1983]) pour mesurer la distribution de la phase secondaire MgO présentant, comme dans le cas des pores,
un contraste aux neutrons par rapport à la matrice $UO_2$.

Les essais que nous avons réalisés sur $UO_2$ et $(U, Pu)O_2$ en vue de l'étude de leur porosité ont été réalisés au Laboratoire Léon Brillouin (LLB) sur la ligne de neutrons PAXE (5G-4) du réacteur Orphée au CEA de Saclay.

Le dispositif expérimental est celui représenté en figure 2.2. La grandeur intéressante, caracté-

<figure><img src="image_57.png" /><figcaption>FIGURE 2.2 – schéma de principe du dispositif DNPA sur la ligne PAXE au LLB (http://www-llb.cea.fr/fr-en/pdf/paxe-llb.pdf)</figcaption></figure>

ristique du matériau étudié est la section efficace différentielle (par unité d'angle solide) et par unité de volume, notée : $\sigma(q)$ où $\vec{q}$ est le vecteur de diffusion des neutrons.

$$ \sigma(q) = \left. \frac{d\Sigma}{d\Omega} \right|_q \qquad (2.1) $$

La grandeur directement mesurée par le détecteur XY est l'intensité diffusée $I_N(q)$ dans l'angle solide $\Delta\Omega$ d'une cellule de détecteur. Elles sont reliées par la relation :

$$ I_N(q) = \Phi \epsilon A \Delta\Omega T e \sigma(q) \qquad (2.2) $$

dans laquelle :

- $\Phi$ est le flux incident de neutrons en $cm^2s^{-1}$,

- $\epsilon$ est l'efficacité du détecteur (sans unité),

- $A$ l'aire de la surface d'une cellule du détecteur en $cm^2$,

- $\Delta\Omega$ l'angle solide sous lequel de l'échantillon (considéré comme ponctuel) on voit une cellule du détecteur,

- $T$ la transmission de l'échantillon (sans unité),

- $e$ est l'épaisseur de l'échantillon en $cm$,
- $\sigma$ est en $cm^{-1}$.

Pour un ensemble d'atomes différents, chaque élément $\alpha$ comportant $n_\alpha$ atomes de longueur de diffusion cohérente $a_\alpha$ et de section efficace de diffusion incohérente $\sigma_{inc}^\alpha$ on a [JPC1999] :

$$ \sigma(q) = \sum_{\alpha} n_{\alpha} \frac{\sigma_{inc}}{4\pi} + \sum_{\alpha,\beta} \langle a_{\alpha} \rangle \langle a_{\beta} \rangle \sum_{i,j}^{n_{\alpha},n_{\beta}} \left\langle e^{i\vec{q} \cdot (\vec{r}_i^{\alpha} - \vec{r}_j^{\beta})} \right\rangle \quad (2.3) $$

où $i$ et $j$ désignent les sites cristallographiques des atomes qui diffusent. Le premier terme qui est une diffusion incohérente (pas de cohérence entre la position et la valeur du spin nucléaire) est constante et ne varie pas avec $q$. Elle n'apporte pas d'information en DNPA, on l'enlève du signal.

La seconde partie n'est pas utilisée sur les atomes eux-mêmes mais sur des groupes d'atomes appelés "diffuseurs". Il est de la forme :

$$ \sigma_{coh}(q) = \frac{1}{V} (\rho_{\alpha} - \rho_{\beta})^2 \left| \int_V e^{i\vec{q} \cdot (\vec{r}^{\alpha} - \vec{r}^{\beta})} dV \right|^2 \quad (2.4) $$

Cette dernière intégrale dépend de l'organisation du solide qui est aussi reliée à ses propriétés thermodynamiques. Elle s'écrit en termes de facteurs de structures partiels $S$ et de facteurs de forme $P$ des particules. Pour un milieu biphasé avec un soluté $\alpha$ et un solvant $l$, $\rho_\alpha = a_\alpha - a_l \frac{v_\alpha}{v_l}$ où $v$ représente le volume.

Au final, pour un milieu bi-phasé on a :

$$ \sigma_{coh}(q) = (\Delta\rho)^2 N_p V_p S(q) P(q) \quad (2.5) $$

$\Delta\rho$ est appelé contraste. On calcule ici le contraste entre $(U_{0,89}Pu_{0,11})O_2$ et l'air. $\bar{b}_{matrice} = 0,89\bar{b}_{UO_2} + 0,11\bar{b}_{PuO_2}$ avec $\bar{b}_{UO_2} = \frac{b_{238U}+2b_O}{3}$ et $b_{238U} = 0,841710^{-12}cm$ et $b_O = 0,5810^{-12}cm$ avec 74% de $^{239}Pu$ et 19% de $^{240}Pu$ on a $b_{Pu} \approx 0,74 \cdot 0,7710^{-12}+0,19 \cdot 0,3510^{-12} \approx 0,68610^{-12}cm$. Par ailleurs, $V_{at} \approx 13,578Å^3$. On arrive à un contraste $(\Delta\rho)^2 = 2.0510^{-3}$

Il reste à traiter les termes $S(q)$ et $P(q)$. Ils dépendent des approximations les plus adaptées au système que l'on étudie quant à la forme et la taille des inclusions (vis-à-vis de la longueur d'onde des neutrons).

Si la taille des pores est petite devant la longueur d'onde des neutrons alors, on est dans le cas de l'approximation de Guinier :

$$ I(q) = I(0)e^{-\frac{1}{3}\xi^2 q^2} \quad (2.6) $$

avec $\xi$ le rayon de l'inclusion, appelé rayon de giration, soit une loi en logarithmes de type : $\ln(I(q)) = \ln(I(0) - \frac{q^2\xi^2}{3})$.

Pour les pores de grande taille devant la longueur d'onde, l'allure de $\ln(I)$ en fonction $\ln(q)$ est une droite décroissante de pente $-4$ (loi de Porod). Dans cette plage, les neutrons voient les interfaces grain/gaz comme s'il s'agissait de plans. L'approximation suivante est alors valable :

$$ I \approx 2\pi (\Delta\rho)^2 \phi (1 - \phi) \frac{4}{\xi} \frac{1}{q^4} \quad (2.7) $$