# Energy conservation of atomistic/continuum coupling

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

► To cite this version:

Pascal Aubertin, Julien Réthoré, René de Borst. Energy conservation of atomistic/continuum coupling. International Journal for Numerical Methods in Engineering, 2009, 78 (11), pp.1365-1386. 10.1002/nme.2542. hal-00402217

HAL Id: hal-00402217

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

Submitted on 6 Apr 2021

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
Energy conservation of atomistic/continuum coupling

Pascal Aubertin¹, Julien Réthoré¹ and René de Borst², *, †

¹Université de Lyon, CNRS INSA-Lyon, LaMCoS UMR 5259, France

²Eindhoven University of Technology, P.O. Box 513, NL-5600 MB Eindhoven, The Netherlands

The efficient and accurate coupling of two dissimilar domains presents a major challenge, especially when wave propagation is considered. Overlap coupling methods are promising in the sense that spurious wave reflections can be avoided and loss of energy due to the coupling scheme can be minimized. However, the conservation properties and the proper physical representation of the forces depend on the precise formulation of the algorithm for coupling such dissimilar models. This is unlike that of coupling similar domains. We will demonstrate this with the help of numerical studies in continuum-to-continuum coupling and continuum-to-discrete coupling.

KEY WORDS: multiscale methods; coupling methods; molecular dynamics; finite elements; energy conservation

# 1. INTRODUCTION

The macroscopic behaviour of engineering materials is to a large extent determined by physical processes which occur at a scale that is one to several orders of magnitude smaller than the macroscopic scale of observation. There are several ways to incorporate this observation in engineering analyses. A first possibility is to ignore this inherent multiscale character and to use phenomenological constitutive relations that have been derived directly at the macroscopic level. Although many successful analyses have been conducted in this manner, their range of validity remains limited and their physical underpinning can be weak. To improve this situation, more recently, constitutive relations have been derived at the macroscopic level which depart from micromechanical considerations at a lower scale, followed by subsequent upscaling. Constitutive relations derived in this manner are more sound from a physics point of view and tend to have a larger domain of validity.

*Correspondence to: René de Borst, Eindhoven University of Technology, P.O. Box 513, NL-5600 MB Eindhoven, The Netherlands.

†E-mail: r.d.borst@tue.nl

Contract/grant sponsor: DGA (Délégation Générale pour l'Armement)
A further approach is to model every detail and to conduct the entire analysis on a fine scale. To circumvent the insurmountable high computational costs inherent in this approach, multiscale methods are being sought for. In these approaches a relatively small part of the domain—typically the part that is of high interest—is modelled at a lower scale, while the remainder of the domain is modelled directly at the macroscopic level. Evidently, refinements exist, where more than two levels are considered, but this is usually outside the current capabilities.

A major issue in the latter class of problems is the coupling of domains where the modelling is done at different scales, see for instance [1–4], or the book by Liu *et al.* [5]. Broadly speaking, two classes of coupling methods have been put forward: methods where the coupling is achieved at a discrete interface, and the methods where a zone of a finite size is employed, often called overlap or zonal coupling. We consider edge-to-edge coupling [6], the quasi-continuum method [7], and macroscopic atomistic *ab initio* dynamics [8] to belong to the first class of methods. The latter class, which includes approaches like the Arlequin method [9], the bridging domain method [10–14], discrete-to-continuum bridging [15], the discontinuous enrichment method [16], and bridging scale decomposition [17, 18], enables a more gradual transition from one domain to the other. It is gaining prominence not only in solid mechanics, but also in certain fluid mechanics applications, e.g. [19]. The ability of a gradual transition becomes especially important for highly dissimilar domains and when wave propagation phenomena are considered, where preservation of the energy and avoiding spurious reflections when a wave exits one domain and enters the other can become an issue. In spite of the importance of assessing the robustness of coupling schemes under these conditions, detailed investigations are relatively scarce in the literature. For this reason we have conducted detailed studies on energy conservation in overlap coupling in this contribution. Inspired by earlier work by Ben Dhia and Rateau [9] and Xiao and Belytschko [10] we have chosen a weak coupling between the models in the two adjacent domains.

A common feature of zonal or overlap coupling methods is that a weighting function is introduced to partition a certain quantity over both models within the domain where they overlap. Herein, we have chosen the energy, its major advantage being that it is the only quantity which is extrinsic to the modelling, and Lagrange multipliers are used to enforce the energy transfer between both models in the coupling domain. The precise way in which the weighting function is introduced can be of importance for the proper energy transfer. This is less critical in case of a continuum-to-continuum transition or a discrete-to-discrete transition, since similar models are then coupled. However, in the case of continuum-to-discrete transition, or vice versa, we have very different quantities: Classical notions like strains and stresses in the continuum model, but in the atomistic domain only quantities like displacements and the bonding energy make sense from a physics point of view. These differences make that the precise form of the algorithm becomes more critical if we wish to properly preserve energy and avoid the appearance of non-physical forces when a wave propagates through the coupling zone. This is of particular importance since the characteristic space and time scales in the atomistic domain are much smaller than those in the continuum domain. For an efficient multiscale analysis, it is therefore imperative that large time steps and/or a coupling zone that contains only a limited number of atoms and/or a limited number of finite elements do not destroy the accuracy and conservation properties of the coupling algorithm.

This contribution is organized as follows. First, the governing equations for a continuum and for an atomistic domain are recalled. Subsequently, the sets of equations are cast in a weak format and the partition-of-unity for the energy distribution over both models is introduced. The coupling conditions are introduced and enforced via Lagrange multipliers. Next, the spatial discretization for the continuum domain is done, the equations are cast into a matrix-vector format for both models
and for the Lagrange multipliers, and the time integration scheme is introduced. Finally, the energy preservation properties are analysed for continuum-to-continuum coupling and for continuum-to-discrete coupling. Differences are highlighted which depend on the precise formulation, and the favourable properties of one formulation in terms of energy conservation for multiscale analyses are indicated.

## 2. PROBLEM STATEMENT

### 2.1. Governing equations

We consider the domain $\Omega$ with the boundaries $\partial\Omega_1$ and $\partial\Omega_2$ on which the displacements and tractions are prescribed, respectively. The domain is thought to consist of two subdomains, $\Omega_M$ and $\Omega_m$, which are modelled using a continuum approach and a molecular dynamics approach, respectively. The boundary between both domains is assumed to be fixed, in order to study and quantify the energy conservation properties of the coupling scheme between both subdomains.

Assuming small deformation gradients for simplicity, the governing equations in the continuum subdomain $\Omega_M$ can be written in a standard manner as

For $\mathbf{x} \in \Omega_M(t)$ and $t \in [0; T]$, given the initial conditions ($\mathbf{u}(\mathbf{x}, 0), \dot{\mathbf{u}}(\mathbf{x}, 0)$)

$$ \text{find } (\mathbf{u}, \boldsymbol{\sigma}) \in \mathcal{U}^{\text{ad}} \times \mathcal{F}^{\text{ad}} \text{ such that} \quad (1) $$

$$ \rho \ddot{\mathbf{u}} = \text{div} \boldsymbol{\sigma} + \mathbf{g}_d $$

with $\mathbf{u}$ as the continuum displacement vector, $\boldsymbol{\sigma}$ the stress tensor, and $\mathbf{g}_d$ the body force vector applied in $\Omega$, subject to the boundary conditions

$$ \mathcal{U}^{\text{ad}} = \{\mathbf{u}=\mathbf{u}(\mathbf{x}, t) \in [\mathcal{H}^1(\Omega_M)]^3; \mathbf{u}=\mathbf{u}_d \text{ on } \partial_1\Omega, \forall t \in [0, T]\} \quad (2) $$

$$ \mathcal{F}^{\text{ad}} = \{\boldsymbol{\sigma}=\mathbb{K}: \nabla \mathbf{u}(\mathbf{x}, t) \in [\mathcal{L}^2(\Omega_M)]^6; \boldsymbol{\sigma} \cdot \mathbf{n}=\mathbf{F}_d \text{ on } \partial_2\Omega, \forall t \in [0, T]\} $$

where $\mathbb{K}$ is the fourth-order stiffness tensor, $\mathbf{n}$ the outward normal vector to $\partial\Omega_2$, and $\mathbf{u}_d$ and $\mathbf{F}_d$ the prescribed displacements and tractions at $\partial\Omega_1$ and $\partial\Omega_2$, respectively.

For the discrete domain, i.e. $\Omega_m$, we build a grid of $N_a$ atoms, and, accordingly, the initial value problem in this domain can be written as

For $1 \le i \le N_a(t)$ and $t \in [0; T]$, given the initial conditions ($\mathbf{d}(0), \dot{\mathbf{d}}(0)$)

$$ \text{find } (\mathbf{d}, \mathbf{f}) \in \mathcal{D}^{\text{ad}} \times \mathcal{F}^{\text{ad}} \text{ such that} \quad (3) $$

$$ m_i \ddot{\mathbf{d}}_i = \mathbf{f}_i $$

with

$$ \mathcal{D}^{\text{ad}} = \{\mathbf{d}=(\mathbf{d}_i(t))_{1 \le i \le N_a}, \forall t \in [0, T]\} \quad (4) $$

$$ \mathcal{F}^{\text{ad}} = \{\mathbf{f}=(\mathbf{f}_i(t)=-\nabla_i \mathcal{U}(\mathbf{d}(t)))_{1 \le i \le N_a}, \forall t \in [0, T]\} $$
from where it transpires that the interatomic forces are derived from a potential energy $\mathcal{U}$. **d** and **f** assemble the discrete displacements and forces of the individual atoms, respectively. The internal energy of the discrete domain can be viewed as the sum of each atomic contribution $\mathcal{U}_j$

$$ \mathcal{U} = \sum_j \mathcal{U}_j(\mathbf{d}) \qquad (5) $$

and the force **f**_i acting on atom i can be written as the sum of all elementary forces

$$ \mathbf{f}_i = -\frac{\partial \mathcal{U}}{\partial \mathbf{d}_i} = \sum_{j \neq i} \mathbf{f}_{ij} \qquad (6) $$

with

$$ \mathbf{f}_{ij} = - \left( \frac{\partial \mathcal{U}_i}{\partial d_{ij}} + \frac{\partial \mathcal{U}_j}{\partial d_{ji}} \right) \frac{\mathbf{d}_{ij}}{d_{ij}} = - \mathbf{f}_{ji} \qquad (7) $$

and

$$ d_{ij} = |\mathbf{d}_{ij}| = |\mathbf{d}_j - \mathbf{d}_i| \qquad (8) $$

We will use the simplest possible descriptions for the potential energy, taking into account only the pair interactions in the total energy. With pair potentials, we can elaborate the potential energy as

$$ \mathcal{U} = \frac{1}{2!} \sum_{i=1}^{N_a} \sum_{j=1}^{N_a} V(\mathbf{d}_i, \mathbf{d}_j) \qquad (9) $$

and the interaction force reads as

$$ \mathbf{f}_{ij} = - \frac{(V_{ij})'}{d_{ij}} \mathbf{d}_{ij} \qquad (10) $$

Examples of pair potential functions are the well-known Lennard-Jones potential

$$ V^{\text{LJ}}(r) = 4a \left[ \left( \frac{b}{r} \right)^{12} - \left( \frac{b}{r} \right)^{6} \right] \qquad (11) $$

which is a combination of attractive and repulsive interactions, and the Morse potential:

$$ V^{\text{Morse}}(r) = D(1 - e^{-\gamma(r-r_e)})^2 - D \qquad (12) $$

The constants *a*, *b*, *D*, and *γ* denote material parameters and *r*<sub>e</sub> is the interatomic equilibrium distance. We note that other potentials have been developed, which allow for a better description of metallic bonds, for example the embedded atoms method potential [20], and the present analysis can rather straightforwardly be extended to such more complicated potentials.
## 2.2. Weak formulations

We now specify the weak formulation, in order to arrive at an energetical framework and to allow for a discretization of the continuum subdomain. Then, for the continuum subdomain

$$ \forall \mathbf{v}^* \in \dot{\mathcal{U}}^{ad,0}, \text{ given the initial conditions } (\mathbf{u}(\mathbf{x}, 0), \dot{\mathbf{u}}(\mathbf{x}, 0)) $$

find $\mathbf{u} \in \mathcal{U}^{ad}$ such that

$$ \int_{\Omega_M} \rho \cdot \ddot{\mathbf{u}} \cdot \mathbf{v}^* \, d\Omega + \int_{\Omega_M} \nabla \mathbf{u} : \mathbb{K} : \nabla \mathbf{v}^* \, d\Omega = \int_{\partial_2 \Omega} \mathbf{F}_d \cdot \mathbf{v}^* \, ds + \int_{\Omega_M} \mathbf{g}_d \cdot \mathbf{v}^* \, d\Omega \quad (13) $$

with $\mathbf{v}^*$ as the test function, or written in a more compact manner as

$$ a_M(\mathbf{u}, \mathbf{v}^*) = l_M(\mathbf{v}^*) \quad (14) $$

where

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

$$ l_M(\mathbf{v}^*) = \int_{\partial_2 \Omega} \mathbf{F}_d \cdot \mathbf{v}^* \, ds + \int_{\Omega_M} \mathbf{g}_d \cdot \mathbf{v}^* \, d\Omega \quad (16) $$

In the subdomain $\Omega_M$ the weak formulation becomes

$$ \forall \mathbf{w}^* \in \dot{\mathcal{D}}^{ad,0}, \text{ given the initial conditions } (\mathbf{d}(0), \dot{\mathbf{d}}(0)) $$

find $\mathbf{d} \in \mathcal{D}^{ad}$ such that

$$ \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^* = 0 \quad (17) $$

with $\mathbf{w}^*$ as the test function, or, again in a more compact fashion

$$ a_m(\mathbf{d}, \mathbf{w}^*) = 0 \quad (18) $$

and

$$ 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^* $$

## 2.3. Coupling models

In order to establish efficient coupling between both subproblems introduced before and to specify the handshaking conditions, we enforce a weak coupling on a common zone $\Omega_c = \Omega_M \cap \Omega_m$. We consider the mechanical energy, which expresses the concept of duality between stress and strain, or loading and displacement, as the most relevant quantity. Indeed, we set out to couple two kinds of models that describe the same physical reality, but each with a different description. Dualizing the formulations and writing them in a weak format make it possible to obtain a global formulation that preserves the descriptive properties of each model, and to focus on the quantity of interest, i.e. the energy, which should not depend on the kind of the model.
<figure><img src="image_2.png" /><figcaption>Figure 1. Partition of unity for the energy distribution.</figcaption></figure>

In the handshaking zone $\Omega_c$, we enforce a velocity coupling in a weak sense, and we distribute the energy between both models via a partition of unity [9]. For this purpose we define the following functions, see Figure 1:

$$ \alpha: \Omega_M \to [0, 1] \qquad (19) $$

$$ \beta: \Omega_m \to [0, 1] \qquad (20) $$

such that

$$ \alpha(\mathbf{x}) = 1 \quad \text{for } \mathbf{x} \in \Omega_M \setminus \Omega_c $$

$$ \beta(\mathbf{x}) = 1 \quad \text{for } \mathbf{x} \in \Omega_m \setminus \Omega_c \qquad (21) $$

$$ \alpha(\mathbf{x}) + \beta(\mathbf{x}) = 1 \quad \text{for } \mathbf{x} \in \Omega_c $$

The weak formulation with the distribution of the energy then reads as:

$$ \begin{gathered} \forall (\mathbf{v}^*, \mathbf{w}^*) \in \dot{\mathcal{W}}^{ad,0} \times \mathcal{D}^{ad,0}, \text{ given the initial conditions } (\mathbf{u}(\mathbf{x}, 0), \dot{\mathbf{u}}(\mathbf{x}, 0), \mathbf{d}(0), \dot{\mathbf{d}}(0)) \\ \text{find } (\mathbf{u}, \mathbf{d}) \in \mathcal{W}^{ad} \times \mathcal{D}^{ad} \text{ such that} \\ a_{\alpha,M}(\mathbf{u}, \mathbf{v}^*) + a_{\beta,M}(\mathbf{d}, \mathbf{w}^*) = l_{\alpha,M}(\mathbf{v}^*) \end{gathered} \qquad (22) $$

### Remark 1

The modified forms $a_{\alpha,M}$, $a_{\beta,m}$, and $l_{\alpha,M}$ take into account the weighting functions $\alpha(\mathbf{x})$ and $\beta(\mathbf{x})$. However, the precise form of the weighting functions has not yet been specified.

### Remark 2

The formulation distributes the energy over both models, but does not enforce a coupling. A coupling condition still has to be introduced.

## 2.4. Coupling condition

Before we rigorously formulate the coupling condition on $\Omega_c$, we recall that:

* The displacement and velocity fields in the domains $\Omega_M$ and $\Omega_m$ have a different nature. In $\Omega_M$ we have a continuum field, whereas in $\Omega_m$ we have a discrete field, which is only defined at the geometrical points corresponding to the atoms. Thus, we have to compare two fields, $\dot{\mathbf{u}}$ and $\mathbf{d}$, on $\Omega_c$ which have a very different nature.
• The velocity coupling between both models is enforced in a weak sense, and the condition therefore has to be written in a weak format, i.e. in a ‘global’ or ‘integral’ way.

For these reasons we construct a new space, denoted by $\mathcal{M}$ and named the ‘mediator space’, on which we will project the fields $\dot{u}$ and $\dot{d}$ in order to be able to compare them. The nature of $\mathcal{M}$ is constrained by the discrete character of the atomistic field. Indeed, its displacements cannot be extrapolated outside the atomic positions if we want to maintain a physical interpretation at the fine scale. Accordingly, $\mathcal{M}$ has to be a subspace of the physical atomistic space. More precisely, we will project the velocities using an operator $\Pi$ on a discrete subset $\bar{\Omega}_c$ of the atomic positions included in $\bar{\Omega}_c$. Considering that $\mathcal{M}$ is built as a Hilbert space, we introduce a scalar product $c$ from $\mathcal{M} \times \mathcal{M}$ to $\mathbb{R}$.

With these definitions we formulate the velocity coupling as

$$ \forall \mu^* \in \mathcal{M}, \quad c(\mu^*, \Pi\dot{u} - \Pi\dot{d}) = \langle \mu^*, \Pi\dot{u} - \Pi\dot{d} \rangle_{\mathcal{M}} = 0 \quad (23) $$

with $c$ as the classical scalar product on $\mathcal{M}$:

$$ \forall (\mathbf{y}, \mathbf{z}) \in \mathcal{M} \times \mathcal{M}, \quad c(\mathbf{y}, \mathbf{z}) = \langle \mathbf{y}, \mathbf{z} \rangle_{\mathcal{M}} = \sum_{\bar{\Omega}_c} (\mathbf{y}(\mathbf{x}_i) \cdot \mathbf{z}(\mathbf{x}_i)) \quad (24) $$

The global equations are coupled via Lagrange multipliers and can subsequently be written as:

$$ \begin{gathered} \forall (\mathbf{v}^*, \mathbf{w}^*, \boldsymbol{\mu}^*) \in \dot{\mathcal{M}}^{ad,0} \times \mathcal{D}^{ad,0} \times \mathcal{M}, \quad \text{given the initial conditions } (\mathbf{u}(\mathbf{x}, 0), \dot{\mathbf{u}}(\mathbf{x}, 0), \mathbf{d}(0), \dot{\mathbf{d}}(0)) \\ \text{find } (\mathbf{u}, \mathbf{d}, \boldsymbol{\lambda}) \in \mathcal{U}^{ad} \times \mathcal{D}^{ad} \times \mathcal{M} \text{ such that} \\ a_{\alpha,M}(\mathbf{u}, \mathbf{v}^*) + a_{\beta,M}(\mathbf{d}, \mathbf{w}^*) + c(\boldsymbol{\lambda}, \Pi\mathbf{v}^* - \Pi\mathbf{w}^*) = l_{\alpha,M}(\mathbf{v}^*) \\ c(\boldsymbol{\mu}^*, \Pi\dot{u} - \Pi\dot{d}) = 0 \end{gathered} \quad (25) $$

### 3. DISCRETIZED PROBLEM

We now introduce the spatio-temporal discretizations. First, the continuum subproblem on $\Omega_M$ is discretized with a finite element method, and subsequently a time discretization is applied in order to solve the dynamic global problem.

#### 3.1. Spatial discretization for the continuum problem

With the shape functions $N_i$ and nodal vectors $\mathbf{u}_i$, which assemble the unknowns, one obtains for the subdomain $\Omega_M$

$$ \forall \mathbf{x} \in \Omega_M, \quad \mathbf{u}_h(\mathbf{x}) = \sum_{i \in \mathcal{N}_h} \mathbf{N}_i(\mathbf{x})\mathbf{u}_i = \mathbf{N}^T\mathbf{U} \quad (26) $$

With the latter, symbolic notation, the bilinear form $a_M$ and the linear form $l_M$, cf. Equations (15) and (16), become

$$ a_M(\mathbf{u}_h, \mathbf{v}_h^*) = \mathbf{V}^{*T}\mathbf{M}\ddot{\mathbf{U}} + \mathbf{V}^{*T}\mathbf{K}\mathbf{U} \quad (27) $$

$$ l_M(\mathbf{v}_h^*) = \mathbf{V}^{*T}\mathbf{F} \quad (28) $$
where the term that represents the body forces has been omitted for simplicity, and

$$ \mathbf{M} = \mathbf{N}^T \mathbf{N} \qquad (29) $$

$$ \mathbf{K} = \nabla \mathbf{N}^T \llbracket \nabla \mathbf{N} \rrbracket \qquad (30) $$

the standard mass and stiffness matrices, respectively. Introducing the weighting function $\alpha$, the bilinear and linear forms as formulated in a standard fashion in Equations (27) and (28) then become

$$ a_{\alpha, M}(\mathbf{u}_h, \mathbf{v}_h^*) = \mathbf{V}^{*T} \mathbf{M}_{\alpha} \ddot{\mathbf{U}} + \mathbf{V}^{*T} \mathbf{K}_{\alpha} \mathbf{U} \qquad (31) $$

$$ l_{\alpha, M}(\mathbf{v}_h^*) = \mathbf{V}^{*T} \mathbf{F}_{\alpha} \qquad (32) $$

where the precise form of the matrices $\mathbf{M}_{\alpha}$, $\mathbf{K}_{\alpha}$, and the vector $\mathbf{F}_{\alpha}$, which governs the energy distribution via the weighting function $\alpha$, is yet left unspecified. With the standard definition of the scalar product, the coupling term in the continuum can be discretized as follows:

$$ c(\lambda, \Pi \mathbf{v}_h^*) = \mathbf{V}^{*T} \mathbf{C}_M \Lambda = \mathbf{V}^{*T} \mathbf{F}_M^L \qquad (33) $$

with $\mathbf{C}_M$ as the continuum coupling matrix. The vector $\Lambda$ contains the Lagrange multipliers and its size equals the $\bar{\Omega}_c$ subset cardinal times the dimension of the space considered. $\mathbf{F}_M^L$ can be regarded as a fictitious force due to the coupling via the Lagrange multipliers. This force has a non-zero value only in the coupling zone $\Omega_c$.

## 3.2. Weighting function for the atomistic problem

In the subdomain $\Omega_m$ the problem is discrete by nature and there is no need for spatial discretization. We can directly use a standard matrix formulation

$$ \mathbf{m} = [m_{i,j}] = [\delta_{i,j} m_i] \qquad (34) $$

$$ \mathbf{f} = [\mathbf{f}_i] = [-\nabla_i \mathcal{U}(\mathbf{d})] \qquad (35) $$

and Equation (17) becomes:

$$ \mathbf{w}^* \cdot \mathbf{m} \ddot{\mathbf{d}} = \mathbf{w}^* \cdot \mathbf{f} \qquad (36) $$

The modified form $a_{\beta,m}$ that takes into account the distribution of the energy then reads as:

$$ a_{\beta,m}(\mathbf{d}, \mathbf{w}^*) = \mathbf{w}^* \cdot \mathbf{m}_{\beta} \ddot{\mathbf{d}} - \mathbf{w}^* \cdot \mathbf{f}_{\beta} \qquad (37) $$

with

$$ \beta_i = \beta(\mathbf{d}_i) \qquad (38) $$

$$ \mathbf{m}_{\beta} = [\beta_i \delta_{i,j} m_i] \qquad (39) $$

$$ \mathbf{f}_{\beta} = [\mathbf{f}_{\beta,i}] \qquad (40) $$

The weighted atomic forces $\mathbf{f}_{\beta,i}$ directly stem from the weighted internal energy $\mathcal{U}_{\beta}$:

$$ \mathbf{f}_{\beta,i} = - \frac{\partial \mathcal{U}_{\beta}}{\partial \mathbf{d}_i} \qquad (41) $$
Using (5) and (9), we write:

$$ \mathcal{U}_{\beta} = \sum_{k=1}^{N_a} \frac{\beta_k}{2} \sum_{j=1}^{N_a} V(\mathbf{d}_k, \mathbf{d}_j) \qquad (42) $$

The expression of the weighted interatomic forces involves the evaluation of the weighting function at the positions of the neighbours of each atom. We elaborate this by combining Equation (41) with (42) to obtain:

$$ \mathbf{f}_{\beta,i} = - \sum_{k=1}^{N_a} \frac{\beta_i + \beta_k}{2} \frac{\partial V_{ik}}{\partial \mathbf{d}_i} = \beta_i \cdot \mathbf{f}_i - \sum_{k=1}^{N_a} \frac{\beta_k - \beta_i}{2} \frac{\partial V_{ik}}{\partial \mathbf{d}_i} \qquad (43) $$

From this formulation we observe that we have two contributions to the weighted force: the first contribution relates to the original force weighted by the function $\beta$ at the current atomic position, while the residual sum stands for the weighted elementary interatomic forces. This residual force

$$ \mathbf{f}_{\beta,i}^{R} = - \sum_{k=1}^{N_a} \frac{\beta_k - \beta_i}{2} \frac{\partial V_{ik}}{\partial \mathbf{d}_i} \qquad (44) $$

can become significant when $\beta_i \neq \beta_k$.

Under some non-restrictive assumptions we will study the behaviour of this residual term. We assume, for simplicity, that the $\beta$ function is linear in the coupling overlapping zone. Furthermore, we assume that we only take into account a finite and reasonable number of neighbours for the computation of the interatomic forces. Put differently, the summation will be carried out on $N_n$ neighbours on a subset of $\Omega_m$, within a cutoff radius $r_c$. Then, the residual force $\mathbf{f}_{\beta,i}^R$ is bounded as follows:

$$ \| \mathbf{f}_{\beta,i}^R \| \le N_n \times \max \left| \frac{\beta_k - \beta_i}{2} \right| \times \max \| \mathbf{f}_{ik} \| \le N_n \times \frac{r_c}{2L_c} \times \max \| \mathbf{f}_{ik} \| \qquad (45) $$

Accordingly, when the ratio $r_c/L_c$ vanishes, which implies that the coupling domain is much larger than the atomic cutoff radius, the residual force will be negligible compared with the ‘classical’ force $\mathbf{f}_i$:

$$ \| \mathbf{f}_{\beta,i}^R \| \ll \| \mathbf{f}_i \| \qquad (46) $$

Henceforth, we will assume that this condition is satisfied and we will use this approximation for the weighted forces:

$$ \mathbf{f}_{\beta,i} = \beta_i \mathbf{f}_i \qquad (47) $$

### 3.3. Coupled system

We next write the coupling force using the Lagrange multipliers in the atomistic domain similar to that in the continuum domain

$$ \mathbf{f}_m^L = \mathbf{C}_m \Lambda \qquad (48) $$

whereupon the weighted and coupled system (25) can be cast in a matrix-vector format:

$$ \mathbf{V}^{*T}(\mathbf{M}_{\alpha}\ddot{\mathbf{U}} + \mathbf{K}_{\alpha}\mathbf{U} + \mathbf{C}_M\Lambda) + \mathbf{W}^{*T}(\mathbf{m}_{\beta}\ddot{\mathbf{d}} - \mathbf{f}_{\beta} - \mathbf{C}_m\Lambda) = \mathbf{V}^{*T}\mathbf{F}_{\alpha} \qquad (49) $$

$$ \mu^{*T}(\mathbf{C}_M^T \dot{\mathbf{U}} - \mathbf{C}_m^T \dot{\mathbf{d}}) = 0 \qquad (50) $$
Since this set must hold for any admissible ($V^*, W^*, \mu^*$), we finally obtain

$$ 
\begin{aligned}
\mathbf{M}_{\alpha} \ddot{\mathbf{U}} + \mathbf{K}_{\alpha} \mathbf{U} &= \mathbf{F}_{\alpha} - \mathbf{C}_{M} \Lambda \\
\mathbf{m}_{\beta} \ddot{\mathbf{d}} &= \mathbf{f}_{\beta} + \mathbf{C}_{m} \Lambda \\
\mathbf{C}_{M}^{\mathrm{T}} \dot{\mathbf{U}} &= \mathbf{C}_{m}^{\mathrm{T}} \dot{\mathbf{d}}
\end{aligned}
\qquad (51) $$

with ($U, d, \Lambda$) the set of unknowns.

### 3.4. Time integration scheme

The time scheme relies on the discretization with a time step $\Delta t$ and has five stages:

* Given the quantities at step $n$, compute the displacements $U_{n+1}$ and $d_{n+1}$.

* Compute the predictive accelerations $\ddot{\mathbf{U}}_{n+1}^*$ and $\ddot{\mathbf{d}}_{n+1}^*$, neglecting the Lagrange forces.

* Compute the predictive velocities $\dot{\mathbf{U}}_{n+1}^*$ and $\dot{\mathbf{d}}_{n+1}^*$.

* Adjust these velocities to give the final velocities $\dot{\mathbf{U}}_{n+1}$ and $\dot{\mathbf{d}}_{n+1}$ by taking into account the coupling terms and Lagrange multipliers $\Lambda_{n+1}$.

* Adjust the predictive accelerations to give the final accelerations $\ddot{\mathbf{U}}_{n+1}$ and $\ddot{\mathbf{d}}_{n+1}$.

Below we specify the different steps of this predictor-corrector scheme:

* With the displacements, velocities and accelerations at step $n$, we compute the displacements at step $n+1$ as follows:

$$ 
\begin{aligned}
\mathbf{U}_{n+1} &= \mathbf{U}_n + \dot{\mathbf{U}}_n \Delta t + \frac{1}{2} \ddot{\mathbf{U}}_n (\Delta t)^2 \\
\mathbf{d}_{n+1} &= \mathbf{d}_n + \dot{\mathbf{d}}_n \Delta t + \frac{1}{2} \ddot{\mathbf{d}}_n (\Delta t)^2
\end{aligned}
\qquad (52) $$

* The predictive accelerations at step $n+1$ are computed with the system (51) but without coupling terms:

$$ 
\begin{aligned}
\ddot{\mathbf{U}}_{n+1}^* &= \tilde{\mathbf{M}}_{\alpha}^{-1} (\mathbf{F}_{\alpha, n+1} - \mathbf{K}_{\alpha} \mathbf{U}_{n+1}) \\
\ddot{\mathbf{d}}_{n+1}^* &= \mathbf{m}_{\beta}^{-1} \mathbf{f}_{\beta, n+1}
\end{aligned}
\qquad (53) $$

Note that we use, for the continuum, a lumped mass matrix $\tilde{\mathbf{M}}_{\alpha}$, which is standard for explicit time integration.

* The predictive velocities are computed with the Newmark scheme:

$$ 
\begin{aligned}
\dot{\mathbf{U}}_{n+1}^* &= \dot{\mathbf{U}}_n + \frac{1}{2}(\ddot{\mathbf{U}}_n + \ddot{\mathbf{U}}_{n+1}^*) \Delta t \\
\dot{\mathbf{d}}_{n+1}^* &= \dot{\mathbf{d}}_n + \frac{1}{2}(\ddot{\mathbf{d}}_n + \ddot{\mathbf{d}}_{n+1}^*) \Delta t
\end{aligned}
\qquad (54) $$

* We next adjust the velocities by introducing the coupling terms:

$$ 
\begin{aligned}
\dot{\mathbf{U}}_{n+1} &= \dot{\mathbf{U}}_{n+1}^* - \frac{1}{2} \tilde{\mathbf{M}}_{\alpha}^{-1} \mathbf{F}_{M, n+1}^{\mathrm{L}} \Delta t \\
\dot{\mathbf{d}}_{n+1} &= \dot{\mathbf{d}}_{n+1}^* + \frac{1}{2} \mathbf{m}_{\beta}^{-1} \mathbf{f}_{m, n+1}^{\mathrm{L}} \Delta t
\end{aligned}
\qquad (55) $$
• Finally, we compute the accelerations:

$$ \ddot{U}_{n+1} = \ddot{U}_{n+1}^* - \tilde{M}_{\alpha}^{-1} \mathbf{F}_{M,n+1}^{\mathrm{L}} \quad (56) $$

$$ \ddot{d}_{n+1} = \ddot{d}_{n+1}^* + \mathbf{m}_{\beta}^{-1} \mathbf{f}_{m,n+1}^{\mathrm{L}} $$

The first three steps are simple and do not need further explanation. The last steps enforce the coupling condition (23).

From Equations (23) and (51), the coupling condition becomes:

$$ C_M^T (\dot{U}_{n+1}^* - \tilde{M}_{\alpha}^{-1} C_M \Lambda_{n+1} \Delta t) = C_m^T (\dot{d}_{n+1}^* + \mathbf{m}_{\beta}^{-1} C_m \Lambda_{n+1} \Delta t) \quad (57) $$

The new values of the Lagrange multipliers $\Lambda_{n+1}$ are subsequently computed by solving

$$ A \Lambda_{n+1} = b_{n+1} \quad (58) $$

with

$$ A = \frac{\Delta t}{2} (C_M^T \tilde{M}_{\alpha}^{-1} C_M + C_m^T \mathbf{m}_{\beta}^{-1} C_m) \quad (59) $$

$$ b_{n+1} = C_M^T \dot{U}_{n+1}^* - C_m^T \dot{d}_{n+1}^* $$

and $b_{n+1}$ stands for the weak coupling condition on the predictive velocities. Thus, this term is a measure of the error compared with the solution that satisfies the system (51).

### 3.5. Energy balance of the resolution scheme

We now focus on the energy conservation in the coupled system. Indeed, a most important aspect of a coupling method is that it conserves the global energy and ensures a proper transfer of the relevant quantities from one domain to the other. We assume, for simplicity, that there is no external force acting on the continuum subdomain. We will use the following notations for the mean value and the jump of a quantity $V$ between $t_n$ and $t_{n+1}$ as:

$$ \langle \mathbf{V} \rangle = \frac{1}{2} (\mathbf{V}_{n+1} + \mathbf{V}_n) \quad (60) $$

$$ [\mathbf{V}] = \mathbf{V}_{n+1} - \mathbf{V}_n $$

where we have the following property:

$$ \langle \mathbf{V} \rangle^T [\mathbf{V}] = \frac{1}{2} [\mathbf{V}^T \mathbf{V}] \quad (61) $$

Thus, the Newmark scheme described in the previous subsection can be written more compactly as:

$$ [\mathbf{U}] = \langle \dot{\mathbf{U}} \rangle \Delta t - \frac{(\Delta t)^2}{4} [\ddot{\mathbf{U}}] \quad (62) $$

$$ [\dot{\mathbf{U}}] = \langle \ddot{\mathbf{U}} \rangle \Delta t \quad (63) $$

$$ [\mathbf{d}] = \langle \dot{\mathbf{d}} \rangle \Delta t - \frac{(\Delta t)^2}{4} [\ddot{\mathbf{d}}] \quad (64) $$

$$ [\dot{\mathbf{d}}] = \langle \ddot{\mathbf{d}} \rangle \Delta t \quad (65) $$
We use the system (51) to compute the mean value of the equilibrium equations between time step $t_n$ and $t_{n+1}$ and we multiply by $[U]$ for the continuum subdomain and by $[d]$ for the atomistic subdomain to obtain the following discrete sets:

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

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

where the following quantities can be identified: (i) the kinetic energy $E_K^M = \frac{1}{2} \dot{\mathbf{U}}^{\mathrm{T}} \tilde{\mathbf{M}}_{\alpha} \dot{\mathbf{U}}$ in the continuum; (ii) the internal energy $E_P^M = \frac{1}{2} \mathbf{U}^{\mathrm{T}} \mathbf{K}_{\alpha} \mathbf{U}$ in the continuum; (iii) the kinetic energy $E_K^m = \frac{1}{2} \dot{\mathbf{d}}^{\mathrm{T}} \mathbf{m}_{\alpha} \dot{\mathbf{d}}$ in the atoms; (iv) the acceleration terms $((\Delta t)^2/(8)) \dot{\mathbf{U}}^{\mathrm{T}} \tilde{\mathbf{M}}_{\alpha} \ddot{\mathbf{U}}$ and $((\Delta t)^2/(8)) \ddot{\mathbf{d}}^{\mathrm{T}} \mathbf{m}_{\alpha} \ddot{\mathbf{d}}$ due to the central difference scheme; (v) the work of the Lagrange multipliers related to the displacements of the continuum and to the displacements of the atoms, $[\mathbf{U}]^{\mathrm{T}} \langle \mathbf{C}_M \mathbf{\Lambda} \rangle$ and $[\mathbf{d}]^{\mathrm{T}} \langle \mathbf{C}_m \mathbf{\Lambda} \rangle$, respectively; and (vi) the work of the interatomic forces within the atomic displacements $[\mathbf{d}]^{\mathrm{T}} \langle \mathbf{f}_{\beta} \rangle$.

### Remark 1

We observe from this discretized energy balance that the relevant quantity for the atomistic domain is the work of the interatomic forces. Indeed, we could have used the total potential energy, but, within this discretized scheme, we then would have had to compute the work of the atomic forces.

### Remark 2

The precise form of the different matrices has still been left open. We will specify and investigate different forms in the remainder.

## 4. ANALYSIS OF THE COUPLING SCHEME

In the present section we will analyse the coupling scheme with respect to energy conservation. First, we will consider the coupling between two continuum subdomains and show that the precise formulation is less critical, in particular with respect to energy conservation. Subsequently, we will demonstrate that this is not the case when two dissimilar subdomains are coupled, e.g. a continuum subdomain and a subdomain that is inherently discrete.

### 4.1. Coupling of two continuum subdomains

At first, we compute the lumped mass matrix, cf. Equation (53), as follows:

$$ \tilde{\mathbf{M}}_{\alpha,I} = \alpha(X_I) \tilde{\mathbf{M}}_I \quad (68) $$

where $X_I$ is the position of the Ith node and $\tilde{\mathbf{M}}$ is the classical lumped mass matrix. Similarly, the stiffness matrix is computed by weighting the elementary terms, thus taking into account the influence of the weighting function $\alpha$ [13]:

$$ \mathbf{K}_{\alpha,I,J} = \int_{\Omega} \alpha(X) \nabla \mathbf{u}(\mathbf{N}_I(X)) : \mathbb{K} : \nabla \mathbf{u}(\mathbf{N}_J(X)) d\Omega \quad (69) $$
<figure><img src="image_3.png" /><figcaption>Figure 2. Propagating traction wave in two continuum subdomains: (a) displacements at $t=0$ and (b) displacements at $t=0.1837244 \times 10^{-12}$ s.</figcaption></figure>

We consider two longitudinal bars, each modelled with finite elements, with a coupling region $L_c$ in between. The left bar is submitted to a traction wave, the initial configuration being displaced on the left-most 50 elements. The right bar is fixed at its right end. The entire domain is $(94.337 \times 10^{-9} \text{ m})$ long and each bar has been modelled using 150 elements. The computation extends over 2500 time steps each of $\Delta t = 1.2248 \times 10^{-15}$ s, which corresponds to 99% of the critical time step. Figure 2(a) shows the initial displacements in the first bar. Figure 2(b) shows how in the second bar the travelling wave is captured in the coupling zone.

We now focus on the mechanical energy transfer when the wave passes through the coupling zone from the first domain to the second domain. We consider three cases in which the coupling length $L_c$ contains 2, 10, and 40 elements, respectively, see Figure 3. In each case, the mechanical energy, calculated as described in (66), passes from one domain to the other without significant loss. We just notice that at each crossing of the coupling domain, the total energy oscillates before attaining its correct value. This phenomenon corresponds to the energy that is temporarily stored
<figure><img src="image_4.png" /><figcaption>Figure 3. Energy plots for the finite element–finite element coupling. In subfigures (a), (c), and (e), the drawn line is the energy in the first (left) bar and the dashed line is the energy in the second (right) bar. The bold drawn line is the total mechanical energy: (a) energy transfer for $L_c$ consisting of 2 elements; (b) oscillations in the mechanical energy for $L_c$ consisting of 2 elements; (c) energy transfer for $L_c$ consisting of 10 elements; (d) oscillations in the mechanical energy for $L_c$ consisting of 10 elements; (e) energy transfer for $L_c$ consisting of 40 elements; and (f) oscillations in the mechanical energy for $L_c$ consisting of 40 elements.</figcaption></figure>

in the Lagrange multipliers. Indeed, as argued in the preceding section, there is an amount of work stored in the Lagrange multipliers when the wave passes the coupling zone. We will return to this issue in more detail in the next subsection.

Next, we investigate the behaviour of the method when we use the standard mass and stiffness matrices in the prediction step, i.e. without multiplying them by the weight function. The second step of the time integration scheme is the computation of the equilibrium based on the prediction step, Equation (53). As we do not use the Lagrange multipliers in the prediction step, we assume that both subdomains are not coupled during this step and, thus, that both subdomains are free.
<figure><img src="image_5.png" /><figcaption>Figure 4. Energy plots with standard matrices, i.e. without weighting in the predictor phase. In subfigure (a), the drawn line is the energy in the first (left) bar and the dashed line represents the energy in the second (right) bar. The bold drawn line is the total mechanical energy: (a) energy transfer for $L_c$ consisting of 10 elements and (b) oscillations in the mechanical energy for $L_c$ consisting of 10 elements.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>Figure 5. Displacements at $t=0$. The left and right sides constitute the continuum domains, while the inside zone covers the atomistic region.</figcaption></figure>

Indeed, we first solve the ‘uncoupled’ equations and subsequently, we link them using the Lagrange multipliers. The predictive accelerations are then computed as

$$ \ddot{U}_{n+1}^* = \tilde{M}^{-1}(F_{n+1} - KU_{n+1}) \qquad (70) $$

This provides the initial solutions that have to be coupled using the Lagrange multipliers. For this subsequent coupling step, it is necessary to use the weighted mass matrix, cf. Equation (59). Now, we obtain the plots for the mechanical energy, which are given in Figure 4(a). The energy conservation is also satisfactory, the oscillations being slightly different, but the transfer is very good and looks similar to that when using the weighted stiffness matrix also in the predictive step, see Figure 4(a), with the same order of magnitude.
<figure><img src="image_7.png" /></figure>

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

<figure><img src="image_9.png" /><figcaption>Figure 6. Wave propagation in a one-dimensional bar for the prediction step with weighted matrices. The left and right sides are continuum domains, while the inside zone represents the atomistic region: (a) displacements at $t=60 \times 10^{-15}$ s; (b) displacements at $t=105 \times 10^{-15}$ s; and (c) displacements at $t=150 \times 10^{-15}$ s.</figcaption></figure>

## 4.2. Coupling atomistics to a continuum

In this subsection we focus on the coupling between a continuum and an atomistic domain, applying the same methodologies as before, and using again a one-dimensional bar to carry out a quantitative assessment. We now consider a longitudinal bar described with finite elements, containing an
atomistic region, with a coupling zone on both sides. The bar is submitted to a traction wave, which is enforced by displacing the left-most 20 elements in the initial configuration, see Figure 5. The right-hand end is free. The whole domain is $59.142528 \times 10^{-9}$ m long and 100 atoms have been put in the atomistic domain. The interatomic distance is $r_e=0.1234708 \times 10^{-9}$ m and the finite element size is $h=r_e$. We use a Lennard-Jones potential as constitutive model for the atoms, with $a=32.043529 \times 10^{-21}$ J and a mass $m=0.0016599 \times 10^{-24}$ g. The elastic material properties for the finite element model have been derived from the atomic properties [21]. The computation continues for 2000 time steps with $\Delta t=1 \times 10^{-15}$ s, which amounts to 95% of the critical time step.

As for the continuum-continuum problem, we first analyse the problem using weighted matrices for the predictive part of the solution strategy, cf. Equations (68) and (69). The coupling length $L_c$ includes 5 elements. The displacements are shown for different times in Figure 6. When the wave passes through the coupling domain, we obtain an amplification of the amplitude and a non-negligible reflection. Considering the energy plots, Figure 7, we observe that, even though the total energy is conserved—which is logical within this integration scheme—the energy transfer between the domains is poor. Many reflections occur, causing information to be lost completely in the end.

An explanation of this phenomenon resides in the discretization of the initial-value problem. The introduction of the weighting in the elementary terms of the internal forces in the continuum is tantamount to considering a porous material where the porosity changes progressively with the distance. At the right-hand side, i.e. near the end of the coupling zone, the left continuum bar has a near-zero Young's modulus and a near-zero mass. This construction does not cause any spurious reflections in the case of a coupling between two continuum domains, but it introduces a non-symmetric problem between a continuum and a domain which is composed of atoms. The atomistic internal forces have been weighted with the weight function, cf. Equation (47), and to solve the problem the continuum internal forces must be constructed such that a full symmetry will exist between both domains.

<figure><img src="image_10.png" /><figcaption>Figure 7. Energy transfer for the finite element—molecular dynamics coupling. The drawn line is the energy in the first (left) continuum domain. The dashed line is the energy in the second (right) continuum domain, and the dash-dotted line represents the energy in the atomistic zone. The bold drawn line is the total mechanical energy.</figcaption></figure>
<figure><img src="image_11.png" /><figcaption>Figure 8. Wave propagation in a one-dimensional bar. Weighted matrices. The left and right sides are continuum domains, while the inside zone is the atomistic region: (a) displacements at $t = 105 \times 10^{-15}$ s; (b) displacements at $t = 150 \times 10^{-15}$ s; and (c) displacements at $t = 240 \times 10^{-15}$ s.</figcaption></figure>

The original continuum internal forces were expressed as:

$$F_{\alpha}=K_{\alpha}U \qquad (71)$$

In order to follow the same procedure as the atomistic one, we write:

$$F_{\alpha,I}=\alpha(X_{I})KU \qquad (72)$$
<figure><img src="image_12.png" /><figcaption>Figure 9. Mechanical energy (drawn line) and the work stored in the Lagrange multipliers (dashed line).</figcaption></figure>

By writing the internal forces in this manner, symmetry has been restored. Moreover, we have obtained a simpler solution scheme. Indeed, by writing the continuum internal forces in this manner, the second step of the procedure is equivalent to solving the ‘free’ problems separately, since

$$ \ddot{\mathbf{U}}_{n+1}^* = -\tilde{\mathbf{M}}_{\alpha}^{-1} \mathbf{K}_{\alpha} \mathbf{U}_{n+1} = -\tilde{\mathbf{M}}^{-1} \alpha^{-1} \alpha \mathbf{K} \mathbf{U}_{n+1} \quad (73) $$

$$ \ddot{\mathbf{d}}_{n+1}^* = \mathbf{m}_{\beta}^{-1} \mathbf{f}_{\beta, n+1} = \mathbf{m}^{-1} \beta^{-1} \beta \mathbf{f}_{n+1} $$

where $\alpha$ and $\beta$ are the diagonal weighting matrices (e.g. $\alpha_I = \alpha(X_I)$). Finally, we obtain:

$$ \ddot{\mathbf{U}}_{n+1}^* = -\tilde{\mathbf{M}}^{-1} \mathbf{K} \mathbf{U}_{n+1} \quad (74) $$

$$ \ddot{\mathbf{d}}_{n+1}^* = \mathbf{m}^{-1} \mathbf{f}_{n+1} $$

The problem is therefore highly simplified. We first solve the ‘free’ uncoupled problems separately, and then, as a second step, we couple them using the Lagrange multipliers. For this step we use the weighted mass matrices.

We now apply this procedure to the example. The displacements are shown for different times in Figure 8.

The simulations do not reveal spurious reflections anymore when the wave passes through the coupling zone, and the information at $t=0$ is preserved during the computation. Figure 9 shows the total mechanical energy during the computation. We observe some fluctuations each time the wave crosses a coupling zone. In fact, work is stored by the Lagrange multipliers, and subsequently put back in the mechanical system when the wave leaves the coupling zone. We clearly observe that the work of the Lagrange multipliers is complementary to the mechanical energy, and the energy balance formulated in Equations (66)–(67) is thus satisfied. Considering the energy plots of Figure 10 we observe that for different coupling lengths (i) the total energy is preserved, and that (ii) the energy correctly passes from one domain to the other when the wave traverses the coupling zones.
<figure><img src="image_13.png" /><figcaption>Figure 10. Energy plots for the finite element—molecular dynamics coupling. In subfigures (a), (c), and (e) the drawn line is the energy in the first (left) continuum, the dashed line is the energy in the second (right) continuum, and the dash-dotted line represents the energy stored in the atoms. The bold drawn line is the total mechanical energy: (a) energy transfer for $L_c$ consisting of 2 elements; (b) energy drift for $L_c$ consisting of 2 elements; (c) energy transfer for $L_c$ consisting of 5 elements; (d) energy drift for $L_c$ consisting of 5 elements; (e) energy transfer for $L_c$ consisting of 10 elements; and (f) energy drift for $L_c$ consisting of 10 elements.</figcaption></figure>
<figure><img src="image_14.png" /><figcaption>Figure 11. Energy plots for the finite element—molecular dynamics coupling in the multiscale example. In subfigure (a), the drawn line is the energy in the first (left) continuum, the dashed line is the energy in the second (right) continuum, and the dash-dotted line is the energy stored in the atoms. The bold drawn line is the total mechanical energy: (a) energy transfer for $L_c$ consisting of 5 elements and (b) energy loss for $L_c$ consisting of 5 elements.</figcaption></figure>

Finally, we present some calculations that demonstrate the potential for truly multiscale computations. Instead of choosing a finite element length $h$ that is of the same order of magnitude as the interatomic distance $r_e$, we put 20 atoms in each finite element. The coupling length covers 5 elements and we put 26 Lagrange multipliers in this zone, at atomic positions. Again, we observe that the energy balance is satisfied, Figure 11(a), and the total energy loss is very small, Figure 11(b).

## 5. CONCLUDING REMARKS

In this contribution, we have discussed methods for coupling two subdomains, which are typical of multiscale analyses. A most important criterion for the accuracy and robustness of any coupling method is that it conserves energy, in particular when a propagating wave passes from one domain to the other. Zonal or overlap coupling methods seem better suited for this task, since, in principle, they have a better capability to avoid spurious wave reflections.

We have proposed a class of overlap coupling methods and have analyzed the energy conservation properties when a wave passes from one subdomain to another. Like other proposals in the literature [9, 10] it enforces a weak kinematic coupling. Herein, the energy is used to enforce the coupling via Lagrange multipliers that are active in the overlapping region.

The energy transfer has been analysed for two different cases: Coupling similar domains (continuum to continuum) and coupling dissimilar domains (continuum to discrete). The former case is rather insensitive to the precise formulation of the predictor-corrector algorithm. It does not make much difference whether the weighting function that distributes the energy in the overlapping domain is included in the predictor step or not. The situation is very different, however, when dissimilar domains are considered. Then, to avoid spurious reflections and a proper energy transfer it is essential that weighting functions are only used in the corrector phase when both subdomains are coupled. This is tantamount to first discretizing the continuum subdomain prior to multiplication by the weighting function, rather than the reverse procedure. The resulting
algorithm has potential for truly multiscale analyses, where the size of the finite elements significantly exceeds the interatomic distance, while coupling zones can be used which contain only a few finite elements.

## ACKNOWLEDGEMENTS

We thank the DGA (Délégation Générale pour l'Armement) for the financial support to the first author.

## REFERENCES

1. Karpov EG, Park HS, Liu WK. A phonon heat bath approach for the atomistic and multiscale simulation of solids. *International Journal for Numerical Methods in Engineering* 2007; **70**:351–378.

2. Xiao SP, Yang WX. A temperature-related homogenization technique and its implementation in the meshfree particle method for nanoscale simulations. *International Journal for Numerical Methods in Engineering* 2007; **69**:2099–2125.

3. Fish J, Nuggehally MA, Shephard MS, Picu CR, Badia S, Parks ML, Gunzburger M. Concurrent ATC coupling based on a blend of the continuum stress and the atomistic force. *Computer Methods in Applied Mechanics and Engineering* 2007; **196**:4548–4560.

4. Nuggehally MA, Shephard MS, Picu CR, Fish J. Adaptive model selection procedure for concurrent multiscale problems. *International Journal of Multiscale Computational Engineering* 2007; **5**:369–386.

5. Liu WK, Karpov EG, Park HS. *Nano Mechanics and Materials*. Wiley: Chichester, 2006.

6. Belytschko T, Xiao SP. Coupling methods for continuum model with molecular model. *International Journal of Multiscale Computational Engineering* 2003; **1**:115–126.

7. Tadmor EB, Ortiz M, Phillips R. Quasicontinuum analysis of defects in solids. *Philosophical Magazine* 1996; **73**:1529–1563.

8. Broughton JQ, Abraham FF, Bernstein N, Kaxiras E. Concurrent coupling of length scales: methodology and application. *Physical Review B* 1999; **60**:2391–2403.

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

10. 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.

11. Zhang SL, Khare R, Lu Q, Belytschko T. A bridging domain and strain computation method for coupled atomistic-continuum modelling of solids. *International Journal for Numerical Methods in Engineering* 2007; **70**:913–933.

12. Guidault PA, Belytschko T. On the L2 and the H1 couplings for an overlapping domain decomposition method using Lagrange multipliers. *International Journal for Numerical Methods in Engineering* 2007; **70**:322–350.

13. Xu M, Belytschko T. Conservation properties of the bridging domain method FOD coupled molecular/continuum dynamics. *International Journal for Numerical Methods in Engineering* 2008; **76**:278–294.

14. Guidault PA, Belytschko T. Bridging domain methods for coupled atomistic-continuum models with L2 or H1 couplings. *International Journal for Numerical Methods in Engineering*, DOI: 10.1002/nme.2461.

15. Fish J, Chen W. Discrete-to-continuum bridging based on multigrid principles. *Computer Methods in Applied Mechanics and Engineering* 2004; **193**:1693–1711.

16. Farhat C, Harari I, Hetmaniuk U. The discontinuous enrichment method for multiscale analysis. *Computer Methods in Applied Mechanics and Engineering* 2003; **192**:3195–3209.

17. Wagner GJ, Liu WK. Coupling of atomistic and continuum simulations using a bridging scale decomposition. *Journal of Computational Physics* 2003; **190**:249–274.

18. Farrell DE, Park HS, Liu WK. Implementation aspects of the bridging scale method and application to intersonic crack propagation. *International Journal for Numerical Methods in Engineering* 2007; **71**:583–605.

19. Akkerman I, Hulshoff SJ, de Borst R. Multiscale overlap coupling for hybrid computations. *Computer Methods in Applied Mechanics and Engineering* 2007; **196**:4294–4304.

20. Daw MS. Model of metallic cohesion: the embedded-atom method. *Physical Review B* 1989; **39**:7441–7452.

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