Fracture toughness and bond trapping of grain boundary cracks

The fracture toughness of materials in which brittle cracks initiate and propagate along grain boundaries depends not only on the energy of the grain boundary, but also on its microscopic structure. How exactly and to which degree brittle grain boundary fracture is aﬀected by the local atomic structure at the crack tip has not yet been studied in detail. Here, we use molecular static simulations to study the atomic-scale fracture behavior of six large angle tilt grain boundaries in tungsten bicrystals. The fracture toughness depends critically on the propagation direction and on the position of the crack tip within the structural units of the grain boundary. Furthermore, the grain boundary fracture toughness can be signiﬁcantly larger than for single crystals in the same orientation. These results cannot be explained by the usual thermodynamic approach in continuum-scale fracture mechanics but can be understood by considering the eﬀect of bond trapping of grain boundary cracks.


Introduction
The resistance against crack propagation is undoubtedly one of the most important properties of structural materials. Whether a stressed component fractures by brittle cleavage or by ductile rupture is determined by the competition between bond-breaking at the crack tip and plastic deformation in the vicinity of the crack [1,2]. Which of these processes dominates at a given temperature and strain rate depends on the nature of the chemical bond between the atoms [3] as well as on the microstructure of the material [2,[4][5][6]. In linear elastic fracture mechanics (LEFM), the fracture toughness (or brittle crack initiation toughness) is defined as the critical stress intensity factor K Ic at which an atomically sharp crack would start to propagate under mode I loading [7]. In the case of negligible plastic deformation, the theoretical value K G is related to the critical energy release rate G c , i.e., the change in elastic strain energy per unit area of crack advance, by where E denotes an appropriate elastic modulus [8]. For completely brittle fracture of a single crystal (SC) the critical energy release rate G c is according to Griffith [9] equal to the surface energy 2γ surf of the two newly created surfaces of unit area: In case the crack propagates along a grain boundary (GB) between two grains 1 and 2, the energy required to create the two surfaces, which can be crystallographically different, has to be considered relative to the GB energy γ GB [7,10]: In brittle materials, fracture along symmetrical GBs should therefore always be more favorable than cleavage on a plane parallel to the GB (γ 1 surf = γ 2 surf = γ surf ). Whether fracture in a brittle polycrystal takes place by transgranular cleavage or intergranular fracture however depends, besides geometrical factors, on the ratio of G GB on oxide ceramics have hinted at the importance of bond trapping in the selection of the transgranular mode of fracture over intergranular fracture [12,19], a systematic study of the influence of bond trapping on GB fracture is currently lacking.
Knowledge about the fracture toughness of GBs is of particular importance for continuum models of intergranular fracture [20][21][22][23][24][25][26][27][28][29][30][31][32], which allow to investigate failure of polycrystalline materials at the macroscale as well as for grain boundary engineering of structural materials [33,34]. The experimental determination of the fracture toughness of individual GBs is however challenging, as it requires the preparation of notched flat and defect-free GBs with defined misorientation. Therefore, only few measurements of GB fracture toughness exist [35,36]. Most experiments instead determine the fracture strength of GBs by measuring the tensile stress at which a GB fails in an uniaxial tensile test. The fracture strength of GBs was however investigated that way only for relatively few materials, e.g., high-purity Mo [37][38][39][40], Ni 3 Al [41,42], Ni-20Cr [43] and Si [44]. In the context of the present work, it is interesting to note that these experiments often show that samples with GBs normal to the tensile direction can have fracture stresses of the order of single crystals [33,39,40,43] or higher [38,41,42,44] and even show intergranular instead of transgranular fracture [40][41][42][43]. In most cases the strong GBs are Σ3, Σ9, Σ17 coincidence site lattice (CSL) GBs [38,[41][42][43], but there are exceptions [40,44]. In addition, low Σ GBs like Σ5 and Σ7 showed lower fracture stresses than the strong GBs and intergranular fracture [41], and the GB fracture stress is not always correlated to the GB energy [42,43].
The determination of the fracture stress in tensile tests involves plastic de-formation, which in turn depends on the crystallographic orientation of both grains. Such experiments therefore cannot provide direct information on the GB strength. However, these experiments show that GBs are not per-se weak and that their strength is not directly predictable [41].
Here we report on atomistic calculations of grain boundary fracture of a series of tilt GBs in tungsten, focusing on the quantitative determination of the fracture toughness and its direction and position dependence.

Methods
In order to study bond trapping in GB fracture, we performed molecular statics (MS) simulations of mode I cracks in five symmetrical and one asymmetrical tilt GBs under variation of their reciprocal density of coincidence sites, Σ, [10], and on single crystals of corresponding orientations.
Moreover, tungsten's near-perfect elastic isotropy simplifies the simulation set-up and analysis.
The atomic interaction is modeled by a Finnis-Sinclair potential [62,63] which has been shown to represent well the experimentally observed fracture behavior of tungsten [15]. The potential properties relevant for this study are summarized in in Tab. 1.

Construction of grain boundaries
The GBs are constructed following the approach detailed in [64]. Two crystal grains are rotated with respect to each other according to the rotation axes and angles provided in Tab. 2 and joined along the GB plane. Periodic boundary conditions (PBC) are used for the x− and z-directions within the GB plane whereas free boundary conditions are used along the GB normal (=y-direction). The GB structure is optimized with respect to the microscopic degrees of freedom by rigidly translating the two crystals with respect to each other and subsequently minimizing the energy while only allowing for motion of the atoms in direction of the GB normal. Using the microscopic translations of the two crystals which lead to the structure with the minimal potential energy as initial condition, the GB is fully relaxed without restriction of atomic motion and allowing for changes of the box size in the x-and z-direction to obtain a stress-free structure.
The exact size of the simulation box in the x-and y-directions is determined by the periodicity of the coincidence site lattice. The width in the z-direction is chosen such that it is at least three times the cut-off radius r cut of the potential, the extension of the GB in y-direction is chosen such that it is larger than 30 nm, as is the sample size in the x-direction.
Throughout this study the FIRE algorithm [65] is used to minimize the potential energy of a system. Samples are considered to be in a minimum energy configuration when the force norm per degree of freedom is smaller than 10 −8 eV/Å. The so obtained GB structures and energies agree well with the literature, see e.g. [48,66].

Simulation set-up
The simulation set-up is schematically shown in Fig. 1. From the relaxed configurations, cylindrical samples with radius R = 15 nm which contain the GB in the middle are cut out. Atoms within a distance R − 2r cut ≤ r ≤ R from the center of the cylinder are held fixed. An atomically sharp crack is introduced at the center by displacing the atoms according to the displacement field u xy (r, θ, K) around a crack tip in an anisotropic body obtained from LEFM under plane strain conditions [8]. For more details on the introduction of cracks see [67]. PBC are used along the cylinder axis which corresponds to the crack front direction as well as the GB tilt axis.
In addition, cracks with the same crack plane and crack front direction were studied using the set-up in SCs with crystallographic orientations corresponding to the rotated grains.

Determination of critical stress intensity factors
The determination of the critical stress intensity factor at which the crack either propagates or crack tip plasticity is initiated requires a well-defined initial structure. I.e., upon energy minimization, the crack tip has to remain stable and at the exact center of the cylindrical set-up. The stress intensity factor at which this is the case, K in , is determined by displacing the atoms according to the displacement field u xy (r, θ, K test ) for various test values K test and subsequently minimizing the energy of the system. The values for K test are chosen to be slightly below the stress intensity factor K G according to the Griffith criterion for the given GB or SC (Eqs. (1-3)). The initial stress intensity factor K in then corresponds to the value of K test at which the crack tip remains stable at the initial position at the center of the cylinder.
From the initial stress intensity K in the crack is loaded by incrementally displacing the atoms according to u xy (r, θ, ∆K I ), which corresponds to an increase in load of ∆K I = 0.007 MPa √ m. The crack structure is subsequently relaxed with the FIRE algorithm, and if no significant change of the bond lengths at the crack tip or of their configuration is detected, the relaxed structure is loaded by a applying a further displacement field increment.
The stress intensity factor at which the bond configuration at the crack tip changes either by brittle fracture or by crack tip plasticity is in the following referred to as K Ic . For more details on the iterative loading procedure and determination of K Ic see [68].
For brittle propagating cracks, K + is identical to K Ic . K − is determined by incremental unloading of the crack tip, i.e., by the application of a displacement field u xy (r, θ, ∆K I ) with ∆K I = -0.007 MPa √ m to the atomistic configuration derived from the previous unloading step. Finally, the stress intensity factor where the crack closes is denoted as K − .
The influence of the boundary conditions were tested by performing simulations on selected GB cracks using four times larger cylinder radii. The observed fracture mechanisms were identical to the ones in the larger set-up, and the values of K Ic varied only within 2 − 3∆K I , the resolution of the stepwise iterative loading procedure.

Cracks in single crystals
Nine different combinations of crack plane and crack front directions were studied in the SC set-up. The theoretical fracture toughness according to  [15]. This value exceeds K SC G , which has been attributed to lattice trapping [55]. The lattice trapping in this case was determined to ∆K SC =0.17, see Tab. 4. Similarly, all other cleavage cracks require higher loads than K SC G to propagate.
With the exception of the (710)[001] crack, which kinked out onto the (110) low-energy plane, all brittle cracks propagated on the initial cleavage plane. An example for brittle crack propagation is shown in Fig. 3

Grain boundary cracks
The grain boundary energies γ GB of the six studied GBs are summarized in Tab. 2, together with the corresponding surface energies γ surf , and the theo- . This table also includes the information of the type of event at K GB Ic , i.e., whether the crack propagated by cleavage or plastic deformation was observed. Note that GB fracture was always studied using at least two opposite crack propagation directions within the GB plane. The direction dependence of GB fracture is further exemplified in Fig. 2b) where the fracture toughness for ten different crack front directions in the Σ3 (111) (112) crack plane. This difference in fracture surfaces is related to the bond breaking process. Fig. 3 shows that the SC crack advances by successive breaking of single bonds whereas quasistatic crack propagation within the GB takes place by simultaneous breaking of two different bonds.
The atomic crack tip configuration for this crack system is shown in more detail in Fig. 4a,b). As expected, the bonding situation is clearly different in the GB compared to the SC. In particular, the two atoms forming the critical bond at the crack tip see a different environment: in the grain boundary, the atoms directly at the crack tip interact with more atoms and have a lower potential energy as compared to the single crystal (see Fig. 4). Following [16], the breaking of atomic bonds at the crack tip is determined by monitoring the bond distances across the crack plane for the two bonds indicated in Figs. 4a,b) as function of the applied load, Fig. 4c). Upon loading, the bond distances remain about constant up to a critical load K + > K G where the bond distance suddenly increases. When the crack is unloaded from this configuration, the distance between the now separated atoms decreases nearly linearly, until the bond forms again at K − which is lower than K G . Further unloading leads to the formation of a new bond with an atom behind the original location of the crack tip (dashed lines in Fig. 4).
In the SC case, the breaking of the first bond (marked by arrows in

Discussion
The usual thermodynamic description of brittle GB fracture in a continuum mechanics framework implies that the fracture toughness of GBs is lower than the fracture toughness of a perfect single crystal in the same orientation, and that the fracture toughness does not depend on the crack tip position within the GB or on the crack propagation direction. Our simulation results show clearly, that at the atomic scale these assumptions can no longer be used.
Direction dependent fracture behavior of GBs has been already reported for experiments [69] and simulations [50] on Cu bicrystals. In these cases, the direction dependence can be explained as a direct consequence of the Rice-Thomson model [69][70][71]: for GB fracture, the angles between the crack propagation direction and the slip systems in the two grains are different for a crack propagating along the +[hkl]-direction from a crack propagating in the −[hkl]-direction. I.e. in one direction emission of dislocations or twins is easier than in the other. Our simulations however show that also for purely brittle fracture K Ic can depend on the crack propagation direction. Fig. 2a)  have significantly higher fracture toughness than the SC cracks on the same crack plane (see also Video 1 and 3 in Supplemental Material at [URL will be inserted by publisher]). It is important to note that also in these cases, fracture takes place in a perfectly brittle manner, and the differences in fracture toughness are not caused by plasticity. These GBs can also be tougher than perfect single crystals oriented for fracture on the primary cleavage plane.

Comparing in
A large number of such tough GBs could therefore favor transgranular over intergranular fracture. In this context it is important to note that the increased toughness, i.e., a significantly higher K Ic than K GB G , is not correlated with the GB energy. E.g., the Σ25 GB has a significantly higher fracture toughness than the lower energy Σ3 and Σ9 GBs.
The deviation of the fracture toughness from the thermodynamic prediction K G has to be rooted in the nature of the atomic bonds at these GBs. Locally, the required energy release rate for crack advance can however be larger in the GB than in the SC case.
Although GBs show regions with low fracture toughness, it is the region with the highest K Ic which determines ultimately the crack initiation toughness of a GB. This is clearly shown in Fig. 5c-f): under quasi-static loading, the crack tip located initially at the weakest spot of the GB just propagates by breaking one bond and is arrested at the next, stronger unit. Yet, if the crack tip is positioned at the strongest bond, the crack propagates through all GB units. This situation might be different for dynamically propagating cracks. The determination of the crack arrest toughness is, however, beyond the scope of this study.
Similar to studies on lattice trapping [16,55], the present study of the effect of bond trapping on grain boundary cracks is limited to static calculations in a quasi two dimensional set-up. Temperature effects and the presence of kinks on the crack front are neglected. However, the 2D quasi-static simulations on lattice trapping were able to explain macroscopically observed fracture behavior, e.g., the occurrence of primary {100} cleavage planes in tungsten [15], or the anisotropy with respect to the propagation direction along different crystallographic directions within the same fracture plane in silicon [16]. Similarly, bond trapping of GB cracks could have macroscopic implications e.g. on the selection of the transgranular mode of fracture over intergranular fracture as speculated in [12,19].
In general, controlled experiments on GB cracks are very difficult to perform [41]. The presence of GB defects like steps or absorbed dislocations, segregated atoms, and the often non-planar nature of GBs and curved crack fronts, as well as the presence of dislocations in the grains, make a direct comparison of our simulation results with experiments nearly impossible. Recent developments in nanomechanical testing however now allow the determination of the fracture toughness of interfaces from in-situ bending tests of notched cantilever beams [74]. Such tests could be performed with well-characterized GBs in polycrystalline materials, avoiding the difficulties in preparing bicrystal specimens. In our simulations, the difference of the GB fracture toughness for different crack propagation directions is rather small (see Fig. 2) and will therefore most probably not be detectable in experiments. However, according to our simulations, the determined fracture toughness can deviate from the thermodynamic prediction according to Eq. (3) by up to 90%, see Fig. 2 and Tab. 3. Comparing the experimentally determined GB fracture toughness in an inherently brittle material with the corresponding theoretical prediction or with a numerical calculation using a cohesive zone model (CZM) [20,75,76] could thus allow to reveal the influence of bond trapping on GB fracture toughness. In fact, bond trapping might be one of the underlying reasons for the reported extraordinary high tensile strengths of some GBs [38,41,42,44,61], in particular when they are not directly correlated with the GB energy [42,43] or otherwise not directly predictable [41].
Currently, GB fracture is typically modeled on the meso-and continuum scale using a cohesive zone approach [20][21][22][23][24][25][26][27][28][29][30][31][32]. Commonly used CZMs [20,75,76] require as parameters the ultimate/maximum tensile strength of the GB and the work of fracture needed to separate the cohesive zone. This information can be determined from atomistic [25,26,30,32] or ab-initio simulations [77]. However, the crack is often not directly modeled in these approaches, and instead these values are derived from tensile straining of GBs [25,77]. Such simulations are averaging over the atomic structure of GBs, and continuum simulations with CZM parameters derived that way failed to reproduce the atomistic simulation results [29]. By their very nature, simple CZMs are furthermore not able to properly describe the direction dependence of the fracture toughness of GBs. Here, theories like quantized fracture mechanics [78], which take into account lattice trapping, could be extended to include bond trapping of GB cracks. Ultimately, however, parameterized GB CZMs would need to be validated using fracture experiments on defined GBs. To our knowledge, such direct comparisons are currently still lacking.

Summary
The present study highlights the role of bond trapping effects on the frac- Excellence 'Engineering of Advanced Materials' at the Friedrich-Alexander-Universität Erlangen-Nürnberg (FAU) which is funded by the DFG within the framework of its 'Excellence Initiative'. Table 1: Summary of the properties of the Finnis-Sinclair potential [62,63] relevant to the simulation of fracture in tungsten in comparison to experimental and DFT data.
The elastic constants are denoted by C ij , and γ surf (hkl) are the surface energies of the corresponding (hkl) planes.

Property
Experimental data DFT data F-S Potential a Elastic constants at 300 K according to Bujard [79] b Surface energies for average orientation taken from Miedema [80] c Elastic constants calculated by Einarsdotter et al. [81] (using DFT with the local density approximation, a plane-wave basis set, and norm-conserving pseudopotentials) d Surface energies calculated by Vitos et al. [82] (using DFT with the generalized gradient approximation, full charge density linear muffin-tin orbitals method)  The crack deviates from its initial (710) crack plane onto low-energy planes.           in crystal lattices (red, according to [7]), GBs (blue) and the corresponding linear continuum limits. The periodicity of bond trapping corresponds to the repeat distance of the structural units, and the magnitude of the bond trapping can be significantly larger than the lattice trapping. The toughness of a GB can therefore locally surpass that of a crystal fracturing along the same plane. The magnitude of the bond trapping effect can furthermore show a backward-forward or direction dependence since the GB structure is generally asymmetric.