# Une méthode variationnelle en temps pour le couplage de schémas hétérogènes en dynamique transitoire

Alain Combescure, Najib Mahjoubi, Anthony Gravouil, Nicolas Greffet

► To cite this version:

Alain Combescure, Najib Mahjoubi, Anthony Gravouil, Nicolas Greffet. Une méthode variationnelle en temps pour le couplage de schémas hétérogènes en dynamique transitoire. 9e Colloque national en calcul des structures, May 2009, Giens, France. hal-01391196

HAL Id: hal-01391196

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

Submitted on 3 Nov 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 variationnelle en temps pour le couplage de schémas hétérogènes en dynamique transitoire

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

¹ Universite de Lyon, Insa de Lyon, LaMCoS UMR CNRS 5259
18-20 Allée des sciences 69621 Villeurbanne Cedex
{alain.combescure,najib.mahjoubi,anthony.gravouil}@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é** — Cet article présente très succinctement les derniers résultats de l’équipe sur le couplage de codes en dynamique. La méthode proposée est très générale et elle est basée sur une vision faible de l’intégration temporelle des équations d’équilibre dynamique. Cette manière de comprendre les algorithmes de calcul en dynamique permet de mettre en place une stratégie de couplage qui assure que le couplage ne perturbe en aucune façon les bilans d’énergie des sous domaines assemblés. Il en résulte que la méthode proposée permet d’assembler des sous domaines ou des codes différenets utilisant leur propre intégrateur temporel, leur propre maillage, et leur propre pas de temps. Un exemple est proposé pour illustrer le propos.

**Mots clés** — Dynamique transitoire, couplages de sous domaines, éléments finis.

## 1 Introduction

Le couplage de sous domaines en dynamique est un sujet de recherche actif ces dernières années. Après une période de doutes sur la faisabilité de ce type de concept, l’arrivée de la méthode FETI ([19],[20]) a permis de comprendre l’importance des méthodes de Shurr duales pour le couplage de domaines. Les travaux d’A Gravouil ([6], [24]) dans les années 90 ont montré qu’un couplage de schémas de Newmark assurait une parfaite conservation de l’énergie mais aussi la stabilité du couplage si les collages étaient faits en vitesse pour des pas de temps identiques dans les sous domaines. Si les pas de temps sont différents entre sous domaines la stabilité est toujours garantie par le schéma proposé mais l’interface consomme un peu d’énergie (quelques pour cent de l’énergie mécanique) si bien que les schémas assemblés perdent un ordre de convergence. Prakash ([33]) a proposé récemment une amélioration du schéma qui n’a plus cet inconvénient mais qui est complexe et lourde de mise en oeuvre car on doit conserver les états de tous les petits pas de temps pendant l’intégration sur un grand pas de temps. De plus la méthode de Prakash ne s’applique qu’aux schémas de Newmark. Ce papier donne un cadre général qui permet de coupler des sous domaines avec des schémas quelconques à un pas. Le papier est organisé en trois parties. La première donne le cadre théorique de la formulation. La seconde expose la mise en oeuvre sur quelques algorithmes et précise les schémas numériques. La troisième développe un exemple.
## 2 Formulation faible de l'équilibre dynamique

Nous décrivons ici la formulation faible du principe des puissance virtuelles pour les systèmes mécaniques.

### 2.1 Un seul domaine

Intégrons sur $[0, t]$ le principe des puissance virtuelles. Il vient :

$$W^*(t) = \int_0^t \left( f_{\Omega} \underline{v^*}(P, \tau)^T \left[ \underline{\underline{\underline{\sigma(P,t)}}} - \underline{\underline{f^{vol}(P,t)}} - \rho \underline{\underline{\ddot{u}(P, \tau)}} \right] d\Omega(P) \right) d\tau = 0 \quad \forall \underline{v^*} \quad (1)$$

Dans l'équation (1) $\Omega$ est le domaine, $\underline{v^*}(P, \tau)$ la vitesse virtuelle au point P, $\tau$ le temps considéré, $\underline{\underline{\sigma}}$ le tenseur de contraintes, $\rho$ la masse volumique, $\ddot{u}$ l'accélération, $\underline{\underline{f^{vol}}}$ les forces de volume. Si on introduit maintenant une discrétisation spatiale (éléments finis par exemple) on trouve :

$$W^*(t) = \int_0^t \underline{v^*}(\tau)^T [\mathbf{F}^{\text{int}}(\tau) - \mathbf{M}\ddot{u}(\tau) - \mathbf{F}^{\text{ext}}(\tau)] d\tau \quad (2)$$

Dans l'équation (2) $\mathbf{F}^{\text{int}}$ est le vecteur des forces internes, $\mathbf{F}^{\text{ext}}$ les forces externes, $\mathbf{M}$ la matrice masse. Si on applique des contraintes cinématiques des forces de liaison $\lambda(\tau)$ sont ajoutées. L'équation(2) devient :

$$W^*(t) = \int_0^t (\mathbf{v}^*(\tau))^T [\mathbf{F}^{\text{int}}(\tau) - \mathbf{M}\ddot{u}(\tau) - \mathbf{F}^{\text{ext}}(\tau) + \lambda(\tau)] + (\mathbf{v}(\tau))^T \lambda^*(\tau) d\tau \quad (3)$$

Observons que les efforts de liaison jouent le même rôle que les efforts extérieurs et qu'ils travaillent avec les vitesses virtuelles. Regardons maintenant l'intégration temporelle sur un pas de temps fini $h$ ($h = t_{n+1} - t_n$) de l'équation(3). Appelons $\Delta(x)$ l'incrément d'une variable $x$ sur le pas de temps $h$. La contribution du pas de temps $h$ à l'équation(3) est :

$$\Delta W^* = \int_{t_n}^{t_{n+1}} (\mathbf{v}^*(\tau))^T [\mathbf{F}^{\text{int}}(\tau) - \mathbf{M}\ddot{u}(\tau) - \mathbf{F}^{\text{ext}}(\tau) + \lambda(\tau)] + \lambda^*(\tau)^T \mathbf{v}(\tau) d\tau \quad (4)$$

Observons que si l'on suppose que les fonctions du temps $\mathbf{F}^{\text{int}}(\tau)$, $\ddot{u}(\tau)$, $\mathbf{F}^{\text{ext}}(\tau)$, $\mathbf{v}(\tau)$ et $\lambda(\tau)$ sont constantes sur le pas de temps, vérifier l'équilibre est équivalent à minimiser l'incrément de puissance virtuelle $\Delta W^*$. Les équations d'équilibre dynamique ne sont pas vérifiées de façon forte sur tout l'intervalle de temps $h$, mais seulement au sens faible. C'est cependant l'hypothèse qui est faite pour tous les schémas d'intégration numérique habituels (Newmark ([32]), HHT ([4]), Simo([35]), Krenk ([30]) et Verlet ([2]) (formulation en vitesse (qui est d'ailleurs identique au schéma de Newmark dynamique explicite))) Pour les efforts de liaison $\lambda$ nous allons faire l'hypothèse très simple qu'ils sont eux aussi constants sur le pas de temps $h$. Le théorème de la valeur moyenne permet d'affirmer que si les fonctions sont continues en temps il existe un instant (inconnu) $t_{n+\theta} \in [t_n, t_{n+1}]$ tel que l'équilibre est exactement vérifié. Cet instant est fonction du schéma numérique choisi. Appelons maintenant $\bar{x}$ la valeur moyenne de la variable $x$. Il vient :

$$\overline{\Delta W^*} = \left[ \int_{t_n}^{t_{n+1}} \mathbf{v}^*(\tau)^T d\tau \right] \left[ \mathbf{F}^{\text{int}} - \mathbf{M}\ddot{u} - \mathbf{F}^{\text{ext}} + \mathbf{L}\bar{\lambda} \right] + \left[ \int_{t_n}^{t_{n+1}} \mathbf{v}(\tau)^T d\tau \right] \left( \mathbf{L}\bar{\lambda}^* \right) \quad (5)$$

### 2.2 Couplage de deux sous domaines

Nous allons maintenant coupler deux sous domaines A et B collés par une interface $\Gamma$. Les forces de liaison imposeront l'égalité des vitesses des sous domaines sur l'interface. En posant :

$$\begin{aligned} \mathbf{v}_{A}^{\gamma} &= \mathbf{L}_{A} \mathbf{v}_{A} \\ \mathbf{v}_{B}^{\gamma} &= \mathbf{L}_{B} \mathbf{v}_{B} \end{aligned} \quad (6)$$
L'incrément de puissance virtuelle sur le pas de temps h s'écrit :

$$ (\overline{\Delta W^*})^{tot} = (\overline{\Delta W^*})^A + (\overline{\Delta W^*})^B + (\overline{\Delta W^*})^{link}_{AB} \quad (7) $$

La variation de puissance associée aux liaisons s'écrit :

$$ (\overline{\Delta W^*})^{link}_{AB} = h(\mathbf{L_A V_A} + \mathbf{L_B V_B})^T \overline{\lambda^*}_{AB} \quad (8) $$

On voit immédiatement sur l'équation (8) que si les vitesses vérifient les liaisons la variation de travail d'interface est nulle et réciproquement. On en conclut que cette formulation assure que l'énergie d'interface est rigoureusement nulle. L'énergie du système est donc exactement égale à la somme des énergies des parties.

Cette idée est la clé de la méthode de couplage temporelle proposée.

- remarque 1 : Aucune hypothèse n'a été introduite de schéma numérique n'a encore été introduite : il en résulte que cette méthode de couplage est valable quel que soit le schéma d'intégration temporel pourvu qu'il permette de passer de l'état $n$ à l'état $n+1$.

- remarque 2 : La stabilité du système est gouvernée par la seule stabilité des parties, car la variation d'énergie d'interface est rigoureusement nulle.

- remarque 3 : Aucune hypothèse n'est faite sur la variation temporelle des forces de liaison. Elles sont constantes par morceaux sur chaque intervalle de temps.

- remarque 4 : La condition de liaison choisie est telle que les histoires de déplacement des deux côtés de l'interface peuvent être différentes mais que leur incrément fini sur le pas de temps sera indentique. En effet on a :

$$ \int_{t_n}^{t_{n+1}} (\mathbf{L_A V_A}(\tau) + \mathbf{L_B V_B}(\tau)) d\tau = 0. \quad (9) $$

- remarque 5 : L'extension de la méthode proposée à un nombre quelconque de sous domaines est évidente.

## 3 Schémas d'intégration temporelle

Nous allons ici présenter une forme générale pour les schémas HHT, Newmark, SIMO, Krenk et Verlet. Nous avons choisi de formuler le problème en incrément de vitesse. L'équation incrémentale dynamique avec liaisons au pas de temps n+1 s'écrit :

$$ \left(\frac{1}{h\gamma}\mathbf{M} + h\frac{\beta}{\gamma}(1+\alpha)\mathbf{K}\right)\Delta\mathbf{v} + \mathbf{L}^T\bar{\lambda} = \eta\mathbf{f}_{n+1} + \mu\mathbf{f}_n \\ -(\eta+\mu)\mathbf{K}\mathbf{u}_n - h(1+\alpha)\mathbf{K}\mathbf{v}_n \\ + (\eta-\mu)\left(\left(\frac{1}{\gamma}-1\right)\mathbf{M} + (1+\alpha)(\beta-\frac{\gamma}{2})h^2\mathbf{K}\right)\mathbf{a}_n \quad (10) \\ \mathbf{L}\Delta\mathbf{v} = \mathbf{0} $$

Les 5 paramètres nécessaires pour tous les modèles sont définis dans la table(1) : Définissons la raideur équivalente $\mathbf{K}_*$ comme :

$$ \mathbf{K}_* = \mathbf{M} \frac{1}{\gamma h} + (1+\alpha) \frac{\beta}{\gamma} h \mathbf{K} \quad (11) $$

le système d'équations (22) s'écrit alors :

$$ \begin{bmatrix} \mathbf{K}_* & \mathbf{L}^T \\ \mathbf{L} & \mathbf{0} \end{bmatrix} \begin{bmatrix} \Delta \mathbf{U} \\ \bar{\lambda} \end{bmatrix} = \begin{bmatrix} \mathbf{F} \\ \mathbf{0} \end{bmatrix} - \begin{bmatrix} \mathbf{N} & \mathbf{0} \\ \mathbf{0} & \mathbf{0} \end{bmatrix} \begin{bmatrix} \mathbf{U}_n \\ \mathbf{0} \end{bmatrix} \quad (12) $$
<table><caption>Tableau 1 – Valeur des paramètres pour l'écriture des algorithmes dans un formalisme unique</caption><thead><tr><th>integrator type</th><th>HHT</th><th>Newmark</th><th>SIMO</th><th>Verlet</th><th>Krenk</th></tr></thead><tbody><tr><td>α</td><td>α</td><td>0.</td><td>0.</td><td>0.</td><td>0.</td></tr><tr><td>β</td><td>β</td><td>β</td><td>1/4</td><td>0.</td><td>1/4</td></tr><tr><td>γ</td><td>γ</td><td>γ</td><td>1/2</td><td>1/2</td><td>1/2</td></tr><tr><td>η</td><td>1 + α</td><td>1</td><td>1</td><td>1</td><td>1</td></tr><tr><td>μ</td><td>-α</td><td>0</td><td>1</td><td>0</td><td>1</td></tr></tbody></table>

avec

$$ \mathbb{K}_* = \begin{bmatrix} \mathbf{K}_* & \mathbf{0} & \mathbf{0} \\ -\frac{\beta h}{\gamma} \mathbf{I} & \mathbf{I} & \mathbf{0} \\ -\frac{1}{\gamma h} \mathbf{I} & \mathbf{0} & \mathbf{I} \end{bmatrix}, \mathbb{L}^T = \begin{bmatrix} \mathbf{L}^T \\ \mathbf{0} \\ \mathbf{0} \end{bmatrix} \qquad (13) $$

$$ \Delta \mathbb{U} = \begin{bmatrix} \Delta \mathbf{v} \\ \Delta \mathbf{u} \\ \Delta \mathbf{a} \end{bmatrix}, \quad \bar{\lambda} = \begin{bmatrix} \bar{\lambda} \\ \mathbf{0} \\ \mathbf{0} \end{bmatrix} \qquad (14) $$

et

$$ \mathbb{F} = \begin{bmatrix} \eta \mathbf{f}_{n+1} + \mu \mathbf{f}_n \\ \mathbf{0} \\ \mathbf{0} \end{bmatrix}, \quad \mathbb{U}_n = \begin{bmatrix} \mathbf{v}_n \\ \mathbf{u}_n \\ \mathbf{a}_n \end{bmatrix} \qquad (15) $$

$$ \mathbb{N} = \begin{bmatrix} -(1+\alpha)h\mathbf{K} & -(\eta+\mu)\mathbf{K} & (\eta-\mu)\left[\left(\frac{1}{\gamma}-1\right)\mathbf{M} + (1+\alpha)(\beta-\frac{\gamma}{2})h^2\mathbf{K}\right] \\ h\mathbf{I} & \mathbf{0} & \left(\frac{\gamma}{2}-\beta\right)h^2\mathbf{I} \\ \mathbf{0} & \mathbf{0} & -\frac{1}{\gamma}\mathbf{I} \end{bmatrix} \qquad (16) $$

L'incrément de déplacement généralisé est la somme d'un incrément “libre” (sans conditions de liaison) et d'un supplément dû aux liaisons :

$$ \Delta \mathbb{U} = \Delta \mathbb{U}_{free} + \Delta \mathbb{U}_{link} \qquad (17) $$

On a besoin des formules de récurrence pour passer de l'éat n à l'état n+1. Pour tous les schémas sauf le schéma de Krenck ces formules s'écrivent :

$$ \mathbf{u}_{n+1} = \mathbf{u}_n + h\mathbf{v}_n + \left(\frac{1}{2} - \beta\right) h^2 \mathbf{a}_n + \beta h^2 \mathbf{a}_{n+1} \qquad (18) $$

$$ \mathbf{v}_{n+1} = \mathbf{v}_n + (1 - \gamma) h \mathbf{a}_n + \gamma h \mathbf{a}_{n+1} \qquad (19) $$

ou γ et β sont les paramètres usuels du schéma de Newmark. Pour le schéma de Krenk les formules de récurrence sont :

$$ \mathbf{u}_{n+1} = \mathbf{u}_n + h\mathbf{v}_n + \frac{h}{2}\Delta\mathbf{v} \qquad (20) $$

$$ \mathbf{v}_{n+1} = \mathbf{v}_n + \Delta\mathbf{v} \qquad (21) $$

Le domaine de stabilité du schéma HHT est une fonction connue de α, β et γ ([3]). On n'a pas de démonstration de la stabilité de ce schéma en non linéaire. La stabilité des schémas de SIMO et VERLET est la même que celle de Newmark en linéaire mais elle est aussi établie en non linéaire. Il en est de même pour le schéma de Krenk.
<figure><img src="image_2.png" /><figcaption>Figure 1 – (a) Collage de sous domaines.</figcaption></figure>

## 4 Coupler les sous domaines

Des forces de liaison sont à ajouter entre les sous domaines pour les “coller” comme sur la figure Figure 1.

Chaque sous domaine (ou code) k est intégré avec son propre schéma d’intégration k. Les forces de liaison sont des multiplicateurs de Lagrange. Les équations du mouvement pour le sous domaine k sont :

$$ 
\left\{ 
\begin{array}{l} 
\left( \frac{1}{h\gamma_k} \mathbf{M}_k + h \frac{\beta_k}{\gamma_k} \mathbf{K}_k \right) \Delta \mathbf{v}^k + \mathbf{L}_k^T \bar{\lambda} = \eta_k \mathbf{f}_n^k + \mu_k \mathbf{f}_n^k \\ 
\qquad - (\eta_k + \mu_k) \mathbf{K}_k \mathbf{u}_n^k - h \mathbf{K}_k \mathbf{v}_n^k \\ 
+ \left( (\eta_k - \mu_k) \left( \frac{1}{\gamma_k} - 1 \right) \mathbf{M}_k - \left( \frac{\gamma_k}{2} - \beta_k \right) h^2 \mathbf{K}_k \right) \mathbf{a}_n^k 
\end{array} 
\right\} 
\quad (22) 
$$

Nous allons maintenant présenter le couplage entre deux sous domaines A et B ayant le même pas de temps h puis généraliser au cas ou les pas de temps sont différents.

### 4.1 Même pas de temps

Les pas de temps sont identiques mais les schémas sont différents. Le système à résoudre est alors :

$$ 
\left\{ 
\begin{array}{l} 
\left( \frac{1}{h\gamma_A} \mathbf{M}_A + h \frac{\beta_A}{\gamma_A} (1+\alpha_A) \mathbf{K}_A \right) \Delta \mathbf{v}^A + \mathbf{L}_A^T \bar{\lambda} = \eta_A \mathbf{f}_n^A + \mu_A \mathbf{f}_n^A \\ 
\qquad - (\eta_A + \mu_A) \mathbf{K}_A \mathbf{u}_n^A - h(1+\alpha_A) \mathbf{K}_A \mathbf{v}_n^A \\ 
+ (\eta_A - \mu_A) \left[ \left( \frac{1}{\gamma_A} - 1 \right) \mathbf{M}_A - (1+\alpha_A) \left( \frac{\gamma_A}{2} - \beta_A \right) h^2 \mathbf{K}_A \right] \mathbf{a}_n^A \\ \\ 
\left( \frac{1}{h\gamma_B} \mathbf{M}_B + h \frac{\beta_B}{\gamma_B} (1+\alpha_B) \mathbf{K}_B \right) \Delta \mathbf{v}^B + \mathbf{L}_B^T \bar{\lambda} = \eta_B \mathbf{f}_n^B + \mu_B \mathbf{f}_n^B \\ 
\qquad - (\eta_B + \mu_B) \mathbf{K}_B \mathbf{u}_n^B - h(1+\alpha_B) \mathbf{K}_B \mathbf{v}_n^B \\ 
+ (\eta_B - \mu_B) \left[ \left( \frac{1}{\gamma_B} - 1 \right) \mathbf{M}_B - (1+\alpha_B) \left( \frac{\gamma_B}{2} - \beta_B \right) h^2 \mathbf{K}_B \right] \mathbf{a}_n^B \\ \\ 
L_A \Delta \mathbf{v}^A + L_B \Delta \mathbf{v}^B = 0 
\end{array} 
\right. 
\quad (23) 
$$

L'écriture matricielle est :

$$ 
\begin{bmatrix} \mathbf{K}_A^* & \mathbf{0} & \mathbf{L}_A^T \\ \mathbf{0} & \mathbf{K}_B^* & \mathbf{L}_B^T \\ \mathbf{L}_A & \mathbf{L}_B & \mathbf{0} \end{bmatrix} \begin{bmatrix} \Delta \mathbf{U}^A \\ \Delta \mathbf{U}^B \\ \bar{\lambda} \end{bmatrix} = \begin{bmatrix} \mathbb{F}^A \\ \mathbb{F}^B \\ \mathbf{0} \end{bmatrix} - \begin{bmatrix} \mathbb{N}_A & \mathbf{0} & \mathbf{0} \\ \mathbf{0} & \mathbb{N}_B & \mathbf{0} \\ \mathbf{L}_A & \mathbf{L}_B & \mathbf{0} \end{bmatrix} \begin{bmatrix} \mathbb{U}_n^A \\ \mathbb{U}_n^B \\ \mathbf{0} \end{bmatrix} 
$$
La même méthode que dans la cas d'un seul sous domaine s'applique. On résoud le problème libre pour chaque sous domaine. On calcule les efforts de liaison qui assurent la continuité cinématique, puis on en déduit les inconnues liées dans chaque sous domaine qu'on ajoute aux inconnues libres pour avoir la solution du problème. Les étapes 1 et 3 sont naturellement paralleles.

## 4.2 Pas de temps différents

On supposera ici que les pas de temps des sous domaines A et B sont multiples. On a :

$$h_A = mh_B \qquad (24)$$

Soit j un des sous pas ($j \in [1, m-1]$). Posons :

$$\overline{\Delta \mathbf{x}} = (\mathbf{x}_{j+1} - \mathbf{x}_j), \quad \overline{\bar{\mathbf{x}}} = \int_{jh_B}^{(j+1)h_B} x(\tau) \, d\tau \qquad (25)$$

Supposons maintenant que le multiplicateur de Lagrange à l'échelle fine $\overline{\bar{\lambda}}_B^j$ est constant sur tout les pas de temps $h_B$ et égal à $\bar{\lambda}_B$ pendant le pas de temps $h_A$. On en déduit :

$$\Delta W_{AB}^{\text{interface}} = \left( (\Delta \mathbf{v}^A)^T \mathbf{L}_A^T + \sum_{j=1}^{m} (\Delta \mathbf{v}_j^B)^T \mathbf{L}_B^T \right) \bar{\lambda} \qquad (26)$$

Les conditions d'égalité des vitesses sont :

$$\mathbf{L}_A \Delta \mathbf{v}^A + \mathbf{L}_B \Delta \mathbf{v}^B = 0 \qquad (27)$$

qui s'écrit aussi :

$$\mathbf{L}_A \Delta \mathbf{v}^A + \sum_{j=1}^{m} L_B \Delta \mathbf{v}_j^B = 0 \qquad (28)$$

On en déduit que le travail d'interface est nul. A chaque grand pas de temps l'équation d'équilibre dynamique s'écrit :

$$K_A^* \Delta \mathbf{U}^A + \mathbf{L}_A^T \bar{\lambda} + N_A \mathbf{U}^A = \mathbf{F}_A \qquad (29)$$

Pour chaque petit pas de temps j l'équation d'équilibre est :

$$K_B^* \Delta (\mathbf{U}_j^B) + \mathbf{L}_B^T \bar{\lambda} + N_B \mathbf{U}_j^B = \mathbf{F}_B \qquad (30)$$

les vecteurs généralisés d'inconnues sur le sous domaine B (petits pas de temps) $\Delta \mathbf{U}_j^B$ et $\mathbf{U}_j^B$ résultent tout simplement de la sommation des vecteurs aux pas de temps précédents :

$$\Delta \mathbf{U}_j^B = \sum_{i=1}^{j-1} \Delta \mathbf{U}_i^B \quad \text{and} \quad \mathbf{U}_j^B = \mathbf{U}_0^B + \sum_{i=1}^{j-2} \Delta \mathbf{U}_i^B \qquad (31)$$

Après quelques calculs un peu longs on trouve que l'on doit résoudre le système matriciel suivant :

$$\begin{bmatrix} K_B^* & \mathbf{0} & \mathbf{0} & \mathbf{0} & \mathbf{0} & \mathbf{0} & \mathbf{L}_B^T \\ N_B & K_B^* & \mathbf{0} & \mathbf{0} & \mathbf{0} & \mathbf{0} & \mathbf{L}_B^T \\ N_B & N_B & K_B^* & \mathbf{0} & \mathbf{0} & \mathbf{0} & \mathbf{L}_B^T \\ \vdots & \vdots & \ddots & \ddots & \ddots & \vdots & \vdots \\ N_B & N_B & N_B & N_B & K_B^* & \mathbf{0} & \mathbf{L}_B^T \\ \mathbf{0} & \mathbf{0} & \mathbf{0} & \mathbf{0} & \mathbf{0} & K_A^* & \mathbf{L}_A^T \\ \mathbf{L}_B & \mathbf{L}_B & \mathbf{L}_B & \mathbf{L}_B & \mathbf{L}_B & \mathbf{L}_A & \mathbf{0} \end{bmatrix} \begin{bmatrix} \Delta \mathbf{U}_1^B \\ \Delta \mathbf{U}_2^B \\ \Delta \mathbf{U}_3^B \\ \vdots \\ \Delta \mathbf{U}_m^B \\ \Delta \mathbf{U}_m^A \\ \bar{\lambda}_* \end{bmatrix} = \begin{bmatrix} \mathbf{F}_1^B - N_B \mathbf{U}_0^B \\ \mathbf{F}_2^B - N_B \mathbf{U}_0^B \\ \mathbf{F}_3^B - N_B \mathbf{U}_0^B \\ \vdots \\ \mathbf{F}_m^B - N_B \mathbf{U}_0^B \\ \mathbf{F}_m^A - N_A \mathbf{U}_0^A \\ -\mathbf{L}_A \mathbf{U}_0^A - \mathbf{L}_B \mathbf{U}_0^B \end{bmatrix} \qquad (32)$$
<figure><img src="image_3.png" /><figcaption>Figure 2 – Partition en 20 sous domaines et 5 schémas numériques (dans les cercles)</figcaption></figure>

On voit sur ce système formel que la méthode est simple à implémenter. En effet dans le sous domaine A on fait un pas pour trouver les inconnues libres, et pour le sous domaine B on doit faire m résolutions consécutives. L'incrément de la solution libre sur le sous domaine B est tout simplement la somme des m sous incrément s. On résoud ensuite le problème permettant de calculer les forces liaisons entre sous domaines. On calcule ensuite les corrections aux inconnues pour les deux sous domaines de la même manière que pour les inconnues libres. Cette stratégie est donc complètement parallelle et il n'y a pas besoin de faire communiquer A et B pendant les petits pas de temps. Cette méthode ne génère aucune énergie d'interface même pour des pas de temps différents, et de plus elle est plus simple que la méthode proposée par Parkash. Enfin elle est aussi plus générale puisqu'elle s'applique quelque soit le schéma d'intégration temporelle choisi.

# 5 Exemple numérique

On regarde la propagation d'une onde dans une plaque plane carrée de coté $L = 1$. Le module d'Young vaut $E = 1$, ainsi que sa masse volumique $\rho = 1$. Le coefficient de Poisson's est $\mu = 0.3$. La plaque est chargée par un créneau d'amplitude 1 appliqué au coin supérieur droit de la plaque. La plaque est encastrée en son coin inférieur gauche. Le maillage est fait avec des éléments quadrangles à 4 noeuds. La Figure (2) montre le sous découpage en 20 sous domaines ayant chacun leur schéma et leur pas de temps. Les paramètres des schémas sont donnés dans la Table 2. 5 algoritmes avec 5 échelles de temps sont utilisées ici : $h_i = m_i h_1 \forall i \in \{1..5\}$ le plus petit pas étant $h_1 = 0.016$. Remarquons aussi que des maillages spatiaux sont incompatibles aux interfaces. On impose ici une condition de continuité cinématique pour coller les sous domaines spatialement. Le calcul est fait avec 100 petits pas de temps. Le travail d'interface (Figure 3) est de l'ordre de grandeur de $10^{-15}$. L'énergie mécanique totale étant de l'ordre de 1, on constate la nullité de
l'énergie d'interface.

<table><caption>Tableau 2 – Paramètres utilisés α, γ and β</caption><thead><tr><th>Schéma n°</th><th>α</th><th>γ</th><th>β</th></tr></thead><tbody><tr><td>1</td><td>0.0</td><td>0.5</td><td>0.25</td></tr><tr><td>2</td><td>0.0</td><td>0.5</td><td>0.66</td></tr><tr><td>3</td><td>0.0</td><td>0.6</td><td>0.3025</td></tr><tr><td>4</td><td>-0.2</td><td>0.7</td><td>0.36</td></tr><tr><td>5</td><td>0.0</td><td>0.55</td><td>0.27565</td></tr></tbody></table>

<table><caption>Tableau 3 – Ratios des pas de temps utilisés</caption><thead><tr><th>m<sub>1</sub></th><th>m<sub>2</sub></th><th>m<sub>3</sub></th><th>m<sub>4</sub></th><th>m<sub>5</sub></th></tr></thead><tbody><tr><td>1</td><td>2</td><td>4</td><td>8</td><td>16</td></tr></tbody></table>

# 6 Conclusion

Ce travail donne une méthode simple permettant de construire des assemblages de schéma numériques qui ne dissipent aucune énergie d'interface. Le concept essentielle est la formulation faible de l'intégration temporelle des équations d'équilibre. Ces schémas laissent cependant une certaine liberté aux variables cinématiques à l'intérieur du pas de temps. Cette méthode s'étend simplement au non linéaire. Il reste cependant à vérifier les effets des hypothèses faites pour le collage temporel sur la qualité des interfaces dans l'espace des fréquences : quelles distorsions fréquentielles sont gérées par ce type de collage temporel ?

# Références

[1] Chung J, Hulbert GM. A time integration algorithm for structural dynamics with improved numerical dissipation : the generalized α-method. *Journal of Applied Mechanics* 1993; **60** : 371-375.

[2] Verlet L., in.; *Phys. Rev.* (1967); **165**, 201 : 67-94.

[3] Gerardin M, Rixen D. *Mechanical Vibrations, Theory and Applications to Structural Dynamics* (2nd edn). Wiley : Chichester, 1997.

[4] 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.

[5] Hughes TJR. *The Finite Element Method : Linear Static and Dynamic Finite Element Analysis*. Prentice-Hall : Englewood Cliffs, NJ, 1987.

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

[7] Herry B, Di Valentin L, Combescure A. An approach to the connection between subdomains with non-matching meshes for transient mechanical analysis *International Journal for Numerical Methods in Engineering* 2002; **55** :973-1003.
<figure><img src="image_4.png" /><figcaption>Figure 3 – Interface Work W<sub>global</sub></figcaption></figure>
[8] Faucher V, Combescure A. Local modal reduction in explicit dynamics with domain decomposition. Part 1 : extension to subdomains undergoing finite rigid rotations *International Journal for Numerical Methods in Engineering* 2004; **60** :2531-2560.

[9] Belytschko T, Mullen R. Mesh partitions of explicit-implicit time integration In *Proceedings of U.S. - Germany Symposium on Formulations and Computational Algorithms in Finite Element Analysis* 1976

[10] Belytschko T, Mullen R. Stability of explicit-implicit mesh partitions in time integration *International Journal for Numerical Methods in Engineering* 1978; **12** :1575-1586.

[11] Ben Dhia, H. Multiscale mechanical problems : the Arlequin method *Comptes Rendus de l'Académie des Sciences - Series IIB - Mechanics-Physics-Astronomy* 1998; **12** :899-904.

[12] Bourel B. Calcul multi-domaines et approches multi-échelles pour la simulation de crashs automobiles *Thèse de doctorat*. Institut National des Sciences Appliquées de Lyon, 2006.

[13] Bourel B, Combescure A, Valentin LD. Handling contact in multi-domain simulation of automobile crashes *Finite Elements in Analysis and Design* 2006; **42** :766-779.

[14] Cavin P, Gravouil AA, Combescure A. Automatic energy conserving space-time refinement for linear dynamic structural problems. *International Journal for Numerical Methods in Engineering* 2005; **50** :304-321.

[15] Combescure A, Gravouil, A. A numerical scheme to couple subdomains with different time-steps for predominantly linear transient analysis *Computer methods in applied mechanics and engineering* 2002; **191** :1129-1157.

[16] Daniel, WJT. A study of the stability of subcycling algorithms in structural dynamics. *Computer methods in applied mechanics and engineering* 1998; **156** :1-13.

[17] Dodds JRH, Lopez LA. Substructuring in linear and nonlinear analysis. *International Journal Numerical Methods in Engineering* 1980; **15** :583-597.

[18] Farhat C, Roux F-X. A method of finite element tearing and interconnecting and its parallel solution algorithm. *International Journal Numerical Methods in Engineering* 1991; **32** :1205-1227.

[19] Farhat C, Crivelli L, Geradin M. On the spectral stability of time integration algorithms for a class of constrained dynamics problems. In *Collection of Technical Papers AIAA/ASME Structures, Structural Dynamics and Materials Conference*. AIAA : Washington DC, U.S.A., 1993 ; 80-97.

[20] Farhat C, Crivelli L, Roux F-X. Transient FETI methodology for large-scale parallel implicit computations in structural mechanics. *International Journal Numerical Methods in Engineering* 1994; **37** :1945-1975.

[21] Faucher V, Combescure A. Local modal reduction in explicit dynamics with domain decomposition. Part 1 : extension to subdomains undergoing finite rigid rotations *International Journal for Numerical Methods in Engineering* 2004; **60** :2531-2560.

[22] Gérardin M, Rixen, D. *Théorie des vibrations - Application à la dynamique des structures (2ème édition)*. Paris : Masson Editions, 1996.

[23] Glowinski R, Le Tallec, P. Augmented Lagrangian interpretation of the non-overlapping Schwarz alternating method In *Domain decomposition method* SIAM : Philadelphia, U.S.A., 1990 ; 224-231.

[24] A Gravouil, A Combescure A Numerical Scheme to Couple Subdomains with Different Time Steps for Predominantly Linear Transient Analysis. *Computer Methods In Applied Mechanics And Engineering, 2002.* **1191-1112** :1129-1157.
[25] Herry B. Développement dune approche multiéchelle parallèle pour la simulation de crash automobible non linéaire des structures. *Thèse de doctorat*. Ecole Normale Supérieure de Cachan, 2002.

[26] Herry B, Di Valentin L, Combescure A. An approach to the connection between subdomains with non-matching meshes for transient mechanical analysis *International Journal for Numerical Methods in Engineering* 2002; **55**:973-1003.

[27] Hughes TJR, Liu WK. Implicit-Explicit finite elements in transient analysis : stability theory. *ASME Journal of Applied Mechanics* 1978; **45**:371-374.

[28] Hughes TJR, Liu WK. Implicit-Explicit finite elements in transient analysis : implementation and numerical examples. *ASME Journal of Applied Mechanics* 1978; **45**:375-378.

[29] Joly P, Rodríguez J. An Error Analysis of Conservative Space-Time Mesh Refinement Methods for the 1D Wave Equation *SIAM J. Numer. Anal.* 2005; **43**:825-859.

[30] Krenk S. Energy conservation and high-frequency damping in numerical time integration. In *Computational Methods in Structural Dynamics and Earthquake Engineering* 2007; Resthymno, Crete, Greece, 2007; 12.

[31] Magoules F, Roux FX. Lagrangian formulation of domain decomposition methods : A unified theory. *Applied Mathematical Modelling* 2006; **30**:593-615.

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

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

[34] Abaqus User's Manual *Theory Manual, Version 5.4., p. 2.3.1-3, 1994.*

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