# Simulation de la propagation de fissures dans les solides elastiques en modes mixtes par la methode des equations integrales duales

Hocine Kebir, Jean Marc Roelandt, Jocelyn Gaudin

► To cite this version:

Hocine Kebir, Jean Marc Roelandt, Jocelyn Gaudin. Simulation de la propagation de fissures dans les solides elastiques en modes mixtes par la methode des equations integrales duales. Revue Européenne des Éléments Finis, 2000, 9 (8), pp.893-916. 10.1080/12506559.2000.10511494. hal-01381733

HAL Id: hal-01381733

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

Submitted on 14 Oct 2016

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
Simulation de la propagation de fissures dans
les solides élastiques en modes mixtes par la
méthode des équations intégrales duales

Hocine Kebir* — Jean Marc Roelandt*

Jocelyn Gaudin**

* Laboratoire de Génie Mécanique pour les Matériaux et les Structures - UPRES A CNRS 6066. Université de Technologie de Compiègne, BP 529 - 60206 Compiègne

** AEROSPATIALE MATRA, Centre Commun de recherches Louis-Blériot
12, rue Pasteur. BP 76, 92152 Suresnes Cedex

RÉSUMÉ Nous présentons dans cet article une démarche numérique basée sur la méthode des éléments de frontière permettant de résoudre différents problèmes complexes de propagation de fissures dans les solides élastiques en modes mixtes Le contour de la structure est discrétisé par des éléments quadratiques, sauf en fond de fissure que nous représentons par des éléments singuliers permettant d'y reproduire la singularité des déformations. Les facteurs d'intensité de contraintes sont évalués à partir des déplacements des nœuds de l'élément singulier. Une comparaison est effectuée avec des résultats expérimentaux de la littérature

ABSTRACT The present paper is concerned with the effective numerical implementation of the two dimensional Dual Boundary element method to analyse the mixed-mode crack growth All the boundaries are discretized with discontinuous quadratic boundary elements and the crack-tip is modeled by singular elements that exactly represent the strain field singularity $1/\sqrt{\tau}$. The Stress Intensity Factors can be computed very accurately from the crack opening displacement at collocation points extremely close to the crack tip. The algorithm is applied to several two dimensional examples and the results obtained are in very good agreement with experimental results.

MOTS-CLÉS méthode des équations intégrales, élément singulier, propagation de fissures en mode mixte, coalescence

KEY WORDS Boundary Element Method, singular element, mixed-mode propagation of a crack, coalescence
# 1. Introduction

Le calcul des structures fissurées a pour objet principal la détermination de différentes caractéristiques de la fissure permettant de répondre aux questions essentielles suivantes

- Est-ce que la fissure va se propager ?

- Si oui, est-ce que la propagation est stable ?

- Si la propagation est stable, avec quelle vitesse la fissure va-t-elle se propager ?

- Pour quelle taille de fissure et à quel moment la propagation devient instable ?

Dans le cadre d'une étude de propagation de fissures en fatigue à grand nombre de cycles, il est admis que la mécanique linéaire élastique de la rupture permet de répondre à ces questions de manière satisfaisante. Elle est basée sur le calcul des facteurs d'intensité de contraintes (FIC). Pour les obtenir, il est nécessaire dans le cas général de déterminer les champs des déplacements et de contraintes en tous points du domaine par une méthode numérique. Les méthodes d'équations intégrales de frontière sont d'un usage assez fréquent en mécanique de la rupture. Outre l'avantage du gain d'une dimension d'espace par rapport aux méthodes d'éléments finis, une meilleure précision est obtenue dans l'évaluation des grandeurs locales au voisinage du front de fissure.

Les procédures numériques de prévision des trajets de propagation de fissures, dans le cadre classique de la mécanique linéaire de la rupture fragile, utilisent toutes la même démarche, chaque pas de calcul comprend trois étapes :

1. La géométrie de la structure fissurée et le chargement étant donnés, on effectue le calcul numérique des FIC. $K_I$ et $K_{II}$.

2. On recherche le trajet de propagation, grâce à l'application d'un critère de direction de propagation qui dépend en général du rapport $K_I/K_{II}$.

3. On prolonge la fissure sur une petite distance (suffisamment grande pour éviter les problèmes numériques) suivant la direction précédemment définie.

Les procédures se distinguent cependant par la manière d'exécuter chaque étape la première peut être réalisée en utilisant, soit la méthode des équations intégrales, soit la méthode des éléments finis. La seconde donne lieu à de multiples possibilités, non seulement suivant le critère employé, mais également suivant que l'on se contente de rechercher l'angle de branchement définissant la direction future de propagation ou que l'on désire également connaître la courbure du trajet futur. Enfin, pour la troisième étape, on doit utiliser le remaillage local ou global dans le cas où l'on utilise la méthode des éléments finis pour l'évaluation des facteurs d'intensité de contraintes, la méthode des équations intégrales duales, celle adoptée dans notre étude, ne nécessite pas de remaillage, la prolongation de la fissure se fait simplement en ajoutant deux éléments singuliers, géométriquement confondus.
Dans cet article, on présente la démarche suivie pour faire propager une fissure dans un milieu homogène. La méthode numérique que nous avons adoptée est la méthode des équations intégrales duales, qui a montré une grande efficacité dans l'étude des structures fissurées en milieux fragiles.

## 2. Formulation

Il existe plusieurs types de formulations intégrales pour des solides élastiques fis-surés. Une de ces formulations, appelée méthode de sous-structuration, consiste à dé-composer le domaine en deux sous-domaines fictifs. Cette technique permet l'utilisation des équations intégrales en déplacement mais engendre deux problèmes élastiques couplés à résoudre. Une autre méthode, appelée méthode des discontinuités de déplacement [COR 79] évite ce couplage mais nécessite l'utilisation d'équations intégrales en traction et l'introduction de variables supplémentaires. L'utilisation des fonctions de Green est une autre formulation des équations intégrales, qui a été essentiellement développée par T A Cruse [CRU 71]. Cette méthode permet d'obtenir une équation intégrale faisant intervenir les tensions et les déplacements de la frontière du domaine, à l'exception de ceux des lèvres de la fissure, ce qui la rend applicable seulement pour quelques situations de géométries particulières.

Dans notre étude, on utilise la méthode des équations intégrales duales qui a été développée en 1992 par Portela et Aliabadi [POR 92]. Cette méthode consiste à écrire l'équation intégrale en déplacement sur les points de collocation d'une lèvre de la fissure et à écrire son équation duale (l'équation intégrale en tension) sur les points de collocation de l'autre lèvre de la fissure.

La formulation de la méthode des équations intégrales en élastostatique est basée sur le théorème de réciprocity de Maxwell-Betti [BRE 89]. En l'absence des forces de volumes et pour un point P n'appartenant pas aux deux lèvres de la fissure, l'équation intégrale en déplacement s'écrit :

$$C_{ij}(P)u_j(P) + \int_{\partial\Omega} T_{ij}(P, Q)u_j(Q)ds(Q) = \int_{\partial\Omega} U_{ij}(P, Q)t_j(Q)ds(Q) \quad [1]$$

où $U_{ij}(P, Q)$ et $T_{ij}(P, Q)$ sont les solutions fondamentales de Kelvin, $C_{ij}$ étant un coefficient qui dépend de la discontinuité de la normale en P $u_j(Q)$ et $t_j(Q)$, déplacement et tension du point Q appartenant au contour du domaine, sont les inconnues du problème.

Si le point P est sur les lèvres de la fissure, cette même équation s'écrit :

$$C_{ij}(P^+)u_i(P^+) + (1 - C_{ij}(P^+))u_i(P^-) \int_{\partial\Omega} T_{ij}(P^+, Q)u_j(Q)ds(Q) = \\ \int_{\partial\Omega} U_{ij}(P^+, Q)t_j(Q)ds(Q) \quad [2]$$
où $u_i(P^+)$ et $u_i(P^-)$ sont les déplacements des lèvres de la fissure au point P. Son équation duale s'écrit :

$$ 
\begin{align*}
C_{ij}(P^{+})t_{i}(P^{+}) &+ (1 - C_{ij}(P^{+}))t_{i}(P^{-}) \\
&+ n_{i}(P^{+}) \int_{\partial\Omega} S_{ijk}(P^{+}, Q)u_{k}(Q)ds(Q) \quad [3] \\
&= n_{i}(P^{+}) \int_{\partial\Omega} D_{ijk}(P^{+}, Q)t_{k}(Q)ds(Q)
\end{align*}
$$

### 3. Calcul des facteurs d'intensité de contraintes

Dans notre étude, nous avons calculé les F.I.C. par :

- l'extrapolation des déplacements au voisinage du fond de la fissure.

- l'estimation de l'intégrale de contour de Rice (Intégrale J).

#### 3.1. Méthode d'extrapolation des déplacements

<figure><img src="image_2.png" /><figcaption>**Figure 1.** Repères cartésien et polaire en fond de fissure</figcaption></figure>

Au voisinage du fond de la fissure, dans le cadre de l'élasticité linéaire, les champs de déplacements et de contraintes sont définis par une série infinie découplée en mode I et II [HAR 69]. En ne considérant que le premier terme de la série, le saut des déplacements des lèvres de la fissure pour les points $P^+(r, \theta = \pi)$ et $P^-(r, \theta = -\pi)$ est donné par :

$$ u_2(\theta = \pi) - u_2(\theta = -\pi) = \frac{\kappa + 1}{\mu} K_I \sqrt{\frac{r}{2\pi}} \quad [4] $$

$$ u_1(\theta = \pi) - u_1(\theta = -\pi) = \frac{\kappa + 1}{\mu} K_{II} \sqrt{\frac{r}{2\pi}} \quad [5] $$
avec $\kappa = \frac{3-\nu}{1+\nu}$ en contraintes planes et $k = 3 - 4\nu$ en déformations planes, $\nu$ étant le coefficient de Poisson et $\mu$ le module de cisaillement.

La connaissance des déplacements des lèvres de la fissure nous permet d'évaluer directement les F.I.C. à partir des formules (4) et (5).

<figure><img src="image_3.png" /><figcaption>**Figure 2.** *Eléments de fond de fissure*</figcaption></figure>

Après résolution du problème, les vecteurs déplacements sont connus en tout point de collocation, notamment pour ceux des deux éléments du fond de la fissure. Donc, d'après les équations (4) et (5), on peut écrire :

$$
\left\{
\begin{array}{l}
K_I^{M_2} = (u_2^{M_2^+} - u_2^{M_2^-}) \frac{2\mu}{\kappa+1} \sqrt{\frac{\pi}{L}} \\[1em]
K_{II}^{M_2} = (u_1^{M_2^+} - u_1^{M_2^-}) \frac{2\mu}{\kappa+1} \sqrt{\frac{\pi}{L}}
\end{array}
\right.
\quad
\left\{
\begin{array}{l}
K_I^{M_3} = (u_2^{M_3^+} - u_2^{M_3^-}) \frac{2\mu}{\kappa+1} \sqrt{\frac{3\pi}{5L}} \\[1em]
K_{II}^{M_3} = (u_1^{M_3^+} - u_1^{M_3^-}) \frac{2\mu}{\kappa+1} \sqrt{\frac{3\pi}{5L}}
\end{array}
\right.
\qquad [6]
$$

Si on effectue une extrapolation linéaire des valeurs des facteurs d'intensité de contraintes entre les points $M_3$ et $M_2$ vers le fond de la fissure on aura :

$$
\left\{
\begin{aligned}
K_I &= \left[ 5 \left( u_2^{M_2+} - u_2^{M_2^-} \right) - \frac{3\sqrt{15}}{5} \left( u_2^{M_3+} - u_2^{M_3^-} \right) \right] \frac{\mu}{\kappa+1} \sqrt{\frac{\pi}{L}} \\
K_{II} &= \left[ 5 \left( u_1^{M_2+} - u_1^{M_2^-} \right) - \frac{3\sqrt{15}}{5} \left( u_1^{M_3+} - u_1^{M_3^-} \right) \right] \frac{\mu}{\kappa+1} \sqrt{\frac{\pi}{L}}
\end{aligned}
\right.
\qquad [7]
$$

## 3.2. Intégrale de Rice (intégrale J)

Considérons un repère cartésien défini en fond de fissure comme le montre la figure 3 (on suppose que les lèvres des fissures ne sont pas chargées). En négligeant les forces de volume, Rice a introduit l'intégrale de contour suivante :

$$J(\bar{u}) = \int_{\Gamma} (W n_1 - t_j u_{j,1}) \, d\Gamma \quad [8]$$
<figure><img src="image_4.png" /><figcaption>Figure 3. Contour Γ de l'intégrale J</figcaption></figure>

où W est la densité de l'énergie de déformation, Γ est un contour orienté ouvert entou-
rant le fond de la fissure dont les extrémités se trouvent sur les lèvres de la fissure, t_j
est la composante du vecteur tension suivant la direction x_j, définie le long du contour
(\vec{x}_1 étant l'axe de la fissure). u_j est la composante du vecteur déplacement suivant la
direction x_j.

Bui [BUI 81] a présenté une méthode élégante de décomposition des modes qui
consiste à décomposer l'intégrale J en la somme de deux intégrales :

$$
J = J(\vec{u}^I) + J(\vec{u}^{II}) \quad \text{avec} \quad J(\vec{u}^I) = \frac{K_I^2}{E'}, \quad J(\vec{u}^{II}) = \frac{K_{II}^2}{E'} \quad [9]
$$

où $E'$ est le module d'élasticité ; $E' = E$ en contraintes planes et $E' = E/(1 - \nu^2)$ en déformations planes.

$$
\left\{ \begin{matrix} u_1^I \\ u_2^I \end{matrix} \right\} = \frac{1}{2} \left\{ \begin{matrix} u_1 + u'_1 \\ u_2 - u'_2 \end{matrix} \right\}, \quad \left\{ \begin{matrix} u_1^{II} \\ u_2^{II} \end{matrix} \right\} = \frac{1}{2} \left\{ \begin{matrix} u_1 - u'_1 \\ u_2 + u'_2 \end{matrix} \right\} \quad [10]
$$

et

$$
u_i'(x_1, x_2) = u_i(x_1, -x_2) \quad [11]
$$

4. Un nouvel élément singulier

Bien que la méthode de l'intégrale de Rice donne de très bons résultats pour l'éva-
luation des facteurs d'intensité de contraintes, elle reste moins performante, voire in-
utilisable dans les cas où la fissure est proche d'un bord ou d'une autre fissure (cas
que nous devons traiter pour les problèmes de la coalescence de deux fissures ou la
coalescence d'une fissure et d'un bord) (Figure 4).
<figure><img src="image_5.png" /><figcaption>**Figure 4.** Intégrale de Rice et les problèmes de coalescence</figcaption></figure>

Un calcul direct des facteurs d'intensité de contraintes à partir des déplacements des nœuds de collocation des éléments de fond de fissure est la solution idéale. Une étude de convergence effectuée avec des éléments quadratiques montre que, quelle que soit la discrétisation, les erreurs relatives en fond de fissure restent importantes (> 8 %). Cette difficulté de représentation des déplacements provient des forts gradients existant en fond de fissure, puisque les déplacements varient en $\sqrt{r}$ ($\frac{\partial u_i}{\partial r}(r = 0) = \infty$) alors que la représentation n'est que quadratique. Il apparaît donc nécessaire d'inclure ces gradients dans un élément spécifique.

Pour ce faire, nous avons développé un nouvel élément singulier, qui a pour particularité de représenter exactement le champ des déplacements en fond de fissure [KEB 99]. Cet élément contient trois nœuds de collocation positionnés comme ceux de l'élément quadratique mais dont les fonctions de forme contiennent un terme linéaire et un terme en $\sqrt{r}$:

$$N_i(s) = a_i + b_i s + c_i \sqrt{s} \quad (i = 1, 3) \quad [12]$$

tel que :

$$N_i(s_j) = \delta_{ij} \quad \text{avec} \quad (s_j)_{j=1,3} = \left(\frac{L_m}{6}, \frac{L_m}{2}, \frac{5L_m}{6}\right) \quad [13]$$

En résolvant les systèmes d'équations (12) et (13) on obtient :

$$\left\{ \begin{aligned} N_1(s) &= \frac{1}{\sqrt{10}+3\sqrt{2}-8} \left( 2(\sqrt{10}-5) + \frac{6(2-\sqrt{10})}{L_m} s + \frac{6\sqrt{3}}{\sqrt{L_m}} \sqrt{s} \right) \\ N_2(s) &= \frac{1}{\sqrt{10}+3\sqrt{2}-8} \left( \sqrt{2}(5-\sqrt{5}) + \frac{6\sqrt{2}(\sqrt{5}-1)}{L_m} s - \frac{8\sqrt{3}}{\sqrt{L_m}} \sqrt{s} \right) \\ N_3(s) &= \frac{1}{\sqrt{10}+3\sqrt{2}-8} \left( 2(1-\sqrt{2}) + \frac{6(\sqrt{2}-2)}{L_m} s + \frac{2\sqrt{3}}{\sqrt{L_m}} \sqrt{s} \right) \end{aligned} \right. \quad [14]$$
<figure><img src="image_6.png" /><figcaption>**Figure 5.** Fonctions de forme de l'élément singulier</figcaption></figure>

**5. Direction de propagation de la fissure**

Lorsqu'une fissure est sollicitée en mode mixte, on observe généralement une bifurcation qui se caractérise par un changement de la direction de propagation. Si l'on se contente de la prédiction de l'angle de bifurcation, le seul élément théorique requis est le critère de propagation fournissant cet angle.

Le critère le plus simple à utiliser est celui d'Erđogan et Sih [ERD 63], qui stipule que la direction de propagation est celle suivant laquelle la contrainte tangentielle $\sigma_{\theta\theta}$ est maximale (le cisaillement $\sigma_{r\theta}$ est nul).

<figure><img src="image_7.png" /><figcaption>**Figure 6.** Repères cartésien et polaire en fond de fissure</figcaption></figure>
Au voisinage du fond de fissure (figure 6), les champs de contraintes sont donnés en coordonnées polaires par :

$$
\begin{cases}
\sigma_{\theta\theta}(M) = \frac{1}{\sqrt{2\pi r}} \cos \frac{\theta}{2} [K_I \cos^2 \frac{\theta}{2} - \frac{3}{2} K_{II} \sin \theta] \\
\sigma_{r\theta}(M) = \frac{1}{\sqrt{2\pi r}} \sin \theta [K_I \sin \theta + K_{II}(3 \cos \theta - 1)]
\end{cases}
\quad [15]
$$

L'angle de bifurcation θt correspond au maximum de la contrainte tangentielle σθθ, il
vérifie

$$
\tan \frac{\theta_t}{2} = \frac{1}{4} \left[ \frac{K_I}{K_{II}} \pm \sqrt{\left(\frac{K_I}{K_{II}}\right)^2 + 8} \right] \quad \text{avec} \quad -\frac{\pi}{2} < \theta_t < \frac{\pi}{2} \qquad [16]
$$

6. Pilotage de la fissure

Il existe deux manières différentes pour le pilotage de la fissure

- la première est le pilotage en longueur de fissure; l'évolution de la fissure s'effectue par pas de longueur (Δa)ᵢ, ce pas peut être constant ou variable durant la propagation. Le nombre de cycles nécessaires pour faire propager la fissure de (Δa)ᵢ se calcule alors à partir de sa vitesse de propagation.

$$
(\Delta N)_i = \frac{(\Delta a)_i}{\|\vec{v}(\Delta K_{eff})\|} \quad [17]
$$

La vitesse de propagation $\vec{v}(\Delta K_{eff})$ de la fissure est obtenue à partir de la connaissance du FIC. effectif $K_{eff} = \sqrt{K_I^2 + 2 K_{II}^2}$ et de la loi de fissuration (loi de Paris, loi de Forman, ou loi expérimentale (point par point)) choisie en fonction de son domaine de validité,

- la seconde est le pilotage en nombre de cycles, l'évolution de la fissure est suivie par saut de $(\Delta N)_i$ cycles. La taille de l'incrément de propagation $(\Delta a)_i$ durant les $(\Delta N)_i$ cycles sera déterminée par

$$
(\Delta a)_i = \| \vec{v}(\Delta K_{eff}) \| (\Delta N)_i \quad [18]
$$

Le pilotage en nombre de cycles crée un problème de mise en œuvre numérique
de la méthode. En effet, la taille de l'incrément de fissure dépend de la vitesse de
propagation et peut varier énormément d'un pas à l'autre, le maillage n'est alors en
général plus régulier, ce qui induit des imprécisions d'ordre numérique, le remaillage
local est alors nécessaire

7. Phénomène de la coalescence

Il y a coalescence ou rupture de ligament quand deux fissures se rejoignent, ou
quand une fissure rejoint le contour (un trou ou un bord libre) et cela sans que la
rupture totale de la structure ne se produise nécessairement. Prévoir à quel moment se produit cette jonction, sous sollicitation de fatigue, revient à déterminer la ou les longueurs des fissures au moment de la rupture du ligament.

## 7.1. Critères de coalescence fissure/fissure

<figure><img src="image_8.png" /><figcaption>**Figure 7.** Coalescence fissure/fissure selon Swift</figcaption></figure>

Un critère de coalescence permet de prévoir dans quelles conditions deux fissures coalecent (rupture instable du ligament). Le critère le plus couramment utilisé est celui de Swift [JEO 95]. Il stipule que la coalescence intervient quand les zones plastiques générées en fond de fissure se rejoignent.

$$ \text{Coalescence des fissures} 2a \text{ et } 2b \Leftrightarrow L = R_p(a) + R_p(b) \quad [19] $$

avec $R_p$ = taille de zone plastique en pointe de fissure

La taille de la zone plastique en pointe de fissure est liée à la valeur du facteur d'intensité de contraintes effectif $K_{eff} = \sqrt{K_1^2 + 2K_{II}^2}$ par la formule d'Irwin.

$$ R_p = \frac{1}{\pi} \left( \frac{K_{eff}}{R_{0.2}} \right)^2 \quad \text{Formule d'Irwin (correction de zone plastique)} \quad [20] $$

où $R_{0.2}$ est la limite d'élasticité du matériau

## 7.2. Critère de coalescence fissure/contour

Le critère de Swift est défini pour prédire la jonction de deux fissures. Aucune approche spécifique n'existe ni pour la coalescence d'une fissure avec un trou ni pour la coalescence d'une fissure avec un bord libre. Nous avons étendu l'approche de Swift à ces deux configurations.

On suppose que la rupture du ligament se produit lorsque la zone plastique en tête de fissure rejoint le contour (figure 8).
<figure><img src="image_9.png" /><figcaption>**Figure 8.** Coalescence fissure/contour</figcaption></figure>

# 8. Mise en œuvre numérique

## 8.1. Discrétisation

Toute la structure est discrétisée par des éléments quadratiques non conformes sauf en fond de fissure qui sont modélisés par des éléments singuliers. Ceci nous a permis de représenter exactement la singularité des déformations en fond de fissure, ce qui par conséquent, nous évite d'y raffiner le maillage et d'utiliser les invariants intégraux tout en obtenant une très bonne précision.

L'algorithme de la méthode peut être résumé comme suit [KEB 96] :

- les contours sont discrétisés par des éléments quadratiques non conformes, sauf en fond de fissure où l'on a utilisé des éléments singuliers ;

- l'équation en déplacement (4) est appliquée aux points de collocation d'une lèvre de la fissure ;

- l'équation en traction (5) est appliquée aux points de collocation de l'autre lèvre de la fissure ;

- l'équation de déplacement (1) est appliquée aux nœuds de collocation situés sur le contour extérieur ;

- les facteurs d'intensité de contraintes sont évalués par la méthode d'extrapolation des déplacements, en exploitant les déplacements des nœuds de l'élément singulier ;

- l'incrémentation de la fissure se fait en ajoutant deux nouveaux éléments singuliers, qui génèreront de nouvelles équations qui seront mises à jour avec celles déjà existantes ;

Au pas *i*, après détermination de la longueur de propagation $(\Delta a)_i$, la mise en œuvre numérique de la propagation se fait en ajoutant deux éléments singuliers géométriquement confondus en fond de fissure (figure 9).
<figure><img src="image_10.png" /></figure>

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

■ Noeud de géométrie

E.Q.N.C. : éléments quadratiques non conformes

• Noeud de collocation

E.S. : éléments singuliers

<caption>**Figure 9.** Incrémentation numérique de la fissure</caption>

## 8.2. Constitution du système d'équations

Avant propagation, un premier calcul élastique nous fournit un système d'équations de la forme :

$$[A_0]\{u_0\} = [B_0]\{t_0\} \qquad [21]$$

$\{t_0\}$ et $\{u_0\}$ désignent respectivement les vecteurs tensions et déplacements aux noeuds de collocation.

A chaque pas de propagation, l'incrémentation de la fissure se fait en ajoutant deux éléments singuliers en fond de fissure, cela se traduit par $2 \times 6$ degrés de liberté supplémentaires. Les tensions sur les lèvres de la fissure étant nulles, seules 6 équations additionnelles sont nécessaires pour la résolution du problème.

Pour les avoir, il suffit d'écrire :

- l'équation intégrale en déplacement pour les trois points de collocation de l'élément singulier $\Gamma_{(NEF+1)}^f$ suivant les deux directions $\vec{e}_i$ ($i = 1, 2$);

- l'équation intégrale en tension pour les trois points de collocation de l'élément singulier $\Gamma_{(NEF+2)}^f$ suivant les deux directions $\vec{e}_i$ ($i = 1, 2$).

D'un point de vue numérique, on ne calcule que les coefficients liés aux nouveaux
<figure><img src="image_12.png" /><figcaption>**Figure 10.** Schéma représentant la constitution du système d'équations</figcaption></figure>

éléments singuliers, qui constituent les 6 dernières colonnes et les 6 dernières lignes des matrices [A] et [B] (figure 10).

D'autre part, les éléments singuliers d'un pas *i* ne seront plus singuliers au pas *i* + 1, puisqu'ils ne seront plus en fond de fissure. Ils redeviennent alors des éléments quadratiques. Ceci a pour conséquence un nouveau calcul des coefficients des matrices [A] et [B] correspondant aux deux anciens éléments singuliers.

## 8.3. Résolution du système d'équations

La résolution du sytème d'équations est la partie la plus coûteuse en temps de calcul lors d'une simulation par la méthode des équations intégrales. Dans le but d'optimiser le temps de résolution durant la propagation de la fissure, un algorithme spécial est mis en œuvre. Il consiste principalement à sauvegarder les résultats de triangularisation du pas précédent, de calculer la partie du système liée aux nouvelles inconnues et de triangulariser partiellement le nouveau système d'équations.

Cet algorithme nous a permis de réduire considérablement le temps de calcul d'un problème de propagation de fissures. A titre indicatif, le temps de calcul est divisé par 10 pour un calcul de propagation de 20 pas (700 degrés de liberté à la fin de la propagation).

## 9. Validation de l'élément singulier

Pour illustrer la précision de l'élément singulier nous avons considéré une plaque rectangulaire contenant une fissure oblique et centrée (Figure 11).

La plaque est soumise à une traction nominale de $\sigma = 100 \text{ MPa}$. La fissure a une longueur de $2a$ et elle est inclinée à $45^\circ$. La géométrie est définie par les rapports
$h/w = 2$ et $a/h = 0.6$. La discrétisation du contour extérieur utilise 34 éléments quadratiques; on a fait varier le maillage de la fissure, discrétisée par $2N$ éléments de longueurs identiques $L = a/N$. Les résultats obtenus sont présentés dans la table 1

<figure><img src="image_13.png" /><figcaption>Figure 11. Maillage et déformée d'une plaque contenant un fissure oblique</figcaption></figure>

<table><caption>Table 1. Facteurs d'intensité de contraintes en modes I et II</caption><thead><tr><th>N</th><th>1</th><th>2</th><th>3</th><th>4</th><th>6</th><th>10</th><th>Ref [MUR 87]</th></tr></thead><tbody><tr><td>K<sub>I</sub></td><td>0.6445</td><td>0.6566</td><td>0.6587</td><td>0.6595</td><td>0.6603</td><td>0.6609</td><td>0.661</td></tr><tr><td>K<sub>II</sub></td><td>0.5482</td><td>0.5582</td><td>0.5615</td><td>0.5631</td><td>0.5647</td><td>0.5660</td><td>0.567</td></tr></tbody></table>

La figure 12 représente l'écart en pourcentage entre la solution obtenue par l'élément singulier et celle de référence (Murakami [MUR 87]) en fonction du rapport $L/a$ (Longueur de l'élément singulier/Longueur de la fissure).

On constate que les résultats obtenus par l'élément singulier sont très précis même pour un faible nombre d'éléments. En effet, la simple discrétisation de la fissure non débouchante par deux éléments ($L/a = 1$) nous donne une erreur relative inférieure à 3.5%, en mode I et en mode II. De plus, on constate que les résultats ne sont pas très dépendants de la taille des éléments en fond de fissure. En effet, l'écart entre deux résultats obtenus par des éléments singuliers de longueurs $L/a = 0.5$ et $L/a = 0.1$ ne dépasse pas 1%.
<figure><img src="image_14.png" /><figcaption>**Figure 12.** Ecart en % entre les F.I.C. calculés et ceux de Murakami en fonction du rapport (longueur de l'élément singulier / longueur de la fissure)</figcaption></figure>

# 10. Exemples de propagation de fissure

Cette démarche a été implantée dans notre code de calcul nommé PRF2D. Dans ce qui suit, on présente deux exemples de simulation de la propagation de fissures.

## 10.1. Poutre trouée et fissurée

Dans ce premier exemple, on simule la propagation d'une fissure dans une poutre en polyméthylmétacrylate (PMMA) et on procède à une comparaison avec les résultats expérimentaux obtenus par Ingraffea et présentés dans la thèse de Autesserre [AUT 95].

La pièce étudiée (figure 13) présente trois trous circulaires disposés selon une ligne verticale : la fissure est initiée par une entaille verticale de longueur *b* pouvant prendre deux positions selon la valeur de *a*.

Le chargement consiste en une "force ponctuelle" appliquée verticalement de haut en bas au milieu de la pièce, équilibrée au moyen d'appuis simples situés aux extrémités.

L'intérêt de la présence de trous est que ceux-ci modifient considérablement le champ de contraintes et engendrent ainsi des trajectoires de fissuration particulières, en semblant "attirer" les fissures (effet du bord).
<figure><img src="image_15.png" /><figcaption>Figure 13. Géométrie initiale de la poutre étudiée (mm)</figcaption></figure>

Une étude expérimentale a montré que le trajet de propagation dépend fortement de la taille de l'entaille et de sa position relative par rapport à l'axe vertical passant par le centre des trous ; en fonction de ces paramètres, la fissure sera attirée par, et pourra éventuellement déboucher dans tel ou tel trou.

La discrétisation du contour extérieur est réalisée à l'aide de 29 éléments quadratiques non conformes, les trous sont discrétisés en 8 éléments quadratiques non conformes chacun, et l'entaille initiale en 3 éléments quadratiques non conformes sur chaque lèvre. La taille de l'incrément de propagation est de 0.1mm. La fissure arrête automatiquement de se propager dès qu'elle atteint un bord.

<figure><img src="image_16.png" /><figcaption>Figure 14. Trajet numérique et expérimental de la première configuration</figcaption></figure>

Nous présentons tout d'abord les résultats relatifs à la première configuration dans laquelle l'entaille est loin de l'axe vertical des trous ($a = 2mm$ et $b = 1mm$). La figure 14 montre d'une part le trajet exprimental de la fissure, qui débouche sur le
trou central, et d'autre part le trajet prédit par $PRF2D$.

La propagation est faite en 42 pas, ce qui nous donne une longueur de propagation de 4.2mm. Le trajet numérique reproduit bien l'attirance de la fissure vers le deuxième trou.

La seconde configuration correspond à la position la plus proche de l'axe vertical des trous ($a = 1mm$ et $b = 1.5mm$). Les trajets numérique et expérimental sont présentés en figure 15. On voit que la fissure se dirige d'abord vers le premier trou, puis s'en écarte pour finalement déboucher dans le deuxième trou. La propagation est faite en 36 pas, ce qui nous donne une longueur de propagation de 3.6mm.

On note la très bonne corrélation entre le trajet expérimental et le trajet numérique prédit par le code de calcul $PRF2D$ dans les deux configurations étudiées.

<figure><img src="image_17.png" /><figcaption>**Figure 15.** Trajet numérique et expérimental de la seconde configuration</figcaption></figure>

## 10.2. Ecaillage des matériaux fragiles

L'écaillage des matériaux fragiles est un exemple intéressant de propagation de fissure. En effet, quand une plaque fragile est sollicitée le long de son bord supérieur, une fissure s'amorce, se propage parallèlement au bord latéral de la plaque, ce qui peut provoquer le détachement d'une partie de la plaque (écaillage) (figure 16).

Les mécanismes de fissuration et d'écaillage qui accompagnent la mise en charge des plaques fragiles sur la tranche ont été étudiés expérimentalement par M. D. Thouless et al [THO 87]. Les expériences réalisées sur du polyméthylmétacrylate (PMMA dont le comportement est élastique linéaire dans les conditions des expériences) ont révélé des tendances systématiques dans la localisation de la fissure et dans le début de l'écaillage. En particulier, les tests ont permis d'identifier un domaine de croissance stationnaire de la fissure dans lequel les fissures se propagent parallèlement à la surface latérale.

La modélisation numérique de la propagation a été effectuée en discrétisant le contour de la plaque en 50 éléments et la fissure initiale en 5 éléments sur chaque lèvre.
<figure><img src="image_18.png" /></figure>

$$d_0 < d \qquad d_0 = d \qquad d_0 > d$$

<caption>**Figure 16.** Configurations d'écaillages d'éprouvettes en PMMA [THO 87]</caption>

La propagation est effectuée par un pas $\Delta a = 1 \text{ mm}$ sur une plaque rectangulaire de dimension $100 \text{ mm} \times 50 \text{ mm}$. Plusieurs calculs ont été effectués en faisant varier la position initiale de la fissure $d_0$ ainsi que la position du chargement $P$ définie par la distance $h$ (figure 17).

<figure><img src="image_19.png" /><figcaption>**Figure 17.** Maillage et déformée de la plaque après propagation</figcaption></figure>
Les calculs de propagation de fissures ont permis de comparer les trajectoires des fissures avec les mesures expérimentales. La tendance générale à la fissuration est compatible avec le comportement prédit. On notera en particulier les points suivants :

- La fissure se propage parallèlement au bord latéral à une distance *d*, de telle manière que le rapport *d*/*h* soit constant (où *h* est la distance entre le chargement *P* et le bord latéral de la plaque).

- Les différents chemins de propagation en fonction de *h* sont présentés en figure 18, la table 2 donne les valeurs numériques et expérimentales du rapport *d*/*h*. La valeur numérique de *d*/*h* est en très bon accord avec celle issue des expériences. En effet, l'écart maximum avec l'expérience ne dépasse pas 2 %.

- Pour un chargement situé à une distance donnée *h*, la fissure a tendance à re-venir vers sa position de propagation stationnaire *d* quelle que soit la valeur de sa distance initiale *d*<sub>0</sub> (figure 19).

- Quand la fissure arrive à l'état stationnaire de propagation, le facteur d'intensité de contraintes *K*<sub>I</sub> prend une valeur quasi constante (figure 20) alors que *K*<sub>II</sub> est approximativement nul (figure 21).

- Le temps CPU pour l'étude de la propagation de la fissure en 50 pas ne dépasse pas 150 secondes sur une machine DEC Alpha 500 (400 MHz).

<table><caption>**Table 2.** Valeurs numériques et expérimentales de *d*/*h* (*h* en mm)</caption><thead><tr><th></th><th>Essais</th><th>calcul</th><th>calcul</th><th>calcul</th><th>calcul</th></tr><tr><th></th><th>Valeur moyenne</th><th>h = 2.0</th><th>h = 2.75</th><th>h = 3.75</th><th>h = 5.00</th></tr></thead><tbody><tr><th>d/h</th><td>4.8</td><td>4.82</td><td>4.81</td><td>4.76</td><td>4.71</td></tr></tbody></table>

## 11. Conclusions

L'utilisation de la méthode duale des équations intégrales associée à une représentation singulière des déformations au fond de fissure s'avère très performante pour l'évaluation des grandeurs de rupture fragile. Nous avons développé et implanté un algorithme permettant de simuler différents problèmes de propagation de fissures en modes mixtes, situations difficiles à reproduire par la méthode des éléments finis.

Les résultats obtenus montrent que l'approche proposée pour l'évaluation des facteurs d'intensité de contraintes est simple, d'un niveau de précision élevé et n'est pas très dépendante de la taille des éléments en fond de fissure. Les résultats obtenus sont en très bon accord avec les résultats expérimentaux. Une approche probabiliste basée sur la méthode de Monte-Carlo a été développée pour prévoir la durée de vie des assemblages boulonnés en utilisant des simulations de propagation par le modèle présenté.
<figure><img src="image_20.png" /><figcaption>Figure 18. Trajectoires de propagation en fonction de h ($d_0$ = 15 mm)</figcaption></figure>

<figure><img src="image_21.png" /><figcaption>Figure 19. Trajectoires de propagation en fonction de $d_0$ ($h$ = 2 mm)</figcaption></figure>
<figure><img src="image_22.png" /><figcaption>**Figure 20.** F.I.C. $K_I$ (en MPa$\sqrt{m}$) au cours de la propagation
($h$ = 2 mm, $d_0$ = 15 mm, $P$ = 33kN)</figcaption></figure>

<figure><img src="image_23.png" /><figcaption>**Figure 21.** F.I.C. $K_{II}$ (en MPa$\sqrt{m}$) au cours de la propagation
($h$ = 2 mm, $d_0$ = 15 mm, $P$ = 33kN)</figcaption></figure>
## 12. Annexe : Equations intégrales duales pour les structures fissurées

<figure><img src="image_24.png" /><figcaption>Figure 22. Structure fissurée ($P \in \Gamma$)</figcaption></figure>

L'équation intégrale de frontière (1) est valable telle quelle pour un solide élastique non fissuré. La présence d'une fissure entraîne des difficultés supplémentaires. Considérons donc un solide élastique $\Omega$, contenant une fissure $\Gamma$; les deux lèvres $\Gamma^+$, $\Gamma^-$, géométriquement confondues et libres de contraintes, sont munies de normales unitaires $n^+$ et $n^-$, opposées et choisies de sorte que $n^-$ soit orientée de $\Gamma^-$ vers $\Gamma^+$ (figure 22).

### Equation intégrale en déplacement pour les points de la fissure

Décomposons par la pensée $\Omega$ en deux sous-domaines $\Omega^+$ et $\Omega^-$ séparées par un contour $C$ et la fissure $\Gamma$ (figure 22). On peut alors écrire, pour un point $P^+$ appartenant à la lèvre $\Gamma^+$ de la fissure, l'équation intégrale en déplacement (1) :

$$ 
\begin{aligned}
& \frac{1}{2} u_i(P^+) + \int_{\partial\Omega^+ \cup C^+ \cup \Gamma^+} T_j^i(P^+, Q) u_j(Q) \, ds(Q) \\
& = \int_{\partial\Omega^+ \cup C^+ \cup \Gamma^+} U_j^i(P^+, Q) t_j(Q) \, ds(Q)
\end{aligned} 
\quad [22] $$

De la même manière, l'équation intégrale en déplacement (1), pour un point $P^-$ appartenant à la lèvre $\Gamma^-$, s'écrit :

$$ 
\begin{aligned}
& \frac{1}{2} u_i(P^-) + \int_{\partial\Omega^- \cup C^- \cup \Gamma^-} T_j^i(P^-, Q) u_j(Q) \, ds(Q) \\
& = \int_{\partial\Omega^- \cup C^- \cup \Gamma^-} U_j^i(P^-, Q) t_j(Q) \, ds(Q)
\end{aligned} 
\quad [23] $$

Tous calculs faits, l'addition membre à membre des deux égalités (22) et (23), prises pour le même point $P = P^+ = P^-$, donne, compte tenu des conditions de raccord et l'opposition des normales $n^+$ et $n^-$ sur $C$ :

$$ \frac{1}{2} (u_i(P^+) + u_i(P^-)) + \int_{\partial\Omega \cup \Gamma^+ \cup \Gamma^-} T_j^i(P, Q) u_j(Q) \, ds(Q) $$
$$ = \int_{\partial\Omega \cup \Gamma^{+} \cup \Gamma^{-}} U_{j}^{i}(P, Q) t_{j}(Q) \, ds(Q) \quad [24] $$

ou encore, sachant que $\Gamma = \Gamma^{+} \cup \Gamma^{-}$ :

$$ \begin{aligned} & \frac{1}{2} (u_i(P^+) + u_i(P^-)) + \int_{\partial\Omega \cup \Gamma} T_j^i(P, Q) u_j(Q) \, ds(Q) \\ & = \int_{\partial\Omega \cup \Gamma} U_j^i(P, Q) t_j(Q) \, ds(Q) \quad [25] \end{aligned} $$

L'équation (25) représente l'équation intégrale en déplacement pour les points de la fissure.

## Equation intégrale en contrainte

L'équation intégrale en contrainte est calculée à partir de la loi de comportement de Hooke (26) :

$$ \sigma_{jk}^{i} = \lambda \delta_{jk} u_{i,l}^{i} + \mu (u_{j,k}^{i} + u_{k,j}^{i}) \quad [26] $$

En dérivant les termes de l'équation (1) selon la loi de Hooke, on obtient :

$$ \frac{1}{2} \sigma_{ij}(P) + \int_{\partial\Omega} S_{jk}^{i}(P, Q) u_{k}(Q) \, ds(Q) = \int_{\partial\Omega} D_{jk}^{i}(P, Q) t_{k}(Q) \, ds(Q) \quad [27] $$

## Equation intégrale en tension pour les points de la fissure

En suivant le même raisonnement que pour l'équation intégrale en déplacement, l'équation intégrale en contrainte pour un point $P = P^{+} = P^{-}$ situé sur les lèvres de la fissure s'écrit :

$$ \begin{aligned} & \frac{1}{2} (\sigma_{ij}(P^{+}) + \sigma_{ij}(P^{-})) + \int_{\partial\Omega \cup \Gamma} S_{jk}^{i}(P, Q) u_{k}(Q) \, ds(Q) \\ & = \int_{\partial\Omega \cup \Gamma} D_{jk}^{i}(P, Q) t_{k}(Q) \, ds(Q) \quad [28] \end{aligned} $$

La multiplication des deux membres de l'équation (28) par la normale $n^{+}$ du point $P^{+}$ ou $n^{-}$ du point $P^{-}$ nous donne l'équation intégrale en tension :

$$ \begin{aligned} & \frac{1}{2} (t_i(P^+) + t_i(P^-)) + n_j(P^+) \int_{\partial\Omega \cup \Gamma} S_{jk}^i(P, Q) u_k(Q) \, ds(Q) \\ & = n_j(P^+) \int_{\partial\Omega \cup \Gamma} D_{jk}^i(P, Q) t_k(Q) \, ds(Q) \quad [29] \end{aligned} $$
Pour une fissure non chargée l'équation intégrale en tension se réduit à :

$$ 
\begin{aligned}
n_j(P^+) & \int_{\partial\Omega \cup \Gamma} S_{jk}^i(P, Q) u_k(Q) ds(Q) \\
& = n_j(P^+) \int_{\partial\Omega} D_{jk}^i(P, Q) t_k(Q) ds(Q)
\end{aligned} 
\quad [30] $$

## 13. Bibliographie

[AUT 95] AUTESSERRE P., «Simulation numérique de trajets de propagation de fissures en mécanique linéaire bidimensionnelle de la rupture», Thèse de doctorat, Université Pierre et Marie Curie, 1995.

[BRE 89] BREBBIA C. A. et DOMINGUEZ J., Boundary elements an introductory course, Mc Graw-Hill, 1989.

[BUI 81] BUI H., *Advances Fracture Research*, Pergamon Press, 1981.

[COR 79] CORNET F. H., «Comparative Analysis by the Displacement-Discontinuity Method of Two Energy Criteria of Fracture», *J. of applied Mechanics*, vol. 46, 1979, p. 349-355.

[CRU 71] CRUSE T. A. et VANBUREN W., «Three dimensional elastic stress analysis of a fracture specimen with an edge crack», *Int. Jour. Fract. Mech.*, vol. 7, 1971.

[ERD 63] ERDOGAN F. et SIH G. C., «On the crack extension in plates under plane loading and transverse shear», *Journal of basic engineering*, vol. 85, 1963, p. 519-527.

[HAR 69] HARTRANFT R. et SIH G., «The use of eigenfunction expansions in the general solution of three-dimensional crack problems», *J. of Math. and Mech.*, vol. 19, 1969, p. 123-138.

[JEO 95] JEONG D. Y. et BREWER J. C., «On the linkup of multiple cracks», *Engineering Fracture Mechanics*, vol. 51, 1995, p. 233-238.

[KEB 96] KEBIR H., ROELANDT J. et FOULQUIER J., «Dual boundary element incremental analysis of crack growth in bolted joints», In *MAT-TEC 96 Improuvement of materials*, Paris, 1996, I.I.T.T. International, p. 223-230.

[KEB 99] KEBIR H., ROELANDT J. et FOULQUIER J., «new singular boundary element for crack problems Application to bolted joints», *Engineering Fracture Mechanics*, vol. 62, 1999, p. 497-510, PERGAMON.

[MUR 87] MURAKAMI Y., *Stress Intensity Factors Handbook*, Pergamon Press, Oxford, 1987.

[OLU 93] OLUKOKO O. et A. B., «A new boundary element approach for contact problems with friction», *Int. J. for Numerical Methods in Eng.*, vol. 36, 1993, p. 2625-2642.

[POR 92] PORTELA A. et ALIABADI M. H., «The dual boundary element method effective implementation for crack problems», *Int. J. for numerical methods in engineering*, vol. 33, 1992, p. 1269 - 1287.

[THO 87] THOULESS M. D., EVANS A. G. ASHBY M. F. et HUTCHINSON J. W., «The edge cracking and spalling of brittle plates», *Acta Metallurgica*, vol. 35, 1987, p. 1333-1341.