# Simulation numérique par FFT de la fissuration des matériaux hétérogènes à inclusions gonflantes

Gérald Feugueur, Lionel Gélébart, Corrado Maurini, Sandrine Miro, Christophe Bourcier

► To cite this version:

Gérald Feugueur, Lionel Gélébart, Corrado Maurini, Sandrine Miro, Christophe Bourcier. Simulation numérique par FFT de la fissuration des matériaux hétérogènes à inclusions gonflantes. 25e Congrès Français de Mécanique, Aug 2022, Nantes, France. hal-04280120

HAL Id: hal-04280120

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

Submitted on 10 Nov 2023

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
# Simulation numérique par FFT de la fissuration des matériaux hétérogènes à inclusions gonflantes

G. FEUGUEUR<sup>a,b</sup>, L. GÉLÉBART<sup>a</sup>, C. MAURINI<sup>b</sup>, S. MIRO<sup>c</sup>, C. BOURCIER<sup>d</sup>

a. Université Paris-Saclay, CEA, Service de Recherches Métallurgiques Appliquées, 91191, Gif/Yvette, France - gerald.feugueur@cea.fr & Lionel.GELEBART@cea.fr

b. CNRS, Institut Jean Le Rond d'Alembert, Sorbonne Université, UMR 7190, 75005, Paris, France - corrado.maurini@sorbonne-universite.fr

c. CEA, DES, ISEC, DE2D, Université de Montpellier, Marcoule 30207, France - sandrine.miro@cea.fr

d. Université Paris-Saclay, CEA, Service de Thermo-hydraulique et de Mécanique des Fluides, 91191, Gif-sur-Yvette, France - Christophe.BOURCIER@cea.fr

...

**Mots clés :** AMITEX_FFTP, Simulation numérique, Fissuration, Matériaux hétérogènes, FFT

## Résumé :

Les verres nucléaires actuellement utilisés pour le confinement des produits de fission et actinides mineurs peuvent inclure un taux de charge allant jusqu'à 18.5%. Les matériaux vitrocéramiques envisagés pour cette application seraient une alternative intéressante qui permettrait de diminuer le volume de colis. Pendant le stockage, les inclusions de phases cristallines, riches en produits de fission, sont soumises à une auto-irradiation α provoquant un gonflement susceptible d'engendrer une fissuration de la matrice vitreuse. La mise en œuvre de simulations de grande taille devra permettre d'évaluer l'effet de la microstructure sur la fissuration du matériau en se rattachant à des quantités d'intérêt comme le temps avant la première fissuration ou la surface fissurée.

## Abstract :

The nuclear glasses currently used for the containment of fission products and minor actinides can include a loading rate of up to 18.5%. The glass-ceramic materials envisaged for this application would be an interesting alternative that would make it possible to reduce the volume of the packages. During storage, the crystalline phase inclusions, rich in fission products, are subjected to α self-irradiation causing swelling that may lead to cracking of the glass matrix. The implementation of large-scale simulations should make it possible to evaluate the effect of the microstructure on the cracking of the material by relating to quantities of interest such as the time to first cracking or the cracked surface.

## 1 Introduction

Les méthodes « FFT » [15, 16], sont particulièrement adaptées à la simulation du comportement mécanique de matériaux hétérogènes. En effet, comparés à l'utilisation de codes Élément Finis standards, les codes FFT sont souvent bien plus performants et très bien adaptés à une implémentation parallèle en mémoire distribuée. Initialement proposée pour des modèles de comportement locaux, linéaires ou
non-linéaires, l'utilisation de ces méthodes s'étend désormais au cadre des modèles non-locaux, tels que la plasticité ou l'endommagement à gradient [7, 13]. Par ailleurs, pour simuler la fissuration, les modèles de champ de phase [10, 4] font l'objet d'un intérêt croissant dans la communauté mécanicienne. L'implémentation du modèle de champ de phase proposé par Miehe, Hofacker et Welschinger [14] a été récemment mis en place au sein du code massivement parallèle AMITEX_FFTP [2, 7]. L'objet de cette communication est donc de présenter les premiers résultats et améliorations mises en place pour l'étude des matériaux hétérogènes à inclusions gonflantes.

## 2 Génération de microstructures et utilisation des voxels composites

Afin d'évaluer l'effet de la microstructure, il est nécessaire de générer de manière aléatoire ces dernières en jouant sur la forme, la densité et la taille des inclusions. Au CEA, des outils génériques ont été développés pour concevoir ces microstructures et les voxeliser. L'utilisation des voxels composites [12, 11] a montré un intérêt pour améliorer la qualité des solutions numériques [5]. Ainsi, dans une dernière sous-section, on cherche à évaluer l'effet de l'utilisation de ces voxels composites dans le cadre de l'application visée : le gonflement d'inclusions au sein d'une matrice.

### 2.1 Présentation des outils génériques développés au CEA

#### 2.1.1 Combs : un outil de génération de microstructures aléatoires

Combs est une bibliothèque Python open source, basée sur Salome. Elle fait appel au modeleur géométrique de Salome pour définir les formes des inclusions et les insérer dans le Volume Élémentaire Représentatif (VER). Le VER est défini par les dimensions de sa boîte, le type de formes à insérer (sphères, sphères multi-couches, ou toute autre forme créée dans Salome GEOM), le nombre de formes à insérer (ou sa fraction volumique), la distance minimale entre chaque inclusion (ou si le chevauchement est autorisé), et une condition de frontière périodique ou non.

Pour le placement, les inclusions sont triées de la plus grande à la plus petite. Pour chaque inclusion, un point est choisi aléatoirement pour définir le centre de l'inclusion. L'inclusion est ensuite tournée de manière aléatoire et si aucun obstacle chevauchant n'est trouvé parmi les obstacles à l'intérieur de ce voxel, l'inclusion est stockée, sinon, l'inclusion est tournée aléatoirement autour de l'obstacle, jusqu'à ce qu'une place vide soit trouvée [19]. Le placement s'arrête si toutes les inclusions ont été placées ou s'il ne parvient pas à en placer une après un certain nombre d'essais.

Une fois le placement effectué, plusieurs méthodes peuvent être utilisées pour générer un maillage du VER, en utilisant des algorithmes de maillage avancés fournis dans Salome : NETGEN (open source) ou Mesh Gems de Distene. La méthode habituelle pour mailler des formes dans Salome est de faire une partition avec tous les objets et de mailler la forme résultante. Cette méthode directe prend beaucoup de temps pour un grand nombre d'inclusions. Une autre méthode plus performante est employée ici (voir [3] pour plus de détails).

#### 2.1.2 Voxelize : un outil de voxelisation

Voxelize utilise VoxSurf, un voxeliseur surfacique (C++), et MEDCoupling, une librairie de manipulation de maillage et de champs (C++, Python), dans le but de générer rapidement des images voxelisées. Cet outil prend en entrée un fichier de maillage avec des groupes volumiques définissant les différents matériaux ou bien un champ par élément avec des numéros définissant les matériaux. En sortie, cet outil, fournit une image 3D (i.e. un maillage régulier d'éléments cubiques) avec un champ de matériaux défini
sur chaque voxel. Concernant les voxels qui sont traversés par une interface, appelés aussi voxels com-
posites, des fichiers décrivant ces voxels sont créés suivant la convention proposée pour leur utilisation
au sein du code AMITEX_FFTP [2]. Pour ces voxels, on détermine la fraction volumique des phases
contenues dans le voxel, et pour chaque portion d'interface présente au sein du voxel : l'aire, les vecteurs
normal et tangent, et le barycentre.

Dans le cadre d'une application sur le béton avec des inclusions polyédriques, la microstructure est dans
un premier temps générée par Combs. Le maillage ainsi obtenu est voxelisé avec une résolution de 100
voxels par côté en prenant en compte ou non l'option permettant d'identifier et de caractériser les voxels
composites qui sont représentés en blanc (cf. Fig. 1).

<figure><img src="image_2.png" /><figcaption>FIG. 1 – Maillage avec deux groupes volumiques (à gauche), image voxélisée obtenue avec une résolution de 100 voxels par direction (au centre) et image contenant des voxels composites (à droite)</figcaption></figure>

## 2.2 Application à la simulation numérique du gonflement des phases cristallines

### 2.2.1 Description de la microstructure et conditions aux limites

On considère ici une microstructure avec une seule inclusion (disque centré) et on se limite à un cadre 2D périodique (en pratique le calcul est 3D avec une seule couche de voxels, la condition de périodicité revient alors à un calcul en déformation plane (ou plane généralisée)). La boîte a pour côté 2 et l'inclusion a pour rayon 0.5. Concernant la discrétisation spatiale, dans cette étude, elle sera supposée régulière et paramétrable par $N_x$, la résolution de $N_x$ voxels par direction allant de 10 à $10^4$ voxels. Au niveau du chargement, on impose une contrainte moyenne $\langle \sigma \rangle$ nulle et une déformation libre $\epsilon^s$, simulant le gonflement de l'inclusion tel que :

$$ \epsilon_{ii}^s(\tilde{x}) = \epsilon^s = 0.1 \quad \forall i = x, y \text{ et } \tilde{x} \in \Omega^I \quad (1) $$

La figure ci-dessous présente à gauche l'image voxélisée de la microstructure pour une simulation sans voxels composites, et une image de la fraction volumique des phases en présence pour une simulation avec voxels composites (cf. Fig. 2). Pour cette étude, l'évaluation des fractions volumiques est réalisée par sous-échantillonnage.

Les propriétés élastiques de la matrice sont fixées. En revanche celles de l'inclusion sont paramétrées par le contraste mécanique C :

$$ \left\{ \begin{array}{ll} \lambda_1 = 0 & \text{et } \mu_1 = 2 \\ \lambda_0 = 0 & \text{et } \mu_0 = C\mu_1 \end{array} \right. \qquad (2) $$
<figure><img src="image_3.png" /><figcaption>FIG. 2 – Champ matériau (à gauche) et champ de fraction volumique (à droite)</figcaption></figure>

Le coefficient de Poisson $\nu$ est pris égal à 0 sur tout le domaine afin d'avoir une contrainte hors-plan nulle ce qui implique d'avoir le premier coefficient de Lamé $\lambda$ nul pour les deux matériaux.

### 2.2.2 Lois de comportements pour les voxels composites

Une fois les voxels composites identifiés et caractérisés, une loi d'homogénéisation doit être choisie pour établir leurs propriétés élastiques effectives. Dans AMITEX_FFTP [2], trois lois d'homogénéisation ont été mises en œuvre : Voigt, Reuss et Laminé. Ces lois sont utilisées pour déterminer le tenseur élastique effectif $\tilde{C}$ à partir des $N$ fractions volumiques de phases $\alpha_j$ et des propriétés élastiques $C_j$ de chaque phase $j$. La loi de Voigt suppose une déformation homogène sur toutes les phases. À l'inverse, la loi de Reuss suppose une contrainte homogène sur toutes les phases. Le modèle de Voigt fournit une borne supérieure en raideur du comportement tandis que le modèle de Reuss fournit une borne inférieur. En plus de ces bornes supérieures et inférieures, la règle d'homogénéisation des stratifiés est proposée. Elle fournit la solution exacte pour un matériau stratifié infini. Dans ce cas, les contraintes et les déformations sont constantes par phase. Les conditions d'équilibre et de compatibilité sont satisfaites par la continuité de la contrainte anti-plane $\sigma^{\bar{P}}$ et de la déformation dans le plan $\varepsilon^p$. L'ensemble d'équations correspondant est donné ci-dessous :

$$ \left\{ \begin{array}{l} \langle \sigma^{\bar{P}} \rangle = \sigma_j^{\bar{P}} \quad \forall j = 0, N-1 \\ \langle \varepsilon^p \rangle = \varepsilon_j^p \quad \forall j = 0, N-1 \\ \sigma_j = C_j : (\varepsilon_j - \varepsilon^s) \quad \forall j = 0, N-1 \end{array} \right. \quad \text{et} \quad \left\{ \begin{array}{l} \langle \sigma^p \rangle = \sum_{j=0}^{N-1} \alpha_j \sigma_j^p \\ \langle \varepsilon^{\bar{P}} \rangle = \sum_{j=0}^{N-1} \alpha_j \varepsilon_j^{\bar{P}} \end{array} \right. \quad (3) $$

La résolution de cet ensemble d'équations permet d'évaluer la contrainte moyenne sur un voxel composite pour une déformation moyenne imposée sur ce même voxel. Les $\varepsilon_j^p$ sont directement déterminées par la continuité de la déformation dans le plan $\varepsilon^p$. La continuité de la contrainte anti-plan $\sigma^{\bar{P}}$ combinée à la loi de comportement fournit un système linéaire où les $\varepsilon_j^{\bar{P}}$ sont les inconnues.

### 2.2.3 Effets du contraste et de l'utilisation des voxels composites

Le choix optimal d'une loi homogénéisée pour les voxels composites dépend fortement du contraste. En effet, le modèle de Voigt semble être, pour ce cas de figure, le plus intéressant pour une inclusion molle ($C < 1$), alors que le modèle de Reuss quant à lui, semble être à privilégie pour une inclusion raide ($C > 1$). En revanche, le modèle Laminé qui permet de prendre en compte l'orientation de l'interface semble être le plus intéressant quelque soit le niveau de contraste.

Pour un niveau de contraste important, l'état de déformation dans la matrice y est plus élevé. Le champ de déformation y est plus dispersée pour un niveau de contraste plus élevé (cf. Fig. 3).
<figure><img src="image_4.png" /><figcaption>FIG. 3 – Déformation axiale moyenne $\overline{\varepsilon}_{xx}$ (à gauche) et écart-type de la déformation axiale $\sigma_{\overline{\varepsilon}_{xx}}$ (à droite) sur la cellule</figcaption></figure>

Pour les voxels standards, on peut remarquer la présence d'oscillations au niveau de l'interface. L'utilisation d'une loi homogénéisée de type Voigt ou Reuss permet d'atténuer ces oscillations et améliore ainsi la qualité de la solution (cf. Fig. 4).

<figure><img src="image_5.png" /><figcaption>FIG. 4 – Champ de contrainte axiale $\sigma_{xx}$ sans la prise en compte des voxels composites (en haut à gauche), en utilisant une loi homogénéisée de type Voigt (en haut à droite), de type laminé (en bas à gauche) et de type Reuss (en bas à droite) pour les voxels composites avec $C = 2$ et $N_x = 25$</figcaption></figure>

# 3 Simulation numérique de la fissuration matricielle par champ de phase

## 3.1 Approche de la rupture fragile par champ de phase

Cette partie illustre comment un problème d'endommagement à gradient peut être résolu de manière alternée. L'approche employée ici est celle par champ de phase dans laquelle deux problèmes mécaniques différents sont impliqués :

- un problème de déplacement (problème en $u$) impliquant un modèle constitutif élastique potentiellement non linéaire à endommagement fixée qui comprend des conditions unilatérales
interdisant l'évolution des fissures dans des états de compression.

— un problème d'endommagement (problème en $\phi$) comprenant un terme de gradient d'endommagement associé à une longueur interne $l_0$ de régularisation qui est à la base de l'approche de la rupture par champ de phase.

## 3.2 Introduction de l'approche de la rupture fragile par champ de phase

L'approche de la rupture fragile par le champ de phase est issue de la régularisation mathématique de l'approche variationnelle de la rupture de Francfort et Marigo [10]. Suite aux travaux mathématiques d'Ambrosio et Tortorelli [1], la régularisation proposée par Bourdin, Francfort et Marigo [4] repose sur l'énergie totale $\mathcal{E}$ qui est la somme des énergies potentielle $\Psi_s$ et de rupture $\Psi_c$ s'écrivant comme suit :

$$ \mathcal{E}(\mathbf{u}, \phi) = \Psi_s(\mathbf{u}, \phi) + \Psi_c(\phi) = \int_{\Omega} \Psi(\varepsilon(\mathbf{u}), \phi) \, dV + \int_{\Omega} \frac{G_c}{c_{\alpha}} \left( \frac{\alpha(\phi)}{l_0} + l_0 \|\nabla\phi\|^2 \right) \, dV \quad (4) $$

Les solutions lors de l'évolution quasi-statique à des pas de temps discrets ($\mathbf{u}^n, \phi^n$) sont obtenues à partir du principe de minimisation de l'énergie global sur les deux variables avec une contrainte d'irréversibilité supplémentaire $\phi^n > \phi^{n-1}$ pour la variable d'endommagement.

Cette formulation initiale ne fait pas de distinction entre les états de contrainte en traction et en compression, ce qui entraîne potentiellement la formation de fissures parasites dans des régions en compression. Un rempart couramment utilisé consiste à séparer la densité d'énergie $\Psi_0$ en une partie positive (associée aux états de traction) et une partie négative (états de compression) et à appliquer la fonction de dégradation $g$ uniquement sur la partie positive, à savoir :

$$ \Psi_s(\mathbf{u}, \phi) = \int_{\Omega} g(\phi) \Psi_0^+(\varepsilon(\mathbf{u})) + \Psi_0^-(\varepsilon(\mathbf{u})) \, dV \quad (5) $$

Les deux parties sont calculées pour les solides isotropes par :

$$ \Psi_0^{\pm} = \frac{\lambda}{2} \langle \mathrm{tr}(\varepsilon) \rangle_{\pm}^2 + \mu (\varepsilon^{\pm} : \varepsilon^{\pm}) \quad (6) $$

où $\lambda$ et $\mu$ sont les coefficients de Lamé, $\langle x \rangle_{\pm} = (x \pm |x|)/2$, et $\varepsilon^{\pm}$ représente la décomposition spectrale introduit par Miehe, Hofacker et Welschinger [14] du tenseur de déformation $\varepsilon^{\pm} = \sum_{i=1}^3 \langle \varepsilon^i \rangle_{\pm} n^i \otimes n^i$, avec $\varepsilon^i$ et $n^i$ étant les valeurs propres et les vecteurs propres du tenseur de déformation.

Un choix classique pour la fonction de dégradation est $g(\phi) = (1 - \phi)^2 + k$ où $k \ll 1$ est une rigidité résiduelle évitant la rigidité singulière dans l'état complètement endommagé. Un choix couramment utilisé pour la densité d'énergie de rupture est le modèle AT2 ($\alpha(\phi) = \phi^2$ et $c_{\alpha} = 2$) [1]. Cette formulation est couramment utilisée, notamment dans les travaux de Miehe, Hofacker et Welschinger [14].

Le problème couplé peut être résolu en minimisant la fonctionnelle d'énergie (Éq. 4), ce qui conduit aux équations locales suivantes :

$$ \frac{G_c}{l_0} (\phi - l_0^2 \Delta \phi) = 2 (1 - \phi) \mathcal{H}(\varepsilon) \quad (7) $$
$$
\left\{
\begin{array}{l}
\mathrm{div}(\boldsymbol{\sigma}(\boldsymbol{u}, \phi)) = 0 \\
\boldsymbol{\sigma} = g(\phi) (\lambda \langle \mathrm{tr}(\varepsilon) \rangle_{+} \mathbf{I} + 2\mu \varepsilon^{+}) + (\lambda \langle \mathrm{tr}(\varepsilon) \rangle_{-} \mathbf{I} + 2\mu \varepsilon^{-}) \\
\varepsilon = (\nabla \boldsymbol{u} + \nabla^{T} \boldsymbol{u}) / 2
\end{array}
\right.
\qquad (8)
$$

où $\boldsymbol{\sigma}(\boldsymbol{u}, \phi)$ est la contrainte locale, qui suit une loi comportement du matériau dépendant du champ de déformation $\varepsilon$ et du champ de phase $\phi$. L'énergie élastique $\mathcal{H}(\varepsilon) = \Psi_0^+(\varepsilon)$ joue le rôle de terme source dans l'équation 7. Pour s'assurer que le champ d'endommagement $\phi$ ne décroît pas, ce terme source est modifié en prenant le maximum de l'énergie élastique de la manière suivante :

$$
\mathcal{H}(\mathbf{x}, t) = \max_{s \in [0, t]} \Psi_0^+(\varepsilon(\mathbf{x}, s)) \quad (9)
$$

3.3 Implémentation du schéma alterné en FFT

Le schéma alterné de Miehe, Hofacker et Welschinger [14] est utilisé pour résoudre séparément le pro-
blème de champ de phase (Éq. 7) et le problème mécanique (Éq. 8). Le schéma alterné consiste à résoudre
le problème du champ de phase en considérant le champ mécanique connu puis à faire évoluer le champ
mécanique à champ de phase donné. Ce processus peut être itératif ou non. Dans ce dernier cas, une
analyse de convergence vis-à-vis du pas de chargement est indispensable et peut conduire à des pas de
temps très faibles. On se propose ci-dessous de résoudre ces problèmes à l'aide de méthodes FFT.

3.3.1 Résolution du problème mécanique

L'inconnue dans le problème mécanique (Éq. 8) est la déformation au lieu du déplacement. Le champ d'endommagement est considéré comme un paramètre à un pas de temps donné. La solution du problème mécanique (Éq. 8) peut être obtenue par résolution itérative de l'équation de Lippmann Schwinger :

$$
\epsilon(\mathbf{x}) = -\mathbf{\Gamma}^0 * \tau(\mathbf{x}) + \mathbf{E} \quad \text{avec} \quad \tau(\mathbf{x}) = \sigma(\mathbf{x}) - \mathbf{C}_0 : \epsilon(\mathbf{x}) \quad (10)
$$

où $\tau$ est le champ de polarisation, et $\Gamma^0$ est l'opérateur de Green d'ordre 4 pour la déformation qui est calculé dans l'espace de Fourier en fonction des constantes élastiques du matériau de référence $C_0$ choisi comme dans [16] permettant d'obtenir un taux de convergence optimal pour la solution mécanique. Le champ de contrainte $\sigma$ est calculé dans l'espace réel en utilisant la loi de comportement de l'équation 8, ce qui permet de prendre en compte les effets de l'évolution de l'endommagement et de la condition unilatérale. Le schéma itératif de point fixe donné par l'équation 10 correspond au schéma de base proposé initialement par Moulinec et Suquet [16]. L'utilisation d'un opérateur de Green modifié [18] s'avère ici strictement équivalent à l'utilisation d'éléments finis hexaédriques à intégration réduite dans un code Élément Finis standard [17]. Pour s'en convaincre, toujours dans le cadre d'une inclusion gonflante, on remarque que les champs de contrainte $\sigma_{xx}$ sont visuellement identiques quelque soit la méthode employée (une analyse quantitative des écarts confirme l'équivalence des deux méthodes, à la précision près des algorithmes itératifs utilisés) (cf. Fig. 5).

Une procédure d'accélération de la convergence [6] est également utilisée pour améliorer les performances de l'algorithme du point fixe.
<figure><img src="image_6.png" /><figcaption>FIG. 5 – Champ de contrainte axiale $\sigma_{xx}$ obtenu avec la méthode FFT (AMITEX_FFTP [2] à gauche), et la méthode EF (FEniCSx [9] à droite) pour $C = 2$ et $N_x = 25$</figcaption></figure>

### 3.3.2 Résolution du problème d'endommagement

Pour résoudre le problème du champ de phase, un algorithme itératif avec un point fixe imbriqué dans le schéma alterné est employé en utilisant la technique FFT [7]. L'objectif de cet algorithme itératif est de trouver le champ de phase $\phi^{n+1}(\boldsymbol{x})$ au pas de temps $t_{n+1}$ sur la base de la variable histoire $\mathcal{H}^n(\boldsymbol{x})$ et du champ de phase $\phi^n(\boldsymbol{x})$ qui sont connus au pas de temps précédent $t_n$.

La prise en compte de l'hétérogénéité a été proposée par Ernesti [8] en introduisant une fonction caractéristique $\theta$ au niveau de la fonction de dégradation $g$, où dans les zones endommageables, cette fonction $\theta$ vaut 1 et à l'inverse, dans les zones non endommageables, cette dernière vaut 0. L'expression proposée par Ernesti a été légèrement modifiée pour se ramener à une forme plus couramment utilisée pour la fonction de dégradation $g$ :

$$g(\boldsymbol{x}, \phi) = \theta(\boldsymbol{x}) (k + (1 - \phi)^2) + (1 - \theta(\boldsymbol{x})) \quad (11)$$

Afin de définir un champ de polarisation et un opérateur de Green pour le problème de champ de phase, il est nécessaire de réécrire l'équation 7 sous la forme :

$$A^n(\boldsymbol{x})\phi(\boldsymbol{x}) - \Delta\phi(\boldsymbol{x}) = B^n(\boldsymbol{x}) \quad (12)$$

En introduisant un paramètre homogène $A_0^n$, l'équation 12 devient :

$$A_0^n \phi(\boldsymbol{x}) - \Delta \phi(\boldsymbol{x}) = \chi(\boldsymbol{x}) \quad (13)$$

avec $\chi(\boldsymbol{x})$, le champ défini ci-dessous :

$$\chi(\boldsymbol{x}) = B^n(\boldsymbol{x}) - (A^n(\boldsymbol{x}) - A_0^n)\phi(\boldsymbol{x}) \quad (14)$$

L'algorithme du point fixe est ensuite utilisé pour résoudre l'équation 13 afin de trouver un nouveau champ d'endommagement $\phi^{i+1}(\boldsymbol{x})$, en supposant que le champ $\chi^i(\boldsymbol{x})$ est connu au pas d'itération précédent $i$ (cf. Algorithme 1). L'opérateur Laplacien est coûteux en calcul dans l'espace réel, mais dans l'espace de Fourier, il peut être facilement calculé par multiplication locale avec le vecteur d'onde $\xi$. Ainsi, la solution de l'équation 13 peut être exprimée dans l'espace de Fourier de la manière suivante :

$$\hat{\phi}(\boldsymbol{\xi}) = \frac{\hat{\chi}(\boldsymbol{\xi})}{A_0^n + \boldsymbol{\xi} \cdot \boldsymbol{\xi}} \quad (15)$$
**Algorithme 1** Résolution du problème d'endommagement par un schéma de type point-fixe :

Entrées : $\phi^n(x)$ et $\mathcal{H}^n(x)$

Sortie : $\phi^{n+1}(x)$

Initialisation : $\phi^i(x) = \phi^n(x)$

Tant que $\epsilon > \epsilon^{\text{tol}}$ faire

Espace réel

$$\chi^i(x) = B^n(x) - (A^n(x) - A_0^n)\phi^i(x)$$

$$\phi^{i+1}(x)$$

$$\chi^{i+1}(x) = B^n(x) - (A^n(x) - A_0^n)\phi^{i+1}(x)$$

$$\epsilon = \frac{\|\chi^{i+1}(x) - \chi^i(x)\|_2}{\|\chi^{i+1}(x)\|_2}$$

$\xrightarrow{\text{FFT}}$

$\xleftarrow{\text{(FFT)}^{-1}}$

Espace de Fourier

$$\hat{\chi}^i(\xi)$$

$$\hat{\phi}^{i+1}(\xi) = \frac{\hat{\chi}^i(\xi)}{A_0^n + \xi \cdot \xi}$$

Fin tant que

Une fois le champ d'endommagement déterminé à un incrément donné, ce champ est ensuite injecté pour résoudre le champ mécanique (en l'occurrence la déformation $\epsilon$). L'algorithme du point fixe est alors utilisé pour résoudre l'équation 10 (cf. Algorithme 2).

**Algorithme 2** Résolution du problème mécanique par un schéma de type point-fixe :

Entrées : $\epsilon^n(x)$, $\phi^{n+1}(x)$ et $\Delta E^{n+1}$

Sortie : $\mathcal{H}^{n+1}(x)$

Initialisation : $\epsilon^i(x) = \epsilon^n(x) + \Delta E^{n+1}$

Tant que $\epsilon > \epsilon^{\text{tol}}$ faire

Espace réel

$$\sigma^i(x) = g(\phi^{n+1}(x), \theta(x))C(x) : \epsilon^i(x)$$

$$\tau^i(x) = \sigma^i(x) - C_0(x) : \epsilon^i(x)$$

$$\epsilon^{i+1}(x)$$

$$\sigma^{i+1}(x) = g(\phi^{n+1}, \theta(x))C(x) : \epsilon^{i+1}(x)$$

$\xrightarrow{\text{FFT}}$

$\xleftarrow{\text{(FFT)}^{-1}}$

Espace de Fourier

$$\hat{\tau}^i(\xi)$$

$$\hat{\epsilon}^{i+1}(\xi) = -\hat{\Gamma}^0(\xi) \cdot \hat{\tau}^i(\xi)$$

$$\epsilon = \frac{(\langle\|\xi \cdot \hat{\sigma}(\xi)\|^2\rangle)^{1/2}}{\|\hat{\sigma}(0)\|}$$

Fin tant que

## 3.4 Simulation numérique de la fissuration matricielle

### 3.4.1 Présentation du cas d'étude

On considère ici la même microstructure que précédemment avec une seule inclusion centrée. La boîte a pour côté 100 µm et l'inclusion a pour rayon 25 µm. La discrétisation spatiale, régulière, est donnée par $N_x = 79$ (avec un voxel dans l'épaisseur). Au niveau du chargement, on impose toujours une contrainte moyenne nulle et une déformation libre simulant le gonflement de l'inclusion en 2D avec un coefficient de Poisson nul toujours pour la même raison.

Les propriétés élastiques sont identiques sur tout le domaine et les propriétés liées au modèle de champ de phase sont appliquées au niveau de l'inclusion gonflante (matériau indicé par 0) :

$$\left\{ \begin{array}{l} \lambda_0 = \lambda_1 = 0 \text{ GPa et } \mu_0 = \mu_1 = 34 \text{ GPa} \\ G_c = 8 \text{ J/m}^2 \text{ et } l_0 = 3 \text{ µm} \end{array} \right. \qquad (16)$$
### 3.4.2 Rupture matricielle sous l'effet du gonflement

Dans ce cas test, on observe ici l'apparition d'une fissure en croix puis une propagation autour de l'inclusion (cf. Fig. 6).

<figure><img src="image_7.png" /><figcaption>FIG. 6 – Champs d'endommagement à différents instants situés à droite de la réponse macroscopique en fonction du chargement</figcaption></figure>

Le schéma de différenciation DF1 qui s'apparente à des éléments finis hexaédriques à intégration réduite avec un décalage d'un 1/2 voxel pour l'une des deux grilles (le champ de phase est aux nœuds et les champs de contrainte/déformation au niveau du centre des voxels) fournit un champ oscillant qui ne sera pas conservé (cf. Fig. 7). Deux autres schémas de différenciation ont été proposés : l'un dont la dérivée seconde par direction est évaluée à partir des deux voisins dans la même direction (DF2) et l'autre qui est le reflet du schéma DF1 sans décalage pour lequel il est donc nécessaire d'interpoler $\phi$ (des noeuds vers les centres de voxels) pour l'utiliser dans la loi de comportement du problème mécanique, et d'interpoler la variables $\Psi_0^+$ (des centres vers les noeuds de voxels) pour l'utiliser dans le problème de champ de phase (DF3).

<figure><img src="image_8.png" /><figcaption>FIG. 7 – Champs d'endommagement pour le schéma de différenciation DF1 (à gauche), DF2 (au centre) et DF3 (à droite)</figcaption></figure>
Dans ce cas test, l'utilisation ou non de la variable histoire a finalement peu d'effet sur la réponse macroscopique et sur le trajet de fissuration. En revanche, sans l'effet histoire il y a une guérison qui apparaît au voisinage de la fissure (cf. Fig. 8).

<figure><img src="image_9.png" /><figcaption>FIG. 8 – Champs d'endommagement pour le schéma de différenciation DF2 avec (au centre) et sans (à droite) la variable histoire et réponse macroscopique (à gauche) en fonction du chargement</figcaption></figure>

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

# 4 Conclusions et perspectives

La génération de microstructures constitue un axe majeur pour évaluer l'effet de ces dernières sur la fissuration matricielle. Les outils génériques développés au CEA pour générer les microstructures, les voxeliser et identifier les voxels composites permettent dorénavant de réaliser ce genre d'évaluation. Il serait par ailleurs judicieux de tendre vers des microstructures réelles issues d'images MEB-FIB 3D après avoir analysé les effets de microstructures aléatoires virtuelles. Les voxels composites quant à eux, permettent d'améliorer la qualité des simulations numériques avec des champs plus lisses au voisinage de l'interface. L'implémentation du modèle de champ de phase proposé par Miehe, Hofacker et Welschinger a été récemment mis en place au sein du code massivement parallèle AMITEX_FFTP. Ce modèle de champ de phase est également en cours d'implémentation au sein du code de calcul par Éléments Finis FEniCSx. En effet, la prise en compte de ces modèles d'endommagement dans les simulations étant relativement délicate, il semble important de pouvoir confronter les résultats des deux codes. Pour ce faire, un premier cas test 2D a été mis en place côté AMITEX_FFTP, et est en cours côté FEniCSx. La prise en compte des voxels composites dans le modèle de champ de phase sera à implémenter et à investiguer. L'évaluation de l'effet de la distribution, la forme et la taille des inclusions sur le comportement à rupture des matériaux vitrocéramiques sera par la suite à analyser et la sensibilité de quantités d'intérêts, comme le temps à rupture et la surface fissurée, aux paramètres microstructuraux devra être quantifiée.

# Références

[1] Luigi Ambrosio and Vincenzo Maria Tortorelli. Approximation of functional depending on jumps by elliptic functional via $\Gamma$-convergence. *Communications on Pure and Applied Mathematics*, 43(8):999–1036, 1990.

[2] AMITEX_FFTP. https://amitexfftp.github.io/AMITEX.
[3] C Bourcier, W Dridi, L Chomat, E Laucoin, B Bary, and E Adam. Combs : open source python library for rve generation. application to microscale diffusion simulations in cementitious materials. In *SNA+ MC 2013-Joint International Conference on Supercomputing in Nuclear Applications+ Monte Carlo*, page 02107. EDP Sciences, 2014.

[4] Blaise Bourdin, Gilles A Francfort, and Jean-Jacques Marigo. Numerical experiments in revisited brittle fracture. *Journal of the Mechanics and Physics of Solids*, 48(4):797-826, 2000.

[5] Romain Charière, Aldo Marano, and Lionel Gélébart. Use of composite voxels in FFT based elastic simulations of hollow glass microspheres/polypropylene composites. *International Journal of Solids and Structures*, 182:1-14, 2020.

[6] Yang Chen, Lionel Gélébart, Camille Chateau, Michel Bornert, Cédric Sauder, and Andrew King. Analysis of the damage initiation in a sic/sic composite tube from a direct comparison between large-scale numerical simulation and synchrotron x-ray micro-computed tomography. *International Journal of Solids and Structures*, 161:111-126, 2019.

[7] Yang Chen, Dmytro Vasiukov, Lionel Gélébart, and Chung Hae Park. A FFT solver for variational phase-field modeling of brittle fracture. *Computer Methods in Applied Mechanics and Engineering*, 349:167-190, 2019.

[8] Felix Ernesti, Matti Schneider, and Thomas Böhlke. Fast implicit solvers for phase-field fracture problems on heterogeneous microstructures. *Computer Methods in Applied Mechanics and Engineering*, 363:112793, 2020.

[9] FEniCSx. https://fenicsproject.org/.

[10] Gilles A Francfort and J-J Marigo. Revisiting brittle fracture as an energy minimization problem. *Journal of the Mechanics and Physics of Solids*, 46(8):1319-1342, 1998.

[11] Lionel Gélébart and Franck Ouaki. Filtering material properties to improve FFT-based methods for numerical homogenization. *Journal of Computational Physics*, 294:90-95, 2015.

[12] Matthias Kabel, Dennis Merkert, and Matti Schneider. Use of composite voxels in FFT-based homogenization. *Computer Methods in Applied Mechanics and Engineering*, 294:168-188, 2015.

[13] Aldo Marano, Lionel Gélébart, and Samuel Forest. FFT-based simulations of slip and kink bands formation in 3d polycrystals : influence of strain gradient crystal plasticity. *Journal of the Mechanics and Physics of Solids*, 149:104295, 2021.

[14] Christian Miehe, Fabian Welschinger, and Martina Hofacker. Thermodynamically consistent phase-field models of fracture : Variational principles and multi-field fe implementations. *International Journal for Numerical Methods in Engineering*, 83(10):1273-1311, 2010.

[15] Hervé Moulinec and Pierre Suquet. A fast numerical method for computing the linear and nonlinear mechanical properties of composites. *Comptes Rendus de l'Académie des sciences. Série II. Mécanique, physique, chimie, astronomie*, 1994.

[16] Hervé Moulinec and Pierre Suquet. A numerical method for computing the overall response of nonlinear composites with complex microstructure. *Computer methods in applied mechanics and engineering*, 157(1-2):69-94, 1998.

[17] Matti Schneider, Dennis Merkert, and Matthias Kabel. FFT-based homogenization for microstructures discretized by linear hexahedral elements. *International journal for numerical methods in engineering*, 109(10):1461-1489, 2017.
[18] François Willot. Fourier-based schemes for computing the mechanical response of composites with accurate local fields. *Comptes Rendus Mécanique*, 343(3):232-245, 2015.

[19] Peter Wriggers and S.O Moftah. Mesoscale models for concrete : Homogenisation and damage behaviour. *Finite elements in analysis and design*, 42(7):623-636, 2006.