Monographie Materiaux Du Nucleaire Outils Modelisation

Les outils de la modélisation
L’approche multi-échelles
de la modélisation des matériaux
L a simulation multi-échelles des matériaux est une disci- sibles directement par l’expérience. D’autre part, le coût et la
pline émergente qui vise à prédire le comportement d’un difficulté de réaliser des expériences sur les matériaux irra-
matériau en couplant des techniques de modélisation opérant diés. Et, enfin, la nécessité de prédire le comportement des
à différentes échelles d’espace et de temps. Ces échelles vont matériaux avant de les tester et de transposer les résultats
du nanomètre (quelques distances interatomiques) au milli- obtenus en vieillissement accéléré par irradiation aux ions,
mètre ou au mètre (échelle des composants), et de la picose- notamment.
conde (période de vibration d’un atome dans un solide) à
quelques années (durée de fonctionnement d’un composant) La très grande majorité des simulations multi-échelles consiste
ou quelques siècles (pour les matériaux de stockage). Cette à coupler deux échelles de simulations, ou plus, mais unique-
approche est liée au caractère intrinsèquement multi-échelles ment par passage d’information d’une échelle à l’autre, non
des phénomènes qui régissent l’évolution des propriétés des en couplant différentes échelles dans une même simulation.
matériaux : c’est l’énergétique des défauts à l’échelle atomique Ces échelles peuvent être grossièrement subdivisées en
qui pilote la cinétique d’évolution de la microstructure du maté- deux : les échelles les plus fines, celles de la physico-chimie
riau et des propriétés macroscopiques qui en découlent. des matériaux et les plus grandes, celles de la mécanique des
matériaux. Ces deux catégories, qui se rejoignent typiquement
Le nucléaire fait historiquement partie des secteurs en pointe à l’échelle de la dynamique des dislocations, obéissent à une
dans le développement et l’utilisation des différents outils de logique un peu différente et sont présentées successivement
simulation des matériaux, pour plusieurs raisons. D’une part, ci-après (fig. 17).
la complexité : la matière est fortement perturbée par l’irradia-
tion neutronique, qui créé des défauts spécifiques, notamment La modélisation des matériaux sous irradiation a toujours été
les interstitiels, dont les propriétés sont difficilement acces- en pointe : une des premières applications de la dynamique
moléculaire a été l’étude de la forma-
tion de défaut par irradiation dans un
cristal [1].
Temps
50 nm
1 an
Les outils de la
modélisation physique
des matériaux
pour le nucléaire
10 µm
Les physiciens ont tendance à partir de
Monte-Carlo sur Objet
Monte-Carlo cinétique (OkMC) l’échelle la plus fine, celle à laquelle la
1s Atomique (AkMC) Dynamique d’amas
physique est la plus robuste. C’est
l’échelle des calculs de structure élec-
tronique. Les calculs dits « ab initio* »
basés sur la Théorie de la
Fonctionnelle de la Densité (DFT)*
ont connu un formidable essor ces der-
1 ps nières années (voir infra, pp. 29-32,
40 nm
« Les calculs électroniques ab initio
Dynamique moléculaire
classique (DM) pour la matière condensée »). Ils per-
Structure électronique
Ab initio (DFT) Distance mettent de prédire, sans paramètres
1 nm 1 µm 1 mm ajustables, la plupart des propriétés
dans la majorité des matériaux avec
Fig. 17. Représentation schématique des principaux outils de la modélisation physique
une précision de quelques pourcents,
des matériaux. mais aussi de nombreuses restrictions
Les matériaux du nucléaire 25

Modélisation et simulation des matériaux de structure
liées, en particulier, à la taille du système simulé (au plus angulaires (comme les métaux des groupes V et VI), et l’ajus-
quelques centaines d’atomes). Les approximations les plus tement des paramètres pour les alliages reste un travail déli-
courantes (LDA ou GGA) permettent en particulier de prédire cat. Combinés à la Dynamique moléculaire*, les potentiels
les énergies de formation ou de migration de défauts dans les empiriques permettent de simuler la trajectoire de systèmes
métaux avec une précision de l’ordre du dixième d’eV. Pour contenant jusqu’à plusieurs millions d’atomes et sur des temps
les non métaux (semi-conducteurs, isolants) et pour les sys- allant jusqu’à la microseconde. Les applications les plus cou-
tèmes à électrons fortement corrélés comme les actinides il rantes dans le domaine des matériaux pour le nucléaire sont
faut aller au-delà de ces approximations, et utiliser des les cascades d’irradiation et les interactions entres disloca-
méthodes, encore en cours de développement, qui nécessi- tions et défauts d’irradiation. Il est également possible d’effec-
tent des temps de calculs plus longs jusqu’à plusieurs ordres tuer des simulations de dynamique moléculaire quantique, en
de grandeurs (GW, RPA, LDA+U, DFMT). Les calculs ab ini- utilisant les énergies et les forces issus de calculs ab initio.
tio dans les matériaux ont fortement bénéficié de l’augmenta- Ces simulations, très coûteuses en temps de calcul, devraient
tion de la puissance de calcul de ces dernières années. Ils se développer de plus en plus, en particulier pour accéder aux
peuvent être avantageusement complétés par des calculs de énergies libres de formation ou de migration de défauts.
structure électronique semi-empiriques (de type liaisons
fortes), moins lourds numériquement, mais qui nécessitent Il existe ensuite une série de modèles cinétiques*, permet-
des paramètres ajustables (voir infra, pp. 29-32, « Les calculs tant de simuler l’évolution de systèmes sous irradiation ou
électroniques ab initio pour la matière condensée »). sous sollicitation thermique sur des temps beaucoup plus
longs qu’en dynamique moléculaire (voir infra, pp. 57-60, « Les
Les potentiels interatomiques* permettent une description modèles cinétiques »). Dans les simulations de Monte-Carlo
moins quantitative des interactions entre atomes, mais beau- cinétique atomique, il n’est pas tenu compte du détail des
coup moins gourmande en temps de calcul (voir infra, pp. 33- vibrations des atomes autour de leur site, seuls les change-
37, « Les potentiels interatomiques »). Ils sont de plus en plus ments de site sont pris en considération. Les données d’entrée
souvent ajustés sur des bases de données obtenues par cal- de ces modèles, les fréquences de sauts, sont maintenant
cul ab initio. Dans les métaux, différents modèles permettent obtenues le plus souvent à partir de calculs ab initio. Pour
de bien rendre compte du caractère non additif des liaisons. accéder à des temps encore plus longs, les méthodes de
C’est le cas par exemple des très populaires potentiels EAM Monte-Carlo cinétique mésoscopique considèrent directement
(Embedded Atom Method). Un certain nombre de métaux les macro-sauts des objets étudiés (Monte-Carlo sur objet) ou
purs sont relativement bien décrits par ces potentiels. Des limi- directement les probabilités de rencontre entre objets (Monte-
tations sont rencontrées dans les métaux avec de fortes forces Carlo sur événements). Pour accéder aux temps encore plus
+ = rien
+
+ =
+
+ =
10 ps ; 10-100 nm3 103 s ; 10-100 µm3

Dommage primaire Évolution de la microstructure Effets sur les propriétés du matériau
Simulations Modèles cinétiques Propriétés mécaniques :

de Dynamique moléculaire Monte-Carlo, Dynamique des dislocations
de cascades d’irradiation dynamique d’amas
Propriétés des dislocations

Propriétés des défauts (plans de glissement, mobilité)
Calculs ab initio Dynamique moléculaire
Fig. 18. Modélisation multi-échelles des matériaux sous irradiation ; les carrés représentent les lacunes, les haltères représentent
des « dumbbells », paires d’interstitiels.
26 L’approche multi-échelles de la modélisation des matériaux

longs, on ne considère plus que les concentrations des diffé- ralement assuré par la simplification et l’homogénéisation des
rentes populations d’objets, en négligeant toutes les corréla- mécanismes décrits.
tions spatiales. Ces méthodes sont complémentaires.
Par exemple, une dislocation est représentée en dynamique
Cette liste de méthodes n’est pas exhaustive. L’évolution de la moléculaire par un défaut de l’empilement des atomes par rap-
microstructure peut par exemple être également décrite par port au cristal parfait, s’étendant linéairement sur quelques
des méthodes de type champ de phase*. Ces méthodes sont centaines de distances interatomiques. Son mouvement sous
actuellement encore peu développées pour les matériaux contrainte est régi par un potentiel interatomique tel que décrit
sous irradiation. plus haut. Typiquement, il est possible de simuler, dans ces
conditions, le mouvement d’une dislocation pour des temps
La manière dont ces différentes méthodes sont utilisées de de l’ordre de la picoseconde. Afin d’accroître la durée et le
manière séquentielle (sequential multiscale modeling) pour volume simulés, cette même dislocation est représentée par
simuler le comportement d’un matériau sous irradiation est une ligne porteuse d’un champ de déformation et dotée d’une
illustrée sur la figure 18. Le dommage primaire, c’est-à-dire la loi de mobilité pouvant être déduite des calculs de dynamique
distribution spatiale des différents types de défauts après la moléculaire. Cette simplification est à la base des codes de
cascade, est obtenue par dynamique moléculaire classique. dynamique des dislocations qui permettent l’étude du com-
Les potentiels interatomiques utilisés sont ajustés sur les cal- portement collectif des dislocations décrivant la formation de
culs ab initio. Les modèles cinétiques, de type Monte-Carlo ou structure de dislocation jusqu’à l’échelle de quelques microns
Dynamique d’amas*, permettent de simuler comment la et pour des temps de l’ordre de la centaine de microsecondes
migration de ces défauts fait évoluer la microstructure du (fig. 18).
matériau. Les données d’entrée de ces modèles (les énergies
de migration et de liaisons des défauts) sont également obte- Pour des temps plus longs, la multiplication des dislocations,
nues à partir de calculs ab initio. Enfin, l’effet de cette micro- sous l’effet d’une sollicitation extérieure, envahit peu à peu le
structure sur les propriétés du matériau est simulé par la cristal et le mouvement de celles-ci se traduit par des déforma-
Dynamique des Dislocations* (voir infra, pp. 63-65, « La tions à l’échelle du grain, généralement incompatibles avec
dynamique des dislocations »), les phénomènes élémentaires celles des grains voisins et créant, par conséquent, des
à prendre en compte étant eux-mêmes tirés de simulations en contraintes internes au polycristal. À cette échelle, il n’est plus
dynamique moléculaire. possible de considérer les dislocations individuellement ; il
peut y avoir, en effet, plusieurs kilomètres de ligne de disloca-
Les cas où plusieurs méthodes de simulations sont utilisées tion par millimètre cube de matériau. Les dislocations sont
simultanément, au cours de la même simulation (concurrent alors représentées par leurs densités dans les modèles de
multiscale modeling), sont en pratique assez rares. Nous pou- plasticité cristalline. Par exemple, on peut admettre que la
vons signaler, par exemple, le cas des fractures : les quelques vitesse de déformation d’un volume cristallin sera proportion-
atomes en tête de fissure sont simulés par des calculs ab ini- nelle à la densité des dislocations mobiles qu’il contient et à la
tio. Les atomes qui les entourent sont décrits par des poten- vitesse de déplacement de celles-ci (loi d’OROWAN). Dans ce
tiels empiriques. La zone la plus éloignée est décrite à l’échelle cadre, il est possible d’accéder à des temps longs et des
continue, par éléments finis. volumes de matière dont le comportement est représentatif
de l’état macroscopique. Cela n’est, toutefois, possible que si
on dispose d’une description de la microstructure, c’est-à-dire,
Les outils de la modélisation de la disposition spatiale relative des grains et de leurs orien-
mécanique des matériaux tation cristallographiques. La représentation numérique de la
pour le nucléaire microstructure et les lois de plasticité cristalline permettent,
alors, via des outils numériques comme les éléments finis ou
La modélisation mécanique des matériaux pour le nucléaire
les transformations de Fourier rapides (sigle anglais FFT ou
a pour objectif de décrire les propriétés mécaniques et surtout
Fast Fourier Transform) de reproduire le comportement macro-
leur évolution en condition de fonctionnement nominale ou
scopique. La figure 19, page suivante, illustre le principe de
incidentelle. Tout comme la modélisation des effets de l’irra-
cette démarche.
diation évoquée précédemment, celle-ci se fonde sur une
approche multi-échelles. Les outils à l’échelle atomique sont
Comme nous le verrons plus loin, dans cette section, l’endom-
identiques : calculs ab initio et dynamique moléculaire. Ces
magement et la rupture des matériaux fait l’objet d’une
simulations extrêmement précises de la physique présentent
démarche similaire, bien que simplifiée. En effet, l’introduction
l’inconvénient de ne pouvoir considérer que des volumes de
explicite de discontinuités dans le milieu simulé complexifie et
matière et des temps petits devant ceux visés par la modéli-
alourdit de façon significative les calculs ; il y a, dans ce
sation du comportement mécanique. Il faut donc disposer
domaine, de nombreux progrès à réaliser.
d’outils qui assurent simultanément le changement d’échelle
d’espace et de temps de simulation. Ce changement est géné-

Représentation numérique
Modélisation de la microstructure
de la microstructure
Microstructure
Orientations cristallographiques
Géométrie et dimensions des grains
Comportement
macroscopique
Dynamique
des dislocations Éléments finis ou FFT
Mécanismes Comportement du monocristal

Règles de mouvement des dislocations en plasticité cristalline
Obstacles pour les dislocations
(autres dislocations, cavités, précipités…)
Fig. 19. Principaux outils de la modélisation de la mécanique des matériaux. Le lien avec les échelles inférieures se fait, principalement,
au travers des mécanismes considérés par la dynamique des dislocations.
François WILLAIME et Bernard MARINI,

Département des matériaux pour le nucléaire
Références
[1] J.B. GIBSON, A.N. GOLAND, G.H. VINEYARD, Phys. Rev., 120 (1960),
p. 1229.
28 L’approche multi-échelles de la modélisation des matériaux

Les calculs électroniques ab initio

pour la matière condensée
L a physico-chimie des phénomènes mis en jeu dans la Évidemment, cette transformation a un prix. Tous les phéno-
matière condensée provient quasiment exclusivement du mènes quantiques ont été cachés dans une expression rece-
comportement des électrons. Ainsi, tous ces phénomènes lant toute la complexité du problème : l’énergie d’échange-cor-
variés tels que la structure des solides, leurs constantes élas- rélation. L’échange provient de la nature indiscernable des
tiques, leur couleur, leur caractère métallique, semi-conduc- électrons qui sont des fermions, tandis que la corrélation
teur ou isolant, leur teneur en défauts, leur déviation de la stœ- désigne, par définition, tout le reste des effets quantiques. Plus
chiométrie, se comprennent par les propriétés des électrons. précisément, le terme de « corrélation » correspond au fait
Les électrons sont, par nature, des particules quantiques ; leur qu’on ne peut décrire la position d’un électron sans tenir
description nécessite donc la résolution de l’équation fonda- compte de la position des autres. Leurs mouvements respec-
mentale de la mécanique quantique, l’équation de tifs sont dits « corrélés ».
SCHRÖDINGER.
Heureusement, des approximations extrêmement simples,
Les calculs électroniques dits « ab initio* », c’est-à-dire à mais pourtant très précises, ont été proposées dès les travaux
partir des premiers principes, proposent donc d’étudier les pro- pionniers des années 60 [2]. Dans la DFT, l’énergie
priétés des électrons sans ajustement spécifique au système, d’échange-corrélation est, en principe, une fonctionnelle de la
directement à partir de la résolution de l’équation de densité électronique (). Cela signifie que cette énergie est
Schrödinger. Cela représente, à première vue, un défi formi- une fonction dont la valeur dépend de la connaissance de la
dable, puisque l’équation de Schrödinger devient déjà très dif- fonction () en tout point de l’espace. L’approximation de la
ficile à résoudre, dès que le nombre de particules quantiques densité locale (LDA) propose de limiter cette dépendance
interagissant dépasse l’unité. Comment donc procéder alors drastiquement : l’énergie d’échange- corrélation n’est fonction
que le nombre d’électrons d’un solide macroscopique est de que de la valeur de la densité au point considéré :
l’ordre de 1023 ?
! = "()ε !(()).
La théorie de la fonctionnelle ε ! désigne la densité d’énergie d’un système de référence

de la densité dont on connaît l’énergie d’échange-corrélation : le gaz d’élec-
trons homogène qui a pu être calculé numériquement de
Inventée dans les années 60 [1,2], mise en pratique dans des
façon exacte, grâce à l’avènement du Monte-Carlo quantique,
programmes informatiques, au cours des années 80-90, la
au début des années 80 [6]. Dans le cadre de la LDA, l’éner-
théorie de la fonctionnelle de la densité (DFT)* a connu un
gie d’échange-corrélation n’est plus une fonctionnelle, mais
essor énorme, à tel point qu’elle règne aujourd’hui sans par-
une simple fonction. Disposant de cette simple évaluation de
tage sur les calculs électroniques dans la matière condensée.
l’énergie d’échange-corrélation, tous les termes sont connus
et calculables.
La DFT, récompensée par le prix Nobel de chimie attribué à
Walter KOHN, en 1998 [3], est un ensemble de théorèmes
La figure 20, page suivante, illustre ces calculs par une étude
mathématiques qui établissent rigoureusement une méthode
du volume atomique du zirconium en fonction de la pression
simplifiant la mécanique quantique des électrons en interac-
pour différentes phases cristallines. Les calculs DFT ont per-
tion. Cette théorie prouve formellement ce qui était pressenti
mis de prédire l’apparition sous pression d’une phase cubique
depuis les origines de la mécanique quantique [4,5], à savoir
centrée pour ce métal.
que les propriétés d’un système d’électrons découlent de la
seule densité des électrons. En d’autres termes, l’objet cen-
On peut améliorer l’approximation locale en incluant une
tral de la mécanique quantique n’est plus la fonction d’onde
dépendance vis-à-vis, non pas seulement de la densité élec-
tronique () elle-même, mais également vis-à-vis de son gra-
des électrons (fonction incroyablement compliquée, car
dient $(). Cette approximation, appelée « approximation du

dépendant de la position de tous les électrons à la fois), mais
plutôt la densité électronique de l’état fondamental (fonction
gradient généralisé (GGA) », reste bien locale (on dit « semi-
beaucoup plus simple, car ne dépendant que de la position).
locale ») et permet une mise en œuvre numérique aisée. Elle
améliore sensiblement les résultats pour les molécules.

que la LDA ou la GGA, mais elles sont indispensables pour
certaines propriétés liées à la bande interdite.
*
Énergie (Ry)
hcp
bcc
Mise en œuvre numérique : pseudopotentiels
0,15 et fonctions de base
La théorie étant posée, il s’agit maintenant de comprendre les
Énergie (mRy)
10
grandes lignes de la mise en œuvre pratique de la DFT. La

0
0,10 *-hcp
DFT, dans sa variante standard KOHN-SHAM [2], requiert le cal-
-10
fcc-hcp
bcc-hcp
cul de fonctions d’onde auxiliaires d’électrons indépendants,
14 16 18 20 22
dont la seule qualité demandée est de donner la bonne den-
Volume (Å3) sité électronique. C’est le calcul et la manipulation de ces fonc-
0,05 tions d’onde auxiliaires qui va constituer le plus grand effort
d’un calcul ab initio.
Les calculs numériques se séparent en deux grandes

0 familles : les calculs « tous électrons » et les calculs à pseu-
dopotentiels*. En effet, la question se pose ainsi : est-il
nécessaire de décrire dans les calculs les électrons proches
des noyaux qui ne participent pas aux liaisons chimiques dans
0 15 20
les solides ? Pour la matière condensée, ce choix penche sou-
Volume (Å3)
vent du côté de la non-prise en compte de ces électrons chi-
Fig. 20. Énergie des différentes phases cristallines du zirconium pur
miquement inactifs, à moins que l’on ne soit intéressé spéci-
en fonction du volume atomique. La stabilité sous pression de la fiquement par ces électrons-là. La technique dite « des
phase cubique centrée (bcc) a d’abord été prédite par les calculs pseudopotentiels » (et ses améliorations ultérieures telles que
avant d’être observée expérimentalement [Physical Review, B 48, la technique Projector Augmented Wave) proposent ainsi de
16269 (1990)].
considérer ensemble le noyau positif et les électrons de cœur
qui en sont proches. Cela permet de traiter moins d’électrons
explicitement et de faire des choix numériques plus adéquats.
Cependant, l’amélioration de la précision est anecdotique pour
la matière condensée. Afin de calculer les fonctions d’onde des électrons, il faut
résoudre des équations différentielles similaires à l’équation
Bien qu’en général très robustes, ces deux approximations, de Schrödinger pour un électron dans un potentiel effectif.
LDA et GGA, échouent néanmoins pour certaines propriétés, Dans l’immense majorité des codes, les fonctions d’ondes
en particulier sur ce qu’on appelle le « problème de la bande Φ&() sont recherchées comme combinaison de fonctions de
interdite ». La bande interdite, lorsqu’elle existe, est la zone base (au sens mathématique) φ&() :
d’énergie qui sépare les états électroniques occupés des états
vides. La connaissance de cette bande interdite est primor- Φ&() = ∑ !)&φ)(),
)
diale, puisque de nombreuses propriétés physiques en décou-
lent. Pour un cristal, son caractère métallique (pas de bande où l’indice & désigne l’ensemble des nombres quantiques per-
interdite), semi-conducteur (bande interdite de l’ordre de l’élec- tinents pour le système considéré : en général, pour les cris-
tron-volt), ou isolant (bande interdite plus large) en est une taux, ce sont un indice de bande (s,p,d,f), un vecteur d’onde
conséquence directe !… Malheureusement, la LDA et la GGA et éventuellement un indice de spin.
échouent lamentablement pour cette propriété essentielle :
toutes deux sous-estiment systématiquement la largeur de la Les fonctions de base les plus communément utilisées pour
bande interdite d’un facteur 2 ou 3. les solides sont illustrées dans la figure 21. D’une part, les
ondes planes, qui sont des oscillations délocalisées dans tout
Aller au-delà en ajoutant une dépendance en termes de déri- le système calculé, offrent un gage de flexibilité et de robus-
vées d’ordre supérieur de la densité électronique complique tesse, puisqu’elles ne dépendent ni des atomes présents ni
énormément la théorie et n’a pas apporté grand-chose jus- de leur position. L’amélioration de la base est systématique :
qu’à maintenant. De nos jours, certaines approches promet- en ajoutant des ondes planes d’oscillations plus rapides, on
teuses, telles que les fonctionnelles hybrides ou l’approxima- permet naturellement une description plus fine des fonctions
tion GW, sortent du cadre strict de la DFT en ajoutant des d’onde Φ&(). Par ailleurs, les ondes planes tirent avantageu-
dépendances autres que la simple densité électronique. Ces sement partie des transformées de Fourier rapides (FFT) qui
approches sont encore peu utilisées pour la matière conden- autorisent de nombreuses astuces de calcul. En revanche,
sée, car leur coût en terme de calcul est infiniment plus grand elles payent le prix de leur simplicité : les ondes planes impo-
30 Les calculs électroniques ab initio pour la matière condensée

sent, par construction, de traiter tout le système avec la même
précision, qui est alors imposée par l’élément du système
induisant les variations les plus rapides.
L’autre grande famille de fonctions de base est constituée par

les bases localisées. La figure 21 représente un exemple de
fonctions localisées. Leur allure précise peut prendre diffé-
rentes formes (gaussienne, orbitale de SLATER, orbitale numé-
rique quelconque, ondelettes...), mais leur trait commun est Ondes planes
que les fonctions de base sont centrées sur les positions des
atomes. Ainsi, les fonctions de base dépendent de la nature
des éléments en présence et de leur position dans l’espace.
Chaque élément requiert un certain travail pour la définition et
la vérification des fonctions de base employées. Mais cela a
l’avantage de traiter différentes parties du système avec une
précision différente de façon flexible. Par exemple, les zones
vides ne sont pas décrites inutilement. Les calculs utilisant une Orbitales localisées
base localisée sont généralement rapides, car il est inutile de
calculer les termes impliquant deux fonctions de base qui ne Fig. 21. Illustration des deux principales familles de fonctions
se recouvrent pas. Les orbitales localisées sont la voie identi-
fiée comme permettant d’arriver à l’« ordre + » : la complexité
de base servant à décrire les fonctions d’ondes : les ondes planes
et les orbitales localisées.
du calcul croîtrait alors linéairement avec la taille du système.
associée à leur légèreté en termes de temps de calcul, les a
rendues incontournables. Bien sûr, certaines propriétés des
Les codes de calcul
électrons sont mal décrites par ces approximations simples,
Les codes informatiques de calcul ab initio pour les solides mais elles restent l’exception plutôt que la règle.
ont connu une transformation énorme depuis la fin des
années 90. La communauté scientifique est passée d’un Les calculs ab initio sont utilisés dans de nombreuses études
modèle fondé sur une multitude de petits codes développés présentées dans cette monographie. Dans ces études, ils sont
au sein de chaque équipe de recherche à un modèle de déve- le plus souvent employés comme un moyen d’obtenir des
loppement de grands codes (libres ou commerciaux) dispo- informations quantitativement correctes sur l’énergie de divers
nibles pour une large communauté d’utilisateurs. Ainsi, il s’est arrangements atomiques. Ainsi, ils sont à la base des modé-
produit une séparation entre les activités de développement lisations thermodynamiques de diagrammes de phases d’al-
de code et les activités applicatives. L’interface des codes se liages sur réseau rigide (voir infra, pp. 49-52, « Les dia-
simplifiant, il n’est désormais plus nécessaire d’être un expert grammes de phase ») qui sont illustrées dans le cas des
dans tous les arcanes et tous les détails numériques de ce alliages fer-chrome (voir infra, pp. 123-126, « La thermodyna-
type de calcul pour obtenir des résultats de simulation ab ini- mique des alliages fer-chrome »). De même, ils permettent de
tio. De l’autre côté de la chaîne, le développement des codes calculer précisément les énergies de formation et de migra-
modernes est supporté par de grandes équipes ou par des tion des défauts ponctuels (voir infra, pp. 39-42, « Le paysage
collaborations entre équipes. énergétique à l’échelle atomique »). Les résultats de ce type
d’étude sont présentés pour plusieurs matériaux : carbure de
Malgré l’accroissement de la complexité des codes ab initio silicium (voir infra, pp. 133-136, « Structure et cinétique des
pour les solides, leur nombre reste néanmoins important. Pour défauts ponctuels dans le carbure de silicium »), métaux de
ne citer qu’un sous-ensemble représentatif de ceux fondés sur transition (voir infra, pp. 137-140, « Structure et cinétique des
la technique des pseudopotentiels (ou ses améliorations), défauts d’irradiation dans le fer »). Ils sont également utilisés
mentionnons ABINIT [7,8], CASTEP, CPMD, Quantum- pour étudier l’évolution microstructurale dans des alliages
Espresso, VASP pour la base des ondes planes, et CP2K, proches des aciers (voir infra, pp. 157-160, « Évolution micro-
Crystal, FHI-AIMS, ONETEP, SIESTA pour les bases locali- structurale dans les alliages modèles d’aciers de cuve ») et
sées. Dans cette liste, la moitié est effectivement utilisée dans pour le cas de l’hélium dans des alliages de fer (voir infra,
les différents services de la DEN. pp. 141-145, « Énergétique et cinétique de l’hélium dans le fer
et les alliages à base de fer »).
En résumé, les calculs ab initio dans les solides sont devenus
quasiment synonymes de calculs DFT. Aujourd’hui, les
approximations les plus simples, telles que la LDA et la GGA, Fabien BRUNEVAL, Jean-Paul CROCOMBETTE
sont mises en œuvre avec une relative facilité dans des codes et François WILLAIME
libres ou commerciaux. La précision de ces approximations, Département des matériaux pour le nucléaire

Références
[1] P. HOHENBERG et W. KOHN, « Inhomogeneous Electron Gas »,
Physical Review, 136 (1964), B864-B871.
[2] W. KOHN et L.J. SHAM, « Self-Consistent Equations Including

Exchange and Correlation Effects », Physical Review, 140 (1965),
A1133-A1138.
[3] W. KOHN, « Nobel Lecture: Electronic structure of matter – wave

functions and density functionals », Review of Modern Physics, 71
(1999), pp. 1253-1266.
[4] L.H. THOMAS, « The calculation of atomic fields », Proc. Cambridge

Phil. Soc., 23 (1927), pp. 542-548.
[5] E. FERMI, « Un Metodo Statistico per la Determinazione di alcune

Prioprietà dell’Atomo », Rend. Accad. Naz. Lincei, 6 (1927), pp. 602-
607.
[6] D. M. CEPERLEY et B. J. ALDER, « Ground State of the Electron Gas

by a Stochastic Method », Physical Review Letters, 45 (1980), pp. 566-
569.
[7] X. GONZE, B. AMADON, P.-M. ANGLADE, J.-M. BEUKEN, F. BOTTIN,

P. BOULANGER et al., « ABINIT : first-principles approach to material and
nanosystem properties », department of chemistry and Institute for
Research in Materials, Dalhousie University, Halifax, Canada
Computer Physics Communications (Impact Factor: 311), 180, 12
(2009), pp. 2582-2615.
[8] M. TORRENT, F. JOLLET, F. BOTTIN, G. ZERAH and X. GONZE,

« Implementation of the Projector Augmented-Wave Method in the
ABINIT code. Application to the study of iron under pressure »,
Comput. Mat. Science, 42 (2008), p. 337.
32 Les calculs électroniques ab initio pour la matière condensée

Les potentiels interatomiques
L es simulations de physique statistique fondées sur une matière condensée »). La contrepartie de la précision des cal-
description atomistique (comme la dynamique moléculaire) culs ab initio est un temps de calcul qui évolue comme au
nous permettent d’accéder à l’évolution de la macrostructure moins le cube N3 du nombre d’atomes considérés N. Les
en partant d’une échelle atomique des matériaux. La robus- limites actuelles des puissances de calcul ne permettent pas
tesse du modèle physique qui décrit la liaison interatomique de traiter plus d’une centaine d’atomes. Cette contrainte
est l’ingrédient clé des simulations et assure leur caractère confine ce type de calcul à l’étude des propriétés physiques de
réaliste, transférable et prédictif. De manière systématique, matériaux massifs ou contenant des défauts de petite taille.
une simplification du modèle physique implique une efficacité Cette contrainte restreint, de plus, les calculs de dynamique
numérique croissante qui va de pair avec une perte importante moléculaire ab initio à l’échelle de la picoseconde.
du caractère prédictif et transférable de la méthode. Cette sim-
plification est l’objet principal des potentiels empiriques. Les Il est légitime d’envisager de repousser les limites spatio-tem-
potentiels empiriques s’attachent à décrire les liaisons entre porelles de ces calculs pour avoir accès à la description de
les atomes à l’aide de fonctionnelles simples, comme celle volumes de l’ordre du nanomètre cube (nm3) pendant des
représentée figure 22. temps de simulation de l’ordre de la nanoseconde. L’extension
spatiale est nécessaire à l’étude de défauts de grande taille,
Les propriétés des matériaux sont pilotées par les électrons et tandis que l’extension temporelle est impérative pour le calcul
leur caractère quantique. Une prédiction exacte de toute pro- de propriétés physiques qui demandent un traitement statis-
priété physique passe par la résolution de l’équation de tique temporel (conductivité thermique, diffusion, énergie libre,
Schrödinger pour tous les électrons du système. Les calculs changements de phase, etc.). L’étude des dommages pri-
de structure électronique ab initio donnent accès aux proprié- maires d’irradiation (voir infra, pp. 127-132, « La modélisation
tés de cohésion des matériaux, quelle que soit la nature des du dommage primaire d’irradiation par dynamique molécu-
liaisons entre les atomes (métallique, covalente ou ionique). laire ») impose à la fois une extension spatiale (quelques 105
Ces calculs ont un grand pouvoir de prédiction car ils décri- atomes) et une extension temporelle (quelques dizaines de
vent explicitement les électrons. Cette description passe par la picosecondes). Notons que ces temps accessibles à ce type
résolution des équations de Schrödinger, via des hypothèses de simulation restent néanmoins courts, comparés à ceux
simplificatrices mais robustes (par exemple, la DFT ; voir atteints par les calculs de Monte-Carlo cinétique sur réseau,
supra, pp. 29-32, « Les calculs électroniques ab initio pour la de champ moyen auto-cohérent ou de dynamique d’amas
(voir infra, pp. 57-60, « Les modèles cinétiques »).
Étendre le domaine spatio-temporel des calculs ab initio

Énergie
Potentiel répulsif requiert une simplification du modèle physique. Une approche

intermédiaire entre les méthodes de type ab initio et celles des
potentiels empiriques s’appuie sur l’observation simplificatrice
que le potentiel effectif ressenti par les électrons est proche
Distance d’équilibre de la superposition des potentiels des atomes neutres.
L’essentiel des modifications électroniques provoquées par les
liaisons chimiques est porté par les électrons des couches
0
électroniques externes ou de valence. Cette observation est à
Potentiel attractif l’origine d’une classe de méthodes illustrée par la méthode
des liaisons fortes. Cette méthode permet d’atteindre des sys-
tèmes de quelques milliers d’atomes. Elle comporte néan-
moins quelques faiblesses, en particulier pour les métaux
Distance interatomique
nobles où il existe une compétition entre les électrons locali-
sés (ceux de la bande d) et délocalisés (ceux de la bande sp).
22. Potentiel d’interaction de type LENNARD-JONES entre deux
atomes. La partie rouge est répulsive, et la verte attractive. La méthode des liaisons fortes constitue néanmoins le point
La distance d’équilibre correspond au minimum d’énergie de départ de la grande majorité des potentiels semi-empi-
du potentiel d’interaction. riques des métaux et des matériaux covalents.

L’objet des potentiels interatomiques – empiriques ou semi- des termes complémentaires à courte distance (entre 1,5 et
empiriques – est de remplacer la description fine de la struc- 12 Å) qui tiennent compte des spécificités iono-covalentes de
ture électronique des liaisons interatomiques par des fonction- chaque liaison. Ces termes complémentaires prennent la
nelles judicieusement choisies pour (a) reproduire ces liaisons forme de potentiels de paire de Buckingham ou de Morse
interatomiques et (b) réduire le temps de calcul. Ces fonction- dans la plupart des oxydes. C’est le cas pour UO2 et l’en-
nelles contiennent un certain nombre de paramètres qu’il semble des composés de structure fluorite (ZrO2, CeO2…),
convient d’ajuster. Bien entendu, la simplification créée par mais aussi pour la zirconolite, les pyrochlores, les spinelles,
l’utilisation de ces fonctionnelles conduit à une perte qualita- les perovskites, etc. Certains oxydes (comme les verres ou les
tive de la description des liaisons. Les potentiels empiriques apatites) contiennent des silicates (SiO4) ou des phosphates
conservent la description des liaisons en fonction de la dis- (PO4). Leurs liaisons sont fortement covalentes et donc direc-
tance entre les atomes, contrairement aux calculs sur réseau tionnelles. La description de la stabilité des angles O-Si-O
de type Monte-Carlo cinétique, champ moyen auto-cohérent nécessite alors de faire appel à des potentiels à trois corps ou
ou dynamique d’amas (voir infra, pp. 57-60, « Les modèles plus qui tiennent compte de la dépendance angulaire.
cinétiques »). En revanche, une des particularités des poten- Finalement, la simulation des dommages primaires (voir infra,
tiels empiriques est de ne pas – ou peu – tenir compte de la pp. 57-60, « Les modèles cinétiques ») par cascade de dépla-
chimie. La conséquence la plus évidente de ce défaut de cements fait intervenir des interactions à très courte distance
construction est, par exemple, l’impossibilité de transférer des (moins de 1 Å). Ces interactions sont décrites par des poten-
potentiels interatomiques dédiés aux métaux à la description tiels de ZIEGLER-BIERSACK-LITTMARK qui sont raccordés aux
d’isolants, et réciproquement. Pour les isolants, les électrons potentiels à courte et longue distances.
sont localisés et conduisent à des liaisons ioniques ou iono-
covalentes, dont la caractéristique principale est la présence Chaque potentiel de paire contient un certain nombre de para-
d’un terme de Coulomb. En revanche, les liaisons métalliques mètres. Le potentiel de Coulomb entre deux ions – par
ne contiennent pas de terme « coulombien », en raison du fort exemple, l’uranium et l’oxygène– dans UO2 – contient deux
écrantage électronique. En général, le transfert d’un potentiel paramètres : les valeurs q de la charge de chaque ion. Les
interatomique d’un composé à l’autre se révèle difficile, même valeurs formelles q = +4 |e| pour l’uranium et q = -2 |e| pour
pour une classe de matériaux donnée. De plus, pour un com- l’oxygène peuvent être considérées comme une première
posé donné, qu’il soit métallique, isolant ou semi-conducteur, hypothèse. Il s’avère que l’analyse des calculs ab initio dément
les potentiels ont un domaine de validité restreint. Enfin, pour cette hypothèse, quelle que soit la méthode d’analyse
un domaine de validité donné, si les propriétés physiques cal- employée. Du point de vue des potentiels empiriques, cela
culées à l’aide de potentiels interatomiques peuvent parfois signifie que les valeurs des charges de l’uranium et de l’oxy-
atteindre une précision de l’ordre du pourcent, par rapport aux gène sont des paramètres. Les valeurs de charges peuvent
valeurs expérimentales ou calculées ab initio, le plus souvent varier, à la condition que la neutralité électronique soit assu-
la précision atteinte est de l’ordre de 10 à 20 %. Les résultats rée, c’est-à-dire que q (uranium) = -2 q (oxygène). Les poten-
déduits des calculs en potentiels empiriques apportent ainsi, tiels à courte distance – de type BUCKINGHAM ou MORSE –
en général, des informations phénoménologiques – c’est-à- contiennent deux ou trois paramètres.
dire, par exemple, l’identification de mécanismes à l’échelle
atomique – et plus rarement des données quantitatives. Les modèles semi-empiriques pour la liaison métallique ont la
caractéristique d’être de courte portée et à plusieurs corps.
Qu’en est-il, par conséquent, de l’intérêt des potentiels empi- L’écrantage des électrons, qui apparaît dans les métaux,
riques, de l’intérêt de les construire et de leur utilisation ? La annule le caractère à longue distance des interactions cou-
réponse est simple : les potentiels empiriques doivent être lombienne. Le terme principal devient un terme d’entourage
construits pour répondre à un objectif. Cela signifie (a) qu’ils à N-corps, complété par un potentiel de paire à courte dis-
sont construits avec un domaine de validité limité comme tance. Ce terme d’entourage ou d’immersion peut se justifier
hypothèse de départ (b) que, dans ce domaine, ils sont soi- à l’aide de la méthode des liaisons fortes ou par analogie avec
gneusement construits et validés de manière indépendante un gaz d’électrons libres. Il représente l’immersion de l’atome
(c’est-à-dire avec des bases de données de construction et dans le gaz des électrons délocalisés de ses voisins et dépend
de validation indépendantes) (c) que les résultats qui en donc de la densité électronique ressentie par chaque atome.
découlent ne sont interprétés que dans ce domaine, et cela Cette classe de potentiels est appelée EAM (pour Embedded
en tenant compte des barres d’erreur. Pour illustrer ces pro- Atom Method ). Ces potentiels sont bien adaptés aux métaux
pos, nous prendrons deux exemples de potentiels empiriques. ayant une structure cubique à faces centrées ou hexagonale
Le premier concerne les solides ioniques et le second les compacte.
alliages métalliques.
En revanche, les métaux de transition de structure cubique
Les interactions dans les solides ioniques contiennent un centrée demandent d’introduire des termes angulaires. Ces
terme principal qui est le terme « coulombien ». Ce terme de termes décrivent le caractère légèrement covalent des liaisons
paires à très grande distance (plus de 12 Å) est modulé par métalliques qui provient de l’occupation partielle des
34 Les potentiels interatomiques

bandes . Les potentiels EAM peuvent être ainsi adaptés en que le choix des fonctionnelles pour décrire le fer, par exemple.
M-EAM (Modified-EAM). Les potentiels BOP (Bond Order Cela conduit à des potentiels ajustés de façon plus précise et
Potential) constituent une version encore plus élaborée et per- transférables qui contiennent des informations difficilement
mettent de tenir compte de manière plus détaillée de la struc- accessibles expérimentalement. Ainsi, l’intégration des forces
ture électronique. Ils ont fait l’objet de nombreux développe- agissant sur les atomes ouvre la voie à la prise en compte non
ments, ces dernières années, et semblent très prometteurs, seulement des configurations de minima ou de points-col,
en particulier pour le tungstène, le fer et le molybdène. mais aussi des configurations hors équilibre (les forces agis-
sant sur les atomes sont alors non nulles).
Comme les potentiels pour les liaisons ioniques, les potentiels
pour les liaisons métalliques contiennent un certain nombre Les paramètres des potentiels empiriques sont ajustés sur la
de paramètres pour chaque terme d’interaction. Pour les première base de données décrite précédemment.
métaux de structure cubique à faces centrées ou hexagonale L’ajustement se fait par minimisation aux moindres carrés de
compacte, par exemple, le terme d’immersion à N-corps la différence entre les valeurs cibles et les valeurs calculées,
dépend du carré de la densité électronique. Celle-ci dépend qui représente une fonction coût. Suivant l’objectif défini, le
de deux paramètres, qui décrivent l’intensité des liaisons en poids relatif de chaque valeur cible est modulé. La minimisa-
fonction de la distance et de l’environnement. Le terme de tion est alors réalisée à l’aide d’algorithmes mathématiques
« paire » contient au moins deux ou trois paramètres. variés. Les méthodes stochastiques issues de la physique sta-
tistique (de type Metropolis, algorithmes génétiques, simplex,
L’ajustement des paramètres des potentiels empiriques – que etc.) sont relativement gourmandes en temps de calcul : la
ce soit pour les potentiels ioniques ou pour les potentiels minimisation d’un potentiel à N-corps pour le fer prend
métalliques – se fait en plusieurs phases. La première étape 24 heures sur les ordinateurs actuels. Les méthodes détermi-
est de définir le domaine de validité des potentiels. Ce nistes qui utilisent le gradient et/ou l’Hessien de la fonction
domaine de validité indique les propriétés physiques que l’on coût (méthodes de Newton, d’Euler, gradient conjugué, etc.)
souhaite décrire au mieux ; elle détermine l’objectif à atteindre. sont plus efficaces : le temps de calcul est réduit à une minute.
Pour les dommages d’irradiation, par exemple, la stabilité rela- L’estimation du gradient de la fonction coût devient alors un
tive des phases susceptibles d’apparaître doit être décrite de problème crucial et très difficile d’un point de vue numérique,
manière fiable – au moins qualitativement. Dans le cas d’UO2, du fait de la forte non-linéarité des fonctions coûts. Cette esti-
il est impératif que la phase amorphe soit moins stable que la mation peut être faite selon les méthodes des différences
phase fluorite (la plus stable à température ambiante). De finies ou, de manière plus efficace, par la méthode de l’adjoint
même, pour les pyrochlores A2B2O7 (A étant un cation triva- du Jacobien de la fonction coût.
lent et B un cation tétravalent) ou pour les spinelles AB2O4,
les phases ordonnées doivent être plus stables que les En pratique, un seul ajustement conduit rarement à une solu-
phases désordonnées (pyrochlore (A1/2B1/2)4O7 ou spinelle tion satisfaisante. Il est, en effet, fréquent que la relaxation des
(A1/3B2/3)3O4), elles-mêmes moins stables que les phases défauts – postérieure à l’ajustement – conduise à des gran-
amorphes [1]. En contrepartie de cette exigence, les proprié- deurs physiques qui s’écartent fortement des valeurs cibles
tés élastiques seront moins bien décrites. dans les alliages métalliques. Il est donc nécessaire de vali-
der les potentiels en recalculant les propriétés et en les com-
Une fois l’objectif déterminé, le moyen d’y parvenir consiste parant à la base de données initiale et à une seconde base
en la construction de deux bases de données physiques. La de données, indépendante de la première. Cette seconde
première contient les données cibles sur lesquelles sont ajus- base de données peut contenir des grandeurs physiques déli-
tés les paramètres. L’approche la plus simple est de choisir cates à intégrer dans l’ajustement. On y trouve fréquemment
des observables expérimentales telles que des propriétés la dilatation thermique – c’est-à-dire l’évolution de la maille en
structurales comme le paramètre de maille et des paramètres fonction de la température ou encore la stabilité relative des
internes, l’énergie de cohésion, la température de DEBYE, les phases ordonnées – désordonnées – amorphes. La validation
propriétés élastiques, etc. Cette approche simple s’est cepen- de l’ajustement conduit donc à un processus itératif de mini-
dant révélée insuffisante à plusieurs reprises, en particulier misation / comparaison en réinjectant les données antérieures
pour la description des défauts dans le fer. En effet, le pay- des configurations relaxées dans la minimisation suivante. Un
sage énergétique des défauts ponctuels est décrit de manière algorithme d’ajustement typique est schématisé dans la
inexacte par rapport aux résultats des calculs ab initio. Cette figure 23a, page suivante.
observation a conduit F. ERCOLESSI et J.B. ADAMS au dévelop-
pement de la méthode « force matching method », à la fin des Les critères de la qualité d’un ajustement, et donc de la qua-
années 90 [2]. Des données plus précises, comme l’énergie lité des potentiels interatomiques, dépendent de l’objectif fixé.
totale et les forces agissant sur les atomes dans différentes Les paramètres structuraux sont, en général, bien décrits, car
configurations, ont été intégrées à la première base de don- on leur attribue un poids important dans la fonction coût. En
nées cibles pour l’ajustement des potentiels. Ces données, au- revanche, la qualification peut apparaître très subjective, et
delà des observables expérimentales, jouent le même rôle clef parfois déroutante, lorsque le potentiel est jugé satisfaisant,

Potentiel i Données expérimentales :
a configuration ab initio
d’équilibre
Configurations ab initio Ajustement

hors équilibre
Potentiel i+1
TESTS Relaxation
POTENTIEL
Potentiel de départ
2
b c Potentiel ajusté
EC15-Epara (eV)
Énergie de formation (eV)
5 DFT
-2
4
-4
-6
0
<110> <111> <100> OCT TED 0 2 4 6 8 10 12 14
Taille de l’amas d’auto-interstitiels
Fig. 23. a) Schéma de principe d’un ajustement de potentiel. b) Comparaison entre les valeurs cibles de la base de données (en orange les
calculs ab initio DFT) et les valeurs ajustées pour l’auto-interstitiel de fer (en vert le potentiel interatomique ajusté sur la base de données DFT)
et en bleu le potentiel de départ. c) Comparaison entre les calculs ab initio et les valeurs prédites par les deux potentiels EAM pour les nouvelles
configurations d’amas d’interstitiels non incluses dans la base d’ajustement. On peut noter qu’un meilleur accord entre les calculs DFT et EAM
pour les configurations des mono-interstitiels entraîne un accord remarquable pour les nouvelles configurations d’amas d’interstitiels.
alors que les erreurs sont de l’ordre de 20 % ou plus pour cer- nouveau potentiel, qui prend la suite du potentiel proposé en
taines propriétés. Il faut alors se rappeler des objectifs qui moti- [4], permet de décrire avec une grande précision non seule-
vent l’utilisation de potentiels empiriques. L’utilisation des ment la stabilité relative des petits amas d’interstitiels dans le
potentiels empiriques permet d’étendre l’investigation spatio- fer (fig. 23b) mais aussi la dilatation thermique, les barrières de
temporelle d’un solide, ainsi que de son espace configuration- migration des mono-interstitiels et mono-lacunes, les proprié-
nel, mais cela, en partie, aux dépens de la précision. Ainsi, les tés vibrationnelles des défauts, les surfaces et les γ -surfaces
résultats attendus aident à identifier des mécanismes, à dans le fer. Il a été utilisé pour explorer de manière extensive
rendre compte et à explorer des phénoménologies, mais don- le paysage énergétique des interstitiels dans le fer. Il a ainsi
nent rarement accès à des valeurs quantitativement fiables conduit à prédire la stabilité d’une structure cristalline tridimen-
des grandeurs physiques. sionnelle de type phase de Laves C15 pour les amas d’auto-
interstitiels dans ce métal [5], en accord remarquable avec les
Il existe néanmoins des exemples où des résultats quantitatifs calculs ab initio pour ce type d’amas (fig. 23c). Cette prédiction
peuvent être obtenus. Citons le développement récent d’un s’oppose à la morphologie classique des boucles bidimension-
potentiel simple (EAM) pour le cas complexe du fer pur [3]. Ce nelles, sur lesquelles le potentiel avait pourtant été ajusté.
36 Les potentiels interatomiques

La qualité d’un tel potentiel pour le fer pur conduit à en déve-
lopper son domaine de validité. Il pourrait ainsi être intéres-
sant d’intégrer l’oxygène pour accéder au rôle des impuretés
d’oxygène et même à l’oxydation du fer. Ce type de dévelop-
pement impose de décrire dans le même potentiel empirique
tout à la fois des liaisons métalliques et des liaisons iono-cova-
lentes, ce qui revient à introduire un certain caractère chimique
aux liaisons. Les potentiels réactifs de type ReaxFF (Reax
Force Field), par exemple, permettent d’aborder ce type
d’étude. Ils ont déjà été développés pour les systèmes Fe-O-
H et donnent une description prometteuse de l’oxydation du fer
par l’eau. Ce type de méthode – plus coûteux en temps de
calcul – constitue une piste intéressante et complémentaire
aux calculs ab initio. Ces potentiels ouvrent la possibilité de
faire des études sur les interfaces métal-oxyde, que ce soit en
corrosion ou pour les alliages ODS (Oxide Dispersion
Strengthened).
Alain CHARTIER,
Département de physico-chimie
Mihai Cosmin MARINICA et Jean-Paul CROCOMBETTE,
Références
[1] A. CHARTIER, G. CATILLON et J.-P. CROCOMBETTE, Phys. Rev. Lett., 102
(2009), p. 155503.
[2] F. ERCOLESSI et J. B. ADAMS, Eur. Lett., 26 (1994), p. 583.
[3] M.-C. MARINICA et F. WILLAIME, Solid State Phenomena, 129 (2007),

p. 67.
[4] M. I. MENDELEV et al., Phil. Mag., 83 (2003), p. 3977.
[5] M.-C. MARINICA, F. WILLAIME et J.-P. CROCOMBETTE, Phys. Rev. Lett.,

108 (2012), p. 025501.

Le paysage énergétique à l’échelle atomique
L’ énergie totale d’un système sans contrainte est la somme tielle, nous pouvons définir une surface de potentiel thermody-
des énergies potentielle et cinétique de ses constituants. Pour namique adéquat, par exemple l’énergie libre pour le système
la physique des matériaux, toutes les propriétés (état d’équi- couplé à un thermostat. Pour d’autres conditions d’investiga-
libre, changements de phases, évolution de la macrostructure tions (pression ou potentiel chimique constant), il faut utiliser
pendant ou après irradiation, recuit de dommage primaire, le potentiel thermodynamique associé (pour plus de détails,
etc.) découlent de l’énergie potentielle qui est une caractéris- voir supra, pp. 27 et 28, « Les outils de la modélisation méca-
tique intrinsèque du système. Ce concept fut rationalisé sous nique des matériaux pour le nucléaire »). Les informations
la forme de « surface d’énergie potentielle » (ou PES, extraites du paysage énergétique sont utilisables pour le para-
« Potential Energy Surface », en anglais) par R. MARCELIN et métrage de techniques multi-échelles très utilisées pour les
H. EYRING, au début des années 30, pour décrire une réaction études des dommages primaires d’irradiation comme les cal-
chimique comme l’évolution d’un point sur une surface d’éner- culs de Monte-Carlo cinétiques* (voir supra, pp. 127-132, le
gie potentielle dans un espace à 3N dimensions, N étant le chapitre intitulé « La modélisation du dommage primaire d’ir-
nombre d’atomes du système considéré. radiation par dynamique moléculaire »), de champ moyen
auto-cohérent (voir supra, pp. 133-136, le chapitre intitulé
Une surface d’énergie potentielle est souvent présentée « Structure et cinétique des défauts ponctuels dans le carbure
comme un paysage montagneux où les minima correspon- de silicium ») ou de dynamique d’amas* (voir supra, pp. 137-
dent au fond des vallées et les points de cols aux points de 140, le chapitre intitulé « Structure et cinétique des défauts
passage les plus bas en hauteur entre les vallées. d’irradiation dans le fer »).
L’information essentielle pour l’évolution d’un système est
extraite à partir de la distribution de minima et de points de La recherche de minima et de points de col d’une fonction d’un
cols de la surface d’énergie. La distribution de minima déter- grand nombre de variables est un sujet largement étudié par
mine la thermodynamique du système et les points de cols, les physiciens et les mathématiciens. Pratiquement, pour
et les chemins qui relient les minima locaux déterminent l’évo- atteindre un minimum, il est seulement nécessaire de suivre
lution cinétique du système. Selon la distribution de une évolution du système dans la direction de la plus grande
BOLTZMANN, la probabilité de trouver le système dans un état descente (méthode du gradient). Pour localiser les minima,
décroît de manière exponentielle avec l’augmentation de nous pouvons utiliser des algorithmes mathématiques spéci-
l’énergie de l’état. Dans la limite de basse température, cette fiques : soit des méthodes déterministes qui utilisent comme
dépendance exponentielle signifie que le système est le plus information le gradient (dérivée première) et/ou le Hessien
susceptible de se retrouver dans le voisinage du minimum glo- (dérivée seconde) de la fonction, soit des méthodes stochas-
bal. Par conséquent, à basse température, le minimum de la tiques de physique statistique sans gradient, pour les situa-
fonction énergie potentielle donne une bonne description de la tions spécifiques, quand les forces interatomiques sont très
structure atomique dans l’état thermodynamique du système. difficilement calculables (Metropolis, algorithmes génétiques,
Dans la même gamme de température, les chemins qui relient simplex, etc.). Les algorithmes sont généralement plus effi-
les minima locaux nous indiquent les probabilités pour évo- caces s’ils disposent du gradient de la fonction traitée.
luer d’une configuration de minimum à l’autre. C’est exacte- L’algorithme déterministe le plus facile est de suivre le gradient
ment le cas des processus élémentaires de diffusion atomis- de la fonction en s’assurant que l’énergie cinétique du sys-
tiques de défauts ponctuels* ou des dislocations* (voir tème est zéro à chaque pas d’itération. L’application directe de
ci-après), les réactions cinétiques de défaut, annihilation, cette méthode peut prendre différents noms suivant le
croissance ou élimination sur des puits à l’origine de toute évo- contexte de son utilisation : méthode d’EULER, de la descente
lution cinétique de la microstructure. Pour des températures la plus rapide (« steepest descent », en anglais) ou de la dyna-
plus grandes, la notion de paysage énergétique n’est plus mique moléculaire trempée. Cette méthode n’est pas la plus
adaptée et des changements quantitatifs et qualitatifs peuvent efficace, mais elle est très utilisée car ses avantages majeurs
apparaître. À une température suffisamment grande, l’état le sont qu’elle est très robuste et très facile à mettre en œuvre.
plus stable n’est donc pas forcément celui qui possède l’éner- Actuellement, les algorithmes les plus efficaces sont les
gie potentielle la plus faible. Une entropie importante peut sta- méthodes de type quasi-Newton qui utilisent une approxima-
biliser une configuration d’énergie élevée. De la même tion numérique de l’inverse du Hessien. Un exemple est la
manière que nous avons défini une surface d’énergie poten- méthode BROYDEN-FLETCHER-GOLDFARB-SHANNO à mémoire

limitée (L-BFGS). Cette méthode est très efficace dans le (voir infra, p. 123-126, « La thermodynamique des alliages fer-
domaine quadratique de la fonction, mais peut parfois être pié- chrome ») ou étendus (voir infra, p. 169-170, « La modélisa-
gée dans des régions très éloignées du domaine quadratique. tion de la mobilité des dislocations », ainsi que l’exemple ci-
après).
Les points cols sont plus difficiles à déterminer pour la simple
raison qu’il s’agit d’états de transition. Les méthodes de Dans la seconde classe de méthodes, la seule nécessité est
recherche des points col peuvent être divisées en deux caté- de connaître un point de départ, minimum ou point col. Tous
gories : a) celles qui recherchent le point col sur la base de les autres minima et points cols sont trouvés par la méthode
l’interpolation de deux minima et (b) celles qui utilisent unique- de façon automatique. On peut citer : Activation Relaxation
ment des informations locales du paysage énergétique. Technique nouveau (ARTn), Eigenvector-Following (EF),
dimer, Gentle Descent. L’évolution du système dans une
Dans la première classe de méthodes, nous devons avoir une méthode de ce type suit deux étapes : a) à partir d’un mini-
connaissance préalable du paysage énergétique et donc une mum local, le système est poussé vers un point col ; b) à par-
estimation raisonnable du point col. Toutes ces méthodes sont tir d’un point col, le système est relaxé dans un nouveau mini-
une mise en œuvre indirecte du principe de moindre action de mum. Ces méthodes sont très coûteuses en temps de calcul.
LAGRANGE : le chemin suivi par le système entre deux points de À cause de la complexité des paysages énergétiques et de
la PES est celui qui conduit à un extremum de l’action. Le pre- leur paramétrage particulier, ces méthodes n’ont été utilisées
mier pas est la construction d’un chemin initial, par intuition que récemment dans le domaine de la physique de matériaux.
ou, le plus souvent, par une interpolation entre les deux Les études avec ARTn ont permis de mettre en évidence la
minima. Dans cette étape, nous procédons à une discrétisa- complexité du spectre énergétique pour les minima ainsi que
tion du chemin, par plusieurs « images » entre les deux pour les points col des petits amas d’interstitiels dans le fer, le
minima. Le deuxième pas consiste dans l’optimisation du che- Si cristallin et le Si amorphe et SiO2.
min. Chaque méthode a ses propres recettes pour minimiser
l’action. Deux méthodes sont très utilisées dans la commu- Dans ce qui suit, nous développons deux exemples d’applica-
nauté de physique de matériaux : la méthode « drag » et la tions des méthodes décrites précédemment.
méthode « Nudge Elastic Band ». Ces deux méthodes ont été
intensivement utilisées dans la physique de matériaux des Nous présentons dans la figure 24 le paysage énergétique
vingt dernières années pour la diffusion de défauts ponctuels des défauts ponctuels interstitiels dans une représentation de
Énergie (eV)
4,5
3,0
4,0
2,5
3,5
2,0
3,0
2,5
1,5
2,0
1,0 1,5
1,0
0,5
0,5
0 0
Fig. 24. Le paysage énergétique d’un amas de quatre auto-interstitiels dans le fer, c’est-à-dire à quatre atomes supplémentaires dans le réseau
cubique centré. Le paysage énergétique a été déterminé en utilisant la méthode « Activation Relaxation Technique nouveau » [1] et visualisé
avec la technique des graphes : chaque terminaison représente un minimum ou configuration, et chaque nœud du graphe représente un point
col. Nous remarquons plusieurs contributions : en rouge celle des amas avec une structure conventionnelle formée de dumbbells <110> et en
rouge celle des structures de type phase de Laves C15. Pour chaque configuration, les cubes bleus représentent des lacunes de la structure
initiale, et les boules de couleur des atomes interstitiels.
40 Le paysage énergétique à l’échelle atomique

a
yz
yz = 0,02 GPa
yz = 0,4 GPa Hkp(yz) = 0,58 eV
d
Barrière énergétique (eV)

0,5
1
yz
z = [110]
y = [111]
x = [112] Coupe du cristal Hkp(yz) = 0,09 eV
-0,5
b c
-1
Dislocation vis à l’arrêt Dislocation vis glissant
Décrochement
-1,5
Décrochement
0
0 0,2 0,4 0,6 0,8 1,0 1,2
Coordonnées de réaction
Fig. 25. Le paysage énergétique d’une dislocation vis avec les deux positions de minima et le chemin de passage entre les deux configurations.
a) Une boîte de simulation avec des surfaces libres (en orange) qui contient une dislocation vis dans une position de minimum. Les atomes sont
colorés en fonction de leur énergie potentielle. Une coupe du cristal permet la visualisation des atomes à proximité du cœur de la dislocation
(les atomes jaunes). b) la dislocation dans la boîte de simulation dans la même position initiale vue de côté avec la même convention de
couleurs. c) La dislocation vis glissant vers un autre minimum avec une paire de décrochements d) Calcul de la barrière énergétique pour la
formation d’une paire de décrochements sur une dislocation vis dans un modèle de fer pour différentes valeurs de contrainte appliquée à la
cellule de simulation. La méthode utilisée est la méthode Nudge Elastic Band [2]. La hauteur de barrière, notée ici Hkp, est un paramètre
déterminant de la loi de mobilité des dislocations (voir, infra, pp. 169-170, « La modélisation de la mobilité des dislocations »).
graphe pour un amas de quatre auto-interstitiels dans un Dans les métaux de symétrie cubique centrée, comme le fer
métal de structure cubique centré, le fer, en utilisant la pur, les dislocations de type vis glissent suivant un processus
méthode ARTn. Le paysage énergétique des interstitiels est dit « de Peierls », c’est-à-dire par activation thermique de
beaucoup plus compliqué que celui des lacunes qui ne pré- double décrochement le long de la ligne de dislocation (voir la
sente que quelques minima bien définis. Nous pouvons remar- figure 25a-c). La méthode Nudge Elastic Band est utilisée afin
quer les milliers de minima et de points cols. Chaque graphe de calculer la barrière énergétique associée au processus de
correspond à une morphologie bien définie, soit de boucle PEIERLS. La fonction enthalpie est utilisée pour analyser la
bidimensionnelle, soit de structure tridimensionnelle pério- dépendance de cette barrière en fonction de la contrainte
dique. Cette nouvelle structure tridimensionnelle a été récem- résolue appliquée sur les surfaces de nos cellules de simula-
ment déterminée et la structure cristalline sous-jacente cor- tion. Deux exemples de calculs de barrière énergétique pour
respond à la phase de Laves C15 (fig. 24, voir également, la formation de paire de décrochement sur une dislocation vis
infra, pp. 123-126, le chapitre intitulé « La thermodynamique dans le fer ont été reportés en figure 25d. Plus la contrainte
des alliages fer-chrome »). Dans le fer-alpha, ces agrégats appliquée augmente, plus la formation de décrochement est
C15 sont très stables, immobiles et se forment directement facilitée par la diminution de la barrière. Ces calculs nous ser-
dans les cascades de déplacements. Ils peuvent croître en vent de paramètres pour déterminer la loi de mobilité des dis-
capturant des auto-interstitiels. Ils constituent ainsi un nouvel locations (voir infra, pp. 169-170, le chapitre intitulé « La modé-
élément important à prendre en compte dans les prévisions lisation de la mobilité des dislocations »).
de l’évolution microstructurale des matériaux à base de fer
sous irradiation. La structure et la mobilité des amas d’auto-
interstitiels sont une question encore largement ouverte (pour Mihai Cosmin MARINICA, François WILLAIME
plus de détails voir, infra, pp. 123-126, le chapitre intitulé « La et Laurent PROVILLE,
thermodynamique des alliages fer-chrome »). Département des matériaux pour le nucléaire

Références
[1] M.-C. MARINICA et al., Phys. Rev. Lett., 108 (2012), p. 025501.
[2] D. RODNEY et L. PROVILLE, Phys. Rev. B, 79 (2009), p. 094108.
Bibliographie
Calcul ab initio de la plasticité : structure de cœur et mécanismes
de glissement des dislocations vis
DEZERALD (L.), VENTELON (L.), CLOUET (E.), DENOUAL (C.), RODNEY (D.)
et WILLAIME (F.), « Ab initio modelling of the two-dimensional energy
landscape of screw dislocations in bcc transition metals », Phys. Rev.
B, 89 (2014), p. 024104.
42 Le paysage énergétique à l’échelle atomique

Les simulations de dynamique moléculaire
C’ est entre 1950 et 1955 que Enrico FERMI et ses col- le calcul des forces exercées sur chacun des atomes, mais
lègues, John R. PASTA et Stanislaw ULAM, de Los Alamos, réa- autorise, par ailleurs, la réalisation d’événements collectifs ;
lisèrent les expériences numériques consistant à intégrer les
équations de la mécanique newtonienne pour une chaîne • d’autre part, la dynamique moléculaire permet de modéliser
d’oscillateurs couplés par des interactions non linéaires des phénomènes hors d’équilibre à l’échelle atomique. Par
(potentiels d’interactions non quadratiques) [1]. Les 18 000 exemple, on peut noter l’étude d’endommagements méca-
lampes de l’ordinateur MANIAC permettaient alors de traiter niques [5], de cascades de déplacements atomiques provo-
des systèmes unidimensionnels comportant quelques quées par irradiation (voir infra, pp. 127-132, le chapitre inti-
dizaines d’atomes. Cette étude heuristique, selon Enrico FERMI tulé « La modélisation du dommage primaire d’irradiation par
lui-même, peut être considérée comme la première simulation dynamique moléculaire ») ou encore la formation de couche
de dynamique moléculaire de solides modèles. Avec ces cal- d’oxyde [6], la diffusion d’atomes isolés sur certaines sur-
culs numériques sur une chaîne unidimensionelle d’oscilla- faces [7], la mobilité des dislocations (voir infra, pp. 169-170,
teurs, l’objectif de FERMI [1] était d’étudier la thermalisation des le chapitre intitulé « La modélisation de la mobilité des dislo-
solides en tenant compte de la non-linéarité des interactions cations »). Cette liste est non exhaustive.
entre atomes. À la grande surprise de FERMI et de ses col-
lègues, la chaîne d’atomes ne thermalisait pas dans les temps L’explosion de la capacité de calcul des ordinateurs a permis
accessibles à la simulation et il fallut de nombreuses études, non seulement d’augmenter considérablement les volumes de
par la suite, afin de comprendre le processus menant à l’équi- matière étudiés, mais aussi d’améliorer la description des
libre thermodynamique [2]. Cette première étude fondait la interactions atomiques.Très tôt, ces interactions furent décrites
méthode dite « de dynamique moléculaire* », qui consiste à par des fonctionnelles de différentes formes analytiques, prin-
considérer, dans le cadre de la mécanique classique, des cipalement les interactions dites « à deux corps » et celles
atomes ou des ions en interaction et de suivre pas à pas leurs dites « à N-corps » (voir supra, pp. 33-37, le chapitre intitulé
mouvements par discrétisation du temps. Une représentation « Les potentiels interatomiques »). La première famille com-
schématique est donnée dans la figure 26. Compte tenu de la porte des potentiels modèles de type LENNARD-JONES, mais
rapidité des vibrations atomiques, le pas en temps des calculs également des potentiels permettant de décrire les interac-
de dynamique moléculaire est extrêmement court (de l’ordre tions dans les solides ioniques. La seconde famille est consti-
de la femto-seconde). De nombreux physiciens et chimistes tuée de potentiels de forme plus complexe (d’où le nom à N-
utilisent cette méthode, désormais standard. On peut classer
les calculs de dynamique moléculaire en deux grandes
familles :
• D’une part, l’étude des systèmes atomiques proches de

l’équilibre thermodynamique et les propriétés associées.
L’idée est alors d’imposer de façon ad hoc cet équilibre.
Différentes techniques permettent de réaliser cette opéra-
tion, ouvrant alors la voie vers l’étude des propriétés thermo-
dynamiques à l’échelle atomique. L’une des toutes premières
applications fut sans doute la mise en évidence théorique de
la transition de phase solide-liquide dans un modèle idéal de
sphères dures [3]. Pour toutes les études portant sur les pro-
priétés thermodynamiques, la dynamique moléculaire est à
comparer avec la méthode de Monte-Carlo [4]. Selon le type
de système étudié, l’une des deux méthodes sera plus effi-
cace. En effet, chaque incrément de l’algorithme de Monte-
Fig. 26. Représentation schématique de la dynamique moléculaire
Carlo nécessite un seul calcul de l’énergie totale, mais n’au- pour un ensemble de particules (cercles gris) interagissant via un
torise la réalisation que d’un seul et unique événement, alors champ de force (symbolisé par les cercles colorés) et se déplaçant
que la dynamique moléculaire requiert à chaque incrément dans les directions (flèches) déterminées par le champ de force.

corps). On peut citer les potentiels de type Embedded-Atom- ments. En ce qui concerne le calcul de processus non thermi-
Method (EAM) pour décrire les métaux, ainsi que les poten- quement activés comme la diffusion de flux de chaleur [12] ou
tiels traduisant les liaisons covalentes (du type potentiel de la viscosité des liquides ou des verres, nous retrouvons, en
Tersoff). Voir supra, pp. 33-38, « Les potentiels interato- quelque sorte, la problématique mise à jour par FERMI, PASTA
miques ». De nos jours, on trouve également de nombreuses et ULAM. Il s’agit alors d’étudier l’évolution d’un système qui est
études réalisées sur la base des calculs de structure électro- maintenu éloigné de son équilibre thermodynamique, La simu-
niques et ab initio permettant de décrire correctement les liai- lation directe de ce type de processus reste particulièrement
sons chimiques au sein des molécules et des solides [8]. difficile et ne permet actuellement d’établir des prédictions
L’intégration des degrés de liberté électroniques diminue de quantitatives que dans de très rares cas [13]. La combinaison
manière considérable les volumes qui peuvent être traités de la théorie statistique de la réponse linéaire [14] et des simu-
mais, en contrepartie, cela ouvre un vaste champ d’analyse lations atomistiques permet également de progresser dans ce
lié aux propriétés électroniques (voir supra, pp. 29-32, « Les domaine [12].
calculs électroniques ab initio pour la matière condensée »).
Les calculs de dynamique moléculaire peuvent être relative- Laurent PROVILLE, Mihai Cosmin MARINICA
ment lourds. Ils sont limités soit par le volume de matière, et Jean-Paul CROCOMBETTE,
c’est-à-dire le nombre d’atomes considérés, soit par le temps Département des matériaux pour le nucléaire
physique que l’on souhaite simuler. Avec une description des

interactions de type EAM, les tailles d’échantillon numérique
qui peuvent être traitées sont en constante progression ;
actuellement, nous pouvons les estimer à plusieurs millions Références
d’atomes sur une durée de quelques microsecondes ou bien [1] E. FERMI, J. R. PASTA et S. M. ULAM, « Studies of nonlinear pro-
à quelques centaines de millions d’atomes sur quelques pico- blems », Los Alamos Scientific Laboratory Rapport No. LA-1940
secondes. Ces capacités, bien qu’étant encore en-deçà des (1955). Réédité dans Collected papers (University of Chicago Press,
Chicago, 1965), p. 978.
échelles macroscopiques d’espace et de temps, permettent
toutefois de couvrir un spectre d’étude relativement large, car [2] A.C. SCOTT, Nonlinear Science : Emergence and Dynamics of
de nombreux processus physiques sont initiés, voire achevés, Coherent Structures, 2e édition, (Oxford University Press, Oxford,
sur des échelles compatibles avec la dynamique moléculaire. 2003).
Outre les études portant sur des systèmes classiques, dont la [3] B.J. ALDER et T.E. WAINWRIGHT, J. Chem. Phys., 27 (1957), p. 1280.
dynamique est dite « newtonienne », nous pouvons également
[4] N. METROPOLIS et al., J. Chem. Phys., 21 (1953), p. 1087.
mentionner l’étude de systèmes quantiques, dont la détermi-
nation des propriétés d’équilibre est rendue possible grâce aux [5] R. K. KALIA et al., Int. J. of Fracture, 121 (2003), p. 71.
intégrales de chemin de FEYNMAN et à l’isomorphisme entre la
[6] J. DALLA TORRE et al., J. Appl. Phys., 92 (2002), p. 1084.
physique statistique d’une particule quantique et celle d’un
polymère classique [9]. On trouvera, par exemple, l’étude des [7] M.-C. MARINICA et al., Phys. Rev. B, 70 (2004), p. 75415.
électrons piégés dans les centres colorés de sels fondus [9]. [8] R. CAR et M. PARRINELLO, Phys. Rev. Lett., 60 (1988), p. 204.
Il est donc possible de couvrir un grand nombre de problèmes
[9] M. PARINELLO et A. RAHMAN, J. Chem. Phys., 80 (1984), p. 860.
physiques dans différents domaines de température et de
pression. [10] A.F. VOTER, Phys. Rev. B, 57 (1998), p. 13985.
[11] N. MOUSSEAU et al., Journal of Atomic, Molecular and Optical

Hormis l’amélioration des théories décrivant les interactions Physics, 2012 (2012), p. 925278.
entre atomes (voir supra, pp. 33-37, « Les potentiels interato-
[12] J.-P. CROCOMBETTE et L. PROVILLE, Appl. Phys. Lett., 98 (2011),
miques »), les développements méthodologiques sont perma-
p. 191905.
nents. Ils portent essentiellement sur les limites temporelles
de la dynamique moléculaire. Nous décrivons ici quelques [13] G. CICCOTTI et G. JACUCCI, Phys. Rev. Lett., 35 (1975), p. 789.
exemples choisis afin d’illustrer ces développements, sans
[14] R. KUBO, J. Phys. Soc. (Jpn), 12 (1957), p. 570.
toutefois être exhaustif. Plusieurs types de méthodes d’accé-
lération sont possibles suivant la nature du problème physique
abordé. Le premier exemple concerne les processus thermi-
quement activés (voir infra, pp. 169-170, « La modélisation de
la mobilité des dislocations »), pour lesquels l’échelle de temps
est beaucoup plus longue que le temps caractéristique des
vibrations du cristal [10,11]. Il est alors possible de traiter les
deux échelles de temps séparément et d’établir des prédic-
tions via des méthodes de physiques statistiques comme la
théorie de l’état de transition et ses nombreux développe-
44 Les simulations de dynamique moléculaire

Potentiels thermodynamiques
et potentiels de force moyenne
L e premier objectif de la simulation atomistique sur ordina- 8(5 –1 , +, /) = -–〈〉/5 = 1 2 34∑ & exp<– & /1 2 5= (2)
teur est de prévoir la stabilité des phases d’un matériau au
cours du temps en fonction de la température, de la pression La somme, appelée « fonction de partition », comprend tous
et des potentiels chimiques de chacun de ses constituants. Le les états & du système. Les termes exponentiels correspon-
second objectif, tout aussi important, est de calculer les dent aux poids de Boltzmann non normalisés. L’ensemble
constantes de réaction associées aux mécanismes atomiques canonique autorise les échanges d’énergie sous forme de
contrôlant l’évolution de la microstructure. Le moyen le plus chaleur entre la boîte de simulation et un thermostat, ce qui se
direct pour résoudre ces deux problèmes consiste à laisser le traduit par des fluctuations de l’énergie du système autour
système évoluer au cours du temps en résolvant l’équation de d’une valeur moyenne 〈〉.
mouvement de Newton, éventuellement couplée à un thermo-
stat et un barostat. Malheureusement, cette approche néces- Notons que l’on aurait également pu effectuer la transformée
site des temps de calcul importants et ne peut être utilisée que de Legendre de l’énergie par rapport à l’entropie du système
dans des cas simples. La raison en est que les processus ther- en considérant cette dernière comme une variable d’état. La
miquement activés, qui contrôlent les mécanismes atomiques fonction thermodynamique correspondante est alors l’énergie
et le transport dans les matériaux, sont très rarement obser- libre :
vés à l’échelle du pas de temps de la dynamique moléculaire,
lequel est typiquement de l’ordre de la femtoseconde (10-15 s). = – 5- (3)
En pratique, il est presque toujours préférable de calculer un Cette fonction, obtenue en utilisant la relation 5 = ∂/∂-, est
potentiel thermodynamique en fonction des paramètres de équivalente à – 58. L’énergie libre est plus communément
contrôle du système. Cela permet, par la suite, de déterminer considérée dans le domaine de la simulation moléculaire. Elle
les conditions de coexistence à l’équilibre thermodynamique possède une analogie évidente avec le principe de stabilité en
en maximisant les fonctions thermodynamiques appropriées. mécanique classique. Mécaniquement, un système est stable
C’est ce qui est fait, en pratique, dans un alliage binaire dans si son énergie potentielle est dans un minimum local.
lequel deux phases coexistent lorsque l’on met en œuvre la Thermodynamiquement, une phase est stable si son énergie
construction de MAXWELL. libre est plus basse que celle de toutes les autres phases pos-
sibles. À température non nulle, la phase la plus stable n’est
Le potentiel thermodynamique le plus simple est celui qui est donc pas forcément celle qui possède l’énergie interne la plus
associé à l’ensemble microcanonique et qui correspond à l’en- faible : une entropie importante peut stabiliser une phase
tropie de Boltzmann du système : d’énergie interne élevée. C’est le cas pour les structures for-
tement dégénérées présentant une multitude de minima éner-
-(, + ,/) = 1 2 34⍀( ,+ ,/), (1) gétiques voisins et des modes de vibration très lents.
où ⍀(, +, /) est le nombre d’états quantiques d’un système Les similarités entre stabilité mécanique et thermodynamique
dans lequel on impose la valeur de l’énergie, le nombre + nécessitent de clarifier ce que l’on entend par « énergie libre
de particules et le volume /. La quantité ⍀(, +, /) est aussi d’une phase » pour des systèmes de taille finie. En simulation
l’inverse de la densité d’état. La constante 1 2 est la constante numérique, on doit tout d’abord être capable de distinguer les
différentes phases possibles dans un système atomique. Pour
cela, on a recours à un paramètre d’ordre >. Il s’agit, le plus
de Boltzmann. Si l’on se restreint au cadre de la mécanique
classique, l’énergie du système est la somme de l’énergie
cinétique et d’un potentiel interatomique empirique. souvent, d’une fonction réelle des positions atomiques dont
Formellement, une ou plusieurs transformées de Legendre de les valeurs varient de façon continue quand on passe d’une
l’entropie de Boltzmann par rapport aux grandeurs extensives phase à l’autre. Pour un alliage présentant une lacune de mis-
(énergie, volume, nombre de particules…) permet d’obtenir cibilité et simulé dans l’ensemble pseudo-grand canonique, la
les potentiels thermodynamiques des autres ensembles. Par concentration jouera parfaitement ce rôle. Pour un alliage
exemple, le potentiel thermodynamique de l’ensemble cano- ordonné, le paramètre d’ordre pourra être une fonction des
nique à la température inverse 5-1 = ∂-/∂ s’écrit comme une occupations atomiques sur les sous-réseaux de la phase
fonction de Massieu : ordonnée. Dans le cas d’une transition structurale, on utilisera

plutôt un paramètre d’ordre basé sur l’orientation des liaisons En examinant les minima d’énergie libre apparaissant sur la
atomiques. On définit l’énergie libre de LANDAU par rapport au figure 27 en fonction de la température, on constate que la
paramètre d’ordre de la façon suivante : phase liquide (L) est la plus stable pour une température
réduite T supérieure à Tls = 0,17 (l’unité des températures est
(>) = –12 534∑ > ∈∂⍀(>)exp<– & /1 2 5=, (4) la profondeur du potentiel Lennard-Jones divisé par 12). À la
&
température Tls, les phases liquide et solide coexistent, la
où ∂⍀(>) englobe toutes les configurations & pour lesquelles phase solide icosaédrique restant la plus stable aux tempéra-
le paramètre >& est compris entre > – >/2 et > + >/2, > tures supérieures à Tss = 0,12. Aux températures inférieures,
étant la largeur de l’intervalle de l’histogramme. L’énergie libre la phase octaédrique, associée au minimum énergétique glo-
de Landau est également appelée « potentiel de force bal du système, devient la plus stable.
moyenne ». Elle s’obtient en évaluant numériquement l’histo-
gramme : La connaissance de la surface d’énergie libre en fonction de
la température et d’un paramètre d’ordre permet non seule-
A(>) = exp<B – (>)C/12 5=, (5) ment de caractériser les stabilités relatives des phases pos-
sibles dans un système mais aussi de calculer, dans certaines
associé à la probabilité d’observer le système atomique avec conditions, les constantes gouvernant la cinétique du système.
une valeur >& du paramètre d’ordre dans l’intervalle défini En effet, pour une réaction d’un bassin initial A vers un bassin
autour de >. d’arrivée B séparés par une barrière d’énergie libre F calcu-
lée par rapport au paramètre d’ordre, la constante de réaction
Calculer un potentiel thermodynamique ou un potentiel de force k est reliée à la barrière d’énergie libre via la formule suivante :
moyenne est une tâche difficile et un domaine de recherche
125 F
1 = к H —— J exp H – —— J
très actif. La difficulté vient du fait qu’il n’est pas possible de cal-
ℎ 125
culer les fonctions de partition par quadrature, étant donné la (6)
forte dimensionnalité du système. La solution universellement
utilisée consiste à effectuer un échantillonnage pondéré : une où ℎ est la constante de Planck, et к est le facteur de transmis-
méthode Monte-Carlo produit les configurations contribuant le sion. Cette quantité est définie comme la fraction des trajec-
plus à la grandeur que l’on souhaite évaluer. Un estimateur per- toires initiées à partir du sommet de la barrière et qui attei-
met ensuite d’obtenir les moyennes en les repondérant judi- gnent le bassin B avant le bassin A. Si cette quantité n’est pas
cieusement. La moyenne statistique utilisée traduit une rela- trop proche de 0 ni de 1, elle peut être estimée relativement
tion de perturbation thermodynamique pour le calcul d’une facilement à partir d’un échantillon de trajectoires. Le para-
énergie libre, ou bien est associée à la construction de l’histo- mètre d’ordre utilisé dans la simulation est alors une bonne
gramme pour le calcul d’une énergie libre de Landau. Chacune
des deux moyennes peut être prise en laissant le système évo-
luer hors ou à l’équilibre thermodynamique. Pour une dyna-
mique moléculaire hors équilibre, une action extérieure force D
I
mécaniquement le système à visiter les régions inexplorées le
O F(Q6, T)
long du paramètre d’ordre, ou bien transforme progressivement
3
les contraintes thermodynamiques qui s’exercent naturellement
0,06 2,5
sur le système à l’équilibre. La variation du potentiel de force
0,08 2
moyenne est directement reliée aux travaux mécaniques utili-
sés pour forcer le système. De même, la variation du potentiel 0,10 1,5
thermodynamique s’obtient à partir des travaux thermodyna- F(Q6, T)

0,12 1
miques récupérés lorsque l’on transforme le système. La 3

0,14 0,5
méthode d’intégration thermodynamique apparaît comme le 2 0,16 0
cas particulier d’une dynamique moléculaire hors-équilibre 1 0,18

0 T
dans laquelle les paramètres de contrôle sont transformés très L 0,20
0
lentement de manière quasi statique. 0,22
0 0,1 0,2 0,3 0,4 0,5 0,6
Q6
À titre d’illustration, nous avons mis en œuvre la méthode
Monte-Carlo fondée sur des dynamiques moléculaires forcées Fig. 27. Énergie libre de LANDAU de l’agrégat LJ38 en fonction
mécaniquement dans un agrégat monoatomique de 38 du paramètre d’ordre (ici noté Q6) et de la température. Les bassins
atomes [1]. Bien que le potentiel interatomique, de type d’attraction correspondant aux phases liquide (L), icosaédrique (I) et
octaédrique (O) apparaissent clairement à toutes les températures
Lennard-Jones, soit très simple, ce système est très riche ther-
pour les valeurs du paramètre d’ordre comprises dans les intervalles
modynamiquement [2], comme le montre le paysage d’éner- respectifs [0,02, 0,06], [0,11, 0,13] et [0,5, 0,56]. D est une structure
gie libre de la figure 27, porté en fonction de la température et
du paramètre d’ordre (ici noté Q6).
octaédrique présentant une faute d’empilement qui est observée
lors des transitions entre les structures octaédrique et icosaédrique.
46 Potentiels thermodynamiques et potentiels de force moyenne

coordonnée de réaction : il permet de bien caractériser le
degré d’avancement d’un mécanisme réactif. Dans le cas En avant
- log (A(ξ))
15 En arrière
contraire, le paramètre d’ordre n’est pas une bonne coordon- Moyen
née de réaction : le sommet de la surface d’énergie libre décrit 12
incorrectement l’état de transition.

9
Ainsi, dans l’exemple de l’agrégat LJ38, les paramètres 6

d’ordre structuraux les plus simples sont connus pour être de
mauvaises coordonnées de réaction pour la transition struc- 3
turale entre les phases icosaédrique et octaédrique. La sur-
M(ξ)
0,30
face d’énergie libre de la figure 27, page précédente, ne per- 0,25
met pas ainsi d’extraire les constantes de réaction associées 0,20
à cette transition structurale. La raison en est que les états 0,15
de transition et certains minima locaux ont les mêmes valeurs 0,10
des paramètres d’ordre structuraux. Dans cette situation, le 0,05
maximum de la barrière d’énergie libre se trouve donc 0,00
ξ
0 5.10-11 10-10 1,5.10-10 2.10-10 2,5.10-10
déplacé dans le bassin d’attraction du minimum local. La
détermination d’un paramètre d’ordre pertinent peut donc
s’avérer difficile. Fig. 28. Migration lacunaire dans le fer-alpha. La position de l’atome
migrant vers le site vacant est représentée par la coordonnée
de réaction ξ (exprimée en mètres). La probabilité de la configuration
le long d’une migration est A(ξ). Le cologarithme de la moyenne
L’étude de la migration d’un défaut lacunaire dans le fer est
de A(ξ) donne accès au profil d’énergie libre. La moyenne du

nettement moins problématique. Dans une boîte de simulation
cologarithme de A(ξ) est notée M(ξ) et correspond à une divergence.
de fer pur de structure cubique centrée comprenant une seule
lacune, la position d’un atome voisin de la lacune est une Elle renseigne sur la précision des calculs : les estimations sont plus
bonne coordonnée de réaction, qui décrit correctement le précises lorsque l’atome de fer est dans une des deux positions
degré d’avancement de la migration de l’atome vers le site ini- stables ou dans la position métastable entre les deux bosses de la
barrière de migration.
tialement vacant. Nous avons calculé le potentiel de force
moyenne pour la migration lacunaire en utilisant la même
méthode que précédemment. Le Fe- est modélisé avec le
0,7
potentiel de ACKLAND ET MENDELEV [3]. Les dynamiques molé-
culaires forcées sont produites en tirant sur un atome voisin
0,6
0,5
de la lacune à l’aide d’un potentiel harmonique, de façon à ce 0,4
que l’atome échange sa position avec celle de la lacune. À F(ξ, T) [eV] 0,3
partir de la force exercée sur l’atome, la probabilité d’équilibre 0,8
(A(ξ)) de sa position ξ le long d’une direction [111] (passant

0,2
0,7
0,6 0,1
0,5
par l’atome et le site vacant) peut être reconstruite. La figure 28 0,4
0,3
illustre cette reconstruction pour des dynamiques moléculaires 0,2
0,1
générées à partir des deux bassins stables de l’énergie libre 0
(forward ou backward).
3.10-10
La moyenne du cologarithme de A(ξ) définit une divergence

1 000 2.10-10
M(ξ) qui indique les fluctuations des mesures numériques en

900
800 2.10-10
fonction de la coordonnée de réaction ξ. Cette quantité donne

700 2.10-10
ξ [m]
600
500 1.10-10
T [K] 400
5.10-11
indirectement accès à la précision de la moyenne directe de 300
A(ξ) (Moyen sur la figure 28) calculée sur un grand nombre de

200 0.100
100
0 -5.10-11
trajectoires (ici 100). Le cologarithme de cette moyenne est

l’énergie libre de Landau représentée sur la figure 29 pour des Fig. 29. Migration lacunaire dans le fer alpha. Paysage d’énergie libre
températures variant de 20 K à 1 000 K. en fonction de la température et de la coordonnée de réaction ξ
représentant la position de l’atome migrant vers le site vacant
(en mètres).
À basse température, le profil de l’énergie libre de la lacune
présente les deux bosses attendues [3]. L’atome voisin de la
lacune doit effectuer deux sauts pour s’échanger avec la
lacune. Lorsque la température augmente, la double bosse À partir des données de la figure 29, les barrières d’énergie
s’estompe puis disparaît complètement. L’atome voisin ne doit libre sont obtenues en fonction de la température et peuvent
plus franchir qu’un seul col pour migrer. être comparées aux valeurs calculées dans le cadre classique
de l’approximation harmonique [1,2] (fig. 30).

rait par un coût en temps de calcul considérable, hors de por-
0,65
Monte-Carlo tée des ordinateurs actuels. Pour l’estimation d’un taux de
Énergie libre de migration (eV)
CHA migration ou d’un coefficient de diffusion, la solution la plus

0,60 judicieuse à l’heure actuelle est d’obtenir la contribution entro-
pique dans le cadre de l’approximation harmonique, à partir
0,55 des spectres de phonons. Une telle prise en compte de l’en-
tropie de migration harmonique introduit un facteur correctif
de l’ordre de 102 à 103 dans le calcul du taux de migration,
0,50
dans le cadre de la théorie de l’état de transition [1], comparé
à l’utilisation souvent faite de la fréquence de DEBYE pour esti-
0,45 mer la contribution entropique [2]. Cette correction est nette-
ment plus importante que celle due aux effets anharmoniques,
0,40 qui doivent cependant être pris en compte à haute tempéra-
0 100 200 300 400 500 600 700 800 900 1 000 ture. Nous préconisons alors la méthode décrite ci-dessus en
Température (K) utilisant des potentiels semi-empiriques ajustés sur les calculs
ab initio.
Fig. 30. Énergie libre de migration de la lacune dans le fer- en
fonction de la température : les résultats des simulations Monte-
Carlo sont comparés à ceux issus de l’approximation harmonique
Manuel ATHÈNES et Mihai Cosmin MARINICA,
classique (CHA).
À la température de service des réacteurs, la prise en compte

Références
de l’anharmonicité (fig. 30) augmente d’un ordre de grandeur
la mobilité de la lacune dans le fer par rapport aux calculs [1] D.J. WALES, Energy Landscapes, Cambridge UP (2003).
usuels effectués dans le cadre de l’approximation harmonique. [2] M. ATHÈNES et M.-C. MARINICA, « Free energy reconstruction from
Cela montre tout l’intérêt de quantifier de la sorte les effets steered dynamics without post-processing », J. Comp. Phys., 229
anharmoniques et de les prendre en compte, par la suite, dans (2010), pp. 7129-7146.
les modèles cinétiques. [3] G. ACKLAND, M. MENDELEV, D. SROLOWITZ, S. HAN et A. BARASHEV,
« Development of an interatomic potential for phosphorous impurities
Le calcul des potentiels thermodynamiques et des potentiels in alpha-Fe », J. Phys. : Cond. Matter, 16 (2004), S2629.
de force moyenne est une composante importante de la simu-
lation des matériaux à l’échelle atomique. Cette approche
reste cependant cantonnée aux systèmes pouvant être décrits
de façon suffisamment réaliste avec des potentiels semi-empi-
riques. En effet, mieux décrire les interactions atomiques
nécessiterait de combiner les codes de calcul de potentiel
thermodynamique et de structure électronique. Cela se tradui-
48 Potentiels thermodynamiques et potentiels de force moyenne

Les diagrammes de phase
A fin de modéliser les propriétés d’équilibre des alliages, en construire la thermodynamique d’un alliage binaire à partir
particulier leur diagramme de phase*, il faut une approche d’un seul paramètre Q d’interaction entre premiers voisins.
permettant de calculer les fonctions thermodynamiques asso- L’énergie d’une configuration est donnée par :
1
ciées, telles que l’énergie libre, et non seulement l’énergie à
= — R QS&S),
2 &,)
0 K. Il est vrai que les calculs de dynamique moléculaire ab ini- (1)
tio peuvent être réalisés à température finie, mais la partie de
où la somme est restreinte aux sites & et ) premiers voisins. La

l’espace des phases qui est ainsi explorée est généralement
variable de pseudo-spin S& prend la valeur +1 (respectivement

beaucoup trop restreinte pour que l’on puisse en déduire les
propriétés d’équilibre du système. Cette limitation, liée au
-1) si le site est occupé par un atome de type A (respective-
temps de calcul, peut être levée par l’emploi de potentiels
ment B). L’alliage est à tendance à l’ordre pour un paramètre
d’interaction Q >0, et à tendance à la démixtion pour Q < 0. Le
empiriques : des méthodes de simulation existent pour calcu-
ler, à partir de ces potentiels empiriques, les différentes fonc-
modèle peut ensuite être facilement généralisé pour inclure
tions thermodynamiques et construire ainsi un diagramme de
des interactions de paire allant au-delà des premiers voisins,
phase. Cependant, peu de potentiels empiriques existent pour
ainsi que des interactions de multiplets mettant en jeu plus de
les alliages, notamment au-delà des alliages binaires, et ces
deux sites cristallins. Ces interactions de multiplets sont par-
potentiels souffrent d’un manque de robustesse quant à leur
ticulièrement importantes pour dissymétriser la thermodyna-
prédictivité. En outre, les approches thermodynamiques utili-
mique de l’alliage : avec seulement des interactions de paire,
sant des potentiels empiriques demeurent très lourdes en
le comportement thermodynamique est le même des côtés
temps de calcul. Nous présentons dans ce chapitre deux
riche en A et riche en B du diagramme de phase (fig. 31).
approches alternatives permettant toutes deux de modéliser
la thermodynamique d’un système tel qu’un alliage à partir de
Enfin, ce modèle ne se limite pas aux alliages binaires et se
données expérimentales ou de calculs atomiques. La pre-
généralise aux alliages avec plus de deux constituants.
mière approche utilise une description atomistique des diffé-
Les paramètres d’interaction Q apparaissant dans le modèle

rentes configurations possibles d’un alliage, via un modèle de
réseau rigide, pour construire ensuite la thermodynamique du
d’Ising [Eq. (1)] peuvent être déduits de calculs ab initio. Pour
système grâce à la physique statistique. La seconde
approche, la méthode CALPHAD, est fondée sur la thermo-
dynamique macroscopique. Paires
Paires + multiplets
T / T0
La thermodynamique du réseau 4
rigide
3 c.f.c. L12 L10 L12 c.f.c.
Nous nous intéressons ici à la modélisation d’un alliage cris-
tallin. À température non nulle, les atomes vibrent autour de
positions moyennes occupant les nœuds d’un réseau pério- 2
dique. Toutes les configurations de l’alliage peuvent donc être
décrites en précisant le type d’atomes occupant chacun des 1
nœuds du réseau. Il devient alors très simple de construire un
modèle prédisant l’énergie de chacune de ces configurations.
0
Ce modèle est ensuite utilisé pour en déduire la thermodyna-
0 0,25 0,5 0,75 1
mique et les propriétés d’équilibre de l’alliage [1,2]. Fraction atomique de soluté
Fig. 31. Diagramme de phase pour un alliage modèle à tendance

Le modèle d’Ising
voisins ( Q1 > 0 et Q2 < 0), puis en y ajoutant des interactions de
à l’ordre sur un réseau cfc avec des interactions de paire 1ers et 2nds
Initialement développé pour les systèmes magnétiques, le
multiplets (ici le triangle et le tétraèdre de 1ers voisins). Le diagramme
modèle d’Ising est le modèle le plus simple permettant de
de phase est calculé en Cluster Variation Method.

U (1–2 )2
( , 5) = — Q —–—— – 15 V 34( ) + (1– )34(1– )X,
cela, les énergies de formation d’un certain nombre de struc-
2 4
(2)
tures ordonnées sont calculées ab initio et les paramètres d’in-
teraction sont ensuite obtenus par ajustement du modèle sur
la base de données ainsi constituée. On effectue alors un où U est le nombre de premiers voisins d’un atome et 1 la
« développement en amas » de l’énergie de l’alliage. constante de Boltzmann. La même approximation de champ
moyen peut être utilisée pour des structures ordonnées à
Un développement en amas n’est utile que si les interactions longue distance en introduisant des sous-réseaux de compo-
sont à courte portée, permettant ainsi de limiter à un nombre sition nominale différente. Il est également possible de traiter
restreint les interactions apparaissant dans le modèle d’ISING des hétérogénéités telles que celles produites par une inter-
généralisé. En présence de forts effets élastiques, cela n’est face ou une surface en considérant une composition variant
pas possible puisque l’interaction élastique entre atomes est spatialement.
à longue portée : elle décroît comme l’inverse du cube de la
distance. Une solution mathématique consiste alors à complé- Néanmoins, cette approximation néglige toutes les corréla-
ter le développement en amas par des interactions calculées tions pouvant exister à courte distance entre les atomes, c’est-
dans l’espace réciproque. Il est également possible de faire à-dire l’ordre à courte distance. Il est possible d’améliorer cette
appel aux fonctions de GREEN du réseau cristallin pour traiter approximation de champ moyen en considérant des probabi-
ces contributions élastiques. lités sur un amas plus gros que le point, afin de traiter de façon
exacte l’ensemble des corrélations existant à l’intérieur de cet
L’entropie de configuration amas. On obtient ainsi la méthode de variation des amas
(CVM pour Cluster Variation Method). Pour un modèle énergé-
Le modèle énergétique ainsi fixé, il est possible de calculer tique donné, cette méthode tend vers la solution exacte, pré-
rapidement l’énergie de toute configuration de l’alliage. On dite, par exemple, par Monte-Carlo, à mesure que la taille de
peut ainsi échantillonner l’espace des configurations du sys- l’amas maximal augmente (fig. 32).
tème pour en déduire l’énergie libre associée à l’entropie de
configuration et calculer ensuite les propriétés d’équilibre du D’autres méthodes peuvent également être utilisées pour
système, telles que le diagramme de phase. modéliser la thermodynamique d’un alliage à partir d’un
modèle énergétique de type Ising comme les développements
La façon la plus simple d’effectuer cet échantillonnage est basse et haute température.
d’utiliser un algorithme de Monte-Carlo thermodynamique.
Dans de telles simulations, des échanges entre des atomes de
différents types sont réalisés avec des taux d’acceptation don-
nés par les poids de Boltzmann calculés à partir des énergies
des différentes configurations. Il est alors possible de calculer
| j\\ |
kBT
le potentiel chimique pour une composition fixée de l’alliage Bragg-Williams

2
(simulation dans l’ensemble canonique), ou bien, à l’inverse,
d’obtenir la composition pour un potentiel chimique fixé
(ensemble semi-grand canonique). Cette méthode, très géné-
rale, présente le désavantage d’être purement numérique. On
obtient en sortie les valeurs des différentes grandeurs ther- 1 A3B AB AB3
modynamiques, mais pas d’expression analytique qu’il serait

possible d’utiliser, par la suite, dans des modèles aux échelles
supérieures.
0
Les approximations de champ moyen, couramment utilisées, 0 0,25 0,50 0,75 1 CB
permettent d’obtenir de telles expressions analytiques.
| j\\ |
kBT
L’approximation la plus courante est celle de BRAGG-WILLIAMS

2,0 CVM
(ou « champ moyen de point »). Elle considère que la proba-
bilité d’avoir une espèce chimique sur un site est indépendante
de son voisinage et donc égale à la composition nominale de 1,5
A3B AB AB3
l’alliage. Cela permet d’obtenir une expression simple de l’en-
tropie de configuration et donc de l’énergie libre de l’alliage. 1,0
CB
Ainsi, pour un alliage binaire dont l’énergie de chaque confi- 0 0,25 0,50 0,75
guration est donnée par l’équation (1), page précédente,

Fig. 32. Diagramme de phase pour un alliage à tendance à l’ordre
l’énergie libre par atome pour une solution solide de compo-
sition nominale et une température 5 sera :
sur un réseau cfc avec des interactions de paire premiers voisins,
calculé soit avec l’approximation de Bragg-Williams, soit avec une
approximation CVM (Cluster Variation Method).
50 Les diagrammes de phase

Les autres sources d’entropie La méthode CALPHAD
Les méthodes précédentes permettent de modéliser la ther- La méthode CALPHAD* (Calculation of Phase Diagrams) [3]
modynamique associée à l’entropie de configuration de l’al- est une méthode semi-empirique de calcul de diagrammes de
liage. D’autres sources d’entropie existent, liées par exemple phases de systèmes multi-éléments. Cette méthode permet
à des excitations électroniques ou magnétiques, ou bien aux d’obtenir une description de l’enthalpie libre de toutes les
vibrations des atomes autour de leur position d’équilibre. Pour phases d’un système en fonction de la composition, de la tem-
les excitations électroniques et les vibrations, l’échelle de pérature et de la pression à travers des modèles thermodyna-
temps mise en jeu est beaucoup plus petite que celle typique miques.
d’un échange entre deux atomes intervenant dans l’entropie
de configuration. Il est généralement possible de faire une La combinaison de systèmes simples (binaires ou ternaires)
approximation adiabatique afin de séparer ces différents permet, par extrapolation, de prévoir des équilibres de phases
degrés de libertés. On peut alors définir pour chaque configu- dans des systèmes d’ordre supérieur : avantage non négli-
ration de l’alliage une énergie libre intégrant les degrés de geable dans le cadre d’études de systèmes multi-éléments,
liberté rapides (excitations électroniques, vibrations…) puis puisque cela permet la rationalisation d’un plan d’expériences.
effectuer un développement en amas. On obtient ainsi un
modèle d’ISING avec des paramètres d’interaction dépendant Le choix des données expérimentales
de la température pour rendre compte des contributions entro-
piques des degrés de liberté internes. Ce modèle d’ISING peut Les données expérimentales utiles pour optimiser un système
ensuite être utilisé dans les approches classiques décrites ci- peuvent provenir de :
dessus pour prendre en compte l’entropie de configuration.
Cette approche, simple à mettre en œuvre, présente tout de • La détermination de diagrammes de phases : mesures de
même des limites, puisqu’elle suppose que les interactions températures d’invariants, de liquidus* ou de solidus*, de
effectives découlant des contributions entropiques correspon- limites de solubilité ;
dant aux degrés de liberté internes peuvent se réduire à un
nombre fini d’interactions à courte portée, ce qui n’est pas tou- • mesures thermodynamiques : activité, potentiel chimique,
jours justifiable physiquement. Pour l’entropie de vibration, il enthalpie de formation, de transformation ou de mélange ;
peut s’avérer, en fait, plus judicieux, mais également plus coû-
teux en temps de calcul, d’effectuer un développement en • mesures physiques : chaleur spécifique, magnétisme ;
amas des matrices de constantes de forces qui sont à courte
portée pour en déduire ensuite l’entropie de vibration des dif- • données calculées : enthalpies de formation ou de mélange
férentes configurations de l’alliage. par l’approche de MIEDEMA ou par calculs ab initio.
À haute température, les effets anharmoniques peuvent deve- Une fois le jeu de données expérimentales sélectionné, les
nir importants. Comme la principale source d’anharmonicité modèles thermodynamiques adaptés aux phases à modéliser
provient de la variation des fréquences de vibration atomique sont choisis.
avec le volume, l’approximation quasi harmonique s’avère
généralement très satisfaisante pour traiter ces effets anhar- Le choix des modèles thermodynamiques
moniques. Cette approximation consiste à appliquer le modèle
Les modèles thermodynamiques associés aux différentes
harmonique pour différents volumes et à minimiser ensuite
phases d’un système donné sont choisis en tenant compte de
l’énergie libre du système par rapport au volume.
leurs caractéristiques cristallographiques.
Quant aux excitations magnétiques, le modèle d’Ising peut
Il est également important de se préoccuper, à ce stade du
être complété par un modèle d’Heisenberg décrivant les inter-
travail, de la compatibilité des modèles choisis pour un sys-
actions entre les moments magnétiques des atomes ainsi
tème donné avec d’autres systèmes mettant en jeu des élé-
qu’une seconde contribution magnétique décrivant la variation
ments et des phases communs. Cela est fait dans la perspec-
de l’énergie d’un atome avec son moment magnétique. Les
tive de créer des bases de données par mutualisation de
propriétés thermodynamiques du système sont alors obte-
différents systèmes optimisés, et donc de pouvoir extrapoler
nues en échantillonnant non seulement les différentes confi-
des équilibres dans des systèmes multi-constitués.
gurations atomiques du système, mais également ses diffé-
rents états magnétiques (amplitude et direction du moment
Selon le type de phases : éléments purs, composés stoechio-
magnétique de chaque atome).
métriques et solutions de substitution, différents formalismes
mathématiques sont utilisés.

Les approches sur réseau rigide et CALPHAD s’avèrent com-
Le principe de l’optimisation
plémentaires. Elles peuvent faire appel aux mêmes données
L’optimisation consiste à ajuster les paramètres des modèles d’entrée pour construire une modélisation thermodynamique,
par une méthode de moindres carrés. Le principe est de déter- mais elles diffèrent par les théories qui les sous-tendent. En
miner les valeurs des termes d’interaction et des enthalpies outre, les modèles thermodynamiques utilisés par la méthode
libres de formation des composés à partir des données expé- CALPHAD se justifient généralement via une description ato-
rimentales. mistique.
On considère que l’optimisation d’un système est terminée

lorsqu’on a réussi à décrire toutes les données expérimen- Emmanuel CLOUET, Frédéric SOISSON
tales sélectionnées, avec un minimum de paramètres. et Caroline TOFFOLON,
Le calcul d’un diagramme de phases

Une fois qu’un système multi-constitué donné est considéré
comme optimisé, on peut procéder à des calculs d’équilibre. Références
Ces calculs consistent en la minimisation de l’énergie de [1] D. DE FONTAINE, « Configurational thermodynamics of solid solu-
Gibbs en prenant en compte toutes les phases du système. tions », Solid State Physics, 34 (1979), pp. 73-274; “Cluster approach
Le logiciel Thermocalc [4] permet de réaliser ces calculs de to order-disorder transformations in alloys”, Solid State Physics, 47
minimisation d’énergie de GIBBS. À partir de ce calcul, il est (1994), pp. 33-176.
possible de tracer les diagrammes de phases (fig. 33), mais [2] F. DUCASTELLE, Order and phase stability in alloys (North-Holland,
aussi des isoplèthes* ou des fonctions thermodynamiques 1991).
du système (Cp, H, S…) en fonction des diverses variables.
[3] P. SPENCER, « A brief history of CALPHAD », CALPHAD, 32 (2008),
Les calculs thermodynamiques permettent de prédire les équi-
pp. 1-8.
libres stables dans un système, mais également des équilibres
métastables. Ils permettent aussi d’effectuer des extrapola- [4] B. SUNDMAN et al., « The Thermo-Calc databank system », CAL-
tions dans des domaines de température et de composition PHAD, 9 (1985), pp. 153-190.
difficilement atteignables expérimentalement.
3 3
G (kJ/mol)
G (kJ/mol)
2 2
Liquide
1 1 Liquide
Solide Solide
0 0
-1 -1
-2 -2
0 0,2 0,4 0,6 0,8 1 0 0,2 0,4 0,6 0,8 1
A Fraction molaire de B B A Fraction molaire de B B
1000 1000
Température (°C)
Température (°C)
800 800
600 Liquide 600 Liquide
400 400
Solide
Solide
Solide
Solide
200 200
0 0
0 0,2 0,4 0,6 0,8 1 0 0,2 0,4 0,6 0,8 1
A Fraction molaire de B B A Fraction molaire de B B
Fig. 33. Exemple de tracé d’un diagramme de phases A-B à partir

des courbes d’enthalpies libres correspondantes à 200 et 480 °C.
52 Les diagrammes de phase

La simulation des dommages primaires

dans les matériaux en réacteur
L a compréhension du comportement des matériaux sou-

mis à des irradiations exige une approche multi-échelles à
même de transposer les résultats des simulations d’endom-
magement effectuées à l’échelle la plus locale (nanométrique)
vers les échelles microscopiques pertinentes. Les simulations
à l’échelle nanométrique, telles que les « boîtes à outils »
exposées préalablement (voir supra, pp. 43-44, « Les simula-
tions de dynamique moléculaire ») et plus avant (voir infra,
pp. 127-132, « La modélisation du dommage primaire d’irra-
diation par dynamique moléculaire »), souffrent de deux limi-
tations importantes. D’une part, elles ne peuvent, en pratique, E
être effectuées que pour une seule particule incidente et ne
sont donc pas adaptées à la description des effets induits par
des particules ayant un spectre énergétique très étendu, tel R(E)
que le spectre énergétique des neutrons dans un réacteur.
D’autre part, elles sont très coûteuses en termes de capacités Fig. 34. Décomposition d’une cascade de déplacements de taille
de calcul et ne peuvent donc pas facilement être étendues à caractéristique R(E), dont l’ordre de grandeur est de quelques
des simulations pour des particules très énergétiques, pro- centaines de nanomètres, en différentes sous-cascades (ellipses
duits de fission ou faisceaux d’ions. Nous présentons dans ce bleues) quasi indépendantes dont la taille est d’une dizaine
de nanomètres.
paragraphe deux approches mises en œuvre pour dépasser
ces limitations, d’une part en permettant une comparaison des
endommagements obtenus dans divers environnements,
réacteurs ou faisceau d’ions, et d’autre part en analysant la Cela définit un sous-volume fortement endommagé (ellipses
structure des cascades* de déplacements produites par des bleues), où seule la Dynamique Moléculaire* (DM) est apte
particules très énergétiques. à décrire le nombre et la nature des défauts formés. La dimen-
sion caractéristique de ces sous-cascades est de l’ordre de la
Des travaux récents ont permis de montrer que la cascade de dizaine de nanomètres. De plus, cette distribution de sous-
déplacements présentait un caractère fractal. En d’autres cascades est distribuée de manière non uniforme au sein de
termes, la densité spatiale de dommages produits au sein la cascade, présentant un taux de recouvrement non nul. Il
d’une cascade de déplacements est loin d’être uniforme, ce apparaît que l’énergie seuil en-deçà de laquelle on peut défi-
que décrit ce caractère fractal [1]. Ainsi, pour une particule inci- nir les sous-cascades dépend uniquement du matériau et du
dente très énergétique, l’endommagement est distribué en type de particule incidente. Sa détermination est complexe
zones de densités de défauts très hétérogènes, qu’il est alors mais fondamentale, car elle permet de faire le lien entre l’ap-
possible de définir comme des sous-cascades indépendantes. proche par collisions binaires et la DM. En effet, lorsque les
sous-cascades ne se recouvrent pas, une simple sommation
La figure 34 présente une telle schématisation de la décom- des effets de chaque sous-cascade décrite par des simula-
position d’une cascade en sous-cascades [2,3]. Une particule tions de DM est suffisante pour calculer le nombre de défauts
d’énergie E de l’ordre de la dizaine de keV va produire des dans une cascade de déplacements. L’intérêt de cette descrip-
chocs dans une zone dont la taille caractéristique R(E), cal- tion est fondamental, en ce sens qu’elle autorise la détermina-
culable par un modèle de collisions binaires, est de l’ordre de tion de l’endommagement par le seul calcul de ces sous-cas-
quelques centaines de nanomètres. Chaque particule mise en cades, nécessitant de considérer des dimensions de cellules
mouvement suit alors une trajectoire calculable en collisions de calcul réduites.
binaires (traits pleins). Quand l’énergie de la particule dépla-
cée devient inférieure à une énergie seuil, la densité de Dans un réacteur, la distribution énergétique des particules
défauts créés localement augmente du fait de la diminution du incidentes constitue une caractéristique fondamentale. La dis-
parcours moyen entre chocs successifs. tribution en énergie cinétique des neutrons dépend du réac-
teur considéré, qui est ainsi caractérisé par son spectre neu-

tronique (neutrons/cm2/s/eV) qui s’étend entre quelques milli- cule, dans le cadre du modèle de collisions nucléaires iso-
électronvolts (neutrons thermalisés) et quelques dizaines de tropes, les dommages primaires dans un réacteur, ainsi que
millions d’électronvolts (neutrons émis par fission) [4]. Dans le nombre d’atomes mis en mouvement dans la phase pri-
un réacteur, cette distribution en constitue une caractéristique maire de la cascade de déplacements. Néanmoins, il ne per-
fondamentale. Lors de la traversée du matériau par les neu- met pas d’avoir ces informations, dans le cas d’un faisceau
trons, de nombreuses collisions se produisent entre ces par- d’ions. Nous citerons également le programme MARLOWE,
ticules et les atomes constituant le matériau. À titre d’exemple, qui permet d’introduire des matériaux anisotropes et se trouve,
un neutron de 1 MeV transmet une énergie moyenne de par conséquent, plutôt adapté à l’irradiation de monocristaux
35 keV à un atome de fer dans un cristal de fer [5]. Cet atome par des faisceaux d’ions.
de fer met en mouvement environ 500 autres atomes de fer
avant de s’arrêter : ces déplacements atomiques, appelés Du fait qu’aucun de ces programmes ne permet de prendre
aussi « dommages primaires », constituent le premier stade en compte à la fois une distribution en énergie des particules
de l’endommagement du matériau [6]. incidentes et différents types de particules, le programme
DART [1] et [11] a été conçu pour permettre une réelle com-
paraison des dommages primaires résultant d’irradiations neu-
Les dommages primaires troniques ou par faisceaux de particules, ions et électrons.
Sur le plan technologique, le calcul des dommages primaires
dans un matériau permet non seulement une estimation de
l’endommagement de ce matériau dans un environnement
Les estimateurs de dommages
donné, mais au-delà une comparaison plus réaliste du com- primaires
portement des matériaux dans différents réacteurs. Ce Pour permettre cette comparaison entre diverses irradiations,
concept de dommages primaires peut ainsi être étendu à tous il apparaît nécessaire de définir différents estimateurs à même
les moyens d’endommagement sous irradiation, en particulier de représenter les endommagements.
les implantations par faisceau d’ions. On conçoit alors tout l’in-
térêt de développer des outils qui permettent de valider l’usage Le taux de production de dommages A correspond au
d’une irradiation par faisceau d’ions pour simuler une irradia- nombre d’atomes déplacés par un flux de particules incidentes
tion neutronique en réacteur. Ces irradiations présentent, en en une seconde. À partir de cet estimateur, on peut ainsi défi-
effet, d’énormes avantages relativement à celles en réacteur : nir le nombre total de dommages + après une durée d’irradia-
choix de la particule, modulation de la cinétique de l’expé- tion ] : + = A]. + est classiquement appelé le « nombre de
rience, choix de la température, faible disponibilité des réac- déplacements par atome » (dpa*). Bien que cet estimateur
teurs expérimentaux, absence de radioactivité induite [7]… soit largement utilisé pour définir l’endommagement dans un
réacteur, il n’est pas suffisant pour comparer les dommages
dans différents réacteurs ou irradiateurs. En effet, pour des
L’approximation des collisions taux de production de dommages égaux, il peut correspondre
binaires des distributions de défauts très diverses, induisant des effets
Une particule pénétrant dans un solide perd son énergie ciné- microstructuraux très différents.
tique, lors de multiples collisions. Lorsqu’une particule pos-
sède une énergie cinétique importante, de l’ordre de 1 MeV, De surcroît, il est d’usage de définir le spectre de primaires*
la distance entre deux collisions est très grande et les chocs (ou spectre de PKA, Primary Knocked-on Atoms) comme
avec les atomes du cristal sont bien distants. Ils peuvent alors étant la distribution en énergie des premiers atomes déplacés.
être décrits dans le cadre des collisions binaires. Ce spectre permet de connaître l’énergie transférée aux diffé-
rents atomes constituant le matériau.
De nombreux travaux débutés dans les années 60 permettent
de simuler cette phase de collisions binaires, à partir de Enfin, si l’on considère non seulement l’énergie des atomes
modèles universels (modèles dits « Binary Collision primaires, mais aussi la totalité des déplacements qu’ils
Approximation » ou « BCA ») [8]. induisent au cours de leur ralentissement et de celui des
atomes successivement déplacés, on peut déterminer le
Afin d’obtenir une estimation des dommages primaires, les spectre pondéré, distribution en énergie des atomes dépla-
programmes BCA s’avèrent des outils fiables et peu coûteux cés lors de la cascade ; cette distribution conditionne la
en temps de calcul. Le programme SRIM [9], par exemple, est réponse du matériau soumis à une irradiation. C’est ce
largement utilisé pour définir la profondeur de pénétration spectre qui doit être pris en compte pour sélectionner la
d’ions dans un matériau. Un tel programme ne permet pas, en masse et l’énergie de l’ion projectile simulant les mêmes
revanche, de décrire l’effet d’un flux de particules ayant une dommages primaires qu’en réacteur.
distribution en énergie, tel le flux de neutrons d’un réacteur.
Un autre programme couramment utilisé, SPECTER [10], cal-
54 La simulation des dommages primaires dans les matériaux en réacteur

Les dommages primaires par des ions Argon [14,15,16] et la microstructure obtenue
dans les ODS dans les deux cas est tout à fait similaire (dissolution partielle
des oxydes et apparition de fins précipités dans la matrice
Cette approche est illustrée ici en considérant un matériau métallique), ce qui n’est pas le cas pour une irradiation avec
développé dans le cadre des études de nouveaux réacteurs des électrons (diminution de la taille des particules d’oxydes).
à neutrons rapides. Il s’agit d’un acier contenant une disper-
sion de particules nanométriques d’oxydes dit « ODS (Oxide
Dispersion Strengthened) * » [13]. Le spectre pondéré de ce La limitation de l’approximation
matériau sous irradiation neutronique a été déterminé avec le des collisions binaires
code DART à partir de la connaissance du spectre des neu-
trons dans le réacteur, des sections efficaces de collision neu- Le ralentissement des particules dans un matériau induit une
tron-isotope, et du nombre de déplacements induits par la cas- distribution hétérogène de défauts, créant ainsi des micro-
cade. Ce spectre pondéré a ensuite été comparé à celui structures complexes. En effet, lorsque l’énergie d’une parti-
obtenu par différentes irradiations dans les ODS. Pour simu- cule diminue lors d’une collision avec un atome du milieu pour
ler la stabilité structurale des ODS dans un réacteur, des expé- atteindre quelques dizaines de keV, la distance entre deux
riences sont menées en testant ces matériaux sous irradia- chocs devient très faible et il se produit une forte concentration
tion ionique ou électronique. de déplacements atomiques localisés dans un faible volume,
appelée « sous-cascade » ou « pic de déplacement » [1] et
La figure 35 compare le spectre pondéré pour un flux de neu- [17]. Lorsque l’énergie cinétique des particules est inférieure
trons rapides, avec les spectres pondérés produits par diffé- à une centaine d’électronvolts, il se produit une quasi-fusion de
rents faisceaux d’ions (He 1 Mev, Ar 600 keV) ou par une irra- cette zone. Cette étape constitue la pointe thermique et corres-
diation électronique (électrons de 1 MeV). Ici, ce sont les ions pond à la fin de la cascade de déplacements. Ainsi, les chocs
Argon de 600 keV qui donnent le spectre pondéré le plus entre les atomes ne peuvent plus se décrire par des collisions
proche du spectre obtenu dans le réacteur.De cette analyse, binaires et échappent à tout traitement analytique. En particu-
il en ressort que la morphologie des cascades de déplace- lier, on observe la recombinaison d’un grand nombre de
ments induites par des ions Ar de 600 keV est similaire à celle défauts ou, au contraire, leur ségrégation en petits amas,
produite par un réacteur nucléaire. conduisant à des nombres effectifs de défauts très différents
de ceux estimés par les approches globales précédentes.
Des observations au microscope électronique ont été faites Seule la simulation numérique permet d’obtenir une descrip-
sur des échantillons irradiés à 30 dpa en réacteur, ainsi que tion représentative de cette étape. La DM permet ainsi de
décrire qualitativement le nombre d’atomes déplacés et de
défauts effectivement produits dans ce pic de déplacement.
Électron 1 MeV
Laurence LUNÉVILLE, Guido BALDINOZZI
Ion Ar 0,6 MeV
Ion He 1 MeV et David SIMEONE,
Neutron (FBR)
LRC CARMEN CEA/CNRS/ECP
Spectre pondéré
0,8
Références
[1] D. SIMEONE et L. LUNEVILLE, Phys. Rev. E, 82 (2010), p. 011122.
0,6
[2] S. JUMEL et J. VAN DUYSEN, Journal of Nuclear Materials, 328 (2004),
0,4 p. 151.
[3] L. LUNEVILLE, D. SIMEONE et W. WEBER, Journal of Nuclear Materials,

0,2 415 (2011), p. 55.
[4] J. BUSSAC et P. REUSS, Traité de neutronique, Hermann (1986).

0
[5] M. ROBINSON, Journal of Nuclear Materials, 216 (1994), p. 1.
101 102 103 104 105 106 107
Énergie (eV)
[6] G. MARTIN, P. BELLON, Solid State Physics, 53-54 (1997), p. 1.
Fig. 35. Évolution du spectre pondéré dans un acier ODS pour [7] L. LUNEVILLE, D. SIMEONE et Y. SERRUYS, Mater. Res. Soc. Symp. Proc.
différents projectiles : électrons d’1 MeV (orange), ions Ar 1215 (2010), V13-03.
de 600 keV (vert), ions He de 1 MeV (bleu) et pour un flux
[8] J. LINDHARD, V. NIELSEN et M. SCHRAFF, Kgl. Dan. Vid. Mat. Fys.
de neutrons typique d’un réacteur à neutrons rapides (FBR) [rouge] :
Medd., 36 (1968), p. 1.
ce sont les ions Argon de 600 keV qui donnent le spectre pondéré
le plus proche du spectre obtenu dans le réacteur.

[9] J. ZIEGLER, J. BIERSACK and U. LITTMARK, The stopping power and
range of ions in solids, Pergamon (1985).
[10] L. GREENWOOD et R. SMITHER, « SPECTER: Neutron damage cal-

culation for materials irradiations », ANL/FPP/TM (1985), p. 187.
[11] L. LUNEVILLE, D. SIMEONE et C. JOUANNE, Journal of Nuclear

Materials, 353 (2006), p. 89.
[12] D. SIMEONE et O. HABLOT, Journal of Nuclear Materials, 246 (1998),

p. 206.
[13] A. HIRATA, T. FUJA, Y. WEN, J. SCHNEBEL, C. LIU et M. CHEN, Nature

Materials, 10 (2011), p. 922.
[14] M. KLIMIANKOU, R. LINDAU, A. MÔSLANG et J. CRYST, Growth, 249

(2003), p. 381.
[15] I. MONNET, PhD (2000).
[16] C. CHEN, J. SUN et Y. XU, Journal of Nuclear Materials, 283 (2000),

p. 1011.
[17] D. SIMEONE et L. LUNEVILLE, Euro. Phys. Let., 83 (2008), p. 56002.
56 La simulation des dommages primaires dans les matériaux en réacteur

Les modèles cinétiques
L a modélisation de la cinétique d’évolution des matériaux gurations stables d’un alliage, énergie de formation et de
en vieillissement thermique ou sous irradiation peut être faite migration des défauts). De plus, le calcul des fréquences de
à différentes échelles et avec différents niveaux d’approxima- saut – en particulier, celui du terme pré-exponentiel qui intègre
tion, selon le degré de précision requis, la complexité des phé- la contribution des phonons – reste très coûteux [2]. Une
nomènes envisagés et les temps physiques simulés. Nous approche alternative consiste à adapter les méthodes d’inter-
décrivons dans cette section quatre modèles cinétiques* et actions effectives sur réseau rigide déjà présentées (voir
discutons leurs champs d’application. supra, pp. 25-28, « L’approche multi-échelles de la modélisa-
tion des matériaux »), en y intégrant la description des confi-
gurations de col. Dans les deux cas, l’ajustement sur des cal-
Monte-Carlo cinétique atomique culs ab initio peut améliorer considérablement la qualité du
(AKMC) modèle.
Dans ces simulations, on suit la position de tous les atomes du

Une fois connues les fréquences de saut, les simulations sont
système, ce qui nécessite une description précise des proprié-
réalisées avec des systèmes de quelques dizaines de milliers
tés thermodynamiques et cinétiques. À l’échelle atomique, ces
à quelques millions d’atomes. Les algorithmes dits « à temps
propriétés interviennent dans les fréquences de saut des
de résidence » sont particulièrement adaptés aux simulations
atomes, sauts qui se produisent par l’intermédiaire de défauts
des cinétiques, se faisant par l’intermédiaire de défauts ponc-
ponctuels. Le point clé est donc la description des concentra-
tuels, et permettent de définir une échelle de temps physique.
tions et des fréquences de saut des défauts ponctuels, avec,
Si les simulations sont réalisées avec un nombre de défauts
en particulier, la dépendance de ces fréquences avec l’envi-
constant, il faut cependant corriger ce temps pour tenir compte
ronnement local.
du fait que, dans un système réel, les concentrations en
défauts évoluent avec la microstructure [12].
Le calcul de toutes les fréquences de
saut par des méthodes ab initio, que ce
soit tout au long de la simulation ou par 3 nm
une tabulation préalable des fré- 50 h 50 h

quences pour tous les environnements
possibles, semble, pour l’instant, hors
100 h 100 h
de portée. Il faut donc employer des
approximations plus simples. Deux
types de méthodes sont utilisés : celles 150 h 150 h
avec relaxation des positions ato-
miques et celles sur réseau rigide.
240 h 240 h
Dans le premier cas, les fréquences de
saut des défauts sont calculées à l’aide
de potentiels interatomiques, tels que 480 h 480 h
ceux déjà décrits (voir supra, pp. 25-28,

« L’approche multi-échelles de la modé- 812 h 812 h
lisation des matériaux »), en prenant en
compte les relaxations et les vibrations
1 067 h 1 067 h
des atomes autour de positions cristal-
lines parfaites, et donc en intégrant
naturellement l’entropie de vibration et Fig. 36. Décomposition d’un alliage Fe-20 % at.Cr entre une phase
les interactions élastiques à longue por- riche en fer () et une phase riche en chrome (⬘), au cours d’un
tée. En pratique, il est toutefois très difficile de construire des vieillissement thermique à 500 °C. Seuls les atomes de Cr situés
potentiels interatomiques offrant une description satisfaisante dans la phase ⬘ sont représentés. À gauche : simulation par Monte-
Carlo cinétique atomique AKMC, à droite : observation à la sonde
de l’ensemble des propriétés pertinentes (énergie des confi-
atomique 3D [10].

Par rapport aux méthodes de champ moyen, le grand avan- sons ne sont pas suffisantes pour estimer l’ensemble des cou-
tage des simulations de Monte-Carlo* atomique est de per- plages de flux mis en jeu dans un alliage donné.
mettre un calcul exact des fréquences de saut pour chaque
environnement, dans le cadre du modèle énergétique choisi, Dès lors, la tendance actuelle est d’affronter la complexité
ainsi que la prise en compte directe des effets de corrélation d’une approche multi-échelles, d’autant plus que les progrès
entre sauts successifs des défauts ponctuels. Ce dernier point actuels des méthodes de calcul de structure électronique per-
est particulièrement important pour décrire les phénomènes mettent d’envisager une estimation complète des fréquences
de couplage qui contrôlent la ségrégation induite par l’irradia- de saut. On est capable de calculer à la fois les énergies d’ac-
tion [11]. Par ailleurs, les fluctuations thermiques sont naturel- tivation et les fréquences d’attaque des défauts ponctuels en
lement présentes, ce qui permet une bonne description des fonction de la concentration locale en atomes de soluté.
processus de germination [12]. Par rapport aux simulations de Quand le système étudié se prête à une modélisation Monte-
Monte-Carlo sur objets ou sur événements (voir ci-après), les Carlo de la diffusion sur un réseau rigide, on peut mesurer
simulations de Monte-Carlo atomique ne requièrent pas la les coefficients Lij en calculant le parcours quadratique moyen
définition d’objets prédéfinis : les éléments de la microstruc- des particules à l’équilibre. Parallèlement à cette approche
ture et leur évolution sont le résultat des sauts individuels des numérique, des théories de champ moyen ont été dévelop-
défauts ponctuels. La contrepartie est un temps de calcul plus pées pour comprendre les effets du mécanisme de diffusion
important, qui limite l’emploi des méthodes AKMC aux pre- sur les flux atomiques (voir [1] pour une revue sur les modèles
miers stades des transformations de phases (fig. 36). de diffusion). Le calcul en champ moyen des coefficients Lij
n’est cependant pas si simple, cela pour plusieurs raisons.
Dans un alliage, la diffusion est assurée par une tierce
Les équations de diffusion espèce, i.e. le défaut ponctuel. Un alliage contient donc, au
minimum, trois espèces chimiques, deux espèces atomiques
en champ moyen
et le défaut ponctuel. Ce ménage à trois exige du défaut ponc-
La complexité des mécanismes de diffusion, et notamment les tuel, quand il s’échange avec un atome, un choix entre un
effets de corrélation, sont bien traités en AKMC. Cependant, atome de l’une ou l’autre espèce chimique. Cette compétition
comme tous les sauts successifs doivent être considérés, le induit de fait des séquences de saut qui s’écartent du chemin
temps de calcul peut être long, et les temps physiques envi- aléatoire, ce qui constitue les effets de corrélation.
sagés restent limités. La modélisation des cinétiques à une L’estimation de ces séquences, et en particulier les écarts au
échelle plus grossière que l’échelle atomique nécessite de chemin aléatoire qui s’accompagnent d’un couplage entre les
savoir comment déduire une équation de diffusion macrosco- séquences de saut des différentes populations, se complique
pique pour les atomes, à partir des fréquences de saut d’un d’autant plus que l’ordre à courte distance est grand.
défaut ponctuel. Rappelons que l’ordre à courte distance est caractérisé par
le fait que la probabilité d’occupation d’un site dépend de l’oc-
Dans les alliages concentrés, cette approche multi-échelles cupation des sites voisins. Il existe principalement deux théo-
est plus difficile à établir que dans le cas dilué, cela en raison ries qui incluent les effets de l’ordre à courte distance sur le
du grand nombre de fréquences de saut mises en jeu du fait calcul des coefficients Lij, la Path Probability Method (PPM)
de la variation du nombre d’atomes de soluté possibles à [8], et celle que nous avons développée, le champ moyen
proximité du défaut ponctuel. Les premiers modèles de diffu- auto-cohérent (SCMF) [9].
sion se sont appuyés sur la Thermodynamique des Processus
Irréversibles (TPI), qui relie les flux macroscopiques des Dans cette dernière approche, la variation temporelle des
atomes aux forces thermodynamiques par une forme généra- moments de la fonction de distribution est dérivée de l’équa-
lisée de la loi de Fick, en introduisant des coefficients de dif- tion maîtresse microscopique. Le point original de cette théo-
fusion phénoménologiques. Ces derniers, appelés également rie est d’expliciter la correction à la fonction de distribution
« coefficients d’Onsager » ou « coefficients Lij », représentent d’équilibre employée pour calculer les moments associés et
les couplages cinétiques entre les flux. La théorie de la TPI de fournir une méthode auto-cohérente pour calculer cette cor-
est couramment utilisée pour interpréter les évolutions conti- rection. Cette théorie a été appliquée avec succès à l’étude
nues d’un champ de concentration, comme celles apparais- de l’effet de l’ordre à courte distance et a été étendue, par la
sant dans une expérience de diffusion, ou encore les phéno- suite, au mécanisme de diffusion par interstitiel dissocié, cela
mènes de couplage de flux, qui sont très importants sous dans différentes structures cristallographiques.
irradiation. Sur le principe, un ensemble d’expériences de dif-
fusion bien choisies, combiné à des mesures thermodyna- Au-delà du calcul des coefficients Lij dans l’état stationnaire,
miques de type calorimétrie, doivent permettre de déterminer les méthodes de champ moyen peuvent aussi être employées
la variation des forces motrices et des coefficients Lij avec la pour simuler des cinétiques complexes sous irradiation [9]. Ce
température et la composition de l’alliage. Cependant, dans la type de modélisation, qui s’appuie sur des équations de diffu-
pratique, ces expériences de diffusion ne peuvent être réali- sion déterministes, permet une description fine des couplages
sées qu’à haute température, et les données dont nous dispo- de flux et de leur effet sur les variations de composition d’un
58 Les modèles cinétiques

20 20 20 20
a b c d
10 10 10 10
0 0 0 0
-10 -10 -10 -10
-20 -20 -20 -20

-20 -10 0 10 20 -20 -10 0 10 20 -20 -10 0 10 20 -20 -10 0 10 20
Fig. 37. Autour d’une dislocation dans l’alliage Ni-Si, lignes de flux (a) des lacunes, des atomes de silicium aux températures (b) T = 1 010 K,
(c) T = 1 060 K, et (d) T = 1 110 K. La distribution initiale des atomes Si et des lacunes est homogène. Les dimensions sont données en unité de
vecteur de Burgers (0,25 nm) pour une dislocation coin dans le nickel placée à l’origine, avec un axe horizontal selon [11^0] et un axe vertical selon
[111] dans le plan (112^). Les couleurs indiquent l’amplitude du cosinus de l’angle entre la direction du flux et le vecteur radial de la dislocation ⃗,
avec en orange les flux qui sortent du cœur de la dislocation, en bleu ceux dirigés vers le cœur et en blanc ceux orthogonaux à ⃗ [13].
alliage à proximité des puits de défauts ponctuels tels que les est identifié par son nombre d’atomes de chaque espèce
joints de grains, les surfaces ou les dislocations (fig. 37). minoritaire et son nombre de défauts ponctuels, et une forme
Toutefois, cette approche ne permet pas d’étudier les sépara- lui est attribuée. Les différentes configurations atomiques des
tions de phases de type germination-croissance susceptibles amas sont prises en compte dans les termes d’entropie de
d’avoir lieu à proximité de ces mêmes puits de défauts. configuration de l’énergie libre de formation [5], qui sont
cependant souvent négligés, en pratique.
Monte-Carlo cinétique sur objets

Comme pour les simulations AKMC, l’évolution du système se
(OKMC) et sur événements (EKMC) fait en choisissant à chaque instant le prochain événement,
Pour éviter de traiter tous les sauts atomiques successifs, dont la fréquence d’apparition est connue, selon un algorithme
comme c’est le cas en AKMC, une approximation consiste à à temps de résidence. Dans l’approche OKMC, un événement
considérer la cinétique comme une cinétique d’évolution peut être le saut d’un amas sur une distance interatomique, la
d’amas constitués d’espèces minoritaires et/ou de défauts dissociation d’un amas, ou la production de défauts et d’amas
ponctuels créés sous irradiation, sans
traiter explicitement l’espèce majori-
taire. Cela se révèle bien adapté au cas
d’un matériau idéal sans impureté ou
d’un alliage dilué. Le système est donc
vu comme un « gaz d’amas » (fig. 38) :
les amas diffusent dans la matrice
constituée de l’espèce majoritaire,
réagissent entre eux et peuvent se dis-
socier thermiquement. Les méca-
nismes élémentaires de diffusion ne
sont pris en compte que de manière
effective, par les valeurs des coeffi-
cients de diffusion attribués aux diffé-
rents amas. Outre les coefficients de
diffusion effectifs, qui peuvent être 77 K 152 K
déterminés par des calculs AKMC, les
autres paramètres de ce modèle sont
Fig. 38. Modélisation en Monte-Carlo cinétique sur événements d’un échantillon de fer irradié
principalement les énergies libres de aux électrons à 4 K et soumis à un recuit isochrone [4]. Le système est représenté à deux
formation des amas. Ces paramètres températures ; les amas interstitiels sont représentés par des tores (boucles de dislocation),
peuvent être déduits de calculs ab ini- les cavités par des cubes. Les atomes de fer situés sur les sites du réseau cubique centré
tio ou à l’aide de potentiels semi-empi- ne sont pas considérés dans la simulation. À basse température, les défauts sont
essentiellement présents sous la forme de paires de Frenkel (sphères bleues et jaunes,
riques pour les petits amas, et de lois
pour les mono-lacunes et mono-interstitiels respectivement), alors qu’à plus haute
analytiques pour les amas de plus température, la migration des interstitiels entraîne la formation d’amas d’interstitiels
grande taille. Généralement, un amas et une diminution du nombre de défauts due aux recombinaisons avec les lacunes.

de défauts sous irradiation. La migration des amas par sauts
AKMC
atomiques successifs pouvant pénaliser lourdement les cal- 1019
DA
d (cm-3)
culs, une approximation supplémentaire consiste à remplacer
les événements de migration par les conséquences de ces 1018
migrations, c’est-à-dire des événements de réaction entre

1017
amas ou d’absorption d’amas par les puits (surfaces, joints de
100 101 102 103 104 105 106 107 108
grains, etc.). Cette approche est connue sous le nom de 4
R (nm)
Monte-Carlo cinétique sur événements (EKMC). Elle a été uti- 3
lisée notamment pour la précipitation homogène en vieillisse- 2
ment thermique [5] et pour les recuits de résistivité ([4] et infra, 1
0
pp. 137-140, « Structure et cinétique des défauts d’irradiation
100 101 102 103 104 105 106 107 108
dans le fer »).
C1
0,012
0,008
La dynamique d’amas 0,004
0
De même que les simulations EKMC, la dynamique d’amas* 100 101 102 103 104 105 106 107 108
t (s)
est une approche de type « gaz d’amas ». Les paramètres
d’entrée sont très proches : ce sont principalement les éner-
Fig. 39. Cinétique de précipitation à 500 °C dans un alliage Fe-
gies libres de formation des amas et leurs coefficients de dif- 1.34 %Cu, simulée par dynamique d’amas (DA) et par Monte-Carlo
fusion. À la différence des simulations EKMC, le système est cinétique atomique (AKMC), en utilisant des paramètres avec
traité en champ moyen. Ainsi, les variables du problème sont lesquels les amas sont immobiles. De haut en bas : densité d’amas,
les concentrations des différents types d’amas, qui sont sup- rayon moyen, concentration de monomères de Cu dans la matrice.
Les différents stades (germination, croissance, coalescence) sont
posés être uniformément répartis dans l’espace. Cette hypo-
clairement visibles. La dynamique d’amas permet de traiter le stade
thèse n’est pas toujours valide ; par exemple, les corrélations de coalescence sur des temps longs.
spatiales entre défauts créés par irradiation sont souvent
fortes, ce qui peut rendre l’utilisation de cette méthode plus
délicate. Dans certains cas, comme pour la simulation d’irra-
diations aux ions qui induisent un dommage fortement dépen-
dant de la profondeur, il est utile de s’affranchir de l’hypothèse
d’homogénéité sur tout le système. On considère alors un sys- alliages de Zr »), ou encore la précipitation homogène du
tème découpé en tranches, couplées entre elles par la diffu- cuivre dans le fer (fig. 39, [6] et infra, pp. 157-160, « L’évolution
sion des espèces mobiles, le système étant supposé homo- microstructurale dans les alliages modèles d’aciers de cuve »).
gène dans chaque tranche.
L’évolution du système est obtenue en résolvant un jeu d’équa- Thomas JOURDAN, Frédéric SOISSON
tions de cinétique chimique sur les concentrations d’amas. La et Maylise NASTAR,
résolution déterministe du système d’équations permet, dans
de nombreux cas, de simuler l’évolution sur des temps phy-
siques longs et d’avoir accès aux fortes doses (plusieurs
dizaines de dpa). Toutes les classes d’amas étant simulées,
Références
du monomère aux amas contenant quelques millions de
[1] A. R. ALLNATT et A. B. LIDIARD, « Atomic Transport in Solids »,
défauts ou de solutés, les processus de germination, crois-
Cambridge University Press (1993).
sance et coalescence d’amas sont directement pris en
compte, comme dans les approches AKMC et EKMC. Les dis- [2] J.-L. BOCQUET, « Contribution of shear strains to the vibrational
tributions des tailles d’amas (boucles interstitielles, cavités, entropy of defect configurations », Philos. Mag., 87 (2007), pp. 3259-
3295.
précipités) obtenues peuvent être directement comparées aux
observations expérimentales, par exemple en microscopie [3] F. CHRISTIEN et A. BARBU, « Cluster dynamics modelling of irradia-
électronique à transmission. tion growth of zirconium single crystals », Journal of Nuclear Materials,
393 (2009), pp. 153-161.
À titre d’exemple, la dynamique d’amas a été utilisée pour [4] C.-C. FU, J. DALLA TORRE, F. WILLAIME, J.-L. BOCQUET et A. BARBU,
simuler les recuits de résistivité dans le fer en présence de « Multiscale modelling of defect kinetics in irradiated iron », Nature
carbone ([7] et infra, pp. 137-140, « Structure et cinétique des Materials, 4 (2005), pp. 68-74.
défauts d’irradiation dans le fer »), la formation d’amas de [5] T. JOURDAN, J.-L. BOCQUET et F. SOISSON, « Modeling homogeneous
défauts ponctuels sous irradiation dans le zirconium ([3] et precipitation with an event-based Monte Carlo method : Application to
infra, pp. 147-150, « L’évolution de la microstructure des the case of Fe-Cu », Acta Mater., 58 (2010), pp. 3295-3302.
60 Les modèles cinétiques

[6] T. JOURDAN, F. SOISSON, E. CLOUET et A. BARBU, « Influence of clus-
ter mobility on Cu precipitation in α-Fe : A cluster dynamics modeling »,
Acta Mater., 58 (2010), pp. 3400-3405.
[7] T. JOURDAN, C.-C. FU, L. JOLY, J.-L. BOCQUET, M. J. CATURLA et

F. WILLAIME, « Direct simulation of resistivity recovery experiments in
carbon-doped alpha-iron », Physica Scripta, T145 (2011), p. 014049.
[8] R. KIKUCHI, « Statistical dynamics of crystalline diffusion », J. Phys.

Chem. Sol., 20 (1961), p. 17.
[9] M. NASTAR, « Diffusion and coupled fluxes in concentrated alloys

under irradiation: a self-consistent mean-field approach », C. R.
Physique, 9 (2008), pp. 362-369.
[10] S. NOVY, P. PAREIGE et C. PAREIGE, « Atomic scale analysis and

phase separation understanding in a thermally aged Fe-20 at.%Cr
Alloy », Journal of Nuclear Materials, 384 (2009), pp. 96-102.
[11] F. SOISSON, « Kinetic Monte-Carlo simulations of radiation induced

segregation and precipitation », Journal of Nuclear Materials, 349
(2006), pp. 235-250.
[12] F. SOISSON et C.-C. FU, « Cu-precipitation kinetics in α-Fe from ato-

mistic simulations: Vacancy-trapping effect and Cu-cluster mobility »
Phys. Rev. B, 76 (2007), p. 214102.
[13] T. GARNIER, V. MANGA, D. R. TRINKLE, M. NASTAR et P. BELLON,

« Stress-induced anisotropic diffusion in alloys: Complex Si solute flow
near a dislocation core in Ni », Phys. Rev. B, 88 (2013), p. 134108.

La dynamique des dislocations
L a dynamique des dislocations* (DD) a pour objet l’étude Le premier code de dynamique des dislocations 3D a été
des mécanismes de déformation plastique à l’échelle du grain développé en France au début des années 90 [1]. Depuis cette
(quelques µm) par une modélisation du comportement de l’en- date, plusieurs codes ont été développés en France et à
semble des dislocations* sous-jacentes19 (voir l’encadré l’étranger. Parmi les plus significatifs, citons les codes MICRO-
supra, p. 20, intitulé « Les dislocations »). Le mouvement de MEGAS (CNRS-ONERA) et TRIDIS (CNRS-Grenoble INP)
ces défauts cristallins linéiques, leurs interactions mutuelles pour la France, les codes PARANOID (IBM T.J. Watson),
ou avec les autres éléments microstructuraux (joints de grains, Parametric Dislocation Dynamics (UCLA) et PARADIS (LLNL,
précipités…) sont, en effet, responsables de nombreuses pro- UCLA) pour les États-Unis, le code développé au Karlsruher
priétés mécaniques du matériau (écrouissage, tenue à la Institute für Technologie en Allemagne ainsi que le code
fatigue, ténacité…). NUMODIS développé conjointement par le CEA, le CNRS et
l’INRIA.
La DD s’inscrit par nature dans une démarche multi-échelle
et vise à faire le lien entre la dynamique moléculaire (voir
supra, pp. 43-44, « Les simulations de dynamique molécu- Le principe de la dynamique
laire ») et la plasticité cristalline (voir infra, pp. 71-74, « La des dislocations
viscoplasticité cristalline »). Intégrant les mécanismes phy-
La modélisation par dynamique des dislocations permet de
siques élémentaires étudiés par la première méthode (voir
suivre l’évolution d’un ensemble de dislocations, représentées
infra, pp. 163-167, « Le calcul ab initio de la plasticité : struc-
individuellement par la ligne de discontinuité du champ de
ture de cœur et mécanismes de glissement des dislocations
déplacement, au cours du temps et sous une sollicitation
vis ». Voir également, infra, pp. 169-170, « La modélisation de
mécanique donnée [2]. Chaque ligne de dislocation est alors
la mobilité des dislocations »), elle a pour objectif de fournir
discrétisée par une succession de segments caractérisés par
des lois de comportement physiquement motivées à la
leur vecteur de Burgers*, leur plan de glissement voire de
deuxième. La DD permet, en outre, l’étude des mécanismes
montée (fig. 40).
d’endommagement lorsque ceux-ci se jouent à l’échelle du
grain (voir infra, pp. 171-173, « La fatigue des aciers austéni-
Cette représentation conduit, par comparaison avec la dyna-
tiques inoydables »).
mique moléculaire, à une réduction drastique du nombre de
19. Les densités typiques sont de l’ordre de 1012 à 1015 m/m3, soit une lon- degrés de liberté nécessaire au suivi d’une dislocation indivi-
gueur de dislocation de l’ordre de 1 mm à 1 m pour un grain de taille 10 µm.
a b c
Fig. 40. a) Une dislocation « coin ». La présence de cette dislocation induit une distorsion du réseau cristallin caractérisée par la ligne
de dislocation (vecteur normal sortant) et par son vecteur de Burgers b. Le plan de glissement (pointillé bleu) de cette dislocation contient
à la fois la ligne de dislocation et son vecteur de Burgers. La figure (b) présente la discrétisation en segments « vis » et « coin » utilisée par
le code TRIDIS pour représenter une ligne de dislocation courbe (noir pointillé). La figure (c) présente la discrétisation de cette même dislocation
par une succession de segments d’orientations quelconques, telle qu’elle est utilisée par le code NUMODIS.

duelle, ce qui permet de simuler des domaines plus vastes sur
La vitesse des dislocations
des temps plus longs. La contrepartie de cette simplification
est de devoir introduire explicitement les règles de comporte- La vitesse de déplacement d’un élément de dislocation s’ex-
ment de ces segments individuels (interactions, mobilité…). prime en fonction de la force de Peach-Koehler qui s’exerce
sur lui à travers une loi de mobilité 2 :
Une fois cette représentation choisie, un calcul de DD consiste
à intégrer au cours du temps le mouvement de la population g Ah = 2(l)
de dislocations initiales. Le déroulement d’un calcul varie d’un
code de DD à l’autre, en fonction des modèles physiques et Cette loi de mobilité contient une part importante de la phy-
des méthodes de résolution numérique choisis. Il est toutefois sique du matériau et peut être déterminée expérimentalement
possible de distinguer les trois étapes suivantes qui se succè- ou par la modélisation (voir infra, pp. 163-167, « Le calcul ab
dent à chaque pas de temps. initio de la plasticité : structure de cœur et mécanismes de glis-
sement des dislocations vis ». Voir également, infra, pp. 169-
170, « La modélisation de la mobilité des dislocations »).
Les forces exercées sur les dislocations
Les interactions entre segments de dislocations résultent de
Les collisions et recombinaisons des dislocations
la distorsion du réseau cristallin induite par leur présence [3].
Elles peuvent être calculées dans le cadre de la théorie élas- Au gré de leur mouvement au sein d’un même grain, les dis-
tique linéaire sur la base de travaux théoriques réalisés dans locations peuvent se multiplier comme dans le cas du méca-
les années 60-70 [4], dont de nombreuses études [5] ont mon- nisme de Frank-Read, s’annihiler ou se rencontrer entre elles
tré la pertinence, même aux échelles très fines. Le champ de formant ainsi des jonctions (fig. 41a-b). Elles peuvent égale-
contrainte en tout point ` de la simulation peut ainsi s’exprimer ment rencontrer des éléments microstructuraux comme des
comme la somme du champ de contrainte extérieure appli- joints de grains ou des précipités, qui vont alors servir d’obs-
quée Sa ] et des champs de contrainte dus à chaque segment tacle à leur mouvement (fig. 41c). Ces mécanismes d’interac-
de dislocation b : tions doivent être introduits explicitement dans la DD et confè-
segments rent aux dislocations un caractère particulièrement
S(`) = Sa ] (`) + f Sb(`) dynamique.
b
La force par unité de longueur gAh ressentie par chaque élé-

Le code NUMODIS
ment de dislocation est fonction de cette contrainte locale, de S’appuyant sur une collaboration de longue date autour du
son vecteur de Burgers et de sa direction donnée par le vec- code TRIDIS (voir infra, pp. 171-173, « La fatigue des aciers
teur unitaire ξ. Elle est donnée par la formule de Peach- austénitiques inoydables »), la DEN, le CNRS-Paris Est-
Koehler : Créteil, le CNRS-Grenoble INP et l’INRIA de Bordeaux ont
entrepris, depuis 2007, le développement d’un code de
gAh = (S ⋅ ) × ξ. deuxième génération appelé NUMODIS dont la première ver-
a b c
Fig. 41. a-b) Formation d’une jonction de LOMER-COTTRELL dans un matériau de structure cubique à faces centrées : le croisement entre deux
dislocations (a) situés dans deux plans de glissement distincts conduit à la formation d’une jonction située (b) dans la direction commune aux
deux plans. (c) Simulation de la plasticité dans un grain de forme cubique (taille 1 µm) contenant initialement quatre sources de FRANK-READ
de taille 0,1 µm situées au centre du grain. On note que les dislocations interagissent avec les joints du grain en venant s’y empiler.
64 La dynamique des dislocations

100
Contrainte de cisaillement (MPa)

80
60
40
20
0
0 0,2 0,4 0,6 0,8 1 1,2
Déformation m (%)
a b
Fig. 42. Simulation de l’interaction d’une dislocation « coin » avec une population de boucles induites par l’irradiation (boucles de type <100>,
densité 6.1022 m-3, taille 2,6 nm). a) La dislocation apparait en bleu et les défauts d’irradiation en rouge. b) courbe contrainte-déformation pour
un cisaillement effectué dans le plan de glissement de la dislocation « coin ». Cette courbe permet de mesurer le durcissement induit par la
présence des boucles et ainsi de faire le lien avec les lois de plasticité cristalline (voir infra, pp. 69-72, « La viscoplasticité cristalline »).
sion20 a été livrée en 2009. Son objectif est de répondre à l’en- Références
semble des besoins en modélisation par DD de la DEN et du [1] L.P. KUBIN, G.R. CANOVA, M. CONDAT, V. PONTIKIS et Y. BRECHET, Sol.
CNRS, en intégrant, d’une part, les spécificités des matériaux Stat. Phen., 23-24 (1992), pp. 455-472.
du nucléaire et en exploitant, d’autre part, le calcul haute per-
[2] V.V. BULATOV, et W. CAI, Computer simulations of dislocations,
formance, afin de permettre la modélisation de volumes et de
Oxford University Press (2006).
temps suffisamment significatifs. Un exemple de résultat
obtenu avec le code NUMODIS est donné dans la figure 42. [3] J. FRIEDEL, Dislocations, Pergamon, Oxford (1964).
[4] J.P. HIRTH et J. LOTHE, Theory of dislocations, Wiley Interscience,

Parmi les spécificités de ce code figurent notamment (a) la Mc-Graw-Hill, New York (1982).
gestion des toutes les structures cristallographiques d’intérêt
[5] M.J. HŸTCH, J.L. PUTAUX, et J.M. PENISSON, Nature, 423 (2003),
(structure cubique centrée pour l’acier de cuve, structure
pp. 270-273.
cubique à faces centrées pour les aciers austénitiques et
structure hexagonale pour le zirconium (voir infra, pp. 147-150, [6] J. DROUET, L. DUPUY, F. MOMPIOU, S. PERUSIN et A. AMBARD, Journal
« L’évolution de la microstructure des alliages de Zr »), (b) le of Nuclear Materials, 449 (2014), pp. 252-262.
traitement explicite des défauts induits par l’irradiation [6] et [7] X.J. SHI, L. DUPUY, B. DEVINCRE, D. TERENTYEV et L. VINCENT, Journal
[7] comme les boucles de Frank ou les tétraèdres de faute of Nuclear Materials, 460 (2015), pp. 37-43.
d’empilement via la prise en compte de leurs vecteurs de
[8] B. BAKO, E. CLOUET, L.M. DUPUY et M. BLÉTRY, Philosophical
Burgers et la gestion des fautes d’empilement, (c) la possibi-
Magazine, 91 (2011), pp. 3173-3191.
lité d’intégrer des lois de mobilité de glissement ou de mon-
tée [8] et (d) son intégration dans la plateforme de simulation
MATIX, ce qui lui permet d’utiliser des microstructures réa-
listes générées par le module MICROGEN ou d’être couplé
au code CAST3M pour permettre des géométries et des sol-
licitations complexes.
Laurent DUPUY,
20. La version 2.0 actuelle a été livrée fin 2014.

La production de microstructures
de matériaux pour la simulation
L es simulations numériques sur microstructures* sont uti- tion cristalline des grains présents à la surface d’un polycris-
lisées dans le but d’évaluer le comportement homogénéisé du tal (fig. 43a) [1] et la tomographie X, qui conduit à une image
matériau (voir infra, pp. 75-79, « Les méthodes d’homogénéi- 3D du coefficient d’absorption des rayons X permettant de
sation en mécanique des milieux continus ») et/ou d’étudier séparer certains constituants ou des porosités au sein d’un
en détail les répartitions des contraintes, des déformations ou matériau (fig. 43b) [2].
de l’endommagement à une échelle locale. Ces calculs s’ap-
puient sur une représentation de la microstructure, sous la
forme d’une ou de plusieurs cellules élémentaires. Ces cel- La génération de microstructures
lules élémentaires peuvent être directement issues de moyens virtuelles
d’imagerie expérimentaux 2D (microscopie optique, électro-
Les outils de génération de microstructures virtuelles présen-
nique, EBSD…) ou 3D (tomographie X). Ces cellules élémen-
tent l’intérêt de pouvoir modifier simplement certaines carac-
taires peuvent également être simulées, on parle alors de
téristiques de la microstructure, mais aussi, dans le cas de
microstructures virtuelles. Dans ce cas, la cellule simulée doit
microstructures aléatoires, de pouvoir modifier la taille des cel-
être la plus représentative possible de la microstructure étu-
lules élémentaires et multiplier les réalisations de configura-
diée.
tions géométriques différentes mais équivalentes d’un point
de vue statistique. Grâce à ces outils, il est alors possible de
définir la taille de la cellule simulée la plus représentative du
Les microstructures réelles matériau étudié (Volume Élémentaire Représentatif). Par
On cherche ici à utiliser la microstructure réelle des matériaux ailleurs, pour les matériaux polycristallins, l’EBSD 3D est
à partir d’images, discrétisées en pixels (ou voxels, pour des encore trop peu répandu pour être utilisé couramment et les
images 3D), selon une grille régulière orthogonale, obtenues microstructures polycristallines sont alors classiquement
par différentes techniques de microscopie expérimentale. Ces représentées par un découpage en cellules de VORONOI
images numériques peuvent être utilisées directement au sein (fig. 44a) [3].
d’un code de calcul mécanique, généralement après un trai-
tement numérique (filtrage, seuillage, etc.). Parmi les tech- Différentes équipes du CEA ont travaillé autour de cette ques-
niques expérimentales d’imagerie largement utilisées à cet tion et un travail de capitalisation est en cours sous la forme
effet, nous citerons l’EBSD* (Electron Back Scattering d’un module (MICROGEN) de la plateforme logicielle MATIX
Diffraction), qui permet d’obtenir une image 2D de l’orienta- (A Multiscale Materials Modeling and Simulation Platform).
a b
Fig. 43. a) Maillage réalisé à partir d’une « image » EBSD d’une éprouvette en alliage de Zirconium à gros grains. b) Maillage cubique régulier
réalisé à partir d’une image 3D obtenue par tomographie X d’un matériau composite SiC/SiC (en rouge, le matériau, en bleu la porosité).

a b c
Fig. 44. Exemples de maillage d’une microstructure polycristalline (a) d’une microstructure biphasée avec croissance de phase aux joints
de grains (b) et d’une microstructure de type « béton » constituée d’inclusions sphériques et prismatiques distribuées aléatoirement au sein
d’une matrice.
Ce découpage en cellules de Voronoi a également été utilisé (voir la description des composites SiC/SiC, infra, pp. 181-183,
pour produire des microstructures représentatives d’alliages « Le comportement macroscopique des aciers déduit de la
de zirconium polycristallins avec croissance de phase aux plasticité cristalline »), microstructures réelles et virtuelles ont
joints de grains (fig. 44b) [4]. On notera également d’autres été confrontées sur la base de différents descripteurs statis-
types de développements réalisés pour introduire différentes tiques. Les figures 45a et 45b, montrent le bon accord quali-
formes d’inclusions afin de reproduire des microstructures tatif, tandis que la figure 45c révèle l’accord quantitatif, ici
représentatives de béton (fig. 44c) [5]. observé sur la fonction de distribution radiale des centres de
fibres.
Si la comparaison des microstructures réelles et simulées
s’est limitée à un accord qualitatif pour les microstructures pré- À l’échelle supérieure, des outils de génération de microstruc-
sentées ci-dessus, une confrontation plus quantitative permet tures tissées ont également été développés pour les compo-
d’asseoir la méthode de changement d’échelles sur une base sites SiC/SiC. À nouveau, les microstructures simulées ont pu
solide. Une telle démarche a notamment été mise en œuvre être comparées aux microstructures expérimentales (fig. 46)
sur les matériaux composites SiC/SiC [6] à l’échelle du toron observées en 3D par tomographie X.
3
Fonction g
Poisson
2,5 Observations
Simulations
1,5
0,5
0
0 10 20 30 40 50
a b c Distance r (µm)
Fig. 45. Exemples de cellules élémentaires réelles (image MEB) (a) et simulées (b) et confrontation des microstructures par la fonction
de distribution radiale des centres de fibres (c) d’un matériau composite SiC/SiC.
68 La production de microstructures de matériaux pour la simulation

0,9 Simulation
Porosité
Tube 1B
0,8 Tube B2
0,7
0,6
0,5
0,4
0,3
0,2
0,1
0
0 4,1 4,2 4,3 4,4 4,5 4,6 4,7
Rayon (mm)
a b c
Fig. 46. Image d’un tube SiC/SiC obtenue par tomographie X (a). Cellule élémentaire virtuelle (b). Comparaison des évolutions de fraction
surfacique de porosité en fonction du rayon pour deux tomographies de tubes et pour la cellule simulée (c).
Nous disposons désormais d’outils pour simuler de manière Références

réaliste la microstructure de matériaux polycristallins et com- [1] C. LEBON, « Étude expérimentale et simulation numérique des
posites, donnée d’entrée des méthodes d’homogénéisation mécanismes de plasticité dans les alliages de zirconium », Thèse,
numériques (voir infra, pp. 75-79, « Les méthodes d’homogé- Université de La Rochelle (2011).
néisation en mécanique des milieux continus »).
[2] L. GÉLÉBART et al., « X-ray tomographic characterization of the
macroscopic porosity of CVI SiC/SiC composites: effects on the elas-
tic behavior », International Journal of Advanced Ceramic Technology,
Lionel GÉLÉBART, Fabien ONIMUS 7 (2010), pp. 348-360.
et Ludovic VINCENT,
Département des matériaux pour le nucléaire [3] J. PACULL, « Modèle numérique micromécanique d’agrégat polycris-
tallin pour les cavités sous pression dans les combustibles oxydes »,
Thèse, Université de Provence, Aix-Marseille I (2011).
[4] G. TRÉGO, « Comportement en fluage à haute température dans le

domaine biphasé (a + b) de l’alliage M5® », Thèse, École des Mines
de Paris (2011).
[5] B. BARY et al., « Numerical and analytical effective elastic proper-

ties of degraded cement pastes », Cement and Concrete Research, 39
(2009), pp. 902-912.
[6] C. CHATEAU, « Analyse expérimentale et modélisation microméca-

nique du comportement élastique et de l’endommagement de compo-
sites SiC/SiC unidirectionnels », Thèse, École Polytechnique (2011).

La viscoplasticité cristalline
L a viscoplasticité* cristalline s’inscrit dans le cadre de la • des lois d’évolution des différentes variables internes du
mécanique des milieux continus, à l’échelle du cristal. Elle per- modèle en fonction des incréments de glissement sur
met de traduire au moyen de variables continues les phéno- chaque système.
mènes discrets que sont le glissement et la multiplication/anni-
hilation des dislocations, ce qui constitue, en un sens, la Pour faciliter le transfert entre la mécanique des dislocations
première étape d’un changement d’échelle plus vaste entre le discrètes (voir supra, pp. 63-65, « La dynamique des disloca-
comportement collectif des dislocations dans un grain et le tions ») et la mécanique des milieux continus cristallins, la plu-
comportement élasto-visco-plastique macroscopique du part des lois de comportement cristallines actuelles intègrent
matériau. Les différentes techniques qui permettent de pas- comme variables internes les densités de dislocations par sys-
ser du comportement d’un grain décrit par la plasticité cristal- tème s, dont l’évolution suit une loi du type [3] :
line à celui du polycristal font l’objet d’un chapitre spécifique
|m|s| 1
|s = —– } — – ~s
(voir infra, pp. 75-79, « Les méthodes d’homogénéisation en
s
(1)
mécanique des milieux continus »).
De nombreuses lois de viscoplasticité cristalline existent selon avec la norme du vecteur de Burgers, s le libre parcours
le type de matériau et de comportement que l’on souhaite moyen ou distance parcourue par un segment de dislocation
représenter : monotone ou cyclique, visqueux, thermiquement avant qu’il ne soit ancré par interaction avec la microstructure
activé ou non. Au-delà de ces différences, les points suivants et ~ une distance maximale d’annihilation entre deux disloca-
restent communs aux lois de viscoplasticité cristalline : tions de signe opposé. Le libre parcours moyen dépend des
densités de dislocations présentes sur tous les systèmes ren-
• Une partition de la déformation plastique en glissements contrés par le système actif, ainsi que de la taille de grain.
plastiques sur les différents systèmes de glissement du
matériau : εnṗ = ∑b r q s = –₂¹ (v4s⨂ y
q sṁs, où r xs + y
x s⨂v4s ) est la L’application aux matériaux cristallins [5] du cadre des grandes
matrice de Taylor permettant d’orienter chaque système de transformations de la mécanique des milieux continus permet
glissement dans l’espace à partir de la normale au plan de
glissement 4 v s et de la direction de glissement, y x s.ṁs étant
enfin de reproduire la rotation cristalline qui apparaît pour
accommoder les forts niveaux de déformation plastique [4].
l’incrément de glissement viscoplastique sur le système de
glissement - ;
Quelques exemples d’application
• le calcul d’une cission résolue* s sur chaque système de de la viscoplasticité cristalline
glissement sous forme de produit contracté entre la matrice
de Taylor rs et le tenseur de contraintes local S
n:
aux matériaux du nucléaire
s = r
q s: S
n; Les matériaux de structure des centrales nucléaires sont sou-
mis à divers chargements ou dommages dont certains méca-
• le calcul d’une cission critique* {s qu’il faut dépasser pour nismes peuvent être avantageusement décrits à l’échelle de la
activer le glissement d’un système à une température don- plasticité cristalline. Trois exemples d’applications réalisées au
née et qui fait intervenir les interactions mutuelles entre les CEA sont présentés dans ce qui suit.
différents systèmes de glissement au moyen de variables
internes. Peuvent aussi intervenir la solution solide, la pré- Fatigue et fluage à haute température
sence de précipités, les vallées de Peierls…
des aciers martensitiques
• une loi d’écoulement reliant l’incrément de glissement visco- Les aciers martensitiques sont des candidats potentiels
plastique d’un système ṁs à la cission résolue s et à la cis- comme matériaux de structures des réacteurs de quatrième
sion critique {s ; génération. On observe expérimentalement que leur micro-
structure très fine, composée de sous-joints, disparaît progres-
sivement en cours de déformation, ce qui conduit à un adou-
cissement notable du matériau. Les sous-joints peuvent être

représentés comme des réseaux de familles de dislocations
[7] portant une désorientation angulaire s dont l’évolution
Le durcissement d’irradiation
dépend du glissement viscoplastique, par exemple dans le cas Le durcissement d’irradiation résulte de l’interaction des dislo-
d’un joint de flexion symétrique composé de dislocations coins cations avec les défauts produits par l’irradiation neutronique.
[2] : La dynamique moléculaire et la dynamique des dislocations
|s = — s|m|s|
–~ peuvent fournir des informations quantitatives sur les forces
d’obstacles α& de ces défauts qui ancrent les dislocations.
Connaissant la densité s& de ces défauts sur chaque sys-
L’équation d’évolution de la densité de dislocations est similaire tème de glissement, le durcissement d’irradiation Fs& peut
¹ s|m| s|, repré-
—

alors être introduit par une relation de Taylor du type :
x
à (1), avec ajout d’un terme supplémentaire, – –
Fs& = α& s&

sentant l’annihilation de dislocations mobiles avec les disloca-
tions des sous-joints. La diminution de la densité de disloca-
tions et la croissance de la taille de sous-grain, , prédites
permettent alors de décrire l’adoucissement observé. avec le module de cisaillement du matériau. À basse tempé-
Dislocations vis et coin sont distinguées du fait de leur ciné- rature dans le fer de structure cubique centrée, la mobilité des
tiques d’annihilation différentes. dislocations vis est fortement réduite par rapport à celle des
dislocations coins, et l’on observe alors que le glissement plas-
À haute température, les mécanismes de montée des disloca- tique dépend directement de la longueur des segments vis for-
tions sous l’effet de la diffusion de lacunes sont actifs et accé- més, laquelle dépend à son tour de la distance moyenne entre
lèrent les mécanismes d’annihilation entre dislocations. Ces les obstacles qui épinglent les dislocations. Le durcissement
mécanismes prenant en compte la contrainte et l’éventuelle d’irradiation peut alors être introduit, dans ce cas, en recalcu-
supersaturation en lacunes, peuvent également être intégrés lant la distance moyenne entre obstacles, réduite par l’ajout
dans une loi de viscoplasticité cristalline [1]. des défauts d’irradiation dans la microstructure.
Des comparaisons entre prévisions et données expérimen- La densité de défauts d’irradiation peut être amenée à diminuer
tales sont détaillées dans le chapitre intitulé « Déformation et sur les systèmes de glissement actifs par un mécanisme de
endommagement des aciers martensitiques revenus à haute nettoyage, lors du passage des dislocations. Ce mécanisme, à
température », infra, pp. 175-179. l’origine, notamment, de la formation de bandes claires dans
les matériaux fortement irradiés (alliages de zirconium ou aciers
inoxydables austénitiques) [fig. 47] peut se traduire par la rela-
s&|m| s| où le coefficient corres-
s = ––
tion suivante [6] : |&
pond à la distance de capture du défaut d’irradiation par la dis-
220
111
0,3 µm
0,25 µm
a b
Fig. 47. Bandes claires dans un acier inoxydable austénitique irradié aux neutrons puis déformé [8]. a) Cisaillement de bandes claires
d’un système primaire par celles d’un système secondaire (flèches rouges) et (b) déformation au voisinage d’un joint de grains (traits continus)
impacté par une bande claire (traits en pointillés).
72 La viscoplasticité cristalline
Bande
S4()
L
45°
t ∑
Grain
Joint de grains
Matrice Matrice
a b
Fig. 48. Introduction de bandes de localisation de la déformation dans un grain lui-même noyé dans une matrice polycristalline. L’épaisseur et la
longueur des bandes sont notées t et L (qui s’identifie à la taille de grain) (a) ; champ de contrainte intergranulaire, S4(), résultant de l’interaction
de cette bande avec un joint de grains, en fonction de la distance au point d’impact de la bande, (b).
location. Cette diminution du nombre de défauts sur les sys-

Localisation de la déformation
tèmes de glissement actifs conduit à un adoucissement qui
favorise la localisation des déformations observée expérimen- La localisation de la déformation sous forme de bandes à
talement à toutes les échelles, depuis le grain jusqu’à, éven- l’échelle du grain est un phénomène que l’on observe égale-
tuellement, l’éprouvette de traction macroscopique. ment en fatigue, ou même en chargement monotone, et qui
peut donc être étudié avec les mêmes types d’outils numé-
La localisation intense de la déformation dans des bandes riques. On observe, par ailleurs, au cours de sollicitations
conduit à la formation d’une microstructure composite à uniaxiales monotones ou cycliques, que la déformation plas-
l’échelle du grain même (fig. 48a). La forte localisation de la tique dans les matériaux de structure cubique à faces cen-
déformation plastique dans ces bandes induit d’importantes trées résulte essentiellement d’un glissement simple dans la
incompatibilités de déformation entre la bande et les grains plupart des grains pour lesquels le rapport entre facteur de
voisins (fig. 48b) qui produisent un écrouissage cinématique SCHMID secondaire et primaire est suffisamment faible (direc-
(back stress) significatif. Ce dernier peut être également pris tions de traction au centre du triangle standard, figure 49a,
en compte dans les lois cristallines par un terme d’écrouis- cristaux dits « bien orientés »). Ce type de comportement n’est
sage cinématique intra-granulaire (de type ARMSTRONG- pas correctement décrit par la plupart des lois de plasticité
FREDERICK) dans la loi d’écoulement [6] : cristalline fondées sur les écrouissages isotrope et cinéma-
tique. Le glissement simple est imposé dans les grains bien
|s– s| – {s
m|s = ⟨—————– ⟩ b&4 (s – s ),
orientés et deux lois de comportement distinctes sont appli-
h quées selon que le cristal est en glissement simple (bandes de
glissement) ou multiple (cellules, labyrinthes…), en accord
s| = m|s – M| s |m|s|, avec les données expérimentales (fig. 49b) [9]. Les paramètres
où h, 4, et M sont des paramètres dépendant du matériau.

identifiés sont une contrainte résolue critique commune à tous
les grains et un jeu de deux paramètres d’écrouissage ciné-
matique non-linéaire d’ARMSTRONG-FREDERICKS par famille
Afin de calculer précisément les champs de contrainte locaux,
d’orientations. Le comportement macroscopique des polycris-
on peut également introduire dans un calcul aux éléments finis
taux est alors convenablement prédit à partir de la seule
utilisant des lois de plasticité cristalline une ou plusieurs
connaissance des courbes expérimentales mesurées sur
bandes de localisation dont le comportement diffère de celui
monocristaux (voir le chapitre suivant).
du grain dans lequel elles sont apparues (fig. 47a).

Contrainte axiale (MPa)
b
200
—– = !&]
Domaine de
glissement multiple
p
<111>
150
Domaine de 100
glissement simple
50
<100>
<-149>
0
<001> <110> 1.00E-05 1.00E-04 1.00E-03 1.00E-02
Déformation axiale
a b
Fig. 49. a) Définition des zones de glissement simple ou multiple en fonction de la direction de chargement dans un grain sollicité en traction-
compression. Le critère est fondé sur la comparaison entre le rapport des contraintes résolues secondaire et primaire avec une valeur critique
identifiée grâce aux observations des microstructures de dislocations. b) Identification du comportement contrainte-déformation en fonction
du type de glissement, simple ou multiple (<-149>, direction pour laquelle on n’active qu’un seul système de glissement ; <100>, direction
pour laquelle on active plusieurs systèmes de glissement). Les simulations (lignes) sont comparées aux valeurs expérimentales de BUQUE
et BOCHWITZ (symboles).
Ludovic VINCENT, Maxime SAUZAY

et Fabien ONIMUS,
Références [9] A. STECKMEYER, M. SAUZAY, A. WEIDNER et E. HIECKMANN,

[1] D. CAILLARD et J. L. MARTIN, Thermally Activated Mechanisms in « Micromechanical modelling of the cyclic stress-strain behaviour of
Crystal Plasticity, Pergamon, (2003). nickel polycrystals », International Journal of Fatigue, 40 (2012),
pp. 154-167.
[2] M. F. GIORDANA, P.-F. GIROUX, I. ALVAREZ-ARMAS, M. SAUZAY et
A. ARMAS, « Micromechanical modelling of the cyclic softening of [10] C. TEODOSIU, J. L. RAPHANEL et L. TABOUROT, « Finite element simu-
EUROFER 97 steel », Mat. Sci. Eng. A, (2012), à paraître. lation of the large elastoplastic deformation of multicrystals », Mecamat
91, eds. Raphanel et Sidoroff, Balkema (1993), pp. 153-160.
[3] U. F. KOCKS et H. MECKING, « Physics and phenomenology of strain
hardening: the FCC case », Progress in Materials Science, 48 (2003),
pp. 171-273.
[4] M. LIBERT, C. REY, L. VINCENT et B. MARINI, « Temperature dependant

polycrystal model application to bainitic steel behavior under tri-axial
loading in the ductile-brittle transition », International Journal of Solids
and Structures, 48 (2011), pp. 2196-2208.
[5] J. MANDEL, « Plasticité classique et viscoplasticité », vol. 97 of CISM

Courses and lectures, Springer Verlag, (1971).
[6] F. ONIMUS et J. L. BECHADE, « A polycrystalline modeling of the

mechanical behavior of neutron irradiated zirconium alloys », Journal
of Nuclear Materials, 384 (2009), pp. 163-174.
[7] W. T. READ et W. SHOCKLEY, « Dislocation Models of Crystal Grain

Boundaries », Physical Review, 78 (1950), pp. 275-289.
[8] M. SAUZAY, K. BAVARD et W. KARLSEN, « TEM observations and finite

element modelling of chennel deformation in pre-irradiated austenitic
stainless steels – Interactions with free surfaces and grains bounda-
ries », Journal of Nuclear Materials, 406 (2010), pp. 152-165.
74 La viscoplasticité cristalline
Les méthodes d’homogénéisation

en mécanique des milieux continus
L es méthodes d’homogénéisation* en mécanique des géométrique particulière au sein du polycristal. De plus, la

milieux continus visent à évaluer le comportement macrosco- phase cristalline est considérée comme homogène, les hété-
pique (ou moyen) d’un matériau hétérogène en s’appuyant sur rogénéités de dimensions inférieures à celle-ci étant négli-
la connaissance de sa microstructure et sur le comportement gées. Les champs de contrainte et de déformation dans
mécanique de chacune des phases constituant le matériau chaque phase cristalline sont donc estimés en moyenne.
(fig. 50). Certaines approches, comme celle de BRENNER et al. [3], per-
mettent néanmoins de tenir compte de façon statistique de
l’hétérogénéité intra-phase. Afin de décrire les interactions
L’approche analytique
entre les grains du polycristal, une approche statistique de
de l’homogénéisation type « champ moyen » est adoptée. Celle-ci consiste à consi-
Les approches analytiques de l’homogénéisation permettent dérer que, compte tenu du grand nombre de grains dans le
de déduire le comportement homogénéisé du matériau à par- VER, les effets de voisinage sur un grain unique peuvent être
tir d’une description statistique de ce dernier [1,2]. Dans le cas vus en moyenne comme équivalents aux interactions méca-
d’un polycristal, cette approche statistique se conçoit naturel- niques que subirait la phase cristalline correspondante (tous
lement, compte tenu du très grand nombre de grains conte- les grains de même orientation) en inclusion dans un milieu
nus dans le Volume Élémentaire Représentatif* (VER). Les homogène dont le comportement est celui du Milieu
grains ne sont pas considérés individuellement, ceux-ci sont Homogène Équivalent (MHE). Pour les polycristaux consti-
remplacés par la notion de « phase cristalline » (désignée par tués de grains équiaxes distribués aléatoirement, l’inclusion
l’indice ) qui regroupe l’ensemble de tous les grains de même représentative est de forme sphérique. Cependant, afin de
orientation cristallographique, caractérisée par les trois angles déterminer la réponse de la phase cristalline en interaction
d’Euler (φ1, 8, φ2) et la fraction volumique de grains corres- avec le MHE, il est nécessaire de connaître le comportement
pondante (g) (fig. 50). Dans cette description, les grains de du MHE, or c’est justement lui qui est l’inconnue du problème.
même orientation sont indiscernables et n’ont pas de position La formalisation de ce schéma auto-cohérent aboutit à une
équation implicite dont la résolution
nécessite la mise en place d’une pro-
cédure spécifique et complexe, dans
les cas les plus généraux (fig. 50) [1].
Au CEA, plusieurs approches de ce

type ont été mises en œuvre dans le
cas d’aciers à 9 % de Cr, d’aciers ren-
forcés par dispersion d’oxyde ou bien
d’alliages de zirconium. Les lois cristal-
lines caractérisant le comportement
50 m
des grains du polycristal sont celles
Polycristal Milieu homogène équivalent
présentées dans le chapitre intitulé

« La viscoplasticité cristalline », supra,
pp. 71-74. Dans le cas des alliages de
MHE MHE MHE zirconium irradiés, le modèle polycris-
tallin développé a pu être validé à la
+ +…= fois, d’une part au niveau macrosco-
pique sur des essais spécifiques de la

littérature et, d’autre part, sur des
observations réalisées en microscopie
électronique en transmission [4]. Ce
Fig. 50. Illustration de la démarche d’homogénéisation, dans le cas d’un polycristal. modèle a notamment permis de prédire
Description du polycristal en phases cristallines et principe du schéma auto-cohérent. le comportement d’alliages de zirco-

nium irradiés dans des conditions d’essais non encore testées Différentes procédures ont été introduites dans le code de cal-
par les techniques expérimentales à disposition. Ces simula- cul par éléments finis CAST3M [10], afin de faciliter la mise en
tions mettant en œuvre une démarche hiérarchique ont éga- œuvre de cette démarche d’homogénéisation numérique,
lement permis de justifier le choix du formalisme des modèles notamment pour la définition des conditions aux limites et
empiriques opérationnels. l’évaluation du comportement moyen d’une cellule élémen-
taire. Nous pouvons observer sur la figure 51 [5,6] (voir détails
dans la légende), d’une part que le comportement moyen tend
L’approche numérique vers le comportement homogénéisé recherché, lorsque la
de l’homogénéisation taille des cellules augmente, et, d’autre part, que la meilleure
convergence est obtenue avec des conditions aux limites
Parallèlement aux approches analytiques, l’approche numé-
périodiques.
rique consiste à estimer le comportement homogénéisé du
matériau à partir d’un calcul mécanique réalisé sur une cel-
Parallèlement à la méthode des éléments-finis, des méthodes
lule élémentaire telle que celles présentées au chapitre inti-
alternatives fondées sur la Transformée de Fourier Rapide
tulé « La production de microstructures de matériaux pour la
(« méthodes FFT » en anglais), en plein essor dans la com-
simulation », supra, pp. 67-69. Cependant, excepté le cas des
munauté « mécanique des matériaux », sont en cours de
matériaux à microstructure périodiques, les matériaux présen-
développement au CEA. Cette méthode repose sur une réso-
tent, en général, une microstructure aléatoire (comme, par
lution itérative du problème élastique linéaire écrit sous la
exemple, pour les polycristaux ou bien pour les fibres de com-
forme de l’équation intégrale de LIPPMANN-SCHWINGER [11, 13] :
posite). Il est alors nécessaire de simuler non pas une cellule
+ 0 * ! – !0: =
élémentaire mais un grand nombre de cellules élémentaires,
afin d’en déduire le comportement moyen, qui est alors une
où désigne la déformation moyenne imposée, le champ

estimation du comportement homogénéisé [5, 9]. Par ailleurs,
de déformation inconnu, ! le champ de propriétés élastiques

il s’avère également indispensable de s’assurer que le résul-
hétérogènes au sein de la cellule, !0 les propriétés élastiques

tat obtenu n’est pas trop dépendant du choix des conditions
d’un matériau de référence homogène sur la cellule et 0 le

aux limites ou de la taille des cellules utilisées pour ces calculs.
La taille des cellules doit, en effet, être suffisamment grande
tenseur de Green relatif à ce milieu. À chaque étape de l’algo-
rithme itératif, le produit de convolution * est ici évalué simple-
par rapport à l’hétérogénéité considérée.
ment dans l’espace de Fourier, tirant ainsi parti des algo-

rithmes performants de FFT.
CLDH
h11 (GPa)
320
CLMI
CLP
310
300
290
280
270
260
250

0 5 10 15 20 25 30 35 40 45
Fig. 51. Homogénéisation de composites SiC/SiC à l’échelle du toron (voir figure 45 du chapitre intitulé « La production de microstructures de
matériaux pour la simulation », supra, pp. 65-67). Évaluation du coefficient K11 du tenseur d’élasticité moyen (1 = direction perpendiculaire aux
fibres) en fonction de la taille des cellules ( = dimension de la cellule / diamètre des fibres) et du choix des Conditions aux Limites Périodiques
(CLP), Mixtes (CLMI) ou Déformation Homogène (CLDH). Sur ce graphique, chaque point est le résultat de la moyenne d’un grand nombre
de calculs réalisés sur différentes cellules élémentaires. Un exemple de champ de contrainte associé à ces calculs est également proposé
(figure de droite).
76 Les méthodes d’homogénéisation en mécanique des milieux continus

16

180
160 14
140 12
120 10
100
8
80
Hill-Hutchinson 6
60
Kröner 4
40
Éléments finis 2
20
0 0
0 0,002 0,004 –0,002 1E17 0,002 0,004
a Déformation axiale b Déformation axiale
Fig. 52. Prédictions des courbes d’écrouissage cyclique obtenues par les modèles à champ moyen ou la méthode des éléments finis : a) nickel ;
b) aluminium. L’identification des paramètres de plasticité cristalline s’effectue à l’échelle des monocristaux (fig. 49 du chapitre « La viscoplasticité
cristalline », supra, p. 72) [8].
La confrontation des méthodes expérimentales obtenues sur monocristaux en chargement

monotone ou cyclique. Ainsi, une large gamme de modules
analytiques et numériques
élastiques, de contraintes critiques et de coefficients d’écrouis-
d’homogénéisation sage a été prise en compte. Quel que soit le matériau consi-
Récemment, des travaux visant à comparer les démarches déré, les courbes macroscopiques monotones ou cycliques
d’homogénéisation numérique (en champ complet) et analy- prédites par le modèle de HILL-HUTCHINSON diffèrent de moins
tique (en champ moyen) ont également été menés au CEA. de 5 % de celles prédites par la méthode des éléments finis
Le cadre élastoplastique a été choisi ici, ce qui a motivé l’uti- [8], et ce pour des déformations plastiques allant jusqu’à 10-2
lisation des modèles classiques : le modèle thermoélastique (fig. 52) [14].
de KRÖNER et le modèle élastoplastique incrémental de HILL-
HUTCHINSON. Afin de donner une généralité à ces résultats, dif- À l’échelle polycristalline, l’accord entre les prédictions et les
férents métaux cubiques à faces centrées ont été considérés : courbes expérimentales est en général satisfaisant (fig. 53) du
aluminium, cuivre, austénite, nickel… Les lois de plasticité cris- fait qu’aucun paramètre cristallin n’est identifié grâces à ces
tallines ont été identifiées uniquement à partir de données courbes. Ces constats sont valides pour l’acier austénitique
160
140
120
100
80
60
300
40
20
200 0
-0,0003 -0,0002 -0,0001 0 0,0001 0,0002 0,0003
-20
-40
100 -60
-80
-100
0
-120
-140
N = 20 000 [Morrison et al.]
-100 -160
Calcul éléments finis
1.00E-05 1.00E-04 1.00E-03 1.00E-02
a Déformation axiale b Déformation axiale
Fig. 53. Comparaison entre courbes d’écrouissages cycliques mesurée et prédite par la méthode des éléments finis (a), et boucles d’hystérésis
(dites stabilisées) sur du nickel à température ambiante [8].

inoxydables 316L(N) [14,15] ainsi que le nickel [8], le cuivre et [8] A. STECKMEYER, M. SAUZAY, A. WEIDNER et E. HIECKMANN,
l’aluminium [15]. À l’échelle, les contraintes internes sont « Micromechanical modeling of the cyclic stress-strain behavior of nic-
proches des valeurs mesurées par diffraction des rayons X [8]. kel polycrystals », Int. J. of Fatigue, 40 (2012), pp. 154-167.
[9] F. BARBE, L. DECKER, D. JEULIN et G. CAILLETAUD, Int. J. Plast., 17

Enfin, des calculs d’homogénéisation numérique ont été réa- (2001), p. 513.
lisés sur un grand nombre de polycristaux représentatifs d’un
[10] http://www-cast3.cea.fr.
acier de cuve, après avoir identifié une loi de comportement
cristalline adaptée sur la réponse expérimentale macrosco- [11] É. CASTELIER, T. HELFER et J. PACULL, Spécifications d’un code de
pique du matériau. Afin de limiter les temps de calcul dans thermomécanique avec résolution par transformées de Fourier
cette phase d’identification de paramètres, une méthode d’ho- rapides, Note Technique SESC/LSC 09-045- Indice 0- Décembre
2009.
mogénéisation analytique a été utilisée. Une fois les para-
mètres de la loi cristalline choisis, la série de calculs d’agré- [12] S. BRISARD et L. DORMIEUX, « Combining Galerkin approximation
gats a pu être lancée et a conduit à un bon accord entre la techniques with the principle of Hashin and Shtrikman to derive a new
réponse macroscopique moyennée sur l’ensemble de ces cal- FFT-based numerical method for the homogenization of composites »,
culs et celle prédite par le modèle d’homogénéisation. Des Comp. Meth. in Appl. Mech. and Eng. (2012), article accepté.
résultats expérimentaux de mesure par diffraction de neutrons [13] L. GÉLÉBART et C. LESTRINGANT, FFT-based homogenization of the
in situ ont ensuite permis de valider à l’échelle locale ce type mechanical behavior of SiC/SiC composites, Note Technique
d’approche numérique de l’homogénéisation en comparant DEN/DANS/DMN/SRMA/LC2M/NT/2011-3255/A.
les valeurs expérimentales et simulées de contrainte moyenne
[14] F. RACHDI et M. SAUZAY, Prediction of the cyclic behavior of FCC
dans certaines phases cristallines [7]. polycrystals using mean-field homogenization, 16th Conference on
Strength of Materials, Bangalore, 2012.
[15] M. SAUZAY, P. EVRARD, A. STECKMEYER et E. FERRIÉ, « Physically-

Lionel GÉLÉBART, Fabien ONIMUS, Maxime SAUZAY
based modeling of the cyclic macroscopic behavior of metals »,
et Ludovic VINCENT,
Procedia Engineering, 2 (2010), pp. 531-540.
Références
[1] M. BORNERT, T. BRETHEAU et P. GILORMINI, Homogénéisation en
mécanique des matériaux 2, Comportements non linéaires et pro-
blèmes ouverts, Hermès Science Publications, Paris (2001).
[2] U. F. KOCKS, C. N. TOMÉ et H.-R. WENK, Texture and Anisotropy, first

Ed., Preferred Orientations in Polycrystals and their Effect on Materials
Properties, Cambridge University (1998).
[3] R. BRENNER, R. MASSON, O. CASTELNAU et A. ZAOUI, « A “quasi-elas-

tic” affine formulation for the homogenised behaviour of nonlinear vis-
coelastic polycrystals and composites », Eur. J. Mech. A/Solids, 21
(2002), p. 943.
[4] F. ONIMUS et J. L. BECHADE, « A polycrystalline modeling of the

mechanical behavior of neutron irradiated zirconium alloys », Journal
of Nuclear Materials, 384 (2009), pp. 163-174.
[5] C. CHATEAU, Analyse expérimentale et modélisation microméca-

niques du comportement élastique et de l’endommagement de com-
posites SiC/SiC unidirectionnels, Thèse, École Polytechnique (2011).
[6] C. CHATEAU, L. GÉLÉBART, M. BORNERT, J. CRÉPIN et D. CALDEMAISON,

« Multiscale approach of mechanical behaviour of SiC/SiC compo-
sites : elastic behaviour at the scale of the tow », Technische Mechanik,
30 (2010), pp. 45-55.
[7] L. VINCENT, L. GÉLÉBART, R. DAKHLAOUI et B. MARINI, « Stress locali-

zation in BCC polycrystals and its implications on the probability of
brittle fracture », Mat. Sci. Eng. A, 528 (2011), pp. 5861-5870.
78 Les méthodes d’homogénéisation en mécanique des milieux continus

la rupture des matériaux
Le terme ultime d’un essai de traction ou de fluage est la rup- Il y a rupture fragile du matériau fissuré si Sy dépasse la
ture de l’éprouvette. Cette rupture peut être ductile* ou fra- contrainte théorique de clivage S , ce qui arrive si la
gile* c’est-à-dire avec ou sans déformation plastique. contrainte macroscopique S dépasse
1 m.. 1/2
Considérons ici le second cas (rupture fragile), où les deux sur-
S = — —-—,
2 .3
faces de la pièce cassée peuvent être remises en contact l’une
avec l’autre, jusqu’à l’échelle atomique.
Or, le rayon de courbure minimal du fond de fissure est de
Si l’on cherche à écarter deux plans cristallins d’un matériau en l’ordre de la distance interatomique , donc
m. 1/2
traction (rupture dite « en mode I »), il faut exercer une
contrainte de traction qui dépend du module de Young et de S ≈ —— (critère de Griffith*).
la variation de l’écartement entre les plans. Cette contrainte 3
croît d’abord linéairement avec , atteint une valeur seuil cor- L’expérience montre qu’effectivement le produit de la racine
respondant à la contrainte théorique de clivage, puis retombe carrée de la longueur de défauts (3) et la contrainte à la rup-
à zéro au-delà d’une distance alpha entre les plans, désormais ture (S) est à peu près constant.
clivés. La forme la plus simple décrivant cette relation
contrainte-déplacement est une fonction sinus : La quantité h ≡ (2..m.)1/2 est une grandeur intrinsèque au
. .
S = —–. b&4 ——,
matériau, appelée « facteur d’intensité de contrainte* », qui
. caractérise sa ténacité* (résistance à la rupture fragile). On la
où est l’espacement entre plans atomiques et la distance

mesure classiquement par des essais de ténacité sur des
éprouvettes entaillées (essai Charpy). Elle s’exprime en
au-delà de laquelle l’interaction entre plans cristallins s’annule. Pa.m1/2. La connaissance de h permet de savoir à partir de
L’énergie de surface du clivage est le travail nécessaire pour quel niveau de contrainte macroscopique une fissure de lon-
séparer les plans, soit : gueur donnée pourra se propager dans le matériau.
.2
m = —– .
2.
Les calculs simplistes, ci-dessus, n’ont guère qu’une valeur
pédagogique permettant d’introduire les notions essentielles
La densité d’énergie surfacique m des matériaux est de l’ordre et les principaux ordres de grandeur. En fait, les choses sont
de 1 J/m2. La contrainte théorique de clivage est la valeur maxi- plus compliquées dans la mesure où la rupture du matériau
male de la sinusoïde, soit : peut aussi être provoquée par des contraintes de cisaillement
. m. 1/2
S = —–. b&4 ——,
(rupture en mode II et III). Par ailleurs, même en rupture dite
.
« fragile », il y a souvent en fond de fissure une zone où le
matériau subit une déformation plastique dans laquelle se pro-
S prend, en général, des valeurs très élevées, de l’ordre de duit une relaxation des contraintes et une dissipation de l’éner-
10 GPa, beaucoup plus grandes que les contraintes de rup- gie de rupture. On entre alors dans le domaine plus complexe
ture observées sur matériaux réels, de l’ordre de 100 MPa. En de la rupture ductile.
fait, la rupture est due à l’existence de défauts microscopiques
sous forme de fissures dans le matériau, qui concentrent la
contrainte en fond de fissure. Si on considère une fissure ellip-
tique, de longueur 2l et de rayon de courbure en fond de fis-
sure, la théorie de l’élasticité dit que la contrainte Sy en fond
de fissure est 2.(3 /)1/2 fois plus forte que la contrainte
moyenne dans le matériau (fig. 54).
Sy
3
Fig. 54. Matériau sous contrainte de traction, avec fissure

elliptique.

Approche locale et approche probabiliste

de la rupture des matériaux
L’approche locale de la rupture l’endommagement. Les champs mécaniques n’étant pas

homogènes, à cause de la géométrie et du chargement du
Le principe de l’approche locale de la rupture est de décrire les volume simulé, différents tirages aléatoires successifs du para-
mécanismes d’endommagement* à l’échelle d’un volume mètre d’endommagement conduiront à différentes valeurs de
dont le comportement peut être supposé comme relevant de l’instant critique d’amorçage de la rupture, reproduisant ainsi
la mécanique des milieux continus et de la plasticité standard. la dispersion de la ductilité du matériau. De plus, la taille du
Cet élément de volume dit « représentatif » (VER*) constitue VER étant constante, toutes les structures homothétiques
la brique élémentaire du modèle dans laquelle l’évolution de auront des distributions différentes d’instants de rupture. En
l’endommagement est décrite. L’endommagement est géné- particulier, à niveau de chargement équivalent, plus le volume
ralement représenté par un scalaire dans une loi d’évolution sollicité sera grand, plus l’amorçage de la rupture sera facile.
fonction de certaines composantes des champs mécaniques
régnant dans l’élément de volume. Cette loi d’évolution obte- Nous voyons ici qu’une telle approche permet potentiellement
nue par modélisation du mécanisme doit faire l’objet d’une vali- une description probabiliste de la rupture ainsi que la prédic-
dation expérimentale. On distingue deux cas : celui où l’en- tion des effets d’échelle. Dans les cas de rupture brutale, la
dommagement ne modifie pas de façon significative le probabilité de rupture macroscopique peut être obtenue en
comportement mécanique du VER et celui où l’endommage- appliquant le principe du maillon faible qui considère que le
ment est couplé au comportement. matériau est rompu quand un des VER qui le constitue est
rompu. De façon complémentaire, et si les événements de
Dans les cas où l’endommagement n’est pas couplé au com- rupture des VER sont indépendants, la probabilité de survie
portement mécanique, l’intégration de la loi d’évolution au d’un volume de matériau / est le produit des probabilités de
cours du chargement et la comparaison avec la valeur critique survie de chacun des VER qui le constituent. Ainsi, la proba-
déterminée expérimentalement permet de définir l’instant de bilité de rupture de / soumis à un chargement sera don-
la rupture. Ce calcul est réalisé en post-traitement d’un calcul née par :
élastoplastique de la structure à analyser.
A (/, ) = 1– Π [1 – A (/0, S)] (1)
Pour décrire les cas où l’endommagement est couplé au com- /0∈ /
portement mécanique, la variable représentant l’endommage-

ment est souvent introduite dans le potentiel plastique servant Lorsque plusieurs mécanismes de rupture sont en compéti-
à l’évaluation des champs mécaniques dans le VER. Dans ce tion, cette expression peut être généralisée ; la probabilité de
cas, l’endommagement peut introduire un comportement survie du volume / étant le produit des probabilités de survie
adoucissant qui entre en compétition avec l’écrouissage du associées à chaque mécanisme [5].
matériau. L’évaluation de l’endommagement se fait donc
simultanément à celle des champs mécaniques et l’instant de L’exemple de modélisation de référence de l’approche locale
la rupture du VER peut être directement déterminé par la probabiliste de la rupture est le modèle de BEREMIN pour la
chute des contraintes dans le VER du fait du développement rupture par clivage transgranulaire des aciers faiblement alliés
de l’endommagement. L’intérêt de cette approche est de simu- [3]. Dans ce modèle, la rupture par clivage transgranulaire d’un
ler un endommagement progressif du matériau et la redistri- volume élémentaire représentatif est obtenue lorsque ce
bution des champs mécaniques associée, permettant ainsi de volume est en train de se déformer plastiquement et que la
décrire la propagation des fissures. plus grande contrainte principale atteint une contrainte critique.
La rupture de ce VER est probabilisée en considérant que les
particules à l’origine des microfissures responsables de la rup-
L’approche probabiliste ture ont une taille qui suit une distribution asymptotique de
de la rupture type loi puissance. Pour un volume constitué de plusieurs
VER, la probabilité de rupture par clivage s’écrit alors :
Un atout majeur de l’approche locale de la rupture est d’être
S¥ () ¦
utilisable pour reproduire les dispersions des propriétés de
A (/, ) = 1 – a p § – ——— ¨,
S
rupture des matériaux. Pour cela, un paramètre physique dis- (2)
tribué statistiquement est introduit dans la loi d’évolution de

où S¥ est appelée « contrainte de WEIBULL » et y et S sont La modélisation des réseaux
des caractéristiques du matériau, qui dépendent, en particu- de fissures
lier, de l’exposant de la loi puissance décrivant le comporte-
ment asymptotique de la taille des particules considérées plus La description probabiliste précédente considère, d’une part,
haut. Cela illustre la capacité de l’approche locale à prendre en une indépendance des évènements et, d’autre part, que la
compte certains aspects microstructuraux des matériaux, éta- rupture est obtenue dès l’apparition du premier VER rompu. Il
blissant ainsi un lien entre la mécanique et la métallurgie. Cette existe certaines situations, comme dans le cas d’un endom-
capacité constitue un avantage décisif de l’approche locale de magement de fatigue, pour lesquelles de nombreuses fissures
la rupture et représente l’assurance d’une bonne transférabi- peuvent apparaître sans pour autant entraîner la ruine de la
lité aux structures industrielles des propriétés des matériaux structure. L’amorçage de ces fissures et leur propagation
déterminées en laboratoire. dépendent alors de l’état d’endommagement existant du
matériau, une schématisation de ce phénomène étant propo-
sée sur la figure suivante [4]. Parmi tous les sites pouvant
Endommagement et rupture potentiellement conduire à l’amorçage d’une nouvelle fissure
dans le volume élémentaire (point rouges sur la figure 55), seuls ceux (points blancs)
situés en dehors des zones de décharge* (oranges) entou-
Nous avons considéré jusqu’ici le volume élémentaire repré- rant les fissures déjà existantes pourront effectivement amor-
sentatif comme un volume homogène en termes de compor- cer une fissure. De même, les plus petites fissures seront arrê-
tement. Le développement récent des techniques de modéli- tées lorsqu’elles tomberont dans la zone de décharge de
sation des microstructures et des comportements cristallins fissures plus grandes.
permet de décrire le VER comme un ensemble de domaines
ayant des propriétés mécaniques distinctes. Pour modéliser ce type d’endommagement progressif, on
considère le voisinage entourant un défaut D susceptible
Pour les matériaux métalliques, il est possible de considérer le d’amorcer une fissure à + cycles et dans lequel il ne faut pas
VER comme un « agrégat de grains » caractérisés par leur rencontrer de fissures. La zone à considérer dépend, en fait,
géométrie et leurs orientations cristallines, dans lequel les du nombre de cycles +© pendant lesquels une fissure peut
champs mécaniques sont hétérogènes. Ces champs peuvent amorcer, ou plus exactement du nombre de cycles pendant
être soit utilisés directement, soit être représentés par des dis- lequel cette fissure va se propager, à savoir (+-+©) cycles. En
tributions statistiques pour décrire l’endommagement et la rup- effet, plus une fissure apparaîtra tôt, plus sa taille aura aug-
ture du VER plus finement que précédemment. Un exemple menté par propagation et, par conséquent, plus il faut que
d’application de cette possibilité est présenté dans le chapitre cette fissure ait amorcé loin de D pour ne pas l’écranter
intitulé « La modélisation physique du comportement méca- (fig. 55a). Partant d’une distribution de Poisson, la probabilité
nique des alliages de zirconium, avant, après irradiation et en élémentaire pª« (+©) de ne trouver aucune fissure nouvelle-
conditions accidentelles », infra, pp. 193-196. ment créée à +© cycles dans une zone ¬ª« (+-+©) s’écrit :
® (+©)
pª« (+©) = a p § – ——–— ¬ª« (+-+©) +©¨ ,
+
——
N=0 N = Ni1 N = Ni2 N = Ni3 N = Ni4
Site obscurci Fissure active entourée de sa zone d’obscurcissement Fissure stoppée entourée de sa zone d’obscurcissement
Fig. 55. Évolution de l’endommagement d’un matériau par fissuration multiple. Les grandes fissures empêchent la germination et la propagation
des fissures situées dans leur zone d’obscurcissement, zone dans laquelle le champ de contraintes est localement relaxé.
82 Approche locale et approche probabiliste de la rupture des matériaux

20 ®± amorcées
®± actives
Densité de fissures (10-3/mm2)

Nombre de cycles
Probabilité d’obscurcissement
¬ª« (0)
1
+ 0,8 15
D
0,6
¬ª« (+-+©)
10
+© 0,4
5
¬ª« (+-+±²³)
0,2
+±²³
0 0
Espace (n)
0 106 107 104 105 106 107
a Zone d’étude b Nombre de cycles c Nombre de cycles
Fig. 56. a) Définition de la zone espace-Nombre de cycles dans laquelle il faut ne trouver aucun défaut pour que le défaut D puisse amorcer une
fissure après N cycles. b) Évolution de la probabilité d’obscurcissement Aª«. c) Évolution des densités de fissures amorcées ®± et actives ®±.
Application au cas d’un chargement de chocs thermiques cycliques produisant une décroissance rapide du chargement et un arrêt des fissures
dans l’épaisseur du matériau (conditions favorables pour le développement d’un réseau important de fissures en surface).
avec ® (+) la densité de fissures amorçables indépendam- Une application de cette modélisation est proposée sur la
ment de tout effet de voisinage. On réutilise alors la notion de figure 56b et c, dans le cas de chocs thermiques cycliques
maillon faible [eq. (1)], en considérant qu’il suffit de trouver au induisant une décroissance rapide du chargement et un arrêt
moins un défaut dans une des zones ¬ª« (+-+©) l’entourant des fissures dans l’épaisseur du matériau. Ce type de char-
pour empêcher un défaut D d’amorcer une fissure à + cycles, gement favorise alors l’amorçage et la propagation de mul-
ce qui s’écrit : tiples fissures en surface d’échantillons, là où le chargement
Aª« (+)=1– ∏ pª« (+©)=1– a p §– ° ——–— ¬ª« (+-+©) +©¨

® (+©)
+ + est encore élevé, d’où la présence de réseaux denses de fis-
+
—— sures macroscopiques.
+© =+y&4
+ y&4
(3) Toutefois, ce type de modèle peut également s’appliquer au

Cette probabilité d’obscurcissement est en fait une variable cas de sollicitations de fatigue uniaxiale isotherme classique
d’endommagement du matériau qui permet également de cal- (type traction/compression), et l’on constate alors que l’on par-
culer la densité de fissures amorcées ou actives (qui se provient à prédire l’évolution des tout premiers stades de dévelop-
pagent encore). pement de réseaux de fissures, avant qu’une fissure macro-
180 MPa Type I Type I

4 190 MPa 220 MPa 80 Type II 80 Type II
Densité de fissures (mm-2)
210 MPa Modèle Type III Type III
3 60 60
2 40 40
1 20 20
0 0 0
10 000 100 000 1 000 000 0 0,2 0,4 0,6 0,8 0 0,2 0,4 0,6 0,8 1
a Nombre de cycles b N/NR c N/NR
Fig. 57. a) Évolution de la densité de fissures amorçables à l’échelle du grain sous chargement de fatigue uniaxial. Comparaison
modèle/expérience sur un acier inoxydable austénitique de type 304L à température ambiante. b) Mesures expérimentales sur un autre type
d’acier inoxydable (316L) de l’évolution de la densité de fissures en fonction du nombre de cycles N rapporté au nombre de cycles à rupture NR
et de la taille des fissures (type I : une taille de grain, type II : deux tailles de grains. Type III : plus de trois tailles de grains), d’après [2].
c) Prévision par le modèle probabiliste de l’évolution des mêmes évolutions de tailles de fissures. D’après [6].

scopique ne l’emporte sur toutes les microfissures voisines
(fig. 57). Ce résultat est en partie lié à la loi probabiliste rete-
nue pour l’amorçage des fissures à l’échelle microscopique
[6].
Enfin, les prévisions de ce type de modèle probabiliste ont été

comparées avec succès à des simulations stochastiques 3D
aux éléments finis faisant intervenir des tirages aléatoires pour
les nombres de cycles à l’amorçage des défauts [1].
Bernard MARINI et Ludovic VINCENT,

Références
[1] J. BARES, L. GÉLÉBART, J. RUPIL et L. VINCENT, International Journal
of Fatigue, à paraître.
[2] A. BATAILLE et T. MAGNIN, Acta Metallurgica et Materialia, 42 (1994),

pp. 3817-3825.
[3] F.M. BEREMIN, Metallurgical Transactions A, 14 (1983), pp. 2277-

2287.
[4] N. MALÉSYS, L. VINCENT et F. HILD, International Journal of Fatigue,

31 (2009), pp. 565-574.
[5] B. MARINI, Journal de Physique IV, 8 (1998), pp. 95-104.
[6] J. RUPIL, L. VINCENT, F. HILD et S. ROUX, International Journal for

Multiscale Computational Engineering, 9 (2011), pp. 445-458.
[7] W. WEIBULL, Ingeniörs Vetenskaps Akademiens, Hanlingar, 151

(1939).
84 Approche locale et approche probabiliste de la rupture des matériaux

Monographie Materiaux Du Nucleaire Outils Modelisation

Transféré par

Droits d'auteur :

Formats disponibles

Monographie Materiaux Du Nucleaire Outils Modelisation

Transféré par

Informations du document

Copyright

Formats disponibles

Partager ce document

Partager ou intégrer le document

Options de partage

Avez-vous trouvé ce document utile ?

Ce contenu est-il inapproprié ?

Droits d'auteur :

Formats disponibles

Monographie Materiaux Du Nucleaire Outils Modelisation

Transféré par

Droits d'auteur :

Formats disponibles

Les outils de la modélisation

Les matériaux du nucléaire 25

10 ps ; 10-100 nm3 103 s ; 10-100 µm3

Simulations Modèles cinétiques Propriétés mécaniques :

Propriétés des dislocations

26 L’approche multi-échelles de la modélisation des matériaux

Les matériaux du nucléaire 27

Mécanismes Comportement du monocristal

François WILLAIME et Bernard MARINI,

28 L’approche multi-échelles de la modélisation des matériaux

Les calculs électroniques ab initio

La théorie de la fonctionnelle ε ! désigne la densité d’énergie d’un système de référence

dient $(). Cette approximation, appelée « approximation du

Les matériaux du nucléaire 29

grandes lignes de la mise en œuvre pratique de la DFT. La

Les calculs numériques se séparent en deux grandes

30 Les calculs électroniques ab initio pour la matière condensée

L’autre grande famille de fonctions de base est constituée par

Les matériaux du nucléaire 31

[2] W. KOHN et L.J. SHAM, « Self-Consistent Equations Including

[3] W. KOHN, « Nobel Lecture: Electronic structure of matter – wave

[4] L.H. THOMAS, « The calculation of atomic fields », Proc. Cambridge

[5] E. FERMI, « Un Metodo Statistico per la Determinazione di alcune

[6] D. M. CEPERLEY et B. J. ALDER, « Ground State of the Electron Gas

[7] X. GONZE, B. AMADON, P.-M. ANGLADE, J.-M. BEUKEN, F. BOTTIN,

[8] M. TORRENT, F. JOLLET, F. BOTTIN, G. ZERAH and X. GONZE,

32 Les calculs électroniques ab initio pour la matière condensée

Les potentiels interatomiques

Étendre le domaine spatio-temporel des calculs ab initio

Potentiel répulsif requiert une simplification du modèle physique. Une approche

Les matériaux du nucléaire 33

34 Les potentiels interatomiques

Les matériaux du nucléaire 35

Configurations ab initio Ajustement

36 Les potentiels interatomiques

[2] F. ERCOLESSI et J. B. ADAMS, Eur. Lett., 26 (1994), p. 583.

[3] M.-C. MARINICA et F. WILLAIME, Solid State Phenomena, 129 (2007),

[4] M. I. MENDELEV et al., Phil. Mag., 83 (2003), p. 3977.

[5] M.-C. MARINICA, F. WILLAIME et J.-P. CROCOMBETTE, Phys. Rev. Lett.,

Les matériaux du nucléaire 37

Le paysage énergétique à l’échelle atomique

Les matériaux du nucléaire 39

40 Le paysage énergétique à l’échelle atomique

Barrière énergétique (eV)

Les matériaux du nucléaire 41

[2] D. RODNEY et L. PROVILLE, Phys. Rev. B, 79 (2009), p. 094108.

42 Le paysage énergétique à l’échelle atomique

Les simulations de dynamique moléculaire

• D’une part, l’étude des systèmes atomiques proches de

Les matériaux du nucléaire 43

physique que l’on souhaite simuler. Avec une description des

[11] N. MOUSSEAU et al., Journal of Atomic, Molecular and Optical

44 Les simulations de dynamique moléculaire

Les matériaux du nucléaire 45

thermodynamique s’obtient à partir des travaux thermodyna- F(Q6, T)

miques récupérés lorsque l’on transforme le système. La 3

méthode d’intégration thermodynamique apparaît comme le 2 0,16 0

cas particulier d’une dynamique moléculaire hors-équilibre 1 0,18

46 Potentiels thermodynamiques et potentiels de force moyenne

incorrectement l’état de transition.

dient $(). Cette approximation, appelée « approximation du