# Conservation de l'énergie et maillages évolutifs en dynamique

Julien Réthoré, Anthony Gravouil, Alain Combescure

► To cite this version:

Julien Réthoré, Anthony Gravouil, Alain Combescure. Conservation de l'énergie et maillages évolutifs en dynamique. Revue Européenne des Éléments Finis, 2004, pp.703-711. 10.3166/reef.13.703-711. hal-00450560

HAL Id: hal-00450560

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

Submitted on 15 Jul 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
# Conservation de l'énergie et maillages évolutifs en dynamique

## Application aux calculs de propagation de fissure avec remaillage

Julien Réthoré — Anthony Gravouil — Alain Combescure

LaMCoS, Laboratoire de Mécanique des Contacts et des Solides
UMR 5514, INSA Lyon

Bat. Jean d'Alembert, 18,20 rue des Sciences
F69621 Villeurbanne

Alain.Combescure@insa-lyon.fr

**RÉSUMÉ.** On présente une étude générale de stabilité des calculs dynamiques à maillages évolutifs par une méthode énergétique. Cette étude, restreinte au cadre des calculs de propagation dynamique de fissures avec remaillage, permet d'établir des conditions pour garantir la stabilité et maîtriser les transferts d'énergie lors des évolutions du maillage. Pour illustrer l'efficacité des méthodes proposées, on analyse les résultats obtenus pour un exemple.

**ABSTRACT.** This paper presents a general study of the stability of variable-mesh dynamic calculations using an energy approach. This study, whose scope is limited to the calculation of dynamic crack propagation with remeshing, enables us to establish the conditions which are necessary to ensure stability and allow control of energy transfers during the evolution of the mesh. The results obtained on an sample problem are analyzed to illustrate the effectiveness of the proposed methods.

**MOTS-CLÉS :** remaillage, dynamique de la rupture, stabilité numérique, bilan énergétique.

**KEYWORDS:** remeshing, dynamic fracture mechanics, numerical stability, energy balance.
# 1. Introduction

Dans de nombreux domaines, on peut être amené à effectuer des calculs au cours desquels, pour diverses raisons, la géométrie de la structure évolue. Il est nécessaire dans certains cas d'utiliser des procédures de remaillage pour construire le maillage de la structure étudiée en prenant en compte sa nouvelle géométrie. Pour une simulation en dynamique, la discrétisation en espace du problème évoluant d'un instant à l'autre et le schéma d'intégration en temps nécessitant de connaître les composantes du vecteur d'état à l'instant $t_n$ pour calculer la solution à l'instant $t_{n+1}$, le vecteur d'état de l'instant $t_n$ doit être projeté sur la discrétisation en espace de l'instant $t_{n+1}$.

Dans les problèmes de simulation dynamique de propagation de fissures en deux dimensions, les méthodes de déboutonnage ne posent pas le problème évoqué plus haut mais présentent l'inconvénient de nécessiter la connaissance a priori du trajet suivi par la fissure. Les autres méthodes font appel à un maillage évolutif. Des instabilités numériques ont été constatées dans [EIN 00] à cause des opérations de projections successives. Le travail des auteurs s'est alors porté sur la qualité de ces projections. Dans [NIS 01], les auteurs décrivent la nécessité des projections de manière formelle.

Dans cet article, le problème est abordé d'un point de vue énergétique. Tout d'abord nous cherchons à écrire une condition de stabilité, puis à établir un bilan d'énergie discrétisée. Ensuite, on propose une méthode, dite de rééquilibrage, permettant de garerantir à la fois la stabilité et la précision du calcul quel que soit le type de projection utilisé. Les résultats sont illustrés au travers d'un exemple pour lequel on compare les résultats obtenus avec et sans rééquilibrage.

# 2. Stabilité et bilan d'énergie discrétisée

## 2.1. Cas général d'un calcul à maillage évolutif

L'intérêt va maintenant être porté au cas d'un calcul avec remaillage. En étendant la méthode exposée dans [HUG 00] et [GRA 00] à ce type de calcul, la stabilité du schéma peut être étudiée. La discrétisation du problème évolue d'un instant à l'autre et le schéma d'intégration en temps nécessite de connaître les composantes du vecteur d'état à l'instant $t_n$ pour calculer la solution à l'instant $t_{n+1}$. Le vecteur d'état de l'instant $t_n$ doit donc être projeté sur la discrétisation en espace de l'instant $t_{n+1}$. L'opérateur de projection $\Pi_{n,n+1}$ est alors noté comme suit :

$$X_n = \Pi_{n,n+1} \bar{X}_n \quad [1]$$

où $X$ représente la discrétisation d'un tenseur du premier ordre, $X_n$ désignant sa valeur calculée à l'instant $t_n$ sur la discrétisation de l'instant $t_n$ et $\bar{X}_n$ sa valeur calculée à l'instant $t_n$ sur la discrétisation de l'instant $t_{n+1}$. Les correspondances entre formes bilinéaires de masse et de raideur pour les deux discrétisations prennent alors, avec des notations similaires, la forme

$$\bar{M}_n = \Pi_{n,n+1}^T M_n \Pi_{n,n+1} \quad [2]$$
$$ \bar{K}_n = \Pi_{n,n+1}^T K_n \Pi_{n,n+1} \quad [3] $$

Le fait de travailler sur deux discrétisations différentes impose aussi de réécrire :

– le schéma de Newmark pour lequel $\gamma \beta$ sont les constantes habituelles :

$$ \left\{ \begin{array}{l} U_{n+1} = \bar{U}_n + \Delta t \bar{U}_n + \Delta t^2 (\frac{1}{2} - \beta) \ddot{U}_n + \Delta t^2 \beta \ddot{U}_{n+1} \\ \dot{U}_{n+1} = \bar{U}_n + \Delta t (1 - \gamma) \ddot{U}_n + \Delta t \gamma \ddot{U}_{n+1} \end{array} \right. \quad [4] $$

– les opérateurs $\langle \rangle$ et $[]$ :

$$ \left\{ \begin{array}{l} \langle X \rangle = \frac{1}{2} (X_{n+1} + \bar{X}_n) \\ [X] = X_{n+1} - \bar{X}_n \end{array} \right. \quad [5] $$

En utilisant la démarche développée dans [HUG 00], et en remarquant que la stabilité du schéma n'est pas influencée par les efforts extérieurs [HUG 78] (le schéma doit être stable quel que soit le problème posé), la condition de stabilité s'écrit :

$$ \langle \ddot{U} \rangle^T A_{n+1}[\ddot{U}] + \langle \dot{U} \rangle^T K_{n+1}[\dot{U}] = -(\gamma - \frac{1}{2})[\ddot{U}]^T A_{n+1}[\ddot{U}] - \frac{1}{\Delta t} [\dot{U}]^T ([M]\ddot{U}_n + [K]\bar{U}_n) \quad [6] $$

avec

$$ A_{n+1} = M_{n+1} + \frac{\Delta t^2}{2} (2\beta - \gamma) K_{n+1} $$

On reconnaît dans le membre de droite le terme en $\gamma - \frac{1}{2}$ dont dépend la stabilité du schéma dans un calcul sans évolution de maillage, les conditions de stabilité étant dans ce cas :

$$ \begin{cases} \frac{1}{2} \le \gamma \le 2\beta & \text{schéma inconditionnellement stable,} \\ \frac{1}{2} \le \gamma & \text{et } 2\beta \le \gamma \text{ schéma stable si } \Delta t \omega_{max} \sqrt{\gamma - 2\beta} \le \sqrt{2} \end{cases} \quad [7] $$

où $\omega_{max}$ est la plus grande fréquence propre de la matrice $A_{n+1}$. Cependant, la pré-sence d'un terme supplémentaire faisant intervenir $[M]$ et $[K]$ remet en cause les conditions de stabilité exprimées dans le cas précédent. En effet, la stabilité du schéma dépend du membre de droite de l'équation (6).

Si on se place dans les cas particuliers des schémas de la différence centrée ($\gamma = \frac{1}{2}$, $\beta = 0$) ou de l'accélération moyenne ($\gamma = \frac{1}{2}$, $\beta = \frac{1}{4}$), le terme $-(\gamma - \frac{1}{2})[\ddot{U}]^T A_{n+1}[\ddot{U}]$ est nul. La stabilité dépend alors de $-\frac{1}{\Delta t} [\dot{U}]^T ([M]\ddot{U}_n + [K]\bar{U}_n)$. Si rien n'est fait pour contrôler ce terme, on ne peut garantir la stabilité du calcul et seule une étude numérique permettra de conclure sur la présence d'instabilité. Ce terme ne fait pas intervenir les paramètres du schéma qui n'est donc pas à remettre en cause en cas d'instabilité.
De la même façon que pour l'étude de stabilité, le bilan énergétique peut aussi être écrit avec les notations définies plus haut. On obtient alors :

$$ [\dot{U}]^T M_{n+1} \langle \dot{U} \rangle + [U^T] K_{n+1} \langle U \rangle = -\frac{\Delta t^2}{2}(2\beta - \gamma)[\ddot{U}]^T M_{n+1} \langle \ddot{U} \rangle \\ -\frac{\Delta t^2}{2}(\gamma - \frac{1}{2})(2\beta - \gamma)[\ddot{U}]^T M_{n+1}[\ddot{U}] \\ -(\gamma - \frac{1}{2})[\ddot{U}]^T K_{n+\frac{1}{2}}[\ddot{U}] \\ -(\gamma - 1)[\ddot{U}]^T([M]\ddot{U}_n + [K]\ddot{U}_n) \quad [8] $$

Ici encore on retrouve les équations d'un problème sans évolution de maillage augmentés d'un terme faisant intervenir $[M]$ et $[K]$.

## 2.2. Application à la mécanique de la rupture : méthode de rééquilibrage

On s'intéresse maintenant à un calcul de propagation dynamique d'une fissure. Dans ce cas particulier, on sait quel phénomène physique provoque le changement de géométrie. En effet, supposons une extension de fissure $da$ entre les instants $t_n$ et $t_{n+1}$. Le vecteur d'état de l'instant $t_n$ ne peut être en équilibre sur cette nouvelle géométrie que si on applique une distribution de force $F_+$ sur l'extension de la fissure pour la refermer. La méthode de rééquilibrage permet via l'annulation du résidu défini par [9], de garantir l'équilibre du vecteur d'état projeté sur la nouvelle discrétisation :

$$ R = M_{n+1}\bar{\ddot{U}}_n + K_{n+1}\bar{U}_n - F_+ \quad [9] $$

On a alors les égalités suivantes :

$$ [U]^T([M]\ddot{U}_n + [K]\bar{U}_n) = [U]^T F_+ \quad [10] $$

$$ [\dot{U}]^T([M]\ddot{U}_n + [K]\bar{U}_n) = [\dot{U}]^T F_+ \quad [11] $$

où $[U]^T F_+$ et $[\dot{U}]^T F_+$ correspondent respectivement à la puissance et au travail de la distribution de force $F_+$. Or, d'un point de vue physique, on sait ([BUI 78]) que ces deux termes sont égaux à, respectivement, $-2Gda$ et $-2G\dot{a}$ ($G$ est le taux de restitution de l'énergie, et $\dot{a}$ la vitesse d'avance de la fissure). Dans l'algorithme utilisé, l'étape de rééquilibrage est effectuée après projection du vecteur d'état de l'instant $t_n$ sur le nouveau maillage. Et elle permet de vérifier les égalités :

$$ [U]^T([M]\ddot{U}_n + [K]\bar{U}_n) = -2Gda \quad [12] $$

$$ [\dot{U}]^T([M]\ddot{U}_n + [K]\bar{U}_n) = -2G\dot{a} \quad [13] $$

Si on revient sur l'étude de stabilité développée plus haut, calculer

$$ I_n = -[\dot{U}]^T([M]\ddot{U}_n + [K]\bar{U}_n) - 2G\dot{a} $$
permet de quantifier l'instabilité effectivement introduite au cours du pas de temps $[t_n, t_{n+1}]$. De la même façon,

$$D_n = -(\gamma - 1)([U]^T([M]\bar{\bar{U}}_n + [K]\bar{U}_n) + 2Gda)$$

mesure la quantité d'énergie injectée ou dissipée dans les opérations effectuées pour passer d'une discrétisation à l'autre.

Le résidu [9] doit être annulé, on peut utiliser pour ce faire une méthode itérative. Cette étape n'est pas habituelle dans les codes de dynamique explicite. On peut donc se poser la question de la stabilité des calculs de dynamique explicite lors de remaillages successifs, voir [EIN 00]. L'utilisation d'une taille de pas de temps petite et de maillages fins minimise les erreurs effectuées lors des projections et permet souvent de maintenir la stabilité des calculs dans ce type de code sans pour autant assurer la conservation de l'énergie. La méthode proposée garantit la stabilité et la conservation de l'énergie lors de remaillages quelles que soient la taille du pas de temps, la finesse de maillage et la qualité des opérateurs de projection.

### 3. Exemple : simulation de la propagation dynamique d'une fissure dans une éprouvette DCB

Pour illustrer l'efficacité de la méthode de rééquilibrage, on effectue un calcul de propagation dynamique d'une fissure dans une éprouvette de type DCB avec remaillage (voir figure 1). La simulation est programmée dans le code de calcul par éléments finis CASTEM 2000 dans lequel on utilise [ATT 97] pour le calcul du taux de restitution de l'énergie. Le schéma utilisé est celui de l'accélération moyenne ($\gamma = \frac{1}{2}$, $\beta = \frac{1}{4}$). Cette éprouvette est constituée d'un matériau considéré élastique homogène et isotrope. Elle est soumise à un chargement à déplacement imposé à ses extrémités. L'évolution de ce chargement est de type échelon. On donne aussi $c = 0,01m$, $b = 0,05m$, $\frac{a}{b} = 0,2$ pour la géométrie et $E = 186GPa$, $\nu = 0,3$, $\rho = 8000kgm^{-3}$, $K_{1c} = 110MPa\sqrt{m}$ pour les propriétés du matériau.

<figure><img src="image_2.png" /><figcaption>**Figure 1.** Définition de la géométrie initiale</figcaption></figure>
La figure 2 présente les évolutions de la longueur de la fissure pour un calcul sans et un calcul avec rééquilibrage. Sans rééquilibrage, le temps à "rupture" (on arrête le calcul lorsque $a/b = 0,7$) est environ deux fois plus court. Pour expliquer ces différences, on peut tracer les bilans énergétiques pour chacun des calculs. Dans les figures 3 et 4, le terme bilan correspond à l'évolution du terme $D_n$ défini plus haut. Si on n'effectue pas d'étape de rééquilibrage, le déséquilibre énergétique est si important que l'échelle utilisée sur le graphique ne permet pas de distinguer les niveaux d'énergie mis en jeu au début du calcul. On compare le déséquilibre cumulé $D = \sum D_n$ sur la figure 5. Sans rééquilibrage, la quantité d'énergie injectée dans le modèle remet en cause la précision du calcul, et amène à s'interroger sur la stabilité de celui-ci. D'un point de vue de la stabilité, de la même façon que pour le bilan énergétique, on peut tracer l'instabilité cumulée $I = \sum I_n$ lors de la propagation de la fissure. Le calcul diverge lentement, on voit sur la figure 6 que l'instabilité est bien présente.

<figure><img src="image_3.png" /><figcaption>**Figure 2.** Comparaison des évolutions de la longueur de la fissure</figcaption></figure>

Cet exemple où le maillage est volontairement grossier et le pas de temps supérieur au pas de temps critique d'un schéma explicite, met en évidence la qualité des résultats obtenus avec la méthode de rééquilibrage et ce, tant sur le plan de la conservation de l'énergie que de la stabilité. Par rapport à un calcul sans rééquilibrage, ces résultats expliquent les différences constatées sur la figure 2 en termes d'évolution de la longueur de la fissure. Cette méthode permet donc d'effectuer un calcul stable au cours duquel la conservation de l'énergie est respectée même quand le maillage est grossier et le pas de temps supérieur au pas de temps critique d'un schéma explicite.
<figure><img src="image_4.png" /><figcaption>**Figure 3.** Bilan énergétique sans rééquilibrage</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>**Figure 4.** Bilan énergétique avec rééquilibrage</figcaption></figure>
<figure><img src="image_6.png" /><figcaption>**Figure 5.** Comparaison du déséquilibre énergétique cumulé lors des projections successives</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>**Figure 6.** Comparaison de l'instabilité cumulée</figcaption></figure>
## 4. Conclusion

Après avoir exposé les difficultés rencontrées lors des simulations dynamiques avec maillages évolutifs, on propose une étude théorique du problème avec une approche énergétique. Cette étude permet de dégager des conditions à respecter afin de garantir la stabilité et la précision de telles simulations. On s’intéresse ensuite au cas particulier de la propagation d’une fissure. Une méthode de rééquilibrage est présentée, ainsi que les résultats fournis pour un exemple. Cette méthode fournit un cadre général pour garantir la stabilité du calcul et contrôler les tranferts d’énergie lors du remaillage. Elle peut être étendue à n’importe quel problème à maillage évolutif et être utilisée dans tous types de code, implicite ou explicite. On peut aussi penser, pour être encore plus général, à son extension à des problèmes où le type et/ou le nombre des degrés de liberté varient au cours du temps. C’est le cas, par exemple, lors de l’utilisation d’éléments finis étendus ou d’éléments d’interface.

## 5. Bibliographie

[ATT 97] ATTIGUI M., PETIT C., « Mixed-mode separation in dynamic fracture mechanics : New path independent integrals », *International Journal of Fracture*, vol. 84, n° 1, 1997, p. 19-36.

[BUI 78] BUI H. D., *Mécanique de la rupture fragile*, Masson, 1978.

[EIN 00] EINSFELD R. A., MARTHA L. F., BITTENCOURT T. N., « Combination of smeared and discrete approaches with the use of interface elements », *European Congress on Computational Methods in Applied Sciences and Engineering*, 2000.

[GRA 00] GRAVOUIL A., Méthode multi-échelles en temps et en espace avec décomposition de domaines pour la dynamique non-linéaire des structures, Thèse, LMT-Cachan, 2000.

[HUG 78] HUGUES T. J. R., LIU W. K., « Implicit-Explicit Finite Element Transient Analysis : Stability Theory », *Journal of Applied Mechanics*, vol. 45, 1978, p. 371-374.

[HUG 00] HUGUES T. J. R., BELYTSCHKO T., *Nonlinear finite element analysis*, ICE Division, Zace Services Ltd., 2000.

[NIS 01] NISHIOKA T., TOKUDOME H., KINOSHITA M., « Dynamic fracture-path prediction in impact fracture phenomena using moving finite element method based on Delaunay automatic mesh generation », *International Journal of Solids and Structures*, vol. 38, 2001, p. 5273-5301.