# Nouvelle méthode semi-analytique pour exprimer les champs auxiliaires des intégrales d'interactions

Gwenaël Edeline, Marc Bonnet

► To cite this version:

Gwenaël Edeline, Marc Bonnet. Nouvelle méthode semi-analytique pour exprimer les champs auxiliaires des intégrales d'interactions. CFM 2007 - 18ème Congrès Français de Mécanique, Aug 2007, Grenoble, France. hal-03361998

HAL Id: hal-03361998

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

Submitted on 1 Oct 2021

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
# Nouvelle méthode semi-analytique pour exprimer les champs auxiliaires des intégrales d'interactions

Gwenaël EDELINE, Marc BONNET

CNRS – Ecole Polytechnique
Laboratoire de Mécanique des Solides
Ecole Polytechnique, 91128 Palaiseau Cedex
edeline@lms.polytechnique.fr

## Résumé :

Les approches globales telles que la méthode G-θ présentent l'inconvénient de mélanger les modes de rupture qu'il est nécessaire de différencier pour l'application de critères de branchements (tels que le Principe de Symétrie Locale) dans la simulation de propagation de fissure. L'extraction des modes de rupture est effectuée par l'intégrale d'interaction qui exploite le terme bilinéaire de la formulation G-θ appliqué à la superposition de deux champs solutions ; l'un du problème considéré ; l'autre d'un problème connu judicieusement choisi (champ auxiliaire). Son succès dépend de la validité de ces champs auxiliaires. Comme il n'existe pas de solutions analytiques pour des fissures 3D quelconques, il est habituel de prendre comme champs auxiliaires les solutions de Westergard en chaque point du front de fissure. Nous proposons une nouvelle méthode de construction de ces champs auxiliaires (permettant de prendre en compte les effets de courbures du front) à l'aide de représentations intégrales.

## Abstract :

The global approaches, such as the G-θ method, present the disadvantage of mixing the modes of fracture. This is necessary to differentiate for the application of kinking criteria (such as the "Principle of local symmetry") in order to have the ability to numerically simulate the propagation of a three-dimensional crack. The extraction of the modes of fracture is carried out by the integral of interaction which exploits the bilinear term of the G-θ formulation applied to the superposition of two fields' solutions; the first field is the solution of the problem considered; the second field is the solution of a judiciously selected known problem. Its success depends on the validity of the auxiliary field (the second). As there are not analytical solutions for unspecified 3D cracks, it is usual to take as auxiliary fields the solutions of Westergard in each point of the crack face. We propose a new method of construction of these auxiliary fields (allowing to take into account the effects of a curved crack face) using integral representations.

## Mots-clefs :

### Intégrales d'interactions ; Méthodes G-θ ; Nouveaux champs auxiliaire 3D

## 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 pour la fissuration tridimensionnelle en mode mixte. Les approches globales (taux de restitution de l'énergie G) telles que la méthode G-θ introduite par Destuynder (1983) ont de plus en plus de succès pour leur grande robustesse comparée aux méthodes locales (facteurs d'intensité de contrainte K) sur un même maillage éléments finis classique. Elles présentent néanmoins l'inconvénient de mélanger les modes de rupture. La différenciation de ces modes
est souvent nécessaire pour simuler numériquement la propagation (par fatigue ou rupture brutale) d'une fissure tridimensionnelle chargée en mode mixte ; principalement, pour l'application de critères de branchements et de déversements du front de fissure.

Une méthode très robuste pour découpler les modes de rupture à partir de l'expression du taux de restitution de l'énergie G est l'intégrale d'interaction. Elle consiste en un post-traitement du calcul éléments finis, comme c'est le cas pour le taux de restitution d'énergie calculé par la méthode G-θ. Elle exploite le terme bilinéaire de l'expression de G appliqué à la superposition de deux champs à l'équilibre, l'un solution du problème considéré, l'autre connu et solution d'un problème auxiliaire judicieusement défini. Plus précisément, le champ auxiliaire est construit de façon à avoir une distribution donnée de facteurs d'intensité de contrainte. Le succès de cette méthode dépend de la validité des champs auxiliaires choisis, qui sont astreints à vérifier les équations de l'équilibre en élasticité linéaire.

Dans ce contexte, nous avons réécris la méthode d'extraction des modes de rupture par l'intégrale d'interaction pour une formulation G-θ. Après avoir étudié les différents choix de champs auxiliaires à notre disposition dans le cas d'une fissure tridimensionnelle quelconque et les problèmes qu'ils posent; nous avons proposé une nouvelle méthode semi-analytique de construction de champs auxiliaires tridimensionnelles pour une fissure quelconque, basé sur la théorie des représentations intégrales.

## 2 Écriture et détails de la méthode des intégrales d'interactions

L'intégrale d'interaction permet le découplage des différents modes de rupture à partir de l'écriture du taux de restitution d'énergie G (ou de l'intégrale de Rice J) en exploitant le terme bilinéaire de l'expression de G appliqué à la superposition de deux champs solutions à l'équilibre. L'un $(u, \sigma[u])$ est la solution recherché du problème considéré, l'autre $(v, \sigma[v])$ est une solution connu d'un problème judicieusement défini (solution d'un même problème à l'équilibre, mais pas forcément avec les mêmes conditions aux limites).

Mialon (1988) a explicité la formulation variationnelle G-θ pour un solide tridimensionnelle élastique $\Omega$ (sans forces de volumes) contenant une fissure quelconque $\Gamma$ (de front $\Gamma_0$) permettant d'obtenir par post-traitement du calcul éléments finis le taux de restitution de l'énergie. Nous avons écrit cette formulation en fonction du tenseur d'Eshelby $\underline{E}$ :

$$ \int_{\Gamma_0} G^u(s) \theta(s) \cdot m(s) \, dS = - \int_{\Omega} \nabla \theta \cdot \underline{E} \, dV = g(u, \cdot; \theta) \quad \forall \theta \in \Theta \quad [1] $$

$$ \Theta = \{ \mu \text{ tel que } \mu \cdot n = 0 \text{ sur } \partial \Omega \} $$

Dans cette écriture, les lèvres de la fissure sont libres de contraintes, c'est-à-dire :

$$ \sigma[u] \cdot n = 0 \quad \text{sur } \Gamma \quad [2] $$

Le tenseur d'Eshelby $\underline{E}$ s'exprime de la manière suivante :

$$ \underline{E} = \frac{1}{2} (\sigma[u] : \nabla u) I - \nabla u \cdot \sigma[u] \quad [3] $$
<figure><img src="image_2.png" /><figcaption>FIG. 1 – Support du champ $\theta$ le long du front de fissure.</figcaption></figure>

En écrivant cette même formulation pour les champs u+v et u-v, et en les combinant, nous avons exprimé l’intégrale d’interaction $M_{\theta}$ sous forme bilinéaire symétrique de la formulation variationnelle $G-\theta$, ainsi qu’en fonction des facteurs d’intensité de contrainte des deux champs solutions $(u, \sigma[u])$ et $(v, \sigma[v])$ :

$$ M_{\theta} = \int_{\Gamma_0} \left[ \frac{1-v^2}{E} (K_I^u K_I^v + K_{II}^u K_{II}^v) + \frac{1}{2\mu} K_{III}^u K_{III}^v \right] (s) \theta(s) \cdot m(s) \, dS = g(u,v;\theta) \quad [4] $$

$$ M_{\theta} = \int_{\Omega} \left[ \frac{\partial B}{\partial \nabla u} \cdot (\nabla v \cdot \nabla \theta) + \frac{\partial B}{\partial \nabla v} \cdot (\nabla u \cdot \nabla \theta) - B(u,v) \, \mathrm{div} \, \theta \right] d\Omega $$

dans laquelle $B(u,v)$ est la forme bilinéaire défini par : $B(u,v) = \nabla u \cdot C : \nabla v$

La formule [4] est valable en l'absence de forces de volumes, et pour une fissure libre de contraintes. L'intégrale d'interaction permet alors l'extraction des facteurs d'intensité de contrainte associés à u.

Le champ auxiliaire choisi doit respecter les conditions suivantes, selon le mode que l'on désire extraire :

$$ \begin{align} K_1^v &= 1, K_{II}^v &= 0, K_{III}^v &= 0 && \text{pour extraire } K_I^u \notag \\ K_1^v &= 0, K_{II}^v &= 1, K_{III}^v &= 0 && \text{pour extraire } K_{II}^u \tag{5} \\ K_1^v &= 0, K_{II}^v &= 0, K_{III}^v &= 1 && \text{pour extraire } K_{III}^u \notag \end{align} $$

En insérant dans l'équation [5], les valeurs des facteurs d'intensité de contrainte correspondants aux champs auxiliaires choisis, nous pouvons extraire chaque mode de rupture par les formules suivantes :

$$ \int_{\Gamma_0} K_I(s) \cdot \theta_n \, dS = \frac{M_{\theta}^a \cdot E}{2(1-\nu)} ; \quad \int_{\Gamma_0} K_{II}(s) \cdot \theta_n \, dS = \frac{M_{\theta}^b \cdot E}{2(1-\nu)} ; \quad \int_{\Gamma_0} K_{III}(s) \cdot \theta_n \, dS = M_{\theta}^c \cdot \mu $$
dans lesquelles, les $M_0^a, M_0^b, M_0^c$ sont les valeurs respectives de l'intégrale $M-\Theta$ calculées par la formulation variationnelle bilinéaire symétrique pour chaque champ auxiliaire correspondant au mode à extraire. Le succès de cette méthode dépend de la validité des champs auxiliaires choisis, qui sont astreints à vérifier les équations de l'équilibre en élasticité linéaire.

En 2D, il est vraiment très simple de choisir des champs auxiliaires, puisqu'il suffit de prendre les champs de Westergard. Cependant, en 3D, il n'existe pas de champ analytique pour des fissures de géométrie arbitraire. Pour résoudre ce problème, plusieurs méthodes de construction de champs auxiliaires ont été proposées. Certaines [Bui and al. (1985), Wanlters and al. (2002), Gosz and al. (2002)] utilisent en chaque point du front de fissure les solutions de Westergard, dans un système de coordonnées curvilignes spécifiques, en rajoutant ou non des termes prenant en compte les courbures du front et des lèvres de fissure. D'autres tels que, Kim and al. (2001), utilisent les solutions analytiques pour une fissure circulaire tridimensionnelle en s'adaptant à la courbure de chaque point du front en chaque point du front. Toutes ces approches présentent un caractère approché, les difficultés étant liées à la courbure du front de fissure ou de la fissure elle-même. Par exemple, les solutions de Westergard écrites en coordonnées curvilignes ne vérifient plus les équations d'équilibre et de compatibilité.

Nous proposons alors une nouvelle méthode de construction de ces champs auxiliaires pour des fissures tridimensionnelles quelconques (qui prennent en compte les différentes courbures du front de fissure, et qui respectent néanmoins les équations de l'équilibre en élasticité linéaire).

### 3 Nouvelle méthode de construction de champs auxiliaires tridimensionnels

Dans ce travail, nous proposons une nouvelle méthode de construction de champs auxiliaires pour l'extraction des facteurs d'intensité de contrainte. Cette méthode garantit (aux erreurs de discrétisation près) la satisfaction des équations de champ (équilibre, compatibilité) pour une géométrie arbitraire de fissure tridimensionnelle. Le principe consiste à poser les champs auxiliaires $v(x)$ sous la forme d'une formule de représentation intégrale d'un champ de déplacement, pour une fissure tridimensionnelle quelconque dans un milieu infini, écrite en termes d'un saut de déplacement $\phi_i$ à travers la fissure choisi de manière à vérifier l'une des conditions [5] :

$$v(x) = \int_{\Gamma} \phi_i(y) \cdot T_i^k(x, y) \, dS_y \quad x \in D \quad [6]$$

dans laquelle, $T_i^k$ est le vecteur contrainte de Kelvin, défini par Bonnet (1995).

Tout champs de déplacement auxiliaire choisis sous cette forme vérifie par construction les équations de l'équilibre en élasticité linéaire (sans forces de volume), quelle que soit la géométrie de la fissure.

Le choix d'une représentation intégrale pour un milieu infini permet de diminuer la taille du problème sans être restrictif puisque le champ auxiliaire est solution à l'équilibre d'un même problème que le problème recherché mais n'ayant pas forcément les mêmes conditions aux limites.

D'autre part, le support $\hat{S} \subset S$ de ce saut de déplacement est limité à une petite région autour du point du front de fissure auquel l'évaluation des facteurs d'intensité de contrainte est souhaité, de façon à minimiser le volume des calculs nécessités par ce post-traitement.

Le choix de la densité du saut de déplacement à travers la fissure permet d'ajuster les distributions de facteurs d'intensité de contrainte souhaitées selon le mode que l'on veut
extraire. En effet, la répartition des facteurs d'intensité de contrainte étant connue sur les lèvres de la fissure [5], on peut alors déterminer une interpolation aux valeurs nodales du saut de déplacement à travers la fissure par la formulation d'extrapolation cinématique (formule reliant directement le saut de déplacement à travers la fissure à la distribution de facteurs d'intensité de contrainte sur celle-ci) écrite par Bonnet (1995) afin d'obtenir la distribution de facteurs d'intensité de contrainte souhaitée.

Cette technique ne donne toutefois de bons résultats que dans le cas d'utilisation d'éléments de Barsoum [Bonnet (2006)] permettant en utilisant des fonctions d'interpolation usuelles de représenter la variation en $d^{1/2}$ (distance au front de fissure) de la répartition du saut de déplacement au voisinage du front de fissure.

Ce choix de représentation intégrale du déplacement développe des contraintes sur les lèvres de la fissure (la seul frontière du domaine est la fissure, l'équilibre élastique du domaine créé par des efforts de contact sur les lèvres de la fissure est donc considéré) ; ces efforts (symétriques) exercés sur les lèvres de la fissure sont reliés par Bonnet (1995) à la répartition du saut de déplacement le long des lèvres de la fissure et s'expriment par :

$$ T^{\pm}[v] = (\sigma[v](x).n(x))^{\pm} = -C_{ijkl}n_j^{\pm}(x) \left\{ \int_{\Gamma^+} [D_{lb} \phi_a(y) - D_{lb} \phi_a(x)] \Sigma_{ab}^k(x, y) dS_y \right\} \quad [7] $$

dans laquelle $D_{lb}$ est un opérateur différentiel tangentiel faisant référence à la normale $n^-$ ; $A_{ab}^k$ et $\Sigma_{ab}^k$ sont des tenseurs définis par la solution de Kelvin pour un milieu infini [Bonnet (1995)].

La forme bilinéaire symétrique de l'intégrale d'interaction écrite à partir de la formulation G-θ écrite précédemment en [5] ne prend pas en compte les efforts exercés sur les lèvres de la fissure. Il est donc nécessaire d'ajuster la formule de l'intégrale d'interaction pour prendre en compte les efforts développés sur les lèvres de la fissure par le champ auxiliaire (sans force de volume).

En réécrivant la formulation θ avec la nouvelle hypothèse que pour le champ auxiliaire v la fissure n'est pas libre de contrainte mais soumis à un effort T ; et en suivant la même méthodologie précédente d'écriture de l'intégrale d'interaction, nous obtenons la nouvelle écriture de M-θ sous forme bilinéaire de la formulation variationnelle G-θ, ainsi qu'en fonction des facteurs d'intensité de contrainte des deux champs (u, σ[u]) et (v, σ[v]) :

$$ M_{\theta} = \int_{\Gamma_0} \left[ \frac{1-v^2}{E} (K_I^u K_I^v + K_{II}^u K_{II}^v) + \frac{1}{2\mu} K_{III}^u K_{III}^v \right] (s) \theta(s) \cdot m(s) ds \quad [8] $$

$$ M_{\theta} = \int_{\Omega} \left[ \frac{\partial B}{\partial \nabla u} \cdot (\nabla v \cdot \nabla \theta) + \frac{\partial B}{\partial \nabla v} \cdot (\nabla u \cdot \nabla \theta) - B(u, v) \right] \mathrm{div} \theta \, d\Omega + \int_{\Gamma} \left[ |\nabla_s u \cdot \theta| \right] T^{-}[v] \, dS $$

La mise en œuvre numérique nécessite :

- le calcul du vecteur contrainte associé à [5] sur la portion $\tilde{S} \subset S$ où le champ de vitesse de transformation $\theta$ est non nul ;
- le calcul d'une nouvelle contribution à l'intégrale d'interaction ayant la forme d'une intégrale de surface sur $\tilde{S}$ ;

- le calcul de $\nabla v(x)$ aux points d'intégration du support $\tilde{\Omega} \subset \Omega$ du champ. Les régions $\tilde{S}$ et $\tilde{\Omega}$ étant locales (voisinage surfacique et volumique du point du front de fissure) cette procédure est économique du point de vue du calcul numérique.

## 4 Conclusions et perspectives

Une nouvelle méthode de construction de champs auxiliaires pour des fissures tridimensionnelles quelconques a été proposée. Elle a nécessité la réécriture de l'intégrale d'interaction permettant ainsi d'extraire les facteurs d'intensité de contrainte pour chaque mode de rupture, à partir d'une écriture globale $G-\theta$ du taux de restitution de l'énergie pour une fissure tridimensionnelle quelconque chargée en mode mixte.

Cet algorithme consiste en un post-traitement du calcul éléments finis, comme c'est le cas pour le taux de restitution de l'énergie calculé par la méthode G-θ.

Nous présenterons lors du colloque des simulations numériques de validations (en précision et stabilité) par comparaison aux solutions analytiques d'une fissure circulaire tridimensionnelle.

## Références

- Destyunder, Ph., Djaoua, M., Lescure, S. 1983 Quelques remarques sur la mécanique de la rupture élastique. *Journal de Mécanique Théorique et Appliquée*. 2. n°1, 113-135

- Mialon, P. 1988 Calcul de la dérivée d'une grandeur par rapport à un fond de fissure par la méthode $\theta$. *E.D.F. Bulletin de la direction des études et recherches* – Série C.

- Bui, H. D., Proix, J. M. 1985 Découplage des modes mixtes de rupture en thermo-élasticité par des intégrales indépendantes du contour. *Actes du Troisième Colloque Tendances Actuelles en Calcul de Strustures*. Novembre, 631-643

- Walters M. C., and al. 2005 Interaction integral procedures for 3D curved cracks including surface tractions. *Engineering Fracture Mechanics*. 72, 1635-1663

- Gosz, M., Moran, B. 2002 An interaction energy integral method for computation of mixed-mode stress intensity factors along non-planar crack fronts in three dimensions. *Engineering Fracture Mechanics*. 69, 299-319

- Kim, Y. J., and al. 2001 Mode decomposition of three-dimensional mixed-mode cracks via two-state integrals. *International Journal of Solids and Structures*. 38, 6405-6426

- Bonnet, M. 1995 Equations intégrales et éléments de frontière : Applications en mécanique des solides et des fluides. CNRS Editions/Eyrolles.

- Bonnet, M., Frangi, A. 2006 Analyse des solides déformables par la méthode des éléments finis. Les éditions de l'école Polytechnique.