# 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é, Alain Combescure, Anthony Gravouil

## ► To cite this version:

Julien Réthoré, Alain Combescure, Anthony Gravouil. 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. 17e Congrès Français de Mécanique, Aug 2005, Troyes, France. 10.3166/remn.16.827-843. hal-00507152

HAL Id: hal-00507152

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

Submitted on 18 Aug 2021

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published 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é, Alain Combescure & Anthony Gravouil

Institut National de Sciences Appliquées de Lyon
Laboratoire de Mécanique des Contacts et des Solides
Bât. Jean d'Alembert, 18-20 rue des Sciences 69621 Villeurbanne Cedex
Julien.Rethore@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ément 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 une application au calcul de la propagation dynamique de fissure

## 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. An application to dynamic crack growth simulation is presented

## Mots-clefs :

éléments finis étendus ; élasto-dynamique ; discontinuité en temps ; propagation dynamique de fissure

## 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 capturer les singularités du problème. On propose ici d'étendre les raisonnements appliqués au problème en espace (voir [1] par exemple) à 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'une application à la simulation de la propagation dynamique de fissure.
## 2 Méthode des éléments finis étendus en espace et en temps

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

$$ \mathbf{u} = \mathbf{u}_d \quad \text{sur } \partial\Omega_1 \qquad (1) $$

$$ \sigma(\mathbf{n}) = \mathbf{F}_d \quad \text{sur } \partial\Omega_2 \qquad (2) $$

$$ \sigma(\mathbf{n}) = \mathbf{0} \quad \text{sur } \Gamma^+ \text{ et } \Gamma^- \qquad (3) $$

$$ \mathrm{div}(\sigma) + \mathbf{f}_d = \rho \ddot{\mathbf{u}} \quad \text{dans } \Omega \qquad (4) $$

On note $\Omega$ le domaine matériel considéré et $\Gamma$ une fissure contenu dans ce domaine. $\partial\Omega_1$ est la partie de la frontière $\partial\Omega$ de $\Omega$ sur laquelle on impose un déplacement $\mathbf{u}_d$. $\partial\Omega_2$ est la partie de la frontière sur laquelle on impose un effort $\mathbf{F}_d$ (les lèvres de la fissure $\Gamma^+$ et $\Gamma^-$ sont considérées comme des surface libres). $\mathbf{u}$, $\dot{\mathbf{u}}$ et $\ddot{\mathbf{u}}$ désignent le déplacement la vitesse et l'accélération. $\rho$ est la masse volumique du matériau homogène qui forme le domaine $\Omega$. Son comportement est élastique linéaire (E est la module de Young et $\nu$ le coefficient de Poisson). $\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 $\Omega$ est discrétisé en un ensemble de noeuds $\mathcal{N}$ sur lesquels s'appuient des fonctions de forme éléments finis $\{N_i(\mathbf{x})\}_{i \in \mathcal{N}}$. L'interpolation en espace utilise la méthode de partition de l'unité [2]. 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(\mathbf{x}) = 1 \quad \forall \mathbf{x} \in \Omega \qquad (5) $$

Si on note $U^h(t)$ l'approximation en espace de $\mathbf{u}(\mathbf{x}, t)$, en utilisant les enrichissements proposés dans [1] on obtient :

$$ U^h(t) = \sum_{i \in \mathcal{N}} N_i(\mathbf{x}) U_i(t) + \sum_{i \in \mathcal{N}_{\mathrm{cut}}(t)} N_i(\mathbf{x}) \mathcal{H}(\mathbf{x}) a_i(t) + \sum_{i \in \mathcal{N}_{\mathrm{front}}(t)} \sum_{\alpha} N_i(\mathbf{x}) B_{\alpha}(\mathbf{x}) b_{i,\alpha}(t) \qquad (6) $$

Les fonctions d'enrichissement sont :

$$ \mathcal{H}(\mathbf{x}) = \begin{cases} +1 & \text{si } \mathbf{x} \text{ est au dessus de } \Gamma \\ -1 & \text{si } \mathbf{x} \text{ est en dessous de } \Gamma \end{cases} \qquad (7) $$

et

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

Leurs supports sont : $\mathcal{N}_{\mathrm{cut}}$ l'ensemble des noeuds dont le support a été entièrement coupé par la fissure et $\mathcal{N}_{\mathrm{front}}$ l'ensemble des noeuds dont le support a contenu le front de la fissure. Ainsi, avec $U_i = 0$, $b_{i,\alpha} = 0$ et $a_i = 1$ on peut représenter de façon exacte la discontinuité $\mathcal{H}$ dans le champ de déplacement discrétisé.

### 2.2 Discrétisation en temps

Comme on l'a fait en espace, on discrétise l'intervalle $I$ en un ensemble de $N$ piquets 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, la vitesse présente sous l'effet de la propagation d'un front d'onde des discontinuités en temps. Par conséquent, on va, en utilisant les propriétés de partition de l'unité des fonctions de forme $\{\lambda_j\}_{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}^{M} \sum_{i \in \mathcal{N}_j} \lambda_i(t) H(t - t_j) V_{i,j}^e \quad (9)$$

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$. 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 (10)$$

avec

$$\lambda_n(t) = \frac{t_{n+1} - t}{t_{n+1} - t_n}; \quad \lambda_{n+1}(t) = \frac{t - t_n}{t_{n+1} - t_n} \quad (11)$$

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} \quad (12)$$

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} \quad (13)$$

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}$ peut être augmentée par une condition de continuité faible de la vitesse à l'instant $t_n$ :

$$\int_{t_n^+}^{t_{n+1}^-} \delta V \cdot [M^{n+1} \dot{V}^{h,\Delta t} + K^{n+1} U^{h,\Delta t}] dt + \delta V(t_n^+).M^{n+1} V_{n+1}^e = \int_{t_n^+}^{t_{n+1}^-} \delta V \cdot F(t) dt \quad (14)$$
où $F$ est le vecteur force généralisée et $\delta V$ une fonction test choisie dans l'ensemble des fonctions de forme. 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

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ément finis étendus de problème en temps ont été étudiées dans [3]. L'ordre de convergence est de trois 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 plus haut pour $N_{cut}$ et $N_{front}$. On choisit de conserver l'intégralité des fonctions de forme enrichies prenant en compte d'anciennes positions du front. 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 [4]. 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

<figure><img src="image_2.png" /><figcaption>FIG. 1 – Modélisation d'une fissure semi-infinie dans une plaque infinie</figcaption></figure>

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 [5] une solution analytique pour le facteur d'intensité des contrainte $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 (15)$$

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

Le modèle numérique utilisé est décrit par la Figure 1. 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 = 210GPa$, $\nu = 0.3$ et $\rho = 8000kgm^{-3}$. La contrainte appliquée est $\sigma_0 = 500MPa$. 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 [4]. On s'intéresse au cas où la fissure démarre à $t = 1.5t_c$ puis se propage à la vitesse $\dot{a} = v_0 = 1500m.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 2 (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 [8, 7, 4]). 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.

<figure><img src="image_3.png" /><figcaption>FIG. 2 – Solutions numériques et analytique $\bar{K}_1$ pour une fissure fixe puis mobile</figcaption></figure>

## 4 Conclusion

La partition de l'unité constituée pas 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'application présentée montre 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.

L'approche présentée diffère de celle développée dans [6] où les auteurs s'intéressent à une discontinuité en espace évoluant dans une cadre élément 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\mathbb{A}\subset\text{couplage}$ entre la $r\mathbb{A}\subset\text{solution}$ en espace et la $r\mathbb{A}\subset\text{solution}$ en temps. On adopte aussi une démarche différente de [7] en restant dans le cadre de la mécanique élastique linéaire de la rupture. La stratégie d'enrichissement utilisée (voir [4]) permet de conserver les fonctions d'enrichissement singulières pour des problèmes dépendant 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 [8].
Références

[1] N. Moës, J. Dolbow, and T. Belytschko. A finite element method for crack growth without remeshing. *International Journal for Numerical Methods in Engineering*, 46(1):133–150, 1999.

[2] I. Babuska and JM. Melenk. The partition of unity method. *International Journal for Numerical Methods in Engineering*, 40:727–758, 1997.

[3] J. Réthoré, A. Gravouil, and A. Combescure. A combined space time extended finite element method. *International Journal for Numerical Methods in Engineering*, to appear 2005.

[4] J. Réthoré, A. Gravouil, and A. Combescure. An energy conserving scheme for dynamic crack growth with the extended finite element method. *International Journal for Numerical Methods in Engineering*, in press 2005.

[5] LB. Freund. *Dynamic fracture mechanics*. Cambridge Monographs on Mechanics and Applied Mathematics, 1990.

[6] J. Chessa and T. Belytschko. Arbitrary discontinuities in space-time finite elements. *International Journal for Numerical Methods in Engineering*, 2004 in press.

[7] T. Belytschko, H. Chen, X. Jingxiao, and Z. Goangseup. Dynamic crack propagation based on loss of hyperbolicity and a new discontinuous enrichment. *International Journal for Numerical Methods in Engineering*, 58:1873–1905, 2003.

[8] T. Belytschko and H. Chen. Singular enrichment finite element method for elastodynamic crack propagation. *International Journal of Computational Methods*, 1(1):1–15, 2001.