# Intégrale invariante pour la propagation de fissure en élastoplasticité sous chargement de fatigue

Vincent Riolo, Vincent Chiaruttini, Joachim Guilié, Frédéric Feyel

► To cite this version:

Vincent Riolo, Vincent Chiaruttini, Joachim Guilié, Frédéric Feyel. Intégrale invariante pour la propagation de fissure en élastoplasticité sous chargement de fatigue. 11e colloque national en calcul des structures, CSMA, May 2013, Giens, France. hal-01717109

HAL Id: hal-01717109

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

Submitted on 25 Feb 2018

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 CC0 1.0 - Universal - International License
# Intégrale invariante pour la propagation de fissure en élastoplasticité sous chargement de fatigue

Vincent Riolo<fnref n="1" />,<fnref n="2" />, Vincent Chiaruttini<fnref n="1" />,*, Joachim Guilié<fnref n="2" />, Frédéric Feyel<fnref n="1" />

**Résumé** — Afin d'étendre les approches issues de la mécanique linéaire de la rupture aux cas des matériaux standards généralisés et prendre en compte des effets liés aux phénomènes de plasticité, on se propose d'étudier une grandeur énergétique globale, calculée par une intégrale invariante. Une première étude numérique du comportement asymptotique de cette grandeur est menée sur des matériaux élastoplastique à écrouissage cinématique ou isotrope. Une stratégie basée sur l'introduction d'une longueur interne est alors proposée de façon à montrer qualitativement l'intérêt d'une telle approche pour alimenter une loi de propagation en fatigue permettant de prendre en compte les effets d'histoire en front de fissure. Enfin une application en propagation 2D montrant l'influence des effets de surcharges est présentée.

**Mots clés** — fissuration, paramètre énergétique de rupture, matériau standard généralisé, méthode G-θ.

## 1 Introduction

Le haut niveau de sophistication des composants mécaniques industriels dans un grands nombre de domaines (aéronautique, automobile, ferroviaire, etc.), conduit à soumettre des pièces critiques à des chargements en service de plus en plus exigeants. Aussi, il n'est plus rare aujourd'hui d'envisager dès le bureau d'étude, l'apparition de fissures dans des composants fortement sollicités dont il va être nécessaire de savoir prédire efficacement l'évolution de façon à déterminer, en toute sécurité, leur durée de vie et la durée des intervalles inter-inspection. Pour atteindre cet objectif, il est essentiel de disposer à la fois de modélisations fines de ces problèmes physiques et d'outils numériques efficaces pour réaliser des études fiables de prédiction de l'évolution des fissures.

Il existe depuis une dizaine d'années des méthodes très performantes pour étudier les phénomènes de fissuration dans le cadre de la mécanique linéaire de la rupture (basées par exemple, sur la méthode X-FEM [1] ou des techniques par maillage conforme [2]). Cependant, accéder à des problématiques hors de ce cadre (plasticité confinée et généralisée, chargement multi-axiaux, contact, etc.), pour des problématiques industrielles, demeure particulièrement complexe.

On se propose ici d'étendre le cadre de l'approche global de la fissuration (basée sur la représentation d'une fissure macroscopique et le calcul d'une grandeur dissipative) au cas des matériaux standards généralisés. Dans une démarche pragmatique, on utilise ici une grandeur invariante basée sur un calcul de type G-θ [3] avec prise en compte de la plasticité comme un état de contrainte résiduel figé lors de la propagation. Cette approche qui possède des similitudes avec l'intégrale GTP [4], a été étudiée dans le contexte d'une véritable fissure (et non plus d'une entaille), pour des chargement de fatigue. Une étude numérique du comportement asymptotique est menée de façon à montrer numériquement l'invariance de cette intégrale pour différents types de comportement élastoplastiques (à écrouissage cinématique ou isotrope). Nous montrerons qu'en utilisant un maillage suffisamment fin et une longueur caractéristique, il est possible de définir une grandeur dissipative dont l'évolution permet qualitativement de représenter des effets d'histoire en pointe de fissure (effet de refermeture et influence des surcharges).

La section suivante présente la grandeur énergétique étudiée dans ces travaux. La section 3 décrit l'étude numérique du comportement asymptotique en pointe de fissure sur un problème élastoplastique 2D. En utilisant une longueur interne et la grandeur calculée, dans la section 4, une loi de propagation de fissure en fatigue de type Paris est utilisée pour rendre compte facilement des effets de surcharge. Enfin,

<footnote n="1">¹ Onera - The French Aerospace Lab, DMSM/MNU. 29 avenue de la Division Leclerc, F-92322 Chatillon Cedex,</footnote>

<footnote n="2">² LMS, École Polytechnique, 91128 Palaiseau Cedex, FRANCE,</footnote>

<footnote>* Auteur correspondant {vincent.chiaruttini@onera.fr}</footnote>
quelques conclusions et perspectives sur ces travaux seront présentées.

## 2 Grandeur caractéristique pour la fissuration élastoplastique

L'étude de la propagation de fissure dans un milieu élastique fragile est largement présentée dans la littérature [6, 5]. En considérant, un domaine soumis à un chargement dans lequel une fissure plane de longueur $l$ est entrain de se propager et dont le comportement est élastique, on peut montrer que la dissipation s'écrit :

$$D = - \frac{\partial W_{ep}}{\partial l} \Big|_{\mathcal{E}} (\mathcal{E}, l) \frac{dl}{dt}. \qquad (1)$$

avec $\mathcal{E}$ une désignation générique du chargement. La force motrice de la propagation de fissure apparaît ainsi comme la dérivée de l'énergie potentielle par rapport à la position de la fissure à chargement fixé. Cette démonstration a été successivement étendue par [6] et par [5] de manière sensiblement équivalente au cas d'une fissure fragile propageant dans un milieu dissipatif. L'hypothèse de fissure fragile suppose que les processus dissipatifs ne sont pas couplés au processus de rupture. Dans ce cas, les auteurs obtiennent formellement

$$D = - \frac{\partial W_{ep}}{\partial l} \Big|_{\mathcal{E},k} (\mathcal{E}, l, k) \frac{dl}{dt} - \frac{\partial \psi_{libre}}{\partial k} \Big|_{\mathcal{E},l} (\mathcal{E}, l, k) \frac{dk}{dt}. \qquad (2)$$

Dans ce cas la force motrice de la propagation apparaît comme la dérivée de l'énergie potentielle par rapport à la position de la fissure à chargement et variables internes fixées dans la configuration variable. Le second terme de la dissipation totale représente la dissipation matériau classique.

On propose de considérer un problème de structure fissurée en petites déformations possédant un comportement rentrant dans le cadre des matériaux standards généralisés. Nous considérons un solide élastique fissuré occupant le domaine $\Omega \in \mathbb{R}^n, n = 2$ ou 3, soumis à un champ de forces surfaciques $f$ appliquées sur une partie $S$ de $\partial\Omega$ et des déplacements imposés sur la partie complémentaire $U$ de $\partial\Omega$.

Le comportement de la structure est tel qu'en tout point de la structure la relation de comportement suivante est vérifiée :

$$\sigma = E : (\varepsilon - \varepsilon^{ae}) \qquad (3)$$

avec $\sigma$ le tenseur des contraintes, $\varepsilon$ le tenseur des déformations totales linéarisé, $\varepsilon^{ea}$ le tenseur des déformations anélastiques (d'origine thermo-élastoplastique, par exemple) et $E$ le tenseur d'élasticité. Notons que pour définir l'intégrale invariante, il suffit de calculer pour n'importe quel état de la structure à l'équilibre le tenseur des déformations anélastiques en utilisant l'inverse de l'opérateur de rigidité, soit par :

$$\varepsilon^{ae} = \varepsilon - E^{-1} : \sigma \qquad (4)$$

Conformément au cadre formulé dans [3] pour la méthode G-θ, un champ de vitesses virtuel d'avancée du front de fissure $\theta$ est introduit. Nous faisons alors l'hypothèse que les champs de forces surfaciques, ainsi que le tenseur des déformations anélastiques, ne sont pas modifiés par la propagation de la fissure. Ce qui se traduit par :

$$f^\cdot = \nabla f \cdot \theta = 0 \text{ sur } S \quad ; \quad \dot{\varepsilon}^{ae} = \nabla \sigma^{ae} \cdot \theta \qquad (5)$$

Dans ces conditions, on peut définir la quantité $G_P(M)$ en chaque point $M$ du front de fissure $\Gamma_0$ par la solution de l'équation variationnelle suivante :

$$\int_{\Gamma_0} G_P(s) \theta(s) \cdot n(s) ds = - \int_{\Omega} \left[ \frac{1}{2} (\sigma : (\varepsilon - \varepsilon^{ae})) \nabla \cdot \theta - \sigma : (\nabla u \nabla \theta) - (\sigma : \nabla \varepsilon^{ae}) \cdot \theta \right] d\Omega, \forall \theta \in \Theta \qquad (6)$$

Cette équation est valable quels que soient les champs $\theta$ respectant les conditions suivantes :

$$\theta \in \Theta = \{\mu \text{ tels que } \mu \cdot n = 0 \text{ sur } \partial\Omega\} ; \quad \theta_{S \cup S_d} = 0 ; \quad \nabla \theta_{S \cup S_d} = 0 \qquad (7)$$

En effet, les transformations à considérer ne doivent modifier que la position du fond de fissure $\Gamma_0$.

Cette formulation, qui produit une intégrale formellement invariante en considérant l'ensemble des déformations anélastiques comme figées dans la transformation, est à rapprocher de la formulation GTP
introduite dans [4]. Cependant, elle s'avère plus simple et générale puisqu'aucune hypothèse n'est re-
quise sur le type de matériau considéré.

Nous proposons dans la suite de valider numériquement l'invariance d'une telle grandeur dont la principale difficulté de calcul tient à l'évaluation du terme contenant le gradient des déformations anélastiques.

L'étude numérique nous amène ensuite à préciser la difficulté principale de ce type de calcul, la valeur en pointe de fissure de l'intégrale des gradients spatiaux de variables internes. L'étude de cette difficulté n'est pas l'objet de cette article et on propose d'introduire une longueur caractéristique pour s'en défaire. On notera dans la suite que si l'on choisit de ne pas tenir compte de ce terme de gradient, la grandeur ainsi calculée sera notée G.

3 Analyse numérique de l'intégrale

La mise en œuvre du calcul de l'intégrale préalablement définie a été réalisée au sein du code Z-
set/Zebulon [7] pour application au module de fissuration conforme Z-cracks [2], co-développé à l'Onera.
Le calcul du gradient des déformations anélastiques est obtenu en utilisant un processus de lissage par
moindre carré mobile (en pratique 8 à 16 points permettent d'obtenir de bons résultats).

Afin de montrer l'intérêt de cette grandeur pour des fissuration en élastoplasticité, une première étude
de convergence est réalisée ici. On considère ainsi un modèle test représentatif des effets de plasticité
confinée agissant en pointe de fissure. Le type de maillage utilisé est présenté sur la figure 1, il s'agit
d'une structure fissurée circulaire 2D sur laquelle est imposée directement sur la frontière extérieure un
chargement en déplacements imposés donné par une solution de Westergaard en mode I :

$$
u_x = \sqrt{r} \cos \frac{\alpha}{2} \left( \kappa - 1 + 2 \sin^2 \frac{\alpha}{2} \right) ; \quad u_y = \sqrt{r} \sin \frac{\alpha}{2} \left( \kappa - 1 - 2 \cos^2 \frac{\alpha}{2} \right) \quad (8)
$$

où (r, α) sont les coordonnées polaires du point considéré et κ = 3 − 4ν avec ν le coefficient de Poisson du matériau.

<figure><img src="image_2.png" /><figcaption>Fig. 1 – Maillage représentatif du problème de plasticité confinée en pointe de fissure (iso-valeurs de plasticité cumulée).</figcaption></figure>

De cette façon, en vue d'une analyse de convergence en maillage, la longueur de la fissure disparaît de la modélisation est seule la taille de la zone plastique est déterminante. En pratique, on calibre les niveaux de chargement de sorte que la zone plastique reste confinée dans le maillage d'étude. Nous avons ainsi, avec un maillage de rayon unitaire, généré des éléments de taille caractéristique jusqu'à 10⁻⁶ pour un rayon plastique de l'ordre de 0, 1. Aussi les courbes présentées pour la valeur des quantités intégrales obtenues dans la suite sont adimensionnées (en pratique Gₚ et G).

3.1 Cycle charge/décharge

Dans cette partie on considère un cycle de chargement simple : un montée en charge en traction
(t ∈ [0; 1]), puis un retour à une frontière extérieure à déplacement imposé nul (t ∈ [1; 2]). Deux types
de matériaux sont étudiés, tous les deux basés sur un comportement isotrope linéaire (module d'Young $E = 200000$ MPa, $v = 0,3$) avec un seuil de plasticité suivant un critère de von Mises de limite initiale $\sigma_0 = 300$ MPa.

On considère un comportement élastoplastique à écrouissage isotrope linéaire de module $H = 50000$ MPa ou à écrouissage cinématique linéaire de module $C = 50000$ MPa. Avec ce type de modèles, la contrainte garde une singularité en pointe de fissure et dans ces conditions, le calcul numérique de la quantité en gradient $\sigma : \nabla\varepsilon^{ae}.\theta$ n'est pas fiable sur les premières couches d'éléments, aussi la limite suivante reste indéterminée :

$$ \lim_{r \to 0} \int_{\Gamma_0^r} (\sigma : \nabla\varepsilon^{ae}.\theta) dV \quad \text{avec} \quad \Gamma_0^r = \{M \in \Omega, d(M, \Gamma_0) < r\} \qquad (9) $$

On choisit arbitrairement qu'en dessous d'un rayon $r_0$, cette quantité sera supposée nulle et évaluée par moindres carrés mobiles uniquement sur l'intervalle $[r_0, r_\theta]$, si $r_\theta$ est le rayon du support de la fonction d'extension $\theta$.

<figure><img src="image_3.png" /><figcaption>Fig. 2 – Quantités énergétiques calculées pour un cycle charge/décharge avec les deux modèles élastoplastiques.</figcaption></figure>

La figure 2 illustre les quantités intégrées pour les deux modèles en utilisant deux types de maillages : l'un très raffiné avec une taille de maille $hmin$ atteignant $1.e-6$ et l'autre plus grossier $1.e-2$. Le calcul est mené avec quatre étapes :

1. un état de charge intermédiaire $t = 0,5$;

2. un état de chargement maximal $t = 1,5$;

3. un état de décharge intermédiaire $t = 1,5$;

4. un état de décharge complète $t = 2$;

Ainsi, plusieurs courbes quasi horizontales sont tracées, elles représentent soit des valeurs de $G$ (courbe verte maillage fin et violette maillage grossier), soit des courbes de $G_P^{r_0}$ (courbe rouge pour maillage fin avec $r_0 = 2.10^{-6}$, bleu foncé et bleu clair pour maillage fin et grossier avec $r_0 = 5.10^{-2}$ - avec dans ce cas deux dernières courbes très proches et invariantes).

Quel que soit le modèle, pour les deux états de montée en charge, on obtient des résultats identiques, avec une quasi-invariance de l'intégrale $G$ pour des très faibles rayons d'intégration (1/1000 du rayon plastique). Cependant cette justification éventuelle d'une valeur tendant vers zéro pour le terme de gradient, en cas de chargement monotone, n'est plus valide pour la décharge avec écrouissage isotrope.
Pour contourner la difficulté de détermination de la valeur du gradient initial, on peut observer que pour les deux types de comportement, d'une part, une fois la valeur $r_0$ choisie la valeur de $G_P$ reste toujours invariante, et d'autre part, le calcul de $G$ mené sur un maillage beaucoup grossier coïncide exactement dans sa zone de validité avec le résultat sur maillage fin. Ainsi, même si la valeur initiale (pour $r < r_0$) de la grandeur en gradient reste une indétermination, il est possible de définir une valeur stable de l'intégrale $G_P$ associée à une valeur de $r_0$ que l'on peut assimiler à une longueur caractéristique représentative, par exemple, de l'échelle à laquelle a lieu le phénomène de refermeture de fissure en décharge.

On propose ainsi, ici de fixer cette longueur à $R_0 = 0,05$ ce qui correspond dans cet exemple à environ 1/5ème du rayon plastique. Il est possible avec cette valeur d'obtenir un résultat très proche de la solution de référence (du maillage fin), avec le maillage grossier beaucoup moins coûteux, tout en restant parfaitement invariant pour $r_0 > R_0$. Sur ce premier exemple, il semble observable que la quantité ainsi définie $G_P^{R_0}$ est sensible à l'histoire du chargement, puisque les valeurs ne sont pas identiques dans les états de sollicitations d'amplitude identique dans les phases de charge et de décharge (avec une valeur qui devient même négative en décharge complète - qui coïncide avec une éventuel état de compression résultant de la plasticité engendrée en traction).

## 3.2 Cycles avec surcharge

Pour valider l'intérêt d'une telle grandeur pour des problématiques de fatigue, nous avons considéré dans cette partie une sollicitation plus complexe agissant sur notre problème test : il s'agit d'une mission de chargement avec surcharge (décrite en trait continu sur la figure 4). Deux charges/décharges sont appliquées avec un premier niveau de sollicitation, une surcharge de 50% est ensuite imposée, avant deux nouveaux cycles d'amplitude initiale.

<figure><img src="image_4.png" /><figcaption>Fig. 3 – $G_P^{R_0=0,05}$ calculé lors d'une mission de chargement avec surcharge, comparaison entre modèles de comportement pour maillages fin ou grossier.</figcaption></figure>

Dans ces conditions, on observe un fort effet de la surcharge en utilisant un comportement élasto-plastique à écrouissage isotrope : une réduction de plus de 25% de la valeur maximale dans le cycle de chargement initial est observé. De plus, il n'y pas de différence significative entre les valeurs obtenues avec les maillages fin et grossier.

# 4 Analyse qualitative d'un modèle de fissuration élastoplastique

Dans cette dernière partie on propose de montrer qualitativement l'intérêt de la grandeur calculée préalablement pour un cas de fissuration en plasticité confinée avec surcharge. On considère un problème
2D sur une éprouvette fine carrée de dimensions 10 mm × 10 mm possédant un trou de rayon 1 mm en son centre. La structure possède un comportement élastoplastique à écrouissage isotrope linéaire décrit dans la section précédente. Un chargement cyclique de traction verticale avec un effort d'amplitude 100 MPa est appliqué sur les parties supérieure et inférieure.

<figure><img src="image_5.png" /><figcaption>Fig. 4 – Eprouvette percée et fissurée soumise à un chargement de traction verticale.</figcaption></figure>

Un fissure de longueur initiale 0,1 mm est présente sur le côté droit de l'éprouvette. On choisit d'utiliser un loi de Paris caractérisée par un coefficient $C = 10^{-11}$ et un exposant $m = 3$. Cette loi de propagation est alimentée par une valeur en variation de $G_p^{R0}$ pour $R_0 = 0,01$ mm transformée en facteur d'intensité de contrainte $K$ équivalent, vérifiant :

$$K = \sqrt{\frac{E G_p^{R0}}{1 - v^2}} \quad ; \quad \frac{\Delta a}{\Delta N} = C (\Delta K)^m \quad (10)$$

où $a$ est la longueur d'avancée de fissure et $N$ le numéro représentatif de la mission, l'unité retenue dans la simulation pour $\Delta K$ est en MPa$\sqrt{mm}$.

<figure><img src="image_6.png" /><figcaption>Fig. 5 – Traînée plastique plastique observée avec un pas d'avancée de fissure de 0,01 mm sans surcharge (iso-valeurs de plasticité cumulée).</figcaption></figure>

Pour réaliser cette simulation nous avons utilisé le module 2D de l'outil Z-cracks de fissuration par remaillage adaptatif. Un maillage avec une taille d'éléments finis jusqu'à 0,005 mm a été utilisée. En pratique un incrément d'avancée de fissure $\Delta a$ est défini, une mission de chargement est effectuée, et à partir du $\Delta K$ calculé, un incrément de mission $\Delta N$ correspondant est déterminé. Le maillage est alors transformé de façon à avancer le front de fissure. Un algorithme de transfert de champ est requis pour correctement représenter les champs de déplacement discontinus et les informations aux points de Gauss lors de la mise à jour du maillage (ici, par collocation nodale et recherche du point de Gauss le plus proche).
La première validation a conduit sur un cas de chargement simple, à étudier la convergence de la simulation en fonction de l'incrément d'avancée de fissure $\Delta a$. La figure 6 montre une très nette surévaluation de la valeur de $\Delta K$ pour un pas d'avancée de fissure de 0,05 mm, qui produit une évolution beaucoup plus rapide de la fissure (cf. figure 7). Il semble, que si l'incrément d'avancée est trop important, le calcul n'est plus sensible à la traînée plastique qui fait chuter la valeur du facteur d'environ 15% en régime établi par rapport au premier cas de chargement (cf. figure 5). Ainsi, seules les courbes avec des incréments d'avancée de fissure de 0,01 mm ou 0,001 mm produisent des valeurs voisines (bien qu'oscillantes) et des vitesses de propagation très proches.

<figure><img src="image_7.png" /><figcaption>Fig. 6 – Comparaison des valeurs de $\Delta K$ en fonction de l'avancée de fissure.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>Fig. 7 – Comparaison des avancées de fissure en fonction du nombre de cycles.</figcaption></figure>

Le cycle de chargement est alors modifié : à chaque fois que la fissure progresse de 0,1 mm, une surcharge de contrainte de 50% est appliquée. En utilisant le paramètre proposé au sein d'une loi de Paris classique et un pas d'avancée de fissure suffisamment fin, il est possible de rendre compte qualitativement d'un ralentissement de la fissure (observable par la diminution de la valeur de $\Delta K$ juste après la surcharge), en accord avec les observations expérimentales [8].
# 5 Conclusion

Ces travaux présentent une grandeur énergétique caractéristique du phénomène de dissipation d'énergie en fissuration et plasticité confinée. En introduisant une longueur caractéristique utile pour la bonne définition de cette grandeur, elle est capable de rendre compte des effets d'histoire du chargement sur structures élastoplastiques soumises à des sollicitations de fatigue avec surcharge. Il conviendrait de poursuivre ces investigations pour établir un véritable modèle de fissuration confronté à des observations expérimentales.

# Références

[1] N. Moës, J. Dolbow, T. Belytschko. A finite element method method for crack growth without remeshing, Int. Journal for Num. Meth. In Engrg, 46, 131-150, 1999.

[2] V. Chiaruttini, F. Feyel, J.L. Chaboche, A robust meshing algorithm for complex 3D crack growth simulation, ECCM 2010, Paris.

[3] Ph Destuynder and M Djaoua and S Lescure.Quelques remarques sur la mécanique de la rupture élastique, Journal de Mécanique Théorique et Appliquée, 2(1), 113-135, 1983.

[4] G. Debruyne, Proposition d'un paramètre énergétique de rupture pour les matériaux dissipatifs, CRAS 328, Iib, 785-791, 2000.

[5] S.T. Nguyen, L. Dormieux, Y. Le Pape, J. Sanahuja. Crack propagation in viscoelastic structures : Theroretical and numerical analyses. Computational Materials Science, 50(1), 83-91.

[6] J. Tillberg, F. Larsson, K. Runesson. On the role of material dissipation for the crack-driving force. International Journal of Plasticity, 26(7), 992-1012.

[7] Frédéric Feyel. Some new technics regarding the parallelisation of zébulon, an object oriented finite element code for structural mechanics. Mathematical Modelling and Numerical Analysis, 36(5):923-935, sep 2002.

[8] P.Y. Decreuse, S. Pommier, L. Gentot, S. Pattofatto. History effect in fatigue crack growth under mixed-mode loading conditions. International Journal of Fatigue, 31(11), 1733-1741.