Concurrent multiscale modeling of three dimensional crack and dislocation propagation

Advances in Engineering Software xxx (2014) xxx–xxx

Contents lists available at ScienceDirect

Advances in Engineering Software journal homepage: www.elsevier.com/locate/advengsoft

Concurrent multiscale modeling of three dimensional crack and dislocation propagation Hossein Talebi a, Mohammad Silani b,c, Timon Rabczuk c,⇑ a

Graduiertenkolleg 1462, Bauhaus-Universität Weimar, 99423 Weimar, Germany Department of Mechanical Engineering, Isfahan University of Technology, Isfahan 84156-83111, Iran c Institute of Structural Mechanics, Bauhaus-Universität Weimar, Marienstr. 15, D-99423 Weimar, Germany b

a r t i c l e

i n f o

Article history: Available online xxxx Keywords: Extended finite element method (XFEM) Molecular dynamics Bridging domain method Multiscale methods LAMMPS Fracture

a b s t r a c t In this manuscript a concurrent coupling scheme is presented to model three dimensional cracks and dislocations at the atomistic level. The scheme couples molecular dynamics to extended finite element method (XFEM) via the Bridging Domain Method (BDM). This method is based on linear weighting of the strain energy over a region (the bridging domain) which conserves the energy in the entire system. To compute the material behavior in the continuum scale, the Cauchy–Born method is used. Many improvements have been made in the implementation to make the method work for the general case of materials and presence of multi-million degrees of freedom. To show the applicability and productivity of the proposed method, two three dimensional crack examples were modeled. The results show that the method and the corresponding implementation are capable of handling dislocation and crack propagation in the three dimensional space. Ó 2014 Elsevier Ltd. All rights reserved.

1. Introduction In the last decade, one major research topic in materials science as well as mechanics was to gain a fundamental understanding of material failure. Earlier approaches were mainly based on the empirical observations. However, material failure depends often on the behavior of the next lower scale. Typical examples are shear bands in metallic materials, micro-cracks in quasi-brittle materials, dislocation and defects in the atomic lattice, etc. In these cases, the response of the body depends on the accuracy of the finer scale modeling. When modeling the fine scale details, the handling of the dislocations are of particular importance. This is because in many metallic systems, the evolution of cracks and other defects are directly influenced by the dislocations and how dislocations interact with each other. There are different methods for modeling cracks and dislocations in materials both in continuum and atomistic scales. Methods which can work on the resolution of the lower scales are too expensive and hardly can be used in real engineering problems. Therefore new methods to simulate material failure are urgently needed in engineering. This was the motivation of a big class of numerical methods called multiscale methods. In these methods, ⇑ Corresponding author at: Chair of Computational Mechanics, Institute of Structural Mechanics, Bauhaus-Universität Weimar, Marienstr. 15, D-99423 Weimar, Germany.

different length and time scales are correlated to keep both accuracy and efficiency. Moreover, whereas the current computational power is capable of solving several billions of degrees of freedom, the multiscale methods are not only capable of solving problems which are impossible to solve with the current computational power; but also they are useful to solve the problems with noticeably reduced sizes to save the cost, energy, computational power and human resources to implement sophisticated parallel software and run huge simulations [1,2]. In the multiscale failure analysis of materials, different length and time scales are involved that cause many complexities to couple these scales. These problems are more severe in case of atomistic–continuum coupling especially due to different treatment of deformation and temperature in these domains. Without suitable coupling methods, non-physical effects such as artificial wave reflections, heating or melting in the atomic system can occur. To date, several concurrent approaches have been developed. Generally, two so-called ‘Interface’ coupling and ‘Handshake’ coupling can be distinguished. The former approach is susceptible to artificial wave reflections on the boundary (however, effective approaches to reduce artificially reflected waves have been developed recently [3]). In ‘Handshake’ couplings approaches, two areas from different scales are gradually transformed into each other. One of the most popular methods in coupling MD to continuum domains is the quasi-continuum method of Tadmor et al. [4]. The coupling can be interpreted as a seamless MD to continuum

http://dx.doi.org/10.1016/j.advengsoft.2014.09.016 0965-9978/Ó 2014 Elsevier Ltd. All rights reserved.

Please cite this article in press as: Talebi H et al. Concurrent multiscale modeling of three dimensional crack and dislocation propagation. Adv Eng Softw (2014), http://dx.doi.org/10.1016/j.advengsoft.2014.09.016

2

H. Talebi et al. / Advances in Engineering Software xxx (2014) xxx–xxx

transition by means of a force-field potential to a continuum potential changing. However, the nodes in the atomic field should be at the same position of the atoms which results in highly distorted FE meshes especially in the ‘coarsening’. Moreover, the propagation of defects (in MD area) requires the conversion of continuum domain to MD field. On top of all, the quasi-continuum method is not capable of modeling the dynamic behavior of systems which is the main focus in this study. The ‘Coupling of Length Scales’ (CLS) method of Abraham et al. [5] and Broughton et al. [6] links the atomic system with the continuum domain which is described by finite elements. At the transition from atomic to the continuum field, atoms are placed at the points of finite elements. A Chebyshev polynomial mapping function was developed by Bayliss and Turkel [7] that can reduce artificial boundary effects. The bridging scale method of Wagner and Liu [8] can eliminate reflection of unauthorized elastic waves on the atom/continuum transition area. In this coupling method the continuum domain is modeled using finite elements or meshfree discretization [9]. Likewise, in the bridge area the nodes should placed in the atom positions. The bridging domain method (BDM) by Xia and Belytschko [10] is based on a L2 link, i.e. displacement compatibility between atomistic and continuum field would accomplished by means of Lagrange multipliers in the weak form. In this paper, the partition of unity enriched methods and the bridging domain method [11] are blended together for crack modeling. Here, a weak coupling between the models in the two adjacent domains is used [12– 14]. It was shown that the spurious wave reflection is minimum in this kind of coupling method [12]. It should be noted that in the XFEM–BDM coupling method, the atomistic region presents only around crack front but not along the whole length of the crack surface. This part of crack is mainly simulated by XFEM in the countinuum domain. The handshake domain allows the crack to be smoothly represented from the continuum to the atomistic domain. The mentioned advantages of XFEM–BDM coupling method over other concurrent multiscale methods such as the standard BDM or the quasi-continuum method of Tadmor et al. [4], Shenoy et al. [15] and Miller and Tadmor [16], makes this method more efficient and applicable for problems involving dynamic fracture. This affords the modeling and simulation of relatively long cracks, by restricting the MD model to a small zone around the crack front. Since the crack front is modeled by molecular dynamics, the creation and propagation of dislocations are easy to capture by this method. This is because, the most straight forward and natural method of modeling dislocations is using MD. Especially near the crack tip, the behavior of dislocations is really complex. There are some methods to handle dislocations using a continuum formulation. For example, Gracie et al. [17] proposed a method to handle dislocation motion by means of XFEM enrichment. However if not impossible, such methods will have many difficulties handling a complex dislocation pattern in a dynamic simulation. Therefore, we will use molecular dynamics to handle dislocation nucleation and propagation as well as crack tip behavior. In comparison to other methods, the bridging domain method has the advantage of natural and robust coupling of the two regions with minimum additional treatments such as artificial damping which makes the method more suitable to apply to different problems. In a nutshell, we are seeking effective methods to model crack and dislocation nucleation and propagation for highly dynamic problems that are also as general as possible, i.e. they are applicable to a wide range of problems. In this path, other advantage of the current study over the previous studies for example by Gracie et al. [18] is that, this study implements a full three dimensional version which can be extended to many potentials and continuum

meshes. At the many stages of the implementation, the generality of the computational and algorithmic methods were treated with care to make it usable for real physical conditions and computational sizes beyond just simple academic examples. This is achieved by coupling our in-house extended finite element code to the LAMMPS molecular simulator [19] which enables usage of many great capabilities of LAMMPS. This three dimensional method enables us to study many complex phenomena which is not even observed in 2D simulations. This paper is based upon the manuscript by Talebi et al. [20]; however the current paper includes several improvements to the original research. First of all, the text and the figures has been improved in many instances to present the method and the results more clear and coherent. Also, the implementation has been improved so it can be used to model materials with other pair potentials such as EAM, MEAM or Stillinger–Weber. An efficient neighbor search method has been also proposed to locate atoms in the finite element region. Moreover, a global unit system is created in the finite element solver to match with the unit system chosen in the molecular dynamic solver. Though similar, the two examples are simulated using two different unit systems. These improvements are essential to propose methods to adaptively refine the coarse scale region when the discontinuities evolve. In the end, a simple but effective method is proposed to locate the discontinuities approaching the bridging domain that finally assist us to refine the atomistic region. This paper is organized as follows. Section 2 outlines the governing equations. Crack modeling in the continuum domain using the XFEM and the coupling method is discussed briefly in this section. Some implementation details about coupling molecular dynamic code to the finite element code are discussed in the Section 3. Section 4 presents two numerical examples to model crack propagation involving crack tip plasticity problems. Finally, the conclusion of this research is presented in the Section 5.

2. Model description and governing equations 2.1. Definitions Consider a domain X where an existing crack with the surface @ Xc is to be simulated. We define a fine scale region Xfs modeled by molecular dynamics (MD) and place its center near the crack front. The rest of the domain is the coarse scale region Xcs where the system is described through a usual continuum formulation.1 The number of atoms in sub-domain Xfs is denoted by nfs . The outer cs cs boundary of Xcs is denoted by @ Xcs with @ Xcs ¼ @ Xcs t [ @ X u [ @ Xc cs cs cs cs cs cs and @ Xt \ @ Xu ¼ £; @ Xc \ @ Xu ¼ £; @ Xt \ @ Xc ¼ £; subscripts u; t and c indicate ‘displacement-’, ‘traction-’ and ‘crack-’, respectively. This molecular treatment of the crack front zone enables automatic description of voids, dislocations and cracks with the motion of atoms. The mentioned defects are known to strongly influence the initiation and propagation of cracks, see for example [21,22,23,24]. In this study, a base crystalline material with a given lattice structure is assumed everywhere. The fine scale region Xfs is composed of the atoms of this base lattice structure and the finite elements in this region are deactivated (Fig. 1). In the continuum region, Xcs , the continuum model is activated where the material behavior is approximated by the Cauchy–Born rule [4]. The idea of the Cauchy–Born rule is to compute the material properties and stresses from the same potential used in the fine scale domain Xfs . In other words at every point of the continuum model, we have 1 In this paper, the superscripts fs and cs denote the fine scale and coarse scale, respectively.

Please cite this article in press as: Talebi H et al. Concurrent multiscale modeling of three dimensional crack and dislocation propagation. Adv Eng Softw (2014), http://dx.doi.org/10.1016/j.advengsoft.2014.09.016

3

H. Talebi et al. / Advances in Engineering Software xxx (2014) xxx–xxx

Assuming the external potential is resulting only from a constant external force, f ext , and a pair-wise potential, the total potential can be expressed as:

X X fs W fs ¼ W fsext þ W fsint ¼  f ext;a da þ W ðxa ; xb Þ a

ð5Þ

a;b>a

Throughout this manuscript, the subscripts ext and int stand for the external and internal properties, forces and so forth. The canonical form of Hamiltonian equations in fine scale can be expressed as:

8 fs fs @H < p_ a ¼  @x ¼  @W @xa a : x_ a ¼ d_ a ¼

@H fs @pa

ð6Þ

pfs

¼ maa

Finally, these two equations can be combined to: fs

@W fs @W fsext @W int ma d€a ¼  ¼  ¼ f ext  f int @xa @da @da

ð7Þ

In this equation, f int is the internal force exerted on atom a. 2.3. Coarse scale formulation Fig. 1. The solution domain including the crack and the bridging domain.

a very small atomistic model with the same lattice structure and force fields as the fine scale. The sub-domain Xcs \ Xfs where the continuum and atomistic descriptions cohabit is denoted by XB that is called the bridging domain, handshake region or the blending domain. In the bridging domain, the continuum and the atomistic strain energies are blended through a weighting by complementary functions. One criterion which should be fulfilled is that the sum of the continuum and the atomistic energies is unchanged by the coupling method. From a kinematic point of view, in this sub-domain the continuum model is coupled to the atomistic model by enforcing compatibility of the displacements of the atoms and that of the nodes. 2.2. Fine scale formulation

@P € þ q0 b ¼ q0 u @X

X 1 pfsa  pfsa þ W fs ðxa ðtÞÞ a 2ma

ð1Þ

ð8Þ

where P is the first Piola–Kirchhoff stress tensor, b is the body force vector per unit mass, q0 is the initial density and u is the displacement vector. The first Piola–Kirchhoff stress tensor can be calculated from a continuum potential:

P¼

@wcs ðFÞ @F

ð9Þ

where wcs is the potential energy per unit volume and F is the deformation gradient tensor. The total potential energy of the coarse model is:

W cs int ¼

In an isolated system, the sum of the potential and kinetic energies is constant in time. This summation is identified as the Hamiltonian and for the atomistic subdomain (Hfs ), it can be expressed as:

Hfs ðxa ðtÞ; pfsa ðtÞÞ ¼

Let the reference and the current configurations of the domain be denoted by X0 and X, respectively. Similar to the fine scale, the material coordinates of a point in Xcs 0 are denoted by X and the spatial coordinates by x. The linear momentum equation is:

Z Xcs 0

wcs ðFÞdXcs 0

ð10Þ

In the coarse scale, the Hamiltonian is given by:

Hcs ¼ K cs þ W cs ¼

Z Xcs 0

1 qv T v dXcs0 þ W cs 2

Z X cs W cs ¼ W cs f ext;I uI þ ext þ W int ¼  I

where xa is the current position vector and ma is mass of atom a. The location of atom a in the reference (denote by X) and spacial configurations (denote by x) can be related by displacement vector d:

xa ¼ X a þ da

ð2Þ

pa is the momentum of atom a and defined by

pfsa ¼ ma x_ a ¼ ma d_ a

ð3Þ fs

where the dot stands for the time derivatives. W ðxÞ is the potential function which can be due to any kind of force field, such as pairwise interactions or three-body potentials:

W fs ðxa Þ ¼

X X W 1 ðxa Þ þ W 2 ðxa ; xb Þ þ . . . a

a;b>a

ð4Þ

wcs ðFÞdXcs 0

ð12Þ

where v is the velocity vector and K cs is the kinetic energy in the coarse scale. In the continuum sub-domain, the displacement field is approximated by the extended finite element method (XFEM). In this approximation, the displacement field is decomposed into a continuous part and a discontinuous part as in [9,25]:

uh ð X Þ ¼

fs

Xcs 0

ð11Þ

X

X NI ðX ÞuI þ NI ðX ÞSðf I ðX ÞÞaI I2ncs b |fflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflffl} |fflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl ffl{zfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl} I2ncs

ucont

ð13Þ

udiscont

where N is the ordinary finite element shape function, ncs is the set of all nodes in the continuum domain and ncs b is the set of nodes that belong to all elements which are completely cut by the crack. The nodal parameters uI and aI are standard and enriched degrees of freedom, respectively and S is the discontinuous enrichment (Heaviside) function:

Please cite this article in press as: Talebi H et al. Concurrent multiscale modeling of three dimensional crack and dislocation propagation. Adv Eng Softw (2014), http://dx.doi.org/10.1016/j.advengsoft.2014.09.016

4

H. Talebi et al. / Advances in Engineering Software xxx (2014) xxx–xxx

Sðf ðX ÞÞ ¼



1 if f ðX Þ > 0 0

if f ðX Þ < 0

ð14Þ

placements. In the Lagrange multiplier method, the total Hamiltonian is written as:

H L ¼ H þ kT g

with

f ðX Þ ¼ sign½n  ðX I  X Þ min kX I  X k X I 2@ Xc

ð15Þ

where n is the normal to the crack plane. In this study since the crack front is modeled by the atomistic region, no crack-front enrichment is needed.

ð19Þ

where k is Lagrange multipliers vector and g is the gap vector between coarse scale displacement and fine scale displacement. To compute the Lagrange multiplier unknowns we use the method described by Xiao et al. [26]. Finally, we obtain the following semidiscrete equations of motion. In the following, the upper case indices indicate the finite element nodes, e.g. I.

€ J ¼ f ext;I  f int;I þ f kcs 8I 2 ncs MIJ u I ;

2.4. Coupling method

fs € fs

fs

As mentioned before, in the bridging domain method, the total energy of the system is a weighted contributions of the fine and coarse models in the bridging domain XB . To implement this concept, a scalar weight function, w was defined which is unity outside the fine scale domain, zero inside the fine scale and smooth in the blending region:

8 8X 2 Xcs n Xfs > <1 wðXÞ ¼ ½0; 1 8X 2 XB > : 0 8X 2 Xfs n Xcs :

ð16Þ

At any point X in XB ; w can be computed by a normalized distance function:

lðXÞ wðXÞ ¼ l0

fs

H ¼ ð1  wÞH þ wH ¼

þ wW

ð20Þ ð21Þ

€ fs are the accelerations of node J and atom a, respec€ J and d where u a tively. Also, the mass matrix is computed by:

8I; J 2 ncs : MIJ ¼

Z Xcs 0

ð1  wÞq0 NI NJ dXcs 0;

ð22Þ

The mass matrix in the coarse scale is diagonalized according to the mass-lumping scheme for XFEM which was proposed by Menouillard et al. [27]. The internal forces in the coarse scale are:

8I 2 ncs ; f int;I ¼

Z Xcs 0

ð1  wÞP

@NI ðXÞ cs dX0 ; @X

ð23Þ

The forces on each atom are determined from the interatomic potential W as:

    @W rab wðX a Þ þ w X b ; fs 2 @db

X1

8a 2 s1 . . . nfs t; : f fsa ¼ 

b

ð24Þ

where b ranges over all atoms within the cutoff radius of atom a. The forces f kcs in the coarse scale and f kfs in the fine scale, due to the coupling are given by:

8I 2 ncs : f kcs ¼ I

X

ka NI ðX a Þ;

ð25Þ

a2XB0

cs

and

X X pfsa  pfsa pcs  pcs I ð1  wðX a ÞÞ þ ð1  wÞW fs þ wðX I Þ I 2ma 2M I a I cs

kfs

8a 2 s1; . . . ; n t; ma da ¼ f a þ f a ;

ð17Þ

where lðXÞ is the orthogonal projection of X on the interior boundary of the coarse domain Xcs and l0 is the length of this orthogonal projection to the boundary of the fine scale Xfs , Fig. 2. The governing equations are derived from the Hamiltonian of the systems, H, which is the sum of the Hamiltonians of each subdomain:

fs

8a 2 s1 . . . nfs t : f kfs a ¼ ka : ð18Þ

where Hfs and Hcs are the Hamiltonians of the fine and coarse subdomains. W cs and W fs are total potential of coarse and fine scales and are defined in Eqs. (12) and (5), respectively. The Coarse and fine scale domains are constrained in the bridging region, XB by the Lagrange multiplier method. In this overlapping domain, the fine scale displacements are required to conform coarse scale dis-

ð26Þ

To find the coupling forces, ka , we need to compute the coupling matrix which is explained in more details in [12]. We also lump the coupling matrix and therefore the presented method to solve the coupled system in time and enforcing the constraints using the Lagrange multipliers, does not introduce any additional number of equations to the system. 2.5. Continuum constitutive modeling There are several methods to compute the material properties of the single crystal when modeled with continuum approximation. In the continuum-atomistic model, the important step is to find a correspondence between an atomistic energy function Ea and a specific strain energy density W 0 [4,28]. By the assumption that the individual atomic contributions to the total energy can be defined and that the energy of each atom a is uniformly distributed over the volume V a of its Voronoi polyhedron in Fig. 3, both energies can be related as follows [29]: at

W0 ¼

Fig. 2. The weighting function in the handshake domain in two dimensions.

n Eint 1 X ¼ Ea ðra1 ; . . . ; r anat Þ ¼ W 0 ðr a1 ; . . . ; ranat Þ nat V a nat V a a¼1

ð27Þ

where nat is the number of atoms in the lattice system. To establish a relation between atomic distance vectors and continuum deformation, 1st-order Cauchy-Born rule is adopted here. The idea of CBR comes from considering a homogeneous deformations of an

Please cite this article in press as: Talebi H et al. Concurrent multiscale modeling of three dimensional crack and dislocation propagation. Adv Eng Softw (2014), http://dx.doi.org/10.1016/j.advengsoft.2014.09.016

H. Talebi et al. / Advances in Engineering Software xxx (2014) xxx–xxx

5

Fig. 3. The 1st-order Cauchy–Born rule for the homogeneous deformation.

infinite representative crystallite. Suppose that r ab and Rab are lattice vectors in spatial and material coordinates. By defining F ¼ rX u as the local deformation gradient,

rab ¼ FRab

ð28Þ

Note that the distance vectors Rab in material coordinate depend only on initial geometrical crystal lattice structure. For the case of homogeneous deformation and small strains where the energies of all atoms are equivalent and using the pair potentials (standard Cauchy–Born rule), the relation between interatomic force, f ab ¼ f ab and the stress, P a , can be defined as:

1 X Pa ¼ f ab  Rab 2V a b–a

finite element library called PERMIX [36] that is based on Fortran 2003 standard. The implementation includes several aspects such as (a) The connection between PERMIX and LAMMPS to exchange data (the F03LAMMPS2 interface), (b) the information about the atoms which need to be saved such as weighting values (USER–BDM), (c) applying the constraints and (d) the Cauchy–Born method to compute continuum stresses from atomistic data. These aspect are explained in the sequel. (a) F03LAMMPS: This is a Fortran 2003 interface using the ISO_C_BINDING module. This interface allows accessing the LAMMPS object directly from PERMIX. Every class of LAMMPS is mapped with Fortran pointers and an object oriented Fortran structure is created using the Fortran derived data types and type bound procedures. For example a two dimensional matrix in C++ to store the atomic coordinates is saved withdouble ⁄⁄x. In Fortran a similar pointer is made by real(kid = C_DOUBLE),pointer :: x(:,:) which is then mapped to the C++ pointer using the C_F_POINTER function having the starting address of the C++ pointer. The first index of the coordinate matrix refers to the X, Y or Z coordinate and the second refers to the atom number. This way a very fast connection is made between the two packages. (b) USER–BDM: To save the information about the atoms which are essential to apply the boundary conditions in LAMMPS, a new atomic vector is implemented inside the LAMMPS package. This atomic vector adds the necessary matrices and vectors to store the weighting values of the atoms, relative position of the atoms with respect to the continuum crack faces and initial position of the atoms. Also, a new potentials such as Lennard–Jones and EAM are implemented in LAMMPS to account for weighting of the atoms in the bridging domain. (c) Applying the constraints: The computation of the tmatrix and applying the coupling constraints is implemented via an Interaction class in PERMIX. The interaction class handles the initial set up as well. The initial setup procedure will be explained later. (d) The Cauchy–Born method: To implement the Cauchy–Born rule, a very small atomistic part is defined at the integration point level of the coarse scale which handles stress and stiffness calculation. This is done via a new material model in PERMIX that handles also communication with LAMMPS. In LAMMPS, a compute procedure is implemented which

ð29Þ

The above CBR can be applied to uncomplicated material behavior in zero temperature. In this paper, we apply it to atomistic behavior modeled with the Lennard–Jones potential which is described as:

"

Uab ¼ 4

r r ab

12

 

r

6 #

r ab

ð30Þ

where  is the depth of the potential well and r is the finite distance at which the inter-particle potential is zero. Remark. there are many studies on the applicability of the CBR [30,31] that show this method is not stable or accurate in the elevated temperature and near the domain boundary. There are also corrections to incorporate the temperatures and surface effects [32–34]. Although the improvements are significant, for the current problem a linear elastic anisotropic material model should suffice too since the majority of the complex material behavior is modeled using the atomistic simulation around the crack tip. The latter method is also used in [35].

3. Implementation details The fact that continuum and atomistic (discrete) domains are coupled, creates certain difficulties in the implementation due to the different nature of the information available in both domains. The implementation becomes more challenging when the final intention is to model different materials (metals, polymers, etc.) in three dimensions. For the atomistic part a code named LAMMPS [19] is adopted here. LAMMPS is a general purpose molecular dynamics code which is written in C++ and performs messagepassing via MPI calls. The LAMMPS code is coupled to an extended

2

The source code of the interface can be download from the PERMIX website.

Please cite this article in press as: Talebi H et al. Concurrent multiscale modeling of three dimensional crack and dislocation propagation. Adv Eng Softw (2014), http://dx.doi.org/10.1016/j.advengsoft.2014.09.016

6

H. Talebi et al. / Advances in Engineering Software xxx (2014) xxx–xxx

deforms the small atomistic structure based on the input deformation gradient passed from PERMIX and computes the homogenized stress at the center atom. In this implementation we have refined Lagrange multiplier mesh to the atomistic spacing. The size of the coupling region is based on the finite element mesh and it is one element thick. To set up the problem, the following steps are done: 1. Build the atomistic crystal with LAMMPS. 2. Minimize the potential energy of the atomistic part to find the stable positions of the atoms. 3. Read the finite element (coarse scale) model definitions. 4. Read coupling information from the input. 5. Build the neighbor lists for the coarse scale. 6. Find all the atoms in all elements, i.e. the index of element which contains every atom. 7. Find active elements, bridging elements and nodes. 8. Compute the weights of the nodes and integration points. 9. Find active atoms, bridging atoms and ghost atoms. 10. Rescale the masses of atoms. 11. Recompute finite element nodal masses according to the weighting. 12. Set up the Lagrange multiplier points. 13. Create the atomistic groups in LAMMPS based on the last step. 14. Compute the coupling matrix. In the current implementation, the fine scale region is located with a cubic box. The elements which are cut by the box are the bridging elements and the element and nodes inside the box are the inactive ones. All the atoms inside the box are the active atoms and those outside are the ghost atoms, i.e. their displacements are extrapolated from the finite element mesh. 3.1. Efficient neighbor searches Similar to LAMMPS, PERMIX has a similar way of domain and neighbor searches based on the bin search method. This is essential

for efficient attribution of nodes and elements to atoms. The method of neighbor searching is illustrated in Fig. 4. Before creating bins, the minimum and maximum size of finite elements and nodal coordinates are computed. The bins are created slightly larger than the maximum element size. Several lists are then created which save the information on the bins containing nodes, bins containing elements and vice versa. To find the bridging elements, Step 7, the elements which intersect with this box are located and the bridging elements are marked. To find the element containing an atom, first, the coarse scale bin that contains the atom is found. Next, all the elements in this bin are checked if they contain the atom. This check is done with computing the local coordinates of the atom according to the element. In our experience, this method is the most efficient way of locating the atoms. This method is also very general and can be used for any type of element and geometry. The computed local coordinate are later used to compute the weighting of the atoms. 3.2. Consistent units Because of extremely small model sizes, in the most molecular dynamics packages the units of quantities are set according to simulation in hand. This saves a lot of computational time because of the floating point operations. The chosen units often introduce additional factors into time integration method, force, velocity and displacement calculation procedure. In most finite element packages on the other hand, the units of input quantities are arbitrary. Since the time integration procedure is done in PERMIX and also to unify the units, we have introduced the same unit system in PERMIX as of the one in LAMMPS. Please refer to the LAMMPS user manual for detailed explanation of the defined units. To solve the example studies in this paper, we have used two unit systems which are called the LJ units and METAL units. The details of these unit systems are listed in Table 1. In the coupling method also, the units have to be accounted to compute correct coupling forces and displacements.

4. The numerical examples 4.1. Edge crack in a finite plate Consider a three dimensional single crystal with a face centered cubic (FCC) lattice which has dimensions of 8000  800  100 Å. In this example a straight crack of length 620 Å is assumed present in the domain across the whole thickness. Along the left and right sides of the specimen in the continuum domain, the movement of nodes in the Y and Z direction is fixed. The continuum model consists of 197,743 hexagonal elements and 692,064 degrees of freedom. The element size is constant over the domain that is about 15 Å. An atomistic domain of size 580  365  100 Å is placed centered around the crack front. Fig. 5 shows a schematic configuration of the system. Since part of the crack falls within the atomistic domain, the crack must be modeled in the atomistic region as well as in the continuum region. We do not follow the generally adopted method of removing rows of atoms along the crack, as this is somewhat arbitrary and introduces extra parameters in the formulation.

Table 1 The LJ and metal units used in this study.

Fig. 4. Illustration of the neighbor search algorithm to locate an atom in the finite element mesh.

LJ units Metal units

Mass

Distance

Time

Force

Energy

m

r

s

=r



gr mole

Å

femtoseconds

Kcal mole Å

Kcal mole

Please cite this article in press as: Talebi H et al. Concurrent multiscale modeling of three dimensional crack and dislocation propagation. Adv Eng Softw (2014), http://dx.doi.org/10.1016/j.advengsoft.2014.09.016

H. Talebi et al. / Advances in Engineering Software xxx (2014) xxx–xxx

7

the bridging atoms. Consequently, the coupling region is one element wide. With this configuration the model has 1,626,240 active atoms, 166,815 bridging atoms and 148,088 ghost atoms. The driving force for the system is introduced through a velocity boundary condition on the right and left face of the continuum Å region. A velocity 0:1 pico-seconds is set on all the nodes belonging to the right and left boundary of the continuum domain, at every time step. The time step is 0.003 pico seconds (ps). Fig. 6 shows the atoms with higher centro-symmetry value, at different time steps. For the current system the centro-symmetry parameter is a powerful measure of the local lattice disorder around an atom and may be used to characterize and visualize whether the atom belongs to a perfect lattice, a defect (e.g. a stacking fault or a dislocation), or a surface [38]. The centro-symmetry indicator, CS, is computed as explained in [39]: Fig. 5. Schematic view of the example, bridging domain, atomistic domain and boundary conditions for an edge crack in finite plate.

Instead, we modify the neighbor list of the atoms to prevent force transmission across the crack faces. The atomistic domain is a three dimensional lattice from an FCC crystal with lattice constant 3.645 Å. Atomic interactions are modeled by the Lennard–Jones potential with parameters r ¼ 2:29621 Å,  ¼ 0:467 eV, and a cut-off radius of 4.0 Å; the mass gr of all atoms is taken as 63.5 mol . In this example the metal units is used which are explained before. In this study, we have not take the temperature into account, since the focus was on the coupling of (extended) finite element method with Molecular Dynamics. However, in our development, we can easily control the temperature in the atomistic domain with Nose–Hoover thermostat method by [37]. To be able to model a realistic three dimensional problem, where periodicity is usually not available, we have not used any periodic boundary conditions in the system. The coupling of the continuum and atomistic part is performed 3 within a cubic box with of dimensions 540  340  100 Å . The elements which are cut by this box are the bridging elements and the atoms which are located inside bridging elements are

Nnrs 2 2  X   CS ¼ Ra þ RaþNnrs  ;

a¼1

2

ð31Þ

where Nnrs is the number of nearest neighbors for each atom in the underlying lattice of atoms. For example here Nnrs ¼ 12 for the FCC lattice. Ra and RaþNnrs are vectors from the atom of interest to a par2 ticular pair of nearest neighbors. The value in the sum is computed for each atom, and the Nnrs =2 smallest quantities are used. For an atom on a lattice site, surrounded by atoms on a perfect lattice, the centro-symmetry parameter will be 0. It will be also near 0 for small thermal perturbations of a perfect lattice. Fig. 7 shows the stress contours at four different time steps which are the same time steps as in Fig. 6. The atomistic stress computed here is the virial stress tensor. The symmetric virial stress tensor for pair potentials such as the one used here is defined in [40,41]:

" N #

nei 1 X 1X b a ab a a a r ¼ R R F m v v V a 2 b¼1 V

ð32Þ

where b 2 s1; . . . ; Nnei t ranges over the Nnei neighbors of atom a; Ra is the coordinate of atom a; F ab is the force on atom a from atom b; V is the total volume, ma is the mass of atom a and v a is the velocity

Fig. 6. The propagation of the crack front and dislocations shown with atoms with high centro-symmetry value at different time snapshots (a) 97.2 ps, (b) 101.4 ps, (c) 149.4 ps, and (d) 190.2 ps.

Please cite this article in press as: Talebi H et al. Concurrent multiscale modeling of three dimensional crack and dislocation propagation. Adv Eng Softw (2014), http://dx.doi.org/10.1016/j.advengsoft.2014.09.016

8

H. Talebi et al. / Advances in Engineering Software xxx (2014) xxx–xxx

Fig. 7. The stress (rxx ) contour of the specimen at different time snapshots (a) 97.2 ps, (b) 101.4 ps, (c) 149.4 ps, and (d) 190.2 ps

of atom a. The definition of virial stress involves the instantaneous velocities only due to thermal fluctuations. Note that to obtain the equivalent continuum Cauchy stress, the virial stress from the molecular dynamics (MD) simulations has to be averaged over time and space, as explained in [42]. In Fig. 7(a) and (b) a stress concentration is visible, that is initially confined at the crack front; subsequently when propagation occurs, the stress waves are emitted from the crack tip. From Fig. 7(c) and (d) the surface effects of the third dimension are also evident. From this figure we also notice the stress distribution around the dislocations and the crack. Such complex mechanisms of crack and dislocations could not be predicted by any classical continuum description of motion. Another noticeable result of this simulation is the huge computational cost reduction. One could easily compute that with full atomistic simulation, the number of degrees of freedom of the system exceeds 187.5 millions. While the simulations in this paper are performed on a multicore shared memory machine, the distributed memory parallel version of the code will enable us to perform simulations which are several times larger than the system here. This means using this multiscale method, one is able to perform simulations in the accuracy of several billion atoms while keeping the number of degrees of freedom in the order of few hundred millions. 4.2. Double-edge-cracked specimen under combination of shear and tension In this example, a double-edge-cracked specimen was considered. Here, we have two cracks under combination of shear and tension loading. Fig. 8 shows the geometry of example, location

of cracks and boundary conditions. To model the continuum domain, 55,778 elements and 254,016 degrees of freedom are used. In this example the units are changed to the LJ units that saves computation time noticeably. The Lennard–Jones potential parameters are however similar to the ones in the last example with converted values. The dimension of specimen was 2000  2000  20 Å. For atomistic domain, a cubic cell of 1100  1100  20 Å was considered in the middle of continuum domain. Two 750 Å length edge cracks were put in the middle of both lateral sides of the domain. The atomistic domain covers both crack tips. The specimen was fixed at the bottom surface and a combination of shear and tension velocity boundary conditions were applied at the top surface, Fig. 8. The velocities in x and y directions were 0:1 0:1

Å pico-seconds

Å pico-seconds

and

respectively. The material parameters were consid-

ered the same as previous example. The coupling of the continuum and atomistic part is performed within a cubic box with of dimensions 1050  1050  20 Å. Same as previous example, the coupling region is one element wide. Finally, the model has 2,094,840 active atoms, 93,804 bridging atoms and 72,366 ghost atoms. Fig. 9 shows the history of the virial stress in atomistic domain. Fig. 9(a) shows the time of dislocation nucleation. Fig. 9(b) and (c) shows the propagation of dislocation and crack nucleation respectively. Fig. 9 (d) represents the crack propagation direction. It can be realized from Fig. 9(d) on the right hand corner where a dislocation has reached the boundary of the full atomistic domain. The presented method in this paper cannot handle such a case to pass the dislocations to the continuum boundary. One method to handle such a case is presented in the next section.

Please cite this article in press as: Talebi H et al. Concurrent multiscale modeling of three dimensional crack and dislocation propagation. Adv Eng Softw (2014), http://dx.doi.org/10.1016/j.advengsoft.2014.09.016

H. Talebi et al. / Advances in Engineering Software xxx (2014) xxx–xxx

9

Fig. 8. Schematic representation of the double-edge-notched specimen under combined loading, including bridging domain, atomistic domain and boundary conditions.

Fig. 9. The virial stress contours (rxy ) of the atomistic domain at different time snapshots (a) 77.4 ps, (b) 94.8 ps, (c) 131.4 ps, and (d) 225.6 ps.

We would also like to remind that the two examples presented here cannot be replicated with a full continuum model. In order to create a comparable full continuum model, one still needs a reliable criterion for phenomenon such as dislocation nucleation, dislocation climb, crack branching and crack coalescence in the

context of finite element methods. For engineering scaled problems, some empirical rules are applied for example in order to decide whether a crack branches. However, for nano-scale systems, those criteria are not applicable. To propagate the defects in the nano-scale, it is conventional to refine the surroundings of the

Please cite this article in press as: Talebi H et al. Concurrent multiscale modeling of three dimensional crack and dislocation propagation. Adv Eng Softw (2014), http://dx.doi.org/10.1016/j.advengsoft.2014.09.016

10

H. Talebi et al. / Advances in Engineering Software xxx (2014) xxx–xxx

defects down to the atomistic scale. Please refer to the manuscripts by Xiong et al. [43] and by Gracie et al. [18] that present methods to pass the dislocation to the continuum region. In the two mentioned manuscripts, the dislocations are always handled in a full atomistic region and after nucleation and propagation the full continuum scale is used to account for their influence in the system. A comparison of the coupled atomistic-continuum methods is presented in [44]. However, this study focuses on the molecular mechanics and the relevant methods are compared which does not apply here. The comparison of this method to the similar ones such as the Bridging Scale Method [45] is out of the scope of this study. 4.3. Arbitrary refinement of the fine scale The extended bridging domain method is capable of easily refining the coarse scale region. In the current implementation, the evolution of the fine scale is done with adding boxed shaped regions to the original fine scale box. This is done with first locating a hotspot which can be a crack tip approaching the handshake region. Below is the detailed procedure to refine the fine scale region. 1. 2. 3. 4.

Input the size of boxes around the hot spot. Find the element set S cs el1 just inside the fine scale region. fs List the atoms in the S cs el1 to make S 2 . At every nth time increment: (a) Loop over S fs2 to check for the target value D. (b) If D fulfills the refinement criterion then add a new box centered by the atom in S fs2 . (c) In case of new boxes, run the set up procedure again.

Fig. 10 shows an arbitrary refinement procedure. The target value, D, stands for a parameter of interest. D shows occurrence of any complex material behavior in the atomistic region. In crystalline materials, D will show presence of dislocations, voids or cracks. In two dimensions, there are several methods to find the dislocations and compute the burgers vector. However, to refine around a dislocation, it is not necessary to compute the burgers vector or locate the dislocation core and slip plane. Therefore we can use the centro-symmetry criterion, CS as explained in Section 4. The CS value is also higher for cracks and voids therefore it can be used for general case of crystalline materials for the refinement

Fig. 10. The refinement procedure with arbitrary boxes.

process. A similar strategy has been pursued by [46] in the context of an adaptive multiscale methods that automatically transfer fine scale domains into coarse scale domains and vice versa. Efficient coarse graining strategy have been developed subsequently in [47]. 5. Conclusion We presented a coupling method to bridge a three dimensional extended finite element treatment of cracks and molecular dynamics. This method is based on an overlapping domain-decomposition scheme where the displacement compatibility conditions in the overlapping sub-domain are enforced by Lagrange multipliers. We showed how the method can be successfully used to simulate the propagation of a cracks and dislocations, where an atomistic domain is placed on top of the three dimensional extended finite element domain, around the crack front. The propagation of dislocations as well as cracks and the corresponding criteria are handled within the atomistic region automatically. We have observed that our three dimensional coupled method is capable of representing the crack and dislocation propagation for much fewer degrees of freedom than a direct numerical simulation. We also computed the centro-symmetry (CS) parameter and virial stress tensor in the atomistic region. The CS parameters can later be used to refine the atomistic region. We have also proposed a method to refine the continuum region where needed. Acknowledgement The authors thank the support of the German Research Foundation (DFG). References [1] Talebi H, Zi G, Silani M, Samaniego E, Rabczuk T. A simple circular cell method for multilevel finite element analysis. J Appl Math 2012. http://dx.doi.org/ 10.1155/2012/526846. [2] Silani Mohammad, Talebi Hossein, Ziaei-Rad Saeed, Kerfriden Pierre, Bordas Stéphane PA, Rabczuk Timon. Stochastic modelling of clay/epoxy nanocomposites. Compos Struct 2014;118:241–9. [3] Saether E, Yamakov V, Glaessgen EH. An embedded statistical method for coupling molecular dynamics and finite element analyses. Int J Numer Meth Eng 2009;78(11):1292–319. [4] Tadmor Ellad B, Ortiz Michael, Phillips Rob. Quasicontinuum analysis of defects in solids. Philos Magaz A 1996;73(6):1529–63. [5] Abraham F, Broughton JQ, Bernstein N, Kaxiras E. Spanning the length scales in dynamic simulation. Comp Phys 1998;12(6):538–46. [6] Broughton JQ, Abraham FF, Bernstein N, Kaxiras E. Concurrent coupling of length scales: methodology and application. Phys Rev B 1999;60:2391–403. [7] Bayliss A, Turkel E. Mappings and accuracy for chebychev pseudospectral approximations. J Comput Phys 1992;101:349–59. [8] Wagner G, Liu WK. Coupling of atomistic and continuum simulations using a bridging scale decomposition. J Comput Phys 2003;190:249–74. [9] Talebi H, Samaniego C, Samaniego E, Rabczuk T. On the numerical stability and mass-lumping schemes for explicit enriched meshfree methods. Int J Numer Meth Eng 2012;89(8):1009–27. [10] Belytschko T, Xiao SP. A bridging domain method for coupling continua with molecular dynamics. Comp Meth Appl Mech Eng 2004;193:1645–69. [11] Rateau G, Ben Dhia H. The Arlequin method as a flexible engineering design tool. Int J Numer Meth Eng 2005;62:1442–62. [12] Belytschko Ted, Xu Mei. Conservation properties of the bridging domain method for coupled molecular/continuum dynamics. Int J Numer Meth Eng 2008;76:278–94. [13] Gracie R, Belytschko T. Concurrently coupled atomistic and XFEM models for dislocations and cracks. Int J Numer Meth Eng 2009;78(3):354–78. [14] Talebi Hossein, Silani Mohammad, Bordas Stéphane PA, Kerfriden Pierre, Rabczuk Timon. Molecular dynamics/XFEM coupling by a three-dimensional extended bridging domain with applications to dynamic brittle fracture. Int J Multisc Comput Eng 2013;11(6). [15] Shenoy VB, Miller R, Tadmor Eb, Rodney D, Phillips R, Ortiz M. An adaptive finite element approach to atomic-scale mechanics–the quasicontinuum method. J Mech Phys Solids 1999;47(3):611–42. [16] Miller RE, Tadmor EB. The quasicontinuum method: overview, applications and current directions. J Comp-Aided Mater Des 2002;9(3):203–39. [17] Gracie R, Oswald J, Belytschko T. On a new extended finite element method for dislocations: core enrichment and nonlinear formulation. J Mech Phys Solids 2008;56(1):200–14.

Please cite this article in press as: Talebi H et al. Concurrent multiscale modeling of three dimensional crack and dislocation propagation. Adv Eng Softw (2014), http://dx.doi.org/10.1016/j.advengsoft.2014.09.016

H. Talebi et al. / Advances in Engineering Software xxx (2014) xxx–xxx [18] Gracie Robert, Belytschko Ted. An adaptive concurrent multiscale method for the dynamic simulation of dislocations. Int J Numer Meth Eng 2011;86(4– 5):575–97. [19] Plimpton S. Fast parallel algorithms for short-range molecular dynamics. J Comput Phys 1995;117(1):1–19. [20] Talebi M, Silani H, Rabczuk T. The three dimensional extended bridging domain method for brittle fracture. In: Topping BHV, editor. Proceedings of the eighth international conference on engineering computational technology. Stirlingshire (UK): Civil-Comp Press; 2012. Paper 19, doi:http:// dx.doi.org/10.4203/ccp.100.19. [21] Cai Y, Zhuang X, Zhu H. A generalized and ecient method for nite cover generation in the numerical manifold method. Int J Comput Meth 2013;10(5):1350028. [22] Zhuang X, Augarde C, Mathisen K. Fracture modelling using meshless methods and level sets in 3D: framework and modelling. Int J Numer Meth Eng 2012;92:969–98. [23] Lange FF. The interaction of a crack front with a second-phase dispersion. Philos Magaz 1970;22(179):983–92. [24] Ravi-Chandar K, Knauss WG. An experimental investigation into dynamic fracture: II. Microstructural aspects. Int J Fract 1984;26(1):65–80. [25] Belytschko T, Black T. Elastic crack growth in finite elements with minimal remeshing. Int J Numer Meth Eng 1999;45(5):601–20. [26] Xiao SP, Belytschko T. Coupling methods for continuum model with molecular model. Int J Multisc Comput Eng 2003;1:115–26. [27] Menouillard T, Réthoré J, Combescure A, Bung H. Efficient explicit time stepping for the extended finite element method (X-FEM). Int J Numer Meth Eng 2006;68(9):911–39. [28] Shenoy VB, Miller R, Rodney D, Phillips R, Ortiz M, et al. An adaptive finite element approach to atomic-scale mechanics–the quasicontinuum method. J Mech Phys Solids 1999;47(3):611–42. [29] Tadmor EB. The quasicontinuum method, Ph.D. thesis. PhD thesis, Brown University; 1996. [30] Aghaei A, Abdolhosseini Qomi MJ, Kazemi MT, Khoei AR. Stability and sizedependency of cauchy-born hypothesis in three-dimensional applications. Int J Solids Struct 2009;46(9):1925–36. [31] Khoei AR, Qomi MJ, Kazemi MT, Aghaei A. An investigation on the validity of cauchy-born hypothesis using Sutton–Chen many-body potential. Comput Mater Sci 2009;44(3):999–1006. [32] Xiao Shaoping, Yang Weixuan. Temperature-related cauchyborn rule for multiscale modeling of crystalline solids. Comput Mater Sci 2006;37(3):374–9. http://dx.doi.org/10.1016/j.commatsci.2005.09.007.

11

[33] Park Harold S, Klein Patrick A. A surface cauchy-born model for silicon nanostructures. Comp Meth Appl Mech Eng 2008;197(4142):3249–60. http:// dx.doi.org/10.1016/j.cma.2007.12.004. [34] Khoei AR, Aramoon A. A multi-scale modeling of surface effect via the modified boundary cauchyborn model. Mater Sci Eng: C 2012;32(7):1993–2000. http:// dx.doi.org/10.1016/j.msec.2012.05.025. [35] Anciaux G, Ramisetti SB, Molinari JF. A finite temperature bridging domain method for Md–Fe coupling and application to a contact problem. Comp Meth Appl Mech Eng 2012;205–208:204–12. [36] Talebi Hossein, Silani Mohammad, Bordas Stéphane PA, Kerfriden Pierre, Rabczuk Timon. A computational library for multiscale modeling of material failure. Comput Mech 2014;53(5):1047–71. [37] Nose S. Constant-temperature molecular dynamics. J Phys: Cond Matter 1990;2:115–9. [38] Mi C, Buttry DA, Sharma P, Kouris DA. Atomistic insights into dislocationbased mechanisms of void growth and coalescence. J Mech Phys Solids 2011;59:1858–71. [39] Kelchner CL, Plimpton SJ, Hamilton JC. Dislocation nucleation and defect structure during surface indentation. Phys Rev B 1998;58(17):11085. [40] Robert JS. Comments on virial theorems for bounded systems. Am J Phys 1983;51:940–2. [41] Subramaniyan AK, Sun CT. Continuum interpretation of virial stress in molecular simulations. Int J Solids Struct 2008;45(14–15):4340–6. [42] Buehler MJ. Atomistic modeling of materials failure. Springer Verlag; 2008. [43] Xiong Liming, Deng Qian, Tucker Garritt, McDowell David L, Chen Youping. A concurrent scheme for passing dislocations from atomistic to continuum domains. Acta Mater 2012;60(3):899–913. [44] Miller RE, Tadmor EB. A unified framework and performance benchmark of fourteen multiscale atomistic/continuum coupling methods. Model Simul Mater Sci Eng 2009;17:053001. [45] Liu Wing Kam, Park Harold S, Qian Dong, Karpov Eduard G, Kadowaki Hiroshi, Wagner Gregory J. Bridging scale methods for nanomechanics and materials. Comp Meth Appl Mech Eng 2006;195(13):1407–21. [46] Budarapu P, Gracie R, Bordas S, Rabczuk T. An adaptive multiscale method for quasi-static crack growth. Comput Mech 2014;53(6):1129–48. [47] Budarapu P, Gracie R, Shih-Wei Y, Zhuang X, Rabczuk T. Efficient Coarse Graining in Multiscale Modeling of Fracture. Theor Appl Fract Mec 2014;69:126–43.

Please cite this article in press as: Talebi H et al. Concurrent multiscale modeling of three dimensional crack and dislocation propagation. Adv Eng Softw (2014), http://dx.doi.org/10.1016/j.advengsoft.2014.09.016