Pascal Aubertin

Julien Réthoré

René De Borst
email: r.d.borst@tue.nl

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

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

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.

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 [START_REF] Parr | Local density functional theory of atoms and molecules[END_REF][START_REF] Car | Unified approach for molecular dynamics and density-functional theory[END_REF] 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. [START_REF] Abraham | Simulating materials failure by using up to one billion atoms and the world's fastest computer: work-hardening[END_REF][START_REF] Zhou | Three-dimensional fracture via large-scale molecular dynamics[END_REF][START_REF] Miller | Quasicontinuum models of fracture and plasticity[END_REF]. A disadvantage of the approach is that it is computationally demanding. For this reason multi-scale approaches have been introduced, in fracture [START_REF] Kolhoff | Crack propagation in bcc crystals studied with a combined finite element and atomistic model[END_REF], as well as in plasticity [START_REF] Shilkrot | Coupled atomistic and discrete dislocation plasticity[END_REF], 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 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 [START_REF] Gracie | Concurrently coupled atomistic and XFEM models for dislocations and cracks[END_REF], but the current paper makes a further advancement in that it includes dynamic crack propagation.

Partition-of-unity based finite element methods [START_REF] Belytschko | Elastic crack growth in finite elements with minimal remeshing[END_REF][START_REF] Moës | A finite element method for crack growth without remeshing[END_REF][START_REF] Réthoré | A two-scale approach for fluid flow in fractured porous media[END_REF] 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. [START_REF] Wells | Discontinuous analysis of softening solids under impact loading[END_REF][START_REF] Wells | A consistent geometrically non-linear approach for delamination[END_REF][START_REF] De Borst | Numerical aspects of cohesive zone models[END_REF]. In particular, it naturally enables crack propagation, also in dynamics [START_REF] Réthoré | An energy conserving scheme for dynamic crack growth with the extended finite element method[END_REF][START_REF] Menouillard | Efficient explicit time stepping for the extended finite element method[END_REF][START_REF] Remmers | A cohesive segments method for the simulation of crack growth[END_REF][START_REF] Remmers | The simulation of dynamic crack propagation using the cohesive segments method[END_REF] and in multi-phase continua [START_REF] Réthoré | A two-scale model for fluid flow in an unsaturated porous medium with cohesive cracks[END_REF].

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 [START_REF] Aubertin | Energy conservation of atomistic/continuum coupling[END_REF]. 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 [START_REF] Remmers | The simulation of dynamic crack propagation using the cohesive segments method[END_REF][START_REF] Xu | Numerical simulations of fast crack growth in brittle solids[END_REF].

MULTISCALE MODEL

Molecular Dynamics

The classical Newtonian equations of motion for an atom i with mass m i and position r i are:

m i ri = - ∂U (r 1 , . . . , r N ) ∂r i = -∇ i U . = f i (1)
The force f i derives from the potential energy function U and can be viewed as the sum of all the interatomic contributions:

f i = - ∂U ∂r i = j =i f ij (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: where N i are standard finite element shape functions supported by the set of nodes N M included in the discretized domain Ω M . Nodes in N cut have their support cut by the discontinuity. In addition to the standard degrees of freedom ūi they also hold degrees of freedom ûi that correspond to the discontinuous function H Γ d defined by:

u = i∈NM N i ūi + i∈Ncut N i H Γ d ûi (3)
H Γ d (x) = x • n Γ d x • n Γ d (4)
with n Γ d the normal to the discontinuity Γ. Symbolically, eq. ( 3) can be written as

u = NU (5)
where the matrix N contains the standard interpolation polynomials N i as well as the discontinuous function H Γ d , and the array U contains the displacement degrees of freedom ūi and ûi . The transition within the domain Ω 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. [START_REF] Fries | A corrected XFEM approximation without problems in blending elements[END_REF].

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, Ω m , is included in the continuum at a given position. Atoms at a distance r ≤ L c from the Molecular Dynamics Box (MD-Box) boundary are then considered to be in the coupling zone Ω 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:

α : Ω M → [0, 1] (6)
β : Ω m → [0, 1] (7)
where

   α(x) = 1 for x ∈ Ω M \Ω c β(x) = 1 for x ∈ Ω m \Ω c α(x) + β(x) = 1 for x ∈ Ω c (8)
The function α has to be computed in the coupling zone, β being its complement to 1, subject to the boundary conditions:

α(x) = 1 for x ∈ Ω M \Ω c α(x) = 0 for x ∈ Ω m \Ω c (9)
The global formulation utilizes Lagrange multipliers in order to ensure a velocity coupling in the domain Ω c , see [START_REF] Aubertin | Energy conservation of atomistic/continuum coupling[END_REF] 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.

DYNAMIC FRACTURE

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 [START_REF] Basinski | Influence of shear stress on screw dislocations in a model sodium lattice[END_REF][START_REF] Lutsko | Stress and elastic constants in anisotropic solids: molecular dynamics techniques[END_REF][START_REF] Cheung | Atomic-level stress in an inhomogeneous system[END_REF], but, as pointed out by Zhou [START_REF] Zhou | A new look at the atomic level virial stress: on continuum-molecular system equivalence[END_REF], 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:

σ σ σ i = 1 2V i rij <rc f ji ⊗ r ij ( 10
)
where V i is the volume of the atom i, r ij = r ir j , and r ij = |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:

σ σ σ i avg = 1 N n,i Nn,i j=1 σ σ σ j (11)

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 U 0 i . We build a discrete field E, which is named the "energy indicator":

E i = 100 × U i -U 0 i U 0 i ( 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 E i :

E i = 0 if E i < 10 E i = 1 if E i ≥ 10 ( 14
)
which corresponds to a 10% variation of the potential energy U i . The discrete field 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 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 :

B n = {A i | E i ≥ 1} ⇒ branches T n = {A i | E i = 2} ⇒ tips (15)
Each tip T i contains the following information, Figure 3:

• The initial position X i 0 , corresponding to the coordinates of the atom that was holding the crack tip at the onset of the corresponding branch; • The reference position X i r , initialized by: X i r = X i 0 , and subsequently updated during propagation; • The current position X i c , which corresponds to the current coordinates of the atom that holds the tip.

During crack evolution, a displacement indicator δ δ δ i X is used at each time step n:

∀T i ∈ T n , δ δ δ i X = X i c -X i r . (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 :

∀T i ∈ T n , if δ δ δ i X ≥ L u ⇒ MD-Box update (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, X i r = X i c is re-initialized, and a lattice is built centered around X i c , in the direction δ δ δ i X . 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.

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 [START_REF] Agrawal | Predicting trends in rate parameters for self-diffusion on FCC metal surfaces[END_REF]: a = 0.415 eV and b = 0.2277 nm. The Young's modulus E and Poisson's ratio ν for the continuum then become: E = 79.334 GP a and ν = 0.25 [START_REF] Aubertin | Coupling of Atomistic and Continuum Models: Dynamic Crack Propagation[END_REF]. The copper atomic mass is taken as m = 0.105520602596× 10 -24 kg (m Cu = 63.546 NA g with N A the Avogadro number: N A = 6.02214179 × 10 23 mol -1 ), which corresponds to a mass density ρ = 1865.250812586 × 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 [START_REF] Nose | Constant-temperature molecular dynamics[END_REF].

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 approximatively 3 nm and 33% of the atoms in this region hold Lagrange multipliers. For computational reasons, the results that are presented from now on have been obtained by only including the first neighbours in the atomistic interactions. Indeed, the equilibration techniques and the updates are expensive when we take into account many neighbours. 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 ms -1 ; • Specimen 2: V p = 31.6 ms -1 ; • Specimen 3: V p = 47.4 ms -1 ; • Specimen 4: V p = 63.2 ms -1 ; • Specimen 5: V p = 126.5 ms -1 . 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 horizonally, 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 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 σ xx , σ xy and σ yy components also corresponds to those predicted by classical fracture mechanics and the highest value is that of σ 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. [START_REF] Coker | Dynamic crack growth along a polymer composite-homalite interface[END_REF]). 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 ∼ 2815 ms -1 . For specimen 3 one observes different regimes. First, the velocity equals V 1 crack,3 ∼ 1967 ms -1 . Then, there is an acceleration at t ∼ 22 ps, resulting in a crack propagation velocity V 2 crack,3 ∼ 2581 ms -1 . Considering the displacement fields from t ∼ 22 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 1
crack,4 ∼ 1892 ms -1 around t = 8 ps, • First acceleration: V 2 crack,4 ∼ 2447 ms -1 around t = 15 ps, • Second acceleration: V 3 crack,4 ∼ 2841 ms -1 around t = 24 ps. From the displacement field in Figure 11) we observe that around these moments, similar changes occur as observed for specimen 3. For t ∼ 15 ps, Figures 11(a) and 11(b), shear bands and dislocations emerge, and the main direction of the crack path changes. Around t ∼ 24 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 = λ+2µ ρ = 7144 ms -1 ,
• Shear wave speed:

c S = µ ρ = 4124 ms -1 , • Rayleigh wave speed: c R ∼ 3785 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 initiatial velocity, i.e. before it accelerates, is around 50% of the Rayleigh wave speed. We finally examine the evolution of the shear band in specimen 3, Figure 12(a). It arises from the crack at t = 26.251 ps and propagates in one of the main lattice directions (-60 • ), see Figure 12(b). A propagation velocity V shear,3 ∼ 5050 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 a interaction between propagations, accelerations and local bifurcations.

!" # $ %!% % &' !""# " ! % ##$ !%%%# &' !""" (b) Specimen 3: Vp = 47.6 ms -1 !" # ! $% !&& !"& # ! && $% !&&' ! ' # ! $% !&&" ' (

CONCLUDING REMARKS

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

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

Figure 1 .

 1 Figure 1. Discretized domain with the coupling region. The dark, "pear-shaped" area is the domain where an MD calculation is carried out (Ωm). The coupling region Ωc consists of elements that are surrounded by a bold line

( a )

 a Energy indicator at t = 9.488 ps. (b) Threshold E at t = 9.488 ps (c) Energy indicator at t = 13.442 ps (d) Threshold E at t = 13.442 ps (e) Energy indicator at t = 17.395 ps. (f) Threshold E at t = 17.395 ps

Figure 2 . 1 Figure 3 . 2 +

 2132 Figure 2. Atomistic potential energy variation during propagation. Vp = 47.4 ms -1

Figure 4 .

 4 Figure 4. Initial configuration of the MD-Box

  (a) Cracked specimen 1 (Vp = 3.16 ms -1 ) at t = 140.744 ps. (b) Cracked specimen 2 (Vp = 31.6 ms -1 ) at t = 39.535 ps (c) Cracked specimen 3 (Vp = 47.4 ms -1 ) at t = 32.102 ps (d) Cracked specimen 4 (Vp = 63.2 ms -1 ) at t = 29.256 ps (e) Cracked specimen 5 (Vp = 126.5 ms -1 ) at t = 16.605 ps

Figure 5 .Figure 6 .

 56 Figure 5. Cracked specimens at the end of the simulations

  (a) σxx at the crack tip at t = 7.986 ps (b) σxy at the crack tip at t = 7.986 ps (c) σyy at the crack tip at t = 7.986 ps.

Figure 7 .

 7 Figure 7. Stress components before crack initiation. Specimen 4: Vp = 63.2 ms -1

  (a) Displacement field at t = 19.767 ps. (b) Displacement field at t = 20.558 ps. (c) Displacement field at t = 21.349 ps. (d) Displacement field at t = 22.139 ps. (e) Displacement field at t = 22.930 ps. (f) Displacement field at t = 23.721 ps.

Figure 8 .

 8 Figure 8. Debonding and displacement around the crack tip. Specimen 2.

Figure 9 . 2 :

 92 Figure 9. Debonding and Von Mises stress around the crack tip. Specimen 2.

c) Specimen 4 :

 4 Vp = 63.2 ms -1 (d) Displacement field at t = 23.721 ps.

Figure 10 .

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

  (a) Displacement field at t = 15.814 ps. (b) Displacement field at t = 18.977 ps. (c) Displacement field at t = 24.512 ps. (d) Displacement field at t = 26.884 ps.

Figure 11 .

 11 Figure 11. Main crack surrounded by shear bands and crack branching. Specimen 4.

  (a) Shear band starting from the crack at t = 31.786 ps Shear band length as function of time

Figure 12 .

 12 Figure 12. Shear band evolution for specimen 3: Vp = 47.4 ms -1