# Étude d'un couplage discret/continu en dynamique explicite

Sylvain Gavoille, Christian Rey, Arnaud Delaplace, Christian Mariotti

► To cite this version:

Sylvain Gavoille, Christian Rey, Arnaud Delaplace, Christian Mariotti. Étude d'un couplage discret/continu en dynamique explicite. 9e Colloque national en calcul des structures, CSMA, May 2009, Giens, France. hal-01421807

HAL Id: hal-01421807

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

Submitted on 22 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
# Etude d'un couplage discret/continu en dynamique explicite

S. Gavoille¹, C. Rey¹, A. Delaplace ¹ et C. Mariotti ²

LMT - ENS Cachan (ENS Cachan/CNRS/UPMC/PRES UniverSud Paris),
61 avenue du Président Wilson, 94235 Cachan FRANCE.
{gavoille, rey, delaplace} @lmt.ens-cachan.fr

² LSEG/ CEA DAM Ile de France,
BP 12 91680 Bruyères-les-Châtel, FRANCE.
christian.mariotti@cea.fr

**Résumé** — Ce papier décrit une stratégie pour coupler une méthode éléments discrets avec une méthode éléments finis. Une formulation basée sur une partition de l’unité en énergie des champs discrets et continus sur un sous-domaine de raccord est utilisée. Cette approche a pour originalité de limiter le raccord des modèles au bord du sous-domaine de raccord, réduisant le coût de calcul sans nuire à la qualité de la résolution. La stratégie développée est illustrée sur un problème de propagation d’onde dans une structure élancée.

**Mots clés** — couplage de modèles, éléments discrets, éléments spectraux.

# 1 Introduction

Le couplage de modèles offre de nouvelles perspectives en terme d’analyse multi-échelle des matériaux et des structures en statique comme en dynamique. La prise en compte de toute la complexité et la richesse des modèles aux toutes petites échelles est en effet, aujourd’hui encore, inaccessible à grande échelle. Le couplage de modèles nous offre ainsi la possibilité de prendre en compte de tels comportements au moins localement.

Nous considérons ici le couplage d'un modèle de type éléments discrets et d'un modèle continu classique. Un tel couplage nous offre notamment la possibilité de modéliser la propagation d'ondes sismiques à travers un milieu de grande taille constitué de géomatériaux pouvant localement avoir un comportement non linéaire dans une zone d'intérêt. Cette zone est modélisée par la méthode des éléments discrets (DEM) qui permet de prendre en compte avec précision les hétérogénéités matériaux, ainsi que les discontinuités engendrées par la rupture [3, 4]. Quant au reste du domaine, on utilise un modèle continu classique discrétisé par une méthode des éléments spectraux (SEM), méthode élément fini d'ordre élevé adaptée à la modélisation de la propagation d'ondes dans un cadre explicite [1, 2].

Le couplage proposé repose sur une partition des énergies sur une zone de recouvrement. Tou-
tefois, plutôt qu’un couplage volumique [5, 6], nous proposons un couplage surfacique limité au
bord de la zone de recouvrement où la continuité des champs est assurée au moyen de multipli-
cateurs de Lagrange. Une telle approche permet de réduire significativement la taille du problème
interface et le coût de calcul des opérateurs de couplage, pour une réduction limitée de la qualité du
couplage. Les performances de l'approche sont illustrées sur un exemple de propagation d'ondes élastiques.

## 2 Description des modèles

Le problème que l'on souhaite résoudre est discrétisé d'une part avec des éléments spectraux recouvrant le domaine $\Omega^{(S)}$ et d'autre part avec des éléments discrets recouvrant le domaine $\Omega^{(D)}$. Le couplage est réalisé à travers une zone de recouvrement entre ces deux domaines $\Omega^{(J)} = \Omega^{(S)} \cap \Omega^{(D)}$.

### 2.1 Le modèle discret

Le modèle discret considère le mouvement de chaque particule obtenue à partir d'un pavage de Voronoi sur l'ensemble du domaine $\Omega^{(D)}$. Chaque particule a une forme polygonale et est reliée à ces voisines par des liens qui traduisent la cohésion du milieu. La figure 1 illustre ce modèle. Dans le cas où les efforts d'interaction au niveau du lien dépasse un seuil critique, dépendant du matériau étudié, la suppression du lien modélise la rupture.

<figure><img src="image_2.png" /><figcaption>Figure 1 – Représentation du modèle discret</figcaption></figure>

Dans ce système constitué de particules, l'énergie totale, i.e. la somme de l'énergie cinétique et de l'énergie potentielle des particules, est constante en fonction du temps et est identifiée à l'Hamiltonien $H^D$ donné par :

$$H_{\Omega^{(D)}}^{D} = T_{\Omega^{(D)}}^{D} + V_{\Omega^{(D)}}^{D} \qquad (1)$$

où $T_{\Omega^{(D)}}^D$ et $V_{\Omega^{(D)}}^D$ correspondent respectivement à l'énergie cinétique et à l'énergie potentielle du modèle discret dans le domaine $\Omega^{(D)}$. Ces énergies ont pour expression :

$$\left\{ \begin{aligned} T_{\Omega^{(D)}}^{D} &= \frac{1}{2}[\dot{\mathbf{u}}^{D}]^{T}\mathbf{M}_{\Omega^{(D)}}^{D}[\dot{\mathbf{u}}^{D}] + \frac{1}{2}[\dot{\theta}^{D}]^{T}\mathbf{I}_{\Omega^{(D)}}^{D}[\dot{\theta}^{D}] \\ V_{\Omega^{(D)}}^{D} &= \frac{1}{2}[\mathbf{u}^{D}]^{T}\mathbf{K}_{\Omega^{(D)}}^{D}[\mathbf{u}^{D}] - [\mathbf{f}_{\Omega^{(D)}}^{D}]^{T}[\mathbf{u}^{D}] + \frac{1}{2}[\mathbf{\theta}^{D}]^{T}\mathbf{G}_{\Omega^{(D)}}^{D}[\mathbf{\theta}^{D}] - [\mathbf{m}_{\Omega^{(D)}}^{D}]^{T}[\mathbf{\theta}^{D}] \end{aligned} \right. \qquad (2)$$

où :

- $\mathbf{u}^D$ et $\mathbf{\theta}^D$ sont les vecteurs rassemblant respectivement les déplacements et rotations de chaque particule.

- $\mathbf{M}_{\Omega^{(D)}}^D$ et $\mathbf{I}_{\Omega^{(D)}}^D$ sont respectivement les matrices assemblant les masses et inerties de chaque particule. Notons que ces deux matrices sont diagonales, ce qui permet de mettre en place naturellement un schéma de résolution explicite.

- $\mathbf{K}_{\Omega^{(D)}}^D$ et $\mathbf{G}_{\Omega^{(D)}}^D$ sont respectivement les matrices assemblant les raideurs en déplacement et en rotation traduisant la cohésion du milieu.

- $\mathbf{f}_{\Omega^{(D)}}^D$ et $\mathbf{m}_{\Omega^{(D)}}^D$ sont les vecteurs traduisant respectivement les efforts et moments extérieurs.
## 2.2 Le modèle continu

Le modèle continu est constitué des équations de base de l'élastodynamique appliquées au domaine d'étude $\Omega^{(S)}$. Ce domaine répond à un équilibre dynamique et à un comportement matériau caractérisé par la loi de Hooke. L'expression de l'Hamiltonien du modèle continu discrétisé avec des éléments spectraux sur le domaine $\Omega^{(S)}$ est la suivante :

$$H_{\Omega^{(S)}}^{S} = T_{\Omega^{(S)}}^{S} + V_{\Omega^{(S)}}^{S} \qquad (3)$$

$T_{\Omega^{(S)}}^{S}$ et $V_{\Omega^{(S)}}^{S}$ correspondent respectivement à l'énergie cinétique et de l'énergie potentielle associé au domaine $\Omega^{(S)}$. Leur expression complète est la suivante :

$$\left\{ \begin{aligned} T_{\Omega^{(S)}}^{S} &= \frac{1}{2} [\dot{\mathbf{u}}^{S}]^{T} \mathbf{M}_{\Omega^{(S)}}^{S} [\dot{\mathbf{u}}^{S}] \\ V_{\Omega^{(S)}}^{S} &= \frac{1}{2} [\mathbf{u}^{S}]^{T} \mathbf{K}_{\Omega^{(S)}}^{S} [\mathbf{u}^{S}] - [\mathbf{f}_{\Omega^{(S)}}^{S}]^{T} [\mathbf{u}^{S}] \end{aligned} \right. \qquad (4)$$

où $\mathbf{u}^S$ est le champ de déplacement recherché. $\mathbf{M}_{\Omega^{(S)}}^S, \mathbf{K}_{\Omega^{(S)}}^S$ et $\mathbf{f}_{\Omega^{(S)}}^S$ correspondent respectivement à la matrice de masse, la matrice de raideur et au vecteur des efforts extérieurs obtenus à partir d'une discrétisation par des éléments spectraux du modèle continu sur le domaine $\Omega^{(S)}$.

# 3 Méthode de couplage

La méthode de couplage proposée reprend certaines considérations de la méthode Arlequin : la décomposition de domaine est réalisée avec recouvrement. A chaque domaine est associé respectivement un modèle discret et continu. Afin d'être conservatif en énergie, on affecte dans la zone de recouvrement $\Omega^{(J)}$ le paramètre $\alpha$ ($\in[0, 1]$) à l'énergie associée au modèle continu et $(1 - \alpha)$ à l'énergie associée au modèle discret.

L'originalité du couplage proposé consiste à limiter le raccord des champs sur le bord de la zone de recouvrement. Cela permet de résoudre à moindre coût le problème de couplage vis-à-vis des méthodes existantes [5, 6]. Pour cela, on introduit un multiplicateur de Lagrange $\lambda$, représentant un effort d'interface, permettant d'imposer la condition de continuité cinématique au niveau de l'interface. Nous assurons seulement la continuité sur les champs de déplacement, sans considérer la rotation des particules.

Dans cette approche, l'Hamiltonien du problème permettant de lier les deux modèles s'écrit :

$$H = H_{\Omega_c^{(S)}}^S + H_{\Omega_c^{(D)}}^D + \alpha H_{\Omega_c^{(J)}}^S + (1 - \alpha) H_{\Omega_c^{(J)}}^D + \lambda^T (\mathbf{C}^{(S)}\mathbf{u}^S - \mathbf{C}^{(D)}\mathbf{u}^D) \qquad (5)$$

avec $\Omega_c^{(S)} = \Omega^{(S)} \setminus \Omega^{(J)}$ et $\Omega_c^{(D)} = \Omega^{(D)} \setminus \Omega^{(J)}$, et $\mathbf{C}^{(S)}$ et $\mathbf{C}^{(D)}$ les opérateurs de couplage permettant le passage respectivement des éléments spectraux et des éléments discrets au bord de la zone de recouvrement.

Le développement canonique de l'Hamiltonien $H$ donne les équations suivantes :

$$\left\{ \begin{aligned} \mathbf{M}^{(S)}\ddot{\mathbf{u}}^{(S)} + \mathbf{K}^{(S)}\mathbf{u}^{(S)} &= \mathbf{f}^{(S)} - \mathbf{C}^{(S)T}\lambda \\ \mathbf{M}^{(D)}\ddot{\mathbf{u}}^{(D)} + \mathbf{I}^{(D)}\ddot{\theta}^{(D)} + \mathbf{K}^{(D)}\mathbf{u}^{(D)} + \mathbf{G}^{(D)}\theta^{(D)} &= \mathbf{f}^{(D)} + \mathbf{m}^{(D)} + \mathbf{C}^{(D)T}\lambda \\ \mathbf{C}^{(S)}\ddot{\mathbf{u}}^{(S)} - \mathbf{C}^{(D)}\ddot{\mathbf{u}}^{(D)} &= \mathbf{0} \end{aligned} \right. $$

où $\bullet^{(S)}$ désigne les termes correspondant au modèle continu discrétisé avec des éléments spectraux et $\bullet^{(D)}$ ceux correspondant au modèle élément discret. La partition en $\alpha$ sur la zone $\Omega^{(J)}$ nous amène à construire les opérateurs de cette façon :
$$ M^{(S)} = M_{\Omega_c^{(S)}}^{S} + \alpha M_{\Omega_c^{(J)}}^{S} \quad ; \quad K^{(S)} = K_{\Omega_c^{(S)}}^{S} + \alpha K_{\Omega_c^{(J)}}^{S} $$

$$ f^{(S)} = f_{\Omega_c^{(S)}}^{S} + \alpha f_{\Omega_c^{(J)}}^{S} $$

$$ M^{(D)} = M_{\Omega_c^{(D)}}^{D} + (1 - \alpha) M_{\Omega_c^{(J)}}^{D} \quad ; \quad K^{(D)} = K_{\Omega_c^{(D)}}^{D} + (1 - \alpha) K_{\Omega_c^{(J)}}^{D} \quad (6) $$

$$ I^{(D)} = I_{\Omega_c^{(D)}}^{D} + (1 - \alpha) I_{\Omega_c^{(J)}}^{D} \quad ; \quad G^{(D)} = G_{\Omega_c^{(D)}}^{D} + (1 - \alpha) G_{\Omega_c^{(J)}}^{D} $$

$$ f^{(D)} = f_{\Omega_c^{(D)}}^{D} + (1 - \alpha) f_{\Omega_c^{(J)}}^{D} \quad ; \quad m^{(D)} = m_{\Omega_c^{(D)}}^{D} + (1 - \alpha) m_{\Omega_c^{(J)}}^{D} $$

# 4 Résultats numériques

## 4.1 Description du cas test

On s'intéresse à la réponse élastodynamique d'un milieu ramené à un problème à deux dimensions. Le matériau étudié correspond à une roche répondant à la loi de Hooke et ayant pour caractéristiques un module d'Young E de 18.8 GPa, un coefficient de Poisson v de 0.25 et une densité massique ρ de 2200 kg/m³. La géométrie étudiée a une forme rectangulaire de 1000 m de longueur et de 50 m de largeur (figure 2).

<figure><img src="image_3.png" /><figcaption>Figure 2 – Problème de référence</figcaption></figure>

Ce milieu est soumis à une pression surfacique uniforme notée $F_d(t)$ sur le bord gauche. La pression $F_d(t)$ est un créneau de période 0.1 s. Le spectre fréquenciel de cette fonction est complet, ce qui permet d'obtenir une analyse complète en fréquence du couplage. Par la suite, on note $u_{\text{ref}}$ la solution obtenue en résolvant le problème sans couplage par une méthode continue.

## 4.2 Construction de la fonctionnelle d'écart de modèle

Pour qualifier l'influence du couplage, nous proposons de définir une fonctionnelle d'écart de modèle, avec pour référence la solution continue $\ddot{u}_{\text{ref}}$. Pour une particule $i$, l'écart de modèle local $e_i$ est donnée par :

$$ e_i(t) = \frac{\| \ddot{u}_i^D(t) - \ddot{u}_{\text{ref}}(\mathbf{x}_i, t) \|}{\text{Max}_{\mathbf{x}_i} \| \ddot{u}_{\text{ref}}(\mathbf{x}_i, t) \|} \quad (7) $$

où $\mathbf{u}_i^D$ est le vecteur déplacement associé à la particule $i$ et $\mathbf{x}_i$ la position de celle-ci. L'écart de modèle global sur le domaine $\Omega$ est défini par :

$$ \Psi = \frac{\sum_{i \in \Omega} e_i(t)}{N} \quad (8) $$

où $O_i$ est le centre de masse de la particule $i$ et $N$ est le nombre de particules contenues dans $\Omega$.
## 4.3 Analyse numérique du problème couplé

On s'intéresse désormais à l'analyse numérique du problème couplé. Pour ce faire, on réalise une décomposition de domaine du problème présenté précédemment. La partie gauche du domaine, occupant les 250 premiers mètres de la structure, est discrétisée avec des éléments spectraux. Le pas de discrétisation spatiale dans le domaine continu est de 25 m. Les éléments spectraux choisis sont d'ordre 5. Quant à la discrétisation éléments discrets, elle occupe les 775 derniers mètres de la structure. De cette façon, la zone de recouvrement a une dimension de 25 m par 50 m. On note *h* le pas de discrétisation dans le domaine discret qui fera l'objet d'une étude paramétrique. Cela nous donne dans le cas d'une discrétisation spatiale moyenne dans le domaine discret de 6.25 m le maillage illustré à la figure 3.

<figure><img src="image_4.png" /><figcaption>Figure 3 – Problème de couplage discrétisé</figcaption></figure>

La solution obtenue avec la décomposition de domaine continu/discret est comparée à la solution obtenue par un domaine purement discret. On utilise l'écart de modèle $\Psi$ sur le domaine $\Omega_C^{(D)}$, à l'instant $t = 150$ ms. A cet instant, les ondes se sont suffisamment propagées dans le domaine discret pour que l'analyse soit pertinente. L'évolution de l'écart pour différents pas de discrétisation du milieu discret est donnée sur la figure 4. L'écart diminue sans surprise avec le raffinement du maillage, la solution discrète tendant alors vers la solution continue prise ici comme référence. Pour le problème couplé, le paramètre $\alpha$ est fixé à 0,5. Le résultat obtenu montre la bonne qualité du couplage proposé, avec un écart de modèle similaire au problème discret seul.

<figure><img src="image_5.png" /><figcaption>Figure 4 – Etude comparative du problème avec et sans décomposition de domaine.</figcaption></figure>

# 5 Conclusion

Le couplage avec recouvrement permet la résolution de problème avec décomposition de domaine utilisant différents modèles et maillages dans les différents domaines. Nous proposons une stratégie de couplage, basée classiquement sur une partition de l'unité dans la zone de couplage,
mais limitant la continuité du champ cinématique au bord du domaine de recouvrement. Le coût numérique du couplage s'en trouve réduit, sans nuire à la qualité du couplage, destinée ici à la représentation de la propagation d'ondes sismiques dans un sol. La stratégie proposée est illustrée sur un problème d'élastodynamique.

Références

[1] A.T. Patera. *A spectral element method for fluid dynamics, laminar flow in a channel expansion*. J. Comput. Phys., vol. 54, pages 468-488, 1984.

[2] D. Komatitsch, J.P. Vilotte, R. Vai, J.M. Castillo-Covarrubias & F.J. Sanchez-Sesma. *The spectral element method for elastic wave equations application to 2-D and 3-D seismic problems*. Int. J. Numer. Methods Eng., vol. 25, pages 1139-1164, 1998.

[3] P.A. Cundall. *A computer model for simulating progressive large scale movements of blocky rock systems*. Proceedings of the symposium of the international society of rock mechanics, vol. 1, pages 132-150, 1971.

[4] A. Delaplace & R. Desmorat. *Discrete 3D model as complimentary numerical testing for anisotropic damage*. International Journal of Fracture, vol. 148, pages 115-128, 2007.

[5] H. Ben Dhia. *Multiscale mechanical problems : the Arlequin method*. CRAS Paris Série IIb, vol. 326, pages 899-904, 1998.

[6] S.P. Xiao & T. Belytschko. *A bridging domain decomposition for coupling continua with molecular dynamics*. Comput. Methods Appl. Mech. Engrg., vol. 193, pages 1645-1669, 2004.