# A molecular dynamics/extended finite element method for dynamic crack propagation

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

► To cite this version:

Pascal Aubertin, Julien Réthoré, René de Borst. A molecular dynamics/extended finite element method for dynamic crack propagation. Vietnam Journal of Mechanics, 2008, 30, pp.205-217. 10.15625/0866-7136/30/4/5626. hal-00938630

HAL Id: hal-00938630

https://hal.science/hal-00938630

Submitted on 27 May 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 Attribution 4.0 International License
# A MOLECULAR DYNAMICS / EXTENDED FINITE ELEMENT METHOD FOR DYNAMIC CRACK PROPAGATION

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

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

²Eindhoven University of Technology, Eindhoven, Netherlands

**Abstract.** A multiscale method is presented which couples a molecular dynamics approach for describing fracture at the crack tip with an extended finite element method for discretizing the remainder of the domain. After recalling the basic equations of molecular dynamics and continuum mechanics the discretization is discussed for the continuum subdomain where the partition-of-unity property of finite element shape functions is used, since in this fashion the crack in the wake of its tip is naturally modelled as a traction-free discontinuity. Next, the zonal coupling method between the atomistic and continuum models is described, including an assessment of the energy transfer between both domains for a one-dimensional problem. Finally, a two-dimensional computation is presented of dynamic fracture using the coupled model.

**Keywords.** multiscale methods, molecular dynamics, extended finite element method, fracture, crack propagation

## 1. INTRODUCTION

Modern research into fracture commences with the seminal work of Griffith [1]. Later, Irwin [2] and Rice [3] established the relation between the stress intensity factors and the energy release rate, and gave linear elastic fracture mechanics a firm basis. However, linear elastic fracture mechanics only applies to crack-like flaws in an otherwise linear elastic solid and when the singularity associated with that flaw is characterized by a non-vanishing energy release rate. The fracture and any dissipative processes must also be confined to a small region in the vicinity of the crack tip.

When the region in which the separation and dissipative process take place is not small compared to a structural dimension, but any nonlinearity is confined to a surface emanating from a classical crack tip, i.e. one with a non-vanishing energy release rate, cohesive zone models as introduced by Barenblatt [4] and Dugdale [5] apply. The cohesive zone approach was extended by Hillerborg et al. [6] and Needleman [7] to circumstances where: (i) an initial crack-like flaw need not be present or, if one is present, it need not be associated with a non-vanishing energy release rate; and (ii) non-linear deformation behaviour may occur over an extended volume. Initially, cohesive-zone models were incorporate in finite element methods via special-purpose interface elements [8, 9], but more
recently, partition-of-unity finite element methods have shown to very amenable to the incorporation of cohesive-zone models, e.g. [10–12]. In particular, they naturally enable crack propagation, also in dynamics [13–17] and in multi-phase continua [18, 19].

In spite of the power of the cohesive-zone approach, and its wide applicability on a range of scales, it remains a phenomenological approach. Probably, quantum mechanics is physically the most appropriate theory to describe fracture, but the difficulties to relate quantum mechanics to continuum mechanics, e.g. via Density Functional Theory [20, 21] presently seem insurmountable. One scale of observation higher is to use Molecular Dynamics to describe fracture processes from a more fundamental physics point of view. Indeed, researchers have recently used this approach to describe fracture, e.g. [22–24]. A disadvantage of the approach is that it is computationally demanding. For this reason multi-scale approaches have been introduced, in fracture [25], as well as in plasticity [26], in which only a part of the body is analysed using molecular dynamics, while the remaining part of the body is modelled using continuum mechanics and discretized using a finite element method. This manuscript furthers along this line and combines molecular dynamics for modelling the fracture process at the crack tip with an extended finite element method (XFEM), where the partition-of-unity property of the polynomial shape functions is exploited to model the crack in the wake of the tip as a traction-free discontinuity. It is noted that recently another approach has been published that couples atomistics and extended finite elements [27], but the current paper makes a further advancement in that it includes dynamic loadings.

A major issue in multi-scale approaches as discussed above is the accurate coupling of both domains, especially when different descriptions are assumed on either domain. While the coupling can, in principle, either be achieved at a discrete interface, or on a zone of a finite size (overlap or zonal coupling), it is believed that zonal approaches, which include the Arlequin method [28], the bridging domain method [29–32], discrete-to-continuum bridging [33], the discontinuous enrichment method [34], and bridging scale decomposition [35, 36] enable a more gradual transition from one domain to the other. The ability of a gradual transition is 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. Inspired by earlier work by Ben Dhia and Rateau [28] and Xiao and Belytschko [29] we have chosen a weak coupling between the models in the two adjacent domains.

This paper is organized as follows. First, we briefly list the equations of molecular dynamics. This is followed by a succinct recapitulating of the governing equations of continuum mechanics, both in the strong and the weak forms. The discretization of the continuum subdomain is carried out using the extended finite element method, where the partition-of-unity property of finite element shape functions is used to model the traction-free discontinuity in the wake of the crack tip. Subsequently, it is discussed how both domains can be coupled, see also [37] and an analysis is presented of the energy conservation properties of the coupling scheme. The paper concludes with a full two-dimensional
coupled analysis of dynamic crack propagation which shows multiple branching and suggests the formation of dislocations, which shows similarities with recent simulations on dynamic fracture using cohesive-zone models [9, 17].

## 2. MOLECULAR DYNAMICS

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))$, find $(\mathbf{d}, \mathbf{f}) \in \mathcal{D}^{ad} \times \mathcal{F}^{ad}$ such that:

$$m_i \ddot{\mathbf{d}}_i = \mathbf{f}_i \qquad (1)$$

with $m_i$ the mass of atom i and:

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

$$\mathcal{F}^{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}$. $\mathbf{d}$ and $\mathbf{f}$ assemble the discrete displacements $\mathbf{d}_i$ and forces $\mathbf{f}_i$ 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 (3)$$

and the force $\mathbf{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 (4)$$

In order to limit the cost of computing such a force, we reduce the summation by only including so-called “nearest” neighbours, within a cut-off radius $r_c$:

$$\mathbf{f}_i \approx \sum_{r_{ij} < r_c} \mathbf{f}_{ij} \qquad (5)$$

where $r_{ij}$ is the interatomic distance.

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)), \\ \text{find } \mathbf{d} \in \mathcal{D}^{ad} \text{ such that:}$$

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

with $\mathbf{w}^*$ the test function, 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^*$$
### 3. CONTINUUM MODEL

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

$$ \text{For } \mathbf{x} \in \Omega_M(t) \text{ and } t \in [0; T], \text{ given the initial conditions } (\mathbf{u}(\mathbf{x}, 0), \dot{\mathbf{u}}(\mathbf{x}, 0)), \\ \text{find } (\mathbf{u}, \boldsymbol{\sigma}) \in \mathcal{U}^{ad} \times \mathcal{S}^{ad} \text{ such that:} \tag{7} $$

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

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

$$ \mathcal{U}^{ad} = \left\{ \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] \right\} \tag{8} $$

$$ \mathcal{S}^{ad} = \left\{ \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] \right\} $$

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.

To allow for a discretization of the continuum subdomain we next specify the weak formulation:

$$ \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)), \\ \text{find } \mathbf{u} \in \mathcal{U}^{ad} \text{ such that:} $$

$$ a_M(\mathbf{u}, \mathbf{v}^*) = l_M(\mathbf{v}^*) \tag{9} $$

with $\mathbf{v}^*$ the test function, and

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

$$ 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 \tag{11} $$

### 4. COUPLING SCHEME

#### 4.1. Coupling functions

In order to enable an efficient coupling between the two domains, a coupling zone is defined and a coupling function is computed which is used to obtain the global energy. First a coupling length $L_c$ is chosen such that will be the characteristic length of the coupling region. Subsequently, the patch of atoms, $\Omega_m$, is included in the continuum at a given position. Atoms at a distance $\ell \le L_c$ from the Molecular Dynamics Box (MD-Box) boundary are considered to be in the coupling zone $\Omega_c$. The finite elements in this zone are named “coupling elements”. Inside the MD-Box, where only the atomistic model applies, elements are removed. The resulting, discretized domain is shown in Figure 1. In principle, each atom must be coupled to the finite element discretization of the underlying continuum. It suffices to couple only a limited number of atoms to the underlying continuum without loss of accuracy. Since this directly affects the size of the coupling matrices that
<figure><img src="image_2.png" /><figcaption>Fig. 1. Discretized domain with the coupling region. The dark, “pear-shaped” area is the domain where an MD calculation is carried out ($\Omega_m$). The coupling region $\Omega_c$ consists of elements that are surrounded by a bold line</figcaption></figure>

will be derived next, such a limited coupling markedly decreases the computational effort, and therefore, a bigger scale transition becomes possible at reasonable computational costs.

In the coupling 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 [28]. For this purpose we define the following functions, see Figure 2:

$$ \alpha : \Omega_M \rightarrow [0, 1] \qquad (12) $$

$$ \beta : \Omega_m \rightarrow [0, 1] \qquad (13) $$

such that:

$$ \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 \\ \\ \alpha(\mathbf{x}) + \beta(\mathbf{x}) = 1 & \text{for } \mathbf{x} \in \Omega_c \end{array} \right. \qquad (14) $$

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, while in $\Omega_m$ we have a discrete field, which is only defined at the geometrical points corresponding to the atoms. Therefore we construct a so-called “mediator space”, denoted by $\mathcal{M}$, on which we project the fields $\dot{\mathbf{u}}$ and $\dot{\mathbf{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 project the velocities using an operator $\Pi$ on a discrete subset $\bar{\Omega}_c$ of the atomic positions included in
<figure><img src="image_3.png" /><figcaption>Fig. 2. Partition of unity for the energy distribution</figcaption></figure>

$\Omega_c$. Considering that $\mathcal{M}$ is built as a Hilbert space, we introduce a scalar product $c$ from $\mathcal{M} \times \mathcal{M}$ onto $\mathbb{R}$. With these definitions we formulate the velocity coupling operator as:

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

with $c$ the classical scalar product on $\mathcal{M}$. 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{U}}^{ad,0} \times \dot{\mathcal{D}}^{ad,0} \times \mathcal{M}, \\ \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:} \end{gathered} $$

$$ a_{\alpha,M}(\mathbf{u}, \mathbf{v}^*) + a_{\beta,m}(\mathbf{d}, \mathbf{w}^*) + c(\boldsymbol{\lambda}, \Pi\mathbf{v}^* - \Pi\mathbf{w}^*) + c(\boldsymbol{\mu}^*, \Pi\dot{\mathbf{u}} - \Pi\dot{\mathbf{d}}) = l_{\alpha,M}(\mathbf{v}^*) \quad (16) $$

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})$, see [37] for details. In the atomistic subdomain the modified form $a_{\beta,m}$ that takes into account the distribution of the energy reads:

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

with:

$$ \beta_i = \beta(\mathbf{d}_i), \quad \mathbf{m}_{\beta} = [\beta_i \delta_{i,j} m_i], \quad \mathbf{f}_{\beta} = [\mathbf{f}_{\beta,i}] \quad (18) $$

## 4.2. Discretized problem

In a manner which is by now standard the interpolation of each component of the displacement field is enriched with discontinuous functions in order to properly capture the traction-free discontinuity in the wake of the crack tip:

$$ \forall \mathbf{x} \in \Omega_M, \quad u_h(\mathbf{x}) = \sum_{i \in \mathcal{N}_M} N_i(\mathbf{x}) \bar{u}_i + \sum_{i \in \mathcal{N}_{cut}} N_i(\mathbf{x}) \mathcal{H}_{\Gamma_d} \hat{u}_i \quad (19) $$

where $N_i$ are standard finite element shape functions supported by the set of nodes $\mathcal{N}_M$ included in the discretized domain $\Omega_M$. Nodes in $\mathcal{N}_{cut}$ have their support cut by the discontinuity. They hold additional degrees of freedom $\hat{u}_i$ corresponding to the discontinuous function $\mathcal{H}_{\Gamma_d}$ defined by:

$$ \mathcal{H}_{\Gamma_d}(\mathbf{x}) = \frac{\mathbf{x} \cdot \mathbf{n}_{\Gamma_d}}{\|\mathbf{x} \cdot \mathbf{n}_{\Gamma_d}\|} \quad (20) $$
with $\mathbf{n}_{\Gamma_d}$ the normal to the discontinuity $\Gamma$. Symbolically, eq. (19) can be written as

$$ \forall \mathbf{x} \in \Omega_M, \quad \mathbf{u}_h = \mathbf{N}^T \mathbf{U} \qquad (21) $$

where the matrix $\mathbf{N}$ contains the standard interpolation polynomials $N_i(\mathbf{x})$ as well as the discontinuous function $\mathcal{H}_{\Gamma_d}$, and the array $\mathbf{U}$ contains the displacement degrees-of-freedom $\bar{u}_i$ and $\hat{u}_i$. The transition within the domain $\Omega_M$ between the subdomain where the nodes are “enriched” and the part which has just the standard formulation does not affect the Molecular Dynamics computation other than through the coupling matrices.

With the latter symbolic notation the bilinear form $a_{\alpha,M}$ and the linear form $l_{\alpha,M}$ 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 (22) $$

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

where the term that represents the body forces has been omitted for simplicity, and

$$ \mathbf{M} = \int_{\Omega_M} \rho \mathbf{N}^T \mathbf{N} d\Omega \qquad (24) $$

$$ \mathbf{K} = \int_{\Omega_M} \nabla \mathbf{N}^T \mathbf{K} \nabla \mathbf{N} d\Omega \qquad (25) $$

the mass and stiffness matrices, respectively. With the standard definition of the scalar product, the coupling term in the continuum can be discretized as follows:

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

with $\mathbf{C}_M$ the continuum coupling matrix. The vector $\mathbf{\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$.

Using the Lagrange multipliers in the atomistic domain similar to that in the continuum domain:

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

the weighted and coupled system (16) 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\mathbf{\Lambda}) + \mathbf{W}^{*T}(\mathbf{m}_{\beta}\ddot{\mathbf{d}} - \mathbf{f}_{\beta} - \mathbf{C}_m\mathbf{\Lambda}) + \mathbf{\mu}^{*T}(\mathbf{C}_M^T\dot{\mathbf{U}} - \mathbf{C}_m^T\dot{\mathbf{d}}) = \mathbf{V}^{*T}\mathbf{F}_{\alpha} \qquad (28) $$

Since this set must hold for any admissible $(\mathbf{V}^*, \mathbf{W}^*, \mathbf{\mu}^*)$ we finally obtain:

$$ \left\{ \begin{array}{l} \mathbf{M}_{\alpha}\ddot{\mathbf{U}} + \mathbf{K}_{\alpha}\mathbf{U} = \mathbf{F}_{\alpha} - \mathbf{C}_M\mathbf{\Lambda} \\ \\ \mathbf{m}_{\beta}\ddot{\mathbf{d}} = \mathbf{f}_{\beta} + \mathbf{C}_m\mathbf{\Lambda} \\ \\ \mathbf{C}_M^T\dot{\mathbf{U}} = \mathbf{C}_m^T\dot{\mathbf{d}} \end{array} \right. \qquad (29) $$

with $(\mathbf{U}, \mathbf{d}, \mathbf{\Lambda})$ the set of unknowns. Details on the time integration scheme associated with this set of coupled ordinary differential equations are given in Ref. [37].
<figure><img src="image_4.png" /><figcaption>Fig. 3. Mechanical energy (drawn line) and the work stored in the Lagrange multipliers (dashed line) when the prediction step is with the full stiffness and mass matrices.</figcaption></figure>

### 4.3. Energy transfer between the atomistics and continuum domains

The energy transfer for the coupling between the continuum and atomistic domains is studied by means of a longitudinal bar discretized with finite elements, and 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. 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 [39]. The computation continues for 2000 time steps with $\Delta t = 1 \times 10^{-15} s$, which amounts to 95% of the critical time step.

Figure 3 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 exits the coupling zone. We observe that the work of the Lagrange multipliers is complementary to the mechanical energy, and the energy balance is therefore satisfied. From energy plots of Figure 4 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.

## 5. DYNAMIC FRACTURE

We now proceed with a two-dimensional simulation of dynamic fracture. A copper single crystal is considered in its (111) plane, so that the two-dimensional lattice is hexagonal. The Lennard-Jones potential is used in the molecular dynamics simulation with parameters from [38]: $a = 0.415$ eV and $b = 0.2277$ nm. The Young's modulus $E$
Mechanical energy (zJ)

<figure><img src="image_5.png" /><figcaption>(a) Energy transfer for $L_c$ consisting of 2 elements</figcaption></figure>

Total energy drift

<figure><img src="image_6.png" /><figcaption>(b) Energy drift for $L_c$ consisting of 2 elements</figcaption></figure>

Mechanical energy (zJ)

<figure><img src="image_7.png" /><figcaption>(c) Energy transfer for $L_c$ consisting of 5 elements</figcaption></figure>

Total energy drift

<figure><img src="image_8.png" /><figcaption>(d) Energy drift for $L_c$ consisting of 5 elements</figcaption></figure>

Mechanical energy (zJ)

<figure><img src="image_9.png" /><figcaption>(e) Energy transfer for $L_c$ consisting of 10 elements</figcaption></figure>

Total energy drift

<figure><img src="image_10.png" /><figcaption>(f) Energy drift for $L_c$ consisting of 10 elements</figcaption></figure>

Fig. 4. 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.
<figure><img src="image_11.png" /><figcaption>Fig. 5. Deformed configuration at $t = 33.209$ ps</figcaption></figure>

and Poisson's ratio $\nu$ for the continuum are obtained following the procedure given in [39]: $E = 79.334$ GPa and $\nu = 0.25$. The copper atomic mass is taken as $m = 0.105520602596 \times 10^{-24}$ kg ($m_{Cu} = \frac{63.546}{N_A}$ g with $N_A$ the Avogadro number: $N_A = 6.02214179 \times 10^{23}$ mol$^{-1}$), which corresponds to a mass density $\rho = 1865.250812586 \times 10^3$ kgm$^{-3}$. In the present study, the temperature has not been taken into account, since the focus is on the coupling of Molecular Dynamics to an (extended) finite element method for crack propagation. When extending the methodology to explicitly include the temperature a “thermal equilibrium” has to be achieved in addition to the mechanical equilibrium. This can for instance be done using the Nose-Hoover thermostat method [40].

The computational domain is 100 nm long and 77.5 nm wide with an initial crack of 10 nm. The finite element mesh consists of 1221 quadrilateral elements and 4868 nodes. The element size is about 10 times the interatomic distance. 36635 atoms are put in the vicinity of the crack tip. The width of the coupling domain is approximately 3 nm and 33% of the atoms in this region hold Lagrange multipliers. For computational reasons, the results that are presented from now on have been obtained by only including the first neighbours in the atomistic interactions. Indeed, the equilibration techniques and the updates are expensive when we take into account many neighbours. Simulations on a smaller scale have indicated that the results are close to those obtained with more neighbours.

The test consists of applying a velocity on the top and bottom edges of the specimen. The prescribed velocity in this example is $V_p = 47.4$ ms$^{-1}$. The time step is 15.811388 fs (10$^{-15}$s). Figure 5 gives the displacement field after the crack has propagated, at $t = 33.209$ ps. Close inspection reveals that, going from the atomistic zone to the continuum domain, no spurious wave reflections occur, which indicates that the coupling algorithm works properly, and is in agreement with results for linear elasticity presented in [37]. We
also see that, even though the prescribed loading corresponds to Mode-I and the atomistic
lattice is a perfect single crystal, the crack path is not that which we would expect to
obtain with classical continuum methods. We observe crack branching, the occurrence of
dislocations, and locally, mixed-mode behaviour.

6. CONCLUDING REMARKS

A numerical approach has been proposed for combining a molecular dynamics
method and a finite element method that exploits the partition-of-unity property of fi-
nite element shape functions (extended finite element method). The aim is to simulate
dynamic fracture in an efficient manner on basis of elementary physical principles. To
this end the zone around the crack tip is modelled using molecular dynamics. Around
this so-called Molecular Dynamics Box a continuum mechanics approach is adopted, with
the finite element method used for discretization. The partition-of-unity property of the
finite element shape functions is exploited to model the crack in the wake of its tip as a
traction-free discontinuity. The coupling between the continuum and molecular dynamics
zones has a zonal character where the energy is partitioned over both models and a weak
velocity coupling is enforced. In this manner, spurious reflections are avoided and energy
is conserved when a wave travels from one zone into another.

REFERENCES

[1] A.A. Griffith, The phenomena of rupture and flow in solids, *Philosophical Transactions of the Royal Society of London A221* (1921) 163–198.

[2] G.R. Irwin, Analysis of stresses and strains near the end of a crack traversing a plate, *Journal of Applied Mechanics* **24** (1957) 361–364.

[3] J.R. Rice, A path independent integral and the approximate analysis of strain concentration by notches and cracks, *Journal of Applied Mechanics* **35** (1968) 379–386.

[4] G.I. Barenblatt, The mathematical theory of equilibrium cracks in brittle fracture, *Advances in Applied Mechanics* **7** (1962) 55–129.

[5] D.S. Dugdale, Yielding of steel sheets containing slits, *Journal of the Mechanics and Physics of Solids* **8** (1960) 100–108.

[6] A. Hillerborg, M. Modeer and P.E. Petersson, Analysis of crack formation and crack growth in concrete by means of fracture mechanics and finite elements, *Cement and Concrete Research* **6** (1976) 773–782.

[7] A. Needleman, A continuum model for void nucleation by inclusion debonding, *Journal of Applied Mechanics* **54** (1987) 525–531.

[8] J.G. Rots, Smeared and discrete representations of localized fracture, *International Journal of Fracture* **51** (1991) 45–59.

[9] X.P. Xu and A. Needleman, Numerical simulations of fast crack growth in brittle solids, *Journal of the Mechanics and Physics of Solids* **42** (1994) 1397–1434.

[10] G.N. Wells and L.J. Sluys, Discontinuous analysis of softening solids under impact loading, *International Journal for Numerical and Analytical Methods in Geomechanics* **25** (2001) 691–709.

[11] G.N. Wells, R. de Borst and L.J. Sluys, A consistent geometrically non-linear approach for delamination, *International Journal for Numerical Methods in Engineering* **54** (2002) 1333–1355.
[12] R. de Borst, Numerical aspects of cohesive zone models, *Engineering Fracture Mechanics* **70** (2003) 1743–1757.

[13] C.A. Duarte, Hamzeh, T.J. Liszka and W.W. Tworzydlo, A generalized finite element method for the simulation of three-dimensional crack propagation, *Computer Methods in Applied Mechanics and Engineering* **bf190** (2001) 2227–2262.

[14] J. Réthoré, A. Gravouil and A. Combescure, An energy conserving scheme for dynamic crack growth with the extended finite element method, *International Journal for Numerical Methods in Engineering* **63** (2005) 631–659.

[15] T. Menouillard, J. Réthoré and A. Combescure, Efficient explicit time stepping for the extended finite element method, *International Journal for Numerical Methods in Engineering* **68** (2006) 911–939.

[16] J.J.C. Remmers, R. de Borst and A. Needleman, A cohesive segments method for the simulation of crack growth, *Computational Mechanics* **31** (2003) 69–77.

[17] J.J.C. Remmers, R. de Borst and A. Needleman, The simulation of dynamic crack propagation using the cohesive segments method, *Journal of the Mechanics and Physics of Solids* **56** (2008) 70–92.

[18] J. Réthoré, R. de Borst and M.A. Abellan, A two-scale approach for fluid flow in fractured porous media, *International Journal for Numerical Methods in Engineering* **71** (2007) 780–800.

[19] J. Réthoré, R. de Borst and M.A. Abellan, A two-scale model for fluid flow in an unsaturated porous medium with cohesive cracks, *Computational Mechanics* **42** (2008) 227–238.

[20] R.G. Parr, S.R. Gadre and L.J. Bartolotti, Local density functional theory of atoms and molecules, *Proceedings of the National Academy of Sciences* **76** (1979) 2522–2526.

[21] R. Car and M. Parrinello, Unified approach for molecular dynamics and density-functional theory, *Physical Review Letters* **55** (1985) 2471–2474.

[22] F.F. Abraham, R. Walkup, H. Gao, M. Duchaineau, T. Diaz De La Rubia and M. Seager, Simulating materials failure by using up to one billion atoms and the world's fastest computer: work-hardening, *Proceedings of the National Academy of Sciences* **99** (2002) 5783–5787.

[23] S.J. Zhou, P.S. Lomdahl, A.F. Voter and B.L. Holian, Three-dimensional fracture via large-scale molecular dynamics, *Engineering Fracture Mechanics* **61** (1998) 173–187.

[24] R. Miller, M. Ortiz, R. Phillips, V. Shenoy and E.B. Tadmor, Quasicontinuum models of fracture and plasticity, *Engineering Fracture Mechanics* **61** (1998) 427–444.

[25] S. Kolhoff, P. Gumbsch and H.F. Frischmeister, Crack propagation in bcc crystals studied with a combined finite element and atomistic model. *Philosophical Magazine A* **64** (1991) 851–878.

[26] L.E. Shilkrot, R.E. Miller and W.A. Curtin, Coupled atomistic and discrete dislocation plasticity, *Physical Review Letters* **89** (2002) 255011–255014.

[27] R. Gracie and T. Belytschko, Concurrently coupled atomistic and XFEM models for dislocations and cracks, *International Journal of Numerical Methods in Engineering* **78** (2009) 354–378.

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

[29] S.P. Xiao and T. Belytschko, A bridging domain method for coupling continua with molecular dynamics, *Computer Methods in Applied Mechanics and Engineering* **193** (2004) 1645–1669.

[30] S.L. Zhang, R. Khare, Q. Lu and T. Belytschko, A bridging domain and strain computation method for coupled atomistic-continuum modelling of solids, *International Journal of Numerical Methods in Engineering* **70** (2007) 913–933.
[31] P.A. Guidault and T. Belytschko, On the L2 and the H1 couplings for an overlapping domain decomposition method using Lagrange multipliers, *International Journal of Numerical Methods in Engineering* **70** (2007) 322–350.

[32] M. Xu and T. Belytschko, Conservation properties of the bridging domain method coupled molecular/continuum dynamics, *International Journal for Numerical Methods in Engineering* **76** (2008) 278–294.

[33] J. Fish and W. Chen, Discrete-to-continuum bridging based on multigrid principles, *Computer Methods in Applied Mechanics and Engineering* **193** (2004) 1693–1711.

[34] C. Farhat, I. Harari and Hetmaniuk, The discontinuous enrichment method for multiscale analysis, *Computer Methods in Applied Mechanics and Engineering* **192** (2003) 3195–3209.

[35] G.J. Wagner and W.K. Liu, Coupling of atomistic and continuum simulations using a bridging scale decomposition, *Journal of Computational Physics* **190** (2003) 249–274

[36] D.E. Farrell, H.S. Park and W.K. Liu, Implementation aspects of the bridging scale method and application to intersonic crack propagation, *International Journal for Numerical Methods in Engineering* **71** (2007) 583–605

[37] P. Aubertin, J. Réthoré and R. de Borst, Energy conservation of atomistic/continuum coupling, *International Journal for Numerical Methods in Engineering* **78** (2009) 1365-1388.

[38] P. Agrawal, B. Rice and D. Thompson, Predicting trends in rate parameters for self-diffusion on FCC metal surfaces, *Surface Science* **515** (2002) 21–35.

[39] P. Aubertin, Coupling of Atomistic and Continuum Models: Dynamic Crack Propagation. *Dissertation*, INSA Lyon (2008).

[40] S. Nose, Constant-temperature molecular dynamics, *Journal of Physics: Condensed Matter* **2** (1990) 115–119.