# Propagation de fissures trois dimensions et en dynamique rapide

Romain Pelee, Vincent Faucher, Thomas Elguedj, Alain Combescure, Benoît Prabel

► To cite this version:

Romain Pelee, Vincent Faucher, Thomas Elguedj, Alain Combescure, Benoît Prabel. Propagation de fissures trois dimensions et en dynamique rapide. 11e colloque national en calcul des structures, CSMA, May 2013, Giens, France. hal-01717068

HAL Id: hal-01717068

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

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
# Propagation de fissures 3 dimensions et en dynamique rapide

R. Pelee<fnref n="1" />,<fnref n="2" />, V. Faucher<fnref n="2" />,<fnref n="3" />, T. Elguedj<fnref n="1" />, A. Combescure<fnref n="1" />, B. Prabel<fnref n="2" />

**Résumé** — Le but de ces travaux a été de développer dans le code de dynamique explicite du CEA, EUROPLEXUS, une méthode de propagation des fissures en 3 dimensions dans des matériaux ductiles. Pour cela, la méthode X-FEM a été utilisée pour éviter les remaillages, couplée aux level-sets pour la représentation implicite de la fissure.

**Mots clés** — X-FEM, dynamique explicite, fissure.

# 1 La méthode X-FEM

## 1.1 Cinématique du déplacement

La méthode X-FEM, ou méthode des éléments finis étendus, a été proposée par Belytschko et al. [1]. Elle permet, dans le cas d'une structure avec une fissure, de pouvoir représenter et faire évoluer la discontinuité indépendamment du maillage. Pour cela, cette méthode se base sur la partition d'unité décrit par Babuska et al. [2], et permet d'approximer le champ de déplacement $\underline{u}$ par la somme de l'interpolation des déplacements standards $U_i$ et de l'interpolation des déplacements enrichis $U_j^e$ :

$$ \underline{u}(\underline{x}) = \sum_{i=1}^{n} N_i(\underline{x})U_i + \sum_{j=1}^{n_d} H(\underline{x})N_j(\underline{x})U_j^e \quad (1) $$

## 1.2 Intégration

Dans notre cas, nous considérons uniquement les enrichissements sauts à l'aide de la fonction de Heaviside $H(\underline{x})$. Du fait de cette discontinuité, le calcul de l'intégrale élémentaire est rendu inexacte par la règle de quadrature de Gauss, donnée par l'équation 2. La stratégie jusque là utilisée consistait à sous-découper l'élément de référence en plusieurs sous-éléments, et à placer un certain nombre de points de Gauss dans chaque sous-élément. Dans le cas d'un découpage régulier, le sous-découpage et le nombre de points de Gauss sont illustrés figure 1(a) dans un cas 2D. Cette méthode est très couteuse en temps de calcul car elle demande un nombre important de points de Gauss (64 en 2D et 512 en 3D).

$$ \int_{-1}^{+1} H(\underline{\xi})N(\underline{\xi})d\underline{\Omega} = \sum_{k=1}^{n_g} w_k H(\underline{\xi}_k)N(\underline{\xi}_k) \quad (2) $$

De plus, pour éviter de projeter les champs de contraintes et de déformations pendant le calcul, les points de Gauss doivent rester fixes dans l'élément de référence. On a donc décidé de fixer les points de Gauss, de calculer le plus exactement possible les intégrales $I$ et de calculer ensuite les poids de Gauss $w_k$ correspondant afin que l'intégrale soit le plus juste possible avec un faible nombre de points de Gauss. Lorsque l'élément est totalement coupé par la fissure, on calcule analytiquement la valeur de l'intégrale $I$. Si l'élément contient le front de fissure, du fait que l'on a que des enrichissements avec les fonctions de Heaviside, l'erreur est minimisée. La figure 1(b) illustre également l'emplacement des nouveaux points de Gauss et les paramètres choisis pour le calcul de l'intégrale analytique. Le système à résoudre pour trouver les poids $w_k$ est donné en équation 3.

$$ \underline{w}_k = N(\underline{\xi}_k)^{-1}I \quad (3) $$

<footnote n="1">¹ INSA-Lyon, LaMCoS UMR5259 - FRANCE, {romains.pelee-de-saint-maurice, thomas.elguedj, alain.combescure} @insa-lyon.fr</footnote>

<footnote n="2">² CEA, DEN, DANS, DM2S, SEMT, DYN - FRANCE, {romains.peleedesaintmaurice, vincent.faucher, benoit.prabel} @cea.fr</footnote>

<footnote n="3">³ LaMSID UMR CNRS-EDF-CEA 8193, F-92141 Clamart - FRANCE, vincent.faucher@cea.fr</footnote>
<figure><img src="image_2.png" /><figcaption>Fig. 1 – Emplacement des points de Gauss pour un sous-découpage auparavant utilisé 1(a) et après simplification 1(b).</figcaption></figure>

En 2D, pour un quadrangle à 4 noeuds, cela revient à résoudre un système de 8 équations (4 fonctions de formes classiques et 4 fonctions discontinues) et en 3D, pour un tétraédre à 8 noeuds, il faut résoudre un système de 16 équations (8 fonctions de formes classiques et 8 fonctions discontinues). Cette méthode permet d'avoir le calcul des forces internes exacte pour les éléments totalement coupés par une fissure droite.

# 2 Discrétisation du front de fissure et propagation des level-sets

## 2.1 Discrétisation du front de fissure

Afin de représenter la fissure, on utilise un second maillage externe régulier (grille auxiliaire), plus fin que le maillage éléments finis, où sont calculées 2 level-sets. La première level-set $\phi_t$ sert à décrire le "plan" de la fissure et la seconde $\phi_n$ pour repérer le front. Une première étape consiste à discretiser le front de fissure afin d'obtenir un certain nombre de points où calculer les grandeurs caractéristiques de la mécanique de la rupture. Pour cela, on selectionne les noeuds du maillage level-sets qui sont le plus proches du front, puis on projète ces points sur le front par un algorithme itératif. L'algorithme 1 et la figure 2 illustrent les différentes étapes de la projection.

<figure><img src="image_3.png" /><figcaption>Fig. 2 – Projection des noeuds sur le front de fissure en 2 dimensions</figcaption></figure>
**Algorithm 1** algorithme de discretisation du front de fissure

**Require:** Noeuds les plus proches du front $P_{nd}$

**Require:** Point du front $P_{fr}$

**Require:** Nombre de noeuds les plus proches du front $n_{nd}$

**for** i = 1 : tous les noeuds du front **do**

Calcul de la valeur des level-sets $\phi_t$ et $\phi_n$ au point $P_{nd}^i$

Calcul du critère : $\text{criterion} = \sqrt{\phi^2 + \psi^2}$

$P_{prec} \leftarrow P_{nd}^i$

**while** ($\text{criterion} \ge \varepsilon$) **do**

Calcul de la valeur des level-sets $\phi_t$ et $\phi_n$ au point $P_{prec}$

Calcul de la valeur des normales aux level-sets $n_{\phi_t}$ et $n_{\phi_n}$
au point $P_{prec}$.

$P_{sui} \leftarrow P_{prec} + \phi_t n_{\phi_n} + \phi_n n_{\phi_t}$

Calcul de la valeur des level-sets $\phi_t$ et $\phi_n$ au point $P_{sui}$

Calcul du critère : $\text{criterion} = \sqrt{\phi^2 + \psi^2}$

$P_{prec} \leftarrow P_{sui}$

**end while**

$P_{fr}^i \leftarrow P_{sui}^i$

**end for**

## 2.2 Propagation des vitesses du front vers tout le maillage level-sets

Une fois le front discrétisé, les vitesses sont calculées par la mécanique de la rupture décrite au paragraphe 3. En chaque point de discrétisation du front, il faut étendre ces vitesses à tout le maillage level-sets. la méthode utilisée jusqu'ici dans la littérature consistait à résoudre les équations d'Hamilton-Jacobi de propagation d'un champs le long d'une interface, décrites par Gravouil et al. [3]. Cette méthode est couteuse en temps de calcul et souffre d'un soucis de robustesse. Pour cette raison, nous avons choisis de propager les vitesses de façon géométrique, comme illustrée figure 3. Les 2 points du front les plus proches du noeud du maillage où l'on souhaite calculer la vitesse sont recherchés, et si la projection du noeud est localisée entre ces 2 points, sa vitesse est calculée linéairement par l'équation 4. Sinon, la vitesse du point le plus proche est donnée directement au noeud. $d_a$ et $d_b$ sont les distances du noeud M au point respectivement a et b du front.

<figure><img src="image_4.png" /><figcaption>Fig. 3 – Propagation des vitesses des points du front vers tout le maillage lvl-set</figcaption></figure>

$$ \dot{a}_M = \frac{d_a}{d_a + d_b} \dot{a}_b + \frac{d_b}{d_a + d_b} \dot{a}_a \quad (4) $$

## 2.3 Propagation des level-sets

Il existe plusieurs méthodes pour faire propager les level-sets, la plus utilisée jusqu'ici est la méthode de résolution des équations d'Hamilton Jacobi, décrites par Gravouil et al. [3]. Cependant cette méthode,
comme dans le chapitre 2.2, souffre d'un problème de robustesse, qui est d'autant plus important avec des
grands angles de propagation. Pour cette raison, nous avons choisi d'utiliser une méthode géométrique
d'actualisation des level-sets. Duflot [4] a comparé différentes méthodes de propagation de level-sets en
se basant sur le calcul de l'intégrale d'interaction. Nous avons étendu en 3D sa méthode $\phi'r\theta\psi$ mais elle
nécessite de connaître explicitement en chaque point du maillage, la distance au front de fissure (notée
r). De plus, notre critère de propagation de fissure présenté dans le chapitre 3 ne nécessite pas d'avoir
une intégrale d'intéraction exacte. Pour cette raison, nous avons choisi d'utiliser une méthode proche
de la méthode $\phi\psi r\theta$ de Duflot. Le domaine complet du maillage est sous découpé en 4 zones, chaque
zone possède son équation d'actualisation. L'ensemble des zones et des équations sont présentées dans
la figure 4 et le tableau 1.

$$
\theta = \arctan \left( \frac{V_{\phi_t}}{V_{\phi_n}} \right) \qquad (5)
$$

$$
\dot{a} = \sqrt{V_{\phi_t}^2 + V_{\phi_n}^2} \qquad (6)
$$

<table><caption>Tableau 1 – La Méthode d’actualisation des level-sets</caption>
  <thead>
    <tr>
      <th>Zones</th>
      <th>Zone 1</th>
      <th>Zone 2</th>
      <th>Zone 3</th>
      <th>Zone 4</th>
    </tr>
  </thead>
  <tbody>
    <tr>
      <th>Conditions</th>
      <td>-tan(θ)φ<sub>t</sub><sup>0</sup> ≥ φ<sub>n</sub><sup>0</sup><br>φ<sub>n</sub><sup>0</sup> ≤ 0</td>
      <td>-tan(θ)φ<sub>t</sub><sup>0</sup> &lt; φ<sub>n</sub><sup>0</sup><br>φ<sub>n</sub><sup>0</sup> ≤ 0</td>
      <td>-tan(θ)φ<sub>t</sub><sup>0</sup> ≤ φ<sub>n</sub><sup>0</sup><br>φ<sub>n</sub><sup>0</sup> &gt; 0</td>
      <td>-tan(θ)φ<sub>t</sub><sup>0</sup> &gt; φ<sub>n</sub><sup>0</sup><br>φ<sub>n</sub><sup>0</sup> &gt; 0</td>
    </tr>
    <tr>
      <th>φ<sub>t</sub></th>
      <td>φ<sub>t</sub><sup>0</sup></td>
      <td>Equation 7</td>
      <td>Equation 8</td>
      <td>Equation 9</td>
    </tr>
    <tr>
      <th>φ<sub>n</sub></th>
      <td>Equation 10</td>
      <td colspan="3">Equation 11</td>
    </tr>
  </tbody>
</table>

$$
\phi_t = \min(\|\langle \phi_t^0 \rangle\|, \|(\phi_t^0 \cos(\theta) - \phi_n^0 \sin(\theta))\|) \cdot \text{sign}(\phi_t^0) \quad (7)
$$

$$
\phi_t = \phi_t^0 \cos(\theta) - \phi_n^0 \sin(\theta) \qquad (8)
$$

$$
\phi_t = \sqrt{(\phi_t^0)^2 + (\phi_n^0)^2} \operatorname{sign}(\phi_t^0) \quad (9)
$$

$$
\phi_n = \phi_n^0 + \phi_t^0 \tan(\theta) - \frac{\dot{a}\Delta t}{\cos(\theta)} \quad (10)
$$

$$
\phi_n = (\phi_n^0 \cos(\theta) + \phi_t^0 \sin(\theta)) - \dot{a}\Delta t \quad (11)
$$

<figure><img src="image_5.png" /><figcaption>Fig. 4 – Illustration des différentes zones</figcaption></figure>
Cette méthode a donné des très bons résultats, comme par exemple, pour un angle de propagation variant de +50° à -50° dans l'épaisseur d'un bloc cubique de 40 × 40 × 10 élément. La vitesse est calculée pour que la fissure parcourt la moitié du bloc en 40 pas de temps. La figure 5 représente les level-sets après 5 pas de temps d'avancée.

<figure><img src="image_6.png" /><figcaption>Fig. 5 – Propagation des level-sets</figcaption></figure>

## 3 Critère de propagation de fissures en 3 dimensions.

Le critère de propagation choisi est un critère local en pointe de fissure utilisé entre autre par Houbussa [5]. Il consiste à moyenner les grandeurs mécaniques (contraintes, déformations...) en pointe de fissure en appliquant une fonction poids donnée en équation 12.

$$w(r) = e^{-\alpha \frac{r^2}{R^2}} \qquad (12)$$

Cette fonction poids de type Gaussienne permettra de donner plus ou moins d'importance aux points d'intégrations situés autour de la pointe de fissure (cf. figure 6(a) et 6(b)), en fonction de la longueur locale R. En pratique, on utilise une valeur R sur une longueur physique (n fois la taille de grains) et pour laquelle on a identifié la contrainte critique de propagation de la fissure $\sigma_{Ic}$. Le maillage doit être alors suffisamment fin (de l'ordre de 3 fois plus fin que R). Cette équation, appliquée au calcul des contraintes pour un domaine $\Omega$ centré sur le front de la fissure, est illustrée équation 13.

$$\tilde{\sigma} = \frac{\int_{\Omega} \underline{\underline{\sigma}} w(r) d\Omega}{\int_{\Omega} w(r) d\Omega} \qquad (13)$$

<figure><img src="image_7.png" /><figcaption>Fig. 6 – Domaine $\Omega$ pour le calcul des contraintes en pointe de fissure en 2D et en 3D.</figcaption></figure>
Un critère est ensuite calculé et comparé à la contrainte critique $\sigma_{Ic}$ afin de savoir si la fissure va propager ou pas. Si elle propage, la vitesse d'avancée de la fissure $\dot{a}$ est donnée par la formule de Kanninen [7], explicité dans l'équation 14, où $c_r$ est la vitesse des ondes de Rayleigh dans le matériau.

$$ \dot{a} = \left(1 - \frac{\sigma^{eq}}{\sigma_{Ic}}\right) c_r \qquad (14) $$

$$ \sigma^{eq} = \max(\tilde{\sigma}_I, \tilde{\sigma}_{II}, \tilde{\sigma}_{III}) \qquad (15) $$

Dans le cas 3D, il faut déterminer 2 angles pour avoir la direction de propagation de la fissure, $\theta$ donne la direction d'avancée du front dans le plan et $\phi$ donne la direction hors plan, ces 2 angles sont représentés figure 6(a). les équations qui permettent de les calculer sont données dans les équations 19 à 21. Elles sont obtenues en recherchant les valeurs d'angles qui minimisent la valeur de la contrainte circonférentielle $\sigma_{\theta\theta}$. Pour cela, Haboussa [6] a utilisé une notation avec des coordonnées barycentriques de $K_I$, $K_{II}$ et $K_{III}$ décrite dans les équations 16 à 18.

$$ K_I^n = \frac{K_I}{K_I + \|K_{II}\| + \|K_{III}\|} \qquad (16) $$

$$ K_{II}^n = \frac{\|K_{II}\|}{K_I + \|K_{II}\| + \|K_{III}\|} \qquad (17) $$

$$ K_{III}^n = \frac{\|K_{III}\|}{K_I + \|K_{II}\| + \|K_{III}\|} \qquad (18) $$

$$ \theta_c = 2\operatorname{sign}(K_{II}) \arctan \left( \frac{1}{4} \left( \hat{K} - \sqrt{(\hat{K})^2 + 8} \right) \right) \qquad (19) $$

$$ \hat{K} = \frac{1 + K_I^n - (1 - K_{III}^n)^{p(v)}}{K_{II}^n} \qquad (20) $$

$$ \phi_c = \frac{1}{2} \operatorname{sign}(K_{III}) \arctan \left( \frac{2\sigma_{\theta z}(\theta_c)}{\sigma_{\theta \theta}(\theta_c) - \sigma_{zz}(\theta_c)} \right) \qquad (21) $$

(22)

avec $p(v) = \frac{1}{4}(\sqrt{\pi} - 5v)$.

Dans notre modèle, on ne calcule pas les valeurs des facteurs d'intensité des contraintes, les valeurs du tenseur des contraintes moyennées en pointe de fissure pilotera donc les directions et les avancées de fissure, en substituant $\tilde{\sigma}_{22}$, $\tilde{\sigma}_{12}$ et $\tilde{\sigma}_{23}$ en lieu et place de $K_I$, $K_{II}$ et $K_{III}$.

# 4 Résultat de propagation de fissure en 3 dimensions.

Un premier cas test a été calculé, un barreau entaillé par une fissure de 22° centré dans la hauteur, la figure 7 représente l'éprouvette et le chargement. Ce barreau est soumis à une pression uniforme sur ses surfaces supérieures et inférieures. Les dimensions du barreau et les caractéristiques matériau sont données dans le tableau 2 et 3. Le résultat de la propagation de la fissure pour un maillage de 50 × 20 × 10 éléments est donné dans les figures 8(a) à 8(b), on remarque une première étape de propagation de la fissure dans le plan tangent de la fissure, puis les bords inclinés de la fissure vont se rejoindre au centre de l'éprouvette, pour former une forme finale hélicoïdale.
<figure><img src="image_8.png" /><figcaption>Fig. 7 – Géométrie de l'éprouvette.</figcaption></figure>

<table><caption>Tableau 2 – Dimension de l'éprouvette et valeur du chargement</caption><tr><td>Longueur</td><td>A</td><td>10.0m</td></tr><tr><td>Hauteur</td><td>B</td><td>4.0m</td></tr><tr><td>Épaisseur</td><td>ep</td><td>2.0m</td></tr><tr><td>Pression</td><td>σ₀</td><td>500.0Mpa</td></tr></table>

<table><caption>Tableau 3 – Caratéristique matériau de l'éprouvette</caption><tr><td>Module de Young</td><td>E</td><td>210.0 Gpa</td></tr><tr><td>Coefficient de Poisson</td><td>ν</td><td>0.3</td></tr><tr><td>Masse volumique</td><td>ρ</td><td>8000.0 Kg.m<sup>-3</sup></td></tr><tr><td>Contrainte de démarrage de la fissure</td><td>σ<sub>Ic</sub></td><td>50 Mpa</td></tr><tr><td>Célérité des ondes de Rayleigh</td><td>c<sub>r</sub></td><td>1200.0 m.s<sup>-1</sup></td></tr><tr><td>Rayon local</td><td>R</td><td>500.0 Mpa</td></tr></table>

<figure><img src="image_9.png" /><figcaption>Fig. 8 – Trajet de la fissure au temps $t = 0\mu s$ et $25\mu s$.</figcaption></figure>

Références

[1] T. BELYTCSCHKO and T. BLACK. Elastic crack growth in finite elements with minimal remeshing. *International journal for numerical methods in engineering*, 45:601–620, 1999.
[2] I. BABUSKA and J. M. MELENK. The partition of unity method. *International Journal for Numerical Methods in Engineering*, 40(4):727-758, 1997.

[3] A. GRAVOUIL, N. MOES, and T. BELYTCSCHKO. Non-planar 3d crack growth by extended finite element and level sets-part 2 : Level set update. *International journal for numerical methods in engineering*, 53 :2569–2586, 2002.

[4] M. DUFLOT. A study of the representation of cracks with level sets. *International journal for numerical methods in engineering*, 70 :1261–1302, 2006.

[5] D. HABOUSSA. *Modélisation de la transition traction-cisaillement des métaux sous choc par la X-FEM.*
PhD thesis, INSA de Lyon, november 2012.

[6] D. HABOUSSA, D. GRÉGOIRE, T. ELGUEDJ, H. MAIGRE, and A. COMBESCURE. X-FEM analysis of the effects of holes or other cracks on dynamic crack propagations. *Int. J. Numer. Meth. Engng*, 86 :618–636, 2011.

[7] M.F. KANNINEN and C.H. POPELAR. *Advanced Fracture Mechanics*. Oxford Engineering Science Series.
Oxford University Press, USA, 1985.