# Modèle de zone cohésive viscoplastique dédié à la surfragmentation à haute température des céramiques nucléaires pressurisées

Clément Tosi, Rodrigue Largenton, Yann Monerie, Pierre-Guy Vincent

► To cite this version:

Clément Tosi, Rodrigue Largenton, Yann Monerie, Pierre-Guy Vincent. Modèle de zone cohésive viscoplastique dédié à la surfragmentation à haute température des céramiques nucléaires pressurisées. 15e colloque national en calcul des structures, Université Polytechnique Hauts-de-France [UPHF], May 2022, Hyères-les-Palmiers, France. hal-03718339

HAL Id: hal-03718339

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

Submitted on 8 Jul 2022

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_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
# Modèle de zone cohésive viscoplastique dédié à la surfragmentation à haute température des céramiques nucléaires pressurisées

C. Tosi<fnref n="1" />,<fnref n="4" />, R. Largenton<fnref n="2" />, Y. Monerie<fnref n="3" />,<fnref n="4" />,
P-G. Vincent<fnref n="1" />,<fnref n="4" />

**Résumé** — Ce travail concerne le comportement d'un fragment de combustible irradié lors d'un hypothétique accident d'insertion de réactivité. Un nouveau modèle de zone cohésive a été développé en suivant l'approche micromécanique de Nkoumbou Kaptchouang et al. [1] et en utilisant le modèle viscoplastique endommageable pour l'UO₂ de Salvo et al. [2].

**Mots clés** — Micromécanique, endommagement, viscoplasticité, zone cohésive, dioxyde d'uranium.

## 1 Contexte

Le combustible en dioxyde d'uranium UO₂ sous forme d'oxyde polycristallin est un combustible utilisé dans les réacteurs à eau sous pression (REP) [3]. L'UO₂ est mis sous forme de pastilles, obtenues par frittage, qui sont ensuite placées dans une gaine cylindrique creuse en alliage de zirconium. L'ensemble gaine et pastille constitue le crayon combustible. Ces crayons sont regroupés pour former un assemblage combustible. Le combustible UO₂ a une microstructure de type fluorine et la taille des grains est d'environ 10 µm [4]. À l'état irradié, des bulles intragranulaires et intergranulaires apparaissent dans le matériau. Leurs tailles et leurs morphologies varient suivant le taux de combustion, le rayon de la pastille et suivant la cote axiale du crayon combustible. Les pores (ou bulles) formés au cours de l'irradiation sont saturés par des gaz de fission (essentiellement du Xénon et Kryton).

Ce travail s'inscrit dans le cadre des recherches menées pour décrire le comportement d'un crayon combustible lors d'un hypothétique accident d'insertion de réactivité (RIA - Reactivity Insertion Accident). Dans le cœur du réacteur, des grappes de commande sont insérées dans les assemblages de crayons combustibles afin de réguler la réaction nucléaire. En cas de rupture du mécanisme de grappe, l'éjection de la grappe résulte de la différence de pression qui existe entre le circuit primaire (155 bars) et l'enceinte de confinement (pression atmosphérique). Cette éjection entraîne un emballement local de la réaction nucléaire pendant quelques dizaines de millisecondes provoquant une augmentation rapide de la température du combustible.

## 2 Développement d'un modèle cohésif

Le modèle de [5] [2] a été retenu pour décrire le comportement viscoplastique endommageable du matériau pour deux raisons d'importance vis à vis de l'application visée en RIA. Premièrement, ce modèle couvre une vaste gamme de taux de déformation (entre 10⁻⁹ et 10⁻<fnref n="1" /> /s), de contraintes (entre 20 et 250 MPa) et de températures (entre 1100 et 1700°C). Deuxièmement, ce modèle introduit de l'endommagement par l'intermédiaire d'une porosité qui peut croître avec le chargement. L'effet d'une pression p des gaz de fission sur le comportement viscoplastique du matériau pourra être pris en compte ultérieurement. Le modèle de [5] [2] est basé sur l'expression d'un potentiel macroscopique en contrainte Ψ,

<footnote n="1">¹ Institut de Radioprotection et de Sûreté Nucléaire (IRSN) PSN-RES/SEMIA/LSMA, (clement.tosi,pierre-guy.vincent)@irsn.fr</footnote>

<footnote n="2">² EDF, rodrigue.largenton@edf.fr</footnote>

<footnote n="3">³ LMGC, Montpellier, yann.monerie@umontpellier.fr</footnote>

<footnote n="4">⁴ Laboratoire MIST, IRSN-CNRS-Université de Montpellier</footnote>
dont dérive le taux de déformation viscoplastique $\dot{\epsilon}$ :

$$ egin{cases} 
\Psi(\sigma_m, \sigma_{eq}, f) = \sigma_0 K \exp\left(\frac{-Q}{RT}\right) \cosh\left(\frac{\Sigma}{\sigma_0}\right) \quad \text{où } \Sigma = \sqrt{\alpha(f)\sigma_m^2 + \beta(f)\sigma_{eq}^2} \\ 
\dot{\epsilon} = \frac{\partial\Psi}{\partial\sigma} = K \exp\left(\frac{-Q}{RT}\right) \sinh\left(\frac{\Sigma}{\sigma_0}\right) \left(\frac{1}{3}\alpha(f)\frac{\sigma_m}{\Sigma}I + \frac{3}{2}\beta(f)\frac{1}{\Sigma}\sigma_d\right) 
\end{cases} 
\quad (1) $$

où $\sigma_m$ est la partie hydrostatique du tenseur des contraintes, $\sigma_{eq}$ la contrainte équivalente de von Mises ($\sigma_{eq} = \sqrt{\frac{3}{2}\sigma_d}$ : $\sigma_d$ où $\sigma_d$ est la partie déviatorique du tenseur des contraintes), f la porosité, $\sigma_0$, K et Q des constantes, R la constante des gaz parfaits, T la température, I le tenseur identité d'ordre 2, $\alpha(f)$ et $\beta(f)$ des fonctions telles que $\alpha(f) = \frac{9}{4}(n(f^{\frac{1}{n}} - 1))^{\frac{-2n}{n+1}}$ et $\beta(f) = (1 + \frac{2}{3}f)(1 - f)^{\frac{-2n}{n+1}}$ et n une constante.

## 2.1 Expression du vecteur contrainte

Un calcul direct permet d'exprimer le potentiel dual de (1) :

$$ egin{cases} 
\Phi(\dot{\epsilon}_m, \dot{\epsilon}_{eq}, f) = \sigma_0 \dot{E} \left( \sinh^{-1} \left( \frac{\dot{E}}{K \exp\left(\frac{-Q}{RT}\right)} \right) - \sqrt{1 + \left( \frac{K \exp\left(\frac{-Q}{RT}\right)}{\dot{E}} \right)^2} \right) \\ 
où \ \dot{E} = \sqrt{\frac{1}{\alpha(f)}(3\dot{\epsilon}_m)^2 + \frac{1}{\beta(f)}(\dot{\epsilon}_{eq})^2} \\ 
\sigma = \frac{\partial \Phi}{\partial \dot{\epsilon}} = \sigma_0 \sinh^{-1} \left( \frac{\dot{E}}{K \exp\left(\frac{-Q}{RT}\right)} \right) \frac{1}{\dot{E}} \left( \frac{3}{\alpha(f)} \dot{\epsilon}_m I + \frac{2}{3\beta(f)} \dot{\epsilon}_d \right) 
\end{cases} 
\quad (2) $$

où $\dot{\epsilon}_m$ représente la composante hydrostatique du tenseur de vitesse de déformation et $\dot{\epsilon}_{eq} = \sqrt{\frac{2}{3}\dot{\epsilon}_d}$ : $\dot{\epsilon}_d$ et $\dot{\epsilon}_d$ est la partie déviatorique du tenseur de vitesse de déformation.

Ce modèle volumique est projeté en un modèle surfacique de type zone cohésive en s'appuyant sur la méthode à cinématique surfacique réduite utilisée dans [1]. La rupture peut être décrite par la localisation des déformations dans des bandes d'épaisseur réduite où les déformations plastiques s'intensifient (voir fig. 1). Le saut de déplacement à travers la bande de localisation des déformations ($S_h$) est noté $[u] =$

<figure><img src="image_2.png" /><figcaption>FIGURE 1 – Modélisation d'une zone de localisation des déformations comme une bande cohésive [1].</figcaption></figure>

$u^+ - u^- = [u_n]\mathbf{n} + [u_t]\mathbf{t} + [u_s]\mathbf{s}$ (avec $(\mathbf{n}, \mathbf{t}, \mathbf{s})$ la base locale normale-tangente représentée sur la fig. 1). Le tenseur de vitesse de déformation au sein de $S_h$ est approché de la manière suivante :

$$ \dot{\epsilon} \approx \begin{bmatrix} \frac{[\dot{u}_n]}{h} & \frac{[u_t]}{2h} & \frac{[\dot{u}_s]}{2h} \\ \frac{[\dot{u}_t]}{2h} & \dot{\gamma} & 0 \\ \frac{[\dot{u}_s]}{2h} & 0 & \dot{\gamma} \end{bmatrix} \quad (3) $$

où $\dot{\gamma}$ est un paramètre permettant d'étendre la plage des états de chargements possibles dans le modèle de zone cohésive.
Le principe de l'approche par zones cohésives est de représenter le comportement de la zone d'amor-
çage de fissure par une relation de traction/séparation, c'est à dire une relation entre le vecteur contrainte
$R = \sigma.n$ et le saut de déplacement $[u]$ à travers $S_h$.

En utilisant (2) et en supposant que $\varepsilon$ est une fonction monotone radiale ($Hyp_1$) telle que $\varepsilon(t) = \dot{\varepsilon}t$,
l'expression suivante du vecteur contrainte est obtenue :

$$
R = \begin{cases}
R_n = \sigma_0 \sinh^{-1} \left( \frac{E}{tK \exp\left(\frac{-Q}{RT}\right)} \right) \frac{1}{E} \left( \frac{1}{\alpha(f)} \left( \frac{[u_n]}{h} + 2\gamma \right) + \frac{4}{9\beta(f)} \left( \frac{[u_n]}{h} - \gamma \right) \right) \\
R_t = \sigma_0 \sinh^{-1} \left( \frac{E}{tK \exp\left(\frac{-Q}{RT}\right)} \right) \frac{1}{E} \left( \frac{1}{3\beta(f)} \frac{[u_t]}{h} \right) \\
R_s = \sigma_0 \sinh^{-1} \left( \frac{E}{tK \exp\left(\frac{-Q}{RT}\right)} \right) \frac{1}{E} \left( \frac{1}{3\beta(f)} \frac{[u_s]}{h} \right) \\
\text{où } E = \sqrt{\frac{1}{\alpha(f)} (3\varepsilon_m)^2 + \frac{1}{\beta(f)} (\varepsilon_{eq})^2}
\end{cases}
$$

La valeur de $\gamma$ est définie de la façon suivante. Le tenseur des contraintes volumique est pris au
voisinage de la zone cohésive dans une approche aux éléments finis cohésifs-volumique et est noté $\sigma^{vol}$.
La triaxialité des contraintes associée s'écrit $\eta_{\sigma^{vol}} = \frac{\sigma_{m}^{vol}}{\sigma_{eq}^{vol}}$. Cette triaxialité est choisie pour déterminer la
triaxialité du taux de déformation à imposer dans la zone cohésive grâce à l'équation suivante : $\eta_{\dot{\epsilon}} =
\frac{\alpha}{3\beta}\eta_{\sigma^{vol}}$. La valeur de $\gamma$ est alors obtenue en résolvant l'équation suivante (d'après (3) et ($Hyp_1$)) :

$$
\eta_{\dot{\epsilon}} = \frac{\dot{\epsilon}_m}{\dot{\epsilon}_{eq}} = \frac{\varepsilon_m}{\varepsilon_{eq}} = \frac{\frac{[u_n]}{h} + 2\gamma}{2\sqrt{\left(\frac{[u_n]}{h} - \gamma\right)^2 + \frac{3}{4}\frac{[u_t]^2}{h^2} + \frac{3}{4}\frac{[u_s]^2}{h^2}}}
$$

2.2 Évolution de la porosité

L'évolution de la porosité (variable d'endommagement) découle du bilan de masse et s'écrit :

$$
\dot{f} = (1 - f) \dot{\epsilon} : \mathbb{I}
$$

L'intégration en temps conduit à :

$$
f = 1 - (1 - f_0) \exp \left( - \frac{[u_n]}{h} - 2\gamma \right)
$$

où $f_0$ est la porosité initiale.

De la même façon que dans les modèles d'endommagement ductile poro-plastiques de type Gurson
[6], une nouvelle variable de porosité effective $f^*$ est introduite, permettant d'accélérer l'évolution de
la porosité à partir d'une porosité critique notée $f_c$, relatif à la coalescence des pores dans le matériau.
Au delà de $f_c$, la triaxialité en contrainte $n_{\sigma}^{vol}$ est fixée et n'évolue plus (valeur notée $n_{\sigma}^{coalescence}$). Cette
porosité effective s'écrit :

$$
f^* = \begin{cases} f & \text{si } f \le f_c \\ f_c + k(f - f_c) & \text{si } f \ge f_c \end{cases}
$$

avec $k$ une constante.

2.3 Transformation en modèle intrinsèque

Le modèle de zone cohésive obtenu est transformé en un modèle intrinsèque via une fonction tendant
rapidement vers 1 lorsque $[u]$ croit :

$$
R^{\text{int}} = \left( 1 - \exp\left(-\alpha_i \frac{|[u_i]|}{h}\right) \right) R_i \quad \text{où } i \in \{n, t, s\}
$$

où $\alpha_n$, $\alpha_t$ et $\alpha_s$ sont des paramètres.
## 2.4 Traitement de la décharge

Afin de décrire le côté irréversible du processus de fissuration, il convient de définir un critère de décharge, ainsi qu’un comportement en décharge.

La décharge survient respectivement dans les directions **n**, **t**, **s** lorsque :

$$|[u_i]^{(t)}| - |[u_i]^{(t-\delta t)}| < 0 \text{ où } i \in \{n,t,s\}$$

où $t$ est le temps et $\delta t$ le pas de temps.

Le comportement est choisi de façon à ce que la décharge s'effectue de manière élastique vers l'origine. Le vecteur contrainte lors de la décharge est écrit de la manière suivante :

$$\bar{R}_i = (R_i^{int})_{\max} \frac{[u_i]}{[u_i]_{\max}} \text{ où } i \in \{n,t,s\} \text{ et } [u_i]_{\max} = \max(|[u_i]^{(\tau)}|, 0 \le \tau \le t)$$

$(R_i^{int})_{\max}$ est la composante $i$ du vecteur contrainte correspondant aux sauts de déplacement ($[u_n]_{\max}$, $[u_t]_{\max}$, $[u_s]_{\max}$).

## 2.5 Résultats

Le modèle décrit par les équations précédentes a été implémenté sous Python.

La réponse du modèle cohésif est illustrée sur les fig. 2, 3 et 4 en faisant varier respectivement la vitesse de déformation $\dot{\epsilon}_{nn}$ ($\dot{\epsilon}_{nt} = \dot{\epsilon}_{ns} = 0$), la température $T$ et la triaxialité en contrainte prise au voisinage de la bande cohésive $\eta_{\sigma^{vol}}$. Des coefficients matériaux issus des thèses de [5] et [7] pour l'UO$_2$ sont utilisés (voir tableau 1).

Dans tous les cas, la contrainte cohésive atteint un seuil d'écoulement, puis décroît lentement, jusqu'à ce que la porosité atteigne la porosité critique de coalescence. La valeur du coefficient $k$ étant élevée, la porosité croît alors rapidement jusqu'à 1 et la contrainte cohésive chute vers 0. De plus, l'augmentation de la vitesse de déformation $\dot{\epsilon}_{nn}$ a pour effet d'augmenter le seuil d'écoulement. Ce phénomène est aussi observé lorsque la température $T$ décroît. L'augmentation de la triaxialité en contrainte $\eta_{\sigma^{vol}}$ agit aussi en augmentant le seuil d'écoulement mais en rendant le matériau plus fragile (la porosité critique $f_c$ est atteinte bien plus rapidement).

<table><caption>TABLE 1 – Coefficients utilisés pour illustrer la réponse du modèle cohésif.</caption><thead><tr><th>K</th><th>Q</th><th>σ<sub>0</sub></th><th>R</th><th>n</th><th>f<sub>0</sub></th><th>f<sub>c</sub></th><th>k</th></tr></thead><tbody><tr><td>29130 s<sup>-1</sup></td><td>482 kJ/mol</td><td>5 MPa</td><td>8.314 J/mol/K</td><td>6</td><td>0.04</td><td>0.041130348</td><td>1000</td></tr></tbody></table>
<figure><img src="image_3.png" /><figcaption>FIGURE 2 – Évolution de la contrainte cohésive $R_n$ (haut) et de la porosité $f$ (bas) en fonction de la déformation $\varepsilon_{nn} : \eta_{\sigma^{vol}} = \frac{1}{3}, T = 1973K, \varepsilon_{nn} = 0.01s^{-1}$ (bleu), $\varepsilon_{nn} = 0.1s^{-1}$ (orange), $\varepsilon_{nn} = 0.5s^{-1}$ (vert).</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>FIGURE 3 – Évolution de la contrainte cohésive $R_n$ (haut) et de la porosité $f$ (bas) en fonction de la déformation $\varepsilon_{nn} : \eta_{\sigma^{vol}} = \frac{1}{3}, \varepsilon_{nn} = 0.01s^{-1}, T = 1373K$ (bleu), $T = 1773K$ (orange), $T = 1973K$ (vert).</figcaption></figure>

### 3 Conclusion

Ce modèle cohésif offre la possibilité, via un endommagement surfacique et une loi de traction-séparation, de résumer la dégradation du combustible jusqu’à création d’une fissure en tout point de la pastille. Pour poursuivre ce travail, il est prévu d’ajouter au modèle existant un modèle de rupture basé sur un modèle plastique de type Drucker-Prager poreux afin de rendre compte de l’asymétrie du comportement en compression et en traction. Cet ajout s’appuiera sur les travaux de [8] et [9] autour de l’écriture d’un modèle micromécanique de type Drucker-Prager poreux et sur les travaux de [10] autour de l’utilisation d’un modèle de Drucker-Prager pour décrire la rupture des joints de grain dans l’UO₂.
<figure><img src="image_5.png" /></figure>

FIGURE 4 – Évolution de la contrainte cohésive $R_n$ (haut) et de la porosité $f$ (bas) en fonction de la déformation $\epsilon_{nn}: T = 1973\text{K}, \epsilon_{nn} = 0.01\text{s}^{-1}, \eta_{\sigma vol} = \frac{1}{3}$ (bleu), $\eta_{\sigma vol} = 1$ (orange), $\eta_{\sigma vol} = 2$ (vert).

Références

[1] N.B. Nkoumbou Kaptchouang, Y. Monerie, F. Perales, and P.G. Vincent. Cohesive GTN model for ductile fracture simulation. *Engineering Fracture Mechanics*, 242:107437, 2021.

[2] M. Salvo, J. Sercombe, J-C. Ménard, J. Julien, T. Helfer, and T. Désoyer. Experimental characterization and modelling of UO2 behavior at high temperatures and high strain rates. *Journal of Nuclear Materials*, 456:54–67, 2015.

[3] J.F. Parisot, B. Bernard, J.L. Guillet, and Y. Guérin. Les combustibles nucléaires. *Editions du Moniteur*, 242:107437, 2008.

[4] F. Dherbey, F. Louchet, A. Mocellin, and S. Leclercq. Elevated temperature creep of polycrystalline uranium dioxide : From microscopic mechanisms to macroscopic behaviour. *Acta Materialia*, 50(6):1495–1505, 2002.

[5] M. Salvo. *Etude expérimentale et modélisation du comportement mécanique du combustible UO₂ en compression à haute température et forte vitesse de sollicitation*. PhD thesis, Université d'Aix Marseille, 2014.

[6] A. L. Gurson. Continuum theory of ductile rupture by void nucleation and growth : Part I - yield criteria. *Journal of Engineering Materials and Technology*, 99:2–15, 1977.

[7] N.B. Nkoumbou Kaptchouang. *Modélisation micromécanique de l'endommagement ductile par une approche cohésive - volumique : application á l'UO2 irradié*. PhD thesis, Université de Montpellier, 2019.

[8] T.F. Guo, J. Faleskog, and C.F. Shih. Continuum modeling of a porous solid with pressure-sensitive dilatant matrix. *Journal of the Mechanics and Physics of Solids*, 56(6):2188–2212, 2008.

[9] W.Q. Shen, J. Zhang, J.F. Shao, and D. Kondo. Approximate macroscopic yield criteria for drucker-prager type solids with spheroidal voids. *International Journal of Plasticity*, 99:221–247, 2017.

[10] M. Salvo, J. Sercombe, T. Helfer, P. Sornay, and T. Désoyer. Experimental characterization and modeling of UO2 grain boundary cracking at high temperatures and high strain rates. *Journal of Nuclear Materials*, 460:184–199, 2015.