# Finite element modelling of unstable crack propagation of nuclear steel at high temperature with cohesive elements

Nicolas Tardif, Alain Combescure, Michel Coret, Philippe Matheron, Vincent Koundy, Grégory Nicaise

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

## To cite this version:

Nicolas Tardif, Alain Combescure, Michel Coret, Philippe Matheron, Vincent Koundy, et al.. Finite element modelling of unstable crack propagation of nuclear steel at high temperature with cohesive elements. 9e Colloque national en calcul des structures, CSMA, May 2009, Giens, France. hal-01412110

HAL Id: hal-01412110

https://hal.science/hal-01412110v1

Submitted on 8 Dec 2016

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

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

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
# Simulation éléments finis de la fissuration d'un acier à haute température à l'aide d'éléments cohésifs

N. Tardif¹², A. Combescure¹, M. Coret¹, P. Matheron³, V. Koundy², G. Nicaise²

¹ Université de Lyon, CNRS, INSA / LAMCOS, UMR/CNRS 5514
Bât. J. d'Alembert, 20 avenue Albert Einstein, 69621 Villeurbanne Cedex
{nicolas.tardif, alain.combescure, michel.coret} @insa-lyon.fr

² IRSN (DSR/SAGR/BPhAG)
BP17, 92262 Fontenay-aux-Roses Cedex
{vincent.koundy, gregory.nicaise} @irsn.fr

³ CEA, DEN, DM2S, SEMT, LISN
F-91191 Gif-sur-Yvette, France.
philippe.matheron@cea.fr

**Résumé** — Une formulation 3D d'élément cohésif en grands déplacements et grandes déformations est proposée pour simuler la fissuration à haute température (900°C) d'une éprouvette CT dans l'acier de cuve français 16MND5. Nous lui associons une loi de traction séparation dont la forme permet à l'élément de dissiper une grande quantité d'énergie pour une faible valeur de CTOD.

**Mots clés** — modèle de zone cohésive, formulation 3D en grand déplacement et grande déformation, rupture à 900°C.

## 1 Introduction

En cas d'accident grave sur un réacteur à eau pressurisée (REP) avec fusion du coeur combustible et des éléments de structures internes, le fond de cuve est soumis à un chargement thermique et mécanique significatif. Pour un scénario hypothétique où la cuve fissurerait par fluage, il est nécessaire de prédire le temps avant l'initiation de la fissure, sa taille ainsi que sa localisation afin de déterminer une stratégie de gestion d'accident hors cuve, dans l'enceinte de confinement.

Nous proposons d'étudier plus particulièrement l'initiation et la propagation quasistatique (ou stable) de fissure en mode I de la nuance *Krakatoa* de l'acier de cuve français 16MND5 à 900°C. Cette étude a fait l'objet d'essais de fissuration sur éprouvette CT entaillée, réalisés au CEA Saclay. Une simulation numérique de ces essais est proposée à l'aide d'un modèle éléments finis 3D dans le code de calcul abaqus. L'initiation et la propagation de la fissure sont modélisées à l'aide d'éléments cohésifs surfaciques en grands déplacements et grandes déformations.

L'application du modèle cohésif à la rupture ductile des métaux a été initiée principalement par Tvergaard et Hutchinson [6] et étendue par Siegmund et Brocks [4]. Les lois de traction-séparation qu'ils proposent sont adaptées au phénomène de nucléation, croissance et coalescence de vides propre à la rupture ductile.
Peu d'études ont été réalisées à l'aide d'éléments cohésifs dans une formulation 3D en grandes déformations [1, 3, 2]. Ce type de formulation est cependant indispensable à la bonne modélisation de notre cas d'étude du fait de la forte ductilité de notre matériau à haute température.

Une loi de traction séparation du type de Tvergaard et Hutchinson modifiée par l'ajout d'une pente est utilisée. La forme de la loi de comportement permet à l'élément de dissiper une grande quantité d'énergie pour une faible valeur de CTOD.

## 2 Formulation des éléments cohésifs surfaciques en grands déplacements et grandes déformations

Initialement sans épaisseur, l'élément isoparamétrique et bilinéaire d'interface comporte 8 noeuds et 4 points d'intégration. Les noeuds sont répartis en deux surfaces de 4 noeuds qui représentent une partie des lèvres inférieures et supérieures de la fissure en cours de décohésion. Le mouvement des parties supérieures et inférieures du volume induit un déplacement des noeuds de l'élément d'interface. Seuls les déplacements relatifs qui représentent les trois modes de ruptures sont pris en compte dans la formulation (cf figure 1).

<figure><img src="image_3.png" /><figcaption>Figure 1 – élément cohésif dans sa configuration initiale (a), sa configuration déformée (b), mode de séparation (c)</figcaption></figure>

Les coordonnées des noeuds de l'élément dans le repère global (U, V, W) peuvent être écrit à l'aide d'un vecteur U de taille 24x1 où l'indice représente le noeud considéré :

$$U_{24x1} = (U_1, V_1, W_1, \dots, U_8, V_8, W_8)^T \qquad (1)$$

Le calcul de la séparation en tout point de l'élément est associé à un repère de référence attaché à son plan médian [2]. Ce repère de référence (u, v, w) est déterminé par rotation $R_{3x3}^T$ du repère global (figure 1). Les coordonnées u et v correspondent à deux vecteurs normés appartenant au plan médian, le vecteur w est le vecteur normé orthogonal au plan médian :

$$(u, v, w)^T = R_{3x3}^T (U, V, W)^T \qquad (2)$$

Nous pouvons alors déterminer la valeur de la séparation relative des surfaces $\delta_{12x1}$ en chacun des 4 points de Gauss G de coordonnées $(\xi_g, \eta_g)$ dans le repère paramétrique de l'élément :

$$\delta_{12x1}(\xi_g, \eta_g) = \underbrace{\begin{pmatrix} -\mathcal{N}_1(\xi_g, \eta_g)\mathbf{R}_{3x3}^T, & -\mathcal{N}_2(\xi_g, \eta_g)\mathbf{R}_{3x3}^T, & -\mathcal{N}_3(\xi_g, \eta_g)\mathbf{R}_{3x3}^T, & -\mathcal{N}_4(\xi_g, \eta_g)\mathbf{R}_{3x3}^T \\ \mathcal{N}_1(\xi_g, \eta_g)\mathbf{R}_{3x3}^T, & \mathcal{N}_2(\xi_g, \eta_g)\mathbf{R}_{3x3}^T, & \mathcal{N}_3(\xi_g, \eta_g)\mathbf{R}_{3x3}^T, & \mathcal{N}_4(\xi_g, \eta_g)\mathbf{R}_{3x3}^T \end{pmatrix}}_{\mathbf{B}_{12x24}(\xi_g, \eta_g)} U_{24x1} \qquad (3)$$
Les fonctions $\mathcal{N}_i$ correspondent aux 4 fonctions d'interpolations bilinéaires polynomiales usuelles exprimées dans le repère paramétrique de l'élément.

Le vecteur traction $\mathbf{t}_{12x1}$ est calculé en chaque point de Gauss à partir de la loi de traction séparation C présentée dans la section suivante :

$$ \mathbf{t}_{12x1} = C \cdot \delta_{12x1} \qquad (4) $$

Le calcul des forces internes est écrit sur la configuration de fin de pas $n+1$. J correspond au jacobien permettant de passer des coordonnées paramétriques propres à l'élément aux coordonnées cartésiennes.

$$ \mathbf{F} = \int_{S_i \text{ coh}} \mathbf{B}_{n+1}^T \mathbf{t}_{n+1} ds = \int_{\square} \mathbf{B}_{n+1}^T \mathbf{t}_{n+1} J_{n+1} d\square \qquad (5) $$

En utilisant une méthode de Gauss pour le calcul numérique de l'intégration, on obtient les forces nodales :

$$ \mathbf{F}_{24x1} = S_{n+1} \sum_g \omega_g \mathbf{B}_{n+1}^T (\xi_g, \eta_g) \mathbf{t}_{n+1} (\xi_g, \eta_g) \qquad (6) $$

Le calcul de la matrice de raideur s'effectue dans les mêmes conditions. Une matrice de raideur tangente est utilisée.

$$ \mathbf{K}_{\mathbf{T}} = \int_{S_i \text{ coh}} \mathbf{B}_{n+1}^T \mathbf{DB}_{n+1} ds = \int_{\square} \mathbf{B}_{n+1}^T \mathbf{DB}_{J_{n+1}} d\square \qquad (7) $$

$$ \mathbf{K}_{\mathbf{T}24x24} = S_{n+1} \sum_g \omega_g \mathbf{B}_{n+1}^T (\xi_g, \eta_g) \frac{d\mathbf{t}}{d\delta_{n+1}} (\xi_g, \eta_g) \mathbf{B}_{n+1} (\xi_g, \eta_g) \qquad (8) $$

Le calcul de $\mathbf{K}_{\mathbf{T}}$ ne garantit pas le principe d'objectivité car il ne prend pas en compte les effets de la rotation de la surface cohésive considérée. Cette simplification est possible car notre modèle ne comprend pas de rotation de cette surface lors du chargement [3].

## 2.1 Loi de traction séparation

La loi de traction séparation représente, à une échelle macroscopique, les phénomènes dissipatifs conduisant à la décohésion totale des lèvres de la fissure dans la process zone en pointe de fissure. Dans le cas de la rupture d'un métal, la loi de traction séparation représente donc les phénomènes de nucléation, croissance et coalescence de vide.

La forme de la loi de traction-séparation n'est a priori pas connue. Tvergaard et Hutchinson [6] ont proposé la loi décrite dans la figure 2(a). En effet, cette loi permet de restreindre la contrainte à rupture à quelques multiples de la limite d'élasticité tout en restituant un taux d'énergie surfacique important pour un CTOD restant assez faible.

Ils ont montré que, dans le cas d'un matériau élasto-plastique, dont la zone plastique reste confinée en pointe de fissure, la forme de la loi n'a que peu d'influence sur la propagation de la fissure. Seuls les paramètres de traction maximale $t_{max}$ et le taux d'énergie dissipée lors de la décohésion des deux surfaces $\Gamma_0$ sont très influents.

Un bon candidat pour montrer cette influence a été proposé par Turner et Kolednik [5]. Ils ont caractérisé la résistance à la rupture R d'un matériau ductile comme étant l'énergie dissipée nécessaire à l'extension de la surface fissurée d'un incrément $Bda$. Cette énergie est séparée en deux contributions : une énergie dissipée en volume par plasticité et une énergie surfacique dissipée lors de la création des deux lèvres de la fissure. Cette dernière contribution correspond au paramètre $\Gamma_0$ de la loi de traction séparation.
$$R = \frac{1}{B} \frac{d(U - E_e)}{da} = \frac{1}{B} \frac{dE_{pl}}{da} + \Gamma_0 \quad (9)$$

Considérons une loi de comportement en volume élastoplastique avec un écrouissage quelconque. A $\Gamma_0$ fixée, plus la traction maximale $t_{max}$ est importante devant la limite d'élasticité, plus le paramètre R augmente car le taux d'énergie dissipée plastiquement dans le volume augmente. A $t_{max}$ fixée, plus l'énergie surfacique $\Gamma_0$ augmente plus R augmente, de façon intrinsèque, mais aussi parce que la taille de la zone du processus de rupture s'agrandit et provoque l'augmentation du taux d'énergie dissipée plastiquement.

Dans notre étude, les effets combinés de la viscoplasticité et de la plasticité étendue nous obligent à rajouter une pente à la loi de traction séparation afin de modéliser proprement la propagation (figure 2(b)).

<figure><img src="image_4.png" /><figcaption>Figure 2 – Lois de traction séparation</figcaption></figure>

# 3 Simulation éléments finis de la fissuration d'une éprouvette CT à 900°C

## 3.1 Maillage et conditions aux limites

L'éprouvette CT a une géométrie hors normes basée sur une CT20. Le maillage de l'éprouvette CT et les conditions aux limites sont représentées sur la figure 3.

<figure><img src="image_5.png" /><figcaption>Figure 3 – Maillage et conditions aux limites</figcaption></figure>
Les éléments constitutifs du volume sont des éléments brique à intégration réduite (8 noeuds et 1 point d'intégration). Les raccords entre taille de maille et les raccords au niveau des congés de l'éprouvette sont éffectués à l'aide d'éléments prisme (6 noeuds et 2 points d'intégrations).

Les éléments de la zone cohésive sont disposés de telle sorte qu'ils suivent le profil de la pré-
entaille. La surface est quadrillée par 43 éléments dans la direction de la propagation et 20 dans la
direction du front.

Les goupilles sont modélisées par des éléments coques rigides. La traction est appliquée sur l'axe des goupilles à une vitesse de 2.5mm.min⁻¹.

3.2 Comportement du volume et de l'interface

A la température étudiée (900°C), le comportement du matériau dépend de la vitesse de déformation. Cette dépendance est prise en compte par la donnée, à différentes vitesses de déformation logarithmique, d'un nuage de courbes élasto-plastiques à écrouissage isotrope (figure 4(a) E = 23466MPa, ν = 0.3). Le code de calcul réalise une interpolation linéaire pour définir le comportement aux vitesses de déformations non connues.

La loi de traction-séparation en mode 1 est représentée sur la figure 4(b). L'interface est considérée comme rompue pour une décohésion totale c'est à dire pour une séparation supérieure à $\delta_c$ (cf figure 2(b)).

<figure><img src="image_6.png" /><figcaption>Figure 4 – Comportement</figcaption></figure>

3.3 Procédure de calcul

Le calcul est effectué en quasi-statique dans le sens où les termes d'inertie et d'accélération sont négligés lors du calcul de l'équilibre mécanique. Le chargement est cependant corrélaté à la base de temps des essais afin de prendre en compte les phénomènes visqueux du matériau.

La sous intégration des éléments du volume provoquent des modes de déplacements sablier au niveau des éléments adjacents à la zone cohésive. Ces modes de déplacement influencent la propagation de la fissure. Ils sont donc régularisés par l'ajout d'une raideur s'opposant à leur développement.

3.4 Résultats

La simulation offre une très bonne adéquation avec les résultats expérimentaux (figure 5(a)) en
terme de propagation de fissure et de courbe effort déplacement. Il est à noter que la détermination
expérimentale de la propagation de fissure résulte d'une interpolation linéaire d'un signal bruité
(réalisée à partir d'une avancée de fissure de 2mm). Cette remarque peut expliquer la différence
entre les courbes en début de propagation.
<figure><img src="image_7.png" /><figcaption>Figure 5 – Résultats numériques</figcaption></figure>

La figure 5(b) montre la répartition, lors de l'avancée de la fissure, de l'énergie dissipée in-
elastiquement dans le volume (89% de l'énergie interne à la fin du calcul) et de l'énergie dissipée
pour la création des lèvres de la fissure dans la zone cohésive (9.3% de l'énergie interne à la fin du
calcul). Par ailleurs, on peut remarquer que l'énergie utilisée pour le contrôle des modes sabliers
est négligeable devant l'énergie dissipée dans la zone cohésive ce qui conforte l'utilisation de ce
mode de contrôle pour notre simulation.

4 Conclusion

Ce papier présente une formulation d'élément cohésif surfacique en grands déplacements et en grandes déformations. La loi de traction séparation, qui lui est associée, permet de simuler la propagation de fissure dans un matériau hautement ductile et viscoplastique. En effet, sa forme induit une dissipation importante d'énergie surfacique pour une faible valeur de CTOD. Par ailleurs, les effets combinés de la plasticité étendue et de la viscoplasticité du matériau engendrent une dépendance plus marquée à la forme de la loi de traction séparation.

Références

[1] C. R. Chen, O. Kolednik, J. Heerens, and F. D. Fischer. Three dimensional modeling of ductile crack growth : Cohesive zone parameters and crack tip triaxiality. *Engng Fract Mech*, 72 :2072-2094, 2005.

[2] M. Ortiz and A. Pandolfi. Finite deformation irreversible cohesive elements for three dimensional crack propagation analysis. *Int J Numer Meth Engng*, 44:1267-1282, 1999.

[3] S. Roychowdhury, Y. D. A. Roy, and R. H. Dodds Jr. Ductile tearing in thin aluminium panels : experiments and analyses using large displacement, 3 D surface cohesive elements. *Engng Fract Mech*, 69:983-1002, 2002.

[4] T. Siegmund and W. Brocks. A numerical study on the correlation between the work of separation and the dissipation rate in ductile fracture. *Engng Fract Mech*, 67:139–154, 2000.

[5] C. E. Turner and O. Kolednik. A micro and macro approach to the energy dissipation rate model of stable ductile crack growth. *Fatigue Fract Engng Mat Struct*, 17:1089–1107, 1994.

[6] V Tvergaard and J. W. Hutchinson. The relation between crack growth resistance and fracture process parameters in elastic-plastic solids. *J Mech Phys Solids*, 6:1377–1397, 1992.