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

HAL
open science

# Dynamic crack propagation using a combined molecular dynamics/extended finite element approach

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

► To cite this version:

Pascal Aubertin, Julien Réthoré, René de Borst. Dynamic crack propagation using a combined molecular dynamics/extended finite element approach. International Journal for Multiscale Computational Engineering, 2010, 8, pp.221-235. 10.1615/IntJMultCompEng.v8.i2.70. hal-00938617

HAL Id: hal-00938617

https://hal.science/hal-00938617

Submitted on 15 Jul 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_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Dynamic crack propagation using a combined molecular dynamics
/ extended finite element approach

Pascal Aubertin¹, Julien Réthore¹, and René de Borst²,\*

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

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

SUMMARY

A method is presented for simulating dynamic crack propagation using a coupled molecular dynamics / extended finite element method. Molecular dynamics is used at the crack tip while the extended finite element method naturally models the crack in the wake of the tip as a traction-free discontinuity. After recalling the basic molecular dynamics equations, the discretization of the continuum and the traction-free discontinuity via the extended finite element method, and the zonal coupling method between both domains, two-dimensional computations of dynamic fracture are presented, including a discussion how to move and/or expand the zone in which molecular dynamics is used upon crack propagation.

KEY WORDS: multiscale methods, molecular dynamics, extended finite element method, fracture, crack propagation

1. INTRODUCTION

Quantum mechanics is probably the most appropriate theory to describe fracture from
a physics point of view, but the difficulties to relate quantum mechanics to continuum
mechanics, e.g. via Density Functional Theory [1, 2] presently seem insurmountable. One
scale of observation higher is to use molecular dynamics to describe fracture processes from
a fundamental physics point of view. Indeed, researchers have recently used this approach to
describe fracture, e.g. [3, 4, 5]. A disadvantage of the approach is that it is computationally
demanding. For this reason multi-scale approaches have been introduced, in fracture [6], as well
as in plasticity [7], 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 contribution 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

<footnote>*Correspondence to: René de Borst, Department of Mechanical Engineering, Eindhoven University of Technology, P.O. Box 513, NL-5600 MB Eindhoven, Netherlands, R.d.Borst@tue.nl</footnote>
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 [8], but the current paper makes a further advancement in that it includes dynamic crack propagation.

Partition-of-unity based finite element methods [9, 10, 11] have been used before in computational methods for analysing and predicting fracture, since they provide elegant solutions to incorporate stress singularities in domain-based discretization methods. More recently, partition-of-unity finite element methods have also been shown to be very versatile regarding the incorporation of cohesive-zone models for fracture, e.g. [12, 13, 14]. In particular, it naturally enables crack propagation, also in dynamics [15, 16, 17, 18] and in multi-phase continua [19].

This paper is organized as follows. First, we briefly recall the equations of molecular dynamics, the main equations that describe the partition-of-unity based finite element method, and indicate how both domains can be coupled. For a full analysis of the coupling method the reader is referred to [20]. Then, the main body of the paper follows which is concerned with a full two-dimensional coupled analysis of dynamic crack propagation. A major problem is that when the crack tip touches the boundary of the zone around the tip that is modelled using molecular dynamics, crack arrest occurs artificially. For this reason a strategy has been devised to expand and move the zone that is modelled using molecular dynamics during the computation. The numerical results, which are given for different loading rates, shows multiple branching, shear banding and local mixed-mode behaviour, and resemble recent simulations on dynamic fracture using cohesive-zone models [18, 21].

## 2. MULTISCALE MODEL

### 2.1. Molecular Dynamics

The classical Newtonian equations of motion for an atom i with mass $m_i$ and position $\mathbf{r}_i$ are:

$$m_i \ddot{\mathbf{r}}_i = -\frac{\partial \mathcal{U}(\mathbf{r}_1, \dots, \mathbf{r}_N)}{\partial \mathbf{r}_i} = -\nabla_i \mathcal{U} \doteq \mathbf{f}_i \quad (1)$$

The force $\mathbf{f}_i$ derives from the potential energy function $\mathcal{U}$ and can be viewed as the sum of all the interatomic contributions:

$$\mathbf{f}_i = -\frac{\partial \mathcal{U}}{\partial \mathbf{r}_i} = \sum_{j \neq i} \mathbf{f}_{ij} \quad (2)$$

### 2.2. Partition-of-Unity Based Finite Element Method

In the part of the domain that is modelled as a continuum a weak formulation is used as point of departure for the finite element discretization. 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:

$$u = \sum_{i \in N_M} N_i \bar{u}_i + \sum_{i \in N_{cut}} N_i \mathcal{H}_{\Gamma_d} \hat{u}_i \quad (3)$$
<figure><img src="image_3.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>

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_{cut}$ have their support cut by the discontinuity. In addition to the standard degrees of freedom $\bar{u}_i$ they also hold degrees of freedom $\hat{u}_i$ that correspond 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 (4) $$

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

$$ \mathbf{u} = \mathbf{N}\mathbf{U} \qquad (5) $$

where the matrix $\mathbf{N}$ contains the standard interpolation polynomials $N_i$ 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 interpolation can be achieved in various ways, e.g. [22].

### 2.3. Coupling Method

In order to enable an efficient coupling between the two domains, a coupling zone is defined and a coupling function is computed which will be used to obtain the global energy.

First a coupling length $L_c$ must be 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 $r \le L_c$ from the Molecular Dynamics Box (MD-Box) boundary are then considered to be in the coupling zone $\Omega_c$. The finite elements in this zone are named “coupling elements”. Inside the MD-Box, where only the atomistic model applies, elements are removed. The resulting, discretized domain is shown in Figure 1.

In order to write a global formulation for the entire domain, we partition the energy between both models using the functions:

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

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

where

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

The function $\alpha$ has to be computed in the coupling zone, $\beta$ being its complement to 1, subject to the boundary conditions:

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

The global formulation utilizes Lagrange multipliers in order to ensure a velocity coupling in the domain $\Omega_c$, see [20] for details. In principle each atom must be coupled to the finite element discretization of the underlying continuum. However, in practice it suffices to take only a limited number of atoms which have to be coupled via Lagrange multipliers.

### 3. DYNAMIC FRACTURE

#### 3.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 [23, 24, 25], but, as pointed out by Zhou [26], 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:

$$ \sigma_i = \frac{1}{2V_i} \sum_{r_{ij} < r_c} f_{ji} \otimes r_{ij} \qquad (10) $$

where $V_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:

$$ \sigma_i^{\text{avg}} = \frac{1}{N_{n,i}} \sum_{j=1}^{N_{n,i}} \sigma_j \qquad (11) $$
<figure><img src="image_4.png" /><figcaption>Figure 2. Atomistic potential energy variation during propagation. $V_p = 47.4 \text{ ms}^{-1}$</figcaption></figure>
<figure><img src="image_5.png" /><figcaption>Figure 3. Crack tips and branches for an evolving fracture pattern</figcaption></figure>

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}^{avg^2} + \sigma_{yz}^{avg^2} + \sigma_{xz}^{avg^2}) \qquad (12)$$

## 3.2. Moving Molecular Dynamics Box

At a generic stage during the loading process the crack tip will touch the boundary of the MD-Box and, at this point, the crack will be artificially arrested. In order to adapt the MD-Box to the dynamic evolution of the cracks, the crack path has to be identified during the computation. Indeed, even in a perfect crystal with a Mode-I loading, crack branching can develop and the crack path cannot easily be predicted a priori.

Fundamentally, a crack opening corresponds to a dissipation of energy. The basic idea for the evolution of the MD-Box is that when the energy stored at the crack tip exceeds a preset value, propagation can occur. Within the atomistic framework, the crack propagates when bonds are broken, which means that some atoms have displaced sufficiently to be out of the local potential well. Thus, the local atomistic potential energy can be considered to be a good indicator for the identification of the location of cracks. In the intact, undeformed zone, a given atom is at equilibrium, surrounded by its neighbours. If crack propagation occurs in this zone, the atom will move, “lose” some of its neighbours and consequently, its potential energy is affected.

At each time step, the local potential energy $U_i$ is computed for each atom $i$ and is compared
with the initial value $\mathcal{U}_i^0$. We build a discrete field $\mathcal{E}$, which is named the “energy indicator”:

$$ \mathcal{E}_i = 100 \times \frac{\mathcal{U}_i - \mathcal{U}_i^0}{\mathcal{U}_i^0} \qquad (13) $$

and is defined for each atom $i$. From Figures 2(a), 2(c) and 2(e) it is observed that this indicator is zero almost everywhere in the atomistic domain, except at the atoms near the crack. Moreover, the highest value of the energy indicator is located around the crack tip. In Figures 2(b), 2(d) and 2(f), we have put a threshold on $\mathcal{E}_i$:

$$ \left\{ \begin{array}{ll} \bar{\mathcal{E}}_i = 0 & \text{if } \mathcal{E}_i < 10 \\ \bar{\mathcal{E}}_i = 1 & \text{if } \mathcal{E}_i \ge 10 \end{array} \right. \qquad (14) $$

which corresponds to a 10% variation of the potential energy $\mathcal{U}_i$. The discrete field $\bar{\mathcal{E}}$ clearly indicates the crack path in the atomistic domain. One can even distinguish between crack openings and shear bands, depending on the displacement field. Accordingly, this indicator provides useful information: At each time step, the crack tips can be located and their paths can be stored. Considering the “cracked” atoms, i.e. those for which $\bar{\mathcal{E}}_i = 1$, the nearest neighbours are counted and the Von Mises atomistic stress is computed in order to pin-point the crack tips. Indeed, the tips can be viewed as the extremities of the discrete domain of the “cracked” atoms. Using a basic graph theory routine, the global crack path can be viewed as a tree, with a clear identification of extremities and branches. Then, branches are stored and the positions of the crack tips are tracked.

The global crack path consists of the successive positions of all crack tips and the number of cracks is assumed to be equal to the number of branches in the tree. At time step $n$, we then define, for atom $A_i$:

$$ \left\{ \begin{array}{ll} \mathcal{B}^n = \{A_i | \bar{\mathcal{E}}_i \ge 1\} & \Rightarrow \text{branches} \\ \mathcal{T}^n = \{A_i | \bar{\mathcal{E}}_i = 2\} & \Rightarrow \text{tips} \end{array} \right. \qquad (15) $$

Each tip $T_i$ contains the following information, Figure 3:

* The initial position $\mathbf{X}_0^i$, corresponding to the coordinates of the atom that was holding the crack tip at the onset of the corresponding branch;

* The reference position $\mathbf{X}_r^i$, initialized by: $\mathbf{X}_r^i = \mathbf{X}_0^i$, and subsequently updated during propagation;

* The current position $\mathbf{X}_c^i$, which corresponds to the current coordinates of the atom that holds the tip.

During crack evolution, a displacement indicator $\delta_X^i$ is used at each time step $n$:

$$ \forall T_i \in \mathcal{T}^n, \quad \delta_X^i = \mathbf{X}_c^i - \mathbf{X}_r^i. \qquad (16) $$

In order to construct an adapted geometry where the cracks and discontinuities are kept within the MD-Box, an update criterion has been used that compares the displacement indicator to a given value $L_u$:

$$ \forall T_i \in \mathcal{T}^n, \quad \text{if } \|\delta_X^i\| \ge L_u \Rightarrow \text{MD-Box update} \qquad (17) $$

This criterion states that after a prescribed amount of crack growth, the geometry of the atomistic region has to be updated in order to keep the crack tip within the MD-Box. If Eq.
(17) results in an update, $\mathbf{X}_r^i = \mathbf{X}_c^i$ is re-initialized, and a lattice is built centered around $\mathbf{X}_c^i$, in the direction $\mathbf{\delta}_X^i$. At each update of the MD-Box, a new configuration is constructed. Thus, elements have to be removed from the atomistic lattice and the coupling region changes. The weighting functions are computed on this new configuration and the global matrices are re-initialized.

<figure><img src="image_6.png" /><figcaption>Figure 4. Initial configuration of the MD-Box</figcaption></figure>

### 3.3. 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 [27]: $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$ [28]. 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 [29].

The domain of interest has been plotted in Figure 4. It is 100 nm long and 77.5 nm wide with an initial crack. 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>Figure 5. Cracked specimens at the end of the simulations</figcaption></figure>
<figure><img src="image_8.png" /><figcaption>(a) Von Mises stress at the crack tip at $t = 9.219$ ps</figcaption></figure>

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

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

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

Figure 6. Crack initiation for specimen 3: $V_p = 47.4 \text{ ms}^{-1}$

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. 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 removed 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}$.
<figure><img src="image_12.png" /><figcaption>Figure 7. Stress components before crack initiation. Specimen 4: $V_p = 63.2 \text{ ms}^{-1}$</figcaption></figure>

For the lowest loading rate we observe from Figure 5 that there is only one main crack, which propagates in a horizontal sense. For the two next higher loading rates secondary branches develop, which have a shear-band character, and correspond to 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 6 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 7 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 shear bands and cracks. Figure 8 shows the atomic bonds at different time steps for specimen 2. While crack propagation is caused by atomic debonding, shear bands 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. Finally, 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. [30]). Figure 9 shows the same structure, but now the Von Mises stress has been plotted. It is interesting to note that the maximum value of the stress is not located at the crack tip (defined in terms of atomic debonding). This is consistent with the observation that we observe void nucleation ahead of the crack tip. Stress concentrations occurs around defects and ahead of the tip, atoms locally debond, and the propagation subsequently follows this path.

Furthermore, the main crack tip has been tracked during propagation for specimens 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 10(a) – 10(c), the evolution of the crack length has been plotted. For specimen 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}^1 \sim 1967 \text{ ms}^{-1}$. Then, there is an acceleration at $t \sim 22 \text{ ps}$, resulting in a crack propagation velocity $V_{crack,3}^2 \sim 2581 \text{ ms}^{-1}$. Considering the displacement fields from $t \sim 22 \text{ ps}$ onwards, e.g. Figure 10(d)), we observe that main crack is surrounded by two shear bands. 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_{crack,4}^1 \sim 1892 \text{ ms}^{-1}$ around $t = 8 \text{ ps}$,

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

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

From the displacement field in Figure 11) we observe that around these moments, similar changes occur as observed for specimen 3. For $t \sim 15 \text{ ps}$, Figures 11(a) and 11(b), shear bands and dislocations emerge, and the main direction of the crack path changes. Around $t \sim 24 \text{ ps}$, Figures 11(c) and 11(d), a crack branching occurs. New free surfaces appear and the local configuration changes. This can be the explanation for the observed acceleration in the crack propagation.

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_13.png" /><figcaption>Figure 8. Debonding and displacement around the crack tip. Specimen 2.</figcaption></figure>
<figure><img src="image_14.png" /><figcaption>Figure 9. Debonding and Von Mises stress around the crack tip. Specimen 2.</figcaption></figure>
<figure><img src="image_15.png" /><figcaption>(a) Specimen 2: $V_p = 31.6 \text{ ms}^{-1}$</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>(b) Specimen 3: $V_p = 47.6 \text{ ms}^{-1}$</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>(c) Specimen 4: $V_p = 63.2 \text{ ms}^{-1}$</figcaption></figure>

<figure><img src="image_18.png" /><figcaption>(d) Displacement field at $t = 23.721 \text{ ps}$.</figcaption></figure>

Figure 10. Crack tip coordinate as a function of time and displacement field of specimen 3 at $t = 23.721 \text{ ps}$.

We finally examine the evolution of the shear band in specimen 3, Figure 12(a). It arises from the crack at $t = 26.251 \text{ ps}$ and propagates in one of the main lattice directions ($-60°$), see Figure 12(b). A propagation velocity $V_{shear,3} \sim 5050 \text{ ms}^{-1}$ 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.

4. 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
<figure><img src="image_19.png" /><figcaption>Figure 11. Main crack surrounded by shear bands and crack branching. Specimen 4.</figcaption></figure>

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

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 adaptative MD-Box has been constructed, where
a potential energy criterion in the atomistic domain is used to automatically resize the MD-Box
<figure><img src="image_20.png" /><figcaption>(a) Shear band starting from the crack at $t = 31.786$ ps</figcaption></figure>

<figure><img src="image_21.png" /><figcaption>(b) Shear band length as function of time</figcaption></figure>

Figure 12. Shear band evolution for specimen 3: $V_p = 47.4 \text{ ms}^{-1}$
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 local shear banding and local mode mixity can also be identified.

REFERENCES

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

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

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

4. Zhou SJ, Lomdahl PS, Voter AF, Holian BL. Three-dimensional fracture via large-scale molecular dynamics. *Engineering Fracture Mechanics* 1998; **61**: 173-187.

5. Miller R, Ortiz M, Phillips R, Shenoy V, Tadmor EB. Quasicontinuum models of fracture and plasticity. *Engineering Fracture Mechanics* 1998; **61**: 427-444.

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

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

8. Gracie R, Belytschko T. Concurrently coupled atomistic and XFEM models for dislocations and cracks. *International Journal for Numerical Methods in Engineering* DOI: 10.1002/nme.2488.

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

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

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

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

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

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

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

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

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

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

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

20. Aubertin P, Réthoré J, de Borst R. Energy conservation of atomistic/continuum coupling. *International Journal for Numerical Methods in Engineering* DOI: 10.1002/nme.2542.

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

22. Fries, TP. A corrected XFEM approximation without problems in blending elements. *International Journal for Numerical Methods in Engineering* 2008; **75**: 503-532.

23. 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.
24. Lutsko JF. Stress and elastic constants in anisotropic solids: molecular dynamics techniques. *Journal of Applied Physics* 1988; **64**: 1152-1154

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

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

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

28. Aubertin P. Coupling of Atomistic and Continuum Models: Dynamic Crack Propagation. *Dissertation*, INSA Lyon, 2008.

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

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