The atomic scale origin of crack resistance in brittle fracture

A. Mattoni and L. Colombo∗

INFM-SLACS Sardinian LAboratory for Computational materials Science
and Dipartimento di Fisica Università di Cagliari,
Cittadella Universitaria, I-09042 Monserrato (Ca), Italy

F. Cleri
ENEA, Unità Materiali e Nuove Tecnologie, and INFM
Centro Ricerche della Casaccia, CP 2400, I-00100 Roma, Italy
(Dated: May 23, 2005)
We investigate the physical meaning of the intrinsic crack resistance in the Griffith theory of brittle
fracture by means of atomic-scale simulations. By taking cubic SiC as a typical brittle material,
we show that the widely accepted identification of intrinsic crack resistance with the free surface
energy underestimates the energy release rate. Strain dependence of the Young modulus and surface
energy, as well as allowance for lattice trapping, improve the estimate of the crack resistance. In
the smallest scale limit, crack resistance can be fitted by an empirical elasto-plastic model.

PACS numbers: 62.20.Mk; 81.40.Jj; 81.40Np

A basic result of fracture mechanics for brittle ma- G = 2γs where the material parameter γs is the integral
terials is represented by the Griffith theory for crack of the stress vs. separation curve for the atomic planes
stability,1 which describes a planar crack in a homoge- undergoing separation during the fracture process. While
neous medium as a reversible thermodynamic system. γs may ideally coincide with the energy γ of the cleavage
The total energy of the system is written as the sum surface, it could include also other atomic level details.
of a mechanical contribution due to the strain energy ab- Moreover, the maximum σf of such a curve defines the
sorbed from the external work, and a crack resistance critical fracture stress, i.e., the value at which a brittle
term, originating from the material resistance to create crack starts propagating. However, the actual maximum
new free surface by breaking bonds. A crack of given stress at the crack tip could overcome σf before propa-
length is stable at the critical value of load for which the gation because of various features, including lattice trap-
total energy of the system is stationary. Once the load ex- ping, non-ideal sharpness of the (finite-size) crack tip,
ceeds such a critical value, the energy-release rate for unit anharmonic and many-body effects in the atomic-level
area of crack advance G becomes larger than the intrinsic forces. While the Griffith theory does not depend on
crack resistance and therefore the crack propagates. In σf but only on the integral γs , other theories tried to
a perfect homogeneous solid in vacuum the crack resis- also include the role of the maximum stress at the crack
tance energy per unit surface is, in fact, identified with tip.8,11,12
the (unrelaxed) free surface energy γ.2,3 The Griffith cri-
Atomistic simulations offer the opportunity to study
terion was extensively verified in glass specimens contain-
the fundamental issues underlying the Griffith theory in
ing cracks of controlled length and it is still adopted to
ideally pure, perfect single-crystal materials. In this Let-
estimate the surface energy of a brittle material.2,4
ter we present an atomic-scale investigation of brittle
At a more fundamental level it is known that, even fracture addressing the relationship among the critical
for a perfectly brittle material, some modification of the load to fracture, intrinsic crack resistance and surface
Griffith theory has to be taken into account in order to energy, in the framework of the Griffith theory. We find
describe the atomistic nature of the interactions.3 As that the identification of crack resistance with the surface
already indicated by Griffith in his original study, this energy provides only a lower limit to the energy release
turns out to be particularly relevant for microcracks of rate. By including the strain dependence of the Young
very short length. First of all, it was shown by means of modulus and surface energy, a large part of the discrep-
a simple lattice-model of brittle crack propagation in a ancy between the atomistic results and Griffith theory
two-dimensional crystal5 that the discretness of the lat- may be recovered. Moreover, we set an upper bound
tice may increase the effective crack-tip force necessary to for the role of lattice trapping in increasing the effective
break a bond across the surface, therefore causing crack crack resistance. Finally, for the shortest microcracks of
arrest for some range of loads above the theoretical Grif- length of a few lattice spacings, we give an empirical fit
fith value. Such a phenomenon was termed lattice trap- to an elasto-plastic model, asymptotically merging into
ping and was subsequently studied in a number of atomic- the Griffith theory at longer crack length. We focus our
scale simulations.6–9 However the results are sometimes work on cubic silicon carbide (β-SiC) since it is the pro-
difficult to generalize. Secondly, the effective crack resis- totype of an ideally brittle material up to extreme values
tance could include terms beyond the mere energy of the of strain, strain rate and temperature, and because of
unrelaxed cleavage surface. Formally, G is defined10 as its technological relevance as a structural and nuclear
 2

0.40

0.30

σ[111] (eV/Å−3)
0.20

0.10

Z [111]
L

0.00
0000000
1111111
1111111
0000000
0000000
1111111
0.170

γ (eV/Å−2)
2c
0.165

0.160
Y[110]
X [112] 0.155
0 0.02 0.04 0.06 0.08 0.1

L ε[111]

FIG. 1: Geometry and orientation of the simulation cell. In FIG. 2: Top panel: stress-strain curve for SiC (present Tersoff
the present simulations 22 nm < L < 88 nm. The shaded model); bottom panel: unrelaxed surface energy dependence
area represents the crack position. on strain.

material. The available data2 for β-SiC are accurate along the z direction according to a given strain value
enough (despite the microstructural heterogeneities of ex-  = [111] δzz , while keeping xx = 0 and yy = 0 (plane-
perimental samples) to suggest that its intrinsic crack re- strain condition). Periodicity is then removed along z
sistance in vacuum is, indeed, higher than the theoretical, and surface tractions are calculated, in order to preserve
ideal-crystal surface energy. the state of deformation. At this stage a microcrack of
We carried out damped-dynamics atomistic simula- given length is introduced by mathematically cutting the
tions with the aim of reproducing the macroscopic condi- interatomic bonds across a segment of a central (111)
tions of a quasi-static (or adiabatic) crack loading process plane. The interbond distance along [112̄] is c0 =2.644
at T=0 K. External loading was represented in terms of Å. The actual minimum-energy atomistic configuration
surface forces, or tractions,13 applied at the non-periodic is eventually obtained by the above specified damped-
borders of the simulation cell (see below). Atomic po- dynamics procedure. After crack opening, interatomic
sitions were relaxed according to the local values of the forces are fully restored and the microcrack reaches its
forces and by constantly damping velocities to zero, until equilibrium shape following stress relaxation.
the maximum force was less than 0.0001 eV/Å. When applying tensile loads, one can reach a strain
Atomic forces were calculated according to the Tersoff state such that the Si-C bond length at the crack tip lies
model.14 Such an empirical interatomic potential has al- in the range of values for which the cutoff function of
ready been applied15 to the study of mechanical proper- the Tersoff potential operates. Notably, the Tersoff po-
ties in β-SiC and it is able to describe the experimentally tential was not originally intended for being used under
observed brittle behavior16 of cubic β-SiC. Furthermore, such extreme deformations. To overcome artifacts due
the same force model has been applied to investigate the to the cutoff, Tang et al.15 shifted the cutoff distance to
static mechanical response of nanostructured β-SiC to larger values as a function of the applied deformation.
uniaxial tensile loading.17 At variance with that work, in the present study the sys-
tem is deformed non-homogeneously since the crack acts
The simulation cell is represented schematically in
as a stress concentrator.17 Therefore, we set a local cri-
Fig.1. The lowest unrelaxed surface energy of β-SiC is
terion to adjust the cutoff distance. In practice, atoms
that of the (111) shuffle plane,16 having the lowest den-
at the crack tip are allowed to preserve their interaction
sity of dangling bonds. As a consequence, (111)-plane
with the original neighbors of the undeformed systems,
cracks are the most likely to form in experimental con-
by shifting the cutoff distance as the crack-tip bonds are
ditions, and we therefore focused our theoretical analysis
stretched. With such a choice, atoms undergoing bond-
on such a crack arrangement. The simulation cell has the
breaking can explore the full Tersoff curve without any
x,y, and z cartesian axes parallel to the [112̄],[1̄10], and
artifact while all the bulk properties are unchanged.
[111] crystallographic directions, respectively. Therefore,
The stress-strain curve obtained with such a modifica-
the crack front lies parallel to the [112̄] direction, i.e., the
tion of the cutoff agrees with that of Tang et al.15 and is
crack arrangement is (111)[112̄].
represented in Fig.2, top panel. The slope at vanishing
We took special care in order to avoid finite-size effects.
deformation (dashed line) is related to the Young modu-
The size of the simulation cell was chosen so as to get in
lus E through the following equation:18
any case a ratio of L/c > 10 (see Fig.1) as indicated
by previous molecular dynamics studies.7 The resulting E (1 − ν)2
number of atoms ranged from 30000 up to ∼ 250000. σ[111] = [111] = E 0 f [111] (1)
1 − ν 2 1 − 2ν
The external load was applied according to the
constant traction method.13 To this end, the three- where ν is the Poisson coefficient, E 0 = E/(1 − ν 2 ). At
dimensional periodic simulation box is initially deformed deformations 0.02 < [111] < 0.12 the Poisson coefficient is
 3

0.1 erage surface traction which preserves the applied strain

0.09 or, equivalently, from the asymptotic value of the atomic-
0.08 level virial stress equation.
0.07 Consistently with the expected brittle behaviour, we
0.06
found that at loads above the critical strain the micro-
εf

0.05
crack extends in a perfectly brittle way, by preserving
atomically smooth (111) cleavage surfaces. Such a re-
0.04
sult is granted only by the above described modification
0.03
to the cutoff function: with the original Tersoff cutoff we
0.02 observed either crack deflection or incipient plasticity, de-
0 10 20 30 40 50 pending on the loading conditions. On the other hand,
2c / c0
no rehealing of the microcrack edges was ever observed in
FIG. 3: Critical strain, f , as a function of the crack size, our simulations at subcritical values of the load, i.e., the
2c, in units of the [112̄] interbond distance c0 . Symbols are microcrack does not receed back to the perfect crystal.
the data from atomistic simulations; the continuous line is Concerning the limits of the Griffith theory, we note
Griffith’s theory with constant material parameters, and the that for microcrack lengths 2c < 10c0 the critical strain
dashed line is the modified Griffith’s theory, see text. The exceeds the value 0.05, i.e., the value at which devia-
horizontal error bars in the figure are due to the lattice spacing tions from linearity start to appear in the Young modu-
orthogonal to the crack front, while the vertical error bars are lus. As a consequence, this is the minimum crack length
due to the steps chosen to vary the strain, ∆ = 0.001. for which the Griffith theory is applicable in our model
of β-SiC.
in the range 0.047 < ν < 0.07 and 1.0024 < f < 1.0057. For microcracks longer than ∼10c0 the calculated criti-
Consequently we can assume f = 1 with 0.5% tolerance. cal strain is systematically higher than the Griffith theory
In our simulations it is possible to calculate E 0 by the prediction. It is worth noting that, while atomistic sim-
numerical derivative of the stress at zero strain: E 0 = ulations on metallic systems with long-range interatomic
E/(1 − ν 2 ) = dσ[111] /d[111] . We obtain E 0 = 567 ± 1.2 potentials reported a substantial agreement with the
Griffith theory,7 Bernstein et al.6 found a similar discrep-
GPa in good agreement with previous results.19 Finally,
ancy in silicon. According to that study such a discrep-
the value of the unrelaxed (111) surface energy obtained
ancy is by definition the lattice trapping R = σfat /σfG .
with the present Tersoff model is γ=0.158 eVÅ−2 at van-
We believe, however, that the difference between atom-
ishing deformation. It is important to stress that the sur-
istic simulations and the Griffith curve deserves further
face energy to use in the Griffith theory is the unrelaxed
investigation before getting to firm conclusions. The as-
one, since the newly formed crack surfaces are infinitesi-
sumptions of the linear elastic fracture mechanics are,
mal portions of cleavage surfaces from the bulk crystal.
strictly speaking, not correct in the case of a realistic
According to Griffith fracture theory, the critical stress
(i.e., anharmonic) force model as is the present case.
σfG for a sharp planar crack of length 2c is a function of
There are, indeed, at least two possible corrections to
the Young modulus E and of the crack resistance γs :
the common interpretation of the crack resistance term
r in the Griffith theory: (i) the surface energy γ depends
G 2γs E 0
σf = . (2) on the state of strain; (ii) the stress-strain curve is not
πc strictly linear over the range of explored loads, therefore
As said above, γs is usually identified with the free surface the Young modulus E 0 is not constant.
energy γ of the cleavage surface. As a matter of fact, the non-linear dependence of the
A series of atomistic simulations was performed with surface energy γ vs. strain [111] (see Fig.2) can also be
microcracks of length 2c0 < 2c < 50c0 . Based on the computed straightforwardly. According, both the surface
Griffith formula the critical load increases with decreas- energy and Young modulus strain dependence can be in-
ing microcrack length. For the Griffith theory to be valid, troduced into the Griffith curve. The result of such a
the limits of applicability of linear elasticity must be re- modified Griffith theory is the curve reported in Fig.3 as
spected. Such a requirement implicitly defines the mini- a dashed line. The agreement between atomistic data and
mum length at which a finite-size microcrack can still be the modified Griffith theory is now much better, within
considered a “Griffith crack”. the reported error bars. This demonstrates that the com-
In Fig.3 we report the values of critical strain f as a mon identification of the crack resistance γs with the free
function of microcrack length, obtained from our atom- energy γ of the cleavage surface provides only a lower
istic simulations (symbols). The full curve represents the bound to the energy release rate G.
critical strain corresponding to the σf (c) predicted by To quantify more accurately the sources of the ob-
the Griffith theory, by using the above defined values at served discrepancy, in Fig.4 we represent our results in a
 = 0 condition of E 0 and γ for the Tersoff potential. new form. The quantity γs = (σfat )2 πc/2E 0 is the intrin-
The corresponding critical value of the atomistic stress sic crack resistance obtained from atomistic simulations.
at failure σfat is obtained either from the value of the av- In the original form of Griffith’s theory γs does not de-
 4

1.5 pendences are included and provide a γs = 1.10γ as the

1.4 asymptotic limit. The remaining discrepancy may be at-
1.3 tributed to a possible lattice trapping (proportional to
1.2
the square root of the discrepancy), the value R = 1.08
1.1
representing an upper bound for such an effect which,
1
however, includes also the statistical error bars of the
γs / γ

0.9
0.8
simulation results.
0.7 For very short microcracks, of length of a few c0 , the
0.6 critical stress σfat is so high that it becomes difficult to
0.5 discriminate between bond breaking and incipient plas-
0.4 ticity. A description of this regime can be attempted by a
0.3
0 10 20 30 40 50
fit to an empirical elasto-plastic law, such as the Dugdale-
2c / c0 Bilby-Cottrell-Swinden model (DBCS).12 In this case,
the model fracture stress can be deduced by inverting
FIG. 4: Crack resistance γs as a function of the microcrack the expression for the (unknown) crack tip displacement
length 2c. Symbols are atomistic simulation data; the con- δ, as:
tinuous line is the original Griffith theory; the long-dashed
  
line is the modified Griffith theory, with strain-dependent sur- 2σM π∆ 1
face energy and Young modulus; the continous line is the σfD = cos−1 exp − . (3)
π 4(1 − ν) c
fit to the DBCS elasto-plastic model. The horizontal short- The lumped length parameter ∆ should be equal to ∆ =
dashed line at γs /γ = 1.25 represents the asymptotic value for
µδ/σM , with µ the shear elastic modulus, and σM the
the infinite-crack intrinsic resistance estimated from atomistic
simulations. ideal cohesive strength, in the original DBCS model.
The best fit of the DBCS model to the atomistic data
is represented in Fig.4 by a continous curve, merging
pend on the crack length and is therefore a constant: an with the atomistically-corrected Griffith theory result
horizontal alignment of the data is expected if Griffith’s at longer crack lengths. It is worth noting that, with
theory prediction is correct. the fitted values of the parameters, σM = 53 GPa and
We exclude for the moment the first three simulation ∆ = 2.3c0 , we obtain an estimate of the crack tip open-
points from the discussion, since they belong to extremely ing δ ∼ 0.7c0 . This means that, for an ideally brittle
small microcracks for which the Griffith theory does not material, the extent of a “plastic” zone in the incipient
apply. In Fig.4 it can be seen that, although the long- microcrack (a “flaw”) is, indeed, vanishingly small. How-
dashed curve corresponding to the modified Griffith the- ever, this analysis also underscores the presence of com-
ory grossly agrees with the atomistic data, some sys- peting instability modes, e.g., originating from a shear re-
tematic discrepancy arises asymptotically for macroscop- sponse, for microcracks of atomic-size length under very
ically long cracks. This gives a 25% departure from the high stress.
classic Griffith theory (horizontal full line at γs = γ). This work has been funded by Italian MIUR un-
However, a substantial part of this discrepancy is due der P.R.I.N.-2002 and under F.I.S.T.-2002 PROMOMAT
to the lack of strain-dependence of the materials param- projects. Two of us (A.M., L.C.) acknowledge useful dis-
eters. In the modified Griffith theory, these strain de- cussions with A. Carpinteri and N. Pugno (Torino, Italy).

∗
Corresponding author: luciano.colombo@dsf.unica.it (2000).
1 10
A. A. Griffith, Phil. Trans. Roy. Soc. London A 221, 163 J. R. Rice, J. Mech. Phys. Solids 40, 239 (1992).
11
(1920). J. R. Rice and J.-S. Wang, Mater. Sci. Eng. A 107, 23
2
B. R. Lawn, Fracture of Brittle Solids (Cambridge Univer- (1989).
12
sity Press, Cambridge, 1975). B. A. Bilby, A. H. Cottrell, and K. H. Swinden, Proc. Roy.
3
K. B. Broberg, Cracks and Fracture (Academic Press, San Soc. London A 272, 304 (1963).
13
Diego, 1999). F. Cleri, Phys. Rev. B 65, 014107 (2002).
4 14
P. P. Gillis, J. Appl. Phys. 36, 1374 (1965). J. Tersoff, Phys. Rev. B 39, 5566 (1989).
5 15
R. Thomson, C. Hsieh, and V. Rana, J. Appl. Phys. 42, M. Tang and S. Yip, Phys. Rev. B 52, 15150 (1995).
16
3154 (1971). M. Tang and S. Yip, J. Appl. Phys. 76, 2719 (1994).
6 17
N. Bernstein and D. W. Hess, Phys. Rev. Lett. 91, 025501 A. Mattoni, L. Colombo, and F. Cleri, Phys. Rev. B 70,
(2003). 094108 (2004).
7 18
F. Cleri, S. R. Phillpot, D. Wolf, and S. Yip, J. Am. Ceram. L. D. Landau, Theory of elasticity (Butterworth-
Soc. 81, 501 (1998). Heinemann, Oxford, 1986).
8 19
W. A. Curtin, J. Mater. Res. 5, 1549 (1990). W. Li and T. Wang, Phys. Rev. B 59, 3993 (1999).
9
R. Perez and P. Gumbsch, Phys. Rev. Lett. 84, 5347