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

HAL
open science

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

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

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

## ► To cite this version:

Pascal Aubertin, Julien Réthoré, René de Borst. A coupled molecular dynamics and extended finite element method for dynamic crack propagation. International Journal for Numerical Methods in Engineering, 2010, 81, pp.72-88. 10.1002/nme.2675. hal-00938621

HAL Id: hal-00938621

https://hal.science/hal-00938621

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_3.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# A coupled molecular dynamics and 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

²Department of Mechanical Engineering, Eindhoven University of Technology, Eindhoven, Netherlands

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 recapitulated. Finally, it is discussed how the stress has been computed in the atomic subdomain, and a two-dimensional computation is presented of dynamic fracture using the coupled model. The result shows multiple branching, which is reminiscent of recent results from simulations on dynamic fracture using cohesive-zone models.

**KEY WORDS:** 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.

<footnote>*Correspondence to: René de Borst, Department of Mechanical Engineering, Eindhoven University of Technology, P.O. Box 513, NL-5600 MB Eindhoven, Netherlands.</footnote>

<footnote>†E-mail: r.d.borst@tue.nl</footnote>

<footnote>‡Current address: Institute for Computational Engineering and Sciences (ICES), The University of Texas at Austin, 201 East 24th Street, 1 University Station, C0200 Austin, TX 78712-0027, U.S.A.</footnote>
Linear elastic fracture mechanics provides a challenge to standard finite element approaches, since the polynomials that are conventionally applied in finite element methods cannot easily capture the stress singularity at the crack tip which is predicted in linear elastic fracture mechanics. However, methods have been developed to overcome this difficulty, e.g. the so-called quarter-point elements [4, 5], and more recently, the advent of meshless methods [6, 7] and partition-of-unity based finite element methods [8–10] have provided elegant solutions to incorporate stress singularities in domain-based discretization methods. On the other hand, boundary integral methods can naturally incorporate such singularities [11].

When the region in which the separation and dissipative process take place is not small compared to a structural dimension, but any non-linearity 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 [12] and Dugdale [13] apply. The cohesive zone approach was extended by Hillerborg et al. [14] and Needleman [15] 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 incorporated in finite element methods via special-purpose interface elements [16, 17], but more recently, partition-of-unity finite element methods have shown to be very amenable to the incorporation of cohesive-zone models, e.g. [18–20]. In particular, it naturally enables crack propagation, also in dynamics [21–25] and in multi-phase continua [26].

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 [27, 28] 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. [29–31]. A disadvantage of the approach is that it is computationally demanding. For this reason multi-scale approaches have been introduced, in fracture [32], as well as in plasticity [33], 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 [34], 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 [35, 36], the bridging domain method [37–41], discrete-to-continuum bridging [42], the discontinuous enrichment method [43], and bridging scale decomposition [44, 45] 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 [35] and Xiao
and Belytschko [37] 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 and ways to ensure equilibrium of the atomistic domain before starting the computation that involves dynamic propagation of an existing, starter crack. 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 recapitulated how both domains can be coupled, see also [46] which also presents a detailed analysis 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 [17, 25].

## 2. MOLECULAR DYNAMICS

### 2.1. Governing equations in the atomic domain

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: (1)

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

with $m_i$ the mass of atom i and:

$$\mathcal{D}^{ad} = \{\mathbf{d}=(\mathbf{d}_i(t))_{1 \le i \le N_a} \quad \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} \quad \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. For details the reader is referred to the Appendix.

In the subdomain $\Omega_m$ the weak formulation becomes:

$$ \forall \mathbf{w}^* \in \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: (6)

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

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{W}(\mathbf{d}) \cdot \mathbf{w}_i^* $$

## 2.2. Equilibrium of the atomic domain

The atomistic domain has to be embedded in a continuum and we cannot consider periodic conditions. Thus, the zone of atoms has an internal boundary with the continuum zone, which has no physical character. Atoms that are on the boundary between the coupling domain and the continuum domain do not have neighbouring atoms in the continuum. This can produce non-physical forces that break the equilibrium. In order to obtain a global stable system that is in equilibrium before loading the specimen, we therefore have to introduce so-called 'ghost' atoms in the continuum, which act as neighbours for the atoms on the boundaries. Their positions are directly derived from the continuum domain, and their velocities and accelerations do not have to be computed. Their creation merely aims to satisfy the equilibrium of the atomic domain.

For crack propagation we introduce a pre-crack or starter crack in the atomic domain. For this purpose we remove the atoms along a line that defines the position of starter crack. This initial crack crosses the continuum as well as the atomistic zone. The removal of atoms causes that the initial lattice is no longer at equilibrium and there is an additional boundary to consider. However, this boundary has a physical character because it stands for a crack surface.

As alluded to above, equilibrium in the atomistic zone has to be found prior to applying any load to the specimen. A widely used technique is the use of dynamic relaxation that stabilizes a given configuration with adaptative damping. As shown in [47] the use of dynamic relaxation in molecular simulations can be effective in obtaining an equilibrium state. The solution of the problem is viewed as the steady-state solution of a damped wave equation. The classical equations of motion are augmented by a damping term to give:

$$ \mathbf{M}\ddot{\mathbf{x}} + \mathbf{C}\dot{\mathbf{x}} - \nabla_x \mathcal{W}(\mathbf{x}) = 0 \qquad (7) $$

When the steady-state solution is found, the acceleration and the velocity are zero. Therefore, the matrices $\mathbf{M}$ and $\mathbf{C}$ have no physical meaning and can be chosen such that the computation is accelerated in an optimal sense. At each time step of the resolution, the local atomistic stiffness is computed, using the potential energy function, and a critical time step is calculated. Then, either the mass or the damping can be adapted to obtain an optimal convergence. This method has been
used for finding equilibrium in the atomistic zone and gives converged results within an acceptable computational time.

### 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:

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)$)

find ($\mathbf{u}, \boldsymbol{\sigma}) \in \mathcal{U}^{ad} \times \mathcal{S}^{ad}$ such that: (8)

$$ \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} = \{\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 (9) $$

$$ \mathcal{S}^{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.

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

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

find $\mathbf{u} \in \mathcal{U}^{ad}$ such that: (10)

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

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} \boldsymbol{\varepsilon}(\mathbf{u}) : \mathbb{K} : \boldsymbol{\varepsilon}(\mathbf{v}^*) \, d\Omega \quad (11) $$

$$ 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 (12) $$

### 4. COUPLING METHOD

#### 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
<figure><img src="image_4.png" /><figcaption>Figure 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>

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

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

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

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 (15) $$

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

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
<figure><img src="image_5.png" /><figcaption>Figure 2. Partition of unity for the energy distribution.</figcaption></figure>

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 $\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}$ onto $\mathbb{R}$. With these definitions we formulate the velocity coupling operator 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}} \quad (16) $$

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{M}}^{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{M}^{ad} \times \mathcal{D}^{ad} \times \mathcal{M} \text{ such that:} \end{gathered} \quad (17) $$

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

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 [46] for details. In the atomic 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 (18) $$

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 (19) $$

## 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 \bar{\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 (20) $$
where $N_i$ are standard finite element shape functions supported by the set of nodes $N_M$ included in the discretized domain $\Omega_M$. Nodes in $N_{\text{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}\|} \qquad (21) $$

with $\mathbf{n}_{\Gamma_d}$ the normal to the discontinuity $\Gamma$. Symbolically, Equation (20) can be written as

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

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 (23) $$

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

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} \mathrm{d}\Omega \qquad (25) $$

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

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(\lambda, \Pi \mathbf{v}_h^*) = \mathbf{V}^{*T} \mathbf{C}_M \mathbf{\Lambda} = \mathbf{V}^{*T} \mathbf{F}_M^L \qquad (27) $$

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 atomic domain similar to that in the continuum domain:

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

the weighted and coupled system (17) 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 (29) $$

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

$$ \begin{aligned} \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} \qquad (30) \\ \mathbf{C}_M^T\dot{\mathbf{U}} &= \mathbf{C}_m^T\dot{\mathbf{d}} \end{aligned} $$
with (U, d, Λ) the set of unknowns. Details on the time integration scheme associated with this set of coupled ordinary differential equations are given in Reference [46].

## 5. DYNAMIC FRACTURE

### 5.1. Mechanical quantities in the atomic domain

In order to extract mechanical quantities from the atomistic domain, we adopt a continuum mechanics point of view to derive classical stress quantities.

The atomic stress tensor at an atom i is a measurement of the interatomic interactions of the atom with its neighbours. A widely used stress quantity defined on the atomistic domain is the virial stress, which takes into account the interactions and a kinetic energy contribution. Many formulations have been derived from this virial stress [48–50], but, as pointed out by Zhou [51], these definitions, even perfectly correct in a statistical and thermodynamical sense, do not correspond to the Cauchy stress or any other mechanical stress. However, it can be shown that the interatomic interactions part of the virial stress reduces to the Cauchy stress with a physical meaning. We therefore adopt this definition for the stress tensor:

$$ \boldsymbol{\sigma}_i = \frac{1}{2\nu_i} \sum_{r_{ij} < r_c} \mathbf{f}_{ji} \otimes \mathbf{r}_{ij} \qquad (31) $$

where $\nu_i$ is the volume of the atom i. Subsequently, the average of this atomistic stress tensor is computed over the volume around i within the cut-off radius. The average atomistic stress at the atom i thus reads:

$$ \boldsymbol{\sigma}_i^{\text{avg}} = \frac{1}{N_{n,i}} \sum_{j=1}^{N_{n,i}} \boldsymbol{\sigma}_j \qquad (32) $$

Finally, the Von Mises atomic stress $\sigma_c$ at an atom i is defined as follows:

$$ 2(\sigma_{c,i})^2 = (\sigma_{xx}^{\text{avg}} - \sigma_{yy}^{\text{avg}})^2 + (\sigma_{yy}^{\text{avg}} - \sigma_{zz}^{\text{avg}})^2 + (\sigma_{xx}^{\text{avg}} - \sigma_{zz}^{\text{avg}})^2 + 6(\sigma_{xy}^{\text{avg}}^2 + \sigma_{yz}^{\text{avg}}^2 + \sigma_{xz}^{\text{avg}}^2) \qquad (33) $$

### 5.2. Two-dimensional simulation

We now proceed with a two-dimensional simulation. 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 [52]: $a=0.415$ eV and $b=0.2277$ nm. The Young’s modulus $E$ and Poisson’s ratio $\nu$ for the continuum are obtained following the procedure given in [53]: $E=79.334$ GPa and $\nu=0.25$. The copper atomic mass is taken as $m=0.105520602596 \times 10^{-24}$ kg ($m_{Cu}=(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$ kg m$^{-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 [54].
<figure><img src="image_6.png" /><figcaption>Figure 3. Initial configuration.</figcaption></figure>

The computational domain has been plotted in Figure 3. It is 100 nm long and 77.5 nm wide with an initial crack of 10 nm. A large MD-Box is considered in order to properly trace the crack propagation during dynamic loading. 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 approximatively 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 \text{ms}^{-1}$. The time step is 15.811388 fs ($10^{-15} \text{ s}$). Figure 4 gives the displacement field after the crack has propagated, at $t=33.209 \text{ 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 [46]. 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.

Figures 5 and 6 show the evolution of the displacements and the Von Mises stress in the MD-Box during propagation. We recall that the stress representation is used to represent the behaviour of the interatomic forces qualitatively, and should not be considered as a quantitative result. The concept of stress at the atomistic level is discutable, and, moreover, as we simulate a two-dimensional lattice, the computed values are not representative of a real three-dimensional copper crystal. Yet, when considering Figure 6, the global behaviour of the stress field is consistent with classical
<figure><img src="image_7.png" /><figcaption>Figure 4. Deformed configuration at $t=33.209$ ps.</figcaption></figure>

results of dynamic crack propagation: A stress concentration forms, which is initially confined at the crack tip, and subsequently propagation occurs and stress waves are emitted from the tips. Moreover, the crack branching that is observed has been demonstrated also in experiments on fast crack propagation [55], and in recent cohesive-zone simulations [17, 25].

## 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 finite 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 [46].

Two-dimensional simulations show realistic crack patterns, including branching and local mixed-mode behaviour. Also the stress concentrations at the crack tips, and the stress wave propagation show qualitatively realistic patterns. Quantitatively, less can be said as the concept of stress as it has been defined can be questioned at an atomistic level and as three-dimensional effects cannot be represented well in this two-dimensional computation. Furthermore, the computation is limited by the size of the MD-Box. When a crack tip reaches the boundary of the atomistic domain and intersects the coupling zone, the cracks are arrested and cannot propagate further, cf.
<figure><img src="image_8.png" /><figcaption>Figure 5. MD-Box: Displacement field around the crack tip during propagation: (a) displacement field at $t=11.070$ ps; (b) displacement field at $t=17.395$ ps; (c) displacement field at $t=23.721$ ps; (d) displacement field at $t=30.046$ ps; (e) displacement field at $t=36.372$ ps; and (f) displacement field at $t=39.535$ ps.</figcaption></figure>
<figure><img src="image_9.png" /><figcaption>Figure 6. MD-Box: Von Mises stress field around the crack tip during propagation: (a) Von Mises stress at $t=11.070$ ps; (b) Von Mises stress at $t=17.395$ ps; (c) Von Mises stress at $t=23.721$ ps; (d) Von Mises stress at $t=30.046$ ps; (e) Von Mises stress at $t=36.372$ ps; and (f) Von Mises stress at $t=39.535$ ps.</figcaption></figure>
Figures 5(e), (f), 6(e) and (f)). A solution is to let the MD-Box expand and/or move with the evolution of the cracks.

APPENDIX A: COMPUTING A NEIGHBOURS LIST

The neighbours of each atom *i* are identified and stored in a table *neighb*<sub>*i*</sub>. This table contains $N_{n,i}$ pointers to the neighbouring atoms. A basic method would be to loop over all the atoms to identify a possible neighbour. This would result in a computational time that is $O(N_a^2)$ for $N_a$ atoms. To decrease the computational cost we therefore have implemented a cell method. A grid of regular cells is built in the wake of the atomistic domain, such that each atom *i* belongs to a cell $C_I$. The cells dimensions are set such that all the neighbours of an atom are included in the neighbours cells of $C_I$, Figure A1.

Such a method improves the computational efficiency regarding the neighbours list, but we have
to determine the optimal number of cells. Indeed, if the grid is either too coarse, or too fine, the
cell method becomes inefficient. In Figure A2, we have 7500 atoms in the lattice and we have
plotted the CPU time needed for the list. In this case the optimal choice is a grid of ~180 cells.
Subsequently, the optimal number of cells has been determined for different lattices. This results
in Figure A3 (left), in which the optimal number of cells has been plotted as a function of the
number of atoms in the lattice. A logarithmic dependence can be identified that will be used in
the sequel to build the grid in the computations. The CPU time spent for the computation of the
neighbours list now shows a clear improvement, see Figure A3 (right). While the computation
time of the basic method is $O(N_a^2)$, that of the cell method is $O(N_a\sqrt{N_a})$.

For each atom we have built a table *neighb*<sub>*i*</sub> that contains its neighbours. During a computation of dynamic fracture, the neighbourhood of an atom, especially in the vicinity of the crack, will change, and a criterion has to be implemented to update the neighbours list. For this purpose, a Verlet list has been implemented. For each atom, more pointers are stored than just the neighbours. In addition to storing the atoms within the cut-off radius *r*<sub>c</sub>, all atoms within a larger radius, *r*<sub>v</sub> are stored. However, the interatomic force is computed only for the neighbours within an inner circle, that is for *r* (< *r*<sub>c</sub>. Accordingly, if an atom which is originally located such that *r*<sub>c</sub> (< *r*<sub>v</sub>, but which for which *r* < *r*<sub>c</sub> at a generic point in the computation, it is considered as an effective neighbour from

<figure><img src="image_10.png" /><figcaption>Figure A1. Atomic lattice and cells: the neighbourhood of an atom is included in the neighbours cells of $C_I$.</figcaption></figure>
<figure><img src="image_11.png" /><figcaption>Figure A2. CPU time needed for neighbours list calculation as a function of the number of cells in the grid. 7500 atoms.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>Figure A3. Left: optimal number of cells versus number of atoms in the lattice; Right: CPU time for neighbours list computation versus $N_a$, number of atoms. The dashed line is the basic method and bold line represents the cell method.</figcaption></figure>

that moment onwards, and is included in the computation of the interatomic forces. The criterion used to update the neighbours list is a maximum displacement criterion. At each time step the two maximum atomic displacements $d_{1,max}$ and $d_{2,max}$ are considered. If $d_{1,max}+d_{2,max} \ge r_v - r_c$, the list is updated, as, in this case, an atom can have displaced more than the distance allowed in the neighbours tables.

## REFERENCES

1. Griffith AA. The phenomena of rupture and flow in solids. *Philosophical Transactions of the Royal Society of London* 1921; **A221**:163–198.

2. Irwin GR. Analysis of stresses and strains near the end of a crack traversing a plate. *Journal of Applied Mechanics* 1957; **24**:361–364.

3. Rice JR. A path independent integral and the approximate analysis of strain concentration by notches and cracks. *Journal of Applied Mechanics* 1968; **35**:379–386.
4. Henshell RD, Shaw KG. Crack tip finite elements are unnecessary. *International Journal for Numerical Methods in Engineering* 1975; **9**:495-507.

5. Barsoum RS. On the use of isoparametric finite elements in linear fracture mechanics. *International Journal for Numerical Methods in Engineering* 1976; **10**:25-37.

6. Fleming M, Chu YA, Moran B, Belytschko T. Enriched element-free Galerkin methods for crack tip fields. *International Journal for Numerical Methods in Engineering* 1997; **40**:1483-1504.

7. Krysl P, Belytschko T. The element-free Galerkin Method for dynamic propagation of arbitrary 3-D cracks. *International Journal for Numerical Methods in Engineering* 1999; **44**:767-780.

8. Belytschko T, Black T. Elastic crack growth in finite elements with minimal remeshing. *International Journal for Numerical Methods in Engineering* 1999; **45**:601-620.

9. Moës N, Dolbow J, Belytschko T. A finite element method for crack growth without remeshing. *International Journal for Numerical Methods in Engineering* 1999; **46**:131-150.

10. Réthoré J, de Borst R, Abellan MA. A two-scale approach for fluid flow in fractured porous media. *International Journal for Numerical Methods in Engineering* 2007; **71**:780-800.

11. Albuquerque EL, Sollero P, Aliabadi MH. Dual boundary element method for anisotropic dynamic fracture mechanics. *International Journal for Numerical Methods in Engineering* 2004; **59**:1187-1205.

12. Barenblatt GI. The mathematical theory of equilibrium cracks in brittle fracture. *Advances in Applied Mechanics* 1962; **7**:55-129.

13. Dugdale DS. Yielding of steel sheets containing slits. *Journal of the Mechanics and Physics of Solids* 1960; **8**:100-108.

14. Hillerborg A, Modeer M, Petersson PE. Analysis of crack formation and crack growth in concrete by means of fracture mechanics and finite elements. *Cement and Concrete Research* 1976; **6**:773-782.

15. Needleman A. A continuum model for void nucleation by inclusion debonding. *Journal of Applied Mechanics* 1987; **54**:525-531.

16. Rots JG. Smeared and discrete representations of localized fracture. *International Journal of Fracture* 1991; **51**:45-59.

17. Xu XP, Needleman A. Numerical simulations of fast crack growth in brittle solids. *Journal of the Mechanics and Physics of Solids* 1994; **42**:1397-1434.

18. Wells GN, Sluys LJ. Discontinuous analysis of softening solids under impact loading. *International Journal for Numerical and Analytical Methods in Geomechanics* 2001; **25**:691-709.

19. Wells GN, de Borst R, Sluys LJ. A consistent geometrically non-linear approach for delamination. *International Journal for Numerical Methods in Engineering* 2002; **54**:1333-1355.

20. de Borst R. Numerical aspects of cohesive zone models. *Engineering Fracture Mechanics* 2003; **70**:1743-1757.

21. Duarte CA, Hamzeh, Lizska TJ, Tworzydło WW. A generalized finite element method for the simulation of three-dimensional crack propagation. *Computer Methods in Applied Mechanics and Engineering* 2001; **190**: 2227-2262.

22. Réthoré J, Gravouil A, Combescure A. An energy conserving scheme for dynamic crack growth with the extended finite element method. *International Journal for Numerical Methods in Engineering* 2005; **63**:631-659.

23. Menouillard T, Réthoré J, Combescure A. Efficient explicit time stepping for the extended finite element method. *International Journal for Numerical Methods in Engineering* 2006; **68**:911-939.

24. Remmers JJC, de Borst R, Needleman A. A cohesive segments method for the simulation of crack growth. *Computational Mechanics* 2003; **31**:69-77.

25. Remmers JJC, de Borst R, Needleman A. The simulation of dynamic crack propagation using the cohesive segments method. *Journal of the Mechanics and Physics of Solids* 2008; **56**:70-92.

26. Réthoré J, de Borst R, Abellan MA. A two-scale model for fluid flow in an unsaturated porous medium with cohesive cracks. *Computational Mechanics* 2008; **42**:227-238.

27. Parr RG, Gadre SR, Bartolotti LJ. Local density functional theory of atoms and molecules. *Proceedings of the National Academy of Sciences* 1979; **76**:2522-2526.

28. Car R, Parrinello M. Unified approach for molecular dynamics and density-functional theory. *Physical Review Letters* 1985; **55**:2471-2474.

29. Abraham FF, Walkup R, Gao H, Duchaineau M, Diaz De La Rubia T, Seager M. 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* 2002; **99**:5783-5787.

30. Zhou SJ, Lomdahl PS, Voter AF, Holian BL. Three-dimensional fracture via large-scale molecular dynamics. *Engineering Fracture Mechanics* 1998; **61**:173-187.
31. Miller R, Ortiz M, Phillips R, Shenoy V, Tadmor EB. Quasicontinuum models of fracture and plasticity. *Engineering Fracture Mechanics* 1998; **61**:427–444.

32. Kolhoff S, Gumbsch P, Frischmeister HF. Crack propagation in bcc crystals studied with a combined finite element and atomistic model. *Philosophical Magazine A* 1991; **64**:851–878.

33. Shilkrot LE, Miller RE, Curtin WA. Coupled atomistic and discrete dislocation plasticity. *Physical Review Letters* 2002; **89**:255011–255014.

34. Gracie R, Belytschko T. Concurrently coupled atomistic and XFEM models for dislocations and cracks. *International Journal of Numerical Methods in Engineering* 2009; **78**:354–378.

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

36. Prudhomme S, Ben Dhia H, Baumann PT, Elkhodja N, Oden JT. Computational analysis of modeling error for the coupling of particle and continuum models by the Arlequin method. *Computer Methods in Applied Mechanics and Engineering* 2008; **197**:3399–3409.

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

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

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

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

41. Guidault PA, Belytschko T. Bridging domain methods for coupled atomistic-continuum models with L2 or H1 couplings. *International Journal for Numerical Methods in Engineering*. 2009; **77**:1566–1592.

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

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

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

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

46. Aubertin P, Réthoré J, de Borst R. Energy conservation of atomistic/continuum coupling. *International Journal for Numerical Methods in Engineering*. 2009; **78**:1365–1386.

47. Pan L, Metzger D, Niewczas M. *The Meshless Dynamic Relaxation Technique for Simulating Atomic Structures of Materials*, vol. PVP 441. ASME Publications: New York, 2002; 15–26.

48. Basinski ZS, Duesberry MS, Taylor R. Influence of shear stress on screw dislocations in a model sodium lattice. *Canadian Journal of Physics* 1971; **49**:2160–2180.

49. Lutsko JF. Stress and elastic constants in anisotropic solids: molecular dynamics techniques. *Journal of Applied Physics* 1988; **64**:1152–1154.

50. Cheung KS, Yip S. Atomic-level stress in an inhomogeneous system. *Journal of Applied Physics* 1991; **70**: 5688–5690.

51. Zhou M. A new look at the atomic level virial stress: on continuum-molecular system equivalence. *Proceedings of the Royal Society of London, Series A* 2003; **459**:2347–2392.

52. Agrawal P, Rice B, Thompson D. Predicting trends in rate parameters for self-diffusion on FCC metal surfaces. *Surface Science* 2002; **515**:21–35.

53. Aubertin P. Coupling of atomistic and continuum models: dynamic crack propagation. *Dissertation*, INSA Lyon, 2008.

54. Nose S. Constant-temperature molecular dynamics. *Journal of Physics: Condensed Matter* 1990; **2**:115–119.

55. Coker D, Rosakis AJ, Needleman A. Dynamic crack growth along a polymer composite-homalite interface. *Journal of the Mechanics and Physics of Solids* 2003; **51**:425–460.