# Simulation de la fissuration par champ de phase : Analyse technique de l'initialisation et la propagation de fissures

Flavien Loiseau, Edgar Zembra, Hervé Henry, Véronique Lazarus

► To cite this version:

Flavien Loiseau, Edgar Zembra, Hervé Henry, Véronique Lazarus. Simulation de la fissuration par champ de phase : Analyse technique de l'initialisation et la propagation de fissures. Congrès Français de Mécanique 2025, Aug 2025, Metz, France. hal-05221510

HAL Id: hal-05221510

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

Submitted on 25 Aug 2025

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
# Simulation de la fissuration par champ de phase :
Analyse technique de l'initialisation et la
propagation de fissures

F. Loiseau<sup>a</sup>, E. Zembra<sup>b</sup>, H. Henry<sup>b</sup>, V. Lazarus<sup>a</sup>

a. IMSIA, ENSTA, EDF, CNRS, IP Paris, Palaiseau, 91120, France.

(flavien.loiseau@ensta.fr, veronique.lazarus@ensta.fr)

b. Laboratoire PMC, École Polytechnique, CNRS, IP Paris, 91120 Palaiseau, France.

(edgar.zembra@polytechnique.edu, herve.henry@polytechnique.edu)

## Résumé :

Cette étude propose d'étudier les bias numériques dans les simulations de la fissuration par champ de phase. Plus spécifiquement, elle examine l'influence du maillage sur la trajectoire de fissure. Un problème de référence basé sur un essai de Pure Shear avec fissure initiale excentrée est utilisé pour évaluer l'impact de la discrétisation spatiale sur la trajectoire de fissure prédite. Les simulations par champ de phase sont comparées à une solution de référence issue de la mécanique linéaire élastique de la rupture, vers laquelle elles convergent théoriquement. Les résultats montrent que les maillages structurés introduisent une anisotropie artificielle qui contraint la fissure à suivre une trajectoire dépendante du maillage. En revanche, les maillages non structurés permettent une prédiction plus fidèle de la trajectoire exponentielle attendue. En effet, le caractère désordonné du maillage réduit les effets d'anisotropie artificiels induits par la discrétisation. Ces observations mettent en évidence l'importance du choix du maillage dans les simulations de propagation de fissures. Des recommandations pratiques sont finalement proposées pour minimiser les biais numériques dans les simulations par champ de phase.

## Abstract :

This study examines numerical biases in phase-field fracture simulations. More specifically, it investigates the influence of meshing on crack path trajectories in phase-field fracture simulations. A reference problem based on a Pure Shear test with an eccentric initial crack is used to evaluate the impact of spatial discretization on predicted crack trajectories. Phase-field simulations are compared to a reference solution from Linear Elastic Fracture Mechanics, towards which they theoretically converge. Results show that structured meshes introduce artificial anisotropy that forces the crack to follow a mesh-dependent trajectory. Conversely, unstructured meshes enable a more accurate prediction of the expected exponential trajectory. Indeed, the disordered nature of the mesh significantly reduces the artificial anisotropy effects induced by discretization. These observations highlight the importance of mesh selection in crack propagation simulations. Practical recommendations are finally proposed to minimize numerical biases in phase-field simulations.

**Mots clefs :** Mécanique Linéaire Élastique de la Rupture, Approche variationnelle par champ de phase, Simulation numérique.
# 1 Introduction

L'étude de la propagation de fissures est cruciale pour garantir la fiabilité et la durabilité des structures mécaniques. La capacité à prédire précisément l'évolution d'une fissure permet d'optimiser le dimensionnement des structures et d'éviter les défaillances catastrophiques. La validation expérimentale des modèles de propagation de fissures sur des cas simples, traitables analytiquement, constitue une première étape essentielle. Cependant, pour les cas réalistes plus complexes, l'implémentation de ces modèles dans des simulations numériques devient incontournable. Les méthodes numériques, bien que puissantes, peuvent introduire des biais plus ou moins importants dans les prédictions du modèle. Il est donc essentiel de mettre en évidence ces biais et de proposer des méthodes pour les atténuer ou les éviter. Bien que des tests de convergence soient couramment réalisés, ils se limitent généralement aux paramètres de discrétisation spatiale et temporelle.

Ainsi, cette étude vise à identifier différents biais numériques pouvant affecter la simulation de la propagation de fissures par champs de phase [1, 2]. À terme, l'objectif est de proposer des solutions pratiques pour atténuer, voire éliminer, ces biais. L'approche variationnelle de la rupture par champ de phase étant une régularisation de la mécanique de la rupture, l'idée est de s'appuyer sur cette dernière pour mettre en évidence les biais numériques dans les simulations par champ de phase. Deux aspects ont été étudiés dans le cadre de ce travail : d'une part, les techniques de représentation de fissures pré-existantes, et d'autre part, l'influence du maillage sur la trajectoire de propagation.

Le premier aspect ayant fait l'objet d'une étude approfondie dans une publication antérieure [3], nous n'en présenterons pas les détails dans ce document. Cette étude a permis de montrer que la manière d'introduire une fissure dans les simulations par champ de phase n'est pas un choix anodin et qu'un mauvais choix peut introduire un fort biais, notamment une large surestimation du seuil de propagation. Ainsi, des recommandations pratiques ont été proposées pour éliminer ces biais : (a) la fissure initiale doit avoir une largeur d'un élément, et (b) le champ de phase doit être initialisé à l'état complètement fissuré le long des lèvres de la fissure.

Le second aspect porte sur l'influence du maillage sur la trajectoire de fissure. Dans cette étude, on considère un test numérique basé sur l'essai de Pure Shear avec une fissure initiale excentrée. Pour ce cas test, les simulations basées sur la mécanique de la rupture montre que la fissure se propage de manière telle que la pointe de fissure se recentre vers la mi-hauteur de l'échantillon de manière exponentielle durant la propagation. Ainsi, on compare les trajets de fissures prédits par champs de phase en utilisant différents maillages à ceux issus de la mécanique de la rupture pour identifier des biais numériques potentiels. On montre ainsi que le choix de la discrétisation spatiale peut significativement modifier la trajectoire de fissure prédite par la simulation.

La Section 2 présente le cadre analytique de cette étude, en introduisant la mécanique élastique linéaire de la rupture ainsi que sa régularisation par champ de phase. La Section 3 porte sur l'analyse de l'impact de la discrétisation spatiale sur les trajectoires de fissuration, mettant en évidence les artefacts numériques associés aux différents types de maillage. Enfin, la Section 4 conclut cette étude avec des recommandations pratiques pour réaliser des simulations de propagation de fissure par champ de phase en limitant les biais numériques.
## 2 Approche variationnelle pour la rupture par champ de phase

Partant de l'idée originale de Griffith [4], on rappelle l'approche variationnelle pour la rupture proposée par Francfort et Marigo [1]. Cette section couvre les bases théoriques de ce travail, en insistant sur le lien entre la théorie de Griffith et les modèles variationnels par champ de phase pour la rupture. Bourdin et al. [5] propose une discussion plus détaillée sur ce sujet.

### 2.1 Théorie de Griffith

Pour un problème élastique statique, le champ de déplacement d'un corps élastique soumis à des charges externes minimise son énergie potentielle $P = \mathcal{E} - W_{\text{ext}}$. Griffith [4] a introduit l'énergie de surface $G_c$, définie comme l'énergie nécessaire pour créer de nouvelles surfaces libres. Ainsi, l'état d'équilibre d'un solide contenant une fissure doit minimiser l'énergie $P + G_c a$. La dérivation de cette énergie par rapport à la longueur de la fissure $a$ conduit à la définition du taux de restitution d'énergie

$$G = -\frac{\partial P}{\partial a} \quad \text{à charge constante.} \qquad (1)$$

Pour assurer que la propagation soit stable, il faut imposer $G \le G_c$. Les conditions de propagation quasi-statique de la fissure sont :

$$(G - G_c)\dot{a} = 0, \quad G \le G_c, \quad \dot{a} \ge 0. \qquad (2)$$

Ainsi, si $G < G_c$, la fissure ne se propage pas, et si $G = G_c$, la propagation devient possible.

Différents critères existent pour déterminer la direction $\varphi$ dans laquelle la fissure se propage. Parmi ces derniers, on considère le critère de maximisation du taux de restitution d'énergie [6]

$$\varphi = \arg \max_{\varphi^*} G^*(\varphi^*) \qquad (3)$$

où $G^*$ est le taux de restitution d'énergie pour une propagation dans la direction $\varphi$. Il peut être obtenu via la formule d'Amestoy-Leblond [7]. Ce critère est une extension directe du critère de Griffith.

### 2.2 Approche variationnelle de la rupture

Francfort et Marigo [1] ont formalisé la théorie de Griffith dans un cadre variationnel. On considère de nouveau un corps élastique $\Omega$ avec une fissure initiale $\Gamma_0$. Des forces volumiques $f$ sont appliquées dans le volume. Le bord du domaine se découpe en deux sous-parties : $\partial\Omega_u$ et $\partial\Omega_t$, où le déplacement $u$ et les forces de surface sont respectivement imposés. Les surfaces fissurées sont libres. Ainsi, on introduit l'espace des déplacements admissibles

$$\mathcal{U}_{\Gamma} = \{u \in H^1(\Omega \setminus \Gamma) \mid \forall x \in \partial\Omega_u, u(x) = u_{\text{imp}}(x)\}. \qquad (4)$$

Dans ce cadre, l'énergie élastique et le travail externe sont

$$\begin{align}
\mathcal{E}(u, \Gamma) &= \int_{\Omega \setminus \Gamma} \frac{1}{2} \varepsilon(u) : \mathbf{E} : \varepsilon(u) \,dx, \nonumber \\
\mathcal{W}_{\text{ext}}(u) &= \int_{\Omega} \mathbf{f} \cdot u \,dx + \int_{\partial\Omega_t} t \cdot u \,dx. \tag{5}
\end{align}$$
Francfort et Marigo [1] ont proposé de réécrire la minimisation d'énergie proposée par Griffith comme le problème de minimisation

$$ (\mathbf{u}, \Gamma) = \arg \min_{\substack{\mathbf{u}' \in \mathcal{U}_{\Gamma} \\ \Gamma' \supseteq \Gamma_0}} \mathcal{E}(\mathbf{u}', \Gamma') - W_{\text{ext}}(\mathbf{u}') + \int_{\Gamma'} G_c \, d\mathbf{x}. \quad (6) $$

La contrainte $\Gamma' \supseteq \Gamma_0$ impose la croissance de la fissure.

## 2.3 Régularisation via un modèle variationnel de champ de phase

Le problème de minimisation (6) est un problème de discontinuité libre. Francfort et Marigo [1] ont proposé de régulariser ce problème en remplaçant la surface de discontinuité $\Gamma$ par un champ de phase $\alpha$. Ainsi, on considère que le matériau au point $\mathbf{x}$ est sain si $\alpha(\mathbf{x}) = 0$ et qu'il est fissuré si $\alpha(\mathbf{x}) = 1$. Le problème de minimisation devient alors

$$ (\mathbf{u}, \alpha) = \arg \min_{\substack{\mathbf{u} \in \mathcal{U} \\ \alpha \in \mathcal{A}}} \mathcal{E}(\mathbf{u}, \alpha) + \mathcal{D}(\alpha) - W_{\text{ext}}(\mathbf{u}) \quad (7) $$

où $\mathcal{U}$ et $\mathcal{A}$ sont les espaces d'admissibilité du champ de déplacement $\mathbf{u}$ et du champ de phase $\alpha$. L'irréversibilité de la fissuration est imposée par une contrainte de croissance locale du champ de phase. Le terme $\mathcal{D}(\alpha)$ correspond à la version régularisée de l'intégrale de dissipation sur la surface de fissure. Ainsi l'énergie élastique et l'énergie dissipée sont réécrites

$$ \begin{aligned} \mathcal{E}(\mathbf{u}, \alpha) &= \int_{\Omega} \frac{1}{2} a(\alpha) \varepsilon(\mathbf{u}) : \mathbf{E} : \varepsilon(\mathbf{u}) \, d\mathbf{x}, \\ \mathcal{D}(\alpha) &= \frac{G_c}{c_w} \int_{\Omega} \frac{w(\alpha)}{\ell} + \ell \nabla \alpha \cdot \nabla \alpha \, d\mathbf{x}. \end{aligned} \quad (8) $$

À noter que la régularisation introduit une longueur $\ell$, ainsi que les fonctions $a(\alpha)$ et $w(\alpha)$. Pour les fonctions $a$ et $w$, on retient le modèle AT2

$$ a(\alpha) = (1 - \alpha)^2, \quad w(\alpha) = \alpha^2, \quad c_w = 2. \quad (9) $$

De plus, la solution régularisée converge vers la solution de Griffith lorsque $\ell \to 0$ [8]. Ainsi, la longueur de régularisation $\ell$ doit être choisie suffisamment faible pour garantir la convergence vers la mécanique de la rupture.

## 3 Influence du maillage sur la trajectoire de fissure

Cette section porte sur l'analyse de l'impact de la discrétisation spatiale sur les trajectoires de fissuration dans les simulations par éléments finis utilisant l'approche champ de phase. Dans un premier temps, le problème mécanique de référence qui servira de base comparative pour notre étude est présenté. Ce problème est ensuite résolu selon deux approches distinctes : d'une part, des simulations de référence fondées sur la mécanique linéaire élastique de la rupture, et d'autre part, des simulations par éléments finis utilisant la méthode champ de phase. Pour les simulations par champs de phase, différentes discrétisations spatiales sont employées.
<figure><img src="image_2.png" /><figcaption>FIGURE 1 – Géométrie et chargement du problème de base.</figcaption></figure>

<table><caption>TABLE 1 – Résumé des paramètres adimensionnés du problème et de leurs valeurs.</caption><thead><tr><th rowspan="2"></th><th colspan="4">Géométrie</th><th colspan="3">Matériau</th><th colspan="3">Chargement</th></tr><tr><th>H</th><th>W</th><th>a<sub>0</sub></th><th>e</th><th>E</th><th>ν</th><th>G<sub>c</sub></th><th>u<sub>min</sub></th><th>Δu</th><th>δ</th></tr></thead><tbody><tr><th>Symbole</th><td>H</td><td>W</td><td>a<sub>0</sub></td><td>e</td><td>E</td><td>ν</td><td>G<sub>c</sub></td><td>u<sub>min</sub></td><td>Δu</td><td>δ</td></tr><tr><th>Valeur</th><td>1</td><td>10H</td><td>H</td><td>H/4</td><td>1</td><td>0.3</td><td>1</td><td>0.5</td><td>1.5</td><td>1</td></tr></tbody></table>

## 3.1 Problème

Dans ce problème, illustré sur la figure 1, on considère une bande de hauteur adimensionnée $H = 1$ et de longueur $W = 10H$. On introduit dans cette bande une fissure initiale excentrée d'une distance $e = H/4$ du centre et de longueur $a_0 = H$. L'idée derrière ce problème est d'imposer un chargement de traction sur la bande afin de propager la fissure. On cherche la trajectoire de la fissure verifiant la théorie de Griffith (équations (2) et (3)). Afin de garantir une propagation stable de la fissure, cette bande est soumise à un échelon spatial (régularisé) de traction qui se propage progressivement le long de l'éprouvette. Le champ de déplacement imposé s'écrit :

$$ \boldsymbol{u}_{\text{imp}} = \pm \left( \Delta u S \left( -\frac{x-\lambda}{\delta} \right) + u_{\text{min}} \right) e_y, \quad S(x) = \frac{1}{1+\exp(-x)}, \qquad (10) $$

où $S(x)$ est une fonction sigmoïde, $\Delta u$ est le saut de déplacement, $u_{\min}$ est le déplacement imposé avant le passage de l'échelon, $\lambda$ est le paramètre de chargement qui correspond à l'avancée de l'échelon et $\delta$ est la longueur caractéristique de la zone de transition. Les paramètres utilisés pour cette étude sont détaillés dans le tableau 1.

La solution de référence pour ce problème est déterminée en utilisant la mécanique de la rupture classique, comme introduite dans la sous-section 2.1. Plus précisément, la solution numérique est obtenue en réalisant des simulations élastiques par éléments finis où la fissure est explicitement intégrée dans le maillage. Les facteurs d'intensité des contraintes sont calculés via la méthode de l'intégrale d'interaction [9], utilisant une formulation en intégrale de domaine. La résolution de l'équation (3) permet ensuite de déterminer la direction potentielle de propagation de la fissure. Lorsque le critère de Griffith est atteint dans cette direction, la fissure est propagée et le domaine est remaillé pour prendre en compte cette évolution.

De plus, le problème est également résolu en utilisant l'approche variationnelle par champ de phase. Cette résolution s'appuie sur un algorithme de minimisation alternée [2]. Les sous-problèmes portant sur le champ de déplacement et le champ de phase sont résolus par éléments finis via la librairie FEniCSx [10].
<figure><img src="image_3.png" /><figcaption>FIGURE 2 – Illustration des maillages utilisés pour les simulations par champ de phase. Pour clarifier les illustrations, ces maillages sont volontairement grossier ($\Delta x = H/16$).</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>FIGURE 3 – Trajectoire de fissure prédite par champ de phase avec un maillage structuré.</figcaption></figure>

## 3.2 Maillages

Afin d'étudier l'influence du maillage sur la trajectoire de fissures dans les simulations de la fissuration par champ de phase, deux types de maillage sont utilisés : structuré et non structuré. Ces deux types de maillage sont illustrés respectivement dans la figure 2a et la figure 2b. En plus de ces différents maillages, les simulations sont répétées pour différents tailles de maillages $\Delta x \in \{H/64, H/96, H/128\}$.

## 3.3 Comparaison des trajectoires

Dans un premier temps, nous analysons l'impact d'un maillage structuré sur la trajectoire de fissure pour différentes tailles de maillage $\Delta x$. La figure 3 compare ces résultats à la solution de référence basée sur la mécanique de la rupture. Contrairement à la référence, les simulations par champ de phase ne reproduisent pas la trajectoire exponentielle attendue. Comme les simulations par champ de phase devraient converger vers la solution de la mécanique de la rupture, nous concluons que ce biais est d'origine numérique. Plus précisément, le maillage structuré introduit une anisotropie artificielle [11, 12]
<figure><img src="image_5.png" /><figcaption>FIGURE 4 – Trajectoire de fissure prédite par champ de phase avec un maillage non structuré.</figcaption></figure>

qui contraint la fissure à suivre une trajectoire rectiligne alignée avec le maillage.

Nous examinons ensuite l'impact d'un maillage non structuré sur la trajectoire de fissure pour différentes tailles de maillage $\Delta x$. La figure 4 présente les trajectoires obtenues par simulation par champ de phase sur des maillages non structurés. Bien que la trajectoire de fissure ne soit pas parfaitement reproduite, les résultats du champ de phase suivent globalement la trajectoire de référence. Le caractère désordonné du maillage réduit l'effet d'anisotropie induite par le maillage, car les orientations variables des éléments lissent globalement cet effet, permettant ainsi une meilleure prédiction de la trajectoire.

# 4 Conclusion

Ce travail propose une analyse complète des biais numériques induit par le maillage dans les simulations de fissuration par champ de phase, en comparaison avec les résultats de la mécanique de la rupture. Notre étude antérieure [3] avait établi l'importance cruciale de l'initialisation des fissures, démontrant que la largeur initiale et le choix de conditions aux limites influencent significativement les résultats. La présente recherche complète cette analyse en examinant systématiquement l'impact de la discrétisation spatiale sur les trajectoires de propagation. Les résultats obtenus confirment que le choix du maillage joue un rôle déterminant dans la précision des simulations. En particulier, ils montrent que les maillages structurés introduisent des artefacts numériques significatifs, conduisant à des prédictions de trajectoire erronées. À l'inverse, les maillages non structurés permettent une modélisation plus fidèle des phénomènes physiques attendus. Ensemble, ces travaux fournissent un cadre méthodologique pour limiter les biais dans les simulations de fissuration par champ de phase, les rendant ainsi plus fiables et précises. Ces résultats ouvrent des perspectives pour l'emploi de l'approche variationnelle par champ de phase dans les applications industrielles.
Références

[1] Francfort, G. A., & Marigo, J.-J. (1998). Revisiting brittle fracture as an energy minimization problem. *Journal of the Mechanics and Physics of Solids*, 46(8), 1319–1342. https://doi.org/10.1016/S0022-5096(98)00034-9

[2] 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. https://doi.org/10.1016/S0022-5096(99)00028-9

[3] Loiseau, F., & Lazarus, V. (in press). How to introduce an initial crack in phase field simulations to accurately predict the linear elastic fracture propagation threshold? *Journal of Theoretical, Computational and Applied Mechanics*. https://doi.org/10.48550/arXiv.2502.03900

[4] Griffith, A. A. (1920). The phenomena of rupture and flow in solids. *Philosophical Transactions of the Royal Society of London. Series A, Containing Papers of a Mathematical or Physical Character*, 221(582–593), 163–198. https://doi.org/10.1098/rsta.1921.0006

[5] Bourdin, B., Francfort, G. A., & Marigo, J.-J. (2008). The Variational Approach to Fracture. *Journal of Elasticity*, 91(1), 5–148. https://doi.org/10.1007/s10659-007-9107-3

[6] Hussain, M. A., Pu, S. L., & Underwood, J. (1973). Strain Energy Release Rate for a Crack Under Combined Mode I and Mode II. *1973 National Symposium on Fracture Mechanics, Part II*, 2, 1–27. https://doi.org/10.1520/STP33130S

[7] Amestoy, M., & Leblond, J. B. (1992). Crack paths in plane situations – II. Detailed form of the expansion of the stress intensity factors. *International Journal of Solids and Structures*, 29(4), 465–501. https://doi.org/10.1016/0020-7683(92)90210-K

[8] Chambolle, A. (2004). An approximation result for special functions with bounded deformation. *Journal de Mathématiques Pures et Appliquées*, 83(7), 929–954. https://doi.org/10.1016/j.matpur.2004.02.004

[9] Yu, H., & Kuna, M. (2021). Interaction integral method for computation of crack parameters K–T – A review. *Engineering Fracture Mechanics*, 249, 107722. https://doi.org/10.1016/j.engfracmech.2021.107722

[10] Baratta, I. A., Dean, J. P., Dokken, J. S., Habera, M., Hale, J. S., Richardson, C. N., Rognes, M. E., Scroggs, M. W., Sime, N., & Wells, G. N. (2023). DOLFINx: The next generation FEniCS problem solving environment. *Zenodo*. https://doi.org/10.5281/zenodo.10447666

[11] Negri, M. (1999). The anisotropy introduced by the mesh in the finite element approximation of the mumford-shah functional. *Numerical Functional Analysis and Optimization*, 20(9–10), 957–982. https://doi.org/10.1080/01630569908816934

[12] Negri, M. (2003). A finite element approximation of the Griffith’s model in fracture mechanics. *Numerische Mathematik*, 95(4), 653–687. https://doi.org/10.1007/s00211-003-0456-y