# Traitement numérique tridimensionnel de la rupture ductile progressive par une méthode XFEM cohésive

Konstantinos Nikolakopoulos, Antonio Kaniadakis, Jean-Philippe Crété, Patrice Longère

► To cite this version:

Konstantinos Nikolakopoulos, Antonio Kaniadakis, Jean-Philippe Crété, Patrice Longère. Traitement numérique tridimensionnel de la rupture ductile progressive par une méthode XFEM cohésive. 15e colloque national en calcul des structures, Université Polytechnique Hauts-de-France [UPHF], May 2022, Hyères-les-Palmiers, France. hal-03717881

HAL Id: hal-03717881

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

Submitted on 8 Jul 2022

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

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

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
# Traitement numérique tridimensionnelle de la rupture ductile progressive par une méthode XFEM cohésive

K. Nikolakopoulos<fnref n="1" />,<fnref n="2" />, A. Kaniadakis<fnref n="1" />,<fnref n="2" />, J.P. Crété<fnref n="2" />, P. Longère<fnref n="1" />

**Résumé** — Ce travail porte sur le traitement numérique de la rupture ductile dans les matériaux méta-talliques. Les trois étapes successives de micro-cavitation, de localisation et d'ouverture de fissure sont reproduites dans une méthodologie unifiée tridimensionnelle qui combine le modèle de plasticité micro-poreuse GTN, la XFEM et un modèle de zone cohésive. Il est montré que la méthodologie est capable de reproduire numériquement la réponse globale de manière indépendante du maillage et de conduire à des surfaces de rupture réalistes.

**Mots clés** — rupture ductile, XFEM, GTN, CZM.

## 1 Introduction

Les métaux et les alliages sont les matériaux les plus utilisés dans l'industrie du transport (aéronautique, navale et automobile) du fait de leur haute résistance et de leur bonne ductilité. Dans le cadre de la prédiction de leur tenue à la mise en forme ou/et à des surcharges accidentelles, il est donc crucial d'être capable de reproduire qualitativement et quantitativement les étapes successives menant potentiellement à leur rupture (Figure 1). Bien qu'il y ait un consensus dans la communauté scientifique concernant les différentes étapes qui conduisent à la rupture ductile (endommagement par micro-cavitation, localisation de la déformation qui donne lieu à une méso-fissure, formation et propagation de la macro-fissure), le traitement numérique de ces phénomènes reste sujet à débat, surtout si on adopte la MEF. Cette méthode présente en effet une dépendance pathologique au maillage dès le régime adoucissant (perte d'unicité de la solution) qui va conduire à une bande de localisation. Il existe dans la littérature des techniques de régularisation pour la plupart basées sur des approches non locales mais elles nécessitent un maillage très fin et donc un coût de calculs élevé qui s'avèrent rédhibitoires lorsqu'on s'intéresse à des structures de grandes dimensions. On présente ici une méthodologie tridimensionnelle dans laquelle la bande de localisation/la fissure sont incorporées dans l'élément fini dans une formulation basée sur la XFEM. En l'occurrence, la bande de localisation est traitée comme une fissure cohésive (XFEM cohésif) dont la perte progressive de cohésion conduit à la fissure proprement dite (XFEM standard).

<figure><img src="image_2.png" /><figcaption>FIGURE 1 – Etapes de la rupture ductile et méthodologie unifiée développée</figcaption></figure>

<footnote n="1">¹Institut Clément Ader, Université de Toulouse, ISAE-SUPAERO, INSA-T, UT3/UPS, IMT Mines Albi, CNRS, 3 rue Caroline Aigle 31400, Toulouse, France, {konstantinos.nikolakopoulos,antonio.kaniadakis,patrice.longere}@isae-supaero.fr</footnote>

<footnote n="2">²Laboratoire Quartz, ISAE-SUPMECA, 3 rue Fernand Hainaut 93407, Saint-Ouen, France, jean-philippe.crete@isae-supmeca.fr</footnote>
# 2 Traitement numerique

## 2.1 Modèle de plasticité microporeuse GTN

On suppose que le matériau considéré est indépendant de la vitesse de déformation et de la température. Le couplage entre plasticité et endommagement ductile est décrit en utilisant le modèle de GTN [7], où le potentiel plastique s'écrit :

$$ \Phi = \left(\frac{\sigma_{eq}}{\sigma_y}\right)^2 + 2q_1 f \cosh\left(-\frac{3}{2}q_2 \frac{p_m}{\sigma_y}\right) - (q_1 f)^2 - 1 = 0 \quad (1) $$

avec $\sigma_{eq}$ la contrainte équivalente de Mises, $p_m = -\sigma_m$ la pression hydrostatique, $\sigma_y$ la contrainte d'écoulement, $f$ la fraction volumique des cavités et $q_1, q_2$ des constantes matériaux. L'écrouissage isotrope est décrit via une loi de type Voce :

$$ \sigma_y = R_0 + R_\infty [1 - \exp(-k\kappa)] \quad (2) $$

avec $R_0$ la taille initiale du domaine d'élasticité, $R_\infty$ la valeur saturante de l'écrouissage, $k$ une constante liée à l'écrouissage et $\kappa$ la déformation plastique cumulée. Le taux de porosité $\dot{f}$ est décomposée en deux parties, $\dot{f}_g$ pour la croissance des vides existants et $\dot{f}_n$ pour la germination des nouveaux vides :

$$ \dot{f} = \dot{f}_g + \dot{f}_n \quad (3) $$

Pour la croissance on a :

$$ \dot{f}_g = (1 - f) \dot{\varepsilon}_v \quad (4) $$

avec $\dot{\varepsilon}_v$ la vitesse de déformation plastique volumique. Pour la germination on a :

$$ \dot{f}_n = A_n \dot{\kappa} \quad (5) $$

où $A_n$ s'exprime sous la forme d'une fonction gaussienne :

$$ A_n = \frac{f_N}{s_N \sqrt{2\pi}} \exp\left(-\frac{1}{2}\left(\frac{\kappa - \kappa_N}{s_N}\right)^2\right) \quad (6) $$

avec $f_N$ la fraction volumique des vides qui peut être finalement germée, $\kappa_N$ la déformation plastique à la vitesse maximale de la germination des micro-vides et $s_N$ l'écart type.

## 2.2 Étape de localisation (XFEM+CZM)

La première condition pour qu'un élément fini soit un candidat potentiel pour la localisation est basée sur le postulat de stabilité de Drucker, tel que :

$$ \det(\mathbf{L}) = 0 \quad (7) $$

avec $\mathbf{L}$ l'opérateur tangent du matériau. Si cette condition est respectée, la localisation est traitée ici comme un phénomène résultant :

1. soit d'une instabilité plastique.

2. soit de la coalescence des microcavités.

De ce fait, deux critères sont considérés dans ces travaux. Pour le cas 1, la localisation induite par instabilité plastique intervient lorsque le tenseur acoustique $\mathbf{Q}$ devient singulier. Le critère revient alors à rechercher le plan pour lequel la normale vérifie $\det(\mathbf{Q}) < 0$ pour la première fois, tout en minimisant le tenseur acoustique :

$$ \min(\det(\mathbf{Q})) < 0 \quad (8) $$

Pour le cas 2, la localisation induite par la coalescence des micro-cavité requiert en fait deux critères : un pour l'amorçage de la bande et un autre pour l'orientation de la bande. La condition pour l'amorçage est que la porosité soit égale à une porosité critique :

$$ f = f_c \quad (9) $$
Dans ce travail, le choix de l'orientation de la bande est dérivé de l'état de contrainte et rend compte de la compétition entre Mode I et Mode II. Les plans de localisation candidats dans l'élément sont respectivement celui qui maximise la contrainte de traction et celui qui maximise la contrainte de cisaillement. Un critère est alors appliqué pour discerner les localisations de type Mode I et Mode II en fonction du taux de triaxialité des contraintes. La région dans laquelle les deux modes interagissent est décrite par une loi des mélanges, voir Figure 2.

<figure><img src="image_3.png" /><figcaption>FIGURE 2 – Angle d'orientation de la bande en fonction du taux de triaxialité des contraintes dans le cas d'une localisation induite par la coalescence des vides</figcaption></figure>

## 2.3 Modèle de zone cohésive

Dans le modèle de zone cohésive considéré ici, le vecteur force de traction est exprimé par :

$$ \mathbf{t} = (1 - D) \mathbf{t}_0 \quad (10) $$

avec $t_0$ le vecteur de traction au début de la localisation au sein de l'élément, D une variable d'endommagement qui décrit la perte progressive de cohésion au sein de la bande de localisation. L'évolution de D suit une loi puissance qui est fonction du déplacement relatif équivalent $\Delta$ :

$$ D = \begin{cases} 0, & \Delta < 0 \\ \left(\frac{\Delta}{\Delta_c}\right)^{\gamma}, & 0 < \Delta \le \Delta_c \\ 1, & \Delta > \Delta_c \end{cases} \quad (11) $$

avec $\Delta_c$ la valeur critique du déplacement relatif équivalent à la rupture ultime ($D = 1$). Le déplacement relatif équivalent $\Delta$ s'écrit :

$$ \Delta = \sqrt{\ll \delta_n \gg^2 + \delta_1^2 + \delta_2^2} \quad (12) $$

avec $\delta_n$ la composante normale du déplacement, $\delta_1$ et $\delta_2$ les composantes tangentielles dans le plan de la bande. La matrice tangente de la loi cohésive est calculée comme suit, où l'indice loc se réfère au repère local de la bande de localisation :

$$ C_{\text{loc}}(i, j) = \frac{\partial t^i}{\partial \delta_j} = - \frac{\gamma}{\Delta_c} \left(\frac{\Delta}{\Delta_c}\right)^{\gamma-1} \frac{\delta_j}{\Delta} t_0^i \quad (13) $$

avec $t_0^i$ les composantes du vecteur de traction initial.

## 2.4 Couplage XFEM-CZM

La zone cohésive est insérée dans la formulation XFEM pour créer une discontinuité forte "cohésive" (voir Wolf et al. [8]), ce qui conduit au système d'équations discrètes suivant :

$$ \mathbf{K} \begin{Bmatrix} d\alpha \\ d\beta \end{Bmatrix} = \begin{Bmatrix} \mathbf{F}_{\text{ext}}^{\alpha} \\ \mathbf{F}_{\text{ext}}^{\beta} \end{Bmatrix} - \begin{Bmatrix} \mathbf{F}_{\text{int}}^{\alpha} \\ \mathbf{F}_{\text{int}}^{\beta} \end{Bmatrix} \quad (14) $$

où $d\alpha, d\beta$ représentent les incréments de déplacement nodaux des degrés de liberté standards et enrichis et $\mathbf{F}_{\text{int}}$ les efforts internes avec l'exposant $(\alpha, \beta)$ indiquant la correspondance avec les degrés de liberté standards ou enrichis.
$$ \mathbf{F}_{int}^{\alpha} = \int_{V_e} \mathbf{B}^T \boldsymbol{\sigma} dV^e \qquad (15) $$

$$ \mathbf{F}_{int}^{\beta} = \int_{V_e} \mathbf{B}^{*T} \boldsymbol{\sigma} dV^e + \mathbf{F}_{coh} \qquad (16) $$

$$ \mathbf{B}^* = \mathbf{p}\mathbf{N}^* \qquad (17) $$

$$ N_j^* = (H(x) - H_j) N_j \qquad (18) $$

avec $\boldsymbol{\sigma}$ la contrainte de Cauchy, $\mathbf{B}$ matrice qui contient les dérives spatiales des fonctions d'interpolations, $\mathbf{p}$ l'opérateur différentiel. Pour la partie XFEM on indique avec un astérisque la matrice des dérives spatiales $\mathbf{B}^*$ et la matrice des fonctions d'interpolations $\mathbf{N}^*$. $H(x)$ est la fonction d'Heaviside et $H_j$ la valeur de la fonction calculé sur chaque noeud. Dans ce travail, on considère un contrôle en déplacement afin qu'il n'y ait pas d'efforts extérieurs appliqués aux éléments finis. La matrice de rigidité s'écrit alors :

$$ \mathbf{K} = \left[ \begin{matrix} \int_{V_e} \mathbf{B}^T \mathbf{L} \mathbf{B} dV^e & \int_{V_e} \mathbf{B}^T \mathbf{L} \mathbf{B}^* dV^e \\ \int_{V_e} \mathbf{B}^{*T} \mathbf{L} \mathbf{B} dV^e & \int_{\Gamma_D} \mathbf{N}^{*T} \mathbf{C}_{loc} \mathbf{N}^* d\Gamma_D \end{matrix} \right] \qquad (19) $$

$$ \mathbf{F}_{coh} = \mathbf{F}_{int}^{0} - \int_{\Gamma_D} \mathbf{N}^{*T} D \mathbf{t}_0 d\Gamma_D \qquad (20) $$

où $\mathbf{N}^* = \mathbf{R}^{-1}\mathbf{N}$ avec $\mathbf{N}$ la matrice des fonctions d'interpolation, $\mathbf{R}$ la matrice de passage du repère global au repère local de la bande et $\mathbf{t}_0$ et $\mathbf{F}_{int}^0$ respectivement le vecteur de traction et le vecteur des efforts intérieurs de l'élément au moment de la localisation.

## 2.5 Étape de fissuration (XFEM)

L'étape de localisation est finalement suivie d'une fissuration dans le même plan. C'est-à-dire qu'une macro-fissure apparaît en lieu et place de la méso-fissure. Cela correspond à une perte totale de cohésion dès que :

$$ D = D_c \qquad (21) $$

A ce stade, les forces cohésives sont réduites à zéro, la bande de localisation devient une fissure sans résistance en traction, et l'élément devient XFEM standard. Dans ce travail on utilise la formulation de la base décalée et on néglige la partie asymptotique du champ de déplacement, voir [2]. Le champ de déplacement prend alors la forme :

$$ u(x) = u_{reg} + u_{dis} = \sum_{i \in I} N_i \alpha_i + \sum_{j \in J} (H(x) - H_j) N_j \beta_j \qquad (22) $$

L'avantage de cette formulation est d'éviter l'introduction d'éléments de transition entre les éléments fissurés et les éléments sains.

## 2.6 Aspects numériques

Dans le présent travail, le schéma d'intégration adopté est l'intégration moyennée sur le volume, voir [5]. Le schéma d'intégration est complété par l'utilisation de la technique B-Bar afin de traiter le verrouillage volumique.

## 2.7 Représentation spatiale de la bande de localisation

Dans la présente méthodologie, la localisation peut s'amorcer au sein de la structure dans un seul élément et se propager aux éléments adjacents. A noter que la porosité au sein d'un élément fini ne peut pas dépasser la porosité critique $f_c$. En effet le phénomène de localisation et la dégradation associée du matériau étant traités par la loi cohésive, on considère ici que l'endommagement ne peut pas évoluer au delà d'une valeur critique. Considérons deux éléments finis adjacents, un premier désigné $m$ contenant une bande de localisation de normale $n^{(m)}$, et un second désigné $j$ ne contenant aucune bande de localisation (voir Figure 3). Le but est de trouver l'angle $\phi$ satisfaisant l'un des deux critères de localisation (d'instabilité plastique où de coalescence des microcavités). Tout au long de l'analyse, certaines conditions qui
<figure><img src="image_4.png" /><figcaption>FIGURE 3 – Orientation de la bande</figcaption></figure>

empêchent artificiellement la bande de localisation de se propager peuvent se produire. Une telle situation est illustrée en Figure 4 où l'élément en rouge adjacent à l'EF contenant le front de bande satisfait le critère de pré-localisation mais ne répond à aucun des deux critères de localisation. Pour surmonter ce problème, il est nécessaire de réaliser une analyse de localisation non seulement dans le voisin le plus proche mais aussi dans un groupe de voisins contenus dans une pastille de rayon R perpendiculaire au front de la bande. L'orientation de la bande déterminée dans l'élément de couleur cyan est ensuite appliquée à l'élément de couleur rouge en respectant la continuité de la bande. Cette méthode n'est appliquée que dans le cas 1, où la localisation et induite par une instabilité plastique. Comme dans les approches de

<figure><img src="image_5.png" /><figcaption>FIGURE 4 – Méthode de patch hémisphérique</figcaption></figure>

Crété et al. [2] et Wolf et al. [8], la propagation de la bande de localisation suit le concept d'épuisement. Si un élément est déjà localisé soit par instabilité plastique où par coalescence des microcavités, la bande est introduite dans les éléments qui satisfont les critères suivants durant le même incrément :

1. Les éléments doivent satisfaire le critère de pré-localisation.

2. Les éléments doivent être connectés à l'élément déjà localisé et partager un côté.

3. La bande adjacente doit croiser ce côté partagé avec l'élément.

4. Les éléments qui satisfont 1,2 et 3 sont contrôlés par ordre de porosité décroissante en mode spirale autour d'un élément localisé.

5. Les éléments satisfaisant au moins l'un des deux critères de localisation sont localisés.

6. Les étapes 1 à 5 sont répétées jusqu'à ce qu'aucun élément nouveau de ne soit localisé durant l'incrément courant.

7. Si le calcul ne converge pas l'étape de propagation qui vient d'avoir lieu est annulée.

8. Si l'incrément de temps est réduit à moins de $10^{-5}$s aucune propagation n'a lieu.

## 3 Application

La méthodologie unifiée a été implémentée en tant qu'élément fini utilisateur (UEL) dans le code de calculs commercial par éléments finis ABAQUS (pour plus d'informations voir [4]). L'hypothèse des
petites déformations et une intégration temporelle implicite sont adoptées. Deux tailles de maillage différentes sont utilisées : à savoir grossier avec des dimensions d'éléments finis $\approx 1 \times 1 \times 1mm^3$ et fin $\approx 0.5 \times 0.5 \times 0.5mm^3$ dans la zone d'intérêt. Les propriétés du matériau sont reportés dans le Tableau 1. Comme souvent lors de l'utilisation de modèles de type GTN, le matériau a une porosité initiale de 0,1% et la porosité critique pour initier la localisation induite par la coalescence des microcavités est fixée à une valeur de 3%. On étudie ici trois géométries d'éprouvette : axisymétrique entaillée, plate entaillée

<table><caption>TABLE 1 – Paramètres liés à l'écrouissage, plasticité et loi cohésive</caption><thead><tr><th>E</th><th>ν</th><th>R<sub>0</sub></th><th>R<sub>∞</sub></th><th>k</th><th>q<sub>1</sub> = q<sub>2</sub></th><th>Δ<sub>c</sub></th><th>D<sub>c</sub></th><th>γ</th></tr></thead><tbody><tr><td>200 GPa</td><td>0.33</td><td>400 MPa</td><td>150 MPa</td><td>20</td><td>1</td><td>3 mm</td><td>0.5</td><td>2</td></tr></tbody></table>

et de cisaillement. La méthodologie unifiée développée est ici appliquée aux éprouvettes et permet la description de comportements matériau présentant un adoucissement post-localisation important. En Figure 5 est représentée la réponse globale en termes de force de réaction en fonction du déplacement pour l'éprouvette axisymétrique entaillée (a) et plate (b) et pour le cisaillement (c). La méthodologie pour les deux maillages est également comparée aux résultats obtenus avec le modèle de GTN. D'après les Figure 5 (a) et (b) on peut voir que pour les deux éprouvettes de traction la taille du maillage a peu ou pas d'effet sur les résultats obtenus. La dépendance pathologique du maillage est en effet surmontée grâce à l'utilisation de la méthodologie proposée combinant XFEM et CZM. Pour l'éprouvette de cisaillement (Figure 5 (c)) il y a une influence de la taille du maillage dans le régime d'écrouissage, c'est-à-dire bien avant le début de la localisation et donc bien avant l'activation de la méthode XFEM-CZM, puis dans le régime d'adoucissement. Ceci s'explique par la présence de trop peu d'éléments dans la région d'intérêt. En ce qui concerne le trajet de fissuration, des résultats presque identiques sont obtenus avec les deux

<figure><img src="image_6.png" /><figcaption>FIGURE 5 – Force de réaction vs Déplacement, a) éprouvette axisymétrique, b) éprouvette plate, c) éprouvette pour le cisaillement</figcaption></figure>

maillages pour tous les échantillons. La bifurcation attendue de la fissure, résultant de la transition du mode I au mode II, est bien reproduite pour l'éprouvette axisymétrique entaillée (Figure 6). A noter que dans le cas des éprouvettes entaillées, c'est le critère où la localisation est induite par la coalescence des vides qui est activé. Concernant l'éprouvette de cisaillement, c'est le critère 1 où la localisation est induite par instabilité plastique qui est activée. La forme de la loi d'adoucissement CZM influence directement la forme de la réponse globale dans le régime d'adoucissement. Divers comportements dans
<figure><img src="image_7.png" /><figcaption>FIGURE 6 – Surface de rupture des trois éprouvettes. Distribution de la porosité</figcaption></figure>

le spectre entre rupture brutale et progressive peuvent être modélisés en modifiant les valeurs des trois paramètres de la loi cohésive utilisée (voir Figure 7).

<figure><img src="image_8.png" /><figcaption>FIGURE 7 – Force de réaction pour différents paramètres de la loi cohésive pour l’éprouvette axisymétrique</figcaption></figure>

## 4 Conclusion

Il a été montré qu'une méthodologie unifiée (Figure 1) qui combine le modèle de GTN pour décrire l'endommagement diffus dans le cadre de la MEF standard, la XFEM pour représenter la fissure et le couplage de la XFEM avec un modèle de zone cohésive (CZM) pour traiter la phase de localisation et la dégradation induite des propriétés mécaniques du matériau est capable de reproduire numériquement la réponse globale de structures jusqu'à rupture, indépendamment de la taille de maille.

## 5 Perspectives

Étant donné que la rupture ductile peut entraîner des grandes déformations et rotations, un des objectifs sera d'étendre la méthode présentée ici aux grandes transformations. Pour ce faire, on peut appliquer l'approche incrémentale proposée par Bathe [1] en choisissant p.ex. la formulation lagrangienne mise à jour. Dans le cas d'un matériau avec un comportement bien défini dans la configuration courante et qui est soumis à des grandes déformations et rotations, il faut utiliser un algorithme d'intégration pour
passer d'un modèle constitutif basé sur le temps à une loi incrémentale (voir [6]) objective. Le principe est de réaliser la discrétisation temporelle dans un repère local qui neutralise les effets de rotation et puis retourner au repère global. On peut ainsi adopter la formulation de Jaumann par exemple. Pour traiter l'incompressibilité dans le matériau qui est responsable du verrouillage volumique, on ajoute l'approche F-Bar (voir [3]).

# 6 Références bibliographiques

## Références

[1] K.J. Bathe. Finite Element Procedures (Second Edition), K.J. Bathe Watertown MA, 2016.

[2] J.P. Crété, P. Longère, and J.M. Cadou. Numerical modelling of crack propagation in ductile materials combining the GTN model and X-FEM, Computer Methods in Applied Mechanics and Engineering, vol. 275, pp. 204-233, 2014.

[3] E.A. De Souza Neto, D. Perić, M. Dutko, D.R.J. Owen. Design of simple low order finite elements for large strain analysis of nearly incompressible solids, Int. J. Solids Structures, 33 : 3277-3296, 1996.

[4] K. Nikolakopoulos, J.P. Crété, P. Longère. Progressive failure of ductile metals : Description via a three-dimensional coupled CZM-XFEM based approach, Engineering Fracture Mechanics, 243 : 1-30, 2020.

[5] K. Nikolakopoulos, J.P. Crété, P. Longère. Volume averaging based integration method in the context of XFEM-cohesive zone model coupling, Mechanics Research Communications, 104 : p. 103485, 2020.

[6] J.C. Simo, T.J.R. Hughes. Computational Inelasticity, Springer-Verlag, 1998.

[7] V. Tvergaard, A. Needleman. Analysis of the cup-cone fracture in a round tensile bar, Acta Metallurgica, vol. 32, no. 1, pp. 157-169, 1984.

[8] J. Wolf, P. Longère, J.M. Cadou, and J.P. Crété. Numerical modeling of strain localization in engineering ductile materials combining cohesive models and X-FEM, International Journal of Mechanics and Materials in Design, vol. 14, no. 2, pp. 177-193, 2018.