## I. INTRODUCTION

Resistance to crack propagation is undoubtedly one of the most important properties of structural materials. However, our current mechanistic understanding of the energy-dissipation processes involved on the various scales and of their interplay is not sufficiently advanced to predict the fracture behavior of materials ab initio.^{Reference Hutchinson and Evans1,Reference Needleman and Van der Giessen2} As fracture is ultimately caused by the breaking of atomic bonds,^{Reference Bitzek, Kermode and Gumbsch3} there is a high demand for atomistically informed fracture criteria. An overview of ongoing attempts to realize them in hierarchical and concurrent schemes can be found, e.g., in Refs. Reference Shiari and Miller4 and Reference Rabczuk5. Here we propose a hierarchical scheme that meets two common challenges: to extract information on the atomic scale in an efficient manner and to coarse-grain this information to obtain practical and meaningful failure criteria for a finite element (FE) implementation on the mesoscale. Our approach is demonstrated by making use of cohesive zones (CZs), but is not limited to this model.

A so-called cohesive zone model is a seemingly natural way to incorporate atomistic information in a continuum model of fracture. In a traditional CZ model, the fracture process zone (i.e., the region in which nonlinear processes occur) ahead of the crack tip is represented by cohesive surfaces (respectively a narrow band of vanishing thickness) which experience cohesive tractions described by a constitutive law, the so-called traction-separation (TS) law. Crack growth occurs when the separation at the tail of the CZ, i.e., at the physical crack tip, reaches a critical value at which the cohesive traction vanishes. The strip yield model of Dugdale,^{Reference Dugdale6} which was originally proposed for estimating the size of the plastic zone around a crack tip, is regarded as one of the first CZ-type models. Barenblatt^{Reference Barenblatt7} also proposed a cohesive fracture concept with the aim to eliminate the crack tip stress singularity in classical linear elastic fracture mechanics (LEFM). For fracture of perfectly brittle materials, the descriptions by LEFM and the CZ model approach are equivalent if the material behavior in the zone around the crack tip is determined by the stress-intensity (*K*) field and if the cohesive energy density is identical to the critical energy release rate.^{Reference Rice and Liebowitz8} The first condition is met when the fracture process zone is smaller than the *K*-dominated zone (the region with a significant gradient in the stress), which is the case in brittle materials.^{Reference Singh, Kermode, De Vita and Zimmerman9,Reference Shimada, Ouchi, Chihara and Kitamura10}

A crucial point to consider in practical applications within a finite element method (FEM) implementation is that CZ-based numerical techniques require cracks to propagate along the element boundaries. Therefore, an inherent nonvanishing amount of mesh dependency arises, which can become dominant for three-dimensional simulations in a homogeneous material.^{Reference Zhou and Molinari11} For instance, no converged solutions from the simulation of progressive delamination were obtained above a critical mesh size.^{Reference Turon, Dávila, Camanho and Costa12} As a general rule, Hillerborg et al.^{Reference Hillerborg, Modéer and Petersson13} suggested a maximum characteristic length ${l_{{\rm{coh}}}} \approx {{{G_{\rm{c}}}} / {\sigma _{\rm{c}}^2}}$ for the CZ under mode I loading, based on the cohesive energy *G* _{c} and the critical stress to initiate crack advance σ_{c}. This characteristic length is always smaller than the radius of the *K*-dominated zone at critical loading, but also larger than the fracture process zone. In our approach, we couple the mesh size to the critical stress, which is the object of a scaling procedure explained in Sec. III. Furthermore, we do not only treat the fracture process zone as a CZ, but sample the complete *K*-dominated region on the cleavage plane with cohesive elements. This way, we obtain mesh-independent results for the fracture toughness under mode I loading.

In principle, the TS law can capture the physics of deformation and failure on the microscale, provided its shape is chosen carefully and that it includes important parameters such as maximum strength, work of separation, and critical displacement.^{Reference Chandra, Li, Shet and Ghonem14} Within the hierarchical modeling approaches, there are basically two philosophies about how to incorporate this information in TS laws, which differ in the way in which dissipative processes enter the modeling framework.

TS laws derived from atomistic simulations of homogeneous deformation of a single crystal, a bicrystal, or even model sections of polycrystalline microstructures can produce a stress–strain curve which reflects both elastic and dissipative processes during fracture.^{Reference Spearot, Jacob and McDowell15–Reference Barrows, Dingreville and Spearot21} However, dissipating the energy which is released during crack initiation and growth can be severely constrained in microstructures^{Reference Kumar and Curtin22} so that the fracture toughness becomes dependent on length scale and geometry. This means that the TS law would have to be adapted for each specific microstructure. Instead, we want to promote an approach in which the process of material separation by breaking of atomic bonds is separated from energy dissipative processes, such as plastic deformation accommodated by generation and motion of dislocations. The description of bond breaking within the CZ can be combined with different levels of description of the plastic deformation in the material adjacent to the CZ. On the smallest scale, plasticity is described based on discrete dislocation dynamics models in two^{Reference Cleveringa, Van der Giessen and Needleman23–Reference Broedling, Hartmaier and Gao25} or three dimensions.^{Reference Déprés, Prasad Reddy, Robertson and Fivel26} Further developments include the formulation of the virtual internal bond model of cohesive fracture and its application to nanomaterials.^{Reference Zhang, Klein, Huang, Gao and Wu27,Reference Ji and Gao28}

Following the failure criterion by Griffith, the TS law for the bond-breaking process will have to capture the work of separation of the material, but to predict reliably the trends in cohesive behavior, the strength of the bonds is required as well.^{Reference Tahir, Janisch and Hartmaier29} The strength of interatomic bonds is reflected in the theoretical strength of a material, which can be determined via ab initio density functional theory (DFT) calculations for single crystals^{Reference Ogata, Li and Yip30,Reference Černý and Pokluda31} as well as interfaces with and without impurities.^{Reference Šob, Friák, Legut, Fiala and Vitek32–Reference Gibson and Schuh34} However, this procedure assumes homogeneous expansion perpendicular to, or cleavage along, an infinite crystallographic plane. During such a process, all bonds across the anticipated cleavage plane are elongated by the same amount and rupture at the same time. By contrast, individual pairs of atoms across the cleavage plane in the vicinity of the crack tip experience different environments, due to the gradient of the elastic displacement ahead of a crack tip, as shown in Fig. 1. Thus, the equivalence of the stresses obtained from such a rigid body separation (RBS)-type calculation with the critical stress needed to advance a crack tip still has to be demonstrated. This is one of the aims of the present paper (see Sec. II).

Embedding an atomistics-based TS law into a fracture model via CZs in a FE scheme has also practical issues because the theoretical strength is of the order of GPa and the critical separation is of the order of Å.^{Reference Tahir, Janisch and Hartmaier35} The stresses on the order of the theoretical strength occur very locally in the extreme stress gradient in front of the crack tip, i.e., usually within the first few layers of atoms. In contrast to this, any continuum modeling is based on describing average stresses in some finite volume such that typical stresses occurring in FE models are much smaller. Therefore, Serebrinsky et al. first used a scaling relation for TS laws^{Reference Serebrinsky, Carter and Ortiz36} derived from the results of ab initio calculations based on the concepts introduced in Refs. Reference Nguyen and Ortiz37 and Reference Hayes, Ortiz and Carter38. It assumes that the critical stress needed to separate two half crystals depends on the amount of elastic energy that is stored and released from the material next to the cleavage plane, which in turn scales with the system size. This can easily be observed in atomistic supercell calculations in which two half crystals are shifted and the atomic positions are relaxed afterwards. If the theoretical strength is simply identified with the maximum of the derivative of the energy–displacement curve, it decreases with an increasing number of planes parallel to the cleavage plane.^{Reference Nguyen and Ortiz37} However, the amount of energy that localizes and is subsequently released during the fracture process is in fact finite,^{Reference Elsner and Müller39} and the seemingly size-dependent fracture stress can be normalized to obtain the unique strength of the material.^{Reference Enrique and Van der Ven40} Thus, a scaling approach via the elastic energy is practical, but not well justified, and a physically meaningful coarse-graining scheme providing an unambiguous connection between the macroscopic and microscopic stresses is still needed.

With the work at hand, we want to promote tensile tests carried out by RBSs and translations in combination with a scaling approach based on pairwise forces. RBS schemes can be used to characterize single crystals and grain boundaries on the same footing,^{Reference Janisch, Ahmed and Hartmaier41} provide a three-dimensional description of single crystal^{Reference Sun, Beltz and Rice42} and grain boundary separation,^{Reference Pang, Janisch and Hartmaier43} and can be extended to include impurity effects.^{Reference Tahir, Janisch and Hartmaier35} Most importantly, the derived critical stresses are size independent without any normalization. Our scaling approach based on interatomic forces will allow us to relate this fundamental atomistic quantity to the required mesh-sensitive model parameters in a rigorous way.

In the following, the results of atomistic simulations of RBS are compared with those of simulations of *K*-dominated cracks to demonstrate the equivalence of the interatomic forces in both situations; see Sec. II. This is the starting point for our novel approach to scale the atomic-level properties into parameters for continuum models in a physically motivated way, which is presented in Sec. III. This novel approach is tested in the framework of CZ modeling of cleavage fracture in an FE model, which demonstrates its advantages and shortcomings. Our findings are summarized in Sec. IV.

## II. COMPARISON OF INTERATOMIC FORCES AT CRACK TIPS AND DURING RIGID-BODY SEPARATION

### A. Interatomic potentials and simulation setups

We chose tungsten as a brittle, elastically isotropic model material to study cleavage fracture under static loading conditions. The atomic interactions were described either in the framework of ab initio DFT calculations using the local density approximation (LDA) or by using different potentials of the embedded atom method (EAM) type. DFT calculations were carried out using the ABINIT electronic structure code,^{Reference Gonze, Beuken, Caracas, Detraux, Fuchs, Rignanese, Sindic, Verstraete, Zerah, Jollet, Torrent, Roy, Mikami, Ghosez, Raty and Allan44–Reference Gonze, Amadon, Anglade, Beuken, Bottin, Boulanger, Bruneval, Caliste, Caracas, Côté, Deutsch, Genovese, Ghosez, Giantomassi, Goedecker, Hamann, Hermet, Jollet, Jomard, Leroux, Mancini, Mazevet, Oliveira, Onida, Pouillon, Rangel, Rignanese, Sangalli, Shaltaf, Torrent, Verstraete, Zerah and Zwanziger46} employing the ultra-soft pseudopotential of Hartwigsen, Goedecker and Hutter.^{Reference Hartwigsen, Goedecker and Hutter47} For the RBSs described below, a *k*-point mesh of 8 × 8 × 2 *k*-points was used, which, after Fourier transformation, corresponds to a maximum *k*-point spacing of 47.45 a.u. (25 Å). The energy cut-off for the plane-wave basis set was 16 Ha (435 eV).

The EAM potentials used are the classical Ackland–Thetford–Finnis–Sinclair (ATFS) potential from the 1980s^{Reference Finnis and Sinclair48,Reference Ackland and Thetford49} and the more recent potential by Wang.^{Reference Wang, Zhou, Li and Hou50} Their fracture-relevant properties and cutoff radii *r* _{cut} are presented in Table I. All atomistic simulations were performed with the molecular dynamics software package IMD^{Reference Stadler, Mikulla and Trebin51,52} using the MIK^{Reference Beeler53} and FIRE^{Reference Bitzek, Koskinen, Gähler, Moseler and Gumbsch54} relaxators.

The simulation setups for the corresponding atomistic simulations are shown in Fig. 1. Figure 1(a) shows the case of rigid-body separations, which is often used to determine the theoretical strength via ab initio tensile tests. In this setup, two crystal halves are iteratively separated without allowing the atoms to relax their positions. As mentioned in the introduction, this strategy prevents size effects, which would occur if the crystallographic layers were allowed to relax, and furthermore enables a unified description of single crystals and grain boundaries.^{Reference Janisch, Ahmed and Hartmaier41} Note, however, that the strength values determined from methods that include atomic relaxations and lateral contraction of the sample can differ significantly, particularly in the presence of defects, such as grain boundaries.^{Reference Černý, Šesták, Řehák, Všianská and Šob55} The implications will be discussed below. Periodic boundary conditions (PBC) are used in all directions. In the case of a grain boundary calculation, this actually requires the inclusion of two cleavage planes as there would be two grain boundaries in the simulation cell. This is not necessary in the single crystal, but to be able to use the same setup, we did so for consistency.

The size of the simulation box in the present study was *L* _{x} × *L* _{y} × *L* _{z} ≈ 10 × 60 × 10 Å^{3}. The separation distance Δ_{y} was varied from −0.5 Å to 2.5 Å in incremental steps of dΔ_{y} = 0.025 Å. The normal stress σ_{at} acting between the atoms of the two crystal halves at a certain separation distance Δ_{y} was then calculated from the average potential energy per atom, *E* _{pot}, as

where *N* _{atoms} is the number of atoms in the configuration and *A* _{xz} is the cross-sectional area of the created surface. The factor of 2 results from the PBC.

Figure 1(b) shows the setup used to determine the fracture toughness and structure of an atomistically sharp crack. Although this setup can also be used to study cracks with DFT,^{Reference Perez and Gumbsch56} it is usually used in classical atomistic simulations with interatomic potentials. In this setup, the linear-elastic solution of the crack-tip displacement field under plane-strain conditions, *u* _{x,y}(*r*, θ) is directly applied to a cylindrical atomistic configuration with PBC.^{Reference Sih, Paris and Irwin57,Reference Sih, Liebowitz and Liebowitz58} The displacements are determined for a given stress intensity factor *K* _{I} and depend on the distance from the central crack tip, *r*, and the inclination angle θ of the direction of *r* and the crack plane. The critical stress intensity factor for the onset of crack propagation, *K* _{Ic}, was determined by iteratively increasing *K* _{I}, the corresponding Δ*u* _{x,y}, and the subsequent relaxation of the atoms in the interior of the configuration (with a minimum distance of 2*r* _{cut} from the outer surface). The value of *K* _{I} at which the crack propagates by brittle bond breaking was then defined as *K* _{Ic} (we will denote this quantity as $K_{{\rm{Ic}}}^{{\rm{at}}}$ hereafter to emphasize that it was obtained by atomistic simulations). For more details on the setup and simulation methodology, see Refs. Reference Möller and Bitzek59 and Reference Möller and Bitzek60. In the present study, a cylinder of radius 150 Å and length *L* _{z} ≈ 10 Å was used to study a crack on the (100) crack plane along the crack-front direction [011], which is the natural crack system in tungsten.^{Reference Gumbsch, Riedle, Hartmaier and Fischmeister61}

### B. Results

#### 1. Simulations of rigid-body separation

The dependence of the binding energy per unit area, i.e., the energy difference of a configuration with respect to the fully separated crystal halves, on the separation distance is shown in Fig. 2 for the two interatomic potentials in comparison to DFT results. From this, the cleavage energy *G* _{c} is obtained as the difference between the energy at infinite and that at zero separation. The values are summarized in Table II. By comparing with Table I, it can be seen that the cleavage energy approximately equals 2γ_{s}(100), which would be the (relaxed) work of separation. The work of separation (WoS) obtained from the DFT calculations (not shown in Table I) equals 8.58 J/m^{2}, which is in good agreement with other DFT values.^{Reference Vitos, Ruban, Skriver and Kollar62,Reference Piazza, Ajmalghan, Ferro and Kolasinski63} The values for the cleavage energies obtained with the EAM potentials are significantly different from the DFT data. Such deviations are not uncommon: see, e.g., Ref. Reference Möller, Mrovec, Bleskov, Neugebauer, Hammerschmidt, Drautz, Elsässer, Hickel and Bitzek64 for a review of the applicability of various classes of potentials for bcc transition metals to fracture-like situations. However, for the aim of the present study, it is sufficient to show that the atomistically determined fracture toughness can be reproduced independently of the mesh size of our CZ model. In addition, we will show how the DFT results from the simpler RBS calculations can be used in the same approach, opening up the possibility of a more accurate determination of *K* _{Ic} without performing the actual DFT crack-tip simulations.

The TS relationships for the different potentials were readily obtained by applying Eq. (1), i.e., by derivation of the binding-energy versus separation-distance data. The corresponding TS plots are shown in Fig. 2(b). Three characteristic values can be obtained from these plots: the threshold stress, σ_{th}, the critical separation distance δ_{crit} (where σ_{th} is reached), and the ‘final’ distance δ_{f} (where the traction across the surface is zero). All characteristic values are summarized in Table II and compared to values obtained by DFT calculations.

Based on the experience with previous studies on RBS,^{Reference Tahir, Janisch and Hartmaier35,Reference Janisch, Ahmed and Hartmaier41} the shape of the TS relationships of the Wang potentials appears to be somewhat unexpected as it shows two local maxima and a local minimum at around Δ_{y} ≈ 0.65 Å. The origin lies in the not-completely-smooth energy–displacement curves [see Fig. 2(a)]. Both potentials were optimized to reproduce several material properties, but an interplanar separation along the [001] axis was not included. However, as we will see in Sec. III, it is sufficient if the TS law includes the peak stress, and this is reasonably well reproduced by both potentials.

#### 2. Crack-tip simulations

The fracture toughness $\left( {K_{{\rm{Ic}}}^{{\rm{at}}}} \right)$ values determined by atomistic simulations are summarized in Table II. According to the Griffith criterion^{Reference Griffith65} ${K_{\rm{G}}} = \sqrt {2{\gamma _{\rm{s}}}\left( {100} \right) \cdot E}$ (with *E* being the appropriate orientation-dependent elastic constant), critical values in the range of 1.61–1.64 MPa m^{1/2} are expected for both potentials. These theoretically expected values using the potential parameters are markedly smaller than the atomistically determined ones; see Sec. II. This well-known behavior^{Reference Gumbsch66} is generally attributed to the lattice-trapping effect.^{Reference Thomson, Hsieh and Rana67} The atomistic crack-tip configurations showing crack propagation by cleavage on the original (100) plane are provided in Fig. 3.

#### 3. Comparison of interatomic forces

There are different ways to calculate the force acting across the cleavage plane from the present simulation results: $F_{{\rm{ct}}}^{{\rm{at}}}$ is the *y*-component of the individual atomic (superscript ‘at’) force vectors calculated during the simulation of crack-containing configurations (subscript ‘ct’). To compare the evolution of the crack-tip forces with the forces during rigid-body separations, the crack-tip distance *r* _{ij} has to be related to the corresponding separation in the [100] direction, Δ_{y}, as follows:

where *a* _{0} denotes the bulk lattice parameter. Strictly, this correction is only valid under the assumption that the distance between next-nearest (NN) neighbors (and for cracks on the (100) plane NN neighbor bonds are the critical bonds) changes with Δ_{y}, as follows:

i.e., when applying a homogeneous uniaxial strain without allowing Poisson contraction, as is the case in the present study.

$F_{{\rm{rb}}}^{{\rm{at}}}$ is the *y*-component of the individual atomic force vectors calculated during the simulation of rigid-body separations (subscript ‘rb’). Since there are *N* _{surface} interacting NN pairs and ${{{N_{{\rm{surface}}}}} \over 2}$ 2nd NN pairs for each (100) surface [Fig. 1(a)], the total force $F_{{\rm{rb}}}^{{\rm{at}}}$ per surface atom is obtained from the individual forces between NN and 2nd NN neighbors, $F_{{\rm{NN}}}^{{\rm{at}}}$ and $F_{2{\rm{NN}}}^{{\rm{at}}}$, by

$F_{{\rm{rb}}}^{{\rm{avg}}}$ is calculated from the change of the average potential energy per atom in rigid-body separations:

which corresponds to the calculation of the tractions in Eq. (1) but normalized to the number of surface atoms rather than the surface area.

*F* ^{eff} is the derivative of the effective pair potential.^{Reference Johnson68} In this case, the same correction regarding *r* _{ij} and Δ_{y} has to be made as for the crack-tip forces; see Eq. (2).

The evolution of these forces as a function of lateral separation Δ_{y} is shown in Figs. 4(a) and 4(b) for the two potentials. The truly remarkable result of this comparison is how well the force laws from the by-far-simpler model systems (rigid-body separations and effective pair potential) agree with the actual forces at the crack tips, considering that the environments in which the forces have been determined (or calculated as in the case of *F* ^{eff}) are completely different.

Second, it should be noted that the unexpected camel-hump shape of the *T*(Δ_{y}) curve for the Wang potential shown in Fig. 2(b), which is essentially the per-area instead of per-atom versions of $F_{{\rm{rb}}}^{{\rm{avg}}}\left( {{\Delta _y}} \right)$ and thereby $F_{{\rm{rb}}}^{{\rm{at}}}\left( {{\Delta _y}} \right)$, can now be understood by the separation of the total acting force into contributions from NN and 2nd NN neighbors.

Third, the individual force–separation relationships ($F_{{\rm{rb}},{\rm{NN}}}^{{\rm{at}}}$ and $F_{{\rm{rb}},2{\rm{NN}}}^{{\rm{at}}}$) are very similar to each other, and their shape can be qualitatively understood by comparing it to the shape of the derivative of the effective pair potential, *F* ^{eff}. This implicitly shows that the correction made in Eq. (2) is applicable even if it does not account for changes in the electron density due to uniaxial straining, which might, however, be the reason for the minor differences in the magnitudes of $F_{{\rm{rb}},{\rm{NN}}}^{{\rm{at}}}$ and *F* ^{eff}.

An important conclusion which can be drawn from these results is that the forces necessary to break the bond between the crack-tip atoms show only a small—if any at all – dependence on the local atomic surrounding at the crack tip. Even if our observations cannot prove that all EAM potentials will show the same behavior or that crack-tip bonds in metals in general do not show a strong dependence on their local bonding environment, this is to our knowledge the first time that it has been shown that the forces to separate crack-tip atoms can indeed be estimated from simulations of rigid-body displacements. In particular, the comparison between $F_{{\rm{ct}}}^{{\rm{at}}}$ and $F_{{\rm{rb}}}^{{\rm{avg}}}$ in Fig. 4 shows that the characteristics of TS laws derived from simple rigid-body separations, i.e., the maximum force and the critical separation distance δ_{crit}, are in very good agreement with the actual crack-tip forces. This has important consequences for multiscale modeling, as it validates the long-standing use of DFT calculations to determine traction separation laws. Whether the correspondence of forces from rigid-body separation and fracture simulations still holds for systems with complex chemical compositions or grain- and interphase boundaries would still need to be shown. However, an averaging of forces over one structural or chemical unit would be a possible way to extend our approach.

## III. ATOMISTICALLY INFORMED TRACTION SEPARATION (TS) LAW

### A. Derivation of a scalable TS law

In this section, a method is derived to convert the atomic forces at the crack tip into meaningful continuum quantities, i.e., TS laws describing the cohesive behavior of the material in front of the crack tip. To accomplish this, we consider a chain of *N* pairs of atoms in front of the crack tip on the crack plane (see Fig. 5). The atom pairs under consideration are those that are affected by the stress concentration of the crack tip and that lie in the crack plane.

Formally, the atomic forces *F* ^{(i)}, defined in Sec. II.B.2, can be converted into atomic stresses by normalization with the atomic area *A* _{at} as

where *x* ^{(i)} = *ia* _{cd} is the distance of the atom pair *i* = 1, 2, …, *N* to the crack tip and *A* _{at} = *a* _{cf}*a* _{cd} is the projected atomic area on the crack plane, given by the interplanar distance *a* _{cf} along the crack front and the interplanar distance *a* _{cd} in crack direction, respectively. The crack itself is defined by the pairs of atoms that are separated by a distance of more than the critical displacement δ_{f}, such that the position of the crack tip is located at the last pair of atoms with a larger separation than δ_{f}. It is noted here that the elastic distortion of the atomic distances along the crack front and in the crack direction is neglected to obtain a closed-form solution for the scaling relation.

From LEFM, it is known that the stress field in front of a crack tip decays as ${1 / {\sqrt {{x^{\left( i \right)}}} }}$. Figure 6(a) shows that this linear-elastic theory holds excellently at the atomic scale, even very close to the crack tip, in agreement with the findings of Ref. Reference Singh, Kermode, De Vita and Zimmerman9. The deviation for the crack-tip atom is partly due to the 20% increased Voronoi volume compared to the atoms in front of the crack tip. In contrast to the atomic stresses, the atomic forces are largest for the crack-tip atom pair. The atomic stresses along the cleavage plane in front of the crack tip can be reasonably expected to follow

with the stress intensity factor

fully defined by atomistic quantities.

The fundamental relation between the atomistic and the continuum description is demonstrated by calculating the fracture toughness *K* _{Ic} from Eq. (8), as

with the critical force *F* _{c} at which the bond between the first pair of atoms breaks; i.e., this pair of atoms is separated beyond the critical distance δ_{f}, and consequently the crack tip advanced by one atomic distance. This force corresponds to the theoretical strength defined in Sec. II.B.1 and can thus be conveniently calculated with atomistic methods. With the values obtained from atomistic calculations with the Wang potential (see Table II), a fracture toughness of *K* _{Ic} = 1.70 MPa m^{1/2} is calculated from the critical force, $F_{{\rm{rb}},{\rm{max}}}^{{\rm{avg}}}$, and *K* _{Ic} = 1.84 MPa m^{1/2} from the theoretical strength σ_{th} with ${a_{{\rm{cd}}}} = {{{a_0}} / {\sqrt 2 }}$. The force-based value corresponds very well to the *K* _{Ic} values obtained from direct atomistic simulation ($K_{{\rm{Ic}}}^{{\rm{at}}} = 1.68$ MPa m^{1/2}) and from the Griffith criterion ($K_{{\rm{Ic}}}^{\rm{G}} = 1.64$ MPa m^{1/2}). The value from the critical strength is roughly 8% higher as can be expected from the difference in $F_{{\rm{eff}}}^{{\rm{max}}}$ and $F_{{\rm{rb}},{\rm{max}}}^{{\rm{avg}}}$ (see Fig. 4).

However, for most continuum fracture models, a fracture energy or WoS is required in addition to the fracture strength ${\tilde{\sigma }_{\rm{c}}}$ and a critical separation ${\tilde{\delta }_{\rm{c}}}$. These parameters represent the main quantities defining the TS law describing the cohesive behavior of a material in front of a crack tip [see Fig. 6(b)]. While the WoS is scale independent, the fracture strength and the critical separation both depend on the spatial resolution of the applied method. The fundamental material constants are obtained from atomistic methods, but a proper scaling must be developed for the scale-dependent model parameters that takes an appropriate length scale or mesh size into account, as outlined in the introduction.

Since all values used in continuum methods represent average values over a finite volume or a FE, we start by defining the critical fracture strength as average stress over the *K* _{I}-dominated region

as the sum over all atomic forces in this region normalized by the corresponding area. From relation Eq. (7), we also conclude that the atomic forces decay as

With Eq. (11), the relation *A* = *NA* _{at} and the condition that *F* ^{(1)} = *F* _{c}, because the onset of crack advance is considered, we obtain

The sum on the right hand side of this equation can be approximated by

Numerical evaluation of the approximation shows that the relative error is less than 10^{−4} for values of *N* between 10^{5} and 10^{6}; for the latter value, the sum amounts to $\tilde{S}\left( N \right) = 0.002$. The precision of the approximation gets better for larger values of *N*, which is the usual case for applications. Since the sum is finite, a precise numerical evaluation is always possible, although the simple approximation function is still easier to handle.

To calculate a properly scaled final separation of atoms ${\tilde{\delta }_{\rm{f}}}$, we consider a simple bilinear TS law [see Fig. 7(b)] such that $\Gamma = {1 / {2\sigma _{\rm{c}}^{{\rm{at}}}}}\delta _{\rm{f}}^{{\rm{at}}} = {1 / 2}{\tilde{\sigma }_{\rm{c}}}{\tilde{\delta }_{\rm{f}}}$. It follows

In the final step, a reasonable estimate for the length of the *K*-dominated region and thus the number *N* of atom pairs to consider is developed. The criterion for the length of this region is the decay of the stress to a value of σ_{0}, which can be either the yield strength of a plastic material or the far-field stress ${\sigma _0} = {{{K_{\rm{I}}}} / {\sqrt {\pi {a_{{\rm{crack}}}}} }}$ considered in linear-elastic fracture mechanics, where *a* _{crack} is the crack length. Hence, *N* is defined by the requirement that

For the critical condition, when the crack starts to propagate (*F* ^{(1)} = *F* _{c}), we obtain

Noting that the critical atomic stresses are typically a 1000 times larger than the yield strength or far-field stress, a reasonable estimate for the number of atom pairs in the *K*-dominated region is *N* _{c} ≈ 10^{6}. This yields a size of *K*-dominated zone *r* _{coh} = *N* _{c}*a* _{cd} ≈ 200 μm.

In summary, a scaling relation for the atomistic strengths and separations has been derived, in the form

with the scaling factor

For practical applications, *N* is replaced by the ratio of the size of each element to the lattice parameter in crack propagation direction, and Eq. (19) becomes

This results in a critical stress and a final displacement, which are both adjusted to the element size, as shown in Fig. 7(b).

It is noted that the directional dependence of the model parameters in Eq. (20) is artificially introduced by the scaling of the FE size with the atomic interplanar distance in the crack propagation direction. A proper directional dependence could be introduced by using a direction-sensitive energy release rate as WoS.

### B. Application and validation

As an application of the novel scaling relation, we conducted FE fracture simulations for a linear elastic material representing tungsten to calculate the fracture toughness *K* _{Ic}. The numerical simulations have been performed using the commercial FE software ABAQUS. All simulations have been performed using an implicit time integration scheme. The surface-based CZ method is used to simulate fracture. The geometry of the model is the same as in the atomistic simulation [see Fig. 1(b)], whereas the ratio between the radius of the sample and the mesh at the crack tip is kept constant at 10,000. This is done to avoid the variation of results due to mesh-size variation in comparison to the sample size. This mesh size is applied to the area around crack tip up to 1% of the total size of the model and then increased gradually toward the outer boundary of the specimen to avoid too much computational cost. The geometry has been meshed with fully integrated 4-node quadratic elements with linear shape functions (CPE4) for the plane strain case. The mechanical boundary conditions have been applied as *K* _{I}-displacement field in a similar way as in the atomistic simulations. This is achieved by selecting the outermost nodes of the FE model and applying the boundary conditions to those nodes. Isotropic elastic constants, Young’s modulus *E* and Poisson ratio ν, were determined according to the elastic constants for the Wang potential in Table I, which yields *E* = 410 GPa and ν = 0.281.

A series of FE simulations has been performed for *K* _{I} loading using scaled parameters for the TS law according to different mesh sizes, i.e., the scaled critical stress ${\tilde{\sigma }_{\rm{c}}}$, and the displacement at failure ${\tilde{\delta }_{\rm{f}}}$ according to Eq. (20). The value for the initial slope of the TS law has been kept constant at 1.634 × 10^{5} GPa for all cases. The fracture toughness *K* _{Ic} is defined as the loading at which the first cohesive surface element breaks (δ ≥ δ_{f}).

In Fig. 8(a), it can be seen that the values of scaled critical stresses decrease over several orders of magnitude with the increasing mesh size. The values of final displacements (not shown) increases correspondingly to keep the cleavage energy (WoS) constant.

Figure 8(b) shows the obtained values of *K* _{Ic} for different mesh sizes. The value *K* _{Ic} varies between 1.35 and 1.49 MPa m^{1/2} over variation of mesh size by 5 orders of magnitude. The average of the *K* _{Ic} values over the range of mesh sizes is 1.45 MPa m^{1/2}, which is close to the atomistic value $K_{{\rm{Ic}}}^{{\rm{at}}} = 1.68$ MPa m^{1/2}. This is remarkable because there is no fitting of parameters, but rather a consistent scaling of the atomistic quantities according to the mesh size over a range of five orders of magnitude. It can be seen that there is no pronounced scatter around the average value in the sense of a pronounced mesh-size dependence. To emphasize the importance of our scaling procedure, the calculation of *K* _{Ic} with a constant TS law was also performed. In this case, the scaling is done only for the 10^{−7} m mesh and the same parameters for the TS law are used for the calculation of *K* _{Ic} with the 10^{−6} m mesh (4.59 MPa m^{1/2}) and the 10^{−4} m mesh (45.82 MPa m^{1/2}). It can be seen in Fig. 8(b) that in this case the critical stress intensity factor scales with the square root of the mesh size (or in other words, the size of the sample), in agreement with the observation of Bazant^{Reference Bažant69} of *K* _{Ic} variation for self-similar geometries but different sizes.

## IV. CONCLUSIONS

In the strong stress gradient ahead of the crack tip, the definition of mechanical stress is necessarily scale-dependent. Thus, only fundamental atomistic methods provide a rigorous way to define stress-related material properties, e.g., the fracture strength.

In this work, it is shown that the widely used approach to derive the critical bond separation stress from ab initio calculations of RBSs is well justified. To accomplish this, we have compared forces occurring in atomistic simulations of RBS with those at crack tips under *K* _{I} loading. The results demonstrate that they are indeed equivalent, which means that they are less environment-dependent than commonly expected. As a consequence, we can derive meaningful TS laws from much simpler model systems than from actual fracture simulations. This enables the use of computationally expensive, but quantitatively reliable methods such as ab initio DFT calculations.

Since the mesh sizes commonly used in continuum methods do not resolve the high atomistic stresses and the strong stress gradient at crack tips, we have introduced a physically meaningful approach to scale atomistic TS laws with the mesh size. It is based on the insight that the atomistically determined stresses ahead of the crack tip agree well with the prediction of LEFM. The approach defined in this work is based on properly averaging atomic forces over cohesive surface elements to convert them in mechanical stresses. This enables the use of fundamental atomistic material constants, such as theoretical strength, theoretical final separation, and work of separation as parameters in continuum fracture models.

It has also been demonstrated that the fracture toughness resulting from the thus-parameterized continuum model is independent of the mesh size, although numerical fluctuations occur. In this way, the atomistic data can be used directly to calculate TS parameters for continuum fracture models, such as CZ models or XFEM.

In our approach, the fracture model describes only the bond-breaking processes during cleavage; any additional (plastic) energy dissipation needs to be taken into account in the constitutive model of the material surrounding the crack. Therefore, the scaling method is ideally used in conjunction with micromechanical plasticity models or discrete dislocation dynamics simulations.

As ab initio tensile tests can be carried out using complex chemical compositions and including grain or interphase boundaries, the scheme enables the material-specific atomistically informed continuum modeling of fracture of engineering materials.

## ACKNOWLEDGMENTS

EB gratefully acknowledges the funding from European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme through the project “microKIc—Microscopic Origins of Fracture Toughness” (grant agreement No. 725483). AH and JJM gratefully acknowledge funding by the Deutsche Forschungsgemeinschaft (DFG) through Projects C3 and C4 of the collaborative research center SFB/Transregio 103 superalloy single crystals under grant number INST 213/747-2.