Modélisation de la propagation des fissures avec Z-set

Introduction

L'étude de la fissuration des composants industriels joue un rôle fondamental dans la garantie de leur sécurité et de leur fiabilité, en particulier dans des domaines sensibles comme l'aéronautique^1–4. Dans ce secteur, où les contraintes mécaniques et les exigences de performance sont extrêmes, la modélisation précise de l’apparition et de l’évolution des défauts est essentielle afin de prévenir les défaillances catastrophiques.

Pour relever ce défi, divers cadres numériques et théoriques ont été mis au point. La méthode dite des éléments finis étendus (XFEM) permet de représenter, grâce à des fonctions d’enrichissement localisées, des discontinuités de fissures sans recourir au remaillage⁵ et s’avère efficace pour capturer les singularités locales et la bifurcation de fissures complexes. Parallèlement, les critères énergétiques fondés sur les intégrales J et G, issus du formalisme de la mécanique linéaire de la rupture (LEFM), quantifient respectivement la libération d’énergie et la ténacité du matériau dans des conditions élastiques quasi-statiques⁶. Ces approches s’appuient principalement sur des théories de la mécanique linéaire de la rupture, qui conviennent bien pour les matériaux fragiles ou présentant de petits défauts. Cependant, lorsqu’il s’agit de matériaux élasto-plastiques, les limites de ces modèles apparaissent rapidement.

En effet, la plasticité induit une dissipation d'énergie significative et une évolution des contraintes résiduelles, compliquant ainsi l'analyse et la prédiction de la propagation des fissures. L'un des défis majeurs réside dans la prédiction précise de la direction de propagation, particulièrement dans des environnements dissipatifs.

Pour pallier ces problèmes, la méthode G-theta propose une solution efficace. Cette méthode combine le calcul d’une grandeur d’énergie caractéristique de la fissuration en milieu dissipatif avec une prédiction efficace du chemin de propagation des fissures.

Dans cette optique, l'outil numérique Z-cracks, développé en collaboration entre l’ONERA⁷ et TRANSVALOR et implémenté dans la suite de simulation par éléments finis Z-set, offre un cadre générique et performant pour l'étude des fissures en 3D. Ce module permet de calculer les facteurs d'intensité des contraintes (FIC) pour des configurations de fissures statiques et en mode mixte, ainsi que de simuler la propagation des fissures, tout en intégrant des techniques robustes de remaillage adaptatif 3D et un post-traitement multi-threadé efficace pour l’extraction des FIC.

Dans cet article, nous explorons l’efficacité de la méthode G-theta intégrée dans Z-cracks. En appliquant cette méthode à divers cas de fissurations, nous analysons sa convergence, notamment dans un contexte élasto-plastique.

Modélisation de la propagation de fissures par la méthode G-theta

Afin de résoudre numériquement un problème de fissuration, l’algorithme de Z-cracks comprend plusieurs étapes. L’étape la plus délicate consiste à générer un maillage conforme autour de la fissure pendant sa propagation, ce qui constitue une tâche complexe en raison de la singularité des champs de contraintes et de déformations au niveau de la pointe de fissure. Tout d’abord, une fissure ayant une forme géométrique prédéfinie est insérée dans un maillage sain préexistant. L’intersection entre ce dernier et le maillage représentant la fissure est générée. Le maillage de cette structure intermédiaire est ensuite adapté pour obtenir une discrétisation de qualité satisfaisante, essentielle pour pouvoir calculer les facteurs d'intensité des contraintes (FIC) avec précision au niveau du front de fissure. Lors de cette étape de remaillage, Z-cracks fait appel à des outils externes de remaillage automatique adaptatif (e.g. le mailleur MeshGems), offrant une flexibilité de calcul par élément finis notable.

La prévision de la propagation des fissures repose sur l’évaluation précise des grandeurs caractéristiques de la fissuration, en particulier les facteurs d’intensité des contraintes et le taux de restitution d’énergie G. Une analyse thermomécanique par éléments finis permet d’estimer l’état de contraintes dans la structure sous un chargement donné. À partir de ce champ de contraintes, les FIC sont calculés en plusieurs points du front de fissure afin de déterminer l’angle de propagation, l’incrément d’avancement, ainsi que la nouvelle position de chaque nœud du front.

La progression du front est ensuite régie par des lois de propagation (e.g., Paris-Elber^8,9), permettant d’estimer l’évolution géométrique de la fissure dans le matériau. Cette évolution implique une mise à jour continue de la topologie du modèle d’éléments finis.

Le module Z-cracks automatise ce processus en extrayant les grandeurs de fissuration, et l’angle de bifurcation à l’aide de la méthode G-theta. Cette méthode proposée par Destuynder et al.¹⁰, repose sur des intégrales de domaine et sur la dérivation Lagrangienne de l'énergie potentielle, par rapport à un champ de vitesse d'extension virtuelle de la fissure, que nous noterons ici θ. Nous considérons un domaine de référence Ω contenant une fissure ayant un front Γ ⁰. Soit M un point du front et S l’abscisse curviligne le long du front. La transformation F ^η de tous les points M du domaine Ω vers un domaine Ω^η, due à la propagation de la fissure est donnée par :

où η est une constante (η = 0 correspond à la configuration de référence) et θ(M) est le champ de déplacement virtuel de la fissure.

Figure 1:  Présentation d’une fissure dans un domaine volumique Ω. 

• Dans le cas d’un comportement thermo-élastique du matériau, il est possible de calculer une intégrale invariante par domaine. Cette approche permet d’évaluer l’énergie potentielle qui serait libérée lors d’un déplacement infinitésimal du front de fissure, en supposant une évolution purement élastique du matériau soumis aux contraintes résiduelles induites par la thermoélasticité. En l’absence de forces volumiques, le taux global de restitution d’énergie associé à cette extension virtuelle de la fissure s’exprime par :

où ε^th désigne la déformation thermique, et ε est la déformation totale.

• De manière analogue, un invariant de domaine peut être évalué dans le cas d’un comportement élasto-plastique du matériau. L’énergie potentielle libérée lors d’un déplacement infinitésimal du front de fissure est alors estimée en supposant une évolution purement élastique du matériau, soumis à un champ de contraintes résiduelles fixé, résultant des déformations thermiques et viscoplastiques. La grandeur intégrée, notée G_p, est définie par l’équation suivante :

où ε^ae désigne la déformation inélastique.

En intégrant le long du front de la fissure Γ ⁰, le front étant discrétisé à l’aide de n_p points de contrôle associés à des fonctions de forme N_j(s), le taux de restitution d’énergie peut être déterminé par la résolution du système d’équations suivant à n_p inconnues :

Validation de la formulation G-theta

La validation de la formulation G-theta est cruciale afin de s'assurer que la modélisation et l'évaluation des fissures dans un contexte élasto-plastique sont fiables et cohérentes avec la physique des matériaux. Pour ce faire, les valeurs du taux de restitution d’énergie G ont été comparées à celles obtenues pour un matériau purement élastique. Cette comparaison permet de quantifier la contribution des contraintes résiduelles inélastiques (thermiques et viscoplastiques) sur la libération d'énergie dans un système réaliste, où des déformations irréversibles jouent un rôle significatif.

Afin de renforcer la crédibilité de l’approche, une phase de validation a été menée à l’aide de solutions de référence analytiques et semi-analytiques en 2D et 3D. Ces validations ont permis de confirmer la capacité de la méthode à reproduire les mécanismes physiques attendus dans des configurations contrôlées. Cette phase a ensuite été complétée par une analyse de la convergence numérique, en évaluant l’effet du maillage et des paramètres de résolution sur la stabilité et la précision des résultats.

Une étude de sensibilité sur le calcul de la valeur de G (cf. équation (2)) a été menée sur une éprouvette bidimensionnelle entaillée soumise à une traction uniaxiale (cf. figure 2(a)), en fonction des rayons d'intégration intérieur r_b et extérieur r_a​. Ces deux rayons délimitent la zone où G est calculé dans les cas élastique et élasto-plastique. L’analyse consiste à fixer l’un des deux rayons et à observer la convergence de G en fonction de l’autre. En fixant r_b = 0 et en faisant varier r_a​, les résultats dans le cas purement élastique montrent que G converge rapidement en appliquant la règle usuelle consistant à choisir r_a​ égal à trois fois la taille minimale des éléments du maillage (cf. figure 2(b)). En revanche, dans le cas élasto-plastique, la convergence de G n’est atteinte que pour des valeurs de r_a​​ supérieures à 0,2 mm. Cela suggère que le rayon externe r_a​​ doit être choisi au minimum trois fois supérieur à la taille estimée de la zone plastique pour assurer la fiabilité du calcul.

gauche : (a) 

droite : (b)

Figure 2: (a) Éprouvette entaillée à un seul bord, soumise à une charge de traction de 50 MPa et le maillage correspondant. (b) Comparaison de l’évolution de valeur de G pour le cas élastique et élasto-plastique en fonction du rayon d’intégration externe r_a. 

En parallèle, quand le rayon externe r_a​​ est fixé, les résultats montrent que le rayon interne r_b n’a pas d’influence significative sur G dès lors que la distance entre les deux rayons​ est suffisante. En revanche, lorsque r_a est inférieure à la taille de la zone plastique, les variations de r_b affectent le calcul de G.

Cette analyse a permis d'affiner le choix optimal de ces rayons, en particulier pour le cas élasto-plastique où la zone plastique s'étend depuis le front de fissure au-delà de trois éléments de maillage, taille auparavant utilisée pour définir r_a dans le cas élastique.

L'influence du gradient proche du front de fissure sur les valeurs de G a également été examinée. Le gradient de déformation viscoplastique varie sur une région composée de quelques éléments à partir du front de fissure, modifiant ainsi les valeurs de G en fonction du paramètre h_min sélectionné, correspondant à la distance minimale par rapport au front de fissure à partir de laquelle le gradient de la déformation inélastique est pris en compte. Le paramètre h_min permet de négliger l’influence de ce gradient dans les zones très proches du front, où il peut présenter des fluctuations importantes susceptibles de perturber l’évaluation du champ mécanique. Ces études de sensibilité ont également été réalisées sur un maillage non structuré en 2D et pour les grandes déformations, afin d'évaluer la convergence des valeurs de G dans les deux cas élastique et élasto-viscoplastique.

Cette première validation a ensuite été étendue à des cas tridimensionnels sur une géométrie plane de dimensions 1 mm × 1 mm × 0,5 mm, soumise à une contrainte σ₂₂ = 110 MPa. Au sein de cette géométrie présentant un plan de symétrie z = 0, nous avons inséré une fissure en bord comme l’illustre la Figure 3 (a) ci-dessous.

Le champ de déformation plastique cumulée (epcum) a ensuite été calculé au niveau de la surface de la géométrie (z = 0.5 mm) ainsi qu’à cœur (z = 0). Les résultats montrent un champ légèrement différent aux deux endroits (cf. Figure 3 (b)). Cette différence observée dans la zone plastique, entre le cœur et la surface, peut être expliquée par les variations des champs de contrainte et de déformation à travers l’épaisseur.

Nous avons ensuite analysé la convergence de G, en surface et à cœur, en faisant varier le rayon externe r_a​ tout en maintenant le rayon interne r_b = 0 constant. Les résultats de cette dernière analyse, présentés dans la Figure 3(c), montrent que la différence de G entre le cœur et la surface est négligeable, tant dans le cas élastique que plastique, et peut donc être ignorée.

(a) - (b)

(c)

Figure 3: (a) La géométrie considérée, dans laquelle une fissure en bord traversante dans l’épaisseur est intégrée, et le maillage utilisé. (b) Le champ de déformation plastique cumulée (epcum) montrant la zone plastique en surface et au cœur de la géométrie. (c) Résultats du calcul de la convergence de G en fonction du rayon extérieur r_a, en surface et au cœur, pour des matériaux élastique et plastique.

Toutefois, pour une fissure traversante dans l’épaisseur, la convergence de G diffère entre la surface et le cœur de la plaque dans le cadre plastique, tel qu’illustré dans la Figure 4 ci-dessous. En effet, celle-ci est plus rapide à cœur qu’en surface, en raison des variations de la zone plastique à travers l’épaisseur. En revanche, dans le cadre élastique, G converge à la même vitesse que ce soit à cœur ou en surface, bien que les valeurs limites atteintes diffèrent, ce qui n’est pas observé en comportement élasto-plastique.

Figure 4: Résultats du calcul de la convergence de G en surface et au cœur, en fonction du rayon extérieur r_a (r_b = 0), dans le cas d’une fissure traversante dans l’épaisseur, pour des comportements plastique et élastique. 

Dans une dernière étape, nous avons analysé l’effet de l’interaction entre deux fronts de fissure sur la valeur de G. L’étude porte sur une plaque de dimensions 2 mm × 2 mm × 0,2 mm, soumise à une traction uniforme σ₂₂ = 150 MPa appliquée sur ses faces supérieure et inférieure. Une fissure traversante de longueur 2a = 0,5 mm y est introduite (cf. Figure 5).

Les calculs de la convergence de G au cœur de la géométrie ont été réalisés en fonction du rayon externe r_a pour les comportements élastique et plastique. Afin de réduire les temps de calcul, des conditions de symétrie sont appliquées dans deux configurations. Dans le premier cas, le modèle est réduit à la moitié de la géométrie en exploitant la symétrie par rapport au plan médian z = 0, sur lequel la condition de déplacement U₃ = 0 est imposée. Dans le second cas, la réduction est étendue à un quart du domaine en introduisant une seconde symétrie selon le plan x = 0, où l’on impose U₁ = 0, en plus de U₃ = 0 au plan z = 0.

Figure 5:  Résultats du calcul de la convergence de G au cœur (z = 0), en fonction du rayon extérieur r_a (r_b = 0, h_min = 0), dans le cas de l’interaction entre deux fronts de fissure, pour des comportements (a) élastique et (b) plastique. 

Les résultats de l’analyse montrent que la convergence de G pour le premier front est influencée par la présence du second. En effet, pour le premier cas, les valeurs de G diminuent significativement lorsque le rayon d'intégration extérieur s'approche du second front de fissure (r_a = 0,5 mm), perturbant ainsi la convergence des résultats (cf. courbes bleues de la Figure 5).

Dans le second cas, la proximité de la frontière géométrique (r_a = 0,25 mm) affecte les résultats (cf. courbes oranges de la Figure 5), entraînant une diminution des valeurs de G après franchissement de cette limite. Bien que les deux configurations soient mécaniquement équivalentes, la zone d’intégration réduite au-delà de r_a = 0,25 mm explique cette divergence.

Conclusion

Cette étude a confirmé la pertinence de la méthode G-theta pour prédire la propagation des fissures dans des matériaux à comportement élasto-plastique. Intégrée dans le module Z-cracks de la suite Z-set, cette approche permet de calculer l’énergie libérée lors de l’avancement des fissures, même dans des contextes industriels complexes.

Grâce à son remaillage adaptatif 3D et à son post-traitement optimisé, Z-cracks offre une solution performante pour l’analyse tridimensionnelle des fissures. Z-set s’impose ainsi comme une plateforme puissante, prête à relever les défis actuels et futurs en mécanique de la rupture.