# Approche énergétique implicite pour la propagation de fissures en 3D. Prise en compte des contraintes résiduelles

Gwénaël Edeline, Vincent Chiaruttini, Marc Bonnet, Frédéric Feyel

► To cite this version:

Gwénaël Edeline, Vincent Chiaruttini, Marc Bonnet, Frédéric Feyel. Approche énergétique implicite pour la propagation de fissures en 3D. Prise en compte des contraintes résiduelles. 8e Colloque national en calcul des structures, CSMA, May 2007, Giens, France. pp.723–728. hal-01500105

HAL Id: hal-01500105

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

Submitted on 2 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.

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

Distributed under a Creative Commons CC0 1.0 - Universal - International License
# Approche énergétique implicite pour la propagation de fissures en 3D

## Prise en compte des contraintes résiduelles

G. Edeline* — V. Chiaruttini** — M. Bonnet* — F. Feyel**

* Laboratoire de Mécanique des Solides (CNRS-X)
Ecole Polytechnique
91128 Palaiseau Cedex
{edeline, bonnet}@lms.polytechnique.fr

** Onera – DMSE/LCME
BP 72, 29 avenue de la Division Leclerc
F-92 322 Châtillon cedex
{vincent.chiaruttini, frederic.feyel}@onera.fr

**RÉSUMÉ.** Nous proposons une nouvelle écriture des dérivées de l'énergie potentielle par la méthode θ prenant en compte un champ de contraintes résiduelles pour traiter de manière explicite ou implicite la propagation de fissure tridimensionnelle, aussi bien dans le cadre de la rupture par fatigue que celui de la rupture brutale. Des premiers éléments de validation numériques sont présentés à travers l'étude d'un problème fissuration en thermo-élasticité.

**ABSTRACT.** A novel way of expressing the derivatives of the potential energy by the θ-method which takes account a residual stress field is presented. It provides the ability to numerically simulate the propagation of a three-dimensional fatigue crack using explicit or implicit algorithms. A first numerical assessment is carried out through the study of a thermo-elastic fracture problem.

**MOTS-CLÉS:** méthode G-θ, fissure 3D, rupture par fatigue, contraintes résiduelles.

**KEYWORDS:** G-θ method, 3D crack, fatigue, residual stresses.
# 1. Introduction

L'étude des propagations de fissures est depuis longtemps un sujet crucial pour les industriels ; il existe un besoin toujours croissant de modèles et d'outils numériques robustes et prédictifs dans ce domaine. Généralement, les méthodes permettant la simulation des phénomènes de fissuration utilisent des approches locales ou globales qui conduisent au calcul de facteurs d'intensité de contrainte et de taux de restitution d'énergie.

Pour obtenir numériquement le taux de restitution d'énergie sur des configurations tridimensionnelles, certaines méthodes sont basées sur une extension du front de fissure (Parks, 1974), d'autres utilisent des intégrales invariantes (Rice, 1968). La méthode globale $G-\theta$, introduite par (Destuynder et al., 1983), est ici retenue : une transformation virtuelle de la structure dans la zone fissurée est utilisée pour obtenir une définition rigoureuse du taux de restitution d'énergie à l'aide de dérivations par rapport à la transformation. Elle présente deux avantages : (i) elle n'est pas limitée à l'élasticité linéaire, et (ii) elle permet un calcul du taux de restitution de l'énergie robuste par rapport au raffinage de maillage et basé sur des éléments finis ordinaires.

Les fissures peuvent être amenées à se propager au sein d'un champ de contraintes résiduelles préexistant à la fissure et dû à divers phénomènes (plastification, contraintes thermiques ou d'élaboration, ...) subis par le matériau dans son histoire. De nombreuses études (particulièrement en fatigue (Jung, 1997)) ont montré que l'effet de ces contraintes résiduelles sur la propagation de fissure est non négligeable.

Dasn ce contexte, nous avons étendu la méthode $\theta$ (Mialon, 1988) et les travaux réalisés à l'ONERA (Ousset, 1999) aux dérivées premières et secondes de l'énergie potentielle à l'équilibre pour un milieu élastique présentant des contraintes initiales. Notre formulation est également applicable aux machines tournantes (par exemple) en introduisant les effets d'inertie sous forme de forces de volume. Pour valider cette approche et mener des simulations numériques sur des cas industriels, une mise en oeuvre de ces expressions des dérivées de l'énergie potentielle a été réalisée au sein du code ZeBuLoN<fnref n="1" />.

# 2. Schéma implicite de simulation de propagation de fissure tridimensionnelle

Le problème de propagation quasi-statique revient à déterminer le champ $\theta$ (champ de vitesse virtuel d'avancée du front de fissure) pour lequel l'énergie totale E de la structure est minimale, ainsi :

$$E^{(1)}(\underline{\theta}) = 0 \Leftrightarrow D^{(1)}(\underline{\theta}) + J^{(1)}(\underline{\theta}) = 0 \quad [1]$$

où $D^{(1)}(\underline{\theta})$ et $J^{(1)}(\underline{\theta})$ sont respectivement les dérivées premières du potentiel de dissipation (énergie de création de surface) et de l'énergie potentielle (calculée par la méthode $\theta$).

<footnote n="1">1. Disponible sur Internet à http://www.nwnumerics.com.</footnote>
Pour résoudre l'équation [1], on utilise une méthode de Newton-Raphson en utilisant la dérivée seconde de l'énergie totale. Chaque itération consiste à déterminer une correction $\Delta\theta$ de la configuration de la fissure en résolvant l'équation linéarisée :

$$E^{(2)}(\Delta\theta, \underline{\theta}) = -E^{(1)}(\underline{\theta}), \forall \underline{\theta} \text{ admissible} \quad [2]$$

Cette base permet de mettre en oeuvre plusieurs approches de la simulation de propagation de fissure : (i) la rupture critique en implicite ou explicite par un critère de Griffith [3]. (ii) la propagation sous-critique en utilisant soit la loi de Paris (traitement explicite), soit une loi de Paris retournée [4] (traitement implicite). Dans ce dernier cas, on obtient une expression du taux de restitution d'énergie critique comparable à un critère de Griffith à respecter lors de la propagation.

$$D^{(1)}(\underline{\theta}) = \int_{\Gamma_0} (G_c(a) \underline{\theta} \cdot \underline{n}) dS \quad [3]$$

$$D^{(1)}(\underline{\theta}) = \int_{\Gamma_0} \left( \left( \frac{\Delta(a \cdot \underline{n})}{C \Delta N} \right)^{\frac{1}{m}} \underline{\theta} \cdot \underline{n} \right) dS \quad [4]$$

### 3. Calcul du taux de restitution de l'énergie par la méthode θ

Nous avons étendu le calcul par la méthode θ (Mialon, 1988) des dérivées premières et secondes de l'énergie potentielle à l'équilibre au cas d'un milieu élastique présentant des contraintes initiales (pour un état initial en équilibre en l'absence de chargement extérieur). Nous considérons un solide tridimensionnel élastique fissuré occupant le domaine $\Omega \in \mathbb{R}^3$, soumis à un champ de forces surfaciques appliquées sur une partie $S$ de $\partial\Omega$, et d'un champ de forces volumiques appliquées sur $\Omega$. Le champ de contraintes initiales $\sigma^0$ pour le système est en équilibre avec des forces de masse nulles et des efforts surfaciques nuls sur tout le contour $\partial\Omega$ (champ d'autocontrainte).

L'équation de comportement élastique linéarisé du matériau s'écrit :

$$\sigma = \sigma^0 + C:\epsilon \quad [5]$$

Le tenseur des déformations sous contraintes nulle $\epsilon^0$ est lié au tenseur $\sigma^0$ de précontrainte par :

$$\epsilon^0 = -S:\sigma^0 \quad [6]$$

Ainsi, l'énergie potentielle d'un champ $u$ cinématiquement admissible s'écrit :

$$W = \int_{\Omega} (\sigma^0 : \nabla u + \frac{1}{2} \nabla u : C : \nabla u - f \cdot u) d\Omega - \int_{S} g \cdot u \, dS \quad [7]$$

Par définition, le taux de restitution de l'énergie $G$ est :

$$G = - \frac{\partial W}{\partial \Omega} \quad [8]$$
Nous faisons l'hypothèse que les champs de forces volumiques et surfaciques, ainsi que le champ de contraintes résiduelles, ne sont pas modifiés par la propagation de la fissure. Cela se traduit par :

$$f' = \nabla f \cdot \theta \quad ; \quad g' = \nabla g \cdot \theta = 0 \text{ sur } S \quad ; \quad \dot{\sigma}^0 = \nabla \sigma^0 \cdot \theta \qquad [9]$$

où [9] tient compte des hypothèses [11] sur $\theta$.

Le taux de restitution de l'énergie $G(s)$ en chaque point du front de fissure est solution de l'équation variationnelle suivante :

$$\begin{align*}
\int_{\Gamma_0} G(s) \Theta(s) m(s) ds &= -J^{(1)}(\underline{\theta}) \\
&= -\int_{\Omega} \left[ \frac{1}{2} (\nabla u : C : \nabla u) \operatorname{div} \theta + (\sigma^0 : \nabla u) \operatorname{div} \theta - \sigma : (\nabla u : \nabla \theta) \right. \\
&\qquad \left. - (\nabla f \cdot \theta) \cdot u - (f \cdot u) \operatorname{div} \theta + (\nabla \sigma^0 \cdot \theta) : \nabla u \right] d\Omega \quad \forall \theta \in \Theta \quad [10]
\end{align*}$$

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

$$\theta \in \Theta = \{\mu \text{ tels que } \mu.n = 0 \text{ sur } \partial\Omega\} ; \quad \theta_{S \cup S_d} = 0 ; \quad \nabla \theta_{S \cup S_d} = 0 \quad [11]$$

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

Nous avons également établi dans le même cadre la dérivée seconde de l'énergie, nécessaire dans l'équation [2] pour déterminer l'incrément d'avancée du front de fissure. Son expression étant très longue, elle n'est pas donnée ici.

## 4. Mise en œuvre numérique et validation

Pour valider la démarche, nous avons étudié un problème de thermo-élasticité, en traitant les effets induits par un chargement thermique sous forme de contraintes résiduelles. On considère une sphère infinie possédant une fissure circulaire sur les lèvres de laquelle une température est imposée. Une solution analytique de ce problème a été obtenue par (Srivastava et al., 1970) et reprise dans la compilation (Murakami, 1987). Le milieu est modélisé par un matériau de coefficient de dilatation thermique $\alpha = 1,2 \cdot 10^{-5} K^{-1}$, au comportement linéaire élastique isotrope ($E = 200 \text{ GPa}$, $\nu = 0,3$). La température est imposée à 173 K sur les lèvres et maintenue à 273 K à l'infini. Pour l'étude de la propagation de fissure, le chargement est imposé de façon cyclique et on suppose que l'évolution suit une loi de Paris. Le maillage 3D, présenté sur la figure 1, possède une dimension 100 fois plus grande que la longueur de fissure initiale $r$. Il contient 3 528 éléments, avec des fonctions d'interpolations quadratiques, soit 16 161 degrés de libertés en thermique et 48 483 degrés de liberté en mécanique.

Pour mener à bien la résolution de ce problème dans ZéBuLoN, un utilise une approche de type couplage faible entre les phénomènes de thermique et de mécanique :
<figure><img src="image_2.png" /><figcaption>**Figure 1.** Maillage d'un quart de sphère fissuré et visualisation du champ de température à proximité du front.</figcaption></figure>

pour chaque configuration du maillage, (i) le champ de température est obtenu en calculant la solution stationnaire discrétisée du problème de thermique, puis, (ii) le problème de mécanique est résolu à l'aide du champ thermique et le taux de restitution d'énergie est calculé. La qualité du champ thermique ainsi obtenu ne nous a pas permis de prendre en considération le terme en $\nabla\sigma^0$ pour le calcul du taux de restitution d'énergie (les intégrations volumiques devenant trop dépendantes du rayon du volume associé au support de la fonction d'extension $\theta$).

<figure><img src="image_3.png" /><figcaption>**Figure 2.** Représentation du taux de restitution d'énergie pour une entaille unitaire, à gauche, et de l'évolution de rayon de la fissure en propagation sous-critique, à droite.</figcaption></figure>

La première validation repose sur la comparaison de la valeur du taux de restitution d'énergie $G$ obtenue par la méthode $G - \theta$ à l'aide de ZéBuLoN avec la valeur analytique (cf. figure 2). La méthode conduit à une erreur maximale de 9, 8% en utilisant un maillage à l'aide d'éléments $P1$ et de 1, 0% à l'aide d'éléments $P2$ (le volume d'intégration pour la méthode $\theta$ étant torique de rayon $0, 25 r$). Ainsi, on peut vérifier que le terme en $\nabla\sigma^0$ est bien négligeable sur cet exemple.

La figure 2 illustre aussi les résultats obtenus en propagation sous-critique, à l'aide de résolutions analytique et éléments finis. On compare une approche explicite classique et une approche semi-implicite : la résolution thermique n'est pas menée à
chaque itération de la résolution de [2], mais seulement à chaque nouveau pas de
temps. La convergence est obtenue en 4 ou 5 itérations (avec un critère d'arrêt rela-
tif de $10^{-5}$ sur le déplacement du front). Un remaillage est mené à chaque avancée
du front mais aucun transfert de champs n'est requis pour ce problème purement élas-
tique. Ce test montre que la méthode implicite permet d'obtenir un résultat plus proche
de la solution de référence par rapport à la résolution explicite menée avec une largeur
d'incrément de chargement identique ($\Delta N = 5\,000$ cycles).

**5. Conclusions et perspectives**

Une nouvelle écriture des dérivées de l'énergie potentielle par la méthode $\theta$ pre-
nant en compte un champ de contraintes résiduelles a été proposée pour traiter de
manière implicite la propagation de fissure 3D en fatigue. L'expression de la dérivée
du champ de contraintes résiduelles est obtenue en considérant son invariance dans
la propagation de fissure. Néanmoins, même si ces premiers résultats numériques ne
le montrent pas, il peut être important de ne pas négliger dans le calcul numérique le
terme contenant la dérivée lagrangienne du champ de contraintes résiduelles. Ce der-
nier est difficile à calculer avec précision sans une procédure spéciale, que l'on prévoit
de mettre au point dans la suite de ce travail.

**Remerciements.** Ce travail est réalisé et financé dans le cadre du projet MAIA PRC-DDV SAFRAN-CNRS-ONERA-DGA.

**6. Bibliographie**

Destuynder P., Djaoua M., Lescure S., « Quelques remarques sur la mécanique de la rupture élastique », *Journal de Mécanique Théorique et Appliquée*, vol. 2, n° 1, p. 113-135, 1983.

Jung O., Influence du grenaillage sur la propagation de fissures, Thèse de doctorat, ENSMA-Poitiers, 1997.

Mialon P., Calcul de la dérivée d'une grandeur par rapport à un fond de fissure par la méthode $\theta$, Bulletin de la direction des études et recherche n° Série C, EDF, 1988.

Murakami Y., *Stress Intensity Factors Handbook*, Pergamon Press, Japan, 1987.

Ousset Y., « Numerical simulation of delamination growth in layered composite plates », *European Journal of Mechanics*, vol. A/Solids, n° 18, p. 291-312, 1999.

Parks D. M., « A stiffness derivative finite element technique for determination of crack tip stress intensity factors », *International Journal of Fracture*, vol. 10, n° 4, p. 487-502, 1974.

Rice J. R., « A path independent integral and the approximate analysis of strain concentrations by notches and cracks », *Journal of Applied Mechanics*, vol. 35, p. 379-386, 1968.

Srivastava K. N., Dwivedi J. P., « Thermal Stresses in an Elastic Sphere Containing a Penny-shaped crack », *Journal of Applied Mathematics and Physics*, vol. 21, p. 864-886, 1970.