# Un schéma numérique efficace pour le traitement de la fissuration fragile des matériaux quasi-fragiles par champ de phase

Thomas Helfer, Benoît Bary, Olivier Fandeur, Lu Ye

► To cite this version:

Thomas Helfer, Benoît Bary, Olivier Fandeur, Lu Ye. Un schéma numérique efficace pour le traitement de la fissuration fragile des matériaux quasi-fragiles par champ de phase. 24e Congrès Français de Mécanique, 2019, Brest, France. hal-03106063

HAL Id: hal-03106063

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

Submitted on 20 Dec 2024

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 Attribution 4.0 International License
# Un schéma numérique efficace pour le traitement de la fissuration fragile des matériaux quasi-fragiles par champ de phase

Y. Luª, T. Helferª, B. Baryᵇ, O. Fandeurᶜ

a. CEA, DEN, DEC, SESC, LSC, Centre de Cadarache, France, {ye.lu,thomas.helfer}@cea.fr

b. CEA, DEN, DANS, DPC, SECR, LECBA, Centre de Saclay, France, benoit.bary@cea.fr

c. CEA, DEN, DANS, DM2S, SEMT, LM2S, Centre de Saclay, France, olivier.fandeur@cea.fr

## Résumé :

L'implémentation numérique de la fissuration fragile par champ de phase se fait soit par une approche monolithique, soit par une approche échelonnée. La première pose des problèmes de robustesse et ne peut être implantée dans tous les codes traditionnels (Cast3M ici) sans développement supplémentaire. La seconde est connue pour converger très lentement. Ce travail se situe dans la continuité du papier [1] où une approche échelonnée avait été proposée. Différentes conditions d'irréversibilité et algorithmes de résolution sont comparés. Enfin, un nouveau schéma de type semi-implicite est proposé pour un traitement efficace de la propagation instable.

## Mots clefs : Conditions d'irréversibilité, Schéma numérique, Semi-implicite, Propagation instable, Champ de phase

## 1 Introduction

La simulation des éléments combustibles est une mission de la Direction de l'Énergie Nucléaire du CEA [2]. Cette simulation fait intervenir de nombreux phénomènes physiques fortement couplés : thermique et mécanique non linéaire, neutronique, thermo-hydraulique, évolutions physico-chimiques (dénсification et gonflement sous irradiation, relâchement de produits de fission gazeux, etc.), changement de phases, etc.

Un réseau complexe de fissures apparaît dans les pastilles combustibles dès la première montée en puissance du fait d'un fort gradient de température, de l'ordre de 150 K.mm⁻¹. Étant données les durées caractéristiques des irradiations des combustibles nucléaires, de l'ordre de l'année, les simulations sont effectuées en quasi-statique : il s'agit donc plus de trouver des états d'équilibre après propagation de fissures que de décrire finement cette propagation.

La modélisation de cette fissuration fragile est actuellement prise en compte par un modèle d'endommagement local régularisé par la taille de maille [3, 4]. Afin de s'affranchir des limites du modèle actuel (dépendance potentielle de maillage, absence de principe énergétique, etc.), une approche par champ de phase [5, 6] est envisagée.
Une première implémentation d'une approche par champ de phase dans Cast3M [1], a montré que l'on pouvait retrouver qualitativement un réseau de fissures proche de ce qui est observé expérimentalement. Cependant, l'implémentation proposée manquait de robustesse et de performance. L'objet de ce papier est d'en montrer une évolution récente permettant d'en améliorer les performances et la robustesse.

# 2 Rappels sur le traitement de la fissuration fragile par champ de phase

## 2.1 Régularisation spatiale du problème

La méthode de champ de phase consiste à approximer des fissures discrètes par une fonction régulière dont la valeur est égale à 1 sur les points de la fissure et tend vers 0 lorsque l'on s'en éloigne (on adopte ici la même définition que dans [5]). Supposons que la fissure apparaît au point $x = 0$, un exemple 1D de ce type de fonction s'écrit alors :

$$d(x) = e^{-|x|/l} \qquad (1)$$

où $d$ est la variable d'endommagement, $l$ définit la taille de zone endommagée et apparaît comme une longueur caractéristique équivalente à la taille d'un volume élémentaire représentatif de la microstructure. Il peut être démontré que l'équation (1) est la solution du problème variationnel :

$$d = \operatorname{Arg} \left\{ \inf_{d \in W} I(d) \right\} = \operatorname{Arg} \left\{ \inf_{d \in W} \int_{\Omega} \frac{1}{2} (d^2 + l^2 d'^2) dV \right\} \qquad (2)$$

avec $W = \{d|d(0) = 1, d(\pm\infty) = 0\}$. On peut définir la densité de surface endommagée $\gamma$ dans le domaine $\Omega$ :

$$\gamma(d, d') = \frac{1}{2l} (d^2 + l^2 d'^2) \qquad (3)$$

Dans le cas multi-dimensionnel, elle s'exprime comme :

$$\gamma(d, \nabla d) = \frac{1}{2l} (d^2 + l^2 |\nabla d|^2) \qquad (4)$$

## 2.2 Formulation du problème

Considérons un corps élastique linéaire $\Omega$, avec une fissure $\Gamma$ et une force extérieure $\bar{t}$ appliquée à la surface $\partial\Omega$, la fonctionnelle potentielle peut s'écrire comme la somme de l'énergie de déformation, du potentiel de dissipation de rupture et du travail des forces extérieures :

$$\Pi(\mathbf{u}, \Gamma) = \Pi_e + \Pi_d - \Pi_{\text{ext}} = \int_{\Omega} \psi_e(\epsilon(\mathbf{u}))dV + \int_{\Gamma} G_c d\Gamma - \int_{\partial\Omega} \bar{t} u dA \qquad (5)$$

où $\psi_e$ est la densité de l'énergie élastique, $\mathbf{u}$ est le déplacement, $G_c$ désigne l'énergie nécessaire pour créer une fissure de longueur d'unité. Avec l'hypothèse des petites perturbations, nous pouvons définir la déformation comme suit :

$$\epsilon = \nabla_s \mathbf{u} \qquad (6)$$

où $\nabla_s$ désigne l'opérateur de gradient symétrique.

Dans la méthode de champ de phase, la fissure discrète est approximée par une zone endommagée de
densité $\gamma$. L'énergie de dissipation peut donc être approximée par :

$$ \Pi_d = \int_{\Gamma} G_c d\Gamma \approx \int_{\Gamma} \frac{G_c}{2l} (d^2 + l^2 |\nabla d|^2) dV \quad (7) $$

Le potentiel $\Pi(\mathbf{u}, \Gamma)$ est donc remplacé par $\Pi(\mathbf{u}, d)$ :

$$ \Pi(\mathbf{u}, d) = \Pi_e + \Pi_d - \Pi_{\text{ext}} \quad (8) $$

## 2.3 Choix de l'énergie libre élastique

Pour éviter l'ouverture de fissure en compression, l'énergie élastique est séparée en deux parties à l'aide d'une décomposition spectrale de la déformation [5] : $\epsilon = \epsilon^+ + \epsilon^-$. Désignant respectivement par $\psi_e^+$ et $\psi_e^-$ les énergies des modes de traction et compression, l'énergie de déformation devient :

$$ \Pi_e = \int_{\Omega} (g(d)\psi_e^+(\epsilon) + \psi_e^-(\epsilon)) dV \quad (9) $$

et la contrainte s'écrit :

$$ \boldsymbol{\sigma} = g(d)(\lambda\langle\text{tr}[\epsilon]\rangle^+\mathbf{I} + 2\mu\epsilon^+) + \lambda\langle\text{tr}[\epsilon]\rangle^-\mathbf{I} + 2\mu\epsilon^- \quad (10) $$

où une fonction de dégradation : $g(d) = (1-d)^2 + k$ apparaît pour prendre en compte l'effet d'endommagement sur la rigidité du matériau. $k$ est une petite valeur introduite pour éviter le problème de singularité. $\mathbf{I}$ est la matrice d'identité. $\langle \cdot \rangle^{\pm} = (\cdot \pm |\cdot|)$.

## 2.4 Équilibre mécanique

L'équation d'équilibre classique s'obtient par minimisation de l'énergie $\Pi(\mathbf{u}, d)$ par rapport au déplacement. Cela conduit aux équations d'équilibre classiques de la mécanique :

$$ \delta_u \Pi = 0, \forall \delta \mathbf{u} \quad \rightarrow \quad \left\{ \begin{array}{ll} \text{div } \boldsymbol{\sigma} = 0 & \text{dans } \Omega \\ \boldsymbol{\sigma} \cdot \mathbf{n} = \bar{\mathbf{t}} & \text{sur } \partial_{\sigma} \Omega \\ \mathbf{u} = \bar{\mathbf{u}} & \text{sur } \partial_u \Omega \end{array} \right. \quad (11) $$

À endommagement constant, l'équilibre mécanique est classiquement résolu soit par une méthode de Newton, soit par une méthode de point fixe accélérée, ces deux méthodes étant à notre disposition. Dans les exemples présentés, la seconde option a été utilisée car elle apparaît plus efficace, malgré la convergence quadratique de la méthode de Newton.

## 2.5 Traitement de l'irréversibilité

Dans le travail original de Francfort et Marigo [6, 7], l'irréversibilité traduisait le fait que les fissures ne pouvaient se guérir. Les auteurs justifiaient cette approche de l'irréversibilité par un principe de stabilité. La manière de transposer cette irréversibilité sur l'évolution de l'endommagement $d$ diffère suivant les auteurs.
Le point de départ est toujours la minimisation du potentiel $\Pi(\mathbf{u}, d)$ par rapport à $d$, ce qui conduit aux équations suivantes :

$$ \delta_d \Pi = 0, \forall \delta_d \rightarrow \begin{cases} \frac{G_c}{l} (d + l^2 \Delta d) = 2(1-d) \psi_e^+ & \text{dans } \Omega \\ \nabla d \cdot \mathbf{n} = 0 & \text{sur } \partial \Omega \end{cases} \quad (12) $$

À déplacement constant, cette équation est généralement résolue par analogie avec un problème de thermique pour pouvoir être utilisée dans la plupart des codes éléments finis. Elle est donc en principe extrêmement peu coûteuse à résoudre.

Cependant, la prise en compte de l'irréversibilité modifie considérablement les choses. Nous avons retenu trois choix proposés dans la littérature.

### 2.5.1 Minimisation sous contrainte

Suivant Marigo et al. [8], il faut en réalité minimiser le potentiel $\Pi(\mathbf{u}, d)$ sous la contrainte d'endommagement croissant :

$$ d = \operatorname{Arg}\{\inf_{d^*} \mathbf{R}^d(\mathbf{u}, d^*)\} \text{ sous la contrainte } d(t) \geq d(s), \forall s < t \quad (13) $$

L'équation (12) est alors remplacée par l'inégalité :

$$ \begin{cases} \frac{G_c}{l} (d + l^2 \Delta d) \geq 2(1-d) \psi_e^+ & \text{dans } \Omega \\ \nabla d \cdot \mathbf{n} = 0 & \text{sur } \partial \Omega \end{cases} \quad (14) $$

En pratique, cette contrainte traduisant la condition d'irréversibilité est implantée dans Cast3M [9] en adjoignant à chaque nœud une condition unilatérale sur la valeur nodale de l'endommagement via l'opérateur 'BLOC' 'MINI'. Les conditions unilatérales sont gérées par une méthode de statuts. La méthode des doubles Lagrange est alors utilisée, l'activation ou la désactivation d'un statut étant alors une modification mineure d'un terme de la matrice de raideur.

Le coût de cette approche est cependant encore très important. Notre implémentation actuelle pourrait être nettement améliorée en conservant les statuts actifs d'une itération sur l'autre.

### 2.5.2 Projection a posteriori

Une autre façon directe pour imposer la condition irréversible est d'utiliser une projection a posteriori [10] :

$$ \begin{cases} d = \operatorname{Arg}\{\inf_{d^*} \mathbf{R}^d(\mathbf{u}, d^*)\} \\ d(t) = \max(d, d(s)), \forall s < t \end{cases} \quad (15) $$

Cette méthode est évidemment beaucoup plus efficace numériquement, car la condition d'inégalité paraît et le problème de champ de phase devient linéaire.
Les solutions finales trouvées par cette méthode sont similaires à la minimisation sous contraintes [8]. C'est également notre constat sur les cas présentés, mais rien ne garantit que ceci se généralise à des cas de propagations complexes et fissurations multiples.

### 2.5.3 Introduction d'une variable d'histoire

Miehe et al. [5] propose d'utiliser une fonction d'histoire de l'énergie maximale à la place de l'énergie élastique dans l'équation de champ de phase :

$$ \mathcal{H}(x, t) = \max_{s \in [0, t]} \psi_e^+(x, s) \qquad (16) $$

L'équation à résoudre *in fine* est donc :

$$ \left\{ \begin{array}{l} \frac{G_c}{l}(d + l^2 \Delta d) = 2(1-d) \mathcal{H}(x,t) \quad \text{dans } \Omega \\ \nabla d \cdot \mathbf{n} = 0 \quad \text{sur } \partial \Omega \end{array} \right. \qquad (17) $$

L'irréversibilité est donc imposée de manière indirecte. Cette approche est aujourd'hui celle ayant le plus grand succès, grâce à son efficacité et sa simplicité d'implémentation. Dans les cas présentés, les solutions obtenues sont également très proches de ce qui est obtenu avec la projection *a posteriori*.

### 2.5.4 Méthode retenue pour ce papier

L'ensemble des méthodes présentées a été implémenté. Dans le cadre de ce papier, on se limitera à présenter et comparer les résultats obtenus avec la méthode de projection *a posteriori* (15) et celle de Miehe (16).

## 3 Traitement numérique des approches par champ de phase

Les approches par champ de phase de la fissuration fragile amènent à considérer un système de deux équations aux dérivées partielles couplées : la première décrit l'équilibre mécanique, la seconde le champ de phase.

Deux méthodes de résolution sont adoptées dans la littérature :

- une résolution monolithique où les deux équations sont résolues en même temps. Cette approche est rapportée comme efficace numériquement mais peu robuste sans l'utilisation de méthode de globalisation [11], en particulier lors des phases de propagation instable. Cette approche ne pouvant être implantée dans la version actuelle de Cast3M, elle n'est pas considérée ici.

- une résolution échelonnée où chaque équation est résolue successivement [6]. En effet, chaque problème pris individuellement se résume à un problème de minimisation d'une fonctionnelle convexe, a priori bien posé. Les résolutions échelonnées se divisent en deux catégories :

- La première, que nous appellerons purement implicite, s'appuie sur une méthode de point fixe (éventuellement accélérée) pour s'assurer qu'à convergence, les deux équations sont sa-tisfaites par la solution de fin de pas.
— L'autre catégorie, que nous appellerons explicite, considère qu'au moins une des équations peut ne pas être satisfaite en fin de pas [5, 12, 13]. Notre expérience montre que les résultats sont alors très dépendants du pas de temps, ce qui rend cette façon de procéder incompatible avec nos applications.

La solution retenue dans le travail précédent est une approche implicite échelonnée dont l'originalité consiste à résoudre l'équation de champ de phase pendant les itérations d'équilibre de la mécanique. Si l'on perd les bonnes propriétés des algorithmes alternés où chaque problème est résolu entièrement, cette approche permet de conserver certains éléments intéressants de l'algorithmes de résolution de Cast3M : l'accélération de convergence [14] et la procédure dite de convergence forcée.

Dans cette étude, nous proposons des pistes d'amélioration pour cette approche. Un schéma, dit semi-implicite, est proposé. L'idée est d'activer une résolution explicite lorsqu'une forte instabilité est détectée au cours des itérations. Cette méthode est basée sur l'observation qu'un grand nombre d'itérations mécaniques est souvent nécessaire après la convergence de l'endommagement. La résolution explicite avec cet endommagement convergé peut évidemment accélérer cette procédure. Pour illustrer la performance de cette méthode, nous allons comparer les différents algorithmes dans les cas test de traction et cisaillement dans un premier temps. Les différentes manières d'imposer la condition d'irréversibilité dans l'équation de champ de phase seront également comparées.

# 4 Implémentation numérique

## 4.1 Discrétisation spatiale

La méthode des éléments finis est utilisée pour discrétiser les problèmes de mécanique et d'endommagement :

$$ \mathbf{u}(x) = \mathbf{N}_u^T(x)\mathbf{U}, \quad d(x) = \mathbf{N}_d^T(x)\mathbf{d} \qquad (18) $$

où $\mathbf{U}$ et $\mathbf{d}$ sont respectivement les déplacements et endommagement nodaux, $\mathbf{N}_u$ et $\mathbf{N}_d$ sont les fonctions de forme associées. Notant $\mathbf{B}$ le gradient de la fonction de forme, la discrétisation éléments finis conduit aux deux équations suivantes (19) et (20) :

$$ \mathbf{R}^u(\mathbf{u}, d) = 0 \quad \leftrightarrow \quad \begin{cases} \mathbf{F}_{\text{ext}} - \mathbf{F}_{\text{int}} = 0 \\ \mathbf{F}_{\text{int}} = \int_{\Omega} \mathbf{B}_u^T \boldsymbol{\sigma} dV, \quad \mathbf{F}_{\text{ext}} = \int_{\partial\Omega} \mathbf{N}_u^T \mathbf{t} \, dA \end{cases} \qquad (19) $$

et

$$ \mathbf{R}^d(\mathbf{u}, d) = 0 \quad \leftrightarrow \quad \begin{cases} \mathbf{K}\mathbf{d} = \mathbf{F}_d \\ \mathbf{K} = \int_{\Omega} \left(2\psi_e^+ + \frac{G_c}{l}\right) \mathbf{N}_d^T \mathbf{N}_d dV + \int_{\partial\Omega} G_c l \mathbf{B}_d^T \mathbf{B}_d dV, \quad \mathbf{F}_d = \int_{\partial\Omega} 2\mathbf{N}_d^T \psi_e^+ dV \end{cases} \qquad (20) $$

où $\psi_e^+$ devient $\mathcal{H}$ si on introduit la fonction d'histoire (16) pour imposer l'irréversibilité.

## 4.2 Implémentation dans Cast3M

L'implémentation des méthodes proposées a été réalisée dans Cast3M à l'aide de la procédure INCREPL [9], couplée avec le générateur de code MFront [15] pour la loi de comportement.

La procédure implicite peut se schématiser comme suit :
* Initialisation à l'instant $t_n : \mathbf{u}^0 = \mathbf{u}(t_{n-1})$, $d^0 = d(t_{n-1})$

* Itération $i \ge 0$ :

1. Prédiction de l'incrément de déplacement : $\Delta \mathbf{u}^{i+1}$

2. Calcul des variables internes avec MFront (loi de comportement) : $\epsilon, \sigma, \psi$

3. Prédiction de l'endommagement (implicite) : $d^{i+1}$

4. Vérification de convergence et accélération éventuelle de la convergence.

La prédiction de l'incrément de déplacement $\Delta \mathbf{u}^{i+1}$ se base :

- soit sur une correction de Newton (la matrice tangente cohérente est calculée à endommagement constant).

- soit sur une correction de point fixe (ce qui équivaut à une méthode de quasi-Newton utilisant la matrice d'élasticité).

En cas de non convergence, la méthode de la convergence forcée est activée, ce qui permet de passer les étapes de propagation instable.

Nous avons observé que cet algorithme pouvait conduire à de très nombreuses itérations, alors même que le champ d'endommagement devenait stationnaire. Cette observation a conduit à la modification suivante :

* Initialisation à l'instant $t_n : \mathbf{u}^0 = \mathbf{u}(t_{n-1})$, $d^0 = d(t_{n-1})$

* Itération implicite $i$ pour : $\mathbf{u}^{i+1}, d^{i+1}$

* Vérification de convergence :

1. Si $\max(\|d^{i+1} - d^i\|) \le \epsilon_d$ et $\|\mathbf{R}^u(\mathbf{u}^{i+1}, d^i)\| \le \epsilon_u$, alors $d(t_n) = d^{i+1}, \mathbf{u}(t_n) = \mathbf{u}^{i+1}$, Fin de boucle.

2. Si $\max(\|d^{i+1} - d^i\|) \le \epsilon_d$ et $\|\mathbf{R}^u(\mathbf{u}^{i+1}, d^i)\| > \epsilon_u$, alors conserver l'endommagement $d(t_n) = d^{i+1}$. Fin de boucle. Résolution mécanique à endommagement constant $d = d^{i+1}$, $\mathbf{u} = \text{Arg}\{\inf_{\mathbf{u}^*} \mathbf{R}^u(\mathbf{u}^*, d)\}$.

3. Sinon, $(\cdot)^i \leftarrow (\cdot)^{i+1}$ et répéter la résolution implicite.

Cette méthode permet d'améliorer notablement la performance de résolution purement implicite ainsi que le démontrent les exemples numériques de la section suivante.

# 5 Exemples numériques

Cette section illustre la performance de la méthode proposée. Les cas-test de Miehe [5] sont utilisés (voir FIGURE 1). Le premier test compare deux méthodes de résolution dans le cas en cisaillement. La condition d'irréversibilité utilisée est la méthode de Miehe (16). Les solutions finales sont similaires, mais le calcul semi-implicite est plus rapide (voir FIGURE 2 et TABLE 1). Le second test compare deux conditions d'irréversibilité. Peu de différence est constatée dans la solution finale (FIGURE 3). Nous avons constaté un effet de propagation instable plus important dans le test de traction : une chute instantanée dans la courbe de force-déplacement. Cela conduit à un temps de calcul implicite beaucoup plus long. Cependant, le temps CPU est significativement réduit (avec un facteur de 8) avec la méthode semi-implicite (voir TABLE 1). Remarquons que dans tous les cas-test, la discrétisation temporelle reste inchangée. Le fait de conserver la convergence d'endommagement dans le schéma semi-implicite permet d'utiliser les même pas de temps que l'implicite : par exemple pour le cas en traction $\Delta u = 6, 1 \times 10^{-4}$ mm pour les premiers 8 pas de temps et $\Delta u = 6, 1 \times 10^{-5}$ mm pour les 20 pas de temps restants, qui ne sont pas acceptables par les schéma explicites.
<figure><img src="image_2.png" /><figcaption>FIGURE 1 – Cas test de Miehe en traction et en cisaillement [5]</figcaption></figure>

<figure><img src="image_3.png" /><figcaption>FIGURE 2 – Comparaison des méthodes implicite et semi-implicite (cas en cisaillement)</figcaption></figure>

## 6 Conclusion

Un schéma semi-implicite robuste et efficace est proposé pour la résolution de champ de phase. La simulation quasi-statique des propagations instables est significativement accélérée. Les résultats sont très encourageants pour son application à des cas complexes, comme la simulation de la fissuration de combustibles. Les solutions finales obtenues avec les différentes conditions d'irréversibilité sont très similaires dans ces premiers exemples. Des applications de ces méthodes proposées dans les cas de combustibles et d'un VER (Volume Elémentaire Représentatif) de béton sont en cours de développement.

## Remerciements

Ce travail a été mené dans le cadre du projet PLEIADES, co-financé par le CEA (Commissariat à l'Énergie Atomique et aux Énergies Alternatives), EDF (Électricité de France) et Framatome, et du projet interdisciplinaire sur la simulation numérique du CEA.

## Références

[1] Helfer, T., Bary, B., Dang, T. T., Fandeur, O., Michel, B. (2017). Modélisation par champ de phase de la fissuration des matériaux fragiles : Aspects numériques et applications au combustible nucléaire oxyde. In 13ème colloque national en calcul des structures.
<figure><img src="image_4.png" /><figcaption>FIGURE 3 – Comparaison des méthodes implicite et semi-implicite et différentes conditions d’irréversibilité (cas en traction)</figcaption></figure>

<table><caption>TABLE 1 – Comparaison de temps de calcul</caption><thead><tr><th>Cas-test</th><th>Miehe impl</th><th>Miehe semi-impl</th><th>Projection semi-impl</th></tr></thead><tbody><tr><td>Cisaillement</td><td>190 min</td><td>150 min</td><td>-</td></tr><tr><td>Traction</td><td>1930 min</td><td>-</td><td>250 min</td></tr></tbody></table>

[2] Michel, B., Helfer, T., Ramière, I., Esnoul, C. (2018). A new numerical methodology for simulation of unstable crack growth in time independent brittle materials. Engineering Fracture Mechanics, 188, 126-150.

[3] Michel, B., Sercombe, J., Thouvenin, G., Chatelet, R. (2008). 3D fuel cracking modelling in pellet cladding mechanical interaction. Engineering Fracture Mechanics, 75(11), 3581-3598.

[4] Michel, B., Helfer, T., Ramière, I., Esnoul, C. (2016). 3D continuum damage approach for simulation of crack initiation and growth in ceramic materials. Key Engineering Materials, 713, 155-158.

[5] Miehe, C., Hofacker, M., Welschinger, F. (2010). A phase field model for rate-independent crack propagation : Robust algorithmic implementation based on operator splits. Computer Methods in Applied Mechanics and Engineering, 199(45-48), 2765-2778.

[6] Bourdin, B., Francfort, G. A., Marigo, J. J. (2000). Numerical experiments in revisited brittle fracture. Journal of the Mechanics and Physics of Solids, 48(4), 797-826.

[7] Francfort, G. A., Marigo, J. J. (1998). Revisiting brittle fracture as an energy minimisation problem. Journal of the Mechanics and Physics of Solids, 46(8), 1319-1342.

[8] Amor, H., Marigo, J. J., Maurini, C. (2009). Regularized formulation of the variational brittle fracture with unilateral contact : Numerical experiments. Journal of the Mechanics and Physics of Solids, 57(8), 1209-1229.

[9] CEA. Site Cast3M. (2016). Disponible à l'adresse : http://www-cast3m.cea.fr/

[10] Lancioni, G., Royer-Carfagni, G. (2009). The variational approach to fracture mechanics. A practical application to the French Panthéon in Paris. Journal of elasticity, 95(1-2), 1-30.

[11] Wick, T. (2017). Modified Newton methods for solving fully monolithic phase-field quasi-static brittle fracture propagation. Computer Methods in Applied Mechanics and Engineering, 325, 577-611.
[12] Nguyen, T. T., Yvonnet, J., Zhu, Q. Z., Bornert, M., Chateau, C. (2016). A phase-field method for computational modeling of interfacial damage interacting with crack propagation in realistic microstructures obtained by microtomography. Computer Methods in Applied Mechanics and Engineering, 312, 567-595.

[13] Molnár, G., Gravouil, A. (2017). 2D and 3D Abaqus implementation of a robust staggered phase-field solution for modeling brittle fracture. Finite Elements in Analysis and Design, 130, 27-38.

[14] Ramière, I., Helfer, T. (2015). Iterative residual-based vector methods to accelerate fixed point iterations. Computers and Mathematics with Applications, 70(9), 2210-2226.

[15] Helfer, T., Michel, B., Proix, J. M., Salvo, M., Sercombe, J., Casella, M. (2015). Introducing the open-source mfront code generator : Application to mechanical behaviours and material knowledge management within the PLEIADES fuel element modelling platform. Computers and Mathematics with Applications, 70(5), 994-1023.