# Étude énergétique du couplage continu/discret en dynamique explicite

Pascal Aubertin, Julien Réthoré, René de Borst

► To cite this version:

Pascal Aubertin, Julien Réthoré, René de Borst. Étude énergétique du couplage continu/discret en dynamique explicite. 9e Colloque national en calcul des structures, CSMA, May 2009, Giens, France. hal-01421800

HAL Id: hal-01421800

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

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 énergétique du couplage continu / discret en dynamique explicite

P. Aubertin¹, J. Réthoré¹, R. de Borst²

¹ LaMCoS (INSA Lyon/CNRS/Université de Lyon)
20 Avenue Albert Einstein, 69 621 Villeurbanne
{pascal.aubertin,julien.rethore}@insa-lyon.fr

² TU Eindhoven
Eindhoven, The Netherlands
R.d.Borst@tue.nl

**Résumé** — L'analyse multi-échelle et/ou multi-physique est une thématique de recherche en pleine expansion. Il s'agit de coupler, des échelles les plus fines aux échelles les plus grossières, les modèles les plus pertinants pour chaque échelle. On propose ici une étude du couplage dynamique moléculaire / élément finis avec application à la propagation dynamique des fissures.

**Mots clés** — couplage énergétique, dynamique explicite, fissure.

## 1 Introduction

L'analyse multi-échelle et/ou multi-physique est une thématique de recherche en pleine expansion. Il s'agit de coupler, des échelles les plus fines aux échelles les plus grossières, les modèles les plus pertinants pour chaque échelle. Pour cela, il est nécessaire de pouvoir coupler des modèles dédiés à des échelles différentes. On s'intéressera ici au couplage entre un modèle de dynamique moléculaire et une description éléments finis du matériau. On s'inspirera des travaux de Ben Dhia et Rateau [4] pour coupler les deux modèles avec une zone de recouvrement. Néanmoins, les deux modèles utilisés ici sont d'une nature complétement différente : l'un est un modèle continu discrétisé, l'autre un modèle intrinsèquement discret. Cela nécessite, comme on peut le conclure à la lecture de [6], quelques adaptations.

## 2 Formulation

Le domaine étudié $\Omega$ est divisé en $\Omega_M$ et $\Omega_m$ qui feront respectivement l'objet d'une discrétisation par éléments finis et par dymanique moléculaire. Dans $\Omega_M$, la forme forte des équations de conservation de la quantité de mouvement donnent en l'absence de force de volume :

$$ \rho \ddot{u} = \text{div } \sigma \qquad (1) $$

où $u$ est le déplacement continu qui sera discrétisé par éléments finis, $\rho$ la masse volumique et $\sigma$ le tenseur des contraintes de Cauchy. Pour le domaine atomique $\Omega_m$, on a :

$$ m_i \ddot{d}_i = f_i \qquad (2) $$
où $m_i$ désigne la masse de l'atome $i$, $\mathbf{d}_i$ son déplacement et $\mathbf{f}_i$ la somme des forces exercées par les atomes voisins. Il est à noter que cette force dérive d'un potentiel (Lennard-Jones par exemple).

On peut associer à ces deux problèmes leur forme faible. Pour le modèle continu, dans le cadre de l'élasticité linéaire en petites perturbations et en l'absence de force volumique, pour toute fonction test $\mathbf{v}^*$ cinématiquement admissible à 0, on a :

$$a_M(\mathbf{u}, \mathbf{v}^*) = l_M(\mathbf{v}^*) \qquad (3)$$

avec

$$a_M(\mathbf{u}, \mathbf{v}^*) = \int_{\Omega_M} \rho \cdot \ddot{\mathbf{u}} \cdot \mathbf{v}^* d\Omega + \int_{\Omega_M} \boldsymbol{\varepsilon}(\mathbf{u}) : \mathbf{K} : \boldsymbol{\varepsilon}(\mathbf{v}^*) d\Omega \qquad (4)$$

et

$$l_M(\mathbf{v}^*) = \int_{\partial_2 \Omega} \mathbf{F}_d \cdot \mathbf{v}^* dS \qquad (5)$$

où $\mathbf{K}$ est l'opérateur de Hooke et $\mathbf{F}_d$ les efforts imposés sur la partie du bord $\partial_2 \Omega$.

Pour le domaine atomique l'écriture d'une forme faible est purement formelle. On peut néan-
moins écrire :

$$a_m(\mathbf{d}, \mathbf{w}^*) = 0 \qquad (6)$$

avec $\mathbf{w}^*$ les fonctions tests et

$$a_m(\mathbf{d}, \mathbf{w}^*) = \sum_{i=1}^{N_a} m_i \ddot{\mathbf{d}}_i \cdot \mathbf{w}_i^* + \sum_{i=1}^{N_a} \nabla_i \mathcal{U}(\mathbf{d}) \cdot \mathbf{w}_i^*$$

où $\mathcal{U}$ est le potentiel dont dérivent les forces inter-atomiques. Trois points restent alors à préciser :

- écrire des problèmes discrétisés à partir de leur forme faible

- former une partition des énergies entre les deux modèles sur la zone de recouvrement $\Omega_c$

- imposer des conditions de compatibilité cinématique entre les deux modèles dans $\Omega_c$

La troisième étape est la plus directe et nous choisissons d'imposer la compatibilité des vitesses grâce à des multiplicateurs de Lagrange portés par des atomes du modèle atomique dans la zone de couplage. Pour cela, on utilise un opérateur $\Pi$ qui projette les vitesses $\dot{\mathbf{u}}$ et $\dot{\mathbf{d}}$ sur un sous-espace de $\Omega_c$. Si $\mu^*$ sont les fonctions tests des multiplicateurs de Lagrange, les conditions de couplage s'écrivent :

$$\langle \mu^*, \Pi \dot{\mathbf{u}} - \Pi \dot{\mathbf{d}} \rangle = 0 \qquad (7)$$

La deuxième s'effectue au moyen d'une partition de l'unité $\alpha(\mathbf{x}) + \beta(\mathbf{x}) = 1$ qui vient pondérer les énergies de chacun des modèles.

$$\left\{ \begin{array}{ll} \alpha(\mathbf{x}) = 1 & \text{for } \mathbf{x} \in \Omega_M \setminus \Omega_c \\ \beta(\mathbf{x}) = 1 & \text{for } \mathbf{x} \in \Omega_m \setminus \Omega_c \end{array} \right. \qquad (8)$$

Le premier point est obtenu en adoptant un schéma de discrétisation éléments finis (éventuellement étendus pour prendre en compte la présence de fissures) pour le modèle continu. La discrétisation du modèle atomique est immédiate. Pour les deux premières points, nous allons montrer que ces deux opérations ne commutent pas et que l'ordre dans lequel on choisit de les effectuer à une importance considérable sur le fonctionnement du couplage. Nous sommes ici en effet dans le cas très spécifique du couplage de deux modèles de nature intrinsèquement différentes.
## 2.1 Pondération + discrétisation

C'est la voie usuelle, et celle suivie dans [6]. Après avoir introduit les pondérations en énergie $\alpha$ et $\beta$ dans les formes linéaire et bilinéaire associées à chacun des problèmes dans leur écriture faible, la discrétisation est adoptée. Il en résulte le système discret suivant :

$$ egin{cases} 
\mathbf{M}_{\alpha} \ddot{\mathbf{U}} + \mathbf{F}_{\alpha} = \mathbf{F}^{\text{ext}} - \mathbf{C}_{M} \Lambda \\ 
\mathbf{m}_{\beta} \ddot{\mathbf{d}} = \mathbf{f}_{\beta} + \mathbf{C}_{m} \Lambda \\ 
\mathbf{C}_{M}^{T} \dot{\mathbf{U}} = \mathbf{C}_{m}^{T} \dot{\mathbf{d}} 
\end{cases} 
\quad (9) $$

avec $(\mathbf{U}, \mathbf{d}, \Lambda)$ l'ensemble des inconnues du problème couplé, $\mathbf{M}_{\alpha}, \mathbf{m}_{\beta}$ les matrice de masse pondérées, $\mathbf{F}_{\alpha}, \mathbf{f}_{\beta}$ les vecteur des forces internes pondérés, $\mathbf{C}_{M}, \mathbf{C}_{m}$ les opérateurs de couplage discrétisés et $\mathbf{F}^{\text{ext}}$ le vecteur des forces externes.

## 2.2 Discrétisation + pondération

On propose ici de pondérer non pas les opérateurs continus mais les opérateurs discretisés. Pour cela, on définit des matrices de pondération $\alpha$ et $\beta$ agissant sur les énergies cinétique et interne discrétisées. Le système obtenu est le suivant :

$$ egin{cases} 
\alpha (\mathbf{M}\ddot{\mathbf{U}}) + \alpha \mathbf{F} = \mathbf{F}^{\text{ext}} - \mathbf{C}_{M}\Lambda \\ 
\beta (\mathbf{m}\ddot{\mathbf{d}}) = \beta \mathbf{f} + \mathbf{C}_{m}\Lambda \\ 
\mathbf{C}_{M}^{T}\dot{\mathbf{U}} = \mathbf{C}_{m}^{T}\dot{\mathbf{d}} 
\end{cases} 
\quad (10) $$

## 2.3 Conservation de l'énergie discrétisée

Une fois les systèmes discrétisés en espace, un schéma d'intégration en temps est adopté : on choisit le schéma explicite de Newmark de la différence centrée. On peut en étudier la conservation de l'énergie en utilisant la méthode proposée dans [5]. On calcule la valeur moyenne (notation $\ll$) des équations du système (9) (ou (10)) entre deux instants $t_n$ et $t_{n+1}$. On multiplie ensuite par l'incrément de déplacement $[\mathbf{U}]$ pour le modèle continu discrétisé et $[\mathbf{d}]$ pour le modèle moléculaire pour obtenir finalement :

$$ \left[ \frac{1}{2} \dot{\mathbf{U}}^T \tilde{\mathbf{M}}_\alpha \dot{\mathbf{U}} - \frac{(\Delta t)^2}{8} \ddot{\mathbf{U}}^T \tilde{\mathbf{M}}_\alpha \ddot{\mathbf{U}} + \frac{1}{2} \mathbf{U}^T \mathbf{K}_\alpha \mathbf{U} \right] = -[\mathbf{U}]^T \langle \mathbf{C}_M \Lambda \rangle \quad (11) $$

$$ \left[ \frac{1}{2} \dot{\mathbf{d}}^T \mathbf{m}_\alpha \dot{\mathbf{d}} - \frac{(\Delta t)^2}{8} \ddot{\mathbf{d}}^T \mathbf{m}_\alpha \ddot{\mathbf{d}} \right] = [\mathbf{d}]^T \langle \mathbf{f}_\beta \rangle + [\mathbf{d}]^T \langle \mathbf{C}_m \Lambda \rangle \quad (12) $$

où l'on identifie les quantités suivantes : (i) l'énergie cinétique discrétisée $E_K^M = \frac{1}{2} \dot{\mathbf{U}}^T \tilde{\mathbf{M}}_\alpha \dot{\mathbf{U}}$ du modèle continu, (ii) l'énergie interne discrétisée $E_P^M = \frac{1}{2} \mathbf{U}^T \mathbf{K}_\alpha \mathbf{U}$ du modèle continu, (iii) l'énergie cinétique $E_K^m = \frac{1}{2} \dot{\mathbf{d}}^T \mathbf{m}_\alpha \dot{\mathbf{d}}$ dans les atomes, (iv) les termes d'accélération $\frac{(\Delta t)^2}{8} \ddot{\mathbf{U}}^T \tilde{\mathbf{M}}_\alpha \ddot{\mathbf{U}}$ et $\frac{(\Delta t)^2}{8} \ddot{\mathbf{d}}^T \mathbf{m}_\alpha \ddot{\mathbf{d}}$ dus au schéma de la différence centrée, (v) le travail des multiplicateurs de L agrange dans le déplacement du modèle continu et dans le déplacement des atomes, $[\mathbf{U}]^T \langle \mathbf{C}_M \Lambda \rangle$ et $[\mathbf{d}]^T \langle \mathbf{C}_m \Lambda \rangle$, respectivement, et (vi) le travail des forces inter-atomiques $[\mathbf{d}]^T \langle \mathbf{f}_\beta \rangle$.

Grâce à ce bilan d'énergie de la formulation discrétisée, nous allons pouvoir évaluer les performances des deux stratégies de couplage d'un point de vue de la conservation de l'énergie.
<figure><img src="image_2.png" /></figure>

<figure><img src="image_3.png" /></figure>

Figure 1 – Propagation d'une onde 1-D dans un barreau en utilisant la stratégie usuelle pondération + discrétisation.

<figure><img src="image_4.png" /></figure>

<figure><img src="image_5.png" /></figure>

Figure 2 – Propagation d'une onde 1-D dans une barreau en utilisant le stratégie proposée discrétisation + pondération.

# 3 Applications

## 3.1 Barreau 1D

Dans le premier exemple on étudie la propagation d'une onde dans un barreau unidimensionel. La taille totale du domaine est 59.142528 × 10⁻⁹ m et 100 atomes sont utilisés dans le domaine moléculaire. La distance inter-atomique est $r_e = 0.1234708 \times 10^{-9}$ m et la taille des éléments finis du modèle continu est $h = r_e$. Pour les atomes, on utilise le potentiel de Lennard-Jones $a = 32.043529 \times 10^{-21}$ J et une masse atomique $m = 0.0016599 \times 10^{-24}$ g. Les propriétés élastiques du modèle continu sont obtenues à partir du modèle moléculaire [2]. Le calcul est effectué sur 2000 pas de temps de $\Delta t = 1 \times 10^{-15}$ s, soit 95% du pas de temps critique.

Les Figures 1 à 3 présentent les résultats. Avec la stratégie usuelle, le transfert d'énergie n'est pas correct bien que le bilan énergétique soit satisfait (Figure 3(a)). Au contraire, la stratégie proposée permet des transferts d'énergie parfait entre les deux modèles (voir [1] pour plus de détails).

## 3.2 Propagation dynamique d'une fissure dans une plaque

On étudie maintenant une plaque formée par un crystal de cuivre dans son plan crystallographique (111) (le réseau bidimensionel est ainsi hexagonal). On utilise le potentiel de Lennard-Jones pour la partie atomique avec les constantes $a = 0.415$ eV et $b = 0.2277$ nm. Les propriétés du
<figure><img src="image_6.png" /><figcaption>Figure 3 – Comparaison des transferts d’énergie pour les deux stratégies : à gauche pondération + discrétisation (9), à droite discrétisation + pondération (10).</figcaption></figure>

modèle continu sont alors obtenues en suivant [2] : $E = 79.334 \text{ GPa}$ et $v = 0.25$. La masse atomique du cuivre est $m_{Cu} = \frac{63.546}{N_A} \text{ g}$ avec $N_A$ le nombre d'Avogadro $N_A = 6.02214179 \times 10^{23} \text{ mol}^{-1}$ ce qui correspond à une densité $\rho = 1865.250812586 \times 10^3 \text{ kgm}^{-3}$. La plaque étudiée mesure 100 nm par 77.5 nm et on impose une vitesse de $47.4 \text{ ms}^{-1}$ orthogonal à la direction de la pré-fissure placée dans le modèle continu. La taille des éléments finis du modèle continus est environ 10 fois la distance inter-atomique. Ainsi, le maillage est constitué de 1221 éléments quadrangle et 8875 atomes sont placés à la pointe de la pré-fissure. L'épaisseur de la zone de couplage est de 3 nm environ et 33% des atomes qui s'y trouvent portent des multiplicateurs de Lagrange. La Figure 4 montre l'évolution du champs de déplacement et la propagation de la fissure principale ainsi que la formation de bande de cisaillement. L'influence de la vitesse de chargement est détaillée dans [3].

## 4 Conclusions

On propose ici une adaptation de la méthode Arlequin au couplage entre une modèle continu éléments finis et un modèle de dynamique moléculaire. La difficulté réside dans la nature discréte du second modèle qui introduit une dissymmétrie dans le problème couplé. Ainsi, on montre que les opérations de pondération des énergies et de discrétisation des formes faibles ne commutent pas. Dans le cas précis d'un couplage entre deux modèles dont un est discret par nature il faut d'abord discrétisé le problème continu et ensuite pondérer les énergies discrétisées. On garantit ainsi des transferts d'énergies corrects entre les deux modèles ce qui est illustré dans le premier exemple proposé. On montre ensuite une application à la propagation dynamique d'une fissure dans une plaque en couplant des éléments finis étendus et un modèle atomique uniquement à la pointe de la fissure.

## Références

[1] Aubertin P, Réthoré J, de Borst R. Energy conservation of atomistic/continuum coupling. *International Journal of Numerical Methods in Engineering* Accepted for publication.

[2] Aubertin P, Réthoré J, de Borst R. Derivation of material continuum properties based on atomistic simulations. Submitted for publication.

[3] Aubertin P, Réthoré J, de Borst R. Dynamic crack propagation using a combined molecular dynamics / extended finite element approach. To be submitted for publication.
<figure><img src="image_7.png" /></figure>

Figure 4 – Champs de déplacement autour de la pointe de la fissure modélisée en utilisant un mod- èle de dynamique moléculaire. Notons que la géomtrie du domaine moléculaire est ré-actualisée de manière à pouvoir suivre l'évolution de la fissure au cours du temps.

[4] Ben Dhia H, Rateau G. The Arlequin method as a flexible engineering design tool. *International Journal of Numerical Methods in Engineering* 2005; **62**: 1442–1462.

[5] Hughes T, Belytschko T, 2000. Nonlinear finite element analysis. ICE Division, Zace Services Ltd.

[6] Xiao SP, Belytschko T, A bridging domain method for coupling continua with molecular dynamics. *Computer Methods in Applied Mechanics and Engineering* 2004; **193**: 1645–1669.