# Co-simulation implicite/explicite multi-échelle en temps pour la dynamique du contact

Yvan Le Nôtre, Michael Brun, Clément Grenat, Anthony Gravouil

► To cite this version:

Yvan Le Nôtre, Michael Brun, Clément Grenat, Anthony Gravouil. Co-simulation implicite/explicite multi-échelle en temps pour la dynamique du contact. 15e colloque national en calcul des structures, Université Polytechnique Hauts-de-France [UPHF], May 2022, Hyères-les-Palmiers, France. hal-03717983

HAL Id: hal-03717983

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

Submitted on 8 Jul 2022

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
# Co-simulation implicite/explicite multi-échelle en temps pour la dynamique du contact

Y. Le Nôtre<fnref n="1" />, C. Grenat<fnref n="2" />, M. Brun<fnref n="3" />,
A. Gravouil<fnref n="1" />

**Résumé** — Le sujet porte sur la résolution numérique de problèmes dynamiques non-linéaires des structures incluant des phénomènes multi-échelles en temps présents pour des structures soumises à des conditions extrêmes. Ces phénomènes peuvent être simulés à l'aide des méthodes HATI et notamment la méthode BGC. La simulation de ces phénomènes induit des erreurs au cours du temps. L'objectif est l'amélioration des résultats obtenus par ces méthodes multi-échelles en temps.

**Mots clés** — multi-échelle, co-simulation, méthode BGC, dynamique, non-linéaire

# 1 État de l'art

## 1.1 Introduction

Il existe très peu de problèmes physiques ayant pour solution un résultat analytique. La majorité des problèmes est résolue de manière approximative en utilisant la méthode des éléments finis permettant une discrétisation spatiale, seule (permettant de résoudre des problèmes statiques), ou en couplant cette méthode avec une discrétisation en temps, pour résoudre des problèmes dynamiques. La discrétisation temporelle implique l'utilisation de schémas temporels. Plusieurs schémas temporels classiquement utilisés comme les schémas d'Euler, de Newmark, de Hilbert Hughes et Taylor (HHT) [1] et de Wood, Bossak et Zienkiewicz (WBZ) [2]. Ces schémas sont caractérisés par plusieurs paramètres tels que la précision, la consistance, l'ordre de convergence, la stabilité. Certains sont définis explicitement, d'autres implicitement, et le schéma doit être judicieusement choisi en fonction du problème à résoudre.

Les schémas explicites sont utilisés pour des phénomènes non-linéaires comme des chocs. Ces schémas effectuent efficacement les calculs à chaque pas de temps mais sont conditionnellement stables. Par conséquent, la sélection du pas de discrétisation doit être effectuée au sein de domaines prédéfinis afin d'assurer la propriété de stabilité de ces schémas. Les schémas implicites nécessitent plus de temps de calcul à chaque pas de temps. En revanche, ils ont l'avantage d'être inconditionnellement stables. Leurs pas de discrétisation ne sont donc pas bornés à un domaine d'utilisation prédéfini.

Un compromis existe : la co-simulation. En effet, cette méthode permet de combiner la résolution d'un problème avec plusieurs types de schémas temporels et de discrétisations spatiales en définissant plusieurs domaines dans l'espace étudié. La méthode de co-simulation proposée par Gravouil, Combescure et Brun [3] permet d'effectuer un calcul dynamique en utilisant, au sein du même calcul, des schémas implicites et explicites avec des pas de temps différents. Le bilan énergétique est proche de zéro et décroissant. Un des objectifs de ce travail est de prospecter des pistes d'amélioration pour annuler cette dissipation d'énergie à l'interface.

## 1.2 Formulation du problème mécanique

La résolution de problèmes est réalisée la plupart du temps par une approximation en espace et en temps. La première étape de la résolution d'un problème mécanique consiste à résoudre une discrétisation en espace. La méthode des éléments finis permet de trouver une solution approchée du problème de référence. Cette méthode conduit au système d'équations du mouvement semi-discrétisé (1) en espace

<footnote n="1">¹ LaMCoS, INSA Lyon, {yvan.le-notre,anthony.gravouil}@insa-lyon.fr</footnote>

<footnote n="2">² Section DTIMRD, Framatome Paris, clement.grenat@framatome.com</footnote>

<footnote n="3">³ LEM3, Université de Lorraine, michael.brun@univ-lorraine.fr</footnote>
pour une structure non-amortie [6] :

$$
\left\{
\begin{array}{l}
\forall t \in [0, T], \quad \mathbf{M}\ddot{u}(t) + f^{\text{int}}(u(t), \dot{u}(t)) = f^{\text{ext}}(t) \\
u(0) = u_0, \quad \dot{u}(0) = u_0
\end{array}
\right.
\qquad (1)
$$

Une fois le problème semi-discrétisé par la discrétisation spatiale (1), avec les matrices de masse **M**, de raideur **K**, d'amortissement **C**, les forces intérieurs $f^{\text{int}}$ et extérieures $f^{\text{ext}}$ ou simplement $f$ et les quantités cinématiques (déplacement $u$, vitesse $\dot{u}$, accélération $\ddot{u}$), il est nécessaire d'effectuer une seconde semi-discrétisation, temporelle, pour résoudre des problèmes de dynamique. Cette discrétisation est réalisée par le schéma de Newmark (2) avec notamment l'introduction des paramètres $\beta$ et $\gamma$ et le pas de temps $h$ utilisé pour la discrétisation temporelle.

$$
\begin{equation}
\left( \mathbf{M} + \gamma h \mathbf{C} + \beta h^2 \mathbf{K} \right) \ddot{u}_{n+1} = f_{n+1} - \underbrace{\mathbf{C} \left( \dot{u}_n + (1-\gamma)h\ddot{u}_n \right)}_{p \dot{u}_n} - \mathbf{K} \underbrace{\left( u_n + h\dot{u}_n + \left( \frac{1}{2} - \beta \right) h^2 \ddot{u}_n \right)}_{p u_n} \tag{2}
\end{equation}
$$

1.3 Méthode multi-échelle en temps

Les systèmes de plus en plus complexes en science de l’ingénieur nécessitent des discrétisations de plus en plus fines de systèmes linéaires ou non-linéaires. Cette augmentation constante, tant en termes de nombre de degrés de liberté résolus et qu’en termes de modélisations de plus en plus complexes, amène les ingénieurs et chercheurs à chercher des méthodes de résolutions innovantes. En présence de phénomènes non-linéaires raides, les méthodes de résolutions directes nécessitent des pas de discrétisation très fins en plus d’une méthode de résolution adaptée à la simulation de phénomènes non-linéaires, quand bien même, la majorité du domaine de résolution se comporte linéairement. Il en résulte une méthode de résolution non adaptée aux temps de calculs très longs. Les méthodes de décomposition de domaine, quant à elles, permettent de découper en plusieurs sous domaines un domaine global. L’intérêt d’une telle méthode réside dans son adaptabilité. Cette méthode permet à chaque sous domaine de posséder sa propre méthode de résolution, son maillage et son pas de temps. Les systèmes linéaires obtenus sont couplés par des conditions de continuité, puis résolus pour obtenir la solution globale. Cette méthode a l’avantage de pouvoir définir des méthodes de résolution adaptées aux domaines comportant de fortes irrégularités ou des comportements complexes.

Cette méthode de décomposition d'un domaine global en plusieurs sous domaines se décline en deux catégories : les méthodes de décomposition avec recouvrement et les méthodes de décomposition sans recouvrement. La démarche est décrite pour seulement deux sous-domaines mais peut être généralisée pour un plus grand nombre. Seule la décomposition sans recouvrement duale est présentée dans ce document, car utilisée dans notre étude.

Les méthode de décomposition sans recouvrement (de Schur) consistent à décomposer le domaine
global en plusieurs sous domaines disjoints. La continuité entre les sous domaines est assurée à l'interface
en utilisant des quantités cinématiques et l'équilibre des forces.

L'approche de Schur duale consiste, quant à elle, à coupler les sous domaines via une équation de
contrainte. Pour résoudre les systèmes couplés aux interfaces, des multiplicateurs de Lagrange y sont
ajoutés. Ces derniers peuvent être interprétés comme des forces internes fictives permettant aux équa-
tions de contraintes imposées aux interfaces d'être satisfaites. Dans notre cas, les équations de contraintes
considérées correspondent à des conditions cinématiques (continuité des déplacements, des vitesses ou
des accélérations). Cette approche duale est utilisée dans la méthode multi-échelle présentée.

Les méthodes multi-échelles requièrent la décomposition d'un système en plusieurs sous-domaines. Différentes discrétisations ou schémas numériques sont appliqués aux sous-domaines. Les équations de la mécanique (3) restent les mêmes (1). Dans le cas que l'on considère, le domaine est décomposé en deux sous-domaines nommés A et B (SDA et SDB).
$$ \left\{ \begin{aligned} \mathbf{M}^A \ddot{u}_m^A + \mathbf{C}^A \dot{u}_m^A + \mathbf{K}^A u_m^A &= f_n^A + \mathbf{L}^{AT} \Lambda_m, & [t_0, t_m] \\ \mathbf{M}^B \ddot{u}_m^B + \mathbf{C}^B \dot{u}_m^B + \mathbf{K}^B u_m^B &= f_n^B + \mathbf{L}^{BT} \Lambda_m, & [t_0, t_m] \end{aligned} \right. \quad (3) $$

$$ L^A \dot{u}_m^A + L^B \ddot{u}_m^B = 0 \quad (4) $$

La continuité est assurée par la résolution sous contrainte des sous-domaines couplés à l'aide des multiplicateurs de Lagrange. La continuité des vitesses (4) correspond à l'équation de contrainte [4] afin d'assurer la stabilité de l'algorithme de couplage implicite/explicite. Les termes $L^{AT} \Lambda_m$ et $L^{BT} \Lambda_m$ représentent les forces d'interfaces entre le sous-domaine A et le sous-domaine B.

Les méthodes multi-échelles utilisent la décomposition de domaine pour définir plusieurs échelles de temps au sein d'un même système. Elles sont souvent utilisées avec un sous-domaine exprimé de manière implicite, couplé avec un deuxième sous-domaine exprimé explicitement comme cela est illustré dans la Figure 1 :

<figure><img src="image_2.png" /><figcaption>FIGURE 1 – Deux sous-domaines avec différentes discrétisations temporelles</figcaption></figure>

$$ \left\{ \begin{aligned} \mathbf{M}^A \ddot{u}_m^A + \mathbf{C}^A \dot{u}_m^A + \mathbf{K}^A u_m^A &= f_m^A + \mathbf{L}^{AT} \Lambda_m^A, & [t_0, t_m] \\ \mathbf{M}^B \ddot{u}_j^B + \mathbf{C}^B \dot{u}_j^B + \mathbf{K}^B u_j^B &= f_j^B + \mathbf{L}^{BT} \Lambda_j^B, & [t_{j-1}, t_j] \end{aligned} \right. \quad (5) $$

Chacun des schémas numériques comporte ses propres paramètres : $\gamma$, $\beta$, $h$. Le rapport m entre le pas de temps grossier $h^A$ (échelle de temps macro) et le pas de temps fin $h^B$ (échelle de temps micro) est défini (6) et utilisé dans plusieurs équations :

$$ m = \frac{h^A}{h^B} \quad (6) $$

Dans l'équation (5), l'indice m représente l'échelle de temps grossière. L'indice j représente l'échelle de temps fine. Une fois la décomposition en plusieurs sous-domaine effectués, plusieurs méthodes multi-échelles en temps peuvent être définies.

La méthode BGC [3] définie par Brun, Gravouil et Combescure, généralise les résultats obtenus par Gravouil et Combescure (GC) [4] aux méthodes G-α (HHT, WBZ, CB). Ces deux méthodes sont basées notamment sur l'interpolation linéaire des vitesses (9) et des multiplicateurs de Lagrange (7) de l'échelle grossière vers l'échelle fine en temps et la continuité des quantités cinématique (8) sur l'échelle fine en temps :

$$ \Lambda_j^A = \left(1 - \frac{j}{m}\right) \Lambda_0^A + \frac{j}{m} \Lambda_m^A \quad (7) $$

$$ \mathbf{L}^A \dot{u}_j^A + \mathbf{L}^B \ddot{u}_j^B = 0 \Rightarrow \mathbf{L}^A \dot{u}_{j,link}^A + \mathbf{L}^B \ddot{u}_{j,link}^B = \mathbf{L}^A \dot{u}_{j,free}^A + \mathbf{L}^B \ddot{u}_{j,free}^B \quad (8) $$

$$ \dot{u}_{j,free}^A = \dot{u}_{m,free}^A \left(\frac{j}{m}\right) + \dot{u}_{0,free}^A \left(1 - \frac{j}{m}\right) \quad (9) $$

Cette relation permet de faire le lien à chaque itération entre les calculs se déroulant sur le sous-domaine B (indice j) et le calcul se déroulant sur le sous-domaine A (indice m). En effet, à chaque itération de la boucle j, doit être résolu l'équation linéaire ci-dessous (10) obtenue à partir du système d'équation du mouvement (5), pour calculer ensuite l'état des variables $u, \dot{u}, \ddot{u}$ :

$$ H \Lambda_j = -b_j \quad (10) $$
$$
\text{avec : } \begin{cases} H = \mathbf{L}^A \gamma^A h^A \tilde{\mathbf{M}}^{A-1} \mathbf{L}^{A^T} + \mathbf{L}^B \gamma^B h^B \tilde{\mathbf{M}}^{B-1} \mathbf{L}^{B^T} \\ b_j = \mathbf{L}^A \dot{u}_{j, \text{free}}^A + \mathbf{L}^B \dot{u}_{j, \text{free}}^B \end{cases}
$$

La matrice effective de masse $\tilde{\mathbf{M}}$ vient de l'équation (5) exprimée en fonction de l'accélération :

$$
\underbrace{(\mathbf{M}^A + \beta^A h^{A2} \mathbf{K}^A)}_{\tilde{\mathbf{M}}^A} \dot{u}_m^A = f_m^A - \mathbf{K}^A p u_m^A
$$

1.4 Bilan énergétique

Le bilan énergétique discrétisé (11) reprend celui utilisé pour les schémas de Newmark [6] en prenant en considération la décomposition en deux domaines (Figure 1). Ce bilan énergétique [3] sans l'amortissement et à chaque pas de temps peut être écrit.

$$
\left\{
\begin{array}{l}
\Delta W^A = \Delta W_{\text{cin},m}^A + \Delta W_{\text{int},m}^A + \Delta W_{\text{comp},m}^A - \Delta W_{\text{ext},m}^A - \Delta W_{\text{diss},m}^A \\
\Delta W^B = \Delta W_{\text{cin},m}^B + \Delta W_{\text{int},m}^B + \Delta W_{\text{comp},m}^B - \sum_{j=1}^{m} \Delta W_{\text{ext},j}^B - \sum_{j=1}^{m} \Delta W_{\text{diss},j}^B - \Delta W_{\text{interface}}
\end{array}
\right.
\qquad (11)
$$

$$
\Delta W^{\text{Total}} = \Delta W^{A} + \Delta W^{B}
$$

$$
\text{avec : } \Delta W_{\text{interface}} = \Delta u_m^{A^T} \left[ \frac{1}{2} \mathbf{L}_A^T (\Lambda_0 + \Lambda_m) + \left( \gamma_A - \frac{1}{2} \right) \mathbf{L}_A^T \Delta \Lambda_m \right] \\
+ \sum_{j=1}^{n} \left[ \Delta u_j^{B^T} \left[ \frac{1}{2} \mathbf{L}_B^T (\Lambda_{j-1} + \Lambda_j) + \left( \gamma_B - \frac{1}{2} \right) \mathbf{L}_B^T \Delta \Lambda_j \right] \right] \tag{12}
$$

Cette énergie représente le travail des forces d'interfaces des deux sous-domaines A et B. Lorsque
ce terme diffère de zéro, le schéma résultant n'est pas conservatif. Si l'énergie d'interface augmente, le
schéma peut alors devenir instable car l'on injecte de l'énergie dans le système. Dans le cas de la méthode
BGC, l'énergie d'interface (12) est strictement négative et globalement décroissante [3] :

2 Nouvelle approche

**2.1 Introduction**

La méthode multi-échelle en temps BGC est bien établie et a été validée pour de nombreuses applications industrielles. Cependant, la méthode actuelle présente toujours une dissipation numérique au niveau des interfaces entre les sous-domaines. Cet amortissement numérique peut notamment poser problèmes lors de la résolution de système présentant un amortissement du même ordre grandeur ou plus faible que l'amortissement généré par la méthode numérique. Afin d'améliorer les propriétés relatives à la dissipation d'énergie à l'interface, l'utilisation de l'accélération à la place de la vitesse, dans les équations de contraintes aux interfaces ou lors des interpolations entre les deux échelles de temps, est explorée.

**2.2 Condition cinématique de continuité en accélération**

La condition cinématique de continuité, ou aussi appelé collage à l'interface des sous domaines peut être exprimé en fonction du déplacement, de la vitesse ou de l'accélération. Les collages déplacement et en vitesse ayant déjà été explorés dans la thèse d'Anthony Gravouil [4], la méthode suivante appelée BLC, explore l'utilisation d'un collage en accélération :

$$
\mathbf{L}^A \ddot{u}_j^A + \mathbf{L}^B \ddot{u}_j^B = 0 \quad (13)
$$

L'équation de condensation (10) voit son expression également modifiée.
$$
\begin{cases}
H = \mathbf{L}^A h^A \tilde{\mathbf{M}}^{A-1} \mathbf{L}^{A^T} + \mathbf{L}^B h^B \tilde{\mathbf{M}}^{B-1} \mathbf{L}^{B^T} \\
b_j = \mathbf{L}^A \ddot{u}_{j, \text{free}}^A + \mathbf{L}^B \ddot{u}_{j, \text{free}}^B
\end{cases}
$$

Comme pour les vitesses (9) à l'échelle micro temps du sous-domaine A dans la méthode BGC, l'accélération au temps j du sous-domaine A n'existe pas. Elle est donc estimée en utilisant une interpolation linéaire :

$$
\ddot{u}_{j,free}^{A} = \ddot{u}_{m,free}^{A} \left(\frac{j}{m}\right) + \ddot{u}_{0,free}^{A} \left(1 - \frac{j}{m}\right) \quad (14)
$$

2.3 Condition cinématique de continuité en vitesse avec interpolation linéaire sur l'accé-
lération

Inspéré de la méthode BGC, un collage en vitesse avec une interpolation en accélération est une autre
voie à explorer. Les bonnes propriétés explorées avec le collage en accélération et les propriétés du col-
lage initialement effectué en vitesse ouvrent une nouvelle voie : coller en vitesse avec une interpolation
linéaire effectuée sur les termes de l'accélération. La condition de continuité des interfaces est assurée
par un collage en vitesse sur l'échelle fine en temps (8). La vitesse à l'interface au temps j est exprimée en
fonction de la précédente itération. Le pas (h), entre chaque instant, utilisé dans le schéma de Newmark
du sous-domaine A (macro) est celui du sous-domaine B (micro), schématisé dans la Figure 2. La vitesse
à l'instant j du sous-domaine A peut être calculée virtuellement par la relation de récurrence 15 qui relie
le terme 0 propre à l'échelle macroscopique et calculé à l'instant précédent au terme j propre au terme à
l'échelle microscopique :

<figure><img src="image_3.png" /><figcaption>FIGURE 2 – Même pas de temps utilisé tout au long du calcul</figcaption></figure>

$$
\mathbf{L}^A \ddot{u}_j^A = \mathbf{L}^A \left( \ddot{u}_0^A + (1-\gamma^A) \frac{h^A}{m} \ddot{u}_0^A + \sum_{k=1}^{j-1} \frac{h^A}{m} \ddot{u}_k^A + \frac{h^A}{m} \gamma^A \ddot{u}_j^A \right) \quad (15)
$$

L'expression de cette vitesse est réinjectée dans l'expression (8). Finalement l'équation (10) devient :

$$
\begin{equation}
\begin{aligned}
& \left[ \mathbf{L}^A \gamma^A \frac{h^A}{m} \tilde{\mathbf{M}}^{A-1} \mathbf{L}^{A^T} + \mathbf{L}^B \gamma^B h^B \tilde{\mathbf{M}}^{B-1} \mathbf{L}^{B^T} \right] \Lambda_j = - \left\{ \mathbf{L}^A \right. \\
& \qquad \left. \left[ \ddot{u}_0^A + (1-\gamma^A) \frac{h^A}{m} \ddot{u}_0^A \right. \right. \\
& \qquad \left. \left. + \sum_{k=1}^{j-1} \frac{h^A}{m} \left[ \ddot{u}_{m,free}^A + \left(1-\frac{k}{m}\right) \ddot{u}_{0,free}^A \right] \right. \right. \\
& \qquad \left. \left. + \gamma^A \frac{h^A}{m} \left[ \frac{j}{m} \ddot{u}_{m,free}^A + \left(1-\frac{j}{m}\right) \ddot{u}_{0,free}^A \right] \right. \right. \\
& \qquad \left. \left. + \sum_{k=1}^{j-1} \frac{h^A}{m} \tilde{\mathbf{M}}^{A-1} \mathbf{L}^{A^T} \Lambda_k \right] \right. \\
& \qquad \left. + \mathbf{L}^B \ddot{u}_{j,free}^B \right\}
\end{aligned}
\tag{16}
\end{equation}
$$

3 Résultats et conclusion

3.1 Condition des tests

La théorie, étudiée dans les sections précédentes, est appliquée dans un cas à une dimension. Le solide est une barre d'une longueur de 1 m, d'une section de 0,1 m², d'une masse volumique de 7800 kg/m³. Cette barre est divisée en deux sous domaines comportant chacun 50 degrés de liberté. Les forces extérieures, déplacements et vitesses sont nulles à l'instant initial. La barre est encastrée en son extrémité
de gauche et est soumise à une force sur son extrémité de droite, (Figure 4). La force évolue linéairement jusqu’à 0.03 s pour attendre une valeur de 1000 N, puis sa valeur est constante jusqu’à la fin de la simulation comme illustré sur la Figure 3.

<figure><img src="image_4.png" /><figcaption>FIGURE 3 – Force appliquée sur la barre</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>FIGURE 4 – Barre encastrée à gauche et soumise à une force à droite</figcaption></figure>

## 3.2 Résultats

La barre 1D, Figure 4, est divisée en deux sous domaines. Le sous-domaine A est discrétisé temporellement par un schéma de Newmark implicite ($\gamma^A = \frac{1}{2}$ et $\beta^A = \frac{1}{4}$) et le sous-domaine B est discrétisé temporellement par un schéma de Newmark explicite ($\gamma^B = \frac{1}{2}$ et $\beta^B = 0$). Les paramètres utilisés pour cette simulation sont : $h^A = 1e-6$s, $h^B = 1e-8$s (m=100), $T = 0,05$s.

Le bilan énergétique et la variation d'énergie à l'interface (12) des sous domaines sont des indicateurs permettant de savoir si les méthodes BGC et la nouvelle méthode sont conservatives. Une étude d'ordre de convergence des erreurs est aussi menée sur ces deux méthodes et tracée en fonction de $\Omega_A$ qui est le produit du pas de temps $h^A$ avec la pulsation $\omega^A$ de SDA pour un système harmonique. Les Figures 5 et 6 comparent les résultats des deux méthodes présentées dans ce document.

<figure><img src="image_6.png" /><figcaption>FIGURE 5 – Comparaison énergie à l'interface des sous-domaines pour les deux méthodes</figcaption></figure>

La convergence pour une approche multi-échelle en temps (m>1) est d'ordre 2 pour le bilan énergétique et la pseudo énergie. Ce résultat est très intéressant, car il n'existe pas actuellement dans la littérature de méthodes de collage sur l'échelle fine en temps d'ordre 2. En ajoutant la discrétisation spatiale, la dissipation de l'énergie est stabilisée à une valeur 10 fois plus petite que celle de la méthode BGC. Néanmoins, cette solution ne permet pas de tirer pleinement partie des méthodes multi-échelles car les incréments de temps utilisables dépendent de l'incrément de temps maximal à utiliser dans le domaine discrétisé spatialement finement (CFL, condition de stabilité dépendant de la discrétisation spa-
<figure><img src="image_7.png" /><figcaption>FIGURE 6 – Comparaison ordre de convergence bilan énergétique et pseudo énergie</figcaption></figure>

tiale et temporelle définissant un pas de temps critique à ne pas dépasser pour les schéma explicite). La conséquence de cette limitation est l'impossibilité d'avoir des incréments de temps supérieurs à la CFL avec un schéma d'intégration temporelle implicite et donc l'impossibilité d'avoir de grands écarts entre les échelles macro temps et micro temps. Or, c'est justement ces grands écarts qui permettent d'avoir des gains de temps sur les simulations mises en place avec un schéma d'intégration temporelle explicite couplé à un schéma implicite.

## 3.3 Conclusion

La méthode de calcul multi-échelle en temps GC fournit des résultats pertinents et validées dans plusieurs applications industrielles. L'amortissement numérique introduit par cette méthode dans le système ne permet pas d'offrir des résultats comparables à des calculs mono-échelles. La nouvelle approche proposée dans ce document explore l'usage de l'accélération dans la définition de la condition cinématique à l'interface et des interpolation temporelles liant les échelles de temps macro et micro. Cette dernière permet d'augmenter de un l'ordre de convergence de l'algorithme de couplage tout en ayant une variation de cette énergie oscillant autour de 0 au cours du temps. Cependant ce nouveau résultat n'est valable que sous certaines conditions de stabilité.

# 4 Remerciement

Je remercie l'Association Nationale de la Recherche et de la Technologie et mon employeur Frama-tome sans qui ce projet ne pourrait se dérouler. Mes remerciements s'adressent également au LaMCoS qui supporte ces travaux de recherche et qui a monté ce projet de recherche.

## Références

[1] Hilber, Hans M. and Hughes, Thomas J. R. and Taylor, Robert L., *Improved numerical dissipation for time integration algorithms in structural dynamics*, Earthquake Engng. Struct. Dyn., 5 : 283-292, 1977.

[2] Wood, W.L., Bossak, M. and Zienkiewicz, O.C., *An alpha modification of Newmark's method*, Int. J. Numer. Meth. Engng., 15 : 1562-156, 1980.

[3] Gravouil, A., Combescure, A. and Brun, M., *Heterogeneous asynchronous time integrators for computational structural dynamics*, Int. J. Numer. Meth. Engng, 102 : 202-232, 2015.

[4] A. Gravouil. *Méthode multi-échelle en temps et en espace avec décomposition de domaines pour la dynamique non-linéaire des structures*, Thèse de doctorat Physique Cachan, Ecole normale supérieure, 2000.
[5] Teddy Chantrait. *Approche multiéchelle en espace et en temps pour la prévision des endommagements dans les structures composites soumises à un impact de faible énergie*, Matériaux et structures en mécanique, INSA de Lyon, 2014.

[6] Hughes, T.J.R.. *The Finite Element Method : Linear Static and Dynamic Finite Element Analysis*, Dover Publications, 2012.