# Simulation numérique de la déchirure ductile par la méthode des éléments de frontière

Gaëtan Hello, Hocine Kebir, Laurent Chambon, Alain Rassineux, Jean-Marc Roelandt

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

To cite this version:

Gaëtan Hello, Hocine Kebir, Laurent Chambon, Alain Rassineux, Jean-Marc Roelandt. Simulation numérique de la déchirure ductile par la méthode des éléments de frontière. 8e Colloque national en calcul des structures, CSMA, May 2007, Giens, France. hal-01504286

HAL Id: hal-01504286

https://hal.science/hal-01504286

Submitted on 9 Apr 2017

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.

Public Domain
# Simulation numérique de la déchirure ductile par la méthode des éléments de frontière

Gaétan Hello* - Hocine Kébir* - Laurent Chambon**
Alain Rassineux* - Jean-Marc Roelandt*

(*) Laboratoire Roberval
Univerité de Technologie de Compiègne
BP 20529 – 60205 COMPIEGNE CEDEX
gaetan.hello@utc.fr – hocine.kebir@utc.fr
alain.rassineux@utc.fr – jean-marc.roelandt@utc.fr

(**) EADS IW
18 rue Marius Terce BP 13050 – 31025 TOULOUSE CEDEX 3
laurent.chambon@eads.net

## RÉSUME :

La méthode des éléments de frontière est utilisée pour simuler numériquement le phénomène de déchirure ductile intervenant dans les structures minces fissurées. La résolution des problèmes mécaniques dans les milieux fissurés à comportement élastoplastique se base sur la formulation intégrale duale du problème. Cette méthode permettant de décrire avec précision les champs mécaniques en fond de fissure, le critère de rupture ductile de Rice et Tracey peut être employé en post-traitement du calcul. La validité du modèle numérique est testée en comparant ses résultats à ceux obtenus expérimentalement sur une éprouvette CCT.

## ABSTRACT :

This work aims at simulating the ductile fracture phenomenon occuring in thin cracked structures thanks to the boundary element method. A process based on the dual integral formulation of the problem has been developed so as to solve elastoplastic cracked problems. Since this method provides accurate mechanical fields at the crack tip, the Rice and Tracey criterion for ductile fracture can be evaluated during the post-processing step. The validity of the method is at last addressed with the comparison between numerical and experimental results for a test on a CCT specimen.

## MOTS-CLÉS :

méthode des éléments de frontière (BEM), mécanique non-linéaire de la rupture, approche locale de la rupture ductile, courbe-R.

## KEYWORDS :

boundary element method (BEM), non-linear fracture mechanics, local approach of ductile fracture, R-curve.
# 1. Introduction

La prise en compte de la tenue à la rupture joue un rôle capital dans le processus de dimensionnement des structures aéronautiques. En effet la recherche d'un gain de masse doit concilier efficacement les aspects économiques avec les impératifs de sécurité. Il convient donc de pouvoir prédire le comportement sous charge d'une structure fissurée.

S'agissant de structures métalliques minces, cette étude s'appuie classiquement sur la notion de courbe-R qui relie le taux de restitution d'énergie à l'avancée du défaut. En effet, contrairement aux structures fissurées épaisses qui rompent dès lors que le facteur K atteint une valeur critique, les structures minces peuvent présenter une phase de propagation stable relativement importante avant que ne survienne la rupture brutale. Cette stabilité est due à la compétition que se livrent en pointe de fissure les mécanismes irréversibles de dissipation plastique et d'avancée du défaut.

Pour une structure donnée, la courbe-R a comme propriété de dépendre principalement des paramètres de matériau et d'épaisseur. Il est donc possible d'extrapoler à une structure complète la courbe obtenue sur une éprouvette de mêmes épaisseur et matériau.

# 2. Contexte et objectifs

Pour l'heure la détermination de la courbe passe par des essais expérimentaux sur des éprouvettes normalisées. Des travaux antérieurs menés par l'UTC et EADS IW (anciennement EADS CCR) ont montré qu'il était néanmoins possible d'utiliser la méthode des éléments finis pour effectuer des calculs sur éprouvettes voire sur structure complète [DEL92]. Les travaux ici présentés entendent tirer parti des qualités de la méthode des éléments de frontière (bonne représentation des champs en fond de fissure, liberté face aux contraintes de maillage, gestion simplifiée de la propagation) pour permettre de déterminer numériquement les courbes-R. L'enjeu industriel immédiat est de réduire le nombre de configurations d'essais qui peuvent s'avérer coûteuses notamment pour celles de grande taille.

La première étape du travail a consisté à développer une méthode de résolution des problèmes de plasticité dans un milieu fissuré. La formulation éléments de frontière du problème se base sur la formulation intégrale duale présentée dans [ALI92]. La non linéarité matérielle est prise en compte au travers d'un algorithme d'intégration du comportement élastoplastique de type retour radial [BON99]. Les champs mécaniques en fond de fissure ainsi obtenus s'avérant particulièrement précis, il est dès lors envisageable de les exploiter pour alimenter le critère de rupture ductile de Rice et Tracey [RIC69] durant une phase de post traitement. L'obtention de la courbe-R recherchée ne requiert plus ensuite que la détermination de la valeur critique du paramètre de rupture au moyen d'un recalage avec des résultats expérimentaux.
### 3. Résolution des problèmes élastoplastiques fissurés

#### 3.1. Formulation intégrale duale

Souhaitant résoudre un problème de plasticité dans un milieu fissuré nous devons étendre la formulation intégrale classique. Celle-ci repose sur deux équations, l'une en déplacement pour les points du contour [1], l'autre en contrainte pour les points strictement à l'intérieur du domaine [2]. Dans ces deux équations interviennent les champs mécaniques recherchés : déplacements u et tensions t sur le contour, déformations plastiques $\varepsilon^p$ dans le domaine ainsi que les noyaux U, T, D, S, $\sigma$ et $\Sigma$ de la solution fondamentale de Kelvin [ALI02]

$$ \frac{1}{2} \cdot u_i(\underline{X}) = \left[ \begin{array}{l} \displaystyle \int_{\partial\Omega} U_{ij}(\underline{X};\underline{Y}) \cdot t_j(\underline{Y}) \cdot d\Gamma(\underline{Y}) - \int_{\partial\Omega} T_{ij}(\underline{X};\underline{Y}) \cdot u_j(\underline{Y}) \cdot d\Gamma(\underline{Y}) \\ \displaystyle + \int_{\Omega} \sigma_{ijk}(\underline{X};\underline{Y}) \cdot \varepsilon_{jk}^p(\underline{Y}) \cdot d\Omega(\underline{Y}) \end{array} \right] \quad [1] $$

$$ \sigma_{ij}(\underline{X}) = \left[ \begin{array}{l} \displaystyle \int_{\partial\Omega} D_{ijk}(\underline{X};\underline{Y}) \cdot t_k(\underline{Y}) \cdot d\Gamma(\underline{Y}) - \int_{\partial\Omega} S_{ijk}(\underline{X};\underline{Y}) \cdot u_k(\underline{Y}) \cdot d\Gamma(\underline{Y}) \\ \displaystyle + \int_{\Omega} \Sigma_{ijkl}(\underline{X};\underline{Y}) \cdot \varepsilon_{kl}^p(\underline{Y}) \cdot d\Omega(\underline{Y}) + f_{ij}(\varepsilon_{kl}^p(\underline{X})) \end{array} \right] \quad [2] $$

Cette formulation s'avère parfaitement opérante lorsque la géométrie est exempte de défaut [HEL05]. Toutefois l'introduction de la fissure, dont les deux lèvres sont géométriquement confondues, aboutit à des formules inexploitables numériquement.

Parmi les méthodes permettant de lever ce blocage, nous avons ici opté pour la formulation duale décrite dans [ALI92]. Elle présente les avantages de découler assez naturellement de la formulation classique et de ne pas nécessiter de traitement numérique trop spécifique. Elle s'appuie sur les équations classiques auxquelles on adjoint deux nouvelles équations intégrales pour les points de la fissure : l'une en déplacement [3] pour la lèvre inférieure, l'autre en tension [4] pour la lèvre supérieure.

$$ \frac{1}{2} \cdot [u_i(\underline{X}^+) + u_i(\underline{X}^-)] = \left[ \begin{array}{l} \displaystyle \int_{\partial\Omega} U_{ij}(\underline{X}^-;\underline{Y}) \cdot t_j(\underline{Y}) \cdot d\Gamma(\underline{Y}) - \int_{\partial\Omega} T_{ij}(\underline{X}^-;\underline{Y}) \cdot u_j(\underline{Y}) \cdot d\Gamma(\underline{Y}) \\ \displaystyle + \int_{\Omega} \sigma_{ijk}(\underline{X}^-;\underline{Y}) \cdot \varepsilon_{jk}^p(\underline{Y}) \cdot d\Omega(\underline{Y}) \end{array} \right] \quad [3] $$

$$ \frac{1}{2} \cdot [t_i(\underline{X}^+) - t_i(\underline{X}^-)] = n_j(\underline{X}^+) \cdot \left\{ \begin{array}{l} \displaystyle \int_{\partial\Omega} D_{ijk}(\underline{X}^+;\underline{Y}) \cdot t_k(\underline{Y}) \cdot d\Gamma(\underline{Y}) - \int_{\partial\Omega} S_{ijk}(\underline{X}^+;\underline{Y}) \cdot u_k(\underline{Y}) \cdot d\Gamma(\underline{Y}) \\ \displaystyle \int_{\Omega} \Sigma_{ijkl}(\underline{X}^+;\underline{Y}) \cdot \varepsilon_{kl}^p(\underline{Y}) \cdot d\Omega(\underline{Y}) + \frac{1}{2} \cdot [f_{ij}(\varepsilon_{kl}^p(\underline{X}^+)) - f_{ij}(\varepsilon_{kl}^p(\underline{X}^-))] \end{array} \right\} \quad [4] $$
## 3.2. Passage aux éléments de frontière et obtention des systèmes linéaires

Partant des équations intégrales [1] à [4], l'obtention de la forme discrète du problème passe dans un premier temps par la décomposition des intégrales sur des sous domaines d'intégration puis par l'approximation élémentaire des champs mécaniques [KEB98]. Les champs de déplacements et de tension sur le contour, pointes de fissure exclues, sont approchés localement par des formes polynomiales quadratiques 1D non-conformes. Ces mêmes champs sont en revanche approximés par des fonctions singulières 1D non-conformes en pointes de fissures [KEB99]. Les déformations plastiques dans le domaine sont quant à elles approchées par des polynômes linéaires 2D non-conforme.

Cette étape de discrétisation de la géométrie et des champs aboutit à deux systèmes linéaires [5] et [6] d'inconnues discrètes {$u$}, {$t$}, {$ε^p$} respectivement les vecteurs déplacements, tensions du contour et déformations plastiques du domaine.

$$[A] \cdot \{u\} = [B] \cdot \{t\} + [D] \cdot \{\epsilon^p\} \qquad [5]$$

$$\{\sigma\} = -[A'] \cdot \{u\} + [B'] \cdot \{t\} + [E] \cdot \{\epsilon^p\} \qquad [6]$$

## 3.3. Algorithme d'intégration du comportement élastoplastique

La résolution concrète du problème repose sur un algorithme élastoplastique implicite faisant intervenir les systèmes [5] et [6]. L'intégration locale de la loi de comportement s'effectue au moyen d'une méthode classique de retour radial [BON99]. L'algorithme présente donc un aspect incrémental lié à l'accroissement progressif de la sollicitation (charge ou déplacement selon la nature du pilotage) ainsi qu'un aspect itératif lié à la recherche de convergence sur chaque pas.

## 3.4. Résultats numériques

La validation de la méthode s'effectuera en comparant, pour un problème donné, les résultats numériques offerts par notre code éléments de frontière KSP avec ceux fournis par le code éléments finis ABAQUS. Le problème retenu est celui d'une plaque mince fissurée représentée sur la figure [1].

Les paramètres pertinents retenus pour la comparaison sont les composantes de déplacement des points de la fissure. Ces quantités présentent le double intérêt d'être à la fois représentatives de la bonne prise en compte de la plasticité mais également d'être des quantités immédiatement accessibles tant dans la formulation FEM que BEM.
En pratique le modèle éléments finis de référence comprend 4350 éléments CPS6 (contraintes planes quadratique) (fig. 2), le modèle éléments de frontières ayant quant à lui 492 cellules et 60 éléments de contour (fig. 3).

<figure><img src="image_2.png" /><figcaption>**Figure 1.**
*Problème*</figcaption></figure>

<figure><img src="image_3.png" /><figcaption>**Figure 2.**
*maillage FEM*</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>**Figure 3.**
*Discrétisations BEM*</figcaption></figure>

La comparaison entre KSP et ABAQUS s'effectue pour des chargements croissants de 100MPa, 120MPa, 140MPa et 160MPa. On notera que pour cette dernière charge la zone plastique n'est absolument plus confinée en pointe de fissure, elle occupe en effet la quasi-totalité de la discrétisation BEM de surface.

Concrètement on observe sur la figure 4 que les fissures émoussées obtenues par KSP coïncident parfaitement avec celles prédites par les calculs de référence ABAQUS pour les différents chargements considérés.

<figure><img src="image_5.png" /><figcaption>**Figure 4. Comparaison des fissures déformées entre KSP et ABAQUS**</figcaption></figure>
## 4. Simulation numérique de la déchirure ductile

Ayant développé un solveur BEM pour les problèmes élastoplastiques fissurés, il est désormais envisageable de lui adjoindre le critère micro-mécanique de rupture de Rice et Tracey pour simuler le phénomène de déchirure ductile. Ce critère [7] se base sur les champs de contraintes et de déformations plastiques en pointes de fissure pour estimer la croissance de micro-cavités aboutissant à la propagation du défaut.

$$ \frac{\dot{R}}{R} = 0.283 \cdot \exp\left(\frac{3}{2} \cdot \frac{\sigma_h}{\sigma_{VM}}\right) \cdot \dot{p} \quad [7] $$

Ce critère évalué à l'issue de chaque incrément de sollicitation permet de statuer sur l'éventuelle avancée du défaut. Il est alors possible de décrire numériquement les phénomènes irréversibles de propagation de fissures et de plasticité.

## 5. Conclusion

Les résultats numériques concluants obtenus dans la résolution de problèmes élastoplastiques fissurés au moyen de la méthode des éléments de frontière ouvrent des perspectives dans l'application de cette méthode à la simulation de la déchirure ductile. L'obtention finale des courbes-R passera par un calage du paramètre utilisé dans le modèle numérique de rupture avec des données expérimentales.

## 6. Bibliographie

[ALI02] Aliabadi M.H., "The boundary element method volume 2 : applications in solid and structures", John Wiley and Sons (2002).

[ALI92] Y. Mi and M. H. Aliabadi, "Dual boundary element method for three-dimensional fracture mechanics analysis", *Engineering Analysis with Boundary Elements*, Volume 10, Issue 2, 1992, Pages 161-171.

[BON99] Bonnet M., "Boundary Integral Equations Methods for Solids and Fluids", John Wiley and Sons (1999).

[DEL92] Delmotte J., « Modélisation numérique de la déchirure ductile en milieu bidimensionnel à l'aide d'une approche locale. Simulation des courbes de résistance à la fissuration », *Thèse de doctorat*, Université de Technologie de Compiègne.

[HEL05] Hello G, « Application de la méthode des éléments de frontière à la résolution de problèmes élastoplastiques », *Actes Giens 2005*, Hermès, pp 445-450.

[KEB99] H. Kebir, J. M. Roelandt and J. Foulquier, "A new singular boundary element for crack problems: Application to bolted joints", *Engineering Fracture Mechanics*, Volume 62, Issue 6,April 1999,Pages 497-510.
[KEB98] H. Kebir, « Approches déterministe et probabiliste de la prévision de durée de vie de structures aéronautiques à l'aide de la méthode des équations intégrales duales », *Thèse de doctorat*, Université de Technologie de Compiègne.

[RIC69] Rice J.R., Tracey D.M., "On the ductile enlargement of voids in triaxial stress fields", *J. Mech. Phys. Solids*, 17, 201.