r5.03.13 Comportement viscoplastique avec endommagement de HAYHURST#
Résumé :
Le modèle visco-plastique couplé à l’endommagement isotrope dit de Hayhurst est particulièrement adapté pour réaliser des calculs de structures en fluage. La partie viscoplastique du modèle a été proposée par Hayhurst et al dans [bib1] tandis que la loi d’endommagement a été proposée par Charles Pétry dans [bib2]. Ce modèle permet une prédiction satisfaisante de la chute de ductilité en fluage via l’application de critères limites sur la déformation et sur l’endommagement. Ce modèle a pour l’instant été principalement utilisé à EDF R&D/MMC pour des prédictions de durée de vie en fluage sur l’acier de grade 92. Via l’identification de paramètres spécifiques et la réalisation de calculs de structure, ce modèle permet également de prédire de façon satisfaisante le comportement en fluage et la durée de vie de jonctions soudées [bib3].
Ce modèle est implanté dans Code_Aster sous le nom de HAYHURST ; les équations en vitesse sont intégrées numériquement par un schéma explicite de Runge-Kutta d’ordre 2 avec découpage automatique en sous-pas locaux en fonction d’une estimation de l’erreur d’intégration (méthode de Runge-Kutta emboîtée, confer [R5.03.14]) ou par une méthode d’intégration implicite de Newton.
Le test SSNV225 valide l’intégration de ce modèle et est présenté dans le document de validation [V6.04.225] qui fournit également des références expérimentales en rapport avec le cas test.
Formulation du modèle#
Cadre théorique#
Dans la formulation initiale proposée dans [bib1], deux variables d’endommagement distinctes possédant chacune une cinétique propre sont proposées. Une variable \(\phi\) est notamment associée aux évolutions de la microstructure dépendant uniquement du temps (i.e. vieillissement statique du matériau), alors que la variable \(\omega\) décrit les mécanismes de cavitation se développant sous l’influence combinée de la déformation viscoplastique et de la triaxialité des contraintes.
Dans la modélisation retenue dans Code_Aster, la variable \(\phi\) est conservée. En pratique, son identification est délicate, et il est possible de ne pas faire intervenir le vieillissement microstructural en mettant à zéro le coefficient \({k}_{c}\) (cf. équations en section 2.2). Par ailleurs, la variable \(\omega\) est renommée \(D\) car sa loi d’évolution est différente de celle proposée par Hayhurst (cf. [bib1]): la loi est ici en sinus hyperbolique. Les avantages de cette formulation sont détaillés dans [bib2].
Équations du modèle#
Les équations du modèles s’écrivent :
\(\lbrace \begin{array}{c}\mathit{élasticité}:\\ \sigma =(1-D)C{\varepsilon}^{e}\text{et}{\varepsilon}^{e}=\varepsilon -{\varepsilon}^{\mathit{th}}-{\varepsilon}^{p}\\ \text{}\\ \mathit{viscoplasticité}:\\ \dot{{\varepsilon}^{p}}=\frac{3}{2}\dot{p}\frac{\tilde{\sigma}}{{\sigma}_{\mathit{eq}}}\text{avec}\dot{p}={\varepsilon}_{0}\sinh(\frac{{\sigma}_{\mathit{eq}}(1-H)}{K(1-D)(1-\phi )})\\ \dot{\phi }=\frac{{k}_{c}}{3}{(1-\phi )}^{4}\\ \text{}\\ \mathit{écrouissage}:\\ H={H}_{1}+{H}_{2}\\ \dot{{H}_{i}}=\frac{{h}_{i}}{{\sigma}_{\mathit{eq}}}({H}_{i}^{\text{*}}-{\delta}_{i}{H}_{i})\dot{p}\text{pour}i=1,2\\ \text{}\\ \mathit{endommagement}:\\ \text{si}{\alpha}_{\sigma}=0\dot{D}={A}_{0}\sinh(\frac{{\alpha}_{D}\text{<}{\sigma}_{I}{\text{>}}_{+}+{\sigma}_{\mathit{eq}}(1-{\alpha}_{D})}{{\sigma}_{0}})\\ \text{si}{\alpha}_{\sigma}=1\dot{D}={A}_{0}\sinh(\frac{{\alpha}_{D}\text{<}\mathit{tr}(\sigma ){\text{>}}_{+}+{\sigma}_{\mathit{eq}}(1-{\alpha}_{D})}{{\sigma}_{0}})\end{array}\)
où:
\(\varepsilon\) , \({\varepsilon}^{\text{e}}\) , \({\varepsilon}^{\text{th}}\) et \({\varepsilon}^{\text{p}}\) |
sont respectivement les déformations totale, élastique, thermique et plastique, |
\(\text{<}x{\text{>}}_{+}\) |
est la partie positive de \(x\) , |
\({\sigma}_{1}\) |
est la contrainte principale maximale, |
\(\tilde{\sigma}=\sigma –\frac{1}{3}\mathit{Tr}(\sigma )I\) |
est la partie déviatorique du tenseur des contraintes, |
\({\sigma}_{\mathit{eq}}=\sqrt{\frac{3}{2}\tilde{{\sigma}_{ij}}\tilde{{\sigma}_{ij}}}\) |
est la contrainte déviatorique de Von-Mises, |
\(C\) |
est le tenseur de rigidité élastique, |
\(p\) |
est la déformation plastique cumulée, |
\(H\) , \({H}_{1}\) , \({H}_{2}\) |
sont les variables d’écrouissage isotrope viscoplastique, |
\(D\) |
est la variable scalaire d’endommagement isotrope, |
\(\phi\) |
est la variable scalaire d’endommagement microstructural, |
\({\alpha}_{\sigma}\) |
est le paramètre permettant de choisir de calculer l’endommagement par rapport à \(({\sigma}_{\mathit{eq}},{\sigma}_{1})\) pour \({\alpha}_{\sigma}=0\) ou \(({\sigma}_{\text{eq}},\text{Tr}(\sigma ))\) pour \({\alpha}_{\sigma}=1\) , |
\({\alpha}_{\text{D}}\) |
est le paramètre permettant d’ajuster la sensibilité à la triaxialité \(({\alpha}_{\text{D}}=1)\) ou la sensibilité à la contrainte principale maximale \(({\alpha}_{\text{D}}=0)\) pour le calcul de l’endommagement, |
\({\delta}_{i}\) |
vaut 0 ou 1 selon que l’on souhaite un écrouissage isotrope linéaire ou non-linéaire, respectivement. |
Nota Bene :
Les paramètres du modèle \(K,{\epsilon}_{0,}{\sigma}_{0,}{h}_{1,}{h}_{2,}{A}_{0,}{\alpha}_{D},\text{et}{k}_{c}\) peuvent être des fonctions de la température (en \(°C\) ). Dans [bib2], les paramètres identifiés varient en fonction de la température suivant une loi d’Arrhénius.
Afin d’empêcher l’évolution de l’endommagement, l’utilisateur peut employer une valeur de \({A}_{0}\) nulle. Dans ce cas, il s’agit malgré tout de renseigner une valeur pour \({\sigma}_{0}\) , qui intervient dans le sinus hyperbolique, valeur non nulle et suffisamment grande.
Remarque:
Le système d’équations précédent peut être réduit : en effet, celles qui sont relatives à l’évolution de l’écrouissage s’intègrent de la façon suivante:
\({H}_{i}=\frac{{H}_{i}^{\text{*}}}{{\delta}_{i}}\left[1-\exp(\frac{-{h}_{i}{\delta}_{i}}{{\sigma}_{\mathit{eq}}}p)\right]\)
et l’équation relative à l’endommagement micro-structural revient à :
\(\varphi =1-\frac{1}{{(1+{k}_{c}t)}^{1/3}}\)
C’est cette expression qui sera utilisée dans la suite du document.
Paramètres de la loi#
Les paramètres matériau requis pour l’utilisation du modèle dans Code_Aster via la commande DEFI_MATERIAU (cf. doc U4.43.01) sont les suivants:
ASTER |
Symbole |
Définition |
EPS0 |
\({\varepsilon}_{0}\) |
Paramètre agissant sur la cinétique de déformation viscoplastique |
K |
\(K\) |
Paramètre régissant le comportement en sinus hyperbolique de la loi viscoplastique |
H1 |
\({h}_{1}\) |
Module d’écrouissage |
H2 |
\({h}_{2}\) |
Module d’écrouissage |
DELTA1 |
\({\delta}_{1}\) |
Choix du type d’écrouissage |
DELTA2 |
\({\delta}_{2}\) |
Choix du type d’écrouissage |
H1ST |
\({H}_{1}^{\text{*}}\) |
Valeur à saturation de l’écrouissage, dans le cas non-linéaire |
H2ST |
\({H}_{2}^{\text{*}}\) |
Valeur à saturation de l’écrouissage, dans le cas non-linéaire |
BIGA |
\({A}_{0}\) |
Paramètre agissant sur la cinétique d’endommagement |
SIG0 |
\({\sigma}_{0}\) |
Paramètre régissant le comportement en sinus hyperbolique de la loi d’endommagement |
ALPHAD |
\({\alpha}_{\text{D}}\) |
Coefficient jouant sur la contrainte effective pour le calcul de l’endommagement |
KC |
\({k}_{c}\) |
Paramètre régissant la cinétique d’endommagement microstructural |
S_EQUI_D |
\({\alpha}_{\sigma}\) |
Choix de contrainte hydrostatique ou principale maximale |
Les paramètres \({\alpha}_{\text{D}}\) , \({k}_{c}\) , et \({\alpha}_{\sigma}\) sont optionnels et ont des valeurs nulles par défaut.
L’identification des coefficients de la partie viscoplastique du modèle est réalisée à partir d’essais de fluage à différents niveaux de contrainte.
Une fois cette identification réalisée, les coefficients pilotant l’endommagement peuvent être identifiés à partir d’essais de fluage longs à faible contrainte pour lesquels une chute de ductilité est observée (chute des déformations à rupture).
Un exemple type d’identification est détaillé dans [bib2].
Implantation dans Code_Aster#
L’utilisation du modèle HAYHURST est possible en modélisations 3D, axisymétrique, et déformation plane ( 3D , AXIS , D_PLAN , respectivement)
Les modèles de déformation possibles sont PETIT , PETIT_REAC , GDEF_HYPO_ELAS , GDEF_LOG.
L’algorithme utilisé est du type global-local. Les itérations globales utilisent la matrice de rigidité élastique calculée à partir de la matrice de Hooke endommagée : \(\underline{\underline{\underline{\underline{\Lambda}}}}=(1-D){\underline{\underline{\underline{\underline{\Lambda}}}}}^{0}\)
Intégration explicite#
Au niveau des itérations locales (c’est-à-dire en chaque point de GAUSS), l’intégration numérique des équations en vitesse est effectuée par un schéma explicite de Runge-Kutta d’ordre 2 avec découpage automatique en sous-pas locaux en fonction d’une estimation de l’erreur d’intégration (méthode de Runge-Kutta emboîtée) (cf. [R5.03.14]). Le test SSNV225A illustre cette méthode.
Intégration implicite#
Pour l’intégration implicite, on emploiera les notations suivantes:
\({A}^{-}\) , \(A\) et \(\Delta A\) représentent respectivement les valeurs d’une quantité au début et à la fin du pas de temps considéré ainsi que son incrément durant le pas. Le système est discrétisé suivant une thêta-méthode: \(\Delta A=\Delta tg({A}^{-}+\theta \Delta A)=\Delta tg({A}^{\theta})0<\theta \le 1\)
Le problème discrétisé à résoudre est alors le suivant: connaissant l’état au temps \({t}^{-}\) ainsi que les incréments de déformation \(\Delta \varepsilon\) (issus de la phase de prédiction , cf. [R5.03.01]) et de température \(\Delta T\) , déterminer l’état des variables internes au temps \(t\) ainsi que les contraintes \(\sigma\) .
Pour bien prendre en compte la variation des paramètres d’élasticité avec la température, il est nécessaire de discrétiser (de façon implicite) la relation contrainte-déformation élastique de la façon suivante (voir par exemple [R5.03.02]):
\({C}^{\text{-1}}\sigma =\frac{(1-D)}{(1-{D}^{-})}{({C}^{-})}^{\text{-1}}{\sigma}^{-}+(1-D)(\Delta \varepsilon -\Delta {\varepsilon}^{\text{th}})-(1-D)\Delta {\varepsilon}^{p}\) avec \(\Delta {\varepsilon}^{p}=\Delta p\frac{3}{2}\frac{\tilde{\sigma}}{{\sigma}_{\mathit{eq}}}=\Delta pn\)
Pour simplifier les expressions et diminuer le nombre d’opérations lors de la résolution, il est aussi possible d’écrire la première équation en fonction de la déformation élastique; cela suppose de stocker les déformations élastiques ou plastiques comme variables internes. On obtient alors:
\(\Delta {\varepsilon}^{e}-(\Delta \varepsilon -\Delta {\varepsilon}^{\text{th}})+\Delta p{n}^{\theta}=0\) avec \({n}^{\theta}=\frac{3}{2}\frac{\tilde{{\sigma}^{\theta}}}{{\sigma}_{\mathit{eq}}^{\theta}}\) (\({F}_{e}\) )
et les contraintes se calculent ensuite à l’aide des déformations plastiques au temps \({t}^{-}\) par:
\(\sigma =(1-D){\sigma}^{\mathit{nd}}=(1-D)C{\varepsilon}^{e}=(1-D)C({(\varepsilon )}^{-}-{({\varepsilon}^{\mathit{th}})}^{-}-{({\varepsilon}^{p})}^{-}+\theta \Delta {\varepsilon}^{e})\)
Les équations suivantes se déduisent des expressions des dérivées des variables internes:
\(\Delta p-\Delta t{\epsilon}_{0}\sinh\left(\frac{{\sigma}_{\mathit{eq}}^{\theta}(1-{H}_{1}^{\theta}-{H}_{2}^{\theta})}{K\left(1-{D}^{\theta}\right)(1-\varphi )}\right)=0\) (\({F}_{p}\) )
\(\Delta {H}_{i}-\frac{{h}_{i}}{{\sigma}_{\mathit{eq}}}({H}_{i}^{\text{*}}-{\delta}_{i}({H}_{i}^{-}+\theta \Delta {H}_{i}))\Delta p=0\) , pour i=1,2 \(({F}_{{H}_{i}})\)
\(\Delta D-\Delta t{A}_{0}\sinh(\frac{{\alpha}_{D}\text{<}{\sigma}_{p}^{\theta}{\text{>}}_{+}+{\sigma}_{\mathit{eq}}^{\theta}(1-{\alpha}_{D})}{{\sigma}_{0}})=0\) avec \({\sigma}_{p}^{\theta}=\underset{I}{max}{\sigma}_{I}^{\theta}\text{ou}\mathit{tr}({\sigma}^{\theta})\) (\({F}_{D}\) )
On peut formellement écrire ce système: \(F(\Delta Y)=0\) , avec \(\Delta Y={(\Delta {\varepsilon}^{\text{e}},\Delta p,\Delta {H}_{1,}\Delta {H}_{2},\Delta D)}^{t}\)
et \(F(\Delta Y)={({F}_{e},{F}_{p},{F}_{{H}_{1}},,{F}_{{H}_{2}},,{F}_{D})}^{\text{t}}\)
Ce système non-linéaire est résolu par la méthode itérative de Newton [R5.03.14]):
\(F(\Delta {Y}_{k})+{(\frac{\partial F}{\partial \Delta Y})}_{k}(\Delta {Y}_{k+1}-\Delta {Y}_{k})\) en itérant en \(k\) jusqu’à convergence.
La matrice jacobienne du système, nécessaire à la résolution par la méthode de Newton, peut être calculée soit numériquement (ALGO_INTE=”NEWTON_PERT”, cf. test SSNV225B), soit analytiquement.
Dans ce dernier cas l’expression des dérivées est:
\((\frac{\partial {F}_{e}}{\partial \Delta {\varepsilon}^{e}})={I}_{d}+\Delta p\frac{\partial {n}^{\theta}}{\partial \Delta {\varepsilon}^{e}}\) avec \(\frac{\partial {n}^{\theta}}{\partial \Delta {\varepsilon}^{e}}=2\mu \theta \frac{(1-{D}^{\theta})}{{\sigma}_{\mathit{eq}}}[{I}_{\mathit{dev}}-n\otimes n]\) et \({I}_{\mathit{dev}}=\frac{3}{2}({I}_{4}-\frac{1}{3}{I}_{2}\otimes {I}_{2})\)
\((\frac{\partial {F}_{e}}{\partial \Delta p})={n}^{\theta}\) \((\frac{\partial {F}_{e}}{\partial \Delta D})=0\)
\((\frac{\partial {F}_{p}}{\partial \Delta {\varepsilon}^{\text{e}}})=-\Delta t{\varepsilon}_{0}\cosh(\frac{{\sigma}_{\mathit{eq}}^{\theta}(1-{H}_{1}^{\theta}-{H}_{2}^{\theta})}{K(1-{D}^{\theta})(1-\phi )})\frac{(1-{H}_{1}^{\theta}-{H}_{2}^{\theta})}{K(1-{D}^{\theta})(1-\phi )}\frac{\partial {\sigma}_{\mathit{eq}}^{\theta}}{\partial \Delta {\varepsilon}^{\text{e}}}\)
\(\left(\frac{\partial {F}_{p}}{\partial \Delta p}\right)=1\)
\((\frac{\partial {F}_{p}}{\partial \Delta {H}_{i}})=\Delta t{\varepsilon}_{0}\cosh(\frac{{\sigma}_{\mathit{eq}}^{\theta}(1-{H}_{1}^{\theta}-{H}_{2}^{\theta})}{K(1-{D}^{\theta})(1-\phi )})\theta \frac{{\sigma}_{\mathit{eq}}^{\mathit{nd}}}{K(1-\phi )}\)
\((\frac{\partial {F}_{p}}{\partial \Delta D})=0\) car: \(\frac{{\sigma}^{\theta}}{1-{D}^{\theta}}={\sigma}^{\mathit{nd}}\) est indépendant de \(\Delta D\)
\((\frac{\partial {F}_{{H}_{i}}}{\partial \Delta {\varepsilon}^{\text{e}}})=\frac{{h}_{i}}{{\sigma}_{\mathit{eq}}^{2}}\Delta p({H}_{i}^{\text{*}}-{\delta}_{i}{H}_{i}^{\theta})\frac{\partial {\sigma}_{\mathit{eq}}^{\theta}}{\partial \Delta {\varepsilon}^{\text{e}}}\)
\((\frac{\partial {F}_{{H}_{i}}}{\partial \Delta p})=-\frac{{h}_{i}}{{\sigma}_{\mathit{eq}}}({H}_{i}^{\text{*}}-{\delta}_{i}{H}_{i}^{\theta})\)
\((\frac{\partial {F}_{{H}_{i}}}{\partial \Delta {H}_{i}})=1+\frac{{h}_{i}}{{\sigma}_{\mathit{eq}}}{\delta}_{i}\theta \Delta p\)
\((\frac{\partial {F}_{{H}_{i}}}{\partial \Delta D})=\frac{-{h}_{i}}{{\sigma}_{\mathit{eq}}^{2}}\Delta p({H}_{i}^{\text{*}}-{\delta}_{i}{H}_{i}^{\theta})\theta {\sigma}_{\mathit{eq}}^{\mathit{nd}}\)
\((\frac{\partial {F}_{D}}{\partial \Delta {\varepsilon}^{\text{e}}})=-\Delta t\frac{{A}_{0}}{{\sigma}_{0}}\cosh(\frac{{\alpha}_{D}\text{<}{\sigma}_{p}^{\theta}{\text{>}}_{+}+{\sigma}_{\mathit{eq}}^{\theta}(1-{\alpha}_{D})}{{\sigma}_{0}})({\alpha}_{D}\frac{\partial \text{<}{\sigma}_{p}^{\theta}{\text{>}}_{+}}{\partial \Delta {\varepsilon}^{\text{e}}}+(1-{\alpha}_{D})\frac{\partial {\sigma}_{\mathit{eq}}^{\theta}}{\partial \Delta {\varepsilon}^{\text{e}}})\)
\((\frac{\partial {F}_{D}}{\partial \Delta D})=1+\frac{\Delta t{A}_{0}\theta }{{\sigma}_{0}}\cosh(\frac{{\alpha}_{D}\text{<}{\sigma}_{p}^{\theta}{\text{>}}_{+}+{\sigma}_{\mathit{eq}}^{\theta}(1-{\alpha}_{D})}{{\sigma}_{\mathit{nd}}})\left[{\alpha}_{D}<{\sigma}_{p}^{0}>+(1-{\alpha}_{D}){\sigma}_{\mathit{eq}}^{\mathit{nd}}\right]\)
avec : \(\frac{\partial {\sigma}_{\mathit{eq}}^{\theta}}{\partial \Delta {\varepsilon}^{\text{e}}}=2\mu \theta (1-{D}^{\theta}){n}^{\theta}\)
et
\(\frac{\partial \text{<}\mathit{tr}{\sigma}^{\theta}\text{>}}{\partial \Delta {\varepsilon}^{\text{e}}}=\frac{\text{<}\mathit{tr}{\sigma}^{\theta}\text{>}}{\mathit{tr}{\sigma}^{\theta}}(3\lambda +2\mu )\theta (1-{D}^{\theta}){I}_{d}\) dans le cas où \({\sigma}_{p}^{\theta}=\mathit{tr}({\sigma}^{\theta})\)
\(\frac{\partial \text{<}{\sigma}_{1}^{\theta}\text{>}}{\partial \Delta {\varepsilon}^{\text{e}}}=\frac{\text{<}{\sigma}_{1}^{\theta}\text{>}}{{\sigma}_{1}^{\theta}}{I}_{H}\) dans le repère principal, dans le cas où \({\sigma}_{p}^{\theta}={\sigma}_{1}=\underset{I}{max}{\sigma}_{I}^{\theta}\)
avec \({I}_{H}=\left[\begin{array}{cccccc}\lambda +2\mu & \lambda & \lambda & 0& 0& 0\\ \lambda & \lambda +2\mu & \lambda & 0& 0& 0\\ \lambda & \lambda & \lambda +2\mu & 0& 0& 0\\ 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0\end{array}\right]\)
Signification des variables internes#
Les variables internes du modèle aux points de Gauss (mot-clé VARI_ELGA) sontaccessibles par :
\(\mathit{V1}={\varepsilon}_{\text{vp}}^{11}\)
\(\mathit{V2}={\varepsilon}_{\text{vp}}^{22}\)
\(\mathit{V3}={\varepsilon}_{\text{vp}}^{33}\)
\(\mathit{V4}={\varepsilon}_{\text{vp}}^{12}\)
\(\mathit{V5}={\varepsilon}_{\text{vp}}^{13}\)
\(\mathit{V6}={\varepsilon}_{\text{vp}}^{23}\)
\(\mathit{V7}=p\) , la déformation plastique cumulée
\(\mathit{V8}={H}_{1}\) la première variable d’écrouissage isotrope viscoplastique
\(\mathit{V9}={H}_{2}\) , la deuxième variable d’écrouissage isotrope viscoplastique
\(\mathit{V10}=\phi\) , la variable d’endommagement microstructural
\(\mathit{V11}=D\) , la variable d’endommagement
\(\mathit{V12}=0\) , variable interne non utilisée en explicite, indicateur de plasticité en implicite.
Bibliographie#
MUSTATA, R. & HAYHURST, D. «Creep constitutive equations for a 0.5Cr 0.5 Mo 0.25V ferritic steel in the temperature range 565°C-675°C», International Journal of Pressure Vessels and Piping, **2005, 82* , 363-372
PETRY, C. «Comportement en fluage de l’acier P92 : caractérisation expérimentale et modélisation», Note H-T24-2009-00594-FR.
PETRY, C. «Comportement en fluage des jonctions soudées de l’acier P92 : caractérisation expérimentale et modélisation», Note H-T24-2010-00225-FR
LATOURTE, F. « Étude du comportement mécanique de l’acier grade 92 : résultats d’essais de caractérisation de l’interaction fatigue-fluage et premières tentatives de modélisation », Note H-T24-2011-02094-FR
Description des versions du document#
Version Aster |
Auteur(s) Organisme(s) |
Description des modifications |
11.2 |
F. LATOURTE EDF R&D/MMC |
Texte initial |
11.4 |
J.M.PROIX EDF R&D/AMA |
Ajout de l’intégration implicite |