Développement d’une méthode Monte Carlo cinétiques à l’échelle cellulaire

Télécharger le fichier pdf d’un mémoire de fin d’études

Simulation Monte Carlo atomique cinétique

La méthode de Monte Carlo fut développée à l’origine pour simuler la diﬀu-sion des neutrons par von Neumann, Ulam et Metropolis [114]. Elle fut ensuite rapidement appliquée à l’étude des dégâts d’irradiation dans les métaux [26, 53].

Une de ses variantes, la méthode de Monte Carlo cinétique à l’échelle atomique (Atomic Kinetic Monte Carlo simulations ou AKMC), est un outil de choix pour étudier les cinétiques de transformation de phase et permet de simuler avec un grand succès la précipitation dans des alliages réels [48, 107, 49]. Son utilisation fut prolifique et le présent paragraphe n’entend pas être une description exhaustive de ces succès. Lorsque les précipités restent cohérents avec la matrice et que les eﬀets élastiques demeurent modérés, une modélisation du cristal par un réseau rigide dont les noeuds sont occupés par les atomes permet une description réaliste de la transformation de phase. Dans un tel cas, les simulations AKMC permettent de simuler la précipitation dans un cristal de quelques millions d’atomes.

Ces simulations permettent de prendre explicitement en compte le mécanisme de diﬀusion, qu’il soit lacunaire [157, 10, 139, 127, 125, 12, 138, 141], interstitiel [60, 73, 72] ou de type balistique comme sous irradiation [139, 55, 154], pour re-produire le chemin de précipitation d’un alliage. Cette considération des spécificités du mécanisme de diﬀusion lui procure un pouvoir prédictif concernant la cinétique de transformation de phase, et conduit à la description d’une large variété de com-portements cinétiques. Notamment, les simulations AKMC permettent de prédire la migration des atomes isolés et de petits amas d’atomes de soluté qui influe sur le taux de germination et le temps d’incubation [142, 98]. Elles permettent aussi une description des interfaces et des eﬀets cinétiques sur ces dernières [133, 134, 107]. Dans les ternaires, des composés transitoires se forment dont la composition ne peut généralement pas être supposée a priori. Ils peuvent par contre être prédit par simulation Monte Carlo [60, 107]. Enfin de tels simulations paramétrées sur des calculs DFT se montrent en très bon accord avec les résultats expérimentaux tels qu’ils sont obtenus par sonde atomique tomographique [48, 107, 124], permettant la construction bottom-up d’une modélisation multiéchelle des matériaux.

Principe

Les méthodes Monte Carlo consistent à utiliser une méthode probabiliste pour calculer une intégrale. En physique statistique et en science des matériaux, elles sont souvent employées pour calculer des valeurs moyennes. Les valeurs moyennes d’une observable thermodynamique peuvent être calculées en eﬀectuant une telle moyenne d’ensemble sur l’espace des phases, tandis qu’un calcul sur l’espace des trajectoires dans l’espace des phases est nécessaire pour étudier les propriétés cinétiques de tels systèmes.

Calculs thermodynamique par la méthode Monte Carlo Les méthodes de Monte Carlo sont extrêmement eﬃcaces pour décrire les propriétés d’équilibre d’un système. Elles permettent ainsi par exemple de tracer le diagramme de phase, d’ex-plorer les profils d’interface (voir [131, 52] par exemple). Ces propriétés d’équilibre sont dictées par la fonction de partition. Or, le calcul exact de cette dernière est hors de portée (A l’exception notable des modèles d’Ising 1D et 2D. Le lecteur intéressé pourra se référer à l’ouvrage [128] pour plus de détails). Cependant, lors de l’éva-luation d’une observable, seul un nombre limité de configurations produisent une contribution notable à la valeur moyenne de cette observable. Les algorithmes de si-mulation Monte Carlo thermodynamique s’appliquent donc à générer une chaîne de configurations explorant plus précisément les configurations de contribution notable, aux dépens de celles de contribution négligeable. Cette exploration sélective se fait en sélectionnant à chaque étape la configuration suivante de la chaîne en se basant sur la distribution de probabilité des états. Ainsi les états les plus probables (i.e. de plus grand poids dans la moyenne) sont explorés préférentiellement, permettant une convergence rapide de l’estimateur de la moyenne. Pour un système thermique, cette distribution de probabilité est typiquement la distribution de Boltzmann.

D’une étape de la chaîne de Monte Carlo à l’autre, des états successifs sont explorés sans s’inquiéter du lien les unissant. Citons l’exemple canonique d’une si-mulation Monte Carlo thermodynamique du moment magnétique du système d’Ising par l’algorithme de Metropolis [114]. A chaque étape,

– un site du réseau est sélectionné aléatoirement,

– la diﬀérence E entre l’énergie du systeme tel qu’il est et après retournement du spin du site sélectionné est calculé,

– un nombre aléatoire R est tiré sur l’intervale [0; 1[,

– si R < e E, le spin est retourné,

– la somme des moments magnetiques est incrémentée du moment du système. La somme est finalement divisée par le nombre d’étapes.

Ainsi, les états les plus probables sont explorés préférentiellement. De plus, un calcul de l’énergie de l’ensemble du système à chaque étape serait extrêmement couteux : le temps CPU d’un tel calcul étant en général proportionnel à la taille du système. La restriction à une diﬀérence d’énergie permet de considérer chaque transition comme un phénomène local dont le temps de calcul est indépendant de la taille de système. D’autre part, on peut noter qu’il importe peu de savoir la raison qui pourrait conduire au retournement du spin selectionné, seule la probabilité d’équilibre de l’état considéré est pertinente pour ce calcul.

Monte Carlo cinétique L’objectif de l’algorithme Monte Carlo cinétique est de décrire l’évolution temporelle du système, et non plus seulement ses propriétés d’équilibre. Il va donc s’intéresser précisément à la manière de passer d’un état à l’autre.

L’utilisation d’un algorithme Monte Carlo est possible à condition que l’évolution puisse être considérée comme un processus de Markov, ou chaque transition est décrite comme aléatoire et indépendante de l’histoire du système, mais pouvant dépendre de l’état présent.

L’évolution de la probabilité P (n; t) de trouver le système dans une configuration donnée n à un temps t est alors donnée par l’équation maîtresse :

dP (n; t)XX = W (n ! n)P (n; t)W (n ! n)P (n; t);(1.3) dt nn où W (n ! n) est la probabilité par unité de temps de passer de la configuration n à la configuration n.

Le calcul des probabilités P (n; t), i.e. la résolution de l’équation maîtresse peut être eﬀectuée formellement [70, 156, 46]. Pour cela l’arbre directionnel décrivant toutes les chaines de configuration k possibles menant à cette configuration n à l’instant t doit être calculé, et chaque chaîne de cet arbre doit être multipliée par le produit de toutes les probabilités de transitions entre les configurations k qui la forment. Cette probabilité doit ensuite être normalisée, ce qui nécessite le calcul d’une fonction de partition étendue dénombrant toutes les probabilités de transition de tous les états. Le calcul exact de la fonction de partition du système d’Ising 3D étant déjà hors de portée, le calcul d’une telle fonction de partition étendue est a fortiori hors d’atteinte.

En conséquence, l’objectif des simulations AKMC n’est pas de calculer cette fonction de partition étendue. Plus modestement, il est de calculer une seule tra-jectoire du système, trajectoire que l’on souhaite la plus représentative possible de la trajectoire “moyenne”.

Ainsi à chaque étape la transition W (n ! n) est sélectionnée de manière probabilistique. Une telle opération est souvent eﬀectuée en utilisant un algorithme à temps de résidence [155, 32].

A chaque étape d’un tel algorithme,

– la liste fng des configurations atteignables avec une seule transition (évène-ment) à partir de l’état présent est établie,

– les taux de probabilité de transition W (n ! n) de chacune de ces transitions sont calculés,

– un nombre aléatoire R de probabilité uniforme sur l’intervale [0; 1[ est tiré

–l’évènement i vérifianti1W (nn)R <iW (nn)

Pnj2fng!Pnj2fng! est choisi,jj

–le temps physique de simulation est incrémenté par l’inverse de la somme des taux de transition t =1,

–n W (n!n)

la transition ainsi choisie est eﬀectuée.

Par exemple, l’ensemble des transitions possibles pour une diﬀusion par mécanisme lacunaire est typiquement l’ensemble des échanges entre chacune des lacunes et chacuns des atomes situés sur un site premier voisin.

L’algorithme de Metropolis peut aussi être utilisé pour de tels calculs, moyennant une définition a priori du temps à chaque étape de calcul. Ces deux algorithmes, temps de résidence et Metropolis, sont en réalité représentatifs de 2 classes d’al-gorithmes, l’un appelé Rejection free pour ceux comme le temps de résidence qui procèdent toujours à une évolution du système, et Null event pour ceux du type de Métropolis. Si leurs résultats sont identiques, leurs performances en temps de calcul varient selon les cas. Pour un système présentant un fort écart entre les probabilités des diﬀérents types d’évènement, un algorithme de type Rejection Free est préfé-rable [43], mais il perd en eﬃcacité lorsque le nombre de transitions augmente [47]. Il est cependant possible de construire un algorithme adaptatif exploitant les forces des deux classes d’algorithme [1, 11].

Description du système à l’échelle atomique

Modèle thermodynamique Le calcul de l’incrément d’énergie joue un rôle cru-cial dans les simulations Monte Carlo thermodynamiques comme il a été rappelé au paragraphe précédent. Dans l’approximation du réseau rigide, cet incrément d’éner-gie doit être calculé à partir des interactions atomiques par rapport aux positions d’équilibre des atomes, correspondant aux sites du réseau cristallin dans l’approxi-mation du réseau rigide. Si des atomes en positions interstitielles sont considérés, ces positions doivent bien entendu être elles aussi prises en compte.

Une configuration est en conséquence définie par le vecteur d’occupation n, dont chacune des composantes représente l’occupation d’un site : n = fnA1; nB1; : : : ; nv1; nA2; nB2; : : : ; nv2; : : :g, où v désigne une lacune, A,B, etc. les espèces chimiques en présence, et où nAi = 1 lorsque le site i est occupé par un atome de type A, et demeure nul sinon.

L’énergie d’une configuration est obtenue grâce à l’Hamiltonien du système E = H(n). Toute fonction des variables d’occupation peut être écrite sous la forme d’un développement en amas. Sans perte aucune de généralité, il est donc possible d’écrire l’Hamiltonien sous la forme d’un développement en amas :

Le modèle d’Ising correspond ainsi à une troncature de l’Hamiltonien à des inter-actions de paire de portée finie. Afin de reproduire les propriétés des matériaux des Hamiltonien plus complexes sont parfois choisis, où des interactions entre amas de trois atomes (ou plus) interviennent [47]. Une telle extension n’entraîne pas de diﬃ-cultés formelles pour les simulations Monte Carlo, mais engendre un ralentissement très important des simulations, en augmentant le temps de calcul de la variation d’énergie.

Modèle cinétique Pour une simulation AKMC dans l’approximation du réseau rigide, chaque transition correspond au déplacement d’un défaut ponctuel d’un site à l’autre. L’hypothèse d’indépendance permettant la résolution de l’équation maîtresse par la méthode de Monte Carlo consiste d’un point de vue physique à découpler les vibrations atomiques, très rapide, et des transitions entre sites. Ces transitions sont liées à un mécanisme, qui peut être lacunaire ou basé sur les atomes en position interstitielle.

Selon la théorie standard de l’état de transition [56, 151], la probabilité de transition prend la forme :

W (n ! n) = n0!nexp( Enact!n); (1.5) où n0!n est la fréquence d’essai (dites aussi fréquences d’attaque) et Enact!n est l’énergie d’activation de la transition. Cette dernière correspond à la diﬀérence d’énergie entre l’état de plus forte énergie du chemin réactionnel, aussi appelé point col, et l’état initial. La fréquence d’essai correspond quant à elle dans l’approximation quasi harmonique au rapport entre la somme des fréquences propres de vibration des atomes dans les configurations initiale et celle en position de col. Pour un système à N atomes, on dénombre alors 3N 3 fréquences propres en position d’équilibre, les fréquences associées au mouvement du centre de masses n’ayant pas d’influence, et 3N 4 en position de col, le mouvement selon l’axe de la transition devant aussi être omis. Cette description a pour conséquence directe de garantir le respect de l’équation de Bilan Détaillé : W (n ! n)P 0(n) = W (n ! n)P 0(n); (1.6) où P 0(n) désigne une probabilité d’équilibre de trouver le système dans la configu-ration n.

L’énergie en position initiale correspond à l’énergie en position d’équilibre telle qu’elle a été décrite au paragraphe précédent. L’énergie en position de col peut être définie de manière similaire à l’énergie en position d’équilibre, en utilisant un développement en amas de l’Hamiltonien qui prenne en compte les positions de col. L’utilisation d’un tel développement en amas restreint aux interactions de paire, tant pour décrire les positions de col que les positions d’équilibre, est souvent résumé par l’expression “modèle de liaisons coupées”. La connaissance de ces interactions ainsi que des fréquences d’essai suﬃt alors à paramétrer un tel modèle.

Le rapport de stage ou le pfe est un document d’analyse, de synthèse et d’évaluation de votre apprentissage, c’est pour cela rapport-gratuit.com propose le téléchargement des modèles complet de projet de fin d’étude, rapport de stage, mémoire, pfe, thèse, pour connaître la méthodologie à avoir et savoir comment construire les parties d’un projet de fin d’étude.

Table des matières

Introduction
1 Simulation de transformations de phase diffusives par la méthode de Monte Carlo cellulaire cinétique
1.1 Méthodes de simulation numérique de cinétiques de transformation de phase diffusives
1.1.1 Simulation Monte Carlo atomique cinétique
1.1.1.1 Principe
1.1.1.2 Description du système à l’échelle atomique
1.1.2 Modèles continus
1.1.2.1 Simulations en champ moyen
1.1.2.2 La Thermodynamique des Processus Irréversibles
1.1.2.3 DICTRA
1.1.2.4 Simulations en Champ de Phase
1.1.3 Simulation Monte Carlo Multi-échelles
1.1.3.1 Simulations Monte Carlo avec transformée en ondelettes de l’Hamiltonien
1.1.3.2 Méthode Coarse Grained Monte Carlo
1.1.3.3 Méthode 􀀀 leap
1.2 Développement d’une méthode Monte Carlo cinétiques à l’échelle cellulaire
1.2.1 Description du système
1.2.2 Echanges entre cellules quasi infinies
1.2.3 Traitement phénoménologique de la discrétisation des lacunes
1.2.3.1 Corrélations cinétiques dans le mécanisme lacunaire
1.2.3.2 Modèles de transition CKMC
1.3 Simulations CKMC
1.3.1 Système de référence
1.3.1.1 Description à l’échelle atomique
1.3.1.2 Description à l’échelle mésoscopique
1.3.2 Reproduction des propriétés cinétiques macroscopiques par simulation CKMC
1.3.2.1 Mesure de la matrice d’Onsager par simulation CKMC
1.3.2.2 Effets cinétiques de taille finie
1.3.3 Naissance et mort d’un précipité
1.3.3.1 Dissolution d’un précipité
1.3.4 Précipitation
1.3.4.1 Précipitation dans un alliage faiblement sursaturé
1.3.4.2 Régime de décomposition spinodale
1.3.4.3 Précipité dans une solution sursaturée
1.4 Conclusion
2 Construction des grandeurs thermodynamiques à l’échelle mésoscopique à partir de l’échelle atomique
2.1 Approche analytique de la thermodynamique : la méthode du champ moyen
2.1.1 Principe du champ moyen et approximation de Bragg Williams
2.1.1.1 Fonction de partition du système d’Ising dans l’approximation de Bragg Williams
2.1.1.2 Approche variationnelle
2.1.2 La méthode de variation d’amas (CVM)
2.1.3 Avantages et limitations
2.2 Paramétrisation numérique des propriétés thermodynamiques des simulations CKMC
2.2.1 Description de l’énergie par une fonctionnelle mésoscopique
2.2.2 Paramétrisation de la fonctionnelle d’énergie et simulation mésoscopique thermodynamique
2.2.3 Article
2.3 Commentaires
3 Construction des grandeurs cinétiques à l’échelle macroscopique à partir de l’échelle atomique
3.1 La matrice d’Onsager
3.1.1 Description macroscopique de la cinétique
3.1.1.1 Thermodynamique des processus irréversibles
3.1.1.2 Diffusion lacunaire dans un alliage sur réseau
3.1.2 Description atomique de la cinétique des alliages
3.1.2.1 Equation Maîtresse
3.1.2.2 Description des alliages : cas limites usuels et description complète
3.1.3 Méthodes de calcul de la matrice d’Onsager
3.1.3.1 Théorie de la réponse linéaire
3.1.3.2 Path Probability Method
3.1.3.3 La méthode du Champ Moyen Auto-Cohérent (SCMF)
3.1.3.4 Simulation Monte Carlo Atomique
3.1.3.5 Commentaires
3.1.4 Conclusion
3.2 Calcul systématique par la méthode SCMF
3.2.1 Généralisation de la formulation de la méthode SCMF
3.2.1.1 Formalisme hors équilibre
3.2.1.2 Equation des moments
3.2.1.3 Résolution formelle pour un système anisotrope
3.2.1.4 Commentaires
3.2.2 Description des fréquences de saut
3.2.2.1 Nomenclature des fréquences
3.2.2.2 Approximation statistique
3.2.3 Performances numériques
3.3 Effet des approximations de la méthode SCMF sur des modèles cinétiques simples.
3.3.1 Approximation statistique
3.3.1.1 Effet sur un alliage dilué
3.3.1.2 Effet sur un alliage concentré
3.3.2 Effet de la portée des interactions effectives
3.3.2.1 Calculs analytiques explicites
3.3.2.2 Applications à des calculs DFT pour des alliages réels.
3.3.2.3 Calcul implicite complet des interactions effectives de paire
3.3.3 Conclusion
3.4 Application de modèles cinétiques avancés à des alliages rééls
3.4.1 Au-delà des interactions entre sites premiers voisins
3.4.1.1 Fréquences mises en jeu
3.4.1.2 Application aux alliages Ni(Cr) et Ni(Fe)
3.4.1.3 Application à l’alliage Fe(Cu)
3.4.1.4 Discussion
3.4.2 Alliage Ni(Si) sous contrainte
3.4.2.1 Modèle cinétique d’une structure FCC soumise à une contrainte biaxiale
3.4.2.2 Effet d’une contrainte biaxiale sur la matrice d’Onsager de l’alliage Ni(Si)
3.4.2.3 Conclusion
Conclusion
Annexes
A Réduction à gros grains de l’équation maîtresse
B Correction des effets cinétiques de taille finie
C Eléments de calcul de la matrice d’Onsager pour un alliage FCC
C.1 Alliage BCC avec interactions entre sites premiers, second et troisièmes voisins
C.2 Alliage FCC avec interactions entre sites premiers voisins