# Étude de la stabilité d'une méthode à raffinement espace-temps automatique (STAR-method)

Pauline Cavin, Anthony Gravouil, Ton Lubrecht, Alain Combescure

► To cite this version:

Pauline Cavin, Anthony Gravouil, Ton Lubrecht, Alain Combescure. Étude de la stabilité d'une méthode à raffinement espace-temps automatique (STAR-method). 7e colloque national en calcul des structures, CSMA, May 2005, Giens, France. hal-01812897

HAL Id: hal-01812897

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

Submitted on 12 Jun 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
Etude de la stabilité d'une méthode à raf-
finement espace-temps automatique (STAR-
method)

Pauline Cavin — Anthony Gravouil — Ton Lubrecht — Alain Com-
bescure

Laboratoire de Mécanique des Contacts et des Solides
UMR 5514
INSA Lyon, Bât. Jean d'Alembert
20 Avenue Albert Einstein
69621 Villeurbanne Cedex
pauline.cavin@insa-lyon.fr

RÉSUMÉ. Ce travail propose un solveur éléments finis permettant de contrôler la précision des résultats au cours du calcul. Pour traiter des problèmes dynamiques, la variable d'espace est couplée à la variable de temps par l'intermédiaire "d'éléments espace-temps". A partir d'un maillage espace-temps grossier, un indicateur d'erreur permet de raffiner, au cours du calcul, automatiquement en espace et en temps dans les zones où la précision souhaitée n'est pas atteinte. Lorsque le raffinement automatique devient local, la stabilité de la méthode est conditionnée par le choix des conditions aux limites sur le domaine espace-temps. De plus, une étude de la stabilité basée sur la méthode énergétique a été menée. Un exemple unidimensionnel en dynamique linéaire illustre cette méthode.

ABSTRACT. This work proposes a finite element solver which allows to control the précision of the calculation. To model time dependent problems, the space and time variables are linked creating "space-time elements". With an initial space time mesh, an error indicator automatically refines during the calculation, whenever the error exceeds the requested precision. When local refinement occurs, the scheme's stability depends on the choice of the boundary conditions. Using the energy method, the scheme's stability is studied.

MOTS-CLÉS : raffinement local de maillage, maillage adaptatif, précision contôlée, indicateur d'erreur, méthode énergétique

KEYWORDS: local mesh refinement, adaptative mesh, controled precision, error indicator, energy method
# 1. Introduction

La capacité croissante des ordinateurs permet aux ingénieurs de développer des modèles éléments finis pouvant atteindre des millions de degrés de liberté. Les techniques numériques de résolution doivent donc s'adapter afin d'effectuer des calculs en des temps raisonnables. De plus, il est fréquent que les utilisateurs de la méthode EF soient plus intéressés par les ordres de grandeurs des quantités qu'ils traitent que par la précision des résultats obtenus. Cependant il est important de connaître cette précision afin de garantir le bon dimensionnement d'une pièce ou la sécurité de fonctionnement d'un mécanisme. C'est pourquoi ce travail s'intéresse à une méthode éléments finis (FEM) multi-échelle à précision contrôlée espace temps pour la dynamique. La STAR méthode (Space Time Automatic Refinement) est un solveur EF particulier qui donne automatiquement, à partir d'un maillage grossier et de la précision souhaitée, la solution et le maillage correspondant pour des problèmes de dynamique linéaire. Ce solveur est basé sur un concept d'élément "espace-temps" qui permet de coupler la discrétisation en espace à celle du temps au cours du calcul.

Différentes approches ont déjà été explorées pour répondre à ce type de questions. Les travaux de Brandt [BRA 77] sur les méthodes multi-grilles permettent d'améliorer considérablement l'efficacité de la résolution. D'autres travaux, notamment ceux de Adams [ADA 02], permettent d'adapter ces méthodes aux résolutions EF. Gravouil et Combescure [GRA 03] proposent de modéliser le problème en un nombre fixe de sous domaines. Chaque sous domaine possède sa propre discrétisation en espace et en temps et peut résoudre les équations de manière explicite ou implicite. Cependant cette approche oblige l'utilisateur à prédéfinir les sous domaines. D'autre part, Combe et al. [COM 02] présentent une approche permettant de contrôler les calculs en dynamique dans le cas d'impact. Le travail s'appuie sur des estimateurs d'erreur. Les travaux de Hughes et Hulbert [HUG 88] développent une méthode éléments finis espace temps adaptive. La méthode présentée ici peut également s'apparenter aux méthodes de pas multiples (Smolinski) [WU 00] et de sous-cyclage (Daniel) [DAN 03].

# 2. Raffinement Local espace-temps (STAR méthode)

Au cours du calcul, des niveaux successifs de maillage espace-temps de plus en plus fins sont créés. Sur le niveau *j* l'élément "espace-temps" $\Omega e^j$ se compose de deux éléments finis $\omega e^j$ et de deux instants $t_n^j$ et $t_{n+1}^j$ où $t_n^j$ représente l'instant *n* au niveau *j* et $t_{n+1}^j = t_n^j + \Delta t$ (figure 1). Pour un élément espace-temps l'axe horizontal *X* représente donc l'espace et l'axe vertical *T* le temps. Les éléments finis $\omega e^j$ peuvent être indifféremment 1, 2 ou 3D. Un élément espace-temps est donc défini sur le domaine $\Omega e^j = \omega e^j * [n, n+1]$. Le domaine espace-temps complet $\Omega^j$ est constitué de l'ensemble des éléments espace-temps $\Omega e^j$, ainsi $\Omega^j = \sum \Omega e^j$.

Par simplicité, le modèle de raffinement espace-temps est détaillé dans le cas d'un problème 1D. La méthode reste identique pour des problèmes 2 ou 3D. Le raffinement espace-temps est similaire au raffinement en espace présenté par les mêmes auteurs
[CAV 04]. Par la suite, la variable $X^j$ désigne le vecteur d'état $\{\underline{u}, \dot{\underline{u}}, \ddot{\underline{u}}\}^j$ représentant les vecteurs de déplacement, vitesse et accélération au niveau j.

<figure><img src="image_2.png" /><figcaption>**Figure 1.** Elément espace-temps $\Omega e^j$ pour un niveau j et ceux correspondants aux niveau suivant j + 1.</figcaption></figure>

Initialement, la solution du problème $X^1$, définie sur le premier maillage grossier du problème complet $\Omega^1$ est calculée. Pour initialiser le processus, un second domaine espace-temps $\Omega^2$, deux fois plus fin en espace et en temps, est créé permettant de calculer $X^2$. Ensuite, par l'intermédiaire d'un indicateur d'erreur, les solutions $X^1$ et $X^2$ sont comparées. Dans les zones où l'indicateur d'erreur est supérieur à la précision demandée le maillage espace-temps est raffiné créant un nouveau domaine espace-temps $\Omega^3$. Ce niveau peut être composé d'un maillage non connexe. En appliquant des CL appropriées à $\Omega^3$, $X^3$ est calculé. Le processus est ensuite répété entre le niveau 2 et 3 puis entre le niveau j et j + 1 jusqu'à ce que la solution atteigne la précision souhaitée sur tout le domaine espace-temps. A la fin du processus le maillage optimal ainsi que la solution correspondant à la précision demandée sont obtenus. Les équations [1] à [4] synthétisent cette démarche.

$$M_{final}\ddot{\underline{u}}_{final} + K_{final}\underline{u}_{final} = \underline{F}_{final} \quad [1]$$

$$\left[ \begin{matrix} [M^1] \\ [M^2] \\ \vdots \\ [M^k] \end{matrix} \right]_{final} \left\{ \begin{matrix} \frac{\ddot{u}^1}{\dot{u}^2} \\ \vdots \\ \frac{\ddot{u}^k}{\dot{u}^k} \end{matrix} \right\}_{final} + \left[ \begin{matrix} [K^1] \\ [K^2] \\ \vdots \\ [K^k] \end{matrix} \right]_{final} \left\{ \begin{matrix} \frac{\dot{u}^1}{u^2} \\ \vdots \\ \frac{\dot{u}^k}{u^k} \end{matrix} \right\}_{final} = \left\{ \begin{matrix} \frac{\underline{F}^1}{\underline{F}^2} \\ \vdots \\ \frac{\underline{F}^k}{\underline{F}^k} \end{matrix} \right\}_{final} \quad [2]$$

$[M^j]$ et $[K^j]$ pour $j = 1...k$ sont les matrices espace-temps de masse et de raideurs définies pour chaque pas de temps $t_n$ avec $n = 0...r$. L'équation [3] donne un exemple de matrice espace-temps.

$$[K^j] = \begin{bmatrix} [K_0^j] & & & \\ & [K_1^j] & & \\ & & \ddots & \\ & & & [K_r^j] \end{bmatrix} \quad [3]$$

La solution finale est obtenue suivant l'équation [4]. La solution $X^{j*}$ est la partie de $X^j$ définie sur le domaine $\Omega^j - \Omega^{j+1}$. Sur le niveau k le plus fin $X^{k*} = X^k$ et la solution du niveau 1 n'est pas prise en compte puisque $\Omega^2 - \Omega^1 = \emptyset$.
$$ \underline{X}_{final} = \underline{X}^k \cup \underline{X}^{k-1*} \cup \dots \cup \underline{X}^{j*} \cup \dots \cup \underline{X}^{2*} \quad [4] $$

Les indicateurs d'erreur en énergie et en déplacement sont utilisés. Pour résoudre numériquement le problème dynamique le schéma incrémental de Newmark de l'ac-célération en moyenne ($\gamma = 0.5$ et $\beta = 0.25$) est employé.

### 3. Stabilité de la STAR méthode

Lorsqu'une partie du maillage est raffinée localement des conditions aux limites (CL) particulières doivent être appliquées autour du nouveau domaine espace-temps. Ces CL doivent être choisies avec précautions puisqu'elles conditionnent la stabilité du calcul. Ainsi l'étude de la stabilité de la STAR méthode, basée sur la méthode énergétique présentée par Hughes et Liu [HUG 78], est menée.

Dans le cas de la STAR méthode il faut s'assurer que l'énergie des forces d'interface entre deux échelles se conserve pendant le calcul. Ainsi l'énergie d'interface entre un élément espace-temps du niveau $j$ ($E_{if}^j$) et les quatre éléments espace-temps correspondants du niveau $j + 1$ ($E_{if+1}^{j+1}$) doit être identique. Ces deux quantités sont décrites par les équations [5] et [6].

$$ E_{if}^j = \frac{1}{\Delta t} [\underline{\dot{u}}^{jT}]_n^{n+1} [\underline{F_i}^j]_n^{n+1} \quad [5] $$

$$ E_{if}^{j+1} = \sum_{i=1}^{2} \frac{2}{\Delta t} [\underline{\dot{u}}^{j+1T}]_{2n+i-1}^{2n+i} [\underline{F_i}^{j+1}]_{2n+i-1}^{2n+i} \quad [6] $$

où $[\underline{q}^{jT}]_n^{n+1} = \underline{q}_{n+1}^{jT} - \underline{q}_n^{jT}$ et $F_i$ représente les forces d'interface.

Pour montrer la conservation d'énergie, il faut développer le terme en vitesse et le terme des forces d'interface de l'équation [6]. Les vitesses à l'interface du niveau $j+1$ sont interpolées linéairement (équation [7]). Les forces d'interface sont développées via les équations d'équilibre aux instants $2n$ à $2n+2$ (équation [8]) et une interpolation linéaire des déplacements et des accélérations semblable à l'équation [7].

$$ \underline{\dot{u}}_{2n+i}^{j+1} = \frac{i}{2} \underline{\dot{u}}_{2n+2}^{j+1} + (1 - \frac{i}{2}) \underline{\dot{u}}_{2n}^{j+1} \quad \text{pour } i = 0, 1, 2 \quad [7] $$

$$ M^{j+1} \underline{\dot{u}}_{2n+i}^{j+1} + K^{j+1} \underline{u}_{2n+i}^{j+1} = \underline{F}_i^{j+1} \quad \text{pour } i = 0, 1, 2 \quad [8] $$

Une dernière hypothèse est à vérifier. L'opérateur de projection $Ps$ permettant de transférer à l'interface les champs du niveau $j$ au niveau $j+1$ doit être compatible avec l'assemblage EF des matrices, par exemple il faut que $M_n^j = Ps^T M_{2n}^{j+1} Ps$.
En utilisant toutes ces hypothèses il est possible de montrer que la quantité [6] équivaut à [5]. Ainsi, en assurant la continuité des déplacements, vitesses et accéléra-tions à l'interface entre deux échelles ainsi qu'une interpolation linéaire de ces quantités, aucune énergie n'est produite ou dissipée à l'interface. Sous ces hypothèses, la STAR méthode est stable.

## 4. Resultats

La figure 2 représente une poutre unidimensionnelle sollicitée en traction. Sa longueur $L$ est de $1 \text{ m}$ et sa section $S$ de $5 \cdot 10^{-3} \text{ m}^2$. Le matériau utilisé est un acier avec les caractéristiques suivantes : $E = 211 \cdot 10^9 \text{ Pa}$, $\rho = 7800 \text{ Kg.m}^{-3}$, $\nu = 0.3$.

<figure><img src="image_3.png" /><figcaption>**Figure 2.** Poutre unidimensionnelle en traction et son chargement.</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>**Figure 3.** Maillage espace-temps final obtenu avec un indicateur d'erreur (a) en énergie (précision demandée $10^{-4}$), (b) en déplacement (précision demandée $5 \cdot 10^{-3}$), (c) convergence pour des indicateurs en déplacement et en énergie</figcaption></figure>

Les figures 3 (a) et (b) représentent le maillage espace-temps final correspondant à la poutre 1D respectivement pour un indicateur d'erreur en énergie et en déplacement.
Il apparaît clairement que le maillage suit la propagation de l'onde dans la poutre. L'utilisation d'indicateur d'erreur permet de localiser le raffinement de maillage uniquement dans les zones sollicitées dynamiquement c'est à dire où l'onde se trouve.

La figure 3 (c) représente l'évolution de l'erreur en fonction des niveaux de maillage pour les deux indicateurs d'erreur. Cette convergence est comparée à un calcul de référence correspondant à un maillage uniforme. Ce graphique permet de conclure que le raffinement espace-temps n'altère pas la convergence du calcul.

## 5. Conclusion

Le couplage de la taille de maille en espace à la taille du pas de temps par l'intermédiaire des éléments espace-temps ainsi que l'utilisation d'indicateur d'erreur permet à la STAR méthode de créer automatiquement un maillage qui suit la propagation de l'onde dans le milieu. Grâce à l'utilisation d'un indicateur d'erreur, le maillage capte cette propagation d'onde tout en contrôlant la précision des résultats. De plus cette étude montre qu'aucune énergie n'est produite ou dissipée à l'interface entre deux niveaux assurant la stabilité de la méthode.

## 6. Bibliographie

* [ADA 02] ADAMS M., « Evaluation of three unstructured multigrid methods on 3D finite element problems in solids mechanics », *International Journal for Numerical Methods in Engineering*, vol. 55, 2002, p. 519-534.

* [BRA 77] BRANDT A., « Multi-level adaptive technique (MLAT) for fast numerical solution to boundary value problems », *Lecture Notes in Physics*, vol. 18, 1977, p. 82-89.

* [CAV 04] CAVIN P., GRAVOUIL A., LUBRECHT A. A., COMBESCURE A., « Efficient FEM calculation with predefined precision through automatic grid refinement », *Finite Element in Analysis and Design*, accepted in 2004.

* [COM 02] COMBE J. P., LADEVÈZE P., PELLE J. P., « Discretization error estimator for transient dynamic simulations », *Advances in Engineering Software*, vol. 33, 2002, p. 553-563.

* [DAN 03] DANIEL W., « A partial velocity approach to subcycling structural dynamics », *Computer Methods in Applied Mechanics and Engineering*, vol. 192, 2003, p. 375-394.

* [GRA 03] GRAVOUIL A., COMBESCURE A., « Multi-time-step and two-scale domain decomposition method for non-linear structural dynamics », *International Journal for Numerical Methods in Engineering*, vol. 58, 2003, p. 1545-1569.

* [HUG 78] HUGHES T. J. R., LIU W. K., « Implicit-explicit finite element in transient analysis : stability theory », *Journal of applied mechanics*, vol. 45, 1978, p. 371-374.

* [HUG 88] HUGHES T. J. R., HULBERT G. M., « Space-time finite element methods for elastodynamics : formulations and error estimates », *Computer Methods in Applied Mechanics and Engineering*, vol. 67, 1988, p. 339-363.

* [WU 00] WU Y., SMOLINSKI P., « A multi-time step integration algorithm for structural dynamics based on the modified trapezoidal rule », *Computer Methods in Applied Mechanics and Engineering*, vol. 187, 2000, p. 641-660.