A comparative study of fracture in Al: quantum mechanical vs. empirical atomistic description
aa r X i v : . [ c ond - m a t . m e s - h a ll ] J u l A comparative study of fracture in Al: quantummechanical vs. empirical atomistic description
Qing Peng and Gang Lu
Department of Physics and Astronomy, California State University Northridge,Northridge, CA, USA
Abstract
A comparative study of fracture in Al is carried out by using quantum me-chanical and empirical atomistic description of atomic interaction at cracktip. The former is accomplished with the density functional theory (DFT)based Quasicontinuum method (QCDFT) and the latter with the originalQuasicontinuum method (EAM-QC). Aside from quantitative differences, thetwo descriptions also yield qualitatively distinctive fracture behavior. WhileEAM-QC predicts a straight crack front and a micro-twinning at the cracktip, QCDFT finds a more rounded crack profile and the absence of twinning.Although many dislocations are emitted from the crack tip in EAM-QC, theyall glide on a single slip plane. In contrast, only two dislocations are nucle-ated under the maximum load applied in QCDFT, and they glide on twoadjacent slip planes. The electron charge density develops “sharp corners”at the crack tip in EAM-QC, while it is smoother in QCDFT. The physicsunderlying these differences is discussed.
Key words:
Plastic Deformation, Dislocations, First-Principles ElectronStructure Theory, Atomistic Simulation, fracture mechanics
PACS:
1. Introduction
Understanding fracture behavior in materials is a challenging undertak-ing. Despite nearly a century of study, several important issues remain un-solved. For example, there is little fundamental understanding of brittle toductile transition as a function of temperature in many materials; there isstill no definitive explanation of how fracture stress is transmitted through
Preprint submitted to JMPS May 29, 2018 lastic zones at crack tips; and there is no complete understanding of thedisagreement between theory and experiment regarding the limiting speed ofcrack propagation. These difficulties to a great extent stem from the fact thatfracture phenomena are governed by processes occurring over a wide range oflength and time scales; these processes are all connected and all contributeto the total fracture energy (Van der Giessen and Needleman, 2002).As emphasized by Van de Giessen and Needleman, although the atomisticinteraction at a crack tip may only account for a small fraction of the totalfracture energy, it can be a controlling factor (Van der Giessen and Needleman,2002) - after all, all fractures take place by breaking atomic bonds. Criticalatomistic information, such as surface energy, stacking fault energy, dislo-cation nucleation/propagation energy and twinning formation energy, etc,has been known to play central roles in fracture. In fact, some of thesequantities are at the heart of fracture mechanics, including the Griffith’s cri-terion for brittle fracture (Griffith, 1920), Rice’s criterion (Rice, 1992) forcrack tip blunting and more recently a criterion for twinning at crack tip(Tadmor and Hai, 2003), to name a few.Because of the inherent multiscale nature of fracture - the process ateach scale depends strongly on what happens at the other scales, the mod-eling and simulation of fracture calls for concurrent multiscale approaches(Lu and Kaxiras, 2005). One of the first concurrent multiscale modeling offracture was based on Macroscopic Atomistic ab initio
Dynamics (MAAD)method for Silicon (Broughton et al., 1999). MAAD couples a quantummechanical description of atoms at crack tip, to an empirical (or classical)atomistic description of atoms at a short distance away from the crack tip,and to the continuum finite-element description of the rest of the system.Since MAAD, several other concurrent multiscale methods have been de-veloped, all involving some level of quantum mechanical modeling at thecrack tip (Cs´anyi et al., 2004; Bernstein and Hess, 2003; Ogata et al., 2001).All these methods were developed/applied for Si owning to the followingtechnical reasons: (1) large-scale electronic structure methods such as linear-scaling algorithms are only applicable to covalently-bonded semiconductorslike Si; general approaches for metals still remain elusive; (2) satisfactoryQM/MM coupling schemes for metals were less well developed until re-cently (Bernstein et al., 2009; Liu et al., 2007; Zhang and Lu, 2007). Onthe other hand, concurrent multiscale approaches that do not involves quan-tum mechanics are readily available (Miller et al., 1998; Kohlhoff et al., 1991;Buehler et al., 2006). Among them, Quasicontinuum (QC) method is par-2icularly promising and it has been widely applied to many materials prob-lems, including fracture in metals (Tadmor and Miller, 2005). QC strives toachieve a “seamless” coupling between atomistic and continuum descriptionsand allows quantum mechanical interactions incorporated in a systematicalmanner. For example, although the original QC was based on classical atomicinteractions, significant progress has been made recently to incorporate quan-tum mechanical interactions in the local QC region (Hayes et al., 2006), non-local QC region (Lu et al., 2006) and entire QC system (Peng et al., 2008).The coarse-graining strategy of QC has also been explored to perform large-scale electronic structure calculations (Gavini et al., 2007).Despite the impressive advance in multiscale methodology development,a crucial question remains unanswered. Although it is clear that an atomisticdescription at a crack tip is indispensable for many purposes, it is not wellestablished whether a quantum mechanical description at the crack tip istruly necessary. This is a poignant point given the continuing improvementof empirical potentials and the still heavy costs for quantum simulations. Itis to address this question that motivates the present study. As a first lookat the problem, we focus on crack tip plasticity in Al and compare resultsreceived from a quantum mechanical description vs. an empirical atomisticdescription at the crack tip, both in the framework of QC. Since the atom-istic resolution is only necessary near the crack tip while the linear elasticfracture mechanics boundary conditions can be applied in the far field, QC iswell poised for such fracture simulations. In addition, QC simulations involvequasi-static energy minimization, thus the unrealistically high strain ratescommon to molecular dynamics simulations are avoided. Unfortunately, as aresult thermally activated processes are precluded in QC. Al is chosen in thisstudy because it is relatively inexpensive for density functional theory (DFT)calculations and an excellent embedded-atom-method (EAM) empirical po-tential exists for Al (Ercolessi and Adams, 1994). The goal of this work is toexamine how and why the results received in the empirical simulations differfrom those in quantum mechanical DFT simulations at the crack tip. To thisend, we employ so-called QCDFT method in which the nonlocal atoms atthe crack tip are treated with DFT. The QCDFT results are compared tothose obtained from the original QC method in which all nonlocal atoms aretreated with EAM empirical potential.The structure of the paper is as follows. The methodology is introducedin Section 2 for both QC and QCDFT methods. A semi-infinite crack un-der mode I loading is set up in Section 3.1. The computational parameters3re described in Section 3.2 and the loading procedure for the crack is sum-marized in Section 3.3. The simulation results and analysis are presentedin Section 4 and discussions are given in Section 5. Finally we conclude inSection 6.
2. Methodology
The QC method (Shenoy et al., 1999; Tadmor et al., 1999) is a multi-scale approach (Lu and Kaxiras, 2005) that combines atomistic models withcontinuum theories, and thus offers an advantage over conventional atom-istic simulations in terms of computational efficiency. The idea underlyingthe QC method is that atomistic processes of interest often occur in verysmall spatial domains (e.g., crack tip) while the vast majority of atoms inthe material behave according to well-established continuum theories. Toexploit this fact, the QC method retains atomic resolution only where neces-sary and coarsens to a continuum finite element description elsewhere. Thisis achieved by replacing the full set of N atoms with a small subset of N r “representative atoms” or repatoms ( N r ≪ N ) that approximate the totalenergy through appropriate weighting. The energies of individual repatomsare computed in two different ways depending on the deformation in theirimmediate vicinity. Atoms experiencing large deformation gradients on anatomic-scale are computed in the same way as in a standard fully-atomisticmethod. In QC these atoms are called nonlocal atoms. In contrast, the en-ergies of atoms experiencing a smooth deformation field on the atomic scaleare computed based on the deformation gradient in their vicinity as befittinga continuum model. These atoms are called local atoms. The total energy E tot (which for a classical system can be written as E tot = P Ni =1 E i , with E i the energy of atom i ) is approximated as E QCtot = N nl X i =1 E i ( { q } ) + N loc X j =1 n j E loc j ( { F } ) . (1)The total energy has been divided into two parts: an atomistic region of N nl nonlocal atoms and a continuum region of N loc local atoms ( N nl + N loc = N r ).The original formulation of QC was limited to classical potentials fordescribing interactions between atoms. However, since many materials prop-erties depend crucially on the behavior of electrons, such as bond break-ing/forming at crack tips or defect cores, chemical reactions with impurities,4urface reactions and reconstructions, electron excitation and magnetism, etc,it is desirable to incorporate appropriate quantum mechanical descriptionsinto the QC formalism. QCDFT is one strategy to fill this role. In specific,QCDFT combines the coarse graining idea of QC and the coupling strategyof the quantum mechanics/molecular mechanics (QM/MM) method. Thismethod can capture the electronic structure at the crack tip within the ac-curacy of DFT and at the same time reach the length-scale that is relevantto experiments(Lu et al., 2006; Peng et al., 2008).The original QC formulation assumes that the total energy can be writ-ten as a sum over individual atom energies. This condition is not satisfiedby quantum mechanical models. To address this limitation, in the presentQCDFT approach the nonlocal region is treated by an EAM-based QM/MMcoupling approach (Lu et al., 2006; Liu et al., 2007): the Kokn-Sham den-sity functional theory (KS-DFT) is coupled to EAM with the interactionenergy calculated also by EAM. The local region, on the other hand, is dealtwith by EAM, which is the same energy formulation used in the MM partof the nonlocal region. This makes the passage from the atomistic to contin-uum seamless since the same underlying material description is used in both.This description enables the model to adapt automatically to changing cir-cumstances (e.g. the nucleation of new defects or the migration of existingdefects). The adaptability is one of its main strengths of QC and QCDFT,which is missing in many other multiscale methods.More specifically, in the present QCDFT approach the material of in-terest is partitioned into three distinct types of domains: (1) a nonlocalquantum mechanical DFT region (region I); (2) a nonlocal classical regionwhere classical EAM potentials are used (region II); and (3) a local region(region III) that employs the same EAM potentials as region II. The cou-pling between regions II and III is achieved via the QC formulation, while thecoupling between regions I and II is accomplished by the QM/MM scheme(Choly et al., 2005; Liu et al., 2007). The total energy of the QCDFT systemis then (Lu et al., 2006) E QCDFTtot = E nl [I + II] + P N loc j =1 n j E loc j ( { F } )= E DFT [I] − E EAM [I] + E EAM [I + II] + P N loc j =1 n j E loc j ( { F } ) , (2)where E nl [I + II] is the total energy of the nonlocal region (I and II combinedwith the assumption that region I is embedded within region II), E DFT [I]is the energy of region I in the absence of region II computed with DFT,5
EAM [II] is the energy of region II in the absence of region I computed withEAM, and E EAM [I + II] is the energy of the nonlocal region computed withEAM.Other types of combination with quantum mechanical and classical atom-istic methods may also be implemented in QCDFT. The great advantage ofthe present implementation is its simplicity; it demands nothing beyond whatis required for a DFT calculation and an EAM-QC calculation. Another im-portant practical advantage of QCDFT method is that, if region I containsmany different atomic species while region II contains only one atom type,there is no need to develop reliable EAM potentials that can describe eachspecies and their interactions. This is because if the various species of atomsare well within region I, the energy contributions of these atoms are canceledout in the total energy calculation. This advantage renders the method par-ticularly useful in dealing with impurities, which is an exceedingly difficulttask for empirical potential simulations.The equilibrium structure of the system is obtained by minimizing thetotal energy in Eq. 2 with respect to all degrees of freedom. Because thetime required to evaluate E DFT [I] is considerably more than that requiredfor computation of the other EAM terms in E QCDFTtot , an alternate relaxationscheme turns out to be useful. The total system can be relaxed by usingconjugate gradient approach on the DFT atoms alone, while fully relaxingthe EAM atoms in region II and the displacement field in region III at eachstep. An auxiliary energy function can be defined as E ′ [ { q I } ] ≡ min { q II } , { q III } E QCDFTtot [ { q } ] , (3)which allows for the following relaxation scheme: (i) minimize E QCDFTtot withrespect to the atoms in regions II ( { q II } ) and the atoms in region III ( { q III } ),while holding the atoms in region I fixed; (ii) calculate E QCDFTtot [ { q } ], and theforces on the region I atoms; (iii) perform one step of conjugate gradientminimization of E ′ ; (iv) repeat until the system is relaxed. In this manner,the number of DFT calculations performed is greatly reduced, albeit at theexpense of more EAM and local QC calculations. A number of tests haveshown that the total number of DFT energy calculations for the relaxationof an entire system is about the same as that required for DFT relaxation ofregion I alone. 6 . Computational details IIIIII [111] (a) (b)(c) region I(d)
Figure 1: (Color online) (a) The overview of the entire crack system with finite-elementmesh; (b) A blown-up view of (a) showing the nonlocal region; (c) region I box and atomicstructure at the crack tip; (d) schematic partition of the system into region I, II and III.The x , y and z axis is along [111],[¯110], and [¯1¯12], respectively. All lengths are in ˚A. A semi-infinite crack in a single Al crystal is studied by both QC andQCDFT for comparisons. The crack is made by removing two layers ofatoms with x < y = 0 & 1 .
41 ˚A. The crack plane is (¯110) and inthis orientation, (111) plane is the only active slip plane for dislocationsemitted from the crack tip; all other { } -type planes lie obliquely to thecrack plane and are thus precluded by the imposed plane strain conditions(Tadmor and Hai, 2003). This configuration was used previously in a MDstudy by Hoagland et al. (Hoagland et al., 1990) and an EAM-QC study byHai et al. (Hai and Tadmor, 2003) although the initial crack opening and the7AM potential used are different in these studies. The initial crack openingin the present work is determined based on two competing considerations: (1)it cannot be too narrow otherwise the crack will close and/or large numberof loading steps is required to observe the onset of plasticity; (2) it cannotbe too wide otherwise the DFT region is too large to render calculationsfeasible. Of course, the DFT region has to be large enough to capture thecrucial plasticity events at the tip. The crack is subject to mode I loadingalong y direction shown schematically in Fig. 1(d).The dimension of the system is 6000 × × .
887 ˚A along the x , y , z directions respectively. The system is periodic in z-direction, and has Dirich-let boundary conditions in the other two directions. The system containsover 11 million Al atoms - a size that is well beyond the reach of any full-blown quantum calculation. The schematic overview of the system is shownin panel (a) of Fig. 1 with finite element meshes. Panel (b) is a zoomed-inview of the nonlocal region and panel (c) displays the atomic detail of thecrack tip. Two comparative calculations are carried out for the crack: EAM-QCvs. QCDFT for Al. EAM-QC is the original QC with EAM potential foratomic interactions. In this work, the EAM potential used is rescaled fromthe original “force-matching” EAM (Ercolessi and Adams, 1994) potentialso that it matches precisely the value of the lattice constant and bulk mod-ulus of Al from the DFT calculations (Choly et al., 2005). Although there-scaling changes very little to the original potential, it eliminates latticeparameter mismatch at the QM/MM interface, and thus reduce QM/MMcoupling errors.DFT calculations are carried out with the Vienna Ab-initio SimulationPackage (VASP) (Kresse et al., 1993; Kresse and Hafner, 1994; Kresse and Furthuller,1996a,b) which is based on Kohn-Shem Density Functional Theory (KS-DFT) with the local density approximation and ultrasoft pseudopotentials.A plane-wave cutoff of 129 eV is used in the calculations and the k -pointsare sampled according to the Monkhorst-Pack method (Monkhorst and Pack,1976) with a 1 × × .3. Loading procedure The simulations are performed quasi-statically with displacement bound-ary conditions where the displacement is prescribed as a function of intendedstress intensity factor (SIF) at each loading step. For EAM-QC, the loadingprocedure is outlined as following:(1) For a small initial stress intensity factor K I , the anisotropic Linear ElasticFracture Mechanics (LEFM) solution (Sih and Liebowitz, 1968) u LEFM ( X , K I )is obtained. Each atom is displaced according to u ( X ) = u LEFM ( X , K I ),where u is the displacement field and X is the position of nodes in themodel.(2) The displacement at the model boundaries is fixed except for the cracksurfaces; the positions for all repatoms are obtained by energy minimizationas explained before.(3) The finite-element mesh, the status of repatoms (local vs. nonlocal) andthe neighbor list are updated.(4) The SIF is increased by a small amount ∆ K I as K I = K I + ∆ K I andrepeat from step (1) until the intended SIF is achieved.In this study the increment ∆ K I = 0 .
001 eV / ˚A . is used. The loadingprocedure adopted here follows closely that of reference (Hai and Tadmor,2003).Because DFT calculations are much more expensive than EAM, we useEAM-QC to load the crack until an incipient plasticity is about to take place,at which point QCDFT is switched on. In other words, a QCDFT relaxationstarts from a configuration that is obtained by EAM-QC for K I . QCDFTthen increases the SIF by ∆ K I . This is a reasonable approximation becauseEAM is known to give accurate results for deformations prior the onset ofcrack tip plasticity or the appearance of lattice defects. As will be shownlater, since the critical load for the onset of plasticity from QCDFT is smallerthan that from EAM-QC, the present loading strategy does not run the riskof missing incipient plasticity of QCDFT.
4. Results and Analysis
To establish the validity of present EAM-QC method, we first performEAM-QC calculations for the same crack studied in reference (Hai and Tadmor,2003). We use the same EAM potential and the same initial crack opening as9 a) (b)
Figure 2: EAM-QC results at SIF K I = 0 .
144 eV / ˚A . . (a)The out-of-plane displacement U z ; the displacement contours range from -0.5 ˚A(darkest) to 0.5 ˚A(lightest). (b)Thezoomed-in view of the crack tip atomic structure and finite-element mesh. The open circlerepresents the atomic position. All distances are in ˚A. that in (Hai and Tadmor, 2003). The results are presented in Fig. 4.1 wheretwo edge dislocations are emitted from the crack tip and they glide at thesame slip plane in a symmetrical manner. The critical SIF is 0.144 eV / ˚A . .All these results are identical to those found in (Hai and Tadmor, 2003) andthus validate the present EAM-QC method.Next, we apply EAM-QC to the crack of interest. The crack is loadedquasi-statically and no crack tip plasticity is observed until the SIF reaches K I = 0 .
180 eV / ˚A . . By comparison, the critical SIF for pure brittle cleavagecomputed from the Griffith criterion for this orientation is K Ic = 0 .
205 eV / ˚A . .In Fig. 3, we present the out-of-plane displacement U z as a function of ap-plied K I values. For K I = 0 .
179 eV / ˚A . , although no dislocation is ob-served, significant deformation at the crack tip is clearly visible (panel a).At K I = 0 .
180 eV / ˚A . , the first dislocation is nucleated and subsequentlymoves away from the crack tip. The dissociated 1 / U z , whose non-zero values indicate that the edge dislocation is dissociated into partials.Interestingly, at the same time, a micro-twin is also nucleated from the cracktip (Fig. 5a). The micro-twn is two layers in both length and width, the10 a) (f)(b) (c) (d) (e) Figure 3: The out-of-plane displacement U z obtained from the EAM-QC calculations atSIF K I of (a) 0.179, (b) 0.180, (c) 0.184, (d) 0.186, (e) 0.191, and (f) 0.198 eV / ˚A . respectively. The displacement contours range from -0.5 ˚A(darkest) to 0.5 ˚A(lightest).The crack surfaces are represented by red curves. All distances are in ˚A. K I value isincreased, more dislocations are emitted on the (111) plane and at the sametime , the micro-twin grows in length but not in width. More specifically, as K I increases to 0.184, 0.186, 0.191, 0.198 eV / ˚A . , two, three, four and fivedislocations are emitted from the crack tip and they glide on the same (111)plane, as shown in the panel of (c), (d), (e) and (f) of Fig. 3 respectively.Correspondingly, the micro-twin grows to three, four, five and six layers inlength, respectively. The micro-twin structures of two, three, five and six lay-ers in length are shown in the panel of (a), (b), (c), (d) of Fig. 5 respectively.The width of the micro-twin is not increasing perhaps due to unfavorablestacking fault energy along the usual twinning plane. The reason that thetwinning was not observed in Hai et al. (Hai and Tadmor, 2003) is proba-bly due to the the different initial crack openings rather than the differentEAM potentials used. We have done additional calculations for the presentcrack opening (2 layers) with the same EAM potential used in reference(Hai and Tadmor, 2003) and found a similar twinning at the crack tip. The (cid:0) A (cid:1) ] D ( x ) dislocation densityfitting Figure 4: The dislocation distribution function D(x) as a function of the inverse of thedistance x. The dashed curve is a fit to the filled circles. emitted dislocations share the following characteristics: (1) They are all edgedislocations of 1 / et al. (Hai and Tadmor, 2003) with the similar EAM potential. Butthis value is too large compared to the experimental splitting distance of 5.5˚A, measured by Mills and Stadelmann(Mills and Stadelmann, 1989). The12iscrepancy may be attributed to still too low stacking fault energy of theEAM potential. (2) They are all on the same { } slip plane whose positionis x = a/
2, where is a is (111) inter-plane distance. The active slip plane isslightly ahead of the crack front position at x = 0. The emitted dislocationsmove away from the crack tip and form a pile-up against the local/nonlocalQC interface. For K I = 0 .
198 eV / ˚A . , the dislocation density D ( x ) definedas the number of dislocations per unit distance along pileup line, is foundto be the square root of the inverse of the distance (Hirth and Lothe, 1970).The distance x refers to the distance between the dislocation center and thelocal/nonlocal interface. The fitted curve (dashed line) in Fig. 4 qualita-tively agrees well with the elastic theory where the fitting function is givenas D ( x ) = 0 . /x ) / . (a) (b)(d)(c) Figure 5: Atomic structure at the crack tip from EAM-QC. Some of atomic planes arehighlighted in pink to indicate the twin. Dashed lines represent the two-layer twin. (a) K I = 0.180, (b) K I = 0.184, (c) K I = 0.191, (d) K I = 0.198 eV / ˚A . . All distances are in ˚A. .2. QCDFT calculation of Al The crack is first loaded by EAM-QC until K I reaches 0.169 eV / ˚A . atwhich point QCDFT is started. This critical loading is determined by trialand error; we launch a number of QCDFT calculations at different K I frompreviously relaxed EAM-QC structures and examine whether any crack tipplasticity is taking palce. The smallest K I value that results in incipientplasticity is the critical loading. In comparison, the critical SIF for purebrittle cleavage computed from the Griffith criterion is K Ic = 0 .
267 eV / ˚A . . (a) (b) (c) (d) (e) Figure 6: The out-of-plane displacement U z obtained from QCDFT calculations at SIF K I of (a) 0.169, (b) 0.170, (c) 0.174, (d) 0.175 (e) 0.178 eV / ˚A . respectively. The dis-placement contours range from -0.5 ˚A(darkest) to 0.5 ˚A(lightest). All distances are in˚A. The crack tip behavior of QCDFT is much different. At K I = 0 . / ˚A . , two dissociated edge dislocations are nucleated - one above thecrack plane and one below it. In contrast to EAM-QC, the two dislocationsglide at two adjacent (111) slip planes, as shown in Fig. 7(d). The positions14f the two slip planes are at x = − a/ x = − a/ K I = 0 .
178 eV / ˚A . , there are only two emitted dislocations gliding at 380 ˚A(up) and 258 ˚A (down) away from the crack plane. Due to the computationalcost of QCDFT, we did not pursue more calculations for even larger loadings.However, it is evident from the present study that the crack tip plasticityobserved in QCDFT is qualitatively different from that in EAM-QC. Thedifferences are more striking given the fact that the QCDFT results arecontinued relaxation of EAM-QC configurations.
5. Discussion
Table 1: Relevant quantities calculated by VASP and EAM for bulk Al; the correspondingexperimental values are extrapolated to T=0 K. a γ γ γ γ sf γ us γ ut τ t (˚A) (J / m ) (J / m ) (J / m ) (J / m ) (J / m ) (J / m )EAM 3.99 0.60 0.98 1.10 0.124 0.134 0.180 0.89DFT 3.99 70.93 1.31 1.18 0.148 0.205 0.262 0.93Exp 4.032 1.14 a a a a Estimates for an “average” orientation.
EAM-QC shows that the crack blunts by emitting dislocations glidingon a single slip plane, and the number of emitted dislocations increases asthe SIF is increased. However, QCDFT predicts that two dislocations arenucleated from the crack tip and they glide at two adjacent slip planes.Moreover, the number of emitted dislocations remains the same (two) for allSIFs. These distinctions are highlighted in Fig. 7 where schematic diagrams,15 a) (b) (c) (d) (e) (f)
A CEA BE D
Figure 7: The schematic diagram, displacement contours and atomic structure with finite-elment mesh showing the dislocation pattern at the crack tip. (a), (b) and (c) are theschematic diagram, displacement contour plot and atomic structure respectively fromEAM-QC calculation for K I = 0 .
198 eV / ˚A . . (d), (e) and (f) are the same from QCDFTcalculation for K I = 0 .
169 eV / ˚A . . x = − / a, − a/ a/
2. Thestacking fault between the Shockley partials is indicated by a short red linesegment. The displacement contour plots are the same as those shown inFig. 3 and Fig. 6; they are reproduced here for convenience to the reader.The open circles in the atomic structure plots represent atomic positions.The sheared finite-element mesh in the atomic structure plots is the result ofpassing-by dislocations (the amount of shear corresponds to the net Burgersvector of the dislocations). The crack tip profile is approximated by blue linesegments in the atomic structure plots. In Fig. 7(c) the crack tip profile isapproximated by rectangular line segments AB + BD + DE to model the factthat the crack tip is blunted by emitting dislocations on a single slip plane.The active slip plane is represented BD whose position is at x = a/
2. In Fig.7(f), however, the crack tip profile of QCDFT is approximated by zigzag linesegments AC + CE because two adjacent slip planes are activated, and theirpositions are at x = − a/ x = − / a . BDE A C a Figure 8: A simple model for crack tip profile. The dotted line represents the relevant slipplanes near the crack tip. The rectangular line segments AB + BD + DE and the zigzagsegments AC + CE approximate EAM-QC and QCDFT crack tip profile respectively. To understand the origin of the differences, we resort to a simple modelthat captures the essential features of the crack configuration shown in Fig. 8.17ne can consult Fig. 7 to understand the correspondence between the modeland the actual crack tip configuration. The reason why the rectangular seg-ments AB + BD + DE are preferred in EAM-QC while the zigzag segments AC + CE are favored in QCDFT can be understood from the following sur-face energy analysis. Without loss of generality, we consider here the casewhere two dislocations are emitted for both EAM-QC and QCDFT. Thelength of BC equals to the magnitude of the Burgers vector because a fulldislocation has been emitted upper-ward from A. As a result, AC repre-sents a { } surface. There are two layers of atoms removed in the initialcrack openning, which is one Burgers vector wide. Therefore in additionto a full dislocation emitted downward from E, the length of CD equals totwice of the Burgers vector magnitude. Hence CE is also a { } surface.Therefore BD , AC (and CE ) and AB (and ED ) represents (111), { } and (¯110) surface, respectively. The total surface energy associated with thezigzag and rectangular segments can be written as ( AC + CE ) γ z and( ABγ + BDγ + DEγ ) z , respectively. z is the repeat distance alongz axis. The various surface energies have been calculated by both DFT andrescaled EAM as summarized in Table 1. For this particular crack configu-ration, the surface energy of the rectangular segments is 0.05 eV lower thanthat of the zigzag segments based on EAM energetics. On the other hand, thesurface energy of the zigzag segments is 1.23 eV lower than that of the rect-angular segments according to the DFT energetics (see note ). Therefore,the zigzag segments are preferred in QCDFT while the rectangular segmentsare favored in EAM-QC. The same conclusion holds for other loadings aswell. According to the Peierls criterion of deformation twinning at a crack tip(Tadmor and Hai, 2003), one can define twinnability which is the likelihoodof a material to twin as opposed to slip at the crack tip. The dimensionlesstwinnability can be expressed as (Tadmor and Bernstein, 2004) τ t = [1 . − . γ sf γ us ] r γ us γ ut (4) a = p / a where a is the lattice constant. AB = a , BC = p / a , BC = √ a , DE = 2 a , AC = p / a , CE = √ a , and z = p / a . Taking a = 2.3036 ˚A, the totalsurface energy of the zigzag configuration is 3.67 eV for EAM and 3.93 eV for DFT. Thetotal surface energy of the rectangular segments is 3.62 eV for EAM and 5.16 eV for DFT. γ sf , γ us , and γ ut are intrinsic stacking fault, unstable stacking fault and unsta-ble twinning energy respectively. A material will emit a dislocation beforetwinning if τ t < τ t >
1. Our DFT and EAM calcula-tions find that τ t is less than 1 as shown in Table I, which suggests that notrue deformation twinning along (111) plane be formed, consistent with theQCDFT and EAM-QC simulation results. However, a two-layer-wide micro-twin along (1¯3¯1) plane is formed in EAM-QC while no such twin is presentin QCDFT. The distinction between the QCDFT and EAM-QC results isdue to the large discrepancy of γ us , which is the energy barrier for a leadingpartial nucleation at a crack tip (Rice, 1992). The large DFT value of γ us renders the nucleation of the leading partial difficult, which in turn prohibitsthe formation of deformation twinning. Figure 9: The electron charge density (˚A − ) at the crack tip for (a) EAM-QC at K I =0 .
198 eV / ˚A . and (b) QCDFT at K I = 0 .
169 eV / ˚A . . The density contours range from0 (blue) to 0.26 (red). The blue sphere stands for atomic position and dashed line is anapproximate profile of the charge density contours. .3. Electron density at the crack tip In Fig. 10, we present the electron charge density around the crack tipfor the EAM-QC configuration of K I = 0 .
198 eV / ˚A . (left) and the QCDFTconfiguration of K I = 0 .
169 eV / ˚A . (right). The charge density for EAM-QC is determined by a superposition of atomic densities (obtained by VASP)centered at each EAM atoms. The distortion of charge density due to thedefects is clearly visible. In EAM-QC the atomic bonding is weakened alongthe twin plane while in QCDFT the atomic bonding is significantly disruptedat the center of the crack. More importantly, the charge density profileof QCDFT is smoother than that of EAM-QC as indicated by the dashedcurves. The “sharp” corners of EAM-QC charge density lead to higher kineticenergy of electrons. Since EAM-QC does not involve quantum mechanics, the“sharp corner” is not penalized energetically and thus permissible. Of course,the electron charge density profile reflects the underlying atomic structure:“sharp corners” correspond to a straight crack front thanks to a single activeslip plane; smooth corners correspond to a more rounded crack front.
6. CONCLUSION
In summary, we have carried out a comparative study of fracture in Al byusing two distinctive atomic interactions: quantum mechanical density func-tional theory and empirical embedded atom method. The DFT descriptionof the crack tip is achieved by QCDFT method while the empirical descrip-tion by EAM-QC method. In addition to quantitative differences, qualita-tively different fracture behavior is also observed between the two methods.EAM-QC predicts a more or less rectangular crack tip configuration whileQCDFT yields a more rounded tip profile. The difference is due to the factthat the emitted dislocations glide on a single slip plane in EAM-QC whiletwo adjacent slip planes are active in QCDFT. As the stress intensity factoris increased, more and more dislocations are emitted from the crack tip inEAM-QC while the number of dislocations remains the same up to the max-imum loading applied in QCDFT calculations. A micro-twin is observed atthe crack tip in EAM-QC, but it is absent in QCDFT. The electron densityprofile at the crack tip is also different between EAM-QC and QCDFT. Allthese differences can be understood in terms of defect energetics, includingsurface energy and stacking fault energy.The different results received suggest that the atomic nature of a cracktip is important and an accurate description of the atomic interaction at the20rack tip is indispensable. Although empirical potentials can be developedby fitting to DFT results, it is unlikely they will reproduce all relevant
DFTenergetics. This is particularly so since one does not know a priori whatare the relevant energetics for a given crack. If several chemical species arepresent in a crack tip, the task of fitting a satisfactory potential becomeseven more daunting. Therefore the solution lies at an explicit quantum me-chanical description of the crack tip, most likely in a form of DFT-basedmultiscale modeling, such as QCDFT. The present paper concerns atomisticaspect of fracture which is important for many purposes. However, there areinteresting fracture phenomena that do not depend on atomistic features andthus are beyond the scope of the present paper. Finally, we have not touchedupon fracture dynamics. The questions - such as will finite temperature dy-namics amplify or diminish the differences that we observed and what arethe best strategies to implement dynamics in a multiscale setting - remainunanswered. We hope that the present paper could spawn more researcheffort in answering these questions.
7. ACKNOWLEDGEMENTS
We thank Ellad Tadmor for his assistance of constructing the crack modeland many helpful discussions. The work at California State UniversityNorthridge was supported by NSF PREM grant DMR-0611562 and DoESciDAC grant DE-FC02-06ER25791.