# Prise en compte de discontinuités en espace et en temps par la méthode des éléments finis étendus

Julien Réthoré, Anthony Gravouil, Alain Combescure

► To cite this version:

Julien Réthoré, Anthony Gravouil, Alain Combescure. Prise en compte de discontinuités en espace et en temps par la méthode des éléments finis étendus. Revue Européenne de Mécanique Numérique/European Journal of Computational Mechanics, 2007, 16 (6-7), pp.827-844. 10.3166/remn.16.827-843. hal-00378154

HAL Id: hal-00378154

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

Submitted on 3 Jun 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
# Prise en compte de discontinuités en espace et en temps par la méthode des éléments finis étendus

## Application à la dynamique de la rupture

Julien Réthoré — Anthony Gravouil — Alain Combescure

LaMCoS, Laboratoire de Mécanique des Contacts et des Structures
UMR 5259, INSA Lyon
Bat. Jean d'Alembert, 18,20 rue des Sciences
F-69621 Villeurbanne
Alain.Combescure@insa-lyon.fr

RÉSUMÉ. On présente une application de la méthode des éléments finis à la dimension temps dans le cadre des problèmes de propagation dynamique de fissure. Ainsi, on exploite les propriétés des éléments finis pour utiliser une base de fonctions de forme enrichie en espace et en temps. Le problème en temps est résolu en utilisant une formulation faible incluant une condition de continuité de la vitesse. Ceci permet de modéliser des problèmes de mécanique présentant des discontinuités à la fois en espace et en temps. On présente des applications au calcul de la propagation dynamique de fissures.

ABSTRACT. This paper is aimed at presenting an application of the extended finite element method to the time variable in the framework of dynamic crack propagation. We take advantage of the partition of the unity properties of finite elements and use an enriched basis of shape functions in space as well as time. To solve the problem in time, we use a weak formulation including a continuity condition for the velocity. It allows to model mechanical problems with space and time discontinuities. Applications to dynamic crack growth simulation are presented.

MOTS-CLÉS : éléments finis étendus ; élasto-dynamique ; discontinuité en temps ; propagation dynamique de fissure.

KEYWORDS: extended finite element ; elasto-dynamics ; time discontinuity ; dynamic crack propagation.
# 1. Introduction

La méthode des éléments finis a été utilisée avec succès dans de nombreux domaines de la simulation des problèmes mécaniques. Cependant il apparaît aujourd'hui un certain nombre de situations difficilement accessibles à cette technique. C'est le cas par exemple des problèmes présentant des discontinuités (propagation de fissure, inclusion, transition de phase...). Ainsi, la méthode des éléments finis étendus utilise les propriétés de partition de l'unité des fonctions de forme éléments finis pour capter les singularités du problème (voir (Babuska et al., 1997; Moës et al., 1999; Remmers et al., 2003) par exemple). Pour des applications quasi statiques, cette technique permet de modéliser les discontinuités d'une manière complètement indépendante du maillage. Si on s'intéresse à la propagation dynamique des fissures, les changements de discrétisation en espace imposés par la propagation introduisent une discontinuité en temps dans le modèle numérique. Cette discontinuité temporelle est une difficulté pour les schémas d'intégration en temps classiques (type Newmark) et la solution est souvent perturbée par des oscillations numériques. Afin d'améliorer le traitement de la discontinuité en temps, on propose ici d'étendre les raisonnements appliqués au problème en espace à la dimension temps. On présente une formulation éléments finis étendus du problème, aussi bien en espace qu'en temps, permettant de modéliser des discontinuités en espace et en temps. Les performances de la méthode sont illustrées au travers d'applications à la simulation de la propagation dynamique de fissures.

# 2. Méthode des éléments finis étendus en espace et en temps

<figure><img src="image_2.png" /><figcaption>**Figure 1. Domaine matériel Ω fissuré**</figcaption></figure>

On se place dans le cadre de l'élasto-dynamique. On résout donc le problème suivant sur l'intervalle $I = [0; T]$ :

$$ 
\begin{align*}
\underline{u} &= \underline{u}_d & \text{sur } \partial\Omega_1 \\
\underline{\sigma}(\underline{n}) &= F_d & \text{sur } \partial\Omega_2 \\
\underline{\sigma}(\underline{n}) &= \underline{0} & \text{sur } \Gamma^+ \text{ et } \Gamma^- \\
\underline{\mathrm{div}}(\underline{\sigma}) + \underline{f}_d &= \rho\underline{\ddot{u}} & \text{dans } \Omega
\end{align*}
\quad [1]
$$
On note Ω le domaine matériel considéré (voir figure 1) et Γ une fissure contenue dans ce domaine. ∂Ω₁ est la partie de la frontière ∂Ω de Ω sur laquelle on impose un déplacement $\underline{u}_d$. ∂Ω₂ est la partie de la frontière sur laquelle on impose un effort $\underline{F}_d$ (les lèvres de la fissure Γ⁺ et Γ⁻ sont considérées comme des surface libres). $\underline{u}$, $\underline{\dot{u}}$ et $\underline{\ddot{u}}$ désignent le déplacement, la vitesse et l'accélération. ρ est la masse volumique du matériau homogène qui forme le domaine Ω. Son comportement est élastique linéaire (E est le module de Young et ν le coefficient de Poisson). $\underline{\sigma}$ est le tenseur symétrique des contraintes de Cauchy dans le cadre de l'hypothèse des petites perturbations.

## 2.1. Discrétisation en espace

Le domaine Ω est discrétisé en un ensemble de nœuds $\mathcal{N}$ sur lesquels s'appuient des fonctions de forme éléments finis {$N_i(\underline{x})$}$_{i \in \mathcal{N}}$. L'interpolation en espace utilise la méthode de partition de l'unité (Babuska et al., 1997). On exploite les propriétés de partition de l'unité des fonctions de forme éléments finis :

$$ \sum_{i \in \mathcal{N}} N_i(\underline{x}) = 1 \quad \forall \underline{x} \in \Omega \quad [2] $$

Si on note $U^h(t)$ l'approximation en espace de $\underline{u}(\underline{x}, t)$, en utilisant les enrichissements proposés dans (Moës et al., 1999) on obtient :

$$ \begin{aligned} U^h(t) = & \sum_{i \in \mathcal{N}} N_i(\underline{x}) U_i(t) + \sum_{i \in \mathcal{N}_{cut}(t)} N_i(\underline{x}) \mathcal{H}(\underline{x}) a_i(t) \\ & + \sum_{i \in \mathcal{N}_{front}(t)} \sum_{\alpha} N_i(\underline{x}) B_{\alpha}(\underline{x}) b_{i,\alpha}(t) \end{aligned} $$

Les fonctions d'enrichissement sont :

$$ \mathcal{H}(\underline{x}) = \begin{cases} +1 & \text{si } \underline{x} \text{ est au dessus de } \Gamma \\ -1 & \text{si } \underline{x} \text{ est en dessous de } \Gamma \end{cases} \quad [3] $$

et

$$ [B_{\alpha}] = \left[ \sqrt{r\sin(\frac{\theta}{2})}, \sqrt{r\cos(\frac{\theta}{2})}, \sqrt{r\sin(\frac{\theta}{2})}\sin(\theta), \sqrt{r\cos(\frac{\theta}{2})}\sin(\theta) \right] \quad [4] $$

Des degrés de libertés enrichis $a_i, b_{i,\alpha}$ sont portés par les nœuds suivants : $\mathcal{N}_{cut}$ l'ensemble des nœuds dont le support a été entièrement coupé par la fissure et $\mathcal{N}_{front}$ l'ensemble des nœuds dont le support a contenu le front de la fissure. Ainsi, avec $U_i = 0, b_{i,\alpha} = 0$ et $a_i = 1$, par exemple, on peut représenter de façon exacte la discontinuité $\mathcal{H}$ dans le champ de déplacement discrétisé.
<figure><img src="image_3.png" /><figcaption>**Figure 2.** Fonction de forme linéaire en temps en fonction enrichie dans $I_{n+1}$</figcaption></figure>

## 2.2. Discrétisation en temps

Comme on l'a fait en espace, on discrétise l'intervalle $I$ en un ensemble de $N$ points de temps $\{t_i\}_{i \in [0..N]}$ régulièrement espacés de $\Delta t$. Sur ces instants discrets s'appuient des fonctions de forme linéaires $\lambda_j(t)$. Dans le cadre de l'élasto-dynamique, le champ des vitesses présente des discontinuités en temps sous l'effet de la propagation d'un front. Par conséquent, on va, en utilisant les propriétés de partition de l'unité des fonctions de forme $\{\lambda_i\}_{i \in [0..N]}$, enrichir l'approximation de la vitesse au moyen de fonctions de Heaviside :

$$V^{h, \Delta t} = \sum_{i=0}^{N} \lambda_i(t) V_i^c + \sum_{j=0}^{N} \sum_{i \in \mathcal{N}_j} \lambda_i(t) H(t - t_j) V_{i,j}^e \quad [5]$$

Les $V^c$ et $V_{i,j}^e$ sont les contributions continue et discontinue des vecteurs vitesse discrétisés sur la base des fonctions de forme en espace. On choisit de limiter le support $\mathcal{N}_j$ des fonctions enrichies $H(t - t_j)$ au seul piquet $t_i$ (voir figure 2). Comme on l'a mentionné dans l'introduction, la prise en compte de la propagation de la fissure introduit une discontinuité dans le modèle numérique. Ces changements dans la discrétisation interviennent aux instants discrets, c'est pourquoi nous choisissons de centrer les fonctions d'enrichissement $H(t - t_j)$ aux piquets de temps. De plus, on a montré dans (Réthoré et al., 2005a) que ce choix conduit à un schéma identique à la méthode de Galerkin discontinue. Néanmoins, ce choix n'est pas une limitation de l'approche présentée.

Comme les $\lambda_j$ ont un support compact $[t_{j-1}; t_{j+1}]$ et que $H(t - t_j)$ est nul pour tout instant antérieur à $t_j$, $H(t - t_n)$ est la seule fonction enrichie active dans $I_{n+1} = ]t_n; t_{n+1}[$. Par conséquent, on peut envisager une approche incrémentale pour résoudre le problème en temps. Dans $I_{n+1}$, l'approximation de la vitesse s'écrit :

$$ \begin{aligned} V^c(t) &= V_n^c \lambda_n(t) + V_{n+1}^c \lambda_{n+1}(t) \\ V^e(t) &= V_{n+1}^e \lambda_n(t) H(t - t_n) \end{aligned} \quad [6] $$
avec

$$ \lambda_n(t) = \frac{t_{n+1} - t}{t_{n+1} - t_n}; \lambda_{n+1}(t) = \frac{t - t_n}{t_{n+1} - t_n} \qquad [7] $$

On note abusivement $V_{n+1}^e = V_{n,n}^e$ car ces degrés de liberté sont obtenus lors de la résolution sur l'intervalle $I_{n+1}$.

Si on décompose les déplacements $U^{h,\Delta t}$ en une partie continue $U^c$ et discontinue $U^e$, en imposant la continuité de $U^{h,\Delta t}$ aux piquets de temps et la relation cinématique $\dot{U}^h = V^h$, on a :

$$ \begin{aligned} U^c(t) &= U_n^c + \int_{t_n}^t V^c(\tau) d\tau \\ U^e(t) &= U_n^e + \int_{t_n^+}^t V^e(\tau) d\tau \end{aligned} \qquad [8] $$

en intégrant analytiquement

$$ \begin{aligned} U_{n+1}^c &= U_n^c + \frac{\Delta t}{2} (V_{n+1}^c + V_n^c) \\ U_{n+1}^e &= U_n^e + \frac{\Delta t}{2} V_{n+1}^e \end{aligned} \qquad [9] $$

ce qui définit la réactualisation du déplacement après résolution de l'équation de conservation de la quantité de mouvement dans $I_{n+1}$.

Pour effectuer cette résolution, on considère $M^{n+1}$ et $K^{n+1}$ les matrices de masse et de rigidité obtenues avec la discrétisation représentant la singularité due à la fissure dans $I_{n+1}$. La formulation de type Galerkin de l'équation de conservation de la quantité de mouvement dans $I_{n+1}$ est augmentée par une condition de continuité faible de la vitesse à l'instant $t_n$ :

$$ \int_{t_n^+}^{t_{n+1}^-} \delta V \cdot \left[ M^{n+1} \dot{V}^{h,\Delta t} + K^{n+1} U^{h,\Delta t} \right] dt + \delta V(t_n^+) \cdot M^{n+1} V_{n+1}^e = \int_{t_n^+}^{t_{n+1}^-} \delta V \cdot F(t) dt \qquad [10] $$

où $F$ est le vecteur force généralisée et $\delta V$ une fonction test choisie dans l'ensemble des fonctions de forme. Le terme $\delta V(t_n^+) \cdot M^{n+1} V_{n+1}^e$ impose de manière faible la continuité en temps de la vitesse à l'instant $t_n$ au sens du produit scalaire défini par l'énergie cinétique. Cette condition permet de controler le saut d'énergie cinétique à $t_n$ et de garantir la stabilité du schéma d'intégration en temps.

Pour les mêmes raisons que précédemment ($\lambda_j$ à support compact, $N_j = t_j$ et $H(t - t_j) = 0 \quad \forall t < t_j$) et parce que $H(t - t_j) = 1 \quad \forall t > t_j$, on choisit les fonctions tests dans $\{\lambda_n, \lambda_{n+1}\}$.
## 2.3. Stabilité, précision, convergence et conservation de l'énergie

<figure><img src="image_4.png" /><figcaption>**Figure 3.** Errur en norme L₂ pour les schémas de Newmark et éléments finis étendus en temps</figcaption></figure>

Les propriétés de stabilité, de précision et de conservation de l'énergie de l'intégrateur en temps obtenu grâce à la formulation éléments finis étendus de problème en temps ont été étudiées dans (Réthoré et al., 2005a). L'ordre de convergence est de trois (voir figure 3) et la stabilité est inconditionnelle.

La stratégie d'enrichissement pour un problème avec une fissure évoluant au cours du calcul est sous-entendue dans les définitions données précédemment pour Ncut et Nfront. On choisit de conserver l'intégralité des fonctions de forme enrichies prenant en compte d'anciennes positions du front (voir figure 4). De plus, on initialise à 0 les nouveaux degrés de liberté correspondant aux fonctions d'enrichissement prenant en compte l'extension de la fissure. Cette stratégie a été étudiée dans (Réthoré et al., 2005b). Elle permet d'assurer la conservation de l'énergie et de maintenir les conditions de stabilité de l'intégrateur en temps lors des évolutions successives de la géométrie de la fissure.

# 3. Applications

## 3.1. Fissure semi-infinie dans une plaque infinie

On s'interesse au problème d'une fissure semi-infinie dans une plaque infinie soumise à une onde de traction. Pour ce problème, on trouve dans (Freund, 1990) une
<figure><img src="image_5.png" /><figcaption>**Figure 4.** Stratégie d'enrichissement en dynamique</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>**Figure 5.** Modélisation d'une fissure semi-infinie dans une plaque infinie</figcaption></figure>
solution analytique pour le facteur d'intensité des contraintes $K_1^{dyn}$ pour une fissure se propageant à une vitesse quelconque :

$$K_1^{dyn}(\dot{a}, t) = \frac{2\sigma_0}{1-\mu} \sqrt{\frac{c_1 t(1-2\mu)}{\pi}} \frac{1-\frac{\dot{a}}{c_r}}{1-\frac{\dot{a}}{2c_r}} \quad [11]$$

où $\dot{a}$ est la vitesse de propagation de la fissure, $c_1$ et $c_r$ les célérités des ondes longitudinales et des ondes de Rayleigh

Le modèle numérique utilisé est décrit par la figure 5. Comme celui-ci est de dimensions finies, on ne pourra comparer les résultats obtenus à la solution analytique que jusqu'au moment où l'onde réfléchie atteint le front de la fissure ($t \le 3t_c = 3H/c_1$). Les dimensions de la plaque sont $H = 2m$, $L = 10m$ et $l = 5m$, les paramètres matériau $E = 210 \text{ GPa}$, $\nu = 0.3$ et $\rho = 8000 \text{ kgm}^{-3}$. La contrainte appliquée est $\sigma_0 = 500 \text{ MPa}$. Le maillage est constitué de 40 × 80 éléments quadrangle linéaires. Pour calculer $K_1^{dyn}$, on utilise une intégrale d'interaction indépendante du domaine calculée autour de la pointe de fissure (Réthoré et al., 2005b). On s'intéresse au cas où la fissure démarre à $t = 1.5t_c$ puis se propage à la vitesse $\dot{a} = v_0 = 1500 \text{ m s}^{-1}$. Les résultats obtenus avec un intégrateur en temps standard (accélération moyenne de Newmark) et la formulation éléments finis étendus en temps sont présentés sur la figure 6 (le facteur d'intensité des contraintes est normalisé par $\sigma_0\sqrt{H}$). Avec un intégrateur en temps classique, la discontinuité en temps provoquée par la propagation de la fissure conduit à des oscillations numériques hautes fréquences qui perturbent la solution (comme dans (Belytschko et al., 2001; Belytschko et al., 2003; Réthoré et al., 2005b)). Avec l'approche présentée ici, la discontinuité est capturée avec précision et la solution obtenue numériquement est très proche de la solution analytique.

## 3.2. Confrontation à l'expérience

Dans la suite de cet article, plusieurs exemples de confrontation expérimentale vont être présentés. Il s'agit de propagation de fissures sous sollicitation mixte donnant lieu à des trajets de fissure complexes. Lors d'une sollicitation en mode mixte, il est nécessaire de déterminer la vitesse et la direction de propagation de la fissure. On considérera que la direction de propagation est gouvernée par l'intensité de contrainte circonférentielle maximum. L'angle de propagation critique $\theta_c$ est alors donné en fonction des facteurs d'intensité des contraintes $K_1^{dyn}$ et $K_2^{dyn}$ :

$$\theta_c = 2\arctan \left[ \frac{1}{4} \left( \frac{K_1^{dyn}}{K_2^{dyn}} - \operatorname{sign}(K_2^{dyn}) \sqrt{8 + \left( \frac{K_1^{dyn}}{K_2^{dyn}} \right)^2} \right) \right] \quad [12]$$

Cette relation, valable en quasi statique et en dynamique pour une fissure fixe, sera aussi utilisée au cours de la propagation. Pour déterminer l'intensité de la sollicitation,
<figure><img src="image_7.png" /><figcaption>**Figure 6.** Solutions numériques et analytique $\bar{K}_1$ pour une fissure fixe puis mobile</figcaption></figure>

on calcule un facteur d'intensité des contraintes équivalent $K_{1eq}^{dyn}$. Dans le cadre de l'utilisation du critère de contrainte circonférentielle maximale, il est égal au facteur d'intensité des contraintes en mode 1 que l'on obtiendrait si la fissure était orientée dans la direction $\theta_c$ :

$$K_{1eq}^{dyn} = \cos^3\left(\frac{\theta_c}{2}\right) K_1^{dyn} - \frac{3}{2}\cos\left(\frac{\theta_c}{2}\right) \sin(\theta_c) K_2^{dyn} \quad [13]$$

Ayant considéré que la propagation était orientée par la direction de contrainte circonférentielle maximale, l'initiation aura alors lieu lorsque le facteur d'intensité des contraintes équivalent atteint une valeur critique $K_{1D}$. On supposera ensuite que la propagation a lieu à intensité de contrainte circonférentielle égale à $K_{1D}$, la vitesse s'adaptant :

$$\begin{cases} \dot{a} = 0 & \text{si } K_{1eq} < K_{1D} \\ K_{1eq}^{dyn} = K_{1D} & \text{si } \dot{a} > 0 \end{cases} \quad [14]$$

On s'attend de plus à ce que $K_{1D}$ dépende de la vitesse de propagation. On trouve dans la littérature (Kanninen et al., 1985) la forme suivante :

$$K_{1D}(\dot{a}) = \begin{cases} K_{1c} & \text{si } \dot{a} = 0 \\ \frac{K_{1A}}{1 - \left(\frac{a}{c_r}\right)^m} & \text{sinon} \end{cases} \quad [15]$$

Par soucis de simplification mais surtout par manque de données expérimentales, on considérera que $m$ vaut 1 et $K_{1A} = K_{1c}$.
<figure><img src="image_8.png" /><figcaption>**Figure 7.** Schématisation des essais Kalthoff</figcaption></figure>

### 3.2.1. Kalthoff (Kalthoff, 2000)

Ce premier cas est dédié à la simulation des essais présentés dans (Kalthoff, 2000). Il s'agit d'essais dynamiques où la sollicitation initiale est en mode 2 pur. La schématisation de ces essais est décrite par la figure 7. Un projectile impacte une plaque dans laquelle deux préfissures ont été pratiquées. La taille du projectile correspond à la distance séparant les deux préfissures. En augmentant la vitesse d'impact $V_0$, Kalthoff a observé un transition d'une rupture de type fragile vers une rupture par propagation de bande de cisaillement. On se placera ici dans un cas $V_0 = 20 \text{ ms}^{-1}$ où la vitesse d'impact provoque une rupture fragile de l'éprouvette par propagation de deux fissures faisant un angle de 65° environ avec les préfissures initiales.

Pour modéliser ces essais, on utilise le modèle décrit par la figure 8. Seule une moitié de l'éprouvette est discrétisée et on impose une condition de symétrie sur l'axe de la symétrie en question. Le projectile est considéré comme infiniment rigide et l'impact est traité en imposant une vitesse horizontale $\dot{u} = V_0$ aux nœuds situés dans la zone impactée. Les dimensions de l'éprouvette sont données par $L = 0.1m$, $l = 0.05m$ et $a = 0.05m$. Le maillage utilisé comporte 80 éléments quadrangle avec une interpolation linéaire comme partition de l'unité. L'intégration en temps est ici effectuée grâce au schéma de Newmark de l'accélération moyenne avec un pas de temps $\Delta t = 1.25\mu s$. On considérera que le matériau constituant l'éprouvette, un acier Maraging 18Ni1900, est homogène de masse volumique $\rho = 8000 \text{ kg m}^{-3}$ et isotrope. On considérera également que son comportement est élastique linéaire. Le module
d'Young $E = 190GPa$, le coefficient de Poisson $\nu = 0.3$ le définissent entièrement. La ténacité de ce matériau est $K_{1c} = 68 MPa\sqrt{m}$.

<figure><img src="image_9.png" /><figcaption>**Figure 8. Modèle numérique pour les essais Kalthoff**</figcaption></figure>

La figure 9 présente les résultats obtenus. On observe la propagation d'une fissure faisant un angle global de 65° environ avec la préfissure, conformément aux observations expérimentales. La propagation s'effectue de manière assez rectiligne avec un vitesse moyenne de $\dot{a} = 0.475c_r$. Néanmoins un léger changement d'orientation est observé à la moitié environ de la propagation. Un examen détaillé de la propagation des ondes dans l'éprouvette permet d'expliquer ce phénomène. En effet, à partir de l'impact un onde de compression se propage dans la partie inférieure de l'éprouvette. Cette onde de compression se réfléchit ensuite sur le bord libre opposé pour se transformer en un onde de traction qui rejoint le front de la fissure dont la propagation a été initiée en mode 2 pur par cisaillement. Ceci donne lieu à une propagation globalement en mode 1 après l'initiation. L'onde de traction se réfléchit ensuite sur le bord libre gauche de l'éprouvette et rejoint le front de la fissure alors qu'il se situe environ à mi-hauteur. La propagation est alors légèrement déviée. Cette explication reste malgré tout assez simpliste dans la mesure où elle ne prend pas en compte la superposition et les interférences des ondes élastiques. C'est pourquoi la fissure n'est que légèrement déviée par l'arrivée du front d'onde.

Les résultats obtenus pour cet exemple sont en accord à la fois avec les observations expérimentales mais aussi avec d'autres résultats numériques de la littérature (Belytschko *et al.*, 2003; de Borst *et al.*, 2006). De plus, la méthode des éléments finis étendus qui permet d'effectuer des calculs de propagation dynamique de fissures sans projection de champs (Réthoré *et al.*, 2005b; Réthoré, 2005), autorise une interprétation avancée par confrontation calcul/expérience.
<figure><img src="image_10.png" /><figcaption>Figure 9. Maillage déformé à t = 32.µs, 44.5µs, 57.0µs, 69.5µs, 82.0µs et 94.5µs</figcaption></figure>
### 3.2.2. Maigre, Rittel (Maigre et al., 1995)

Dans ce deuxième exemple, on propose de simuler les essais décrits dans (Maigre et al., 1995). Il s'agit d'essais effectués sur une géométrie d'éprouvette de type Compact Compression. Ces éprouvettes sont sollicitées grâce à un dispositif expérimental de type bar d'Hoptkinson. Les auteurs ont testé différentes configurations dont celle que nous avons choisie. La modélisation numérique est présentée par la figure 10. Dans la configuration choisie, seul le pied gauche de l'éprouvette est relié à la barre entrante. On impose une force correspondant aux relevés expérimentaux dans cette zone.

<figure><img src="image_11.png" /><figcaption>**Figure 10.** Modèle numérique pour les essais sur éprouvettes CC</figcaption></figure>

La géométrie de ces éprouvettes est définie par la hauteur $0.07 \ m$, la largeur $0.06 \ m$, la hauteur de l'arche $0.035 \ m$ et la distance entre les pieds $0.02 \ m$. Une pré-fissure est effectuée sur l'axe de symétrie de l'éprouvette à partir de la face supérieure et sur une longueur de $0.0183 \ m$. Le maillage de l'éprouvette est également décrit par la figure 10. Le calcul s'effectue avec une taille de pas de temps de $1 \ \mu s$. Les éprouvettes sont en PMMA, matériau qui sera considéré homogène de masse volumique $\rho = 1180 \ kg \ m^{-3}$ et isotrope. Le module d'Young $E = 5.76 \ GPa$, le coefficient de Poisson $\nu = 0.42$ définissent son comportement supposé élastique linéaire. Sa ténacité est $K_{1c} = 4 \ MPa\sqrt{m}$.
<figure><img src="image_12.png" /><figcaption>**Figure 11.** Maillage déformé à $t = 110\mu s, 120\mu s, 130\mu s$ et $140\mu s$</figcaption></figure>

Les résultats sont donnés par la figure 11. Le comportement global de l'éprouvette est de type « pince ». La mise en flexion du pied relié à la barre entrante provoque l'ouverture de la fissure et l'initiation de la propagation dans la direction de ce pied. Ensuite la fissure décrit un trajet curviligne pour rejoindre l'axe de l'éprouvette en fin de propagation. La figure 11 montre la distribution de la contrainte de cisaillement $\sigma_{xy}$ à l'instant $t = 100\mu s$. On observe à cet instant une zone de cisaillement intense sur la gauche de la fissure. L'arrivée du front de fissure dans zone va provoquer la décharge de cette contrainte de cisaillement et le changement d'orientation de la fissure (voir figure 11). Dans cette exemple également, les résultats numériques permettent une interprétation avancée des observations expérimentales et une meilleure compréhension des mécanismes gouvernant la propagation d'une fissure à l'échelle de la structure.
<figure><img src="image_13.png" /><figcaption>**Figure 12.** Carte de $\sigma_{xy}$ à l'instant $t = 100\mu s$</figcaption></figure>

**4. Conclusion**

La partition de l'unité constituée par les fonctions de forme éléments finis est utilisée pour enrichir l'approximation en espace puis en temps. Dans le cadre de la mécanique de la rupture, on utilise des enrichissements singuliers et discontinus dans l'approximation en espace afin de prendre en compte la présence d'une fissure. Comme les problèmes de dynamique de la rupture sont aussi fortement discontinus en temps, l'interpolation en temps est enrichie par des fonctions discontinues. Une formulation faible du problème semi-discret permet d'ajouter une condition de continuité faible en temps au champ de vitesse.

L'approche présentée diffère de celle développée dans (Chessa *et al.*, 2004) où les auteurs s'intéressent à une discontinuité en espace évoluant dans un cadre éléments finis espace-temps. Ici, la discontinuité modélisée concerne à la fois l'espace et le temps mais la formulation permet de conserver un découplage entre la résolution en espace et la résolution en temps. On adopte aussi une démarche différente de (Belytschko *et al.*, 2003) en restant dans le cadre de la mécanique élastique linéaire de la rupture. La stratégie d'enrichissement utilisée (voir (Réthoré *et al.*, 2005b)) permet de conserver les fonctions d'enrichissement singulières pour des problèmes dépendants
du temps. On peut ainsi modéliser la propagation dynamique d'une fissure de géométrie quelconque sans pour autant avoir à supposer le trajet de la fissure comme dans (Belytschko et al., 2001).

Les applications présentées montrent la robustesse des résultats obtenus grâce à la méthode des éléments finis étendus appliquée à la fois au problème en espace et au problème en temps. La confrontation des résultats numériques avec les observations expérimentales dans des cas complexes (mode mixte, trajet de fissure curviligne) apporte des éléments de compréhension de la propagation dynamique de fissure dans les structures. Bien que le critère de propagation utilisé soit relativement simple, il permet de bien rendre compte de l'influence de la propagation des ondes sur le trajet de la fissure, phénomène prépondérant à l'échelle de la structure. Dans des travaux plus récents (Grégoire et al., 2007), les performances de la méthode présentée ici sont illustrées sur un cas où la propagation de la fissure se fait selon un trajet curviligne mais aussi avec un arrêt et un redémarrage. Ce cas particulièrement discriminant renforce la crédibilité des résultats obtenus par la méthode proposée dont le principal avantage concerne la conservation de l'énergie discrète.

## 5. Bibliographie

Babuska I., Melenk J., "The Partition of unity method", *International Journal for Numerical Methods in Engineering*, vol. 40, 1997, p. 727-758.

Belytschko T., Chen H., "Singular enrichment finite element method for elastodynamic crack propagation", *International Journal of Computational Methods*, vol. 1, n° 1, 2001, p. 1-15.

Belytschko T., Chen H., Jingxiao X., Goangseup Z., "Dynamic crack propagation based on loss of hyperbolicity and a new discontinuous enrichment", *International Journal for Numerical Methods in Engineering*, vol. 58, 2003, p. 1873-1905.

de Borst R., Remmers J., Needleman A., "Mesh-independent numerical representations of cohesive-zone models", *Engineering Fracture Mechanics*, vol. 173, n° 2, 2006, p. 160-177.

Chessa J., Belytschko T., "Arbitrary discontinuities in space-time finite elements", *International Journal for Numerical Methods in Engineering*, vol. 61, n° 15, 2004, p. 2595-2614.

Freund L., *Dynamic fracture mechanics*, Cambridge Monographs on Mechanics and Applied Mathematics, 1990.

Grégoire D., Maigre H., Réthoré J., Combescure A., "Dynamic crack propagation under mixed-mode loading - Comparison between experiments and X-FEM simulations", *International Journal of Solids and Structures*, vol. 44, 2007, p. 6517-6534.

Kalthoff J., "Modes of dynamic shear failure in solids", *International Journal of Fracture*, vol. 101, 2000, p. 1-31.

Kanninen M., Popelar C., *Advanced fracture mechanics*, Oxford University Press, New York, 1985.

Maigre H., Rittel D., "Dynamic fracture detection using the force displacement reciprocity : application to the compact compression specimen", *International Journal of Fracture*, vol. 73, 1995, p. 67-79.
Moës N., Dolbow J., Belytschko T., "A finite element method for crack growth without reme-shing", *International Journal for Numerical Methods in Engineering*, vol. 46, n° 1, 1999, p. 133-150.

Remmers J., de Borst R., Needleman A., "A cohesive segments method for the simulation of crack growth", *Computational Mechanics*, vol. 31, 2003, p. 69-77.

Réthoré J., Méthode élément finis étendus en espace et en temps : application à la propagation dynamique des fissures, Thèse, Institut National de Sciences Appliquées de Lyon, 2005.

Réthoré J., Gravouil A., Combescure A., "A combined space time eXtended Finite Element Method", *International Journal for Numerical Methods in Engineering*, vol. 64, 2005a, p. 260-284.

Réthoré J., Gravouil A., Combescure A., "An Energy Conserving Scheme for Dynamic Crack Growth with the eXtended Finite Element Method", *International Journal for Numerical Methods in Engineering*, vol. 63, 2005b, p. 631-659.