r5.03.38 Comportement mécanique isotrope élasto-viscoplastique multi-échelles sous irradiation: applications aux aciers de cuve irradiés#
Résumé :
Ce document présente l’écriture de la loi de comportement VISC_ISOT_PLAS. Il s’agit d’un modèle d’écoulement isotrope élasto-viscoplastique modélisant les effets multi-échelles de la plasticité ainsi que les effets d’irradiation. Une modélisation des grandes déformations est accessible via la cinématique GDEF_LOG.
Les paramètres du modèle comprennent les caractéristiques mécaniques élastiques du matériau (le module de Young et le coefficient de Poisson) ainsi que les caractéristiques mécaniques de la microstructure du matériau considéré: taille de grains et densité de dislocations pour l’état non-irradié, ainsi que concentration et taille des amas de solutés pour l’état irradié. Ces derniers peuvent être directement obtenus expérimentalement à partir d’essais de traction standards.
La loi fournit la contrainte de von Mises en fonction de la déformation (visco)plastique cumulée. L’écrouissage est entièrement isotrope, ce qui restreint l’utilisation de la loi aux chargements monotones avec faibles décharges.
Cette modélisation, principalement développée pour des besoins spécifiques à EDF, est avant tout conçue pour décrire le comportement mécanique monotone des aciers des cuves du parc nucléaire français de type 16MND5 (ou de type A508 classe 3) à l’état non-irradié et irradié.
Le modèle théorique d’écoulement#
La description théorique générale de l’écoulement cristallin de la ferrite est présentée en détails dans [bib1]. Cette dernière a été ensuite simplifiée pour décrire le cas des aciers de type 16MND5 dans [bib2]. Par ailleurs, les équations constitutives du comportement cristallin ont été intégrées analytiquement suivant l’hypothèse d’homogénéisation de Taylor (déformation homogène) présentée dans [bib3] pour aboutir à une description isotrope polycristalline du comportement de l’acier de cuve dans [bib4]. C’est ce dernier article qui constitue la référence de la loi décrite dans le présent document. Nous nous contentons ici de présenter les équations constitutives finales et de décrire les améliorations apportées au modèle depuis la publication de cet article.
Notations et rappels de viscoplasticité#
On rappelle que la viscoplasticité désigne la théorie de la mécanique des milieux continus qui permet de décrire le comportement anélastique d’un matériau dont la déformation dépend de la vitesse à laquelle sont appliqués les chargements sur le solide considéré. Il s’agit donc d’un cas particulier de comportement plastique, le matériau subissant des déformations irréversibles permanentes quand un certain niveau de chargement est atteint. Cependant, le critère de plasticité est ici également dépendant de l’histoire du chargement (et de sa vitesse), généralement au travers d’une variable d’état interne scalaire appelée déformation équivalente. Elle permet de mesurer la capacité du matériau à se durcir sous l’effet de sa déformation plastique (permanente) sans se rompre, indépendamment de la direction choisie. On parle d’écrouissage isotrope et le matériau est qualifié de ductile, comme c’est le cas pour de nombreux métaux et alliages.
Le modèle fait intervenir explicitement la température, notée \(T \left( ^{\circ} C \right)\) lorsque celle-ci est exprimée en degré Celcius ou \(T \left( K \right)\) en Kelvin. Dans le modèle final, on impose que toute température est exprimée en degré Celsius.
On se place ici sous l’hypothèse des petites déformations. En particulier, les déformations totales \(\tensTwo{\varepsilon}^{tot}\) utilisées dans la relation de comportement VISC_ISOT_PLAS sont les déformations linéarisées. Celles-ci se décomposent en une partie élastique \(\tensTwo{\varepsilon}^{elas}\) et une partie viscoplastique \(\tensTwo{\varepsilon}^{vp}\):
La partie élastique suit une loi de Hooke isotrope linéaire standard dépendant du module de Young \(E > 0\) et du coefficient de Poisson \(\nu \in \left] -1; \frac{1}{2} \right[\) asssociés au matériau utilisé. En introduisant le module de cisaillement \(\mu:= \frac{E}{2 \left( 1+\nu \right)}\), la contrainte \(\tensTwo{\sigma}\) s’exprime donc sous la forme:
où \(\text{Tr}\) désigne l’opérateur trace et \(\tensTwo{I}\) le tenseur identité d’ordre deux. On introduit aussi sa partie déviatorique \(\tensTwo{\tilde{\sigma}} : = \tensTwo{\sigma} - \frac{1}{3} \text{Tr} \left( \tensTwo{\sigma} \right) \tensTwo{I}\) ainsi que la contrainte équivalente (de type de von Mises) \(\sigma_{eq} := \sqrt{\frac{3}{2} \tensTwo{\tilde{\sigma}} : \tensTwo{\tilde{\sigma}} }\). Notons que le facteur de normalisation \(\frac{3}{2}\) permet de s’assurer que \(\sigma\) est bien égal à \(\sigma_{eq}\) dans le cadre d’une sollicitation uniaxiale (1D). On définit alors la fonction de charge \(f\) suivante:
où la fonction d’écrouissage isotrope \(R\) (présentée dans la section suivante) dépend à la fois de la vitesse de déformation viscoplastique cumulée \(\dot{p} : = \sqrt{\frac{2}{3} \tensTwo{\varepsilon}^{vp} : \tensTwo{\varepsilon}^{vp} }\) et de la déformation viscoplastique cumulée \(p(t) := \int_{0}^{t} \dot{p} \left( \tau \right) \, \text{d} \tau\), également appelée déformation équivalente. Remarquons que le facteur de normalisation \(\frac{2}{3}\) permet de s’assurer que \(\varepsilon^{vp}\) est bien égal à \(p\) dans le cadre uniaxial (1D). On mentionne aussi qu’une des spécificités de ce modèle d’écoulement est que la loi d’écrouissage isotrope \(R\) dépend de \(\dot{p}\) et pas seulement de \(p\), comme c’est le cas pour les lois d’écrouissage isotrope plus standards. On ne considère pas d’écrouissage cinématique dans ce modèle, ce qui restreint l’utilisation de la loi aux chargements monotones avec faibles décharges.
On suppose maintenant que la règle d’écoulement provient d’un potentiel viscoplastique et qu’elle est associative. Autrement dit, l’évolution de la déformation viscoplastique \(\tensTwo{\varepsilon}^{vp}\) est gouvernée l’équation d’écoulement et la direction d’écoulement est normale à la surface de charge:
où \(\tensTwo{n} := \frac{\partial f}{\partial \tensTwo{\sigma}} = \frac{3 \tensTwo{\tilde{\sigma}}}{2\sigma_{eq}}\) désigne la normale à la surface de charge et \(\lambda\) correspond au multiplicateur viscoplastique. Son évolution est donnée par les conditions de Karush-Kuhn-Tucker, aussi appelées conditions de cohérence:
Nous ne rentrons pas ici dans des considérations énergétiques qui dépassent la portée de ce document. Nous mentionnons simplement que sous certaines hypothèses de régularité, cette formulation est équivalente à un problème de maximisation de la dissipation intrinsèque, généralisant le principe du travail plastique maximal de Hill au cas de l’écrouissage isotrope. En particulier, on déduit de (2200) et du fait que \(\tensTwo{n}: \tensTwo{n} = \frac{3}{2}\) la relation suivante:
On vérifie aussi que \(\text{Tr} \left( \dot{\tensTwo{\varepsilon}}^{vp} \right) = 0\), c’est-à-dire que la déformation viscoplastique est isochore i.e. sans variation de volume. Il ne reste plus qu’à décrire la fonction d’écrouissage \(R\) utilisée dans ce modèle.
Composantes définissant la fonction d’écrouissage isotrope#
La fonction d’écrouissage isotrope \(R\) introduite dans la section précédente est constituée de plusieurs composantes, chacune étant associée à une composante de la microstructure ou à un mécanisme de la plasticité:
La signification des différentes composantes est rapidement donnée dans le tableau ci-dessous.
\(\tensTwo{\sigma}_{HP}\) |
Composante représentant l’effet de la taille de grain (également appelé effet Hall-Petch) |
\(\tensTwo{\sigma}_{VD}\) |
|
\(\tensTwo{\sigma}_{VS}\) |
|
\(\tensTwo{\sigma}_{D}\) |
|
\(\tensTwo{\sigma}_{E}\) |
|
\(\tensTwo{\sigma}_{I}\) |
|
\(\tensTwo{\sigma}_{auto}\) |
Résistance induite par l’intéraction des dislocations du même système |
\(\tensTwo{\sigma}_{LT}\) |
Résistance de la forêt (interactions avec les autres dislocations vis) |
Dans ce bilan, nous omettons la contribution des carbures intragranulaires. D’une part, leur densité est variable (bainite inférieure vs supérieure) et d’autre part, leur contribution reste faible comparée à celle des dislocations. Ainsi, leur éventuelle contribution peut être prise en compte en augmentant légèrement la valeur de la densité de dislocations initiale.
Les composantes \(\tensTwo{\sigma}_{D}\), \(\tensTwo{\sigma}_{VS}\) et \(\tensTwo{\sigma}_{VD}\) n’existaient pas dans le modèle initial décrit dans [bib4]. Elles apportent des contributions de deuxième ordre comparées aux contributions initiales.
Nous décrivons à partir de maintenant l’ensemble des expressions retenues pour les différentes contributions. Leurs versions finales sont exprimées dans le système international d’unités (SI) et une adaptation est effectuée lorsque ce n’est pas le cas pour les unités de longueur et contrainte.
L’effet de taille de grain \(\tensTwo{\sigma}_{HP}\) (aussi appelé effet Hall-Petch)#
Le terme \(\tensTwo{\sigma}_{HP}\) ne découle pas de la modélisation multi-échelles car les effets de taille de grains restent très mal compris. L’effet Hall-Petch (HP) est donc introduit à l’aide des résultats expérimentaux. L’effet a été confirmé et caractérisé dans les matériaux ferritiques [bib5]. Il peut être mis sous la forme:
où \(K_{HP}\) désigne la constante de Hall-Petch (fixée à \(0.45\, MPa\sqrt{m}\)) et \(d\) représente la taille moyenne des paquets de lattes exprimée en mètres, aussi appelé TAILLE_GRAIN dans les paramètres matériaux. La normalisation par le module de cisaillement \(\mu\) est introduite pour rendre compte de la dépendance à la température. Afin de garder une certaine cohérence avec [bib6], on utilise les valeurs de \(E \left( T^{\circ}C \right) := \left( 2.6 MPa \right) \left( 78821 - 23.721 \, T \left( ^{\circ}C \right) \right)\) et \(\nu = 0.3\) données par le RCC-M pour l’expression de \(\mu \left( 25 ^{\circ}C \right): = \frac{E \left( 25^{\circ}C \right)}{2 \left(1 + \nu \right)}\). On en déduit alors l’expression suivante pour la contribution dûe à l’effet Hall-Petch:
où la valeur de la contante \(K_{HP0}\) est donnée dans le tableau ci-dessous. Une adaptation des unités est effectuée lorsque ce ne sont pas les unités SI qui sont choisies pour \(\mu\) et TAILLE_GRAIN.
\(K_{HP0}\) |
\(3.3090 \, 10^{-11} \, m\) |
Les vieillissements dynamique \(\tensTwo{\sigma}_{VD}\) statique \(\tensTwo{\sigma}_{VS}\)#
Ces termes empiriques ont été introduits pour rendre compte du crochet de traction ainsi que de la légère consolidation de la limite d’élasticité au-delà de \(150\,^{\circ}C\). Leur contribution au comportement total est faible mais ils influencent la détermination de la limite d’élasticité. En particulier, le vieillissement statique \(\tensTwo{\sigma}_{VS}\) génère un phénomène d’adoucissement lors de la transition élasto-viscoplastique. Les deux expressions qui rendent compte des résultats expérimentaux sont:
et
où \(\left\langle \bullet \right\rangle_{+}\) désigne la partie positive de la grandeur et les valeurs des différentes contantes sont données dans les tableaux ci-dessous. Une adaptation des unités est effectuée lorsque ce ne sont pas les unités SI qui sont choisies pour \(\tensTwo{\sigma}_{VD}\) et \(\tensTwo{\sigma}_{VS}\).
\(A_{VD}\) |
\(2.5 \, 10^{7} \, Pa\) |
\(B_{VD}\) |
\(10^{-2} \, \left( ^{\circ}C \right)^{-1}\) |
\(C_{VD}\) |
\(2.0\) |
\(A_{VS}\) |
\(8.0 \, 10^{7} \, Pa\) |
\(B_{VS}\) |
\(2.0 \, 10^{5} \, Pa \, \left( ^{\circ}C \right)^{-1}\) |
\(C_{VS}\) |
\(150.0 \, ^{\circ}C\) |
\(D_{VS}\) |
\(200.0\) |
La résistance des crans \(\tensTwo{\sigma}_{D}\)#
Cette nouvelle composante traduit la basse mobilité des crans sur les dislocations vis engendrés par l’intersection avec les dislocations non-coplanaires. Leur mouvement n’est pas conservatif, c’est-à-dire qu’il requiert l’émission thermo-activée d’interstitiels ou de lacunes, ce qui rajoute une résistance supplémentaire au mouvement des dislocations vis déjà soumises à la friction du réseau (i.e le terme \(\tensTwo{\sigma}_{E}\)). L’expression initiale de \(\tensTwo{\sigma}_{D}\) en exponentielle d’Arrhenius est souvent simplifiée en une expression de type puissance:
où les valeurs des différentes constantes sont données dans le tableau ci-dessous. Une adaptation des unités est effectuée lorsque ce ne sont pas les unités SI qui sont choisies pour \(\tensTwo{\sigma}_{D}\).
\(\sigma_{D0}\) |
\(5.0 \, 10^{6} \, Pa\) |
\(\gamma_{D0}\) |
\(10^{-6} \, s^{-1}\) |
\(N_{D}\) |
\(10\) |
- Remarque
Il convient de préciser ici que l’unité de temps considérée est la seconde. Cette dépendance du modèle est apparue tardivement dans le processus d’intégration de la loi de comportement VISC_ISOT_PLAS et fera sûrement l’objet d’évolution future. Par conséquent, il est important de noter que pour le moment, le modèle impose la seconde pour l’unité gérant l’évolution des chargements. Dans le cas contraire, cela peut amener à des résultats incorrects.
La contrainte effective \(\tensTwo{\sigma}_{E}\)#
Suivant la discussion dans [bib4], la contrainte effective dans l’acier de cuve (caractérisée par une grande densité de dislocations initiale) peut être donnée par l’équation:
avec \(\left\langle \bullet \right\rangle_{+}\) la partie positive de la grandeur, \(T \left( K \right) = T \left( ^{\circ}C \right) + 273.15\) la température exprimée en Kelvin et la température d’activation
où les constantes sont ici l’énergie d’activation \(\Delta G_{0} = 0.84 \, eV\), la constante de Boltzmann \(k_{B} = 8.62 \, 10^{-5} eV \, K^{-1}\), la densité de dislocations mobiles \(\rho_{m} = 10^{13} \, m^{-2}\), la norme du vecteur de Burgers \(b = 2.48 \, 10^{-10} \, m\), le facteur de fréquence \(F = 10^{13} \, s^{-1}\) et la taille critique \(l_{c} = 10^{-8} \, m\). C’est cette composante \(\tensTwo{\sigma}_{E}\) qui porte la sensibilité principale à la vitesse de déformation et à la température. On l’exprime ici sous la forme suivante:
où les valeurs des différentes contantes sont données dans le tableau ci-dessous. Une adaptation des unités est effectuée lorsque ce ne sont pas les unités SI qui sont choisies pour \(\tensTwo{\sigma}_{E}\).
\(\sigma_{E0}\) |
\(9.0 \, 10^{8} \, Pa\) |
\(T_{0}\) |
\(273.15 \, ^{\circ} C\) (i.e le zéro absolu) |
\(A_{E}\) |
\(9744.78 \, ^{\circ} C\) |
\(B_{E}\) |
\(19.3289\) |
\(C_{E}\) |
\(1.0 \, s\) |
- Remarque
Il convient de préciser ici que l’unité de temps considérée est la seconde. Cette dépendance du modèle est apparue tardivement dans le processus d’intégration de la loi de comportement VISC_ISOT_PLAS et fera sûrement l’objet d’évolution future. Par conséquent, il est important de noter que pour le moment, le modèle impose la seconde pour l’unité gérant l’évolution des chargements. Dans le cas contraire, cela peut amener à des résultats incorrects.
Contribution des défauts d’irradiation \(\tensTwo{\sigma}_{I}\)#
En plus des composantes de la microstructure, l’irradiation neutronique génère d’autres défauts provoquant le durcissement de l’irradiation. Les principaux défauts identifiés sont les amas de solutés, révélés par la sonde atomique tomographique et les boucles de dislocations (souvent confondues avec le dommage de la matrice). Les analyses des résultats de la sonde et de la microscopie électronique [bib7] montrent que le dommage de la matrice augmente rapidement et sature à une fluence inférieure à \(10^{23} \, n \, m^{-2}\), alors que la concentration des amas de solutés C_AMAS, dont le diamètre TAILLE_AMAS est d’environ \(3 \, nm\), augmente proportionnellement avec la fluence.
Pour simplifier, nous considérons que seuls les amas de soluté existent dans le matériau et que le durcissement d’irradiation est seulement fonction des paramètres matériaux C_AMAS et TAILLE_AMAS. La contribution des amas a été déterminée par les simulations de dynamique des dislocations et la dynamique moléculaire, ce qui a permis d’établir une équation analytique dans [bib7]. Cependant, l’équation est relativement complexe et nécessite des précautions numériques. Or, une simple expression empirique s’est avérée suffisamment robuste pour reproduire les valeurs de l’équation analytique. C’est donc cette expression numérique qui a été finalement retenue pour la loi:
lorsque \(C_{L} \, C\_AMAS > A_{I}\) et sinon \(\tensTwo{\sigma}_{I} = 0\). Les valeurs des différentes contantes sont données dans le tableau ci-dessous. Une adaptation des unités est effectuée lorsque ce ne sont pas les unités SI qui sont choisies pour \(\mu\), TAILLE_AMAS et C_AMAS.
\(C_{L}\) |
\(10^{-22} \, m^{3}\) |
\(A_{I}\) |
\(10^{-2}\) |
\(B_{I}\) |
\(3.35 \, 10^{-8}\) |
\(C_{I}\) |
\(10^{9} \, m^{-1}\) |
\(D_{I}\) |
\(2.3\) |
\(E_{I}\) |
\(8\) |
\(F_{I}\) |
\(7\) |
L’auto-interaction \(\tensTwo{\sigma}_{auto}\) liée à l’évolution de la densité de dislocations#
Une des originalités du modèle consiste à diviser la population de dislocations en quatre groupes formés de dislocations vis, au lieu des douze systèmes usuels dans la plasticité cristalline. En effet, comme la seule variabilité pour une dislocation vis est sa direction de ligne (i.e. le vecteur de Burgers \(\vector{b}\)), il n’existe que quatre directions différentes de \(\vector{b}\): \(\left[ 1 1 1 \right]\), \(\left[ \bar{1} 1 1 \right]\), \(\left[ 1 \bar{1} 1 \right]\) et \(\left[ 1 1 \bar{1} \right]\).
L’interaction entre les dislocations vis parallèles est l’auto-interaction (de l’anglais self-interaction). Elle correspond donc à l’interaction avec un quart de la population, d’où la formule:
où \(\mu\) désigne toujours le module de cisaillement, \(b\) la norme du vecteur de Burgers \(\vector{b}\) (ici fixée à \(2.48 \, 10^{-10} \, m\)) et où \(\rho\) représente la densité (linéique) de dislocations dans le matériau dont l’évolution convient maintenant d’être décrite. Cette description diffère du modèle initial [bib4] et permet de mieux prendre compte les phénomènes liés à de grandes déformations (seulement disponible au formalisme “GDEF_LOG” pour ce comportement). C’est pourquoi nous la détaillons.
Ainsi, comme on vient de le voir, les interactions entre dislocations dépendent de la densité de dislocations \(\rho\). En raison du stockage et de la nécessité de créer de nouvelles dislocations mobiles, cette densité augmente avec la déformation plastique cumulée \(p\) conduisant à l’écrouissage du matériau.
Une importante évolution du modèle initial [bib4] consiste à introduire un nouveau terme représentant le stockage de dislocations non-vis, qui s’ajoute au stockage initial des dislocations vis.
Pendant la déformation plastique, en plus du stockage classique, les dislocations mobiles coupent les autres dislocations en formant des crans (jogs en anglais) non-vis immobiles. A faible déformation, la densité de ces crans est très faible et c’est pourquoi elle a été négligée dans le modèle initial [bib4]. Cependant, à grande déformation, leur densité devient prépondérante car ils ne peuvent pas s’annihiler facilement: on parle d’annihilation thermiquement activée. On note cette densité \(\rho_{j}\) (\(j\) pour jog), à distinguer de la densité des dislocations piégées et stockées après leur mouvement notée \(\rho_{s}\). L’évolution des deux densités donne celle de la densité totale \(\rho\):
Selon la théorie de Kocks-Mecking, les densités de dislocations \(\rho_{s}\) et \(\rho_{j}\) s’immobilisent après avoir parcouru une distance notée \(\lambda_{s}\) et \(\lambda_{j}\), respectivement. Néanmoins, les deux densités saturent lorsque la distance moyenne entre les deux types de dislocations devienne inférieure à \(y_{s}\) et \(y_{j}\), respectivement. Par simplicité, nous normalisons toutes les distances par la norme du vecteur de Burgers \(b\). Nous effectuons ainsi les normalisations suivantes: \(L_{s}:= \frac{\lambda_{s}}{b}\), \(L_{j}:= \frac{\lambda_{j}}{b}\), \(f_{s}:= \frac{y_{s}}{b}\) et \(f_{j}:= \frac{y_{j}}{b}\). L’évolution des deux populations étant indépendantes, nous pouvons donc écrire:
avec \(M\) le facteur de Taylor (que nous considérons égal à \(3\)). L’intégration analytique de chaque équation d’évolution fournit les deux expressions suivantes pour les densités:
où les constantes \(A_{s}\) et \(A_{j}\) doivent être déterminées à partir des conditions initiales. D’une part, on fait l’hypothèse qu’au début de la déformation, la densité de crans \(\rho_{j}\) est nulle, ce qui permet de déterminer \(A_{j}\):
D’autre part, on suppose que la densité initiale de dislocations \(\rho\) est donnée par le paramètre matériaux D_DISLOC, ce qui fournit la valeur de \(A_{s}\):
Par ailleurs, comme le mouvement des crans est restreint, nous considérons qu’ils ne peuvent se déplacer pour s’annihiler. Par simplicité, nous supposons donc que \(M \, f_{j} \approx 1\). On en déduit que la densité totale s’écrit:
On simplifie encore davantage cette expression en faisant l’hypothèse qu’à très grande déformation, la densité de crans sur chaque système ne peut pas dépasser la densité de saturation des dislocations vis. Comme il y a trois types de crans pour chaque dislocation vis (intersections avec les trois autres dislocations vis), nous considérons que la densité finale de cran \(\frac{M}{b^{2} L_{j}}\) est approximativement égale à trois fois la densité des dislocations vis. On en déduit que \(L_{j} = L_{s} f_{s}\) et on réduit l’expression précédente à sa forme finale:
On peut désormais omettre les indices en renotant \(L\) la distance libre moyenne \(L_{s}\) (normalisée par \(b\)) ainsi que \(f\) la distance d’annihilation \(f_{s}\) (aussi normalisée par \(b\)) dont les expressions seront bientôt spécifiées. En renotant \(A_{auto}\) le facteur de Taylor \(M\) et en introduisant \(B_{auto} : = \frac{b^{2}}{4}\), il en résulte l’expression finale de l’auto-intéraction \(\tensTwo{\sigma}_{auto} = \mu \, b \, \sqrt{\frac{\rho}{4}}\) utilisée pour la loi VISC_ISOT_PLAS:
avec les expressions suivantes retenues pour la distance libre moyenne \(L\) et d’annihilation \(f\) (toutes deux normalisées)
Les valeurs des différentes contantes sont données dans le tableau ci-dessous. Une adaptation des unités est effectuée lorsque ce ne sont pas les unités SI qui sont choisies pour \(\mu\), D_DISLOC et C_AMAS.
\(A_{auto}\) |
\(3\) |
\(B_{auto}\) |
\(1.5376 \, 10^{-20} \, m^{2}\) |
\(A_{L}\) |
\(2.222222222 \, 10^{-4}\) |
\(B_{L}\) |
\(6.7925e \, 10^{-6}\) |
\(C_{L}\) |
\(10^{-22} \, m^{3}\) |
\(A_{f}\) |
\(6.8\) |
\(B_{f}\) |
\(1.178\) |
\(C_{f}\) |
\(0.019 \, ^{\circ} C^{-1}\) |
\(T_{0}\) |
\(273.15 \, ^{\circ} C\) (i.e. le zéro absolu) |
- Remarque
On remarquera de ce qui précède que \(A_{auto}\) n’est rien d’autre que le facteur de Taylor et que \(B_{auto} = \frac{b^{2}}{4}\), c’est-à-dire le quart du carré de la norme du vecteur de Burgers. Enfin, on précise également que \(A_{L} = \frac{b}{\lambda_{0}}\) avec \(\lambda_{0} = 1.116 \, \mu m\) le libre parcours moyen dans le matériau non-irradié.
Contribution de la tension de ligne \(\tensTwo{\sigma}_{LT}\) aussi appelée résistance de la forêt#
L’interaction avec les dislocations de la forêt (i.e. différents vecteurs de Burgers \(\vector{b}\)) doit rendre compte de l’interaction avec toutes les autres dislocations. La contribution maximale est donc de \(\mu b \sqrt{\frac{3 \rho}{4}} = \sqrt{3} \tensTwo{\sigma}_{auto}\). Cependant, comme expliqué dans [bib1], la contribution réelle diminue avec la contrainte effective \(\tensTwo{\sigma}_{E}\) en raison de la courbure des dislocations non-vis qui diminue la résistance des jonctions. Cette interaction \(\tensTwo{\sigma}_{LT}\) (long-term en anglais) peut alors être exprimée par la formule suivante:
où on rappelle que \(\left\langle \bullet \right\rangle_{+}\) désigne la partie positive de la grandeur. On décrit dans la section suivante l’implémentation numérique de ce modèle continu.
Implémentation numérique#
L’implémentation numérique du modèle continu décrit dans la section précédente est effectuée en utilisant le module d’intégration de lois de comportement MFront [bib8]. On se referra à la notice d’utilisation [U2.10.02] pour davantage de détails concernant son couplage avec code_aster ainsi qu’au fichier VISC_ISOT_PLAS.mfront (dispnible dans les sources) pour les détails sur l’écriture de la loi de comportement au formalisme MFront.
Tout d’abord, on décrit la formulation retenue pour discrétiser le problème. Au sein de la fonction d’écrouissage isotrope \(R\) définie dans (2201), on commence par isoler le terme \(\tensTwo{\sigma}_{D}\) associé à la résistance des crans en introduisant un pseudo-écrouissage \(\tensTwo{\sigma}_{C} := R - \tensTwo{\sigma}_{D}\) puis en lui associant une pseudo-surface de charge:
On inverse alors la relation puissance (2202) définissant \(\tensTwo{\sigma}_{D}\) pour obtenir une expression de la vitesse de déformation viscoplastique équivalente \(\dot{p}\) sous la forme d’une formulation de type Perzyna:
où \(\left\langle \bullet \right\rangle_{+}\) désigne la partie positive de la grandeur. Ainsi, dans l’implémentation finale de la loi, c’est la différence \(\sigma_{eq} - \tensTwo{\sigma}_{C}\) qui va induire la vitesse de déformation plastique, bien que cette vitesse intervienne également dans la contrainte \(\tensTwo{\sigma}_{C}\).
La directive @DSL (pour Domain Specific Language en anglais) utilisée dans MFront pour décrire la loi de comportement VISC_ISOT_PLAS est le DSL ImplicitII. Celui-ci est identique au DSL Implicit sauf qu’il ne prédéfinit pas le tenseur des déformations élastiques comme première variable d’état interne. En effet, comme on va le montrer, la seule réelle inconnue ici à calculer sera l’incrément de déformation viscoplastique équivalente \(\Delta p\).
- Remarque
Du fait de la dépendance de notre pseudo-écrouissage \(\tensTwo{\sigma}_{C}\) à la vitesse de déformation viscoplastique équivalente \(\dot{p}\), nous n’avons pas été en mesure de choisir directement: ni le DSL IsotropicMisesCreep qui traite le cas \(\dot{p} = \phi \left( \sigma_{eq} \right)\); ni le DSL IsotropicStrainHardeningMisesCreep qui traite le cas \(\dot{p} = \phi \left( \sigma_{eq}, p \right)\); ni le DSL IsotropicPlasticMisesFlow qui traite le cas \(\phi \left( \sigma_{eq}, p \right) \leq 0\) avec \(\dot{p} \geq 0\) et \(\dot{p} \, \phi \left( \sigma_{eq}, p \right) = 0\); ni même la StandardElastoViscoplasticity brick.
On utilise une formulation totalement implicite i.e que le paramètre MFront prédéfini \(\theta\) est fixé à \(1\). De plus, l’algorithme de résolution choisi est la méthode de Newton-Raphson et la jacobienne est calculée analytiquement. Une vérification des types d’unités physiques utilisées pour décrire la loi de comportement est réalisée à la compilation du fichier MFront grâce à l’activation de la directive @Qt (pour Quantity type en anglais).
La stratégie de résolution de la loi de comportement sur un intervalle de temps \(\left[ t, t + \Delta t \right]\) est celle classiquement effectuée dans MFront et décrite dans [bib9]. On adapte juste celle-ci à notre cas. On suppose donc que l’incrément de déformations totales \(\Delta \tensTwo{\varepsilon}^{tot} : = \tensTwo{\varepsilon}^{tot} \left( t + \Delta t \right) - \tensTwo{\varepsilon}^{tot} \left( t \right)\) est connu.
On commence par initialiser les variables locales grâce à la directive @InitLocalVariables. Puis, avec la directive @Predictor, on fait une prédiction élastique \(\sigma_{eq}^{pred}\) et \(\tensTwo{n}^{pred} := \frac{3 \tensTwo{\tilde{\sigma}}^{pred}}{2 \sigma_{eq}^{pred}}\) en supposant que l’incrément de déformations élastiques \(\Delta \tensTwo{\varepsilon}^{elas}\) est égal à \(\Delta \tensTwo{\varepsilon}^{tot}\) (ou autrement dit que l’incrément de déformations viscoplastiques \(\Delta \tensTwo{\varepsilon}^{vp}\) est nul):
Si l’on vérifie \(\sigma_{eq}^{pred} - \tensTwo{\sigma}_C \left( p \left( t \right), 0 \right) < 0\), alors on est bien dans le régime élastique et il n’y a rien à faire de plus: \(\Delta p = 0\), \(\Delta \tensTwo{\varepsilon}^{vp} = 0\) et la loi de Hooke fournit l’incrément de contrainte \(\Delta \tensTwo{\sigma}\).
Dans le cas contraire, il faut corriger l’incrément de déformations élastiques \(\Delta \tensTwo{\varepsilon}^{elas}\) par un incrément de déformation viscoplastique:
la dernière égalité étant une approximation dont on fait l’hypothèse à partir de maintenant. En suivant le raisonnement de [bib9], on peut montrer qu’on a alors les relations suivantes: \(\tensTwo{n} \left( t + \Delta t \right) = \tensTwo{n}^{pred}\) et \(\sigma_{eq} \left( t + \Delta t \right) = \sigma_{eq}^{pred} - 3 \mu \Delta p\). On en déduit que la relation (2203) prise au temps \(t + \Delta t\) se discrétise sous la forme suivante:
En dénotant \(F\) le membre de gauche précédent, la relation précédente se réécrit sous la forme \(F \left( \Delta p \right) = 0\). On est donc amené à résoudre une équation scalaire d’une inconnue scalaire \(\Delta p\) en utilisant la méthode de Newton-Raphson via la directive @Integrator et que l’on détaillera dans la section suivante. Une fois que le zéro de la fonction \(F\) a été trouvé, on peut alors calculer l’incrément de contrainte grâce à la loi de Hooke et via la directive @ComputeFinalStress:
où \(\vector{\tensTwo{I}}\) désigne le tenseur identité d’ordre quatre. Finalement, grâce à la directive @TangentOperator, on calcule analytiquement l’opérateur tangent \(D_{t} : = \frac{ \partial \left( \Delta \tensTwo{\sigma} \right) }{ \partial \left( \Delta \tensTwo{\varepsilon}_{tot} \right)}\) en dérivant l’expression précédente. En effet, le fait d’avoir une seule inconnue scalaire fournit directement l’inversion de la jacobienne
tandis que le calcul des quantités \(\frac{ \partial \sigma_{eq}^{pred}}{\partial \left( \Delta \tensTwo{\varepsilon}^{tot} \right)}\) et \(\frac{ \partial \tensTwo{n}^{pred} }{ \partial \left( \Delta \tensTwo{\varepsilon}^{tot} \right)}\) est standard et détaillé dans [bib9].
On rappelle que la méthode de Newton-Raphson est une méthode itérative permettant de trouver le zéro d’une fonction \(F\) à partir d’un point de départ \(\Delta p^{0}\) et en posant successivement:
où on rappelle que \(F\) est ici explicitement défini par le membre de gauche de (2204) et que sa dérivée peut se calculer analytiquement par dérivation. On vérifie au passage qu’il s’agit d’une fonction continue mais seulement \(C^{1}\) par morceaux du fait des parties positives intervenant dans de nombreux termes de \(\tensTwo{\sigma}_{C}\).
Bien qu’ici, le zéro semble correctement défini de manière unique, il est connu que la méthode de Newton-Raphson peut ne pas converger dans ce genre de situation. En particulier, le choix d’un bon point de départ est nécessaire. De plus, il s’avère que \(\Delta p^{0} = 0\) est le pire puisqu’il conduit forcément à une non-convergence de l’algorithme.
Une méthode de la sécante aurait été plus adaptée dans notre cas (au prix de faire passer l’ordre de convergence quadratique à \(1.6\) environ) mais on a souhaité rester dans le formalisme MFront pour résoudre cette équation et tordre un peu le cou à la méthode de Newton-Raphson. C’est d’ailleurs exactement ce que font certains algorithmes internes de Code_Aster. On va donc choisir notre point de départ à partir d’un intervalle encadrant le zéro qu’on cherche.
On pose d’abord \(\Delta p^{try} = 10^{-11}\). Puis, on distingue deux cas dans la directive @Predictor. D’une part, si \(F \left( \Delta p^{try} \right) < 0\), alors on pose \(\Delta p^{left} = 10^{n_{0}-1} \Delta p^{try}\) et \(\Delta p^{right} = 10^{n_{0}} \Delta p^{try}\) avec \(n_{0}\) le plus petit entier strictement positif pour lequel \(F \left( \Delta p^{right} \right) > 0\). D’autre part, si \(F \left( \Delta p^{try} \right) > 0\), alors on pose \(\Delta p^{left} = 10^{-n_{0}} \Delta p^{try}\) et \(\Delta p^{right} = 10^{1-n_{0}} \Delta p^{try}\) avec \(n_{0}\) le plus petit entier strictement positif pour lequel \(F \left( \Delta p^{left} \right) < 0\).
A la fin de ce processus qui ne nécessite que des évaluations d’une fonction scalaire, on dispose d’un intervalle \(\left[ \Delta p^{left}, \Delta p^{right} \right]\) qui encadre le zéro que l’on cherche. On définit alors notre point de départ pour la méthode de Newton-Raphson \(\Delta p^{0} : = \Delta p^{left}\).
Cependant, on va demander également à l’intervalle d’être stable pour la méthode de Newton-Raphson, c’est-à-dire que le \(\Delta p^{n+1}\) défini par (2205) doit rester dans l’intervalle \(\left[ \Delta p^{left}, \Delta p^{right} \right]\). Durant la directive @Predictor, tant que cette condition n’est pas vérifiée, on utilise une dichotomie \(\Delta p^{n+1}=\frac{\Delta p^{left} + \Delta p^{right}}{2}\) à la place en réduisant la taille de l’intervalle jusqu’à pouvoir basculer sur une méthode de Newton-Raphson. Cependant, même lorsqu’on a basculé sur la méthode de Newton-Raphson, on continue de réduire la taille de l’intervalle à chaque itération en distinguant les cas \(F \left( \Delta p^{n+1} \right)\) positif et négatif. Autrement dit, on pose \(\Delta p^{left} = \Delta p^{n+1}\) si \(F \left( \Delta p^{n+1} \right) < 0\) ou bien \(\Delta p^{right} = \Delta p^{n+1}\) si \(F \left( \Delta p^{n+1} \right) > 0\) de sorte que l’intervalle \(\left[ \Delta p^{left}, \Delta p^{right} \right]\) encadre toujours le zéro que l’on cherche mais que sa taille diminue au fur et à mesure des itérations successives de l’algorithme.
Le critère d’arrêt du DSL ImplicitII est que la condition :math:` vert F left( Delta p^{n} right) vert < epsilon` où le paramètre MFront prédéfini \(epsilon\) est fixé à \(5.0 \, 10^{-14}\) et correspond au mot-clé RESI_INTE. De même, le nombre maximal d’itérations dans l’algorithme de Newton-Raphson est fixé à \(100\) via le paramètre MFront prédéfini \(iterMax\) correspondant au mot-clé ITER_INTE_MAXI. Dans quelques rares cas, l’extrême sensibilité de la fonction \(F\) vis-à-vis de sa variable \(\Delta p\) fait que le critère d’arrêt n’est pas vérifié mais qu’on a atteint l’epsilon-machine pour \(\Delta p\) (et donc qu’on a convergé). Le fait de réduire à chaque itération la taille de l’intervalle permet de repérer et d’éliminer ces cas sans ajouter de critère d’arrêt supplémentaire. De même, il arrive dans de très rares cas que l’algorithme oscille entre deux valeurs très proches, dont la différence est de l’ordre de l’espilon-machine. De façon similaire, la taille de l’intervalle encadrant le zéro permet de gérer ce genre de situations. Pour finir, le choix du critère d’arrêt à \(5.0 \, 10^{-14}\) est établi par une heuristique issue de l’utilisation de cette loi et qui exclut de manière significative les deux phénomènes d’espilon-machine décrits précédemment.
Bibliographie#
Monnet G., Vincent L. and Devincre B., Dislocation-dynamics based crystal plasticity law for the low- and high-temperature deformation regimes of bcc crystal. Acta Materialia, 61(16) (2013): 6178‑6190.
Monnet G., Vincent L., and Gélébart L., Multiscale modeling of crystal plasticity in Reactor Pressure Vessel steels: Prediction of irradiation hardening, Journal of Nuclear Materials, 514 (2019): 128.
Kocks U. F., Tomé C. N., and Wenk H.-R., Texture and Anisotropy: Preferred Orientations in Polycrystals and Their Effect on Materials Properties. Cambridge University Press, 2000.
Monnet G., Analytical flow equation for irradiated low-alloy steels established by multiscale modeling. Journal of Nuclear Materials, 586 (2023).
Nishino K., and Honma H., A Study on Temperature and Strain Rate Dependences of Yield Stress at Low Temperatures in Pure Iron and Fe-Cr Alloys. Tetsu-to-Hagane, 57(6) (1971): 954‑964.
RCC-M : Regles de Conception et de Construction des materiels mecaniques des ilots nucleaires PWR. Edition juin 2000, modificatif juin 2007. Edite par l’AFCEN : Association francaise pour les regles de conception et de construction des materiels des chaudieres electro-nucleaires.
Monnet G., Multiscale modeling of irradiation hardening: Application to important nuclear materials. Journal of Nuclear Materials (2018).
CEA et EDF, MFront web site https://thelfer.github.io/tfel/web/index.html.
Helfer, E. Castelier, V. Blanc, J. Julien, Le générateur de code mfront: écriture de lois de comportement mécanique (2013). https://thelfer.github.io/tfel/web/documents/mfront/behaviours.pdf