# A Multiscale Molecular Dynamics/Extended Finite Element Method for Dynamic Fracture

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

► To cite this version:

Pascal Aubertin, Julien Réthoré, René de Borst. A Multiscale Molecular Dynamics/Extended Finite Element Method for Dynamic Fracture. 18th International Conference on Computer Methods in Mechanics, May 2009, Zielona Góra, Poland. pp.211-237, 10.1007/978-3-642-05241-5_12. hal-04711618

HAL Id: hal-04711618

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

Submitted on 26 May 2025

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 Multiscale Molecular Dynamics / Extended Finite Element Method for Dynamic Fracture

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

**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. It is discussed how the stress has been computed in the atomistic 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.

## 12.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

<footnote>Pascal Aubertin
Université de Lyon, CNRS INSA-Lyon, LaMCoS UMR 5259, France
e-mail: paubertin@gmail.com</footnote>

<footnote>Julien Réthoré
Université de Lyon, CNRS INSA-Lyon, LaMCoS UMR 5259, France
e-mail: Julien.Rethore@insa-lyon.fr</footnote>

<footnote>René de Borst
Eindhoven University of Technology, Department of Mechanical Engineering, P.O. Box 513, 5600
MB Eindhoven, Netherlands
e-mail: R.d.Borst@tue.nl</footnote>
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.

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, 9, 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 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 [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, 19, 20]. In particular, they naturally enable crack propagation, also in dynamics [21, 22, 23, 24, 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 [30, 31] 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. [27, 28, 29]. 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 advance-
ment in that it includes dynamic loadings.

A major issue in multi-scale approaches as discussed above is the accurate cou-
pling of both domains, especially when different descriptions are assumed on either
domain. While the coupling can, in principle, either be achieved at a discrete in-
terface, 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, 38, 39, 40], discrete-to-continuum bridging [42], the discontinuous en-
richment method [41], and bridging scale decomposition [43, 44] 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 molec-
ular dynamics and ways to ensure equilibrium of the atomistic domain before start-
ing the computation that involves dynamic propagation of an existing, starter crack.
This is followed by a succinct recapitulation of the governing equations of con-
tinuum 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 [45] and an analysis is pre-
sented of the energy conservation properties of the coupling scheme. The paper con-
cludes 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].

12.2 Molecular Dynamics

12.2.1 Governing Equations in the Atomistic Domain

For the discrete domain, i.e. Ω<sub>m</sub>, we build a grid of N<sub>a</sub> atoms, and, accordingly, the
initial value problem in this domain can be written as:

For 1 ≤ i ≤ N_a(t) and t ∈ [0; T], given the initial conditions (d(0), ḋ(0)),
find (d, f) ∈ D^{ad} × F^{ad} such that:

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

with $m_i$ the mass of atom $i$ and:
$$ 
\begin{aligned}
\mathcal{D}^{ad} &= \left\{ \mathbf{d} = (\mathbf{d}_i(t))_{1 \le i \le N_a}, \quad \forall t \in [0, T] \right\} \\
\mathcal{F}^{ad} &= \left\{ \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] \right\}
\end{aligned} 
\quad (12.2) $$

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}) \quad (12.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} \quad (12.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} \quad (12.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 \quad (12.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^* $$

### 12.2.2 Equilibrium of the Atomistic 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 boudaries. 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 atomistic domain.

For crack propagation we introduce a pre-crack or starter crack in the atomistic domain. For this purpose we remove the atoms along a line that defines the position of the 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 adaptive damping. As shown in [46] 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:

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

When the steady-state solution is found, the acceleration and the velocity are zero. Therefore, the matrices **M** and **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.

## 12.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{W}^{ad} \times \mathcal{S}^{ad}$ such that:

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

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] \} $$

$$ \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\} \quad (12.9) $$

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 \mathcal{W}^{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:

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

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

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

$$ l_M(\mathbf{v}^*) = \int_{\partial_2 \Omega} \mathbf{F}_d \cdot \mathbf{v}^* \mathrm{d}S + \int_{\Omega_M} \mathbf{g}_d \cdot \mathbf{v}^* \mathrm{d}\Omega \quad (12.12) $$

<figure><img src="image_2.png" /><figcaption>**Fig. 12.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>
## 12.4 Coupling Scheme

### 12.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 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 12.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 12.2:

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

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

such that:

$$ \begin{cases} \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{cases} \qquad (12.15) $$

<figure><img src="image_3.png" /><figcaption>**Fig. 12.2** Partition of unity for the energy distribution</figcaption></figure>

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 $\mathbf{u}$ and $\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 $\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 \mathbf{u} - \Pi \dot{\mathbf{d}}) = \langle \mu^*, \Pi \mathbf{u} - \Pi \dot{\mathbf{d}} \rangle_{\mathcal{M}} \quad (12.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 \mathcal{W}^{ad,0} \times \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{W}^{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 \mathbf{u} - \Pi \dot{\mathbf{d}}) = l_{\alpha,M}(\mathbf{v}^*) \quad (12.17) $$

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 [45] 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 (12.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 (12.19) $$

### 12.4.2 Discretized Problem

In a manner which is 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}_{\bar{f}_d} \hat{u}_i \quad (12.20) $$

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}_{\bar{f}_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 (12.21) $$

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

$$ \forall \mathbf{x} \in \Omega_M, \quad \mathbf{u}_h = \mathbf{N}^T \mathbf{U} \qquad (12.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 (12.23) $$

$$ l_{\alpha,M}(\mathbf{v}_h^*) = \mathbf{V}^{*T}\mathbf{F}_{\alpha} \qquad (12.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} d\Omega \qquad (12.25) $$

$$ \mathbf{K} = \int_{\Omega_M} \nabla \mathbf{N}^T \mathbf{K} \nabla \mathbf{N} d\Omega \qquad (12.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(\boldsymbol{\lambda}, \Pi \mathbf{v}_h^*) = \mathbf{V}^{*T} \mathbf{C}_M \boldsymbol{\Lambda} = \mathbf{V}^{*T} \mathbf{F}_M^L \qquad (12.27) $$

with $\mathbf{C}_M$ the continuum coupling matrix. The vector $\boldsymbol{\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 \boldsymbol{\Lambda} \qquad (12.28) $$

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

Since this set must hold for any admissible $(\mathbf{V}^*, \mathbf{W}^*, \boldsymbol{\mu}^*)$ we finally obtain:
$$ \left\{ \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} \\ \mathbf{C}_{M}^{T} \dot{\mathbf{U}} &= \mathbf{C}_{m}^{T} \dot{\mathbf{d}} \end{aligned} \right. \qquad (12.30) $$

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

### 12.4.3 Time Integration Scheme

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

* Given the quantities at step $n$, compute the displacements $\mathbf{U}_{n+1}$ and $\mathbf{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 $\ddot{\mathbf{U}}_{n+1}$ and $\ddot{\mathbf{d}}_{n+1}$ by taking into account the coupling terms and Lagrange multipliers $\mathbf{\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:

$$ \left\{ \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} \right. \qquad (12.31) $$

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

$$ \left\{ \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} \right. \qquad (12.32) $$

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:

$$ \left\{ \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} \right. \qquad (12.33) $$
* We next adjust the velocities by introducing the coupling terms:

$$ \left\{ \begin{aligned} \dot{\mathbf{U}}_{n+1} &= \mathbf{U}_{n+1}^* - \frac{1}{2} \tilde{\mathbf{M}}_{\alpha}^{-1} \mathbf{F}_{M,n+1}^L \Delta t \\ \dot{\mathbf{d}}_{n+1} &= \mathbf{d}_{n+1}^* + \frac{1}{2} \mathbf{m}_{\beta}^{-1} \mathbf{f}_{m,n+1}^L \Delta t \end{aligned} \right. \qquad (12.34) $$

* Finally we compute the accelerations:

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

The last steps enforce the coupling condition (12.16).

From Eqs (12.16) and (12.30), the coupling condition becomes:

$$ \mathbf{C}_M^T (\dot{\mathbf{U}}_{n+1}^* - \tilde{\mathbf{M}}_{\alpha}^{-1} \mathbf{C}_M \mathbf{\Lambda}_{n+1} \Delta t) = \mathbf{C}_m^T (\dot{\mathbf{d}}_{n+1}^* + \mathbf{m}_{\beta}^{-1} \mathbf{C}_m \mathbf{\Lambda}_{n+1} \Delta t) \qquad (12.36) $$

The new values of the Lagrange multipliers $\mathbf{\Lambda}_{n+1}$ are subsequently computed by solving:

$$ \mathbf{A}\mathbf{\Lambda}_{n+1} = \mathbf{b}_{n+1} \qquad (12.37) $$

with:

$$ \left\{ \begin{aligned} \mathbf{A} &= \frac{\Delta t}{2} (\mathbf{C}_M^T \tilde{\mathbf{M}}_{\alpha}^{-1} \mathbf{C}_M + \mathbf{C}_m^T \mathbf{m}_{\beta}^{-1} \mathbf{C}_m) \\ \mathbf{b}_{n+1} &= \mathbf{C}_M^T \dot{\mathbf{U}}_{n+1}^* - \mathbf{C}_m^T \dot{\mathbf{d}}_{n+1}^* \end{aligned} \right. \qquad (12.38) $$

and $\mathbf{b}_{n+1}$ stands for the weak coupling condition on the predictive velocities. Thus, this term is a measure of the error compared to the solution that satisfies the system (12.30).

<figure><img src="image_4.png" /><figcaption>Fig. 12.3 Wave propagation in a one-dimensional bar where the left and right sides are continuum domains, and the inside zone represents the atomistic region. Displacements at $t = 150 \times 10^{-15}$ s.</figcaption></figure>
<figure><img src="image_5.png" /><figcaption>**Fig. 12.4** Energy transfer for the finite element - Molecular Dynamics coupling in case of prediction step with the weighted matrices. The drawn line is the energy in the left continuum domain. The dashed line is the energy in the 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_6.png" /><figcaption>**Fig. 12.5** 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>

### 12.4.4 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 [52]. The computation continues for 2000 time steps with $\Delta t = 1 \times 10^{-15} s$, which amounts to 95% of the critical time step.

We first analyse the problem using weighted matrices for the predictive part of the solution strategy. We compute the lumped mass matrix, cf. Eq.(12.32), as follows:

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

where $X_I$ is the position of the $I^{th}$ 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$ [40]:

$$ \mathbf{K}_{\alpha,I,J} = \int_{\Omega} \alpha(X) \nabla \mathbf{u}(\mathbf{N}_I(X)) : \mathbf{K} : \nabla \mathbf{u}(\mathbf{N}_J(X)) \, d\Omega \quad (12.40) $$

The displacements shown in Figure 12.3(a) are for a coupling length $L_c$ that includes 5 elements. When the wave passes through the coupling domain, we observe an amplification of the amplitude and a non-negligible reflection. The energy plots of Figure 12.4 show that, even though the total energy is conserved, 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 that is composed of atoms. The atomistic internal forces have been weighted with the weight function, and to solve the problem the continuum internal forces must be constructed such that full symmetry exists between both domains.

The original continuum internal forces were expressed as:

$$ \mathbf{F}_{\alpha} = \mathbf{K}_{\alpha} \mathbf{U} \quad (12.41) $$

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

$$ \mathbf{F}_{\alpha,I} = \alpha(X_I) \mathbf{K} \mathbf{U} \quad (12.42) $$
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:

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

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

$$ \left\{ \begin{aligned} \ddot{\mathbf{U}}_{n+1}^* &= -\tilde{\mathbf{M}}^{-1} \mathbf{K} \mathbf{U}_{n+1} \\ \ddot{\mathbf{d}}_{n+1}^* &= \mathbf{m}^{-1} \mathbf{f}_{n+1} \end{aligned} \right. \qquad (12.44) $$

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 12.3(b). The simulations do not reveal spurious reflections anymore when the wave passes through the coupling zone, and the formation at $t = 0$ is preserved during the computation. Figure 12.5 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. Considering the energy plots of Figure 12.6 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.

# 12.5 Dynamic Fracture

## 12.5.1 Mechanical Quantities in the Atomistic 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 atomistic 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 [47, 48, 49], but, as pointed out by Zhou [50], these definitions, even perfectly correct in a statistical and thermodynamical sense, do not correspond to the Cauchy stress or to 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\mathcal{Y}_i} \sum_{r_{ij} < r_c} \mathbf{f}_{ji} \otimes \mathbf{r}_{ij} \qquad (12.45) $$

where $\mathcal{Y}_i$ is the volume of the atom $i$, $\mathbf{r}_{ij} = \mathbf{r}_i - \mathbf{r}_j$, and $r_{ij} = |\mathbf{r}_{ij}|$. Subsequently, the average of this atomistic stress tensor is computed over the volume around $i$ within the cut-off radius $r_c$. The average atomistic stress at the atom $i$ thus reads: $\boldsymbol{\sigma}_i^{avg} = \frac{1}{N_{n,i}} \sum_{j=1}^{N_{n,i}} \boldsymbol{\sigma}_j$. Finally, the Von Mises atomistic stress $\sigma_c$ at an atom $i$ is defined as follows:

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

### 12.5.2 Examples: Crack Propagation under Velocity Loading

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 [51]: $a = 0.415$ eV and $b = 0.2277$ nm. The Young's modulus $E$ and Poisson's ratio $\nu$ for the continuum then become: $E = 79.334$ GPa and $\nu = 0.25$ [52]. 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 [53].

The domain of interest 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. 8875 atoms are put in the initial MD-Box. The width of the coupling domain is approximately 3 nm and 33% of the atoms in this region hold Lagrange multipliers.
<figure><img src="image_7.png" /><figcaption>(a) Energy transfer for $L_c$ consisting of 2 elements</figcaption></figure>

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

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

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

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

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

**Fig. 12.6** 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_13.png" /><figcaption>(a) Cracked specimen 1 ($V_p = 3.16 \text{ ms}^{-1}$) at $t = 140.744 \text{ ps}$.</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>(b) Cracked specimen 2 ($V_p = 31.6 \text{ ms}^{-1}$) at $t = 39.535 \text{ ps}$.</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>(c) Cracked specimen 3 ($V_p = 47.4 \text{ ms}^{-1}$) at $t = 32.102 \text{ ps}$.</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>(d) Cracked specimen 4 ($V_p = 63.2 \text{ ms}^{-1}$) at $t = 29.256 \text{ ps}$.</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>(e) Cracked specimen 5 ($V_p = 126.5 \text{ ms}^{-1}$) at $t = 16.605 \text{ ps}$.</figcaption></figure>

<caption>**Fig. 12.7** Cracked specimens at the end of the simulations</caption>

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 ac-
count many neighbours. However, simulations on a smaller scale have indicated
that at least qualitatively, the results are rather similar to those obtained with more
neighbours. In order to simulate defects in the lattice, 0.5% of the atoms are re-
moved in a random manner. Example calculations have been carried out for five
different loading rates, applied to the top and bottom edges of the specimen:

* Specimen 1: $V_p = 3.16 \text{ ms}^{-1}$;

• Specimen 2: $V_p = 31.6 \text{ ms}^{-1}$;

• Specimen 3: $V_p = 47.4 \text{ ms}^{-1}$;

• Specimen 4: $V_p = 63.2 \text{ ms}^{-1}$;

• Specimen 5: $V_p = 126.5 \text{ ms}^{-1}$.

For the lowest loading rate we observe from Figure 12.7 that there is only one main
crack, which propagates in a horizontal sense. For the two next higher loading rates
secondary branches develop, which suggest the appearance of evolving dislocations.
For the two highest loading rates multiple crack branching is observed and the main
crack no longer runs horizontally, but shows a tortuous pattern.

Figure 12.8 shows the Von Mises stress just before and after the crack initiation for the specimen that is subjected to the medium loading rate. The right plots in this figure give the stress distributions along the crack – along a horizontal line. The stress concentration clearly localizes at the crack tip before opening. Upon exceeding a threshold value, propagation occurs and the stress decreases rapidly which is accompanied by stress waves that are emitted radially. These results correspond to those obtained in classical fracture mechanics. Indeed, Figure 12.9 shows the stress components at the onset of crack initiation for specimen 4. The shape of the $\sigma_{xx}$, $\sigma_{xy}$ and $\sigma_{yy}$ components also corresponds to those predicted by classical fracture mechanics and the highest value is that of $\sigma_{yy}$, which corresponds to the fact that we have a Mode-I loading in the vertical direction.

Considering the atomistic structure during crack propagation, a distinction can
be made between cracks and dislocations. Figure 12.10 shows the atomic bonds at
different time steps for specimen 2. While crack propagation is caused by atomic
debonding, dislocations are characterized by a local re-arrangement of the lattice,
rather than debonding. Furthermore, one can observe that, even if the main crack
propagates in a horizontal sense, locally its path is governed by lattice defects. Fi-
nally, it is noted that the crack does not propagate in a smooth manner. Instead,
void nucleations are observed ahead of the crack tip. The local process of crack
propagation thus seems a succession of void nucleation, debonding and crack tip
propagation. This phenomenon is known at the atomistic level, and has also been
observed in macroscopic experiments (usually for high velocities, e.g. [54]).

Furthermore, the main crack tip has been tracked during propagation for spec-
imens 2, 3, and 4, in order to plot the curvilinear coordinate of the crack tip as a
function of time. This enables the computation of the crack propagation velocity. In
Figures 11(a) – 11(c), the evolution of the crack length has been plotted. For speci-
men 2 a crack velocity can straightforwardly be identified as $V_{crack,2} \sim 2815 \text{ ms}^{-1}$.
For specimen 3 one observes different regimes. First, the velocity equals $V_{crack,3} \sim$
<figure><img src="image_18.png" /><figcaption>(a) Von Mises stress at the crack tip at $t = 9.219$ ps</figcaption></figure>

<figure><img src="image_19.png" /><figcaption>(b) Stress distribution along the crack line at $t = 9.219$ ps</figcaption></figure>

<figure><img src="image_20.png" /><figcaption>(c) Von Mises stress at the crack tip at $t = 9.805$ ps</figcaption></figure>

<figure><img src="image_21.png" /><figcaption>(d) Stress distribution along the crack line at $t = 9.805$ ps</figcaption></figure>

**Fig. 12.8** Crack initiation for specimen 3: $V_p = 47.4 \text{ ms}^{-1}$

1967 $\text{ms}^{-1}$. Then, there is an acceleration at $t \sim 22$ ps, resulting in a crack propagation velocity $V_{\text{crack},3}^2 \sim 2581 \text{ ms}^{-1}$. Considering the displacement fields from $t \sim 22$ ps onwards, e.g. Figure 11(d)), the main crack seems to be surrounded by two zones of evolving dislocations. This causes wave reflections, which can affect the crack propagation and its velocity. Moreover, as the velocity increases, more dislocations seem to arise. Finally, for specimen 4, three regimes can be identified:

* Initiation: $V_{\text{crack},4}^1 \sim 1892 \text{ ms}^{-1}$ around $t = 8$ ps,

* First acceleration: $V_{\text{crack},4}^2 \sim 2447 \text{ ms}^{-1}$ around $t = 15$ ps,

* Second acceleration: $V_{\text{crack},4}^3 \sim 2841 \text{ ms}^{-1}$ around $t = 24$ ps.

From the displacement field in Figure 12.8 we observe that around these moments, similar changes occur as observed for specimen 3. For $t \sim 15$ ps, Figures 12(a) and 12(b), dislocations emerge, and the main direction of the crack path changes. Around $t \sim 24$ ps, Figures 12(c) and 12(d), a crack branching occurs. New free surfaces appear and the local configuration changes. This can explain the observed acceleration in crack propagation.
<figure><img src="image_22.png" /><figcaption>**Fig. 12.9** Stress components before crack initiation. Specimen 4: $V_p = 63.2 \text{ ms}^{-1}$</figcaption></figure>

From the elastic material properties computed above, the following wave speeds
can be deduced:

* Longitudinal wave speed: $c_L = \sqrt{\frac{\lambda+2\mu}{\rho}} = 7144 \text{ ms}^{-1}$,

* Shear wave speed: $c_S = \sqrt{\frac{\mu}{\rho}} = 4124 \text{ ms}^{-1}$,

* Rayleigh wave speed: $c_R \sim 3785 \text{ ms}^{-1}$.

Comparing these values with the numerically obtained velocities, the final crack
speed appears to be around 75% of the Rayleigh wave speed, and for specimens 3
and 4, the initialial velocity, i.e. before it accelerates, is around 50% of the Rayleigh
wave speed.
<figure><img src="image_23.png" /><figcaption>**Fig. 12.10** Debonding and displacement around the crack tip. Specimen 2.</figcaption></figure>
<figure><img src="image_24.png" /><figcaption>**Fig. 12.11** Crack tip coordinate as a function of time and displacement field of specimen 3 at $t = 23.721$ ps</figcaption></figure>

We finally examine the evolution of the dislocation in specimen 3, Figure 13(a). It arises from the crack at $t = 26.251$ ps and propagates in one of the main lattice directions (−60°), see Figure 13(b). A propagation velocity $V_{shear,3} \sim 5050$ ms⁻¹ can be identified, which, compared with the shear wave velocity $c_S$ seems somewhat high, but it is in the same order of magnitude. The graph also shows that the band propagation is influenced by defects along its path. Indeed, the global velocity is the result of an interaction between propagations, accelerations and local bifurcations.

## 12.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
<figure><img src="image_25.png" /><figcaption>(a) Displacement field at $t = 15.814$ ps.</figcaption></figure>

<figure><img src="image_26.png" /><figcaption>(b) Displacement field at $t = 18.977$ ps.</figcaption></figure>

<figure><img src="image_27.png" /><figcaption>(c) Displacement field at $t = 24.512$ ps.</figcaption></figure>

<figure><img src="image_28.png" /><figcaption>(d) Displacement field at $t = 26.884$ ps.</figcaption></figure>

**Fig. 12.12** Main crack surrounded by crack branching and dislocations. Specimen 4.

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

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 crack is arrested and cannot propagate further. For this reason an adaptive MD-Box has been constructed, the atomistic domain is automatically resized and/or moved when the crack evolves during the computation. In this manner discontinuities always remain inside the atomistic domain. Two-dimensional examples are shown for different loading rates. Depending on the loading rate, different failure modes can be distinguished, ranging from stable propagation of a single crack to multiple branching, while phenomena like evolving dislocations and local mode mixity can also be identified.
<figure><img src="image_29.png" /><figcaption>(a) Evolving dislocation band starting from the crack at $t = 31.786$ ps</figcaption></figure>

<figure><img src="image_30.png" /><figcaption>(b) Evolving dislocation: Length as function of time</figcaption></figure>

**Fig. 12.13** Evolving dislocation band for specimen 3: $V_p = 47.4 \text{ ms}^{-1}$

## References

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

[2] Irwin, G.R.: Analysis of stresses and strains near the end of a crack traversing a plate. Journal of Applied Mechanics 24, 361–364 (1957)
[3] Rice, J.R.: A path independent integral and the approximate analysis of strain concentration by notches and cracks. Journal of Applied Mechanics 35, 379–386 (1968)

[4] Henshell, R.D., Shaw, K.G.: Crack tip finite elements are unnecessary. International Journal for Numerical Methods in Engineering 9, 495–507 (1975)

[5] Barsoum, R.S.: On the use of isoparametric finite elements in linear fracture mechanics. International Journal for Numerical Methods in Engineering 10, 25–37 (1976)

[6] Fleming, M., Chu, Y.A., Moran, B., Belytschko, T.: Enriched element-free Galerkin methods for crack tip fields. International Journal for Numerical Methods in Engineering 40, 1483–1504 (1997)

[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 44, 767–780 (1999)

[8] Belytschko, T., Black, T.: Elastic crack growth in finite elements with minimal remeshing. International Journal for Numerical Methods in Engineering 45, 601–620 (1999)

[9] Moës, N., Dolbow, J., Belytschko, T.: A finite element method for crack growth without remeshing. International Journal for Numerical Methods in Engineering 46, 131–150 (1999)

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

[11] Albuquerque, E.L., Sollero, P., Aliabadi, M.H.: Dual boundary element method for anisotropic dynamic fracture mechanics. International Journal for Numerical Methods in Engineering 59, 1187–1205 (2004)

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

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

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

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

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

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

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

[19] Wells, G.N., de Borst, R., Sluys, L.J.: A consistent geometrically non-linear approach for delamination. International Journal for Numerical Methods in Engineering 54, 1333–1355 (2002)

[20] de Borst, R.: Numerical aspects of cohesive zone models. Engineering Fracture Mechanics 70, 1743–1757 (2003)

[21] Duarte, C.A., Hamzeh, L.T.J., Tworzydło, W.W.: A generalized finite element method for the simulation of three-dimensional crack propagation. Computer Methods in Applied Mechanics and Engineering 190, 2227–2262 (2001)
[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 63, 631–659 (2005)

[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 68, 911–939 (2006)

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

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

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

[27] Abraham, F.F., 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 99, 5783–5787 (2002)

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

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

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

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

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

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

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

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

[36] Prudhomme, S., Ben Dhia, H., Baumann, P.T., Elkhodja, N., Oden, J.T.: Computational analysis of modeling error for the coupling of particle and continuum models by the Arlequin method. Computer Methods in Applied Mechanics and Engineering 197, 3399–3409 (2008)

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

[38] Zhang, S.L., Khare, R., Lu, Q., Belytschko, T.: A bridging domain and strain computation method for coupled atomistic-continuum modelling of solids. International Journal of Numerical Methods in Engineering 70, 913–933 (2007)

[39] Guidault, P.A., Belytschko, T.: On the L2 and the H1 couplings for an overlapping domain decomposition method using Lagrange multipliers. International Journal of Numerical Methods in Engineering 70, 322–350 (2007)
[40] Xu, M., Belytschko, T.: Conservation properties of the bridging domain method for coupled molecular/continuum dynamics. International Journal for Numerical Methods in Engineering 76, 278-294 (2008)

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

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

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

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

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

[46] Pan, L., Metzger, D., Niewczas, M.: The meshless dynamic relaxation technique for simulating atomic structures of materials. ASME Publications PVP 441, 15-26 (2002)

[47] Basinski, Z.S., Duesberry, M.S., Taylor, R.: Influence of shear stress on screw dislocations in a model sodium lattice. Canadian Journal of Physics 49, 2160-2180 (1971)

[48] Lutsko, J.F.: Stress and elastic constants in anisotropic solids: molecular dynamics techniques. Journal of Applied Physics 64, 1152-1154 (1988)

[49] Cheung, K.S., Yip, S.: Atomic-level stress in an inhomogeneous system. Journal of Applied Physics 70, 5688-5690 (1991)

[50] Zhou, M.: A new look at the atomic level virial stress: on continuum-molecular system equivalence. Proceedings of the Royal Society. A. London 459, 2347-2392 (2003)

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

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

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

[54] Coker, D., Rosakis, A.J., Needleman, A.: Dynamic crack growth along a polymer composite-homalite interface. Journal of the Mechanics and Physics of Solids 51, 425-460 (2003)