# Méthode de calcul des facteurs d'intensité des contraintes en mode mixte pour des matériaux à gradient de propriétés mécaniques

Thomas Menouillard, Thomas Elguedj, Alain Combescure

► To cite this version:

Thomas Menouillard, Thomas Elguedj, Alain Combescure. Méthode de calcul des facteurs d'intensité des contraintes en mode mixte pour des matériaux à gradient de propriétés mécaniques. 7e colloque national en calcul des structures, CSMA, May 2005, Giens, France. hal-01813047

HAL Id: hal-01813047

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

Submitted on 12 Jun 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>
# Méthode de calcul des facteurs d'intensité des contraintes en mode mixte pour des matériaux à gradient de propriétés mécaniques

Thomas Menouillard — Thomas Elguedj — Alain Combescure

INSA/LaMCoS, UMR/CNRS 5514
Bât. J. d'Alembert
20 avenue Albert Einstein
69621 Villeurbanne Cedex
Alain.Combescure@insa-lyon.fr

RÉSUMÉ. Nous présentons dans cet article une méthode générale pour le calcul des différents facteurs d'intensité des contraintes dans des matériaux dont la loi de comportement est élastique linéaire continûment variable dans l'espace. La méthode qui permet de prédire les facteurs d'intensité des contraintes est une généralisation de la méthode de l'intégrale d'interaction [KIM 03a] [KIM 03b] [RAO 03]. Une méthode d'obtention systématique des champs asymptotiques en déplacement pour ce type de matériau est aussi développée, et servira à montrer une plus grande qualité des résultats. Les applications numériques se focaliseront sur des matériaux dont le module d'Young dépend de l'espace.

ABSTRACT. In this paper, we present a general method for the calculation of the various stress intensity factors in a material whose constitutive law is elastic, linear and varies continuously in space. The approach used to predict the stress intensity factors is an extension of the interaction integral method [KIM 03a] [KIM 03b] [RAO 03]. For this type of material, we also develop a systematic method to derive the asymptotic displacement fields and use it to achieve better-quality results. Numerical examples will focus on materials with space-dependent Young modulus.

MOTS-CLÉS : Taux de restitution d'énergie, matériaux variables, FGM, facteur d'intensité des contraintes, mode mixte de rupture, champs asymptotiques en pointe de fissure, X-FEM.

KEYWORDS: Energy release rate, non homogeneous material, functionally graded materials, stress intensity factor, mixed-mode fracture, crack tip asymptotic fields, X-FEM.
# 1. Champ de déplacement en pointe de fissure

## 1.1. Forme du champ de déplacement

A partir de la forme du champ de contrainte $\sigma$ obtenu en résolvant $div\sigma = 0$, la relation de comportement ($\mathbb{C}$ linéaire et isotrope) permet de mettre en place la forme du champ de déplacement $U$ en pointe de fissure :

$$ \sigma = K_I g_I + K_{II} g_{II} \qquad (1) $$

$$ \epsilon(U) = \mathbb{C}^{-1} \sigma \qquad (2) $$

Les résultats concernant les champs asymptotiques (notés $u_i^h$) pour les matériaux aux caractéristiques mécaniques constantes dans l'espace (loi de comportement notée $\mathbb{C}^h$) seront un point de départ pour déterminer les champs de matériaux variables dans l'espace (loi de comportement notée $\mathbb{C}^*$, champs asymptotiques notés $U_i^*$). On peut mettre en relation les déformations et lois de comportement de ces deux matériaux :

$$ \text{pour } i = I \text{ et } II : \quad \epsilon(U_i^*) = \mathbb{C}^{*-1} \cdot (\mathbb{C}^h \cdot \epsilon(u_i^h)) $$

On note que le champ de déplacement auxiliaire $U_i^*$ n'est pas compatible avec le champ de contrainte. On applique ce résultat au cas d'un matériau dont le coefficient de Poisson est constant, et le module d'Young variable ; pour i=I et II :

$$ \epsilon \left( \vec{U}_i^h \right) = \frac{E_0}{E(M)} \epsilon \left( \vec{u}_i^h \right) \qquad (3) $$

Les solutions asymptotiques en déplacement en pointe de fissure, pour le cas d'un matériau constant, sont totalement connues. Mais peut-on en déduire les champs asymptotiques des matériaux dont le module d'Young est continûment variable dans l'espace ?

## 1.2. Champs asymptotiques en déplacement : cas d'un matériau au module d'Young variable dans l'espace $E(M) = E_0 e^{\alpha r}$, et au coefficient de Poisson constant

$$ U_{Ir}^*(r, \theta) = \frac{1 + \nu}{E_0} \sqrt{\frac{r}{2\pi}} (k - \cos\theta) \cos\frac{\theta}{2} \frac{1}{\sqrt{\alpha r}} \int_{0}^{\sqrt{\alpha r}} e^{-t^2} dt \qquad (4) $$

$$ \begin{aligned} U_{I\theta}^*(r, \theta) = \frac{1 + \nu}{6E_0} \sqrt{\frac{r}{2\pi}} \left[ \left( (6k - 9)\sin\frac{\theta}{2} + \sin\frac{3\theta}{2} \right) e^{-\alpha r} \right. \\ \left. + \left( (6 - 12k)\sin\frac{\theta}{2} + 2\sin\frac{3\theta}{2} \right) \frac{1}{\sqrt{\alpha r}} \int_{0}^{\sqrt{\alpha r}} e^{-t^2} dt \right] \end{aligned} \qquad (5) $$
$$U_{IIr}^*(r, \theta) = \frac{1+\nu}{E_0} \sqrt{\frac{r}{2\pi}} \left[ -\sin\left(\frac{\theta}{2}\right)(k+\cos\theta) + 2\sin\frac{3\theta}{2} \right] \frac{1}{\sqrt{\alpha r}} \int_{0}^{\sqrt{\alpha r}} e^{-t^2} dt \quad (6)$$

$$U_{II\theta}^*(r, \theta) = \frac{1+\nu}{2E_0} \sqrt{\frac{r}{2\pi}} \left[ \left( (2k-3)\cos\frac{\theta}{2} + \cos\frac{3\theta}{2} - 1 - 3k \right) e^{-\alpha r} + \right. \quad (7)$$

$$ \left( (2 - 4k)\cos\frac{\theta}{2} + 2\cos\frac{3\theta}{2} + 3k + 1 - 2(k + 3)\alpha r \right) \frac{1}{\sqrt{\alpha r}} \int_{0}^{\sqrt{\alpha r}} e^{-t^2} dt] $$

Ces solutions étant complètement analytiques, nous pouvons les utiliser dans les codes de calculs sans avoir à les calculer numériquement.

## 2. Formulation du taux de restitution de l'énergie

M. SUO et COMBESCURE [SUO 92b] [SUO 92a] montrent que le taux de restitution de l'énergie s'écrit, en utilisant le champ virtuel $\Theta$ de la façon suivante (où il n'y a pas de thermique ni de force volumique) :

$$G = \int_{\Omega} \mathrm{tr}(\sigma \cdot \nabla U \cdot \nabla \Theta) \cdot d\Omega - \int_{\Omega} w \cdot \mathrm{div}(\Theta) d\Omega - \frac{1}{2} \int_{\Omega} \mathrm{tr} [\nabla C \cdot \Theta \cdot \epsilon(U) \cdot \epsilon(U)] d\Omega \quad (8)$$

## 3. Analyse en mode mixte pour un matériau variable

La méthode présentée s'appuie essentiellement sur la forme des déplacements en pointe de fissure. On a mis en évidence la relation déplacement - facteurs d'intensité des contraintes $K_I$ et $K_{II}$. La méthode suivante va permettre de calculer les facteurs d'intensité des contraintes.

Il faut introduire la forme bilinéaire symétrique $a(...)$, définie comme suit :

$$a(u, v) = \int_{\Omega} \mathrm{tr} [C\epsilon(u)\nabla v\nabla \Theta + C\epsilon(v)\nabla u\nabla \Theta - (\mathrm{div}(\Theta)C + \nabla C\Theta)\epsilon(u)\epsilon(v)] d\Omega$$

En écrivant $a(U, u_i^s) = K_I a(u_I^s, u_i^s) + K_{II} a(u_{II}^s, u_i^s)$ (pour i=I et II) alors on en déduit directement :

$$K_I = \frac{a(u_{II}^s, u_{II}^s) \cdot a(U, u_I^s) - a(u_I^s, u_{II}^s) \cdot a(U, u_{II}^s)}{a(u_I^s, u_I^s) \cdot a(u_{II}^s, u_{II}^s) - a(u_I^s, u_{II}^s)^2} \quad (9)$$

$$K_{II} = \frac{a(u_{I}^{s}, u_{I}^{s}) \cdot a(U, u_{II}^{s}) - a(u_{I}^{s}, u_{II}^{s}) \cdot a(U, u_{I}^{s})}{a(u_{I}^{s}, u_{I}^{s}) \cdot a(u_{II}^{s}, u_{II}^{s}) - a(u_{I}^{s}, u_{II}^{s})^{2}} \quad (10)$$

On précise que $U$ est le champ de déplacement, $u_I^s$ et $u_{II}^s$ les champs auxiliaires (correspondant aux modes I et II) choisis pour le découplage.
<figure><img src="image_2.png" /><figcaption>**Figure 1.** *Validation numérique : géométrie, chargement, module d'Young*</figcaption></figure>

**4. Applications numériques**

Dans un premier temps, la validation de la formulation des facteurs d'intensité des contraintes permettra de s'interesser au choix de $u_i^s$, par rapport à la dépendance des résultats par rapport au rayon du champ $\Theta$. Les champs $U_i^*$ étant plus proche de la réalité plus loin de la pointe de fissure que les champs $u_i^h$, ils vont permettre d'obtenir des résultats plus indépendants de la taille du champ $\Theta$. On va observer cette différence dans un code éléments finis (CAST3M [SUO 92b]) et dans un code (X-FEM) utilisant la méthode des éléments finis étendus (X-FEM eXtended Finite Element Method [MOE 99] [MOE 02]). Il est important de noter que, dans la méthode X-FEM ici, on découple avec les champs $U_i^*$ ou $u_i^h$, et on enrichit avec $u_i^h$.

Les caractéristiques du **premier cas de calcul** sont présentées sur la figure 1. L'article [KIM 03a] sera considéré comme la référence pour valider la méthode explicitée dans les parties précédentes. La forme bilinéaire $a(\cdot, \cdot)$ est la partie centrale de la validation. Les SIFs sont définis de la façon suivante : $SIF = \frac{K_I}{\overline{\epsilon E}(-0.4)\sqrt{\pi a}}$

Comment choisir les champs auxiliaires qui serviront au calcul de $K_I$. Choisissons d'utiliser les champs asymptotiques d'un matériau dont les caractéristiques mécaniques sont constantes (valeurs en pointe de fissure). Ce seront donc des champs approchés. On va s'interesser à l'influence du rayon du champ $\Theta$ sur les résultats.

La figure 2 présente les calculs (points noirs) des SIFs en fonction de la couronne du champ $\Theta$ virtuel utilisé pour différentes variations du module d'Young. Sont aussi représentés par des droites horizontales, les résultats issus des travaux de Kim et Pau lino [KIM 03a]. Ils proposent diverses méthodes qui sont résumées dans le tableau 1 et comparés aux résultats de la figure 2.

L'intérêt de la méthode développée réside sur le choix des solutions asymptotiques que l'on peut utiliser. Enrichir l'information dans ces solutions permettra d'utiliser des couronnes THETA plus grande, et donc des éléments dans la couronne qui soient plus grands aussi.
<figure><img src="image_3.png" /><figcaption>**Figure 2.** Résultats de la validation numérique</figcaption></figure>

<table><caption>**Tableau 1.** Tableau des résultats de l'expérience comparés aux résultats de KIM et PAULINO [KIM 03a]</caption><thead><tr><th>β</th><th>βa</th><th>Direct</th><th>Method I</th><th>Method II</th><th>Eischen</th><th>cette méthode</th></tr></thead><tbody><tr><td>0</td><td>0</td><td>2.109</td><td>2.133</td><td>2.133</td><td>2.112</td><td>2.1118</td></tr><tr><td>5</td><td>2</td><td>2.289</td><td>2.304</td><td>2.348</td><td>2.295</td><td>2.300</td></tr><tr><td>10</td><td>4</td><td>2.549</td><td>2.589</td><td>2.670</td><td>2.571</td><td>2.586</td></tr><tr><td>15</td><td>6</td><td>2.729</td><td>2.769</td><td>2.879</td><td>2.733</td><td>2.765</td></tr><tr><td>50</td><td>20</td><td>3.050</td><td>3.314</td><td>3.579</td><td>3.228</td><td>3.207</td></tr></tbody></table>

Le second cas de calcul diffère du premier par le chargement (qui sollicite l'éprouvette en mode mixte), et par le module d'Young ($E(r, θ) = E_0e^{αr}$). La figure 3 présente les résultats pour les 4 calculs : découplage dans CAST3M/X-FEM avec les champs auxiliaires du matériau constant/variable dans l'espace. Pour la méthode X-FEM utilisée, l'enrichissement se fait avec les solutions du matériau constant. Ce

<figure><img src="image_4.png" /><figcaption>**Figure 3.** Résultats des calculs correspondant à un chargement en mode mixte pour le matériau $E(r, θ) = 2e^{0,4r}$.</figcaption></figure>
graphique montre que l'utilisation des solutions asymptotiques du matériau variable pour le découplage, permet d'obtenir des erreurs plus faibles qu'en utilisant les solutions asymptotiques du matériau constant.

**5. Conclusion et perspectives**

Le découplage des modes de rupture fait intervenir une couronne virtuelle (champ Θ par exemple). L'intérêt pour le calcul numérique serait d'avoir indépendance par rapport au rayon de ce paramètre. C'est effectivement vrai pour le cas d'un matériau constant dans l'espace. En revanche, lorsque le matériau est variable dans l'espace, les caractéristiques matériau ne sont pas les mêmes partout dans la couronne. Donc l'utilisation des champs de déplacement analytiques du cas matériau constant ne permettent pas d'obtenir l'indépendance par rapport au rayon de la couronne. On a bien montré une plus grande stabilité en utilisant la loi de comportement variable, et les champs asymptotiques du matériau variable. Cependant dans la majorité des cas, il n'y a pas de solution asymptotique pour le matériau E(X) (par exemple FGM). Il faudra réfléchir par la suite au calcul numérique des solutions asymptotiques pour des matériaux variables, qui pourra ensuite être implanté dans X-FEM.

**6. Bibliographie**

[KIM 03a] KIM J.-H., PAULINO G. H., « An accurate scheme for mixed-mode fracture analysis of functionally graded materials using the interaction integral and micromechanics models », *International Journal for Numerical Methods in Engineering*, vol. 58, 2003, p. 1457-1497.

[KIM 03b] KIM J.-H., PAULINO G. H., « Mixed mode J-integral formulation and implementation using graded elements for fracture analysis of nonhomogeneous orthotropic materials », *Mechanics of Materials*, vol. 35, 2003, p. 107-128.

[MOE 99] MOES N., DOLBOW J., BELYTSCHKO T., « A finite element method for crack growth without remeshing », *International Journal of Numerical Methods in Engineering*, vol. 46, 1999, p. 131-150.

[MOE 02] MOES N., GRAVOUIL A., BELYTSCHKO T., « Non-planar 3D crack growth by the extended finite element level sets », *International Journal of Numerical Methods in Engineering*, vol. 53, 2002, p. 2549-2568.

[RAO 03] RAO B., RAHMAN S., « Mesh-free analysis of cracks in isotropic functionally graded materials », *Engineering Fracture Mechanics*, vol. 70, 2003, p. 1-27.

[SUO 92a] SUO X., COMBESCURE A., « Energy release rate and integral for any non homogeneous material », *The American Society of Mechanical Engineers*, vol. 233, 1992, p. 173-179.

[SUO 92b] SUO X., COMBESCURE A., « On the application of G(θ) method and its comparison with De Lorenzi's approach », *Nuclear Engineering and Design*, vol. 135, 1992, p. 207-224.