Bridging crack propagation at the atomistic and mesoscopic scale for BCC-Fe with hybrid multiscale methods

Engineering Fracture Mechanics 155 (2016) 166–182

Contents lists available at ScienceDirect

Engineering Fracture Mechanics journal homepage: www.elsevier.com/locate/engfracmech

Bridging crack propagation at the atomistic and mesoscopic scale for BCC-Fe with hybrid multiscale methods Taolong Xu a,b,c,⇑, Ross Stewart b, Jinghong Fan b, Xiangguo Zeng c, Anlin Yao a a

School of Oil and Natural Gas Engineering, Southwest Petroleum University, Chengdu 610500, China Mechanical Engineering, Alfred University, New York 14802, USA c College of Architecture and Environment, Sichuan University, Chengdu 610065, China b

a r t i c l e

i n f o

Article history: Received 1 September 2015 Received in revised form 16 December 2015 Accepted 17 December 2015 Available online 12 January 2016 Keywords: Cohesive zone model BCC-Fe TS curve Crack propagation Hybrid multiscale method

a b s t r a c t The hybrid multiscale method bridging the atomistic and mesoscopic scale is proposed in the current study, which combines the concurrent generalized particle (GP) dynamics method and the bottom-up hierarchical cohesive zone model (CZM) with embedded traction–separation law. The primary purpose is to transfer the GP-obtained physical parameters to the upper mesoscale finite element method (FEM) to investigate the meso- and microscopic crack propagation. The local fracture energy, stress and opening relationships under tension loading during steady-state crack propagation are extracted from the threescale GP model. In this procedure, the scale duality technique is conducted using the GP analog, which can allow material to exist as particles via a lumping process and allows them to decompose into atoms at crack tips and interfaces. Crack extension resistance is detected by coordination vector (CV) snapshots, and energy release rates in the subdomain are used to evaluate the material behavior against the crack propagation when the current crack tip grows. Using the unique cohesive element length, four-scale tension specimen FE models are designed to reveal the accuracy of the intrinsic correlation between the atomistic and mesoscopic scale under a stress intensity factor to study the brittle body-centered cubic (BCC)-Fe fracture behavior. The result appears reasonable and encouraging. Ó 2015 Elsevier Ltd. All rights reserved.

1. Introduction The macroscopic properties of materials are determined by the atomic arrangement and microstructure, which is a fundamental property–structure concept of material science. A detailed investigation of the fracture parameters of engineering materials on the atomistic scale considering the physical sources of damage is non-trivial. These parameters are supplementary in nature for macro testing and for the visualization of defects nucleation and evolution, as evaluated using open-source software, i.e., visual molecular dynamics (VMD) [23] and AtomEye [28]. However, there is a limitation of establishing a strong linkage that separates the atomic scale and the continuum due to the size and time scale disparity. Thus, the multiscale analysis is proposed to more qualitatively and quantitatively establish a link to integrate related physical and geometrical variables at different scales and to considerably reduce the required computational time. From a methodological perspective, multiscale analysis can be classified into hierarchical (sequential) and concurrent (parallel) methods [18]. The former is one-way coupling from a lower scale to an upper scale, which is described in the present study. The latter ⇑ Corresponding author at: School of Oil and Natural Gas Engineering, Southwest Petroleum University, Chengdu 610500, China. E-mail address: [email protected] (T. Xu). http://dx.doi.org/10.1016/j.engfracmech.2015.12.015 0013-7944/Ó 2015 Elsevier Ltd. All rights reserved.

T. Xu et al. / Engineering Fracture Mechanics 155 (2016) 166–182

167

Nomenclature G Go GC GI

strain energy release rate in an elastic–brittle material intrinsic fracture energy cohesive fracture energy energy release rate of mode I w visco-elastic and/or plastic energy term d separation r traction df ; dfail separation at final failure n separation at failure initiate d0 ; dcr ; dfail n rcr critical opening traction a pre-crack length P displacement loading W crack path domain length B thickness of CT model lpz length of the process zone M a parameter that is dependent on various models
t 0 interface strength Ne cohesive element number along the process zone le the cohesive element length E Young’s modulus t Poisson’s ratio K, Kn, Ks, Kt the required stiffness value in FEA code T0 the user-specified initial constitutive thickness of the cohesive element rn the nominal stress in the pure normal mode rs the nominal stress in the first shear direction rt the nominal stress in the second shear direction CN coordination number CV coordination vector S1 scale one of GP model, atomic domain S2 scale two of GP model S3 scale three of GP model

requires full coupling among different scales. Typically, multi-simultaneous algebraic approaches are used in cracking analysis, and quantum mechanics (QM) and molecular dynamics (MD) should be applied to the crack tip with defect nucleation, twin formation, and dislocation emission. Then, the coupling continuum mechanical method, the so-called finite element method, will be embedded at a remote location where the related atoms and element nodes are related in a one-to-one manner. The direct coupling (DC) [8,36] method and quasicontinuum (QC) method [31,38], which are widely cited, have a good implementation regarding the concurrent multiscale method. However, the incompatibilities between non-continuum (atomics) and continuum (FE) domains result in a statically non-equilibrium ‘‘ghost force”, which makes it difficult to avoid precision loss. The coauthor proposed a GP (generalized particle dynamics) method based on the concurrent multiscale analysis scheme that maintains the basic material structure for different GP scale levels, and the strong iteration algorithm is performed. The conclusive validation for the seamless transition of force, deformation, dislocation process, and slip patterns across the interface between the atomistic and particle domains was studied, all of which showed good comparability to the MD prediction [17]. To bridge the crack analysis using the GP model, the cohesive zone model (CZM) for hierarchical multiscale is recast by the commercial finite code ABAQUS. The theory of CZM can be traced back to the early works by Dugdale [16] and Barenblatt [5] in which the concepts of atomistic de-cohesion and the defect process zone are established. Xu and Needleman [41] first related the CZM to the finite element framework; recently, the CZM approach has gained widespread use with its numerical implementation, especially for crack propagation problems. The CZM defines the interfacial traction–separation relationship, and an inverse identification procedure remains the primary experimental approach that requires one to empirically presuppose the shape of the interface separation law according to the expected types of behavior [24]. Currently, improved detection methods [20,32] are used for crack interface separation without any shape assumption. The direct experimental quantification is difficult, and specific cohesive zone laws are often assumed rather than predicted [44]. In addition, the results of the atomisticbased study will augment the findings of different experimental techniques used to characterize the structure, strength, and fracture behavior of bi-material interfaces [19]. By extending the research of atomistic fracture studies [26,34,35], the

168

T. Xu et al. / Engineering Fracture Mechanics 155 (2016) 166–182

atom-to-continuum CZM approaches from the past 15 years are examined. Spearot et al. [37] extracted interface separation constitutive laws based on MD simulations of a planar, tilt nanograin boundary interface subjected to tension and shear deformations between Cu grain boundaries, although that study did not result in an explicit traction and separation relationship at the local fractured surfaces. Yamakov et al. [42] developed an MD model with a center pre-crack under steady-state conditions to analyze the intergranular fracture along a flat symmetric tilt grain boundary of Al and to define the regions for extracting the parameters for the cohesive zone interface elements in an FE simulation. However, the far field tensile (mode I) loading conditions lacked an in-depth implementation of crack growth on the continuum levels via the cohesive zone law. By extending the work by Yamakov et al. [42] and Zhou et al. [45] developed an MD to examine a mixed mode loading (mode I and mode II), leading to interfacial delamination between two brittle materials and derived direct analytical functions relating the local traction, local displacement, and local loading mode mixing. The effect of elastic constant mismatch between adjacent materials was also explored [44]. Based on Zhou’s work, Lloyd et al. [29] implemented the combination of the traction–separation relationship and the FE/CSE (control systems engineering) framework (referred to as a finite element approach or FEA), and a direct comparison was performed regarding its performance by examining the original geometry and loading rate used in MD. Dandekar and Shin [9,10] conducted the parameterization and obtained a traction–separation law via an atomistic simulation for a ductile–brittle interface (Al–SiC, Al–Al2O3) under tensile and shear loadings, and then used a finite element model to predict the stress–strain response of the metal matrix composite under a high strain rate loading, which was then compared with experimental results. Others have continued the atomistic extracted fracture behavior study, for example, Choi and Kim [7] identified the crack-tip cohesive zone constitutive relationships in an isotropic elastic solid. Krull and Yuan [27] studied the failure mechanism using the cohesive traction–separation law in pre-cracked fcc single-crystal aluminum under the mode I loading condition. Zeng and Li [43] attempted to use the cohesive zone model to determine the coarse grained depletion potential in a multiscale fracture simulation. In this work, we devote more attention to the parameterization of a reasonable assessment of local traction–separation (TS) domains that bridge the crack propagation simulation in MD and the FE code. A three-scale atomistic level model for BCC-Fe with an edge pre-crack is developed in which symmetric local domains are designed along the crack path. Here, the high-accuracy of the atomistic crack simulation is essential and is affirmed by the generalized particle dynamics (GP) method. The TS curves with energy release rates for local domains are extracted under tensile loading, and then mesoscopic scale continuum CT models with a specific element length in the process zone are embedded with the bi-linear TS law. The results show supportive evidence for this hybrid multiscale method that can provide an accessible approximate method for fracture analysis of brittle material.

2. Unique atomistic scale model developed by the GP concurrent multiscale method 2.1. Model design The GP concurrent multiscale method proposed by Fan [17] presents effective degree of freedom (DOF) reduction with particles, as well as a natural scale interface between the domains of different scales, which performs calculations of any particle domain in the corresponding atomistic domain. The main assumption of the GP method is that if the basic material structure can be kept the same as the real material at all scales, then the multiscale analysis will be more accurate and more effective. Based on the advantages of the GP method, a three-scale GP model for pre-crack tensile specimens is designed, as shown in Fig. 1. The entire dimensions are 793.893 Å
1190.105 Å
45.721 Å. The crack is located in two GP scales, and its tip with an elliptical shape is at scale one, which is the atomistic scale and is designated as S1; the lengths of the semi-major axis and the minor semi-axis of the elliptical crack are 5.3 Å and 2.7 Å, respectively. Fig. 2 shows the pre-crack GP model for BCC-Fe, which refers to Fig. 1. With sufficiently large system dimensions that are divided into three scale levels, one can clearly detect the particle intensity (Fig. 2(a)). The total number of particles is 291,008, which includes 109,536 iron atoms, represented by S1. The atom number is greatly reduce when the scale duality procedure is used in the GP analog, which can allow the material to exist as particles via a lumping process and allows them to decompose into atoms for domains near a crack tip (Fig. 2(a)). A ladder layer distribution of atoms around the crack tip is due to the elliptical shape design mentioned in Fig. 2(b). Moreover, the crystallographic orientation of the crack plane is
1 0), which is a 45 degree clockwise rotation around the Z-axis, thereby changing the crack front along the [0 0 1] orien(1 tation, resulting in the cleavage plane along (1 1 0). All of these special pre-established steps on BCC-Fe can promote the observability and efficiency during crack initiation and propagation.

2.2. Simulation procedure Based on the embedded-atom method (EAM) [11], the FS-EAM potential [30] is selected for BCC-Fe simulation. Before loading, the GP model is equilibrated using the NPT method for the Berendsen barostat with 0.000101 GPa for 100 ps at 300 K or 50,000 time steps of 2 fs each, thus maintaining a certain temperature (300 K is applied in this simulation) for each GP scale particle and maintaining a constant overall size. Non-PBCs in the X and Y directions and PBCs in the Z direction are set to avoid the thickness effect. Uniaxial strain loading along the Y direction is conducted after the equilibration, with 50

T. Xu et al. / Engineering Fracture Mechanics 155 (2016) 166–182

L: 1190.105 Å W: 793.893 Å T: 45.721 Å l1: 80 Å w1: 524 Å l2: 198.316 Å Lc: 78.216 Å wc: 10.8 Å lc: 58.216 Å a: 5.3 Å b: 2.7 Å S1: Scale one of GP model, atomic domain S2: Scale two of GP model S3: Scale three of GP model Fig. 1. Schematic diagram.

Fig. 2. (a) Three-scale GP model; (b) crystallographic orientation and atomic configuration at the crack tip.

169

170

T. Xu et al. / Engineering Fracture Mechanics 155 (2016) 166–182

Fig. 3. Traction–separation domains along the cracking path.

loading steps and an application of 5% strain. For each loading step, 5000 NVT relaxation processes are carried out to ensure that the global model is sufficiently stable at 300 K circumstance before the next loading step. 2.3. Local traction–separation domain defined For the CZM simulation procedure, the most important property for characteristic crack evolution behavior is the traction–separation law along the crack growth path. Fig. 3 shows 25 local domains symmetrically distributed along the cracking path, and each brick size is 80 Å
20 Å
45.721 Å, as shown in the figure. To obtain detailed physical data in terms of the crack propagation distance, Da, from the initial tip of the pre-setting edge crack of length, 20 rectangle atomistic subdomains I (I = 1, . . . , 25) are divided for the area in front of the origin of the crack tip. The right edge X coordinate of each generic subdomain is XI = (I  1)
20 Å where 20 Å is the width of the subdomain. We take the subdomain I as the fundamental material object to observe the material behavior against the crack propagation when the current crack tip propagates to the subdomain, i.e., when Da = XI. The traction can be extracted as the Y component under Y-strain loading, which is the general practice to obey the traction separation law (TSL) in an MD simulation [42]. 3. Cracking simulation results at the atomistic scale 3.1. Cracking configuration in S1 Fig. 4(a)–(f) records the configurations of crack initiation and propagation, which includes brittle cracking, crack blunting, void nucleation and the growth process. For each local domain, a different color is used to facilitate the identification. Fracture mechanics define the crack initiation as the appearance of new crack surfaces due to the breaking of atomistic bonds at the crack tip where great significance is identified for energy release in the local domain. Cracking begins at 83 ps in ADD02, as shown in Fig. 4(a) when 3% strain is imposed. From Fig. 4(a) to (b), a rapidly extending cracking behavior with a sharp corner (so-called brittle cracking) emerges until reaching 4.5% strain, along with crack tip crossing of ADD02–ADD08 at 101 ps. Subsequently, crack tip blunting is detected from ADD08 to ADD10 until reaching void nucleation between ADD11 and ADD 12, as shown in Fig. 4(c) and (d). Due to the increasing strain load, the void grows quickly in ADD12 and ultimately connects to the leading crack. Then, the cracking path, which is deviating from the symmetry cracking direction, is redefined. In conclusion, the traction–separation laws performed by domains ADD02–ADD15 should be the primary concern rather than ADD16–ADD25. Furthermore, crack characteristics along S1 are also revealed in Fig. 5, with respect to the tip location and the extending speed, which includes four key crack stages: initiation, brittle cracking, crack blunting, void nucleation and connection. These stages exhibited speeds of 1.99 Å/ps, 14.21 Å/ps, 7.03 Å/ps and 2.22 Å/ps, respectively. The brittle cracking stage has the highest growing speed, which is over 2 times faster than the blunting stage. However, the cracking speed in any of the steps is too fast to be identified from a macro viewpoint. It is obvious that crack initiation and the extending process is so uncertain that it is difficult to predict the cracking behavior based on the macro test results, in addition to predicting the cracking speed and displacement simulated by the atomistic approach. Therefore, the energy concept during cracking is tested under a sufficient atomistic simulation using the three-scale GP model, which is represented by traction–separation laws carried out in local domains. The TS law will be transferred to a high scale by embedding the cohesive element in the finite element model. 3.2. Traction–separation (rI–d) curves extraction Two TSL modes can be classified into Fig. 6(a) and (b) based on the brittleness and blunt sharpness of the crack tip, respectively. The distinction can be easily identified from the curve trade, especially during unloading in which the brittleness leads to a faster unloading. In addition, the extreme stress values for blunting (ADD09–ADD15) are slightly greater than that for brittle cracking (ADD02–ADD08), which can be reasonably explained by the greater energy consumption for new crack

T. Xu et al. / Engineering Fracture Mechanics 155 (2016) 166–182

171

(a) t=83ps, ε = 3%

(b) t=101ps, ε = 4.5%

(c) t=106ps, ε = 5%

(d) t=107ps, ε = 5%

(e) t=109ps, ε = 5.5%

(f) t=115ps, ε = 6% Fig. 4. Crack extension process snapshot under Y-strain loading, which includes brittle cracking (a)–(b), crack blunting (c) and void nucleation and growth (d)–(f).

T. Xu et al. / Engineering Fracture Mechanics 155 (2016) 166–182

Location of crack tip (Angstrom)

172

200 Location of crack tip Fit curve of crack speed

150 100

Void nucleation

Crack blunting

50

2.22 Ang/ps 7.03 Ang/ps

0

Brittle cracking

-50

-150

14.21 Ang/ps

Initiation

-100

1.99 Ang/ps

80

85

90

95

100

105

110

115

120

Time (ps) Fig. 5. Tip location and extending speed of the crack.

surfaces in the blunting period. However, the TS law is typically influenced by the void revolution; for example, in ADD12, the cracking resistance decreases faster than in any other blunting domain, which is ascribed to the growth of defects. 3.3. Energy release rate The energy release rate G, introduced by Rice and Wang [33], plays a significant role in the cracking procedure and is defined as the area under the TS curve. The area under the curve was assumed to be as follows:

Z

1

G¼ 0

rðdÞdd ¼ 2cint ¼ A

ð1Þ

where G is the strain energy release rate in an elastic–brittle material. Rice proposed an important concept regarding the area under the TS curve. It can be interpreted as the adhesion energy or work of fracture, i.e., the energy consumed in opening up two new surfaces from a single one [6]. For a mode I fracture, Tvergaard and Hutchinson [40] proposed a similar CZM concept with different TS curves and noted that the shape parameters (d1/dc; d2/dc) of the TS curve are relatively unimportant, revealing that only the maximum stress (rmax) and the area under the curve (C0 ) play a key role. Under the plane strain tensile condition, the integration result for the TS curves (the so-called energy release rate) GI (I = 1, . . . , 15) of ADD01–ADD15 along the crack path is shown in Fig. 7 in which coordination vector plots are embedded into the typical domains. This vector is defined by the sum of all of the atomic distances of the neighbors that were counted to equal the coordination number (CN); thus, a measure of the local density near a crack tip can be easily identified as an energy concept. For atom i, CV is related to its CN within its near-neighbor (NN) radius, rnn. It is defined as a non-dimensional vector by an average covering all of its NN atoms (i.e., j = 1, 2 . . . CN). The CV’s geometric and vector expressions of atom i are shown in Fig. 8 in terms of the position of the vacancy, neighbor atom j, and rnn. Its mathematical definition is given as follows: 20 ADD02 ADD03 ADD04 ADD05 ADD06 ADD07 ADD08

Local Y stress (GPa)

16 14 12 10 8 6 4 2

16 14 12 10 8 6 4 2 0

0 -2

ADD09 ADD10 ADD11 ADD12 ADD13 ADD14 ADD15

18

Local Y stress (GPa)

18

-2 0

10

20

30

40

50

60

Local Y displacement (Angstrom)

70

0

10

20

30

40

50

Local Y displacement (Angstrom)

(a) Brittle cracking Fig. 6. Traction–separation (rI–d) law of cracking domains.

(b) Blunting

60

T. Xu et al. / Engineering Fracture Mechanics 155 (2016) 166–182

173

Fig. 7. Energy release rate GI (I = 1, . . . , 15) for a local domain and the coordinate vector diagram around a crack tip (CV > 0.4).

Fig. 8. CV concept for atom i with a vacancy.

CV 

CN X 1 ð
rj 
ri Þ rnn  CN j¼1

ð2Þ

There is no vacancy when CV ¼ 0; an atom with a larger CV value indicates more vacancies on its one side. Thus, it may have a greater probability to be close to a pore or to be a site for tiny pore nucleation and may have a greater probability to tend toward cracking based on the new surface formation. CV > 0.4 is applied to display the brittle, blunt and defects nucleation process along local domains. The energy release rate GI increases with crack propagation, except that the void appears at ADD12. The maximum energy to be represented at ADD14 indicates that the crack tip has much more defect affordability during the splitting process. The CV technique provides a meaningful method to identify the atomic-affected zone around the crack tip, which maintains a close relationship with the magnitude of the energy release rate, as shown in Fig. 7. A phenomenological method intelligibly describes the energetics concept for the crack tip regarding the atomistic scale, which can lead to a new concept based on the traditional fracture mechanics. A detailed analysis for the energy and traction laws in the subdomains under strain plane (G_PE and S_PE) and plane stress (G_PS and S_PS) is shown in Fig. 9. Using the plane strain condition as an example, the domains from 2 to 15 were selected as the candidate local areas for fracture parameterized in this work because the crack completely propagates across these domains, otherwise the crack will deviate into the following domains. The initial crack tip is located on the left edge of ADD02 (Fig. 4(a)); thus, ADD01 should not be the destination area for crack initiation. The jumping magnitude of the energy and traction predicts the ultimate capability for fracture initiation, i.e., r2-max = 11.92 GPa and G2-max = 11.97 N/m under the plane strain condition. In the propagation process, two different crack growth patterns consisting of brittle cracking (ADD02– ADD08) and crack blunting (ADD09–ADD15), along with the increased rate of crack lengthening occur. The maximum fracture energy and maximum normal traction are larger in brittle fracture than in the blunting process. Blunting starts at

174

T. Xu et al. / Engineering Fracture Mechanics 155 (2016) 166–182

Fig. 9. Critical energy release rate and maximum normal stress for subdomains. Under the plain strain condition, subdomains 2–15 are the valid range for parameter extraction because of the pre-crack cross domain 1 (Fig. 4(a)), and the crack propagation will deviate after domain 15. There are several different propagation stages of the crack tip along the subdomains: (1) brittle cracking (Fig. 4(b)) in domains 2–8; (2) crack blunting (Fig. 4(c)) appears in domains 10–15, during the blunting period, and (3) void nucleation begins in domain 12 (Fig. 4(d) and (e)).

ADD09 in which normal traction values are r8-max = 16.86 GPa > r9-max = 15.88 GPa < r10-max = 16.73 GPa, but no obvious change is observed for fracture energy. For void nucleation in ADD12, the energy release rate G11-max = 39.28 N/m > G12max = 31.03 N/m < G13-max = 39.24 N/m, and no obvious reduction in normal traction results. In conclusion, blunting on a crack tip can be characterized by the sudden decrease in the maximum normal traction in the subdomain, and void nucleation can be identified when the maximum fracture energy suddenly decreases in the subdomain.

4. Critical fracture parameter calculation in the CZM 4.1. CZM description A continuum model under mechanical loading is developed to predict the interfacial delamination against the selected failure criterion, cohesive element and surface, which is contributed by the cohesive zone model and is subsequently used to simulate the fracture process in different types of composites under different loading conditions. From the schematic of the CZM concept shown in Fig. 10(a), the continual material X is divided into two subdomains X1 and X2 by the new crack surface, named cohesive surface Cc , which is partially closed. If the body forces are neglected, then the stress field rij can be related to the external loading f i on the surface boundary Cf and to the closing tractions T i in the material discontinuity by following relationships:

8 @r ij > ¼0 > > @xj > < rij nj ¼ f j >
i > ui ¼ u > > : rij nj ¼ T i

in X ¼ X1 [ X2 on Cf on Cu

ð3Þ

on Cc


i is the prescribed displacement on the boundary Cu . where nj is the outward normal, and u The nonlinearity model embedded in the CZM along crack path Cc is shown in Fig. 10(b), and the entire body volume remains elastic. The peak value rcr of the traction–separation curve is considered to be the adhesive strength of the material, while the area under the curve is associated with the cohesive fracture energy GC that was previously mentioned. This value involves the intrinsic fracture energy (Go) required to break the intrinsic bonding forces and a visco-elastic and/or plastic energy term (w) that accounts for energy dissipation in the surrounding adhesive layer, as follows [25]:

GC ¼ Go þ w

ð4Þ

The fracture process based on the CZM can be summarized in four steps, as illustrated in Fig. 10(b). Before the defect begins, remote area (1) is immersed in the elastic response, and as the load increases, the crack is initiated (2), and then þ propagated to complete failure (3) with the appearance of new traction free crack surfaces C c and Cc (4), with an opening displacement of d ¼ df . Therefore, the cracking constitutive behavior between (2) and (4) should be defined as the Traction– Separation Law (TSL) derived from cohesive surfaces, rather than a linear stress–strain relationship for the bulk material.

T. Xu et al. / Engineering Fracture Mechanics 155 (2016) 166–182

175

Fig. 10. Schematic for the CZM [2].

Three TS laws that include the exponential model, bilinear model and trapezoidal model (Fig. 10(c)) are widely used to characterize the cracking resistance of elastic materials, and the feasibility for implementation in the FE framework is the most studied topic for molecular-dynamics-based multiscale analysis. 4.2. FE simulation 4.2.1. Model development Because of the consistent results between 2D and 3D models [22], a 2D standard compact tension (CT) specimen was designed based on the dimension rule in ASTM E399-05 [1], as shown in Fig. 11(a). The size of the CT model is 1.25 W
1.2 W, the length of the potential crack path and pre-crack is W and a = 0.5 W, respectively, and the cohesive element is embedded as the CZM in the pre-crack domain. Two displacement loadings with a magnitude of P are applied the in the center of the holes along opposite directions. Four CT model sizes listed in Table 1 are simulated in this study to detect the size effect on the fracture toughness. Iron CT for FEA has a modulus of E = 198 GPa and a Poisson’s ratio of t = 0.3. CPE4R and COH2D4 are selected as the plane strain element type for the entire model and the cohesive element, respectively. The meshed CT model with a size of 1000 lm
960 lm is shown in Fig. 11(b). 4.2.2. Cohesive element size definition in the process zone The cohesive element technique was implemented using the commercial FE code ABAQUS. The cohesive element size is typically determined experimentally in which 3–5 cohesive elements are placed per adjacent continuum element [14,15]. However, based on the interfacial strength derived from the literature or via experimentation, the size of the cohesive elements in the mesh can be estimated, which is based on embedding a sufficient number of cohesive elements within the pro-

(a) Dimensions Fig. 11. Compact tension specimen for FE.

(b) Meshed geometrics

176

T. Xu et al. / Engineering Fracture Mechanics 155 (2016) 166–182

Table 1 List of CT models with different dimensions. NO.

Size (lm)

Pre-crack length a (lm)

Displacement loading P (lm)

Element type

Element number

1 2 3 4

250
240 500
480 1000
960 2000
1920

100 200 400 800

50

COH2D4, CPE4R

54221 75092 126805 214436

cess zone (PZ) that develops in front of a crack tip. The length of the PZ is a material property for non-slender bodies and is given by the following equation:

lpz ¼ ME

GC 2 ð
t0 Þ

ð5Þ

where GC is the interface fracture toughness,
t 0 is the interface strength, E is the Young’s modulus of cracking material, and M is a parameter that is dependent on various models, and its value range is from 0.21 to 1.0 [39]; in this study, M = 1.0. Define N e as the cohesive element number along PZ, and the cohesive element length can be given as follows:

le ¼

lpz Ne

ð6Þ

Typically, N e = 3–5 is sufficient for an accurate simulation for mode I crack propagation [39]. Here, consider that Ne = 5, then the cohesive element size for BCC-Fe plane strain crack propagation can be calculated by Eqs. (5) and (6). From the molecular dynamics perspective, the scale of interface strength for BCC-Fe is
t0 2 ½11:919 GPa; 17:99 GPa, where
t0 = 11.919 GPa, M = 1.0, and E = 198 GPa and t = 0.3 was chosen as the length calculation procedure. In addition, the literature value for the interface fracture strength factor of pure iron is K I 2 ½6 MPa m1=2 ; 20 MPa m1=2  [3]. The size of the cohesive element can be estimated as follows:

le ¼

lpz GC K 2 =Eð1  t2 Þ 0:182K 2I ¼ ME ¼E I ¼ 2 ½0:05 lm; 0:5 lm 2 2 2 Ne Ne ð
t 0 Þ ð
t 0 Þ Ne ð
t 0 Þ

ð7Þ

where le ¼ 0:5 lm in this study is based on the calculation precision and guaranteed efficiency. 4.2.3. Traction–separation parameterization for the CZM In the case of cohesive elements with traction–separation behaviors, the parameters characterizing the traction–separation relationship must be specified, including the initial stiffness, damage initiation threshold, and damage evolution properties; these parameters control the deformation of cohesive elements. The stiffness should not be an infinite value in reality because it would otherwise lead to the non-convergence in ABAQUS/Standard and small time increments in Abaqus/Explicit. Therefore, the uncoupled initial elastic stiffness is set as the relationship between tractions and strains, and the strain is defined as the ratio of separation and the constitutive thickness. Diehl [12–15] has proposed an alternative approach to estimate the elastic stiffness of the cohesive element based on classical fracture mechanics and numerical stability considerations. For bilinear CZM as shown in Fig. 10(c), the elastic stiffness K 0 is given in terms of fracture energy GC, separation at final failure df, and the damage initiation ratio dratio as follows:

K 2GC ¼ T 0 dratio d2f

K0  K¼

2GC dratio d2f

ð8Þ ð9Þ

T0

where dratio ¼ d0 =df , and 0 < dratio < 1, and d0 denotes the separation at failure initiate, K is the required stiffness value in FEA code, and T0 is the user-specified initial constitutive thickness of the cohesive element (typically specified as 1.0). It is assumed that the isotropic materials have equal normal and tangential stiffness components for a mode I fracture, and the area of the triangle in Fig. 9(c) represents the fracture energy GC, where the separation at defect initiation and final failure are dcr ¼ dinit and df ¼ dfail n n , respectively. According to Eqs. (8) and (9), the stiffness can be derived as follows:

Kn ¼ Ks ¼ Kt ¼ K ¼

2GC dratio d2f

T0 ¼ 2

dfail n
N max fail fail 2 2ðdinit n =dn Þðdn Þ

T0 ¼

Nmax dinit n

ð10Þ

Maximum nominal stress criterion in terms of TSL is selected as the damage initiation criteria, as follows:



MAX

hrn i rs rt ; ; Nmax Smax T max



¼1

ð11Þ

177

T. Xu et al. / Engineering Fracture Mechanics 155 (2016) 166–182

where rn is the nominal stress in the pure normal mode, rs is the nominal stress in the first shear direction, and rt is the nominal stress in the second shear direction; in addition, rn ¼ rs ¼ rt . The third important parameterization step is damage evolution based on TSL. Linear response is set as the post damageinitiation softening response of bilinear CZM, which can be divided into displacement-based evolution and energy-based evolution. Mode-independent with linear energy failure evolution criteria are used in this study, which depend on the fracture energy GC according to the curvilinear integral, and the viscosity coefficient is set as 1
105 for the stabilization requirement of the cohesive element. Fig. 12 shows the TSL development from each local domain. Table 2 shows the key parameters derived from the GP method as the input values of bilinear CZM, including initial stiffness Kn, Ks, Kt of the cohesive element, maximum damage normal traction rn and the energy release rate GC.

4.2.4. Results and discussion 4.2.4.1. Force–displacement response. The reaction force in the center point of the hole on the CT specimen, as shown in Fig. 11 (a), is extracted as the traction response of CZM, which reaches an extreme value under crack initiation, as shown in Fig. 13 (b) of reaction force vs. displacement response. The fracture toughness calculation depends on the critical load Pc applied to the hole’s center point at the moment of crack initiation. STATUS outputs as 1 or 0 are designated to ensure recording of the time point when the tip cohesive element defect occurs. The crack of the model case 1 will grow at the point of element 52221 deletion, as shown in Fig. 13(c); meanwhile, the STATUS will be 0. Therefore, the accurate reaction force can perform real-time tracking and can detect the procedure of crack initiation and propagation. According to the TSL parameters for each local domain listed in Table 2, the reaction–displacement relationships for various model sizes are evaluated by FEA, as shown in Fig. 14. Domain ADD14 has the highest critical force to meet the crack initiation, showing the maximum local Y stress of 17.99074 GPa and a high energy release rate of 39.24035 lN/lm, leading to the greater cracking resistance.

20

Local Y stress (GPa)

ADD14 Bilinear law (ADD14)

15

10

5

0 0

10

20

30

40

50

Local Y displacement (Angstrom)

(a) ADD14 Bilinear law (ADD02) Bilinear law (ADD03) Bilinear law (ADD04) Bilinear law (ADD05) Bilinear law (ADD06) Bilinear law (ADD07) Bilinear law (ADD08)

15

10

5

20

Local Y Traction (GPa)

Local Y Traction (GPa)

20

Bilinear law (ADD09) Bilinear law (ADD10) Bilinear law (ADD11) Bilinear law (ADD12) Bilinear law (ADD13) Bilinear law (ADD14) Bilinear law (ADD15)

15

10

5

0

0

0

4

8

12

16

20

24

28

Local Y separation (Angstrom)

(b) Brittle cracking (ADD02-ADD08)

32

0

5

10

15

20

25

30

35

40

45

50

Local Y separation (Angstrom)

(c) Crack blunting (ADD09-ADD15)

Fig. 12. Bilinear CZM for local domains ADD01–ADD15.

178

T. Xu et al. / Engineering Fracture Mechanics 155 (2016) 166–182

Table 2 The simulation parameters for CZM from the GP method. Domain ADD02 ADD03 ADD04 ADD05 ADD06 ADD07 ADD08 ADD09 ADD10 ADD11 ADD12 ADD13 ADD14 ADD15

rn (GPa)

Maximum local Y stress

Local strain at maximum rn

Displacement at maximum rn (Å)

Energy release rate GC (lN/lm)

Stiffness coefficients Kn, Ks, Kt (lN/lm3)

11.91907 12.75061 15.15916 14.9497 15.93879 16.45041 16.85554 15.87659 16.7348 17.18728 17.51299 17.60306 17.99074 17.78654

0.063872 0.054329 0.097933 0.073474 0.085025 0.076009 0.075332 0.14363 0.088021 0.099966 0.092053 0.083418 0.08519 0.11434

5.10975 4.34635 7.83467 5.87791 6.80203 6.08069 6.02658 11.49036 7.04165 7.9973 7.3642 6.67341 6.81521 9.14718

11.97013 14.83577 16.8198 19.91512 21.92361 25.8058 28.91036 34.62545 38.31021 39.27735 31.03211 39.24035 41.10038 38.2293

2.332613E+7 2.933636E+7 1.934882E+7 2.54337E+7 2.34324E+7 2.705353E+7 2.796867E+7 1.381731E+7 2.376545E+7 2.149135E+7 2.378125E+7 2.637791E+7 2.639792E+7 1.944483E+7

Fig. 13. FEA results of the CZM. (a) Stress distribution of the CT specimen of model case 1 (250
240 lm) with TSL from ADD014 and the cohesive element configuration at the crack tip before the failure; (b) plot of the reaction force vs. displacement response; (c) the opening cohesive element 52221 at the crack tip in figure, which will be removed once the load exceeds the threshold value.

4.2.4.2. Fracture toughness. In materials science, fracture toughness is a property that describes the ability of a material containing a crack to resist fracture and is one of the most important properties of any material for many design applications. Brittle fracture is characteristic of materials with less fracture toughness [21]. The linear-elastic fracture toughness of a material is determined from the stress intensity factor (K) at which a thin crack in the material begins to grow. Based on the linear-elastic fracture theory, the KI of a standard CT specimen under normal tensile mode can be calculated using the following equation:

pffiffiffiffi ! a  a 2  a 3  a 4

2 þ Wa P K I ¼ pffiffiffiffiffiffi  þ 14:72  5:60  13:32 0:886 þ 4:64  W W W W B W 1  Wa 3=2

ð12Þ

where P is the applied tensile load on the center nodes, as shown in Fig. 11(a). Fig. 15 depicts the critical tensile load P along the local domains of 4 mode cases; a, W are the pre-crack length and crack path domain length of the CT model, respectively; and B is the thickness of the 2D CT model with plain strain (let B = 1.0 to avoid the thickness effect).

179

T. Xu et al. / Engineering Fracture Mechanics 155 (2016) 166–182

50000

44000

45000

42000

40000

2.7

3.0

Force (μN)

36000

25000

3.3

20000 15000 10000

0

1

2

3

4

5

3.6 3.9 4.2 4.5

30000 20000

6

7

8

9

0

10

0

1

2

3

4

5

6

7

Displacement (μm)

Displacement (μm)

(a) 250×240μ m

(b) 500×480μ m

120000

8

9

10

120000

82500

ADD14

100000

100000

80000 5.6

60000

6.3

Force (μN)

75000

80000 60000

7.2

40000

114000

20000

20000

108000

0

0

40000

0

1

2

3

ADD02 ADD09

4

5

6

7

8

9

10 11

8.0

102000 0

1

2

3

4

5

6

7

8

Displacement (μm)

(c) 1000×960μ m

(d) 2000×1920μ m

ADD03 ADD10

ADD04 ADD11

ADD05 ADD12

ADD06 ADD13

ADD07 ADD14

9

110000 100000

250*240 500*480 1000*960 2000*1920

90000 80000 70000 60000 50000 40000 30000

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16

Domain Fig. 15. Critical force at crack initiation along local domains of 4 model cases.

10 11

ADD08 ADD15

Fig. 14. Reaction force vs. displacement for CT specimen of 4 model cases with CZM related to each local domain.

120000

8.8

ADD14

Displacement (μm)

Critical force at crack initiation (μN)

Force (μN)

40000

10000

5000 0

ADD14

54000

50000

38000

30000

58500

60000

40000

35000

Force (μN)

70000

ADD14

T. Xu et al. / Engineering Fracture Mechanics 155 (2016) 166–182

Stress intensity factor (MPa (μm)1/2)

180

40000

32772 31283 29973

35000

29778

30000 25000 20000 15000 10000

250*240 500*480 1000*960 2000*1920

5000 0 0

1

2

3

4

5

6

7

8

9

10

Displacement (μm)

Critical stress intensity factor (MPa ( m)1/2)

Fig. 16. Size effect on KIC (ADD14).

40000 35000

32772 30000 25000 20000 15000

250*240 500*480 1000*960 2000*1920

20616

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16

Domain Fig. 17. Fracture toughness for different micro-models.

In addition, the size effect is detected according to the KIC calculation on the ADD14 domain for four models. Fig. 16 shows the critical stress intensity factor KIC based on Eq. (12). The evaluated values 32,772 MPa (lm)1/2, 31,283 MPa (lm)1/2, 29,973 MPa (lm)1/2, and 29,778 MPa (lm)1/2 are very close, thus revealing the very slight size effect on the results. Generally, the tested fracture toughness of cast iron is 20,000–60,000 MPa (lm)1/2, as described in the literature [4], which shows that the range of simulated KIC values for the micro CT model is 20,616–32,772 MPa (lm)1/2 appears reasonable, although, there is an existing error between different size models. In Fig. 17, the error is reduced, and the curves are of a larger size 1000
960 lm and 2000
1920 lm, from which a conclusion can be drawn that the larger size of the mode has a more stable fracture toughness. 5. Conclusions The hybrid multiscale methodology is analyzed in detail for fracture parameterization that bridges the molecular dynamics method and the finite element method via the TSL embedded in the cohesive element. The adhesion laws in terms of traction–displacement response during crack initiation and propagation are extracted from the atomic forces and motions around the crack tip using the FS-EAM potential, during which the hierarchical GP method with the scale duality technique can preserve the degree of freedom in large quantities. Then, concurrent use of the CZM successfully transfers the fracture parameters from the non-local atomistic model to the continuum standard pre-crack model so that a significant hybrid multiscale bridge is established. In the scale-spanning analysis of this study, the atomistic energy release rate and the maximum traction force play the key role in guiding the cracking behavior. For each subdomain I (I = 1, . . . , 25), three important physical variables are obtained, as follows:

T. Xu et al. / Engineering Fracture Mechanics 155 (2016) 166–182

181

(1) The maximum resistance stress, rI-max, of the material against the crack passing through the subdomain I. (2) The rate of the resistance increases to reach rI-max as measured by the initial slope of the stress-crack surface separation diagram (i.e., rI–d curve); this parameter indicates the material hardening behavior toward crack propagation. (3) The area under the rI–d curve gives the energy release rate GI, which is an important energy variable for crack propagation. The GP results by drawing curves connecting all of the points of rI-max and GI for the 25 subdomains show the following features of the subdomain resistance. The crack starts propagation at the stress of r2-max = 11.92 GPa and G2 = 11.97 J/m2, then the resistance quickly increase to subdomain 3 of the crack propagation length Da = 40 Å with r3-max = 12.75 GPa and G3 = 14.84 J/m2 with a sharp slope in the stress and energy curves. Then, the slope gradually becomes flat until subdomain 11 of Da = 200 Å with r11-max = 17.19 GPa and G11 = 39.28 J/m2. This maximum resistant stress maintains nearly the same value until the crack reaches subdomain 15 with Da = 280 Å. However, GI only maintains a value close to 38–41 J/m2 in subdomains 10–15 where the crack blunting starts and the GI increases until domain 14 with G14 = 41.10 J/m2, except for the void nucleation in subdomain 12, where Da = 220 Å with G12 = 31.03 J/m2. The important conclusion is that the ultimate fracture energy and stress values can be the criteria for void nucleation and crack blunting, respectively. (4) The successful attempt of using the hybrid multiscale covered Å ? nm ? lm ? cm on BCC-Fe under mode I, which directly transferred the GP-determined fracture properties by three data sets for a cohesive finite element analysis of a compact fracture specimen. The credible parameters derived from the atomistic level should be the guarantee, which enable further improvements in view of the GP method. The result appears reasonable and encouraging. Additional studies related to using the hybrid method for the mixed mode loading and rate effects should be conducted in the future. The FE results from a data set at different subdomains and a discussion related to the commonly used cohesive constitutive law and methods should also be performed. Looking back to the past 20 years, several direct coupling (DC) methods between atoms and FE nodes of continuum in the concurrent multiscale analysis have developed to solve defects problems at particle scale level, which is assured the precision of behavior transfer between different scales. Even though, the MD based cohesive element method shows significant improvement on scale bridging, which is considered as the most challenge fact for DC methods. Authors has made some further meaningful works that will be published so called ‘‘GP-FEA” concurrent multiscale method which is based on the preliminary results in this article. Acknowledgement The paper is supported by ‘‘the Young Scholars’ Development Fund of Southwest Petroleum University of China”. References [1] Test method for linear-elastic plane-strain fracture toughness KIc of metallic materials. ASTM International; 2005. [2] Alfano M, Furgiuele F, Leonardi A, Maletta C, Paulino GH. Mode I fracture of adhesive joints using tailored cohesive zone models. Int J Fract 2009;157:193–204. [3] Ashby F, Jones DRH. Engineering materials 1: an introduction to properties, applications and design. Butterworth-Heinemann; 2012. [4] Ashby MF, Shercliff H, Cebon D. Materials: engineering, science, processing and design. Butterworth-Heinemann; 2013. [5] Barenblatt GI. The mathematical theory of equilibrium cracks in brittle fracture. In: Dryden HL, Kármán Tv, Kuerti G, Dungen FHvd, Howarth L, editors. Advances in applied mechanics. Elsevier; 1962. p. 55–129. [6] Chandra N, Li H, Shet C, Ghonem H. Some issues in the application of cohesive zone models for metal–ceramic interfaces. Int J Solids Struct 2002;39:2827–55. [7] Choi ST, Kim KS. Nanoscale planar field projections of atomic decohesion and slip in crystalline solids. Part I. A crack-tip cohesive zone. Philos Mag 2007;87:1889–919. [8] Curtin WA, Miller RE. Atomistic/continuum coupling in computational materials science. Modell Simul Mater Sci Engng 2003;11:R33. [9] Dandekar CR, Shin YC. Effect of porosity on the interface behavior of an Al2O3–aluminum composite: a molecular dynamics study. Compos Sci Technol 2011;71:350–6. [10] Dandekar CR, Shin YC. Molecular dynamics based cohesive zone law for describing Al–SiC interface mechanics. Compos A Appl Sci Manuf 2011;42:355–63. [11] Daw MS, Baskes MI. Embedded-atom method: derivation and application to impurities, surfaces, and other defects in metals. Phys Rev B 1984;29:6443–53. [12] Diehl T. Modeling surface-bonded structures with ABAQUS cohesive elements: beam-type solutions; 2004. [13] Diehl T. Using ABAQUS cohesive elements to model peeling of an epoxy-bonded aluminum strip: a benchmark study for inelastic peel arms. In: 2006 ABAQUS users’ conference; 2006. [14] Diehl T. On using a penalty-based cohesive-zone finite element approach, Part I: elastic solution benchmarks. Int J Adhes Adhes 2008;28:237–55. [15] Diehl T. On using a penalty-based cohesive-zone finite element approach, Part II: inelastic peeling of an epoxy-bonded aluminum strip. Int J Adhes Adhes 2008;28:256–65. [16] Dugdale DS. Yielding of steel sheets containing slits. J Mech Phys Solids 1960;8:100–4. [17] Fan J. Multiscale analysis across atoms/continuum by a generalized particle dynamics method. Multiscale Modell Simul 2009;8:228–53. [18] Fan J. Multiscale analysis of deformation and failure of materials. Chichester, UK: John Wiley & Sons Ltd; 2010. [19] Gall K, Horstemeyer MF, Van Schilfgaarde M, Baskes MI. Atomistic simulations on the tensile debonding of an aluminum–silicon interface. J Mech Phys Solids 2000;48:2183–212. [20] Ghiassi B, Xavier J, Oliveira DV, Lourenço PB. Application of digital image correlation in investigating the bond between FRP and masonry. Compos Struct 2013;106:340–9. [21] Hertzberg RW. Deformation and fracture mechanics of engineering materials. Wiley; 1996. [22] Hibbitt K, Sorensen. ABAQUS/standard user’s manual. Hibbitt, Karlsson & Sorensen; 2001. [23] Humphrey W, Dalke A, Schulten K. VMD: visual molecular dynamics. J Mol Graph 1996;14:33–8.

182

T. Xu et al. / Engineering Fracture Mechanics 155 (2016) 166–182

[24] Jumel J, Ben Salem N, Budzik MK, Shanahan MER. Measurement of interface cohesive stresses and strains evolutions with combined mixed mode crack propagation test and Backface Strain Monitoring measurements. Int J Solids Struct 2014. [25] Kinloch AJ. Adhesion and adhesives: science and technology. Netherlands: Springer; 1987. [26] Komanduri R, Chandrasekaran N, Raff L. Molecular dynamics (MD) simulation of uniaxial tension of some single-crystal cubic metals at nanolevel. Int J Mech Sci 2001;43:2237–60. [27] Krull H, Yuan H. Suggestions to the cohesive traction–separation law from atomistic simulations. Engng Fract Mech 2011;78:525–33. [28] Li J. AtomEye: an efficient atomistic configuration viewer. Modell Simul Mater Sci Engng 2003;11:173. [29] Lloyd JT, Zimmerman JA, Jones RE, Zhou XW, McDowell DL. Finite element analysis of an atomistically derived cohesive model for brittle fracture. Modell Simul Mater Sci Engng 2011;19:065007. [30] Mendelev M, Han S, Srolovitz D, Ackland G, Sun D, Asta M. Development of new interatomic potentials appropriate for crystalline and liquid iron. Philos Mag 2003;83:3977–94. [31] Miller R, Tadmor EB. The quasicontinuum method: overview, applications and current directions. J Comput Aided Mater Des 2002;9:203–39. [32] Réthoré J, Estevez R. Identification of a cohesive zone model from digital images at the micron-scale. J Mech Phys Solids 2013;61:1407–20. [33] Rice JR, Wang J-S. Embrittlement of interfaces by solute segregation. Mater Sci Engng, A 1989;107:23–40. [34] Rottler J, Barsky S, Robbins MO. Cracks and crazes: on calculating the macroscopic fracture energy of glassy polymers from molecular simulations. Phys Rev Lett 2002;89:148304. [35] Rottler J, Robbins MO. Molecular simulations of deformation and failure in bonds formed by glassy polymer adhesives. J Adhes Sci Technol 2003;17:369–81. [36] Saether E, Yamakov V, Glaessgen EH. An embedded statistical method for coupling molecular dynamics and finite element analyses. Int J Numer Methods Engng 2009;78:1292–319. [37] Spearot DE, Jacob KI, McDowell DL. Non-local separation constitutive laws for interfaces and their relation to nanoscale simulations. Mech Mater 2004;36:825–47. [38] Tadmor EB, Ortiz M, Phillips R. Quasicontinuum analysis of defects in solids. Philos Mag A 1996;73:1529–63. [39] Turon A, Dávila CG, Camanho PP, Costa J. An engineering solution for mesh size effects in the simulation of delamination using cohesive zone models. Engng Fract Mech 2007;74:1665–82. [40] Tvergaard V, Hutchinson JW. The relation between crack growth resistance and fracture process parameters in elastic–plastic solids. J Mech Phys Solids 1992;40:1377–97. [41] Xu XP, Needleman A. Numerical simulations of fast crack growth in brittle solids. J Mech Phys Solids 1994;42:1397–434. [42] Yamakov V, Saether E, Phillips DR, Glaessgen EH. Molecular-dynamics simulation-based cohesive zone representation of intergranular fracture processes in aluminum. J Mech Phys Solids 2006;54:1899–928. [43] Zeng X, Li S. A multiscale cohesive zone model and simulations of fractures. Comput Methods Appl Mech Engng 2010;199:547–56. [44] Zhou XW, Moody NR, Jones RE, Zimmerman JA, Reedy ED. Molecular-dynamics-based cohesive zone law for brittle interfacial fracture under mixed loading conditions: effects of elastic constant mismatch. Acta Mater 2009;57:4671–86. [45] Zhou XW, Zimmerman JA, Reedy ED, Moody NR. Molecular dynamics simulation based cohesive surface representation of mixed mode fracture. Mech Mater 2008;40:832–45.