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

A multiscale method is presented which couples a molecular dynamics approach for describing fracture at the crack tip with an extended ﬁnite 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 ﬁnite element shape functions is used, since in this fashion the crack in the wake of its tip is naturally modelled as a traction-free discontinuity. Next, the zonal coupling method between the atomistic and continuum models is recapitulated. Finally, it is discussed how the stress has been computed in the atomic subdomain, and a two-dimensional computation is presented of dynamic fracture using the coupled model. The result shows multiple branching, which is reminiscent of recent results from simulations on dynamic fracture using cohesive-zone models.


INTRODUCTION
Modern research into fracture commences with the seminal work of Griffith [1].Later, Irwin [2] and Rice [3] established the relation between the stress intensity factors and the energy release rate, and gave linear elastic fracture mechanics a firm basis.However, linear elastic fracture mechanics only applies to crack-like flaws in an otherwise linear elastic solid and when the singularity associated with that flaw is characterized by a non-vanishing energy release rate.The fracture and any dissipative processes must also be confined to a small region in the vicinity of the crack tip.
1 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 quarterpoint elements [4,5], and more recently, the advent of meshless methods [6,7] and partition-ofunity 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 non-linearity is confined to a surface emanating from a classical crack tip, i.e. one with a non-vanishing energy release rate, cohesive zone models as introduced by Barenblatt [12] and Dugdale [13] apply.The cohesive zone approach was extended by Hillerborg et al. [14] and Needleman [15] to circumstances where: (i) an initial crack-like flaw need not be present or, if one is present, it need not be associated with a non-vanishing energy release rate; and (ii) non-linear deformation behaviour may occur over an extended volume.Initially, cohesive-zone models were incorporated in finite element methods via special-purpose interface elements [16,17], but more recently, partition-of-unity finite element methods have shown to be very amenable to the incorporation of cohesive-zone models, e.g.[18][19][20].In particular, it naturally enables 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 [27,28] presently seem insurmountable.One scale of observation higher is to use Molecular Dynamics to describe fracture processes from a more fundamental physics point of view.Indeed, researchers have recently used this approach to describe fracture, e.g.[29][30][31].A disadvantage of the approach is that it is computationally demanding.For this reason multi-scale approaches have been introduced, in fracture [32], as well as in plasticity [33], in which only a part of the body is analysed using molecular dynamics, while the remaining part of the body is modelled using continuum mechanics and discretized using a finite element method.This manuscript furthers along this line and combines molecular dynamics for modelling the fracture process at the crack tip with an extended finite element method (XFEM), where the partition-of-unity property of the polynomial shape functions is exploited to model the crack in the wake of the tip as a traction-free discontinuity.It is noted that recently another approach has been published that couples atomistics and extended finite elements [34], but the current paper makes a further advancement in that it includes dynamic loadings.
A major issue in multi-scale approaches as discussed above is the accurate coupling of both domains, especially when different descriptions are assumed on either domain.While the coupling can, in principle, either be achieved at a discrete interface, or on a zone of a finite size (overlap or zonal coupling), it is believed that zonal approaches, which include the Arlequin method [35,36], the bridging domain method [37][38][39][40][41], discrete-to-continuum bridging [42], the discontinuous enrichment method [43], and bridging scale decomposition [44,45] enable a more gradual transition from one domain to the other.The ability of a gradual transition is especially important for highly dissimilar domains and when wave propagation phenomena are considered, where preservation of the energy and avoiding spurious reflections when a wave exits one domain and enters the other can become an issue.Inspired by earlier work by Ben Dhia and Rateau [35] and Xiao and Belytschko [37] we have chosen a weak coupling between the models in the two adjacent domains.
This paper is organized as follows.First, we briefly list the equations of molecular dynamics and ways to ensure equilibrium of the atomistic domain before starting the computation that involves dynamic propagation of an existing, starter crack.This is followed by a succinct recapitulating of the governing equations of continuum mechanics, both in the strong and the weak forms.The discretization of the continuum subdomain is carried out using the extended finite element method, where the partition-of-unity property of finite element shape functions is used to model the traction-free discontinuity in the wake of the crack tip.Subsequently, it is recapitulated how both domains can be coupled, see also [46] which also presents a detailed analysis of the energy conservation properties of the coupling scheme.The paper concludes with a full two-dimensional coupled analysis of dynamic crack propagation which shows multiple branching and suggests the formation of dislocations, which shows similarities with recent simulations on dynamic fracture using cohesive-zone models [17,25].

Governing equations in the atomic domain
For the discrete domain, i.e. m , we build a grid of N a 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: with m i the mass of atom i and: from where it transpires that the interatomic forces are derived from a potential energy U. d and f assemble the discrete displacements d i and forces f i of the individual atoms, respectively.The internal energy of the discrete domain can be viewed as the sum of each atomic contribution U j : and the force f i acting on atom i can be written as the sum of all elementary forces: 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 : where r ij is the interatomic distance.For details the reader is referred to the Appendix.
In the subdomain m the weak formulation becomes: with w * the test function, and

Equilibrium of the atomic domain
The atomistic domain has to be embedded in a continuum and we cannot consider periodic conditions.Thus, the zone of atoms has an internal boundary with the continuum zone, which has no physical character.Atoms that are on the boundary between the coupling domain and the continuum domain do not have neighbouring atoms in the continuum.This can produce nonphysical forces that break the equilibrium.In order to obtain a global stable system that is in equilibrium before loading the specimen, we therefore have to introduce so-called 'ghost' atoms in the continuum, which act as neighbours for the atoms on the boundaries.Their positions are directly derived from the continuum domain, and their velocities and accelerations do not have to be computed.Their creation merely aims to satisfy the equilibrium of the atomic domain.
For crack propagation we introduce a pre-crack or starter crack in the atomic domain.For this purpose we remove the atoms along a line that defines the position of starter crack.This initial crack crosses the continuum as well as the atomistic zone.The removal of atoms causes that the initial lattice is no longer at equilibrium and there is an additional boundary to consider.However, this boundary has a physical character because it stands for a crack surface.
As alluded to above, equilibrium in the atomistic zone has to be found prior to applying any load to the specimen.A widely used technique is the use of dynamic relaxation that stabilizes a given configuration with adaptative damping.As shown in [47] the use of dynamic relaxation in molecular simulations can be effective in obtaining an equilibrium state.The solution of the problem is viewed as the steady-state solution of a damped wave equation.The classical equations of motion are augmented by a damping term to give: 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.

CONTINUUM MODEL
Assuming small deformation gradients for simplicity, the governing equations in the continuum subdomain M can be written in a standard manner as: For x ∈ M (t) and t ∈[0; T ], given the initial conditions (u(x, 0), u(x, 0)) find (u, r) ∈ U ad ×S ad such that: ü = div r+g d (8) with the mass density and u the continuum displacement vector, r the stress tensor, and g d the body force vector applied in M , subject to the boundary conditions where K is the fourth-order stiffness tensor, n the outward normal vector to * 2 ,a n du d and F d the prescribed displacements and tractions at * 1 and * 2 , respectively.
To allow for a discretization of the continuum subdomain we next specify the weak formulation: ∀v * ∈ Uad,0 , given the initial conditions (u(x, 0), u(x, 0)) find u ∈ U ad such that: with v * the test function, and 4. COUPLING METHOD

Coupling functions
In order to enable an efficient coupling between the two domains, a coupling zone is defined and a coupling function is computed which is used to obtain the global energy.First a coupling length L c is chosen such that will be the characteristic length of the coupling region.Subsequently, the patch of atoms, m , is included in the continuum at a given position.Atoms at a distance ℓ L c from the Molecular Dynamics Box (MD-Box) boundary are considered to be in the coupling where only the atomistic model applies, elements are removed.The resulting, discretized domain is shown in Figure 1.In principle, each atom must be coupled to the finite element discretization of the underlying continuum.It suffices to couple only a limited number of atoms to the underlying continuum without loss of accuracy.Since this directly affects the size of the coupling matrices that will be derived next, such a limited coupling markedly decreases the computational effort, and therefore, a bigger scale transition becomes possible at reasonable computational costs.
In the coupling zone c , we enforce a velocity coupling in a weak sense, and we distribute the energy between both models via a partition of unity [35].For this purpose we define the following functions, see Figure 2: such that: The displacement and velocity fields in the domains M and m have a different nature.In M we have a continuum field, while in 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 M, on which we project the fields u and ḋ in order to be able to compare them.The nature of 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, M has to be a subspace of the physical atomistic space.More precisely, we project the velocities using an operator on a discrete subset c of the atomic positions included in c .Considering that M is built as a Hilbert space, we introduce a scalar product c from M×M onto R. With these definitions we formulate the velocity coupling operator as: with c the classical scalar product on M. The global equations are coupled via Lagrange multipliers and can subsequently be written as: ∀(v * , w * , l * ) ∈ Uad,0 × Ḋad,0 ×M, given the initial conditions (u(x, 0), u(x, 0), d(0), ḋ(0)) find (u, d, k) ∈ U ad ×D ad ×M such that: The modified forms a ,M , a ,m and l ,M take into account the weighting functions (x) and (x),see [46] for details.In the atomic subdomain the modified form a ,m that takes into account the distribution of the energy reads:

Discretized problem
In a manner which is by now standard the interpolation of each component of the displacement field is enriched with discontinuous functions in order to properly capture the traction-free discontinuity in the wake of the crack tip: 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.They hold additional degrees of freedom ûi corresponding to the discontinuous function H d defined by: with n d the normal to the discontinuity .Symbolically, Equation ( 20) can be written as where the matrix N contains the standard interpolation polynomials N i (x) 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 formulation does not affect the Molecular Dynamics computation other than through the coupling matrices.
With the latter symbolic notation the bilinear form a ,M and the linear form l ,M become: where the term that represents the body forces has been omitted for simplicity, and 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: with C M the continuum coupling matrix.The vector K contains the Lagrange multipliers and its size equals the c subset cardinal times the dimension of the space considered.F L M 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 c .
Using the Lagrange multipliers in the atomic domain similar to that in the continuum domain: the weighted and coupled system ( 17) can be cast in a matrix-vector format: Since this set must hold for any admissible (V * , W * , l * ) we finally obtain: with (U, d, K) the set of unknowns.Details on the time integration scheme associated with this set of coupled ordinary differential equations are given in Reference [46].

Mechanical quantities in the atomic domain
In order to extract mechanical quantities from the atomistic domain, we adopt a continuum mechanics point of view to derive classical stress quantities.The atomic stress tensor at an atom i is a measurement of the interatomic interactions of the atom with its neighbours.A widely used stress quantity defined on the atomistic domain is the virial stress, which takes into account the interactions and a kinetic energy contribution.Many formulations have been derived from this virial stress [48][49][50], but, as pointed out by Zhou [51], these definitions, even perfectly correct in a statistical and thermodynamical sense, do not correspond to the Cauchy stress or any other mechanical stress.However, it can be shown that the interatomic interactions part of the virial stress reduces to the Cauchy stress with a physical meaning.We therefore adopt this definition for the stress tensor: where V i is the volume of the atom i.Subsequently, the average of this atomistic stress tensor is computed over the volume around i within the cut-off radius.The average atomistic stress at the atom i thus reads: Finally, the Von Mises atomic stress c at an atom i is defined as follows:

Two-dimensional simulation
We now proceed with a two-dimensional simulation.A copper single crystal is considered in its (111) plane, so that the two-dimensional lattice is hexagonal.The Lennard-Jones potential is used in the molecular dynamics simulation with parameters from [52]: a = 0.415 eV and b = 0.2277 nm.The Young's modulus E and Poisson's ratio for the continuum are obtained following the procedure given in [53]: E = 79.334GPa and = 0.25.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 3kg m −3 .
In the present study, the temperature has not been taken into account, since the focus is on the coupling of Molecular Dynamics to an (extended) finite element method for crack propagation.When extending the methodology to explicitly include the temperature a 'thermal equilibrium' has to be achieved in addition to the mechanical equilibrium.This can for instance be done using the Nose-Hoover thermostat method [54].The computational domain has been plotted in Figure 3.It is 100 nm long and 77.5nm wide with an initial crack of 10 nm.A large MD-Box is considered in order to properly trace the crack propagation during dynamic loading.The finite element mesh consists of 1221 quadrilateral elements and 4868 nodes.The element size is about 10 times the interatomic distance.36635 atoms are put in the vicinity of the crack tip.The width of the coupling domain is approximatively 3 nm and 33% of the atoms in this region hold Lagrange multipliers.For computational reasons, the results that are presented from now on have been obtained by only including the first neighbours in the atomistic interactions.Indeed, the equilibration techniques and the updates are expensive when we take into account many neighbours.Simulations on a smaller scale have indicated that the results are close to those obtained with more neighbours.
The test consists of applying a velocity on the top and bottom edges of the specimen.The prescribed velocity in this example is V p = 47.4ms−1 .The time step is 15.811388 fs (10 −15 s). Figure 4 gives the displacement field after the crack has propagated, at t = 33.209ps.Close inspection reveals that, going from the atomistic zone to the continuum domain, no spurious wave reflections occur, which indicates that the coupling algorithm works properly, and is in agreement with results for linear elasticity presented in [46].We also see that, even though the prescribed loading corresponds to Mode-I and the atomistic lattice is a perfect single crystal, the crack path is not that which we would expect to obtain with classical continuum methods.We observe crack branching, the occurrence of dislocations, and locally, mixed-mode behaviour.
Figures 5 and 6 show the evolution of the displacements and the Von Mises stress in the MD-Box during propagation.We recall that the stress representation is used to represent the behaviour of the interatomic forces qualitatively, and should not be considered as a quantitative result.The concept of stress at the atomistic level is discutable, and, moreover, as we simulate a two-dimensional lattice, the computed values are not representative of a real three-dimensional copper crystal.Yet, when considering Figure 6, the global behaviour of the stress field is consistent with classical results of dynamic crack propagation: A stress concentration forms, which is initially confined at the crack tip, and subsequently propagation occurs and stress waves are emitted from the tips.Moreover, the crack branching that is observed has been demonstrated also in experiments on fast crack propagation [55], and in recent cohesive-zone simulations [17,25].

CONCLUDING REMARKS
A numerical approach has been proposed for combining a molecular dynamics method and a finite element method that exploits the partition-of-unity property of finite element shape functions (extended finite element method).The aim is to simulate dynamic fracture in an efficient manner on basis of elementary physical principles.To this end the zone around the crack tip is modelled using molecular dynamics.Around this so-called Molecular Dynamics Box a continuum mechanics approach is adopted, with the finite element method used for discretization.The partition-of-unity property of the finite element shape functions is exploited to model the crack in the wake of its tip as a traction-free discontinuity.The coupling between the continuum and molecular dynamics zones has a zonal character where the energy is partitioned over both models and a weak velocity coupling is enforced.In this manner, spurious reflections are avoided and energy is conserved when a wave travels from one zone into another [46].
Two-dimensional simulations show realistic crack patterns, including branching and local mixedmode behaviour.Also the stress concentrations at the crack tips, and the stress wave propagation show qualitatively realistic patterns.Quantitatively, less can be said as the concept of stress as it has been defined can be questioned at an atomistic level and as three-dimensional effects cannot be represented well in this two-dimensional computation.Furthermore, the computation is limited by the size of the MD-Box.When a crack tip reaches the boundary of the atomistic domain and intersects the coupling zone, the cracks are arrested and cannot propagate further, cf.The neighbours of each atom i are identified and stored in a table neighb i .This table contains N n,i pointers to the neighbouring atoms.A basic method would be to loop over all the atoms to identify a possible neighbour.This would result in a computational time that is O(N a 2 ) for N a atoms.To decrease the computational cost we therefore have implemented a cell method.A grid of regular cells is built in the wake of the atomistic domain, such that each atom i belongs to a cell C I .The cells dimensions are set such that all the neighbours of an atom are included in the neighbours cells of C I , Figure A1.
Such a method improves the computational efficiency regarding the neighbours list, but we have to determine the optimal number of cells.Indeed, if the grid is either too coarse, or too fine, the cell method becomes inefficient.In Figure A2, we have 7500 atoms in the lattice and we have plotted the CPU time needed for the list.In this case the optimal choice is a grid of ∼ 180 cells.Subsequently, the optimal number of cells has been determined for different lattices.This results in Figure A3 (left), in which the optimal number of cells has been plotted as a function of the number of atoms in the lattice.A logarithmic dependence can be identified that will be used in the sequel to build the grid in the computations.The CPU time spent for the computation of the neighbours list now shows a clear improvement, see Figure A3 (right).While the computation time of the basic method is O(N a 2 ), that of the cell method is O(N a √ N a ).For each atom we have built a table neighb i that contains its neighbours.During a computation of dynamic fracture, the neighbourhood of an atom, especially in the vicinity of the crack, will change, and a criterion has to be implemented to update the neighbours list.For this purpose, a Verlet list has been implemented.For each atom, more pointers are stored than just the neighbours.In addition to storing the atoms within the cut-off radius r c , all atoms within a larger radius, r v are stored.However, the interatomic force is computed only for the neighbours within an inner circle, that is for r <r c .Accordingly, if an atom which is originally located such that r c <r <r v , but which for which r <r c at a generic point in the computation, it is considered as an effective neighbour from   that moment onwards, and is included in the computation of the interatomic forces.The criterion used to update the neighbours list is a maximum displacement criterion.At each time step the two maximum atomic displacements d 1,max and d 2,max are considered.If d 1,max +d 2,max r v −r c ,t h e list is updated, as, in this case, an atom can have displaced more than the distance allowed in the neighbours tables.

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

Figure 2 .
Figure 2. Partition of unity for the energy distribution.

Figure A1 .
Figure A1.Atomic lattice and cells: the neighbourhood of an atom is included in the neighbours cells of C I .

Figure A2 .
Figure A2.CPU time needed for neighbours list calculation as a function of the number of cells in the grid.7500 atoms.

Figure
Figure A3.Left: optimal number of cells versus number of atoms in the lattice; Right: CPU time for neighbours list computation versus N a , number of atoms.The dashed line is the basic method and bold line represents the cell method.