# Une méthode de couplage de schémas numériques en temps hétérogènes et multi-échelle en temps

Najib Mahjoubi, Anthony Gravouil, Alain Combescure, Nicolas Greffet

► To cite this version:

Najib Mahjoubi, Anthony Gravouil, Alain Combescure, Nicolas Greffet. Une méthode de couplage de schémas numériques en temps hétérogènes et multi-échelle en temps. 9e Colloque national en calcul des structures, CSMA, May 2009, Giens, France. hal-01408401

HAL Id: hal-01408401

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

Submitted on 4 Dec 2016

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
# Une méthode de couplage de schémas numériques en temps hétérogènes et multi-échelle en temps

N. Mahjoubi¹,², A. Gravouil¹, A. Combescure¹, N. Greffet²

¹ Université de Lyon, INSA-Lyon, LaMCoS, CNRS UMR5259, 18-20 rue des sciences, F-69621 Villeurbanne, France {najib.mahjoubi, anthony.gravouil, alain.combescure} @insa-lyon.fr

² EDF R&D (AMA/T62), LaMSID, EDF-CNRS, UMR2832, 1, avenue du Général de Gaulle F-92141 Clamart Cedex, France nicolas.greffet@edf.fr

**Résumé** — Le couplage de modèle apparaît comme un défi dans la compréhension des différents phénomènes physique mis en jeu. Les modèles numériques doivent tenir compte des différentes échelles spatiales et temporelles d'un problème. Dans ce contexte, sous réserve que les modèles en espace sont formalisés dans le cadre d'une méthode de décomposition de domaine de Schur duale (sans recouvrement), nous proposons un cadre unifié pour coupler de grandes classes de schémas numériques en temps avec des échelles de temps éventuellement très différentes.

**Mots clés** — Méthode de sous-domaines, Multi-échelle en temps, Schéma numérique, Décomposition de Schur Dual.

## 1 Introduction

Au cours de ces dernières années, le couplage de shéma d'intégration en temps a été un domaine de recherche actif. Belytchko et Mullen [1] ont été parmi les premiers à proposer une méthode de sous-cyclages qui soit multi-échelle en temps. Les méthodes de sous-cyclages ne permettent pas d'utiliser de grand ratio entre les différentes échelles de temps. Plus récemment, nous pouvons citer le travail de Casadei [2] dans le cadre de la dynamique transitoire explicite avec des échelles de temps locales éléments par éléments.

La méthode FETI proposée par Farhat [3] offre un cadre intéressant et relativement général pour exprimer le couplage dans un formalisme de Schur dual. Cette méthode a été initialement étendue à la dynamique transitoire des structures par Gravouil et Combescure (GC méthode) [4]. Ils ont proposé une méthode de couplage multi-échelles en temps des schémas de Newmark et montré une dissipation d'énergie à l'interface pour de grand ratio entre les échelles de temps. Prakash et Heljsmtad [5] ont proposé par la suite une méthode de couplage inspiré de la GC méthode qui permet d'assurer un travail de l'interface nul pour les schémas de Newmark.

Nous proposons dans ce travail une méthode générale de couplage de sous-domaine avec des schémas d'intégration en temps différents [6]. Cette méthode autorisera l'emploi d'échelle de temps différente sous-domaine par sous-domaine avec de large rapport entre les échelles. Cette méthode est élaborée en assurant un travail des interfaces entre les sous-domaines qui soit nul. Nous proposons pour cela un formalisme de couplage de schémas multi-échelles en temps qui permet de choisir le schéma d'intégration en temps et l'échelle de temps approprié sous-domaine par sous-domaine.
## 2 Problème de référence

On considère un corps occupant un domaine $\Omega$ de frontière $\partial\Omega$ au temps $t \in [0, T]$. Le matériau constitutif est supposé homogène et isotrope de masse volumique $\rho$. Les conditions aux limites sont constituées de vitesses $\mathbf{v}_v$ et d'efforts $\mathbf{t}_v$ imposés respectivement sur les frontières $\partial\Omega_v$ et $\partial\Omega_t$ tel que $\partial\Omega_t \cup \partial\Omega_v = \partial\Omega$, $\partial\Omega_t \cap \partial\Omega_v = \emptyset$. Le problème de référence est de trouver le couple $(\mathbf{v}, \underline{\underline{\sigma}})$ définis sur $\Omega \times [0, T]$. L'intégration sur l'intervalle de temps $[0, T]$ du Principe des Puissances Virtuelles sur l'intervalle de temps $[0, T]$ s'écrit :

$$ \int_0^T \left[ \int_{\Omega_*} \mathbf{v}_* \left[ \underline{\underline{div}}(\underline{\underline{\sigma}}) - \mathbf{f} - \rho \frac{\partial \mathbf{v}_*}{\partial t} \right] d\Omega \right] dt = \mathbf{0} \quad \forall \mathbf{v}_* \in \mathcal{V}^* \quad (1) $$

Dans cet article, les espaces fonctionnels ne sont pas détaillés. On suppose les quantités correspondantes suffisamment régulieres et à énergie finie. Il est possible de faire apparaître explicitement les conditions aux limites sur la frontière $\partial\Omega_v$ dans le formulation faible finale (1). On considère par la suite des conditions aux frontières linéaires de la forme $\mathbf{g}(\mathbf{v}) = \mathbf{v} - \mathbf{v}_d$.

$$ \int_0^T \left[ \int_{\Omega_*} \mathbf{v}_* \left[ \underline{\underline{div}}(\underline{\underline{\sigma}}) - \mathbf{f} - \rho \frac{\partial \mathbf{v}_*}{\partial t} + \boldsymbol{\lambda} \right] d\Omega \right] dt + \int_0^T \left[ \int_{\Omega_v} \mathbf{g}(\mathbf{v}_*) \boldsymbol{\lambda}_* d\Omega \right] dt = \mathbf{0} \quad (2) \\ \forall \mathbf{v}_* \in \mathcal{V}^* \text{ et } \forall \boldsymbol{\lambda}_* \in \Lambda^* $$

Le champ $\boldsymbol{\lambda}$ est un champ de multiplicateurs de Lagrange. L'introduction de façon usuelle de la méthode des éléments finis dans la formulation (2) ainsi que la discrétisation de l'intervalle de temps $[0, T]$ permet d'écrire :

$$ \int_{t_{n+1}}^{t_n} \mathbf{v}_* [\mathbf{M}\mathbf{a} + \mathbf{f}^{int} - \mathbf{f}^{ext} + \mathbf{L}^T \boldsymbol{\lambda}] dt + \int_{t_{n+1}}^{t_n} \boldsymbol{\lambda}_* [\mathbf{L}\mathbf{v}_* - \mathbf{v}_d] dt = \mathbf{0} \quad \forall \mathbf{v}_* \in \mathcal{V}^* \quad \forall \boldsymbol{\lambda}_* \in \Lambda^* \quad (3) $$

où $\mathbf{L}$ est l'opérateur de restriction sur l'interface du domaine. La formulation (3) donne une formulation au sens faible de l'équilibre sur l'intervalle de temps $[t_n, t_{n+1}]$. En conséquence, il existe un instant $t_{n+\xi} \in [t_n, t_{n+1}]$ avec $\xi \in [0, 1]$ tel que l'équilibre est vérifié au sens fort (équation (4)a). Nous supposons que les multiplicateurs de Lagrange sont constants par intervalle de temps $[t_n, t_{n+1}]$ ce qui permet d'écrire une condition sur le contour $\Omega_v$ (équation (4)b). Ainsi étant donné que la formulation faible (3) est vérifiée $\forall \mathbf{v}_*$ et $\forall \boldsymbol{\lambda}_*$, on peut écrire :

$$ \left\{ \begin{array}{l} \mathbf{M}\mathbf{a}_{n+\xi} + \mathbf{f}_{n+\xi}^{int} + \mathbf{L}^T \boldsymbol{\lambda} = \mathbf{f}_{n+\xi}^{ext} \\ \int_{t_{n+1}}^{t_n} [\mathbf{L}\mathbf{v} - \mathbf{v}_d] dt = \mathbf{0} \end{array} \right. \quad (4) $$

Ce paragraphe a permis de définir l'écriture au sens faible du problème de référence en dynamique. A partir de là, on propose au paragraphe 3 de donner la formulation faible d'un problème dynamique transitoire décomposé en sous-domaines avec leur propre schéma numérique en temps et leur propre échelle de temps.

## 3 Formulation faible du problème de sous-domaine

Dans ce paragraphe, on considère un domaine de référence $\Omega$ divisé en deux sous-domaines $\Omega_\Lambda$ et $\Omega_B$. Cette partition crée une interface entre les sous-domaines. On considère que cette interface est maillée sur chaque sous-domaine de façon compatible en espace. Nous proposons une méthode de décomposition de domaine qui permet d'utiliser des schémas d'intégration en temps différents
<figure><img src="image_2.png" /><figcaption>Figure 1 – (a) Même échelle de temps (b) Différentes échelles de temps</figcaption></figure>

sur chaque sous-domaine et/ou des échelles de temps différentes. La condition de continuité à l'interface correspondant à l'équation (4)b s'écrit comme :

$$ \int_{t_{n+1}}^{t_n} [\mathbf{L}_A \mathbf{v}^A + \mathbf{L}_B \mathbf{v}^B] \, dt = \mathbf{0} \qquad (5) $$

La formulation faible de ce problème peut se présenter sous la forme :

$$ \left\{ \begin{array}{l} [\mathbf{M}_A \mathbf{a}^A + \mathbf{f}_{int}^A - \mathbf{f}_{ext}^A + \mathbf{L}_A^T \mathbf{\lambda}]_{n+\xi_A} = \mathbf{0} \\ [\mathbf{M}_B \mathbf{a}^B + \mathbf{f}_{int}^B - \mathbf{f}_{ext}^B + \mathbf{L}_B^T \mathbf{\lambda}]_{n+\xi_B} = \mathbf{0} \\ \int_{t_{n+1}}^{t_n} [\mathbf{L}_A \mathbf{v}^A + \mathbf{L}_B \mathbf{v}^B] \, dt = \mathbf{0} \end{array} \right. \qquad (6) $$

Cette méthode propose un formalisme général de couplage des schémas numériques en temps. En effet, les équations (6)a et (6)b ont été établies sans émettre d'hypothèse relative à une quelconque famille de schéma numérique. Ainsi le choix du schéma adéquat est laissée à l'appréciation de l'utilisateur.

Il est également possible d'envisager des échelles de temps différentes sur chaque sous-domaine (Figure 3) en assurant un travail d'interface nul sur l'échelle grossière en temps. En assurant la condition (6)c, la méthode de couplage proposée assure que le travail des efforts d'interface soit nulle. Dans les paragraphes (3.1) et (3.2), la discrétisation de la relation (6)c est présentée respectivement dans le cas mono-échelle et multi-échelle.

On pose la notation $\bar{\mathbf{x}}$ pour définir la moyenne sur l'intervalle de temps d'une quantité vectorielle $\mathbf{x}$ par :

$$ \bar{\mathbf{x}} = \frac{1}{h} \int_{t_n}^{t_{n+1}} \mathbf{x} \, dt \qquad (7) $$

## 3.1 Même échelle de temps

L'intégration en temps est réalisée avec des pas de temps égaux $h_A = h_B$. Nous considérons que les multiplicateurs de Lagrange sont constants sur un intervalle de temps $[t_n, t_{n+1}]$. Cela permet d'assurer que l'intégrale de l'équation (6c) soit vérifiée. En utilisant la notation (7), l'équation (6c) peut être discrétisée comme suit :

$$ \int_{t_{n+1}}^{t_n} [\mathbf{L}_A \mathbf{v}^A + \mathbf{L}_B \mathbf{v}^B] \, dt = h_A \mathbf{L}_A \bar{\mathbf{v}}^A + h_B \mathbf{L}_B \bar{\mathbf{v}}^B = \mathbf{0} \qquad (8) $$

## 3.2 Différentes échelles de temps

Dans le cas où des échelles de temps différentes sont utilisées, on considère que les multiplicateurs de Lagrange sont constants et égaux sur l'échelle fine en temps du sous-domaine B
(Figure 3b) entre les instants correspondants à ceux de l'échelle macro en temps du sous-domaine A. Nous supposons également que les multiplicateurs de Lagrange de l'échelle fine sont égaux à ceux de l'échelle grossière en temps. Ainsi l'équation (6c) est vérifiée et peut s'exprimer en utilisant la notation (7) :

$$ \int_{t_{n+1}}^{t_n} [\mathbf{L}_A \mathbf{v}^A + \mathbf{L}_B \mathbf{v}^B] dt = h_A \mathbf{L}_A \bar{\mathbf{v}}^A + h_B \sum_{j=1}^{m} \mathbf{L}_B \bar{\mathbf{v}}_j^B = \mathbf{0} \quad (9) $$

Pour plus de détails, le lecteur peut se référer à [6].

## 4 Résultats

Afin d'illustrer cette méthode de couplage, le calcul dynamique d'une barre en traction est présenté. La barre est divisée en quatre sous-domaines (A, B, C, et D) de longueur $l = 1m$ (Figure 2). La structure est constituée d'un matériau homogène, isotrope, linéaire élastique avec un module de Young $E = 210 \cdot 10^9 Pa$, une densité volumique de $\rho = 7800 kg \cdot m^{-3}$ et un coefficient de Poisson de $\nu = 0.3$. La barre est soumise à un effort de traction $F_{max} = 1 \cdot 10^8 N$ dépendant du temps dont l'évolution est donnée Figure 3. Les quatre sous-domaines sont intégrés avec

<figure><img src="image_3.png" /><figcaption>Figure 2 – Modèle de la barre soumise à un effort de traction</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>Figure 3 – Evolution temporelle de l'effort de traction</figcaption></figure>

quatre schémas numériques en temps différents (Newmark [7], HHT [8], Simo [9], Krenk [10]) et des échelles de temps différentes. Les pas de temps sur chacun des sous-domaines sont tels que $h_j = \frac{H}{m_j}, j \in \{A, B, C, D\}$ et $H = 10^{-5} s$. Les paramètres des différents schémas numériques ainsi que les ratios entre échelles de temps sont rappelés dans le Tableau 1.

- Sous-domaine A : Simo

- Sous-domaine B : Krenk $\alpha = 0.176$

- Sous-domaine C : HHT $\alpha = -0.1$
– Sous-domaine D : Newmark $γ = 0.5, β = 0$

<table><caption>Tableau 1 – Ratio entre les différentes échelles de temps</caption><thead><tr><th>Sous-domaine</th><th>A</th><th>B</th><th>C</th><th>D</th></tr></thead><tbody><tr><td>$m_j$</td><td>1</td><td>250</td><td>500</td><td>1000</td></tr></tbody></table>

Les déplacements, vitesses et accélération d'un noeud à l'interface entre les sous-domaines B et C sont présentés respectivement sur les Figure 4a, 4b 4c. Les solutions sont comparées à un calcul de référence obtenu à partir d'un calcul mono-domaine. On peut noter la bonne corrélation des résultats avec un calcul de référence. Un bilan d'énergie est présenté à la Figure 4d afin de vérifier que la travail des efforts d'interfaces soit nul. Conformément à la théorie, nous obtenons un travail d'interface nulle d'un point de vue numérique de l'ordre de $10^{-10} J$ c'est à dire quinze ordres de grandeurs inférieurs à la valeur maximale des énergies cinétique et potentielle.

<figure><img src="image_5.png" /><figcaption>Figure 4 – (a) Déplacement - (b) Vitesse - (c) Accélération</figcaption></figure>

## 5 Conclusion

Nous proposons ici une méthode générale de couplage de sous-domaines autorisant l'emploi de schémas d'intégrations en temps différents et des échelles de temps différentes. Cette méthode a été développée à partir d'une formulation faible de l'équilibre afin d'assurer une dissipation d'énergie nulle à l'interfaces entre les sous-domaines. Ainsi, la consistance et la stabilité globale sont gouvernées uniquement par la consistance et la stabilité de chaque sous-domaines sans être perturbées par le travail des efforts d'interfaces. Cette méthode est adaptée aux problèmes de simulation d'impact (crash, séisme) lorsque différentes échelles de temps sont mises en jeu localement. Cette méthode peut être étendue aux problèmes de couplage de modèles, atomistique-continu par exemple, dans lesquels des échelles de temps très différentes doivent être couplées. L'extension de cette méthode aux couplages de problèmes multi-physiques comme le couplage fluide-structure dans une stratégie de couplage de code [12] présente également d'importantes perspectives.
Références

[1] Belytschko T, Mullen R. Explicit integration of structural problems. *Finite Elements in Non-linear Mechanics*, Bergen P, et al (eds), vol. 2, 1977 ; 697-720.

[2] Casadei F, Halleux JP. Binary spatial partitioning of the central-difference time integration scheme for explicit fast transient dynamics. *International Journal Numerical Methods in Engineering* 2008 ; early view

[3] Farhat C, Mandel J. The two-level FETI method for static and dynamic plate problems Part I : An optimal iterative solver for biharmonic systems. *Computer methods in applied mechanics and engineering* 1998 ; **155**:129-151.

[4] A. Gravouil, A. Combescure. A multi-time-step explicit-implicit method for non-linear structural dynamics. *International Journal for Numerical Methods in Engineering* 2001 ; **50**:199-225.

[5] A. Prakash, KD. Hjelmstad. A FETI-based multi-time-step coupling method for Newmark schemes in structural dynamics. *International Journal Numerical Methods in Engineering* 2004 ; **61**:2183-2204.

[6] N. Mahjoubi, A. Gravouil, A. Combescure, N. Greffet A FETI-based multi-time-step coupling method for Newmark schemes in structural dynamics. *Computer methods in applied mechanics and engineering* submitted.

[7] Newmark NM. A method of computation for structural dynamics. *Journal of the Engineering Mechanics Division* (ASCE) 1959 ; **85**: 67-94.

[8] Hilber HM, Hughes TJR, Taylor RL. Improved numerical dissipation for time integration algorithms in structural dynamics. *Earthquake Engineering and Structural Dynamics* 1977 ; **5**: 283-292.

[9] Simo J.C., Wong K.K. Unconditionally stable algorithms for rigid body dynamics that exactly preserve energy and momentum. *International Journal for Numerical Methods in Engineering* 1991 ; **31**:19-52.

[10] Krenk S. Extended state-space time integration with high-frequency energy dissipation. *International Journal for Numerical Methods in Engineering* 2008 ; **73**:1767-1787.

[11] Verlet L., Computer "Experiments" on Classical Fluids. I. Thermodynamical Properties of Lennard-Jones Molecules. *Phys. Rev.* 1967 ; **159**: 201-214

[12] Code_Aster. Reference book. EDF R&D. www.code-aster.org