Notes On Fracture

Notes on Fracture Michael P. Marder Center for Nonlinear Dynamics and Department of Physics The University of Texas at Austin Austin Texas, 78712, USA February 7, 2003

Contents 1 Continuum Theories of Fracture: Scaling 1.1 Introduction . . . . . . . . . . . . . . 1.2 Scaling theory I: Things fall apart . . 1.3 Scaling Theory II: All together . . . . 1.4 Scaling theory III: Cracks at last . . . 1.5 The scaling of dynamic fracture . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

2 Continuum fracture mechanics 2.1 Structure of fracture mechanics. . . . . . . . . . . . . 2.2 Dissipation and the process zone . . . . . . . . . . . . 2.3 Conventional fracture modes . . . . . . . . . . . . . . 2.4 Stresses Around an Elliptical Hole . . . . . . . . . . . 2.5 Stress Intensity Factor . . . . . . . . . . . . . . . . . . 2.6 Universal singularities near the crack tip . . . . . . . . 2.7 The relation of the stress intensity factor to energy flux 2.8 Elasticity in Two Dimensions . . . . . . . . . . . . . . 2.9 Steady States in Plane Stress . . . . . . . . . . . . . . 2.10 Energy flux and limiting crack speed . . . . . . . . . . 2.11 Crack paths . . . . . . . . . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . .

. . . . . . . . . . . . . . . .

. . . . .

2 2 3 3 4 6

. . . . . . . . . . .

9 9 9 11 12 14 15 16 16 19 22 23

3 Incorporating Plasticity

25

4 Atomic Theories of Fracture 4.1 Physical argument for velocity gap . . . . . . . . . . . . . 4.2 Dynamic fracture of a lattice in anti-plane shear, Mode III . 4.3 Definition and energetics of the model . . . . . . . . . . . 4.4 Symmetries: . . . . . . . . . . . . . . . . . . . . . . . . . 4.5 Eliminating the spatial index m: . . . . . . . . . . . . . . 1

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

. . . . .

26 28 31 31 34 35

4.6 4.7 4.8 4.9 4.10 4.11 4.12

Equations solved in terms of a single mass point on crack line: Relation between and
. . . . . . . . . . . . . . . . . . . Phonon emission: . . . . . . . . . . . . . . . . . . . . . . . . Forbidden velocities . . . . . . . . . . . . . . . . . . . . . . . Nonlinear Instabilities . . . . . . . . . . . . . . . . . . . . . . The connection to the Yoffe instability . . . . . . . . . . . . . Generalizing to Mode I: . . . . . . . . . . . . . . . . . . . . .

5 Scaling ideas and molecular dynamics 5.1 Molecular dynamics . . . . . . . . . . . . . . . . 5.2 Interatomic potentials . . . . . . . . . . . . . . . 5.3 Realistic potentials for silicon? . . . . . . . . . . 5.4 Results of zero Kelvin calculations in silicon . . . 5.5 Along  110  . . . . . . . . . . . . . . . . . . . . 5.6 Along  111  : Crackons . . . . . . . . . . . . . . 5.7 Crack behavior at nonzero temperatures in silicon 5.8 Comparisons with experiment . . . . . . . . . . 5.9 Final thoughts . . . . . . . . . . . . . . . . . . .

. . . . . . . . .

. . . . . . . . .

. . . . . . . . .

. . . . . . . . .

. . . . . . . . .

. . . . . . . . .

. . . . . . . . .

. . . . . . .

. . . . . . .

. . . . . . .

. . . . . . .

. . . . . . .

. . . . . . .

. . . . . . .

35 38 40 43 43 45 46

. . . . . . . . .

. . . . . . . . .

. . . . . . . . .

. . . . . . . . .

. . . . . . . . .

. . . . . . . . .

. . . . . . . . .

48 49 49 50 53 55 55 57 58 61

1 Continuum Theories of Fracture: Scaling 1.1 Introduction The world is full of phenomena that are completely familiar but become more and more unexpected the more one thinks about them. Fracture is one of them. Who would guess that throwing a brick at a pane of glass would lead a few atomically thin regions of separation to run at the speed of sound from end to end? It is common to say at the outset of papers by physicists on fracture that little is understood. In fact, few problems have been studied more exhaustively, or solved for practical purposes more completely. Still, it is true that much remains to be solved. The problem has many features that make it particularly challenging. Cracks extract energy from the macroscopic reaches of a solid and funnel it spontaneously to regions as small as the atomic scale. Cracks can creep at rates less than an atom per second, or run at the speed of sound. They travel in solids well described by linear elasticity, but near the tip, where the action is, they are far from equilibrium. The crack connects atomic and macroscopic scales, and connot be understood without including both scales at once. The goal of these lectures is to lay out some of the elements of the physics of fracture. For the purposes of engineering, the most important concern is to prevent objects from breaking. From the perspective of physics, other issues are also important. In particular, before a physics audience, I will want to ask a set of questions that is not settled

 



How can one calculate the energy needed for a crack to propagate? When is the motion stable or unstable? How does a crack choose its path? 2

The lectures will begin by establishing the basic stability criteria that determine when fractures to begin moving. They will proceed to simple scaling theories for rapid crack motion, and then to a brief description of linear elastic fracture mechanics, an ingenious framework that allows a huge number of practical questions about fractures to be answered without needing to resolve how they behave near the tip. However, linear elastic fracture mechanics is not complete, and the lectures will continue with a discussion of current attempts to provide the missing information, including continuum theories of crack tips, atomic theories of crack tips, and techniques to consider cracks simultaneously at many scales.

1.2 Scaling theory I: Things fall apart Consider a piece of rock, of area A and height h. According to equilibrium principles the rock should not be able to sustain its own weight under the force of gravity if it becomes tall enough. I begin with a simple estimate of what the critical height should be. The gravitational potential energy of the rock is  Ah2 g  2 where  is the density. By cutting the rock into two equal blocks of height h  2 and setting them side by side, this energy can be reduced to  Ah2 g  4, for an energy gain of  Ah2 g  4. The cost of the cut is the cost of creating new rock surface. Estimate that there is a atom per square angstrom at the surface, and that the cost of cutting a bond at the surface is 1 eV. The the cost of creating new surface ˚ 2 , or G  10 J/m2 . Taking the density to be   2000kg/m3, the critical height is G  1 eV/A at which it pays to divide the rock in two is h



4G g



14cm

(1)

So every block of stone more than 20 cm tall is unstable under its own weight. A similar estimate applies to steel, concrete, or the bones in our legs. Almost every feature in the solid world around us is out of mechanical equilibrium. Mountains, buildings, and even our own bodies should fall apart in the long run under the influence of gravity. Yet they are apparently very stable, something which can only be explained if strong energy barriers keep them intact. The next task is to estimate the size of these barriers.

1.3 Scaling Theory II: All together It appears that an easy way to obtain a rough value would be by imagining what happens to the atoms of a solid as one pulls it uniformly at two ends. At first, the forces between the atoms increase, but eventually they reach a maximum value, and the solid falls into pieces, as shown in Figure 1. Interatomic forces vary greatly between different elements and compounds, but the forces typically reach their maximum value when the distance between atoms increases by around 20% of their original separation. The force needed to stretch a solid slightly is (Figure 1) F



EA L 3

(2)

F

L

L

dL



A

F

Figure 1: Mechanically stable configurations are often far from their lowest energy state. For example, a solid completely free of flaws would only pull apart when all bonds as long a plane snapped simultaneously, despite the fact that it would generally be energetically advantageous for the solid to be split. Material Young’s Modulus GPa Iron 195-205 Copper 110-130 Titanium 110 Silicon 110-160 Glass 70 Plexiglas 3.6

Theoretical Strength GPa 43-56 24-55 31 45 37 3

Practical Strength GPa 0.3 0.2 0.3 0.7 0.4 0.05

Practical/Theoretical 0.006 0.005 0.009 0.01 0.01 0.01

Table 1: The experimental strength of a number of materials in polycrystalline or amorphous form, compared to their theoretical strength, from [1] and [2] where E is the material’s Young’s modulus, so the force per area needed to reach the break ing point is around (3)  FA  E5 As shown in Table 1, this estimate is in error by orders of magnitude. The problem does not lie in the crude estimates used to obtain the forces at which bonds separate, but in the conception of the calculation. The first scaling argument greatly underestimated the practical resistance of solid bodies to separation, while this one greatly overestimates it. The only way to uncover correct orders of magnitude is to account for the actual dynamical mode by which brittle solids fail, which is by cracking.

1.4 Scaling theory III: Cracks at last Imagining that solid objects quickly reach states of lowest energy leads to the conclusion that nothing over 1 m high is stable. Imagining that solid objects fail through the simul4

taneous failure of a macroscopic plane’s worth of bonds leads to a huge overestimation of practical strength. To obtain even the correct order of magnitude in describing the resistance of solids to failure, one must consider the presence of cracks. A crack is a region of consecutive broken bonds. A crack in an otherwise perfect material creates a stress singularity at its tip. For an atomically sharp crack tip, a single crack a few microns long is sufficient to explain the large gaps between the theoretical and experimental material strengths as shown in Table 1. The theoretical strength of glass fibers, for example, can be approached by means of acid etching of the fiber surface. The etching process serves to remove any initial microscopic flaws along the glass surface. By removing its initial flaws, a 1mm2 glass fiber can lift a piano (but during the lifting process, I wouldn’t advise standing under the piano). The stress singularity that develops at the tip of a crack serves to focus the energy that is stored in the surrounding material and efficiently use it to break one bond after another. Thus, the continuous advance of the crack tip, or crack propagation, provides an efficient means to overcome the energy barrier between two equilibrium states of the system having different amounts of mechanical energy.

 L

 



dl L dl



y x

Figure 2: The upper surface of a strip of height L and thickness  is rigidly displaced upwards by distance . A crack is cut through the center of the strip, and relieves all stresses in its wake. When the crack moves distance dl from (A) to (B), the net effect is to transfer length dl of strained material into a length dl of unstrained material. Consider the crack moving in a thin strip depicted in Figure 2. Far to the right of the crack tip the material is under tension, and stores an elastic energy  right per area 



right



1 2E 2 L 



(4)

Far to the left of the crack, the material is completely relaxed, and therefore if the boundaries of the system are rigid as the crack proceeds in a steady fashion, the energy dissipated by crack motion must be 1 2E (5) G 2 L per unit area crack advance. The reason for this assertion is that when the crack moves ahead in a steady fashion the stress and strain fields around it do not change; as shown in Figure 2, the fields translate in the direction of crack motion, and the change in energy 5

comes from the transformation of stressed material ahead of the crack to unstressed material behind. This conclusion rests upon symmetry, and does not even demand that strains ahead of the crack tip to be so small that linear mechanics be applicable. Steady velocity
 and a corresponding energy flux G are achieved in the long time limit, the natural scale in experiments as well as in analytical calculations, where the crack tip reaches dynamic equilibrium with waves reflecting from top and bottom boundaries. That is, the return of acoustic waves from the system boundaries actually simplifies the task of understanding energy balance in the thin strip geometry. According to fracture mechanics, the relationship between energy flowing to a crack tip and crack velocity is, for a given lattice direction, universal. Having obtained the relation in a strip, one knows it for any of the vast range of geometries to which fracture mechanics is applicable, such as a long crack in a large plate [3].

1.5 The scaling of dynamic fracture





l

 Figure 3: As a crack of length l expands at velocity
in an infinite plate, it disturbs the surrounding medium up to a distance on the order of l. The first analysis of rapid fracture was carried out by Mott [4], and then slightly improved by Dulaney and Brace [5]. It is a dimensional analysis which clarifies the basic physical processes, despite being wrong in many details, and consists of writing down an energy balance equation for crack motion.  Consider a crack of length l  t  growing at rate
 t  in a very large plate where a stress is applied at the far boundaries of the system, as shown in Figure 3. When the crack extends, its faces separate, causing the plate to relax within a circular region centered on 6

the middle of the crack and with diameter of order l. The kinetic energy involved in moving a region of this size is M
2  2, where M is the total mass, and
is a characteristic velocity. Since the mass of material that moves is proportional to l 2 , the kinetic energy is guessed to be of the form KE  CK l 2
2 (6)



The region of material that moves is also the region from which elastic potential energy is being released as the crack opens. Therefore, the potential energy gained in releasing stress is guessed to be of the form (7) PE "! CP l 2



These guesses are correct for slowly moving cracks, but fail qualitatively as the crack velocity approaches the speed of sound, in which case both kinetic and potential energies diverge. This divergence will be demonstrated later, but for the moment, let us proceed fearlessly. The final process contributing to the energy balance equation is the creation of new crack surfaces, which takes energy # l. The fracture energy # accounts not only for the minimal energy needed to snap bonds, but also for any other dissipative processes that may be needed in order for the fracture to progress; it is often orders of magnitude greater than the thermodynamic surface energy. However, for the moment, the only important fact is that creating new surface scales as the length l of the crack. So the total energy of the system containing a crack is given by

$ with

CK l 2


2

%  qs  l 

CP l 2 %

 qs  l &"!

# l

(8) (9)

Consider first the problem of quasi-static crack propagation. If a crack moves forward only slowly, its kinetic energy will be negligible, so only the quasi-static part of the energy,  qs, will be important. For cracks where l is sufficiently small, the linear cost of fracture energy is always greater than the quadratic gain of potential energy, and in fact such cracks would heal and travel backwards if it were not for irreversible processes, such as oxidation of the crack surface, which typically prevent it from happening. That the crack grows at all is due to additional irreversible processes, sometimes chemical attack on the crack tip, sometimes vibration or other irregular mechanical stress. It should be emphasized that the system energy  increases as a result of these processes. Eventually, at length l 0 , the energy gained by relieving elastic stresses in the body exceeds the cost of creating new surface, and the crack becomes able to extend spontaneously. One sees that at l 0 , the energy functional  qs  l  has a quadratic maximum. The Griffith criterion for the onset of fracture is that fracture occurs when the potential energy released per unit crack extension equals the fracture energy, # . Thus fracture in this system will occur at the crack length l 0 where d  qs dl Using Eq. (9) I find that l0



0



  2C#   P 7

(10)

(11)

Quasistatic energy Eqs

l0 Crack length l Figure 4: The energy of a plate with a crack as a function of length. In the first part of its history, the crack grows quasi-statically, and its energy increases. At l 0 the crack begins to move rapidly, and energy is conserved. Eq. (9) now becomes

 qs  l &' qs  l0 !

CP  l ! l0 

2

(

(12)

This function is depicted in Figure 4. Much of engineering fracture mechanics boils down to calculating l 0 , given things such as external stresses, which in the present case have all been condensed into the constant CP . Dynamic fracture starts in the next instant, and because it is so rapid, the energy of the system is conserved, remaining at  qs  l0  . Using Eq. (8) and Eq. (12), with )* qs  l0  gives


+ t , CCP  1 ! ll0 &-
max  1 ! ll0   (13) K

This equation predicts that the crack will accelerate until it approaches the speed
max . The maximum speed cannot be deduced from these arguments, but Stroh [6] correctly argued that
max should be the Rayleigh wave speed, the speed at which sound travels over a free surface. A crack is a particularly severe distortion of a free surface, but assuming that it is legitimate to represent a crack in this way, the Rayleigh wave speed is the limiting speed to expect, and this result was implicit in the calculations of Yoffe [7]. In this system one needs only to know the length at which a crack begins to propagate in order to predict all the following dynamics. As we will see presently, Eq. (13) comes astonishingly close to anticipating the results of a sophisticated mathematical analysis developed over fifteen years. This result is especially puzzling since the forms (6) and (7) for the kinetic and potential energy are incorrect; these energies should actually diverge as the crack begins to approach the Rayleigh wave speed. The reason that the dimensional analysis succeeds anyway is that to find the crack velocity one takes the ratio of kinetic and potential energy; their divergences are of exactly the same form, and cancel out.

8

2 Continuum fracture mechanics In this section I will attempt to review briefly the background, basic formalism and underlying assumptions that form the body of continuum fracture mechanics. There are many textbooks on the subject from many different points of view [8–14] , and I will only cover a small subset of topics in them. I will first discuss the reduction of linear elasticity to two dimensions. Some basic concepts common to both static and dynamic cracks will then be discussed. These will include the creation of a singular stress field at the tip of a static crack together with criteria for the onset and growth of a moving crack. I then turn to a description of moving cracks. I will first describe the formalism used to describe a crack moving at a given velocity and, using this formalism, look at the dependence of the near tip stress fields as a function of the crack’s velocity. I then discuss criteria for determination of a crack’s path and conclude the section by discussing a number of points that the theory can not address.

2.1 Structure of fracture mechanics. The structure of fracture mechanics follows the basic ideas used in the scaling arguments described in the previous section. The general strategy is to solve for the displacement fields in the medium subject to both the boundary conditions and the externally applied stresses. The elastic energy transported by these fields is then matched to the amount of energy dissipated throughout the system, and an equation of motion is obtained. For a single moving crack, as in the scaling argument above, the only energy sink existing in the system is at the tip of the crack itself. Thus, an equation of motion can be obtained for a moving crack if one possesses detailed knowledge of the dissipative mechanisms in the vicinity of the tip.

2.2 Dissipation and the process zone Unfortunately, the processes that lead to dissipation in the tip vicinity are far from simple. Depending on the type of material, there is a large number of complex dissipative processes ranging from dislocation formation and emission in crystalline materials to the complex unraveling and fracture of intertangled polymer strands in amorphous polymers. At first glance, the many different (and in many cases, poorly understood) dissipative processes that are observed would appear to preclude a universal description of fracture. A way around this problem was proposed by Irwin and Orowan [15, 16]. Fracture, together with the complex dissipative processes occurring in the vicinity of the tip, occurs due to intense values of the stress field that occur as one approaches the tip. As I will show shortly, if the material surrounding the tip were to remain linearly elastic until fracture, a singularity of the stress field would result at the mathematical point associated with the crack tip. Since a real material cannot support singular stresses, in this vicinity the assumption of linearly elastic behavior must break down and material dependent, dissipative processes must come into play. Irwin and Orowan independently proposed that the medium around the crack tip be divided into three separate region as follows.

9

?@5

Figure 5: Structure of fracture mechanics. The crack tip is surrounded by a region in which the physics is unknown. Outside this process zone is a region in which elastic solutions adopt a universal form. The process zone: In the region immediately surrounding the crack tip, called the process zone (or cohesive zone), all of the nonlinear dissipative processes that ultimately allow a crack to move forward, are assumed to occur. Fracture mechanics avoids any sort of detailed description of this zone, and simply posits that it will consume some energy # per unit area of crack extension. The size of the process zone is material dependent, ranging from nanometers in glass to microns in brittle polymers, to centimeters or larger in ice near its melting point. The typical size of the process zone can be estimated by using the radius at which an assumed linear elastic stress field surrounding the crack tip would equal the yield stress of the material. The universal elastic region: Everywhere outside of the process zone the response of the material can be described by continuum linear elasticity. In the vicinity of the tip, but outside of the process zone, the stress and strain fields adopt universal singular forms which solely depend on the symmetry of the externally applied loads. In two dimensions the singular fields surrounding the process zone are entirely described by three constants, called stress intensity factors. The stress intensity factors incorporate all of the information regarding the loading of the material and are related, as we shall see, to the energy flux into the process zone. The larger the overall size of the body in which the crack lives, the larger this region becomes. In rough terms, for given values of the stress intensity factors, the size of this universal elastic region scales as \ L, where L is the macroscopic scale on which forces are applied to the body. Thus the assumptions of fracture mechanics become progressively better as samples 10

become larger and larger. Outer elastic region: Far from the crack tip, stresses and strains are described by linear elasticity. There is nothing more general to relate; details of the solution in this region depend upon the locations and strengths of the loads, and the shape of the body. In special cases, analytical solutions are available, while in general one can resort to numerical solution. At first glance, the precise linear problem that must be solved might seem inordinately complex. How can one avoid needing explicit knowledge of complicated boundary conditions on some complicated loop running outside the outer rim of the process zone? The answer is that so far as linear elasticity is concerned, viewed on macroscopic scales the process zone shrinks to a point at the crack tip, and the crack becomes a branch cut. Replacing the complicated domain in which linear elasticity actually holds with an approximate one that needs no detailed knowledge of the process zone is another approximation that becomes increasingly accurate as the dimensions of the sample, hence the size of the universal elastic region, increase. The dissipative processes within the process zone determine the fracture energy, # , defined as the amount of energy required to form a unit area of fracture surface. In the simplest case, where no dissipative processes other than the direct breaking of bonds take place, # is a constant, depending on the bond energy. In the general case # may well be a complicated function of both the crack velocity and history and differ by orders of magnitude from the surface energy defined as the amount of energy required to sever a unit area of atomic bonds. No general first principles description of the process zone exists, although numerous models have been proposed (see e.g. [13]).

2.3 Conventional fracture modes It is conventional to focus upon three symmetrical ways of loading a solid body with a crack. These are known as modes, and are illustrated in Figure 6. A generic loading situation produced by some combination of forces without any particular symmetry is referred to as mixed mode fracture. Understanding mixed-mode fracture is obviously of practical importance, but since our focus will be upon the physics of crack propagation rather than upon engineering applications, we will restrict attention to the special cases in which the loading has a high degree of symmetry. The fracture mode that I will mainly deal with here is Mode I, where the crack faces, under tension, are displaced in a direction normal to the fracture plane. In Mode II, the motion of the crack faces is that of shear along the fracture plane. Mode III fracture corresponds to an out of plane tearing motion where the direction of the stresses at the crack faces is normal to the plane of the sample. One experimental difficulty of Modes II and III is that the crack faces are not pulled away from one another. It is unavoidable that contact along the crack faces will occur. The resultant friction between the crack faces contributes to the forces acting on the crack, but is difficult to measure precisely. For these reasons, of the three fracture modes, Mode I corresponds most closely to the conditions used in most experimental and much theoretical work. In two–dimensional isotropic materials Mode II fracture cannot easily be observed, since slowly propagating 11

Mode I

Mode II

y z

x

Mode III

Figure 6: Illustration of the three conventional fracture modes, which are characterized by the symmetry of the applied forces about the crack plane cracks spontaneously orient themselves so as to make the Mode II component of the loading vanish near the crack tip [17], as I will discuss in Section 2.11. Mode II fracture is however observed in cases where material is strongly anisotropic. Both friction and earthquakes along a predefined fault are examples of Mode II fracture where the binding across the fracture interface is considerably weaker than the strength of the material that comprises the bulk material. Pure Mode III fracture, although experimentally difficult to achieve, is sometimes used as a model system for theoretical study since, in this case, the equations of elasticity simplify considerably. Analytical solutions, obtained in this mode, have provided considerable insight to the fracture process.

2.4 Stresses Around an Elliptical Hole It is possible to find the stresses surrounding an elliptical hole sitting within a plate that has been placed under remote tension, as illustrated in Figure 7. To simplify matters as much as possible, assume that the stresses applied to the plate are such as only to create a displacement u  uz in the z direction, which obeys Laplace’s equation,

]

2

u 0



(14)

The theory of complex variables can once again be brought to bear in order to find solutions. For the boundary problem at hand, conformal mapping is the appropriate technique. Because u is a solution of Laplace’s equation, it can be represented by u _^

a`4 % ^ a `b

2 12

(15)

where is analytic, and `c x % iy. ^ The asymptotic behavior of is easy to determine. Far from the hole, the stress ^ goes to a constant value dfe , where d is dimensionless. So



 ehg uy  e2 i i^j  x % iyY ! i^  x % j l m`b,nog ! id for c ` nqp  ^j yz

iy Rk

 yz

(16) (17)

How does the presence of the hole affect the stress field? Because the edges of the hole are free, all stresses normal to the edge must vanish. Let t be a variable that parameterizes the edge of the hole, so that (18)  x  t  y  t r

travels around the boundary of the hole as t moves along the real axis. Then T

ts g xt g yt u g g

and

 N

ts>! g yt g xt u It is easy to verify that Nv wxTv =0, and because v must be normal. Tv is tangent N g g

(19)

are tangent vector and normal vectors along the edge of the hole, so requiring normal stress to vanish means that



l l l l



 xz yz Yy N  0 ezs g ux g uy u yBsB! g yt g xt u  0 l g uy g xt ! g ux g yt  g g g g g g g g x y % s !{g ix^ g ix^ u g t  s g iy^ !|g iy^ u g t Insert (15) for u. g g g g g g g ^ t }g ^ t This equality holds only on the boundary. g m`b+ g m`b Because ~ is arbitrary up to a constant anyway, there is no worry about dropping a con^ ^ stant of integration.

(20) 0

(21) (22) (23) (24)

when ` lies on the boundary. Equation (24) appears innocent, but is in fact a powerful relation that can lead to rapid solution of complicated boundary-value problems. Suppose, to be definite, that the hole is elliptical and thus can be described by

`€ % € p

(25)

p is a number lying between 0 and 1.

with € lying on the unit circle,

€‚

ei ƒ

(26)

and „ real. When p  0, the boundary is circular, and when p along the real axis. Considering as a function of € , one has

^ ^ …€f, ^ …€f+‡^ †  € 1 

~Šˆ ‰Œ‹=



1, the boundary is a cut

‹ ~8‰Ž‹=

The notation means that if is expanded as a power series in , one should take the complex conjugate of all the coefficients in the expansion.

13

(27)

because € †  1 8€ on the unit circle. Equation (27) is a relation between two analytic functions that holds over the whole unit circle. The difference between €f and †  1 Š€f is an ^ ^ analytic function; the Taylor series of the difference vanishes on the unit circle, and by analytic continuation the difference vanishes everywhere, so the two functions must be equal everywhere in the complex plane. Now …€f can be determined completely by analyzing its asymptotic behavior. Outside ^ of the hole must be completely regular, except for the fact that it diverges as ! i d‘` for ^ large ` . The relation between € and ` is

€

` %“’ ` 2 !

4p

2



If the other sign of the square root were chosen here, the argument below would have to be modified in several places, but the final answer would be the same.

(28)

Therefore, when ` is large, large, €‚”` , and in accord with Eq. (17) one obtains

^ €f,n ! i •d €

for €‚n–p

(29)



Consulting Eq. (27) one concludes also that

^ †  1Š€f,n—! i•d €

for €{nqp

(30)

which means that

†^ €f>n ! € id ^ €fBn i€ d

for €n

0

for €n

0

(31) (32)



However, † …€f can have no other singularities within the unit circle, or else €f would ^ ^ have corresponding singularities outside the unit circle, which is forbidden if u is to be smooth away from the hole. Having determined all the possible singularities of within and without the unit circle, ^ it is determined up to a constant. It must be given by

^ …€fB-! i •d € % i d€ l m`b4-! i d `  1 % ’ 2 ^

(33)

Add Eqs. (32) and (29).

`  1! ’ 1!  % id 2p Use Eq. (28) and 1 ˜m‹0™‘šF‰ 1 ›•œ 1 › 4p ˜š D ˜ 2p. 1 ! 4p =`

4p =`

2

2



(34)

2

 Displacements and stresses can now be determined from Eqs. (15) and (16). Figure 7 shows a plot of yz along the line y  0 for a narrow ellipse with p  0 9. 

2.5 Stress Intensity Factor The limit p n 1 is particularly interesting. The elliptical hole closes down and becomes a thin crack reaching from x ž! 2 to x  2. The function acquires a branch cut over

^

14

e¡d

e¡d

e¡d

Radius of curvature R l

e¡d

yz

15

along y

x

10

Ÿ

5

Stress

z

 

0

y

e¡d¢

20

yz ,

e&d



0 1

3 2 Distance x

4

 in a solid are much greater than those Figure 7: The stresses at the tips of an elliptical hole applied off at infinity. The plot shows a solution of yz derived from Eq. (35) for d{ 1 and m  0 9. Note the 20-fold increase in stress near the crack tip. 



exactly the same region. The displacement u is finite everywhere, but the stress singular, diverging on the x axis as



x as x n 2 n e¡d %  x ! 2\ x 2 \ x! 2 The stress intensity factor K is a coefficient of the diverging stress, defined by  yz

-ehg uy  \ g

e¡d

r ™ x›  rlim £ 0 \ 2¤ r yz  In the case of Eq. (35), K  \ 2¤&& e d .

K

2 is the distance from the crack tip.

yz

becomes (35)

(36)

2.6 Universal singularities near the crack tip The calculations outlined in the previous section demonstrate that stresses become singular near the tip of a long thin hole. The computations were restricted to the case of Mode III, antiplane shear, for which the algebra is relatively simple. The cases of Mode I and Mode II, which are of greater practical importance, require lengthier expressions, but are not radically different. As one approaches the tip of a crack in a linearly elastic material, the stress field surrounding the tip develops a square root singularity. As first noted by Irwin [18], the stress field at a point  r „= near the crack tip, measured in polar coordinates

with the crack line corresponding to „@ 0, takes the form  1 (37) y fi¦ ¥ j §
„E i ¥ j  K¦ \ 2¤ r 15

where
is the instantaneous crack velocity, and ¨ is an index running through Modes I, II, ¦ and III. For each of these three symmetrical loading configurations, f i ¥ j a
„E in Eq. (37) is

a known universal function. The coefficient, K¦ , called the Stress Intensity Factor, contains all of the detailed information regarding sample loading and history. K¦ will, of course, be determined by the elastic fields that are set up throughout the medium, but the stress that locally drives the crack is that which is present at its tip. Thus, this single quantity will entirely determine the behavior of a crack, and much of the study of fracture comes down to either the calculation or measurement of this quantity. One of the main precepts of fracture mechanics in brittle materials is that the stress intensity factor provides a universal description of the fracture process. In other words, no matter what the history or the external conditions in a given system, if the stress intensity factor in any two systems has the same value, the crack tip that they describe will behave in the same way. The universal form of the stress intensity factor allows a complete description of the behavior of the tip of a crack where one need only carry out the analysis of a given problem within the universal elastic region (see 2.2). For arbitrary loading configurations, the stress field around the crack tip is given by three stress intensity factors, K¦ , which lead to a stress field that is a linear combination of the pure Modes:  ©3 1 y fi¦ ¥ j a
„= K¦ (38) i¥ j  ¦«ª 1 \ 2¤ r

2.7 The relation of the stress intensity factor to energy flux How are the stress intensity factors related to the flow of energy into the crack tip? Since one may view a crack as a means of dissipating built-up energy in a material, the amount of energy flowing into its tip must influence its behavior. Irwin [19] showed that the stress intensity factor is related to the energy release rate, G, which is defined as the quantity of energy flowing into the crack tip per unit fracture surface formed. The relation between the two quantities has the form © 3 1 !®­ 2 G A¦¡a
4 K¦ 2 (39) E ¦¬ª 1

where ­ is the Poisson ratio of the material and the three functions A ¦¯§
B depend only upon the crack velocity
. This relation between the stress intensity factor and the fracture energy is accurate whenever the stress field near the tip of a crack can be accurately described by Eq. (37). The near field approximation of the stress fields surrounding the crack tip embodied in Eq. (38) becomes increasingly more accurate as the dimensions of the sample increase.

2.8 Elasticity in Two Dimensions The complete formal development of fracture mechanics is very elaborate. For static cracks, there is a lovely formalism developed by Muskhelishvili [20] that allows one to use conformal mapping to obtain the stress fields around two–dimensional voids in a solid. I will present here a variant of this technique that allows one to find fields surrounding a 16

moving crack. I will not take it very far; just far enough to find the universal parts of the fields in the vicinity of the crack tip. The starting point is the equation of motion for an isotropic elastic body in the continuum limit, 2u  g t 2 °x± % e¡ ]  ] y6u  % e ] 2u (40) g originally found by Navier. u is a field describing the displacement of each mass point from its original location in an unstrained body,  is the density, and the constants e and ± are called the Lam´e constants, have dimensions of energy per volume, and are typically of order 1010 ergs/cm3 . The tensor structure of elasticity makes it particularly difficult to solve fully threedimensional problems, and it is difficult to carry out controlled three-dimensional experiments as well. Fortunately, there are cases in which the theory naturally reduces to two dimensions, where most of the analytical results have been obtained. A first case is called anti-plane shear. Imagine one is tearing a telephone book, with one hand gripping on the left, the other gripping on the right, one pushing up and the other pulling down. The only non-zero displacement is uz , and it is a function of x and y alone. The only non-vanishing stresses in this case are



xz

²e g

yz

²e‘g



and

g g

The equation of motion for uz is 1 2 uz g c2 t 2

g

where

c

uz x

(41)

uz y

(42)

 ]

e

2

uz

(43)

(44)

‘

Therefore, the vertical displacement obeys the ordinary wave equation. A second case corresponds to pulling on a thin plate, and is called plane stress. Let the z direction be the direction that goes through the plate. If the scale over which stresses are varying in x and y is large compared with the thickness of the plate, then one might expect that the displacements in the z direction will come quickly into equilibrium with the local x and y stresses. When the material is being stretched, (think of pulling on a balloon), the plate will contract in the z direction, and when it is being compressed, the plate will thicken. Therefore, one guesses that uz



z f  u x uy 

(45) 

and that ux and uy are independent of z. One can deduce the function f by noticing that

must vanish on the face of the plate. This means that

±³ g g

ux x

% g g

uy y ´

% x± % 2¡e  g 17

g

uz z



0

zz

(46)

at the surface of the plate, which implies that uz z

f  ux uy &žg

So

and one can write

where

g

g

'! ± % ± 2e ³ g uxx % g uyy ´  g g

(47)

 ± % 2e 2e ³ g uxx % g uyy ´ g g g g g  u ¦=µ ± ¶ I¦ µ g ux·· % ezs g uxµ¦ % g x¦µ u g g g ¶  ± 2% e¡2± e ±¸ ux x

% g

uy y

% g

uz z

(48)

(49)

(50)

and ¨ and ¹ now range only over x and y. Therefore, a thin plate obeys the equations of two-dimensional elasticity, with an effective constant ± ¶ , so long as uz is dependent upon ux and uy according to Eq. (47). In the following discussion, the tilde over ± ¶ will usually be dropped, with the understanding that the relation to three-dimensional materials properties is given by Eq. (49). The equation of motion is still Navier’s equation, Eq. (40), but restricted to two dimensions. A few random useful facts: materials are frequently described by the Young’s modulus E and Poisson ratio ­ . In terms of these constants,

±$  1 % ­>Eº ­1 ! 2>­ b ±$ ¶   1 E!®­ ­ 2 b e»

E 2 1 %

­B«

(51)

The following relation will be useful in discussing two-dimensional static problems. First  note that ] y u ¼x± % 2e¡ ¦=¦ (52)



Second, taking the divergence of Eq. (40), one finds that



± % 2e g



2

t

¦=¦ ] 2 



2

¦E¦

g ] Therfore, yºu obeys the wave equation, with the longitudinal wave speed

% ± 2e cl    ]ž½ u also obeys the wave equation, but with the shear wave speed Similarly,

e ct   

(53)

(54)

(55)

The fact that the local thickness of the plate is tied to stresses in the x and y directions leads to two optical methods to determine stress fields experimentally. The first method relies on the fact that when light reflects off a curved surface, the reflected intensity becomes 18

Figure 8: Interference fringes captured during rapid motion of crack though brittle plastic. K Ravi–Chandar. singular at certain points that depend on the details of the geometry. In practice, when light is shined on a crack tip, a sharp dark spot surrounds the tip, and its shape and size can be used to deduce the stresses. This technique is known as the method of caustics. A second method, older and more reliable, relies upon the fact that materials under stress will typically rotate the plane of polarization of transmitted light. It is possible to determine the basic structure of this rotation without detailed calculation. The rotation must depend upon features of the stress tensor which are rotationally invariant, and therefore can depend only upon the two principal stresses, which are the the diagonal elements in a reference frame where the stress tensor is diagonal. In addition, there should be no rotation of polarization when the material is stretched uniformly in all directions, in which case the two principal stresses are equal. So the angular rotation of the plane of polarization must be of the form







K



1

!



2



(56)

where 1 and 2 are the principal stresses (eigenvalues of the stress tensor), and K is a constant that will have to be determined experimentally. Whenever stresses of a twodimensional problem are calculated analytically, the results can be placed into Eq. (56), and compared with experimental fringe patterns. Fast optical systems have been developed to carry out this procedure for rapidly moving cracks. Although one might wonder about the extent to which Eq. (56) is obeyed when cracks move at speeds on the order of the speed of sound, experimental results match rather well with predictions for crack tip fields.

2.9 Steady States in Plane Stress We now proceed to find the form of the singularity around the tip of a Mode I crack. Begin with the dynamical equation for the strain field u of a steady state in the moving frame, 2 x± % e¡ ]  ] yºu  % e ] 2u -¬
2 g xu2  (57) g Divide u into transverse and longitudinal parts so that u "u t % u l 19

(58)

with

 ]


ts g
yt ! g
xt u  (59) g g It follows immediately that ¾ ¾ 2 2 ] % 2 2 m± 2e¡ !¿¬
g x2 À u l  f "! e ] 2 !¿¬
2 g x2 À u t (60) g g for some function f which must be harmonic ( f x ! i fy is a function of x % iy). We have then that ¾ 2 2 (61) ¨ 2 g x2 % g y2 À ] 2
l  0 ¾ g 2 g 2 (62) ¹ 2 g x2 % g y2 À ] 2
t  0 g g u l

where

l

and u t

¬
2 ± % 2e 2 ¹ 2  1 ! ¬e


¨ 2

1!

(63a) (63b)

Therefore the general form of the potentials is


l² 
l0  z  %
l0  z %
l1  x % i¨ y  %
l1  x % i¨ y (64)
t ²
t0  z  %
t0  z  %
t1  x % i¹ y  %
t1  x % i¹ y (65) subject to the constraint of Eq. (60), which gives a relation between
l0 and
t0 . In fact, the purely harmonic pieces
l0 and
t0 disappear entirely from the expressions for u They result from the freedom one has to add a harmonic function to
l and
t simul taneously, and can be neglected. Defining  z ¡
l1  z2 g z and Ác z¯ g
t1  z2 g z I have g ^ for u ux  z¦ % z¦Y % i¹{ÂÃÁÄ zµBÅ ! Á‚Ä zµ4ÅIÆ (66a)   ^ ^  uy  i ¨ÈÇ  z¦Y! z¦É&!|ÂPÁ‚Ä zµ4Å % ÁÄ zµBÅÊÆ (66b)

^ ^ where

z¦U x % i ¨ y

zµ



x % i¹ y



(67)

Equation Eq. (66) gives a general form for elastic problems which are in steady state moving at velocity
. Define also ˓  z2 g z and ÌÍ g ÁÎ z r g z  Then the stresses are g ^ given by  



% ¡e ¡ÇPË¿ z¦ % Ë¿ z¦YÏÉÐÄ 1 !®¨ 2 Å e Ñ Ä 1 % ¨ 2 Å Ç Ë¿ z¦Ó % Ëz z¦Ó É % 2i¹ Â Ì Ä zµ4Å ! Ì Ä zµ4Å=Æ>Ô xx !  yy  2Ò 2 xy  2e Ñ 2i ¨ÈÇÕËz z¦Y! Ëz z¦É&!”Ä2¹ 2 % 1Å KÂ Ì”Ä zµ4Å % ” Ì Ä zµ4ÅIÆ>Ô xx

%



yy



2 x±

20

(68a) (68b) (68c)

we will also need the rotation, which means

žg g

uy x

!Ög g

ux y

"!×Ä 1 !¿¹ 2 Å ÂKÌ²Ä zµ4Å % Ì²Ä zµ4Å=Æ 

(69)

It is worth writing down the stresses directly as well:



yy



'!e Ä 1 % ¹ 2 Å Ç Ë¿ z¦Y % Ëz z¦ É !

2i¹¡e

Â Ì Ä zµ4Å ! Ì Ä zµ4ÅIÆ

(70)

(71) ²eØÄ 1 % 2¨ 2 ¿ ! ¹ 2 Å ÇÙËz z¦Ó % Ë¿ z¦ÏÉ % 2i¹¡eÈÂPÌ²Ä zµ4Å ! Ì²Ä zµ4Å=Æ  The definitions of ¨ and ¹ in Eq. (63) have been used to simplify the expressions. To solve a general problem, one has to find the functions and Á which match boundary ^ conditions. It is interesting to notice that when
¸n 0, the right hand side of Eq. (68a) xx

goes to zero as well. Since one will be finding the potentials from given stresses at the boundaries, Ë must diverge as 1 I
, and the right hand side of Eq  68  will turn into a  derivative of Ë with respect to ¨ . That is why the static theory has a different structure than the dynamic theory. In fact, the dynamic theory is more straightforward. As a first application, I will show that a moving crack under symmetric loading becomes unstable at a certain speed. Assume the crack to lie along the negative x axis, terminating at x  0 and moving forward. The problem is assumed symmetric under reflection about

the x axis, but no other assumption is needed. This instability was first found in a particular case by Yoffe (1951). We know that in the static case, the stress fields have a square root singularity at the crack tip. I will assume the same to be true in this case (the assumption is verified in all cases that can be worked explicitly.) Near the crack tip, I assume that

% ^  z    Br Ác z   Dr %

iBi  z1 Ú

iDi  z1 Ú

2 2

We first appeal to symmetry. Observe that ux J! y ¡ ux  y 

(72) (73)



uy J! y ¡"! u  y 

(74)



Placing Eq. (73) into Eq. (66) and using Eq. (74), we find immediately that B i  Dr  0. Thus i Ë z  B1Ú r2 ÌÛ z   iD (75) 2 1 Ú z z  We also observe that the square roots in Eq. (73) must be interpreted as having their cuts  the crack.  On the crack surface, we have two along the negative x axis, corresponding to  yy vanish. Upon substituting Eq. (75) boundary conditions, which require that xy and into Eq. (71) we find that the condition upon yy is satisfied identically for x Ü 0 y  0.

However, substituting into Eq. (69) with y  0 we find that



xy

²e i Ý 2¨

Br !

Ä ¹ 2 % 1Å 21

Di Þ

³ \

1 x

! \

1 x´



(76)

Thus

2¨ Di  (77) 2 Br ¹ % 1  This relation is enough to find the maximum velocity at which a crack can proceed stably along the x axis. Using Eq. (77) we find that



  with

xx



yy



xy



\ \ \

K Ä ¹ 2 % 1Å Ä 1 % 2¨ 2 !¿¹ 2¤ D i 2 K !i Ä 1 % ¹ 2 Å ³ 1z¦ % 2¤ D \ K 1 ! 2i ¨ Ä ¹ 2 % 1Å ³ 2¤ D \ z¦

D  4¨ß¹$!à 1 %

1 % 1 4ß ! ¨ ¹ k ³ \ ´ \ \ zµ \ z† µ ´ 1 % 4ߨ ¹ ³ 1 % 1 k \ z† ¦Î´ \ zµ \ z† µ ´ 1 1 % 1 ! \ z† ¦ \ zµ \ z† µ ´ 2

1 z¦

Å ³

%

1 z† ¦

¹ 2 2 

(78a) (78b) (78c) (78d)

The constant K is again called the stress intensity factor, and is given by K



lim x n

0á

\ 2¤



x

yy

(79)



In order to find the direction of maximum stress, one must approach the tip of the crack along a line at angle „ to the x axis, and compute the stress perpendicular to that line. So one wants to choose

%

zµ  r cos „ ri¹ sin „ and to evaluate the stress     ƒ2ƒ  cos2 „ yy % sin2 „ xx ! sin  2„E xy  z¦U r cos „

%

ri ¨ sin „

(80)

(81)

When this is carried out, one finds that above a certain velocity (for Poisson ratio ­â 1  3, the critical velocity is about .61) the direction of maximum tearing stress points away from the axis, and the crack would presumably become unstable. The point at which this happens is referred to as the Yoffe instability.

2.10 Energy flux and limiting crack speed Energy flux may be found from the time derivative of the total energy. We have d K % Pkb dt i

¾

d dt

ã

dxdy

 u ¦ u¦ % 2ä ä

1 u¦ g 2 xµ

g



¦Iµ À 

(82)

The spatial integral is taken over an region which is static in the laboratory frame. So d K % Pkb dt i

ã

¾

dxdy

22

æuå ¦ uä ¦ % g g

u¦ ä xµ



¦=µ À

(83)

where the symmetry of the stress tensor under interchange of indices is used for the last term. Using the equation of motion



0uå ¦¡ x&

(84) g xµ ¦=µ  x g  we have ¾  ¦ u % (85) ã dxdy g g xµ ¦=µ uä ¦  gg xä µ ¦Iµ À (86)  ã dxdy gxµ Ç ¦Iµ uä ¦ É g  (87)  ã;ç uä ¦ ¦=µ nµ S  system, and nè is an outward unit normal. where the integral is now over the boundary of the By using the asymptotic forms Eq. (78) for yy and the corresponding expressions for uy from Eq. (66a), one finds that the total energy flowing into the crack tip per unit time is J tot

²
 1 !¿¹ 2  2¨ e 4ߨ ¹Ò!“ 11 % ¹ 2  2 K 2

(88)

The factor 4 ¨ß¹$!à 1 % ¹ 2  2 is a function of crack velocity
whose roots happen to define the Rayleigh wave speed, which is the speed at which waves travel over a free surface. Because this factor vanishes at the Rayleigh wave speed and becomes negative for larger values of crack velocity, the Rayleigh wave speed is conventionally thought to provide an upper limiting speed for cracks. The argument seems fairly secure for this loading geometry. However, cracks loaded in Mode II, or moving along interfaces can go faster [21, 22]. I think Petersan, Deegan, and I have some experimental evidence for cracks traveling faster than the shear wave speed under Mode I loading in rubber, but we are not yet 100% certain.

2.11 Crack paths I now briefly discuss the path chosen by a moving crack. Energy balance provides an equation of motion for the tip of a crack only when the crack path or propagation direction is assumed. Although criteria for a crack’s path have been established for a slowly moving cracks, no such criterion has been proven to exist for a crack moving at high speeds. A slow crack is one whose velocity
is much less than the Rayleigh wave speed c R . The path followed by such cracks obeys the “ the principle of local symmetry,” first proposed by Goldstein and Salganik [23]. This criterion states that a crack extends so as to set the component of Mode II loading to zero. One consequence is that if a stationary crack is loaded in such a way as to experience Mode II loading, upon extension it forms a sharp kink and moves at a new angle. A simple explanation for this rule is that it means the crack is moving perpendicular to the direction in which tensile stresses are the greatest. Cotterell and Rice [17] have shown that a crack obeying this principle of local symmetry is also choosing a direction so as to maximize the energy release rate. The distance over which a crack needs to move so as to set KII to zero is on the order of the size of the process 23

zone; Hodgdon and Sethna [24] show how to arrive at this conclusion using little more than symmetry principles. Cotterell and Rice also demonstrated that the condition KII  0 has the following con sequences for crack motion. Consider an  initially straight crack, propagating along the x axis. The stress field components xx and yy have the following form





xx



yy



KI  2¤ r  1 Ú KI  2¤ r  1 Ú

2 2

%

T

%àé  r1Ú 2 

%àé  r1Ú 2 

(89a) (89b)

The constant stress T is parallel to the crack at its tip. Cotterell and Rice showed that if T ê 0, any small fluctuations from straightness cause the the crack to diverge from the x direction, while if T Ü 0 the crack is stable and continues to propagate along the x axis. They also discuss experimental verification of this prediction. Yuse and Sano [25] and Ronsin, Heslot, and Perrin [26] have conducted experiments by slowly pulling a glass plate from a hot region to a cold one across a constant thermal gradient. The velocity of the crack, driven by the stresses induced by the nonuniform thermal expansion of the material, follows that of the glass plate. At a critical pulling velocity, the crack’s path deviates from straight-line propagation and transverse oscillations develop. This instability is completely consistent with the principle of local symmetry [27,28]; the crack deviates from straightness when T in Eq. (89) rises above 0. Adda-Bedia and Pomeau and Ben-Amar [29] have extended the analysis to calculate the wavelength of the ensuing oscillations. Hodgdon and Sethna [24] have generalized the principle of local symmetry to three dimensions. They show that an equation of motion for a crack line involves in principle nine different constants. It would be interesting for an experimental study to follow upon their work and try to measure the many constants they have described, but we are not aware of such experiments. Larralde and Ball [30, 31] analyzed what such equations would imply for a corrugated crack, and found that the corrugations should decay exponentially. They also performed some simple experiments and verified the predictions. Thus the principle of local symmetry is consistent with all experimental tests that have been performed so far on slowly moving cracks. Nevertheless, it does not rest upon a particularly solid foundation. There is no basic principle from which it follows that a crack must extend perpendicular to the maximum tensile stress, or that it must maximize energy release. The logic of the principle of local symmetry says that bonds under the greatest tension break first, and therefore cracks loaded in Mode I move straight ahead, at least until the velocity identified by Yoffe when a crack is predicted to spontaneously break the symmetry inherent in straight-line propagation. This logic is called into question by a very simple  [32]. calculation, first described by Rice Let us look at the ratio xx  yy right on the crack line. From Eq. (78) it is





% 2% 2 2  a¹ 14º  ¯¨ 1¹$!à2 m1¨ % !¿¹ ¹2  2! 4¯¨ ¹ yy xx

24

(90)



2 §¹ 2 % 1 ;m¨ 4¨ß¹$!à 1 %

2

!¿¹ 2  ! ¹ 2 2

1

(91)

The Taylor expansion of Eq. (91) for low velocities
is


2  ct4 % c4l  % (92) r22 2c2l ct2  cl ! ct º cl % ct  and in fact Eq. (91) greater than unity for all
. This result is surprising because it states 1%

that, in fact, the greatest tensile forces are perpendicular to a crack tip and not parallel to it, as soon as the crack begins to move. Therefore it is hard to understand why a crack is ever supposed to move in a straight line. These remarks do not do justice to the calculations performed with the cohesive zone models. In their most elaborate versions, the crack is allowed to pursue an oscillating path, and the cohesive zone contains both tensile and shear components. In most, although not all of these models, crack propagation is violently unstable to very short-length oscillations of the tip. A summary of this work has been provided by Langer and Lobkovsky [33]. They consider a large class of models, and correct subtle errors in previous analyses. They do find some models in which a crack undergoes a Hopf bifurcation to an oscillation at a critical velocity. However, their “general conclusion is that these cohesive-zone models are inherently unsatisfactory for use in dynamical studies. They are extremely difficult mathematically and they seem to be highly sensitive to details that ought to be physically unimportant.” Attempts to find models in this general class continue, with recent attempts by Karma, Kessler and Levine [34] and Aranson, Kalatsky Vinokur [35]. Both of these papers model Mode III cracks, and it is not clear how to generalize them to Mode I. The more recent paper cures some problematic features of the first one; I am not convinced that either provides a really useful starting point to address outstanding problems in fracture. One possibility is that cohesive zone models must be replaced by models in which plastic yielding is distributed across an area, and not restricted to a line. Dynamic elasticplastic fracture has not, to our knowledge, considered cracks moving away from straight paths. Another possibility is that these calculations signal a failure of continuum theory, and that the resolution must be sought at the atomic or molecular scale. It is not possible right now to decide conclusively between these two possibilities. However, the experimental observation that the dynamic instabilities consist of repeated frustrated branching seems difficult to capture in a continuum description of the process zone.

3 Incorporating Plasticity The lectures included an improvised segment on elastic–plastic fracture mechanics. There was no written material prepared, but here is a brief outline: Elastic–plastic fracture mechanics is a large and well–developed subject within engineering, resting on a large scientific base that considers interactions of cracks and dislocations. Thomson Physics of Fracture (1986) ; Weertman Dislocation–Based Fracture Mechanics (1996) ; Janssen, Zuidema, and Wanhill Fracture Mechanics (2002) 25

Basic results:

 

When plastic deformation at tip is significant, cracks move slowly, and bond–breaking is irrelevant. Total plastic dissipation can be estimated from irreversible separation of crack faces (Crack Tip Opening Displacement):

#) 



where

c



c

(93)

is the shear stress at which plastic flow begins.

Distinction between ductile and brittle materials captured by resistance to shear. Roughly, brittle materials are those where

#Èë



10a

c

(94)

(Rice and Thomson, 1974; Armstrong; Orowan)

4 Atomic Theories of Fracture In contrast to cohesive zone and other continuum models, where the correct starting equations still are not yet known with certainty, and instabilities in qualitative accord with experiment are difficult to find, calculations in crystals provide a context where the starting point is unambiguous, and instabilities resembling those seen in experiment arise quite naturally. In this section I will record some theoretical results relating to the instability. I will first focus on a description of brittle fracture introduced by Slepyan [9,36], and Marder and Liu [37, 38]. The aim of this approach to fracture is to find a case where it is possible to study the motion of a crack in a macroscopic sample, but describing the motion of every atom in detail. In this way, questions about the behavior of the process zone and the precise nature of crack motion can be resolved without any additional assumptions. This task can be accomplished by arranging atoms in a crystal, and adopting a simple force law between them, one in which forces rise linearly up to a critical separation, and then abruptly drop to zero. I propose to call a solid built of atoms of this type an ideal brittle crystal. A force law of this type is not, of course, completely realistic, but has long been thought a sensible approximation in brittle ceramics [13]. It is more realistic in brittle materials than, for example, Lennard-Jones or Morse potentials. A surprising fact is that this force law makes it possible to obtain a large variety of analytical results for fracture in arbitrarily large systems. Furthermore, the qualitative lessons following from these calculations seem also to be quite general. A summary of results from the ideal brittle crystal is Lattice trapping: For a range of loads above the Griffith point, a crack can be trapped by the crystal, and does not move, although it is energetically possible [39, 40]. Steady states: Steady state crack motion exists, and is a stable attractor for a range of energy flux. 26

Phonons: Steadily moving cracks emit phonons whose frequencies can be computed from a simple conservation law. Fracture energy The relation between the fracture energy and velocity can be computed. Velocity gap: The slowest steady state runs at around 20% of the Rayleigh wave speed; no slower-moving steady state crack exists. Instability: At an upper critical energy flux, steady state cracks become unstable, and generate frustrated branching events in a fashion reminiscent of experiments in amorphous materials. One may wonder about the motivation for formulating a theory for the dynamic behavior of a crack in a crystal. Is the lattice essential or can one make the underlying lattice go away by taking a continuum limit? So far as I know, all attempts to describe the process zone of brittle materials in a continuum framework have run into severe difficulties [33]. These problems do not arise in an atomic-scale description. The simplicity of the ideal brittle crystal is somewhat misleading in a number of respects. Therefore, before introducing the mathematics, I will comment on a number of natural questions regarding the generality of the predictions it makes.





Does the simple force law employed between atoms neglect essential aspects of the dynamics? I will demonstrate that the same qualitative results observed in the ideal brittle crystal occur in extensive molecular dynamic simulations using realistic potentials.





The calculations are in a strip. Is this geometry too restrictive? To this question the general formalism of fracture mechanics provides an answer. Fracture mechanics tells us that as long as the conditions of small-scale yielding are satisfied, the behavior of a crack is entirely governed by the structure of the stress fields in the near vicinity of the crack tip. These fields are solely determined by the flux of energy to the crack tip. A given energy flux can be provided by an infinite number of different loading configurations, but the resultant dynamics of the crack will be the same. As a result, the specific geometry used to load the system is irrelevant and no generality is lost by the use of a specific loading configuration. This fact is borne out by the experimental work described in the previous section. Have lattice trapping and the velocity gap ever been seen experimentally? Maybe. Robert Deegan and I have obtained some results, but the data are not yet conclusive. Molecular dynamics simulations indicate that lattice trapping disappears at room temperature. New experiments are needed to obtain a detailed description of dynamic fracture in crystals at low temperatures. Are results in a crystal relevant for amorphous materials? This is still an open question. However, the results of the lattice calculations seem to be remarkably robust. Adding quenched noise to the crystal has little qualitative effect. Effects of topological disorder are not known. However, a certain amount is known about the effects 27

of increasing temperature. When the temperature of a brittle crystal increases above zero, the velocity gap closes [41], and its behavior is reminiscent of that seen in experiments performed on amorphous materials. Kulakhmetova, Saraikin and Slepyan [9, 42] first showed that it is possible to find exact analytical expressions for Mode I cracks moving in a square lattice. They found exact relationships between the energy flux to a crack tip and crack tip velocity. They also observed that phonons must be emitted by moving cracks, and calculated their frequencies and amplitudes. Later calculations extended these results to other crystal geometries, allowed for a general Poisson ratio, showed that there is a minimum allowed crack velocity, found when steady crack motion is linearly stable, calculated the point at which steady motion becomes unstable to a branching instability, and estimated the spacing between branches [38].

Figure 9: A sketch showing steady state motion of crack moving in an ideal brittle crystalline strip loaded in Mode I. These calculations are extremely elaborate. The analytical expressions are too lengthy to fit easily on printed pages, and most of the results were obtained with the aid of symbolic algebra. For this reason, I will first give a qualitative argument for the most surprising of the dynamical results, the velocity gap. Then, after summarizing the results of the technical calculations, I will proceed to describe in detail how the calculations work in the case of a simple one–dimensional model where the algebra is much less demanding but most of the ideas are the same.

4.1 Physical argument for velocity gap The velocity gap is a natural consequence of rapidly snapping bonds. It emerges both from the simple model depicted in Figure 12, and also in more elaborate models, whose solution requires greater effort, that describe cracks moving through fully two-dimensional 28

crystals. In the analytically solvable models, the bond breaking is instantaneous. In any more realistic situation there will instead be a brief period in which the force between atoms falls to zero, but all that is required for the velocity gap is that this time be short compared to the vibrational period of an atom. Dynamic fracture is a cascade of bonds breaking, one giving way after another like a toppling line of dominos. Examining Figure 10, watch what happens as a crack moves ahead. In frame 2, the bond between two atoms has just broken. There is no guarantee that the next bond to the right will break. The crack could fall into a static lattice trapped state. The best chance to avoid this fate is for the atom marked in green to deliver enough of a blow to its right-hand neighbor that the bond on that neighbor is broken in turn. This process had better take place within the first half of the first vibrational period of the green atom. For the longer it vibrates, the more energy will be distributed to its neighbors in all directions in the form of traveling waves, and the smaller become the chances that there will be enough concentrated energy available to snap the next bond down the line. To make this idea slightly more quantitative, let ì b be the effective spring constant acting on an atom after a bond connected to it has snapped. Its oscillatory period is

í 

2¤



M

(95)

ì b

where M is its mass. Assuming that the next bond, at distance a, breaks during some fraction ¨‚Ü 1 of half this oscillatory period, the speed
of the crack will be



a


 ¨¡¤

ì

b

M



(96)

By this same rough logic, the speed of sound is given by solving the wave equation M uå ²ì a

2

g g

2u

x2

(97)

where now ì is an effective spring constant appropriate for atoms surrounded by unbroken bonds. From Eq. (97) follows a wave speed c of c a



ì

M

(98)

One should expect c to be larger than the crack speed
. The main reason is that when one of the springs connected to an atom is snapped, the effective spring constant describing its vibrations should decrease, so ìîê-ì b . There is also an extra factor of ¤ in Eq. (96). The crack speed
estimated here is a lower bound because there is no way that pulling more gently on remote system boundaries can reduce the vibrational period of two atoms once the bond connecting them breaks. Pulling harder on external boundaries can however increase the speed of a crack because it can supply enough energy that the atoms are already moving quickly when the bond between them snaps. These arguments presume that the only source of energy available for snapping bonds is contained in strained material ahead of the tip. This assumption is never completely 29

1

5

2

6 Snaps

Snaps

3

7

4

8

Figure 10: Sequence of eight snapshots from motion of a crack tip. Bonds break in frames two and six. The green atom snaps backwards each time a bond connected to it breaks, providing a crucial portion of the energy for the next bond to snap. This process must occur within the first half of the first vibrational period of the green atom, or it is unlikely ever to occur.

30

correct. At any nonzero temperature, thermal fluctuations can bring extra energy in. For this reason, cracks strained above c always move ahead at any nonzero temperature, but the rate is very small when the thermal fluctuations are rare. Similarly, chemical agents in the environment can help catalyze bond breaking, and crack speeds are sometimes controlled by the rate at which they arrive. Dynamic fracture refers to the motion of cracks when these sorts of effects can be neglected. dè 1 ð0ñ dè 5 ð

dè 2 ð

dè 1 ï dè 2 ï dè 6 ð dè 5 ï dè èd3 ï dè 4 ï 6 ï dè 4 ð dè 3 ð

u 1

u 5

u 3

u 0

u 2

u 6

u 4

Figure 11: The geometry of the lattice used for fracture calculations.

4.2 Dynamic fracture of a lattice in anti-plane shear, Mode III 4.3 Definition and energetics of the model I now work in detail through some of the analytical results available for a crack moving through an ideal brittle crystal loaded in Mode III. The main results are the relation between loading and crack velocity, the prediction of a velocity gap, and the calculation of the point at which steady crack motion becomes unstable. Linear stability of the steady states will not be considered here; the result obtained in [38] is that the steady states, when they exist, are always linearly stable. The techniques used to solve problems of this type were first found by Slepyan [9, 36, 42]. There are differences between details of his solution and mine because Eq. (99) describes motion in a strip rather than an infinite plate, and in a triangular rather than a square lattice. The strip is preferable to the infinite plate when it comes time to compare with numerical simulations, while reducing to the simpler infinite plate results in certain natural limits. The velocity gap and nonlinear instabilities were first found in [37, 38]. Slepyan’s plots of crack speed versus driving contain ranges of solutions that turn out to be inconsistent. Consider a crack moving in a strip composed of 2  N % 1  rows of mass points, shown in Figure 13. All of the bonds between lattice points are brittle-elastic, behaving as perfect 31

ú

1

÷ø ù ú

1 N

õö

xm

÷¬ø ù

1 N

1

ú

1 N

xm ü

õö

÷¬ø ù

1

xm û

òoó

1

òoó

ý

þ

ÿ

2N

ú

1

÷ø ù

ú

1

÷Eø ù

1 N

1 N

ú

÷ø ù

òoó

1

1 N

1

1

xm ü

õö

xm û

1




õö

xm

òôó

1





1

1

Figure 12: In this simplest of all solvable atomic crack models, two lines of atoms are connected by weak springs to a ceiling and a floor, and to horizontal and vertical neighbors by strong springs. The vertical strong springs break if they exceed a certain extension. 32

linear springs until the instant they snap at a separation of 2u f , from which point on they exert no force. The location of each mass point is described by a single spatial coordinate u  m n  , which can be interpreted as the height of mass point  m n  into or out of the page.

The force between adjacent masses is determined by the difference in height between them. The index m takes integer values, while n takes values of the form 1  2, 3  2, N % 1  2. 2r The model is described by the equation

Figure 13: Dynamic fracture of a triangular crystal in anti-plane shear. The crystal has 2  N % 1  rows of atoms, with N  4, and 2N % 1 rows of bonds. Heights of spheres indicate displacements u  m n  of mass points out of the page once displacements are imposed at the

boundaries. The top line of spheres is displaced out of the page by amount u  m N % 1  2  UN  \ 2N % 1, and the bottom line into the page by amount u  m ! N ! 1 2 ¯¼! UN  ! \ 2N % 1. Lines connecting mass points indicate whether the displacement between them has exceeded the critical value of 2u f ; see Eq. (99b)).

uå  m n Y! bu %

with

ä

1 2

© nearest neighbors m

P¥  n

f u  m n ! u  m n Rk

i

j j



(99a)

 

f  u & u „> 2u f ! u  (99b) describing ideal brittle springs, „ the step function, and b the coefficient of a small dissipative term. There is no difficulty involved in choosing alternative forms of the dissipation, if desired. The boundary condition which drives the motion of the crack is u  m N % 1  2k> UN (99c)

Øi



It is important to find the value of UN for which there is just enough energy stored per length to the right of the crack to snap the pair of bonds connected to each lattice site on 33

the crack line. For m



0 one has u  m n &

nUN N % 1 2

(100)

and the height difference between mass points with adjacent values of n is Uright

UN



N % 1 2

(101)

Therefore, the energy stored per unit length in the 2N % 1 rows of bonds is

 

1½ 2 Upper Bonds  Sitek 2 i 1 1 UN 2 2  2N % 1   % 2 2 N 1 2 2Q0  UN  2

½ i

Rowsk

½

½

i

2 Spring Constantk Uright

(102) (103) (104)

with Q0

  2N1% 1 

(105)

The energy required to snap two bonds each time the crack advances by a unit length is



1 2 1 2

½

½ i

2 Bonds  Sitek

2

½

1 2

½

½  2u  2  f

i

Spring Constantk 2u2f

½

i Separation at fracturek

2

(106) (107)



Therefore, equating Eqs. (104) and (107) the proper dimensionless measure of the external driving is (108)  UNu\ Q0 f a quantity which reaches 1 as soon as there is enough energy to the right of the crack to snap the bonds along the crack line, and which is linearly related to the displacements imposed at the edges of the strip.

4.4 Symmetries: Assume that a crack moves in steady state, so that one by one, the bonds connecting u  m 1  2  with u  m % 1 ! 1  2  or u  m ! 1  2  break. They break because the distance

between these points exceeds the limit set in Eq. (99b) and as a consequence of the driving force described by Eq. (99c). Assuming that the times at which bonds break are known, the original nonlinear problem is immediately transformed into a linear problem. However, one has to come back at the end of the calculation to verify that 1. Bonds break at the time they are supposed to. Imposing this condition determines a relation between crack velocity
and loading . 34

2. No bonds break when they are not supposed to. Imposing this condition leads to the velocity gap on the low-velocity end, and to crack tip instabilities above a critical energy flux. Steady states in a crystal are more complicated than steady states in a continuum. In a continuum, steady state acts as u  x !z
t  . The closest one can come in a triangular crystal is by having the symmetry u  m n t 4 u  m % 1 n t % 1 I
4



and also

(109a)



u  m n t 4²! u  m



which implies in particular that u  m 1  2 t b²! u  m

We have defined

gn





 0

! n t ! i 1 2 ! gn k…Ê
B

! 1 2 t ! 1 2
4 

1 mod  n ! 1  2 2 

(109b) (109c)

if n  1  2 5  2

r2 if n  3  2 7  2

r2 in general

(110)

4.5 Eliminating the spatial index m: Assuming that a crack is in steady state, we can therefore eliminate the variable m entirely from the equation of motion, by defining un  t , u  0 n t 

(111)



and write the equations of motion in steady state as uå n  t B

1 2

%  %%

if n ê 1  2, and

%

uå 1 Ú

2

 t B

1 2

un 1  t !à gn á á un  t % 1 I
4 un 1  t !à gn



1

! 12 Ê
B ! 12 Ê
B

!

% %

6un  t 

%

%

u n 1  t ! gn á á un  t ! 1 I
4 u n 1  t ! gn



 t ! 1 Ê
B % % %  t 1I
4! 4u  t   t ! 1 Ê
B %  21Ú 2  t ! u1Ú 2  t R1kN„>Ú 2R! t  %  21Ú 2  t ! 1I
4Y! i i u3 Ú u1 Ú u

2

 t



1

u3 Ú u1 Ú u



1 1

I
4 I
4

2

u1 Ú



c!

bun

ä

(112a)

 

(112b) 2  t RkN„> 1 > 2
4Y! t  ! buä 1Ú 2

if n  1  2 The time at which the bond between u  0 1  2 t  and u  0 ! 1  2 t  breaks has been cho

sen to be t  0, so that by symmetry the time the bond between u  0 1  2 t  and u  1 ! 1  2 t 



breaks is 1  2
.

4.6 Equations solved in terms of a single mass point on crack line: Above the crack line, the equations of motion (112a) are completely linear, so it is simple to find the motion of every atom with n ê 1  2 in terms of the behavior of an atom with 35

n  1  2. Fourier transforming Eq. (112a) in time gives

h! € un €f, ib€ %



 %

2

Let

%

un €f& u1 Ú

2

       … Ú k        k          …Ú k

i gn 1 1 …Ú % i i eei>Ú ! 6 % e i>eÚ i ei gn 1 1…Ú % ei

á … €f €f  …€f

1 2 un 1 1 2 un 1 2 un 1

i

…€f ek  n  1Ú 2 

 !Ú  !

i gn 2

gn

1

gn

1

(113)



(114)



Substituting Eq. (114) into Eq. (113), and noticing that gn % gn

á 1  1 gives  1 k ei gn  1 gn  2 …Ú" 2! % ei gn  1 gn …Ú! 2! k    f €  u e á á

2 1Ú 2 !@€ 2u1Ú 2 €f, ib€ u1Ú 2 …€f %  % 12 u1Ú 2 €f ii ei>Ú ! 6 % e  i>Ú  k % 1 u1Ú 2 €f e  k ei gn  1 á gn  2 …"Ú  2! % ei# gn  1 á gn …Ú" 2! k  2 i (115) l € 2 % ib€ % 2 cosh  k  cos …€f> 2
42 % cos …€fI
4Y! 3  0 (116)  Defining 3 ! cos …€fI
4Y! € 2 ! ib€ z (117) 2 cos …€f 2
B



one has equivalently that

’

y  z%

with

z2 ! 1

y  ek

(118)

(119)



One can construct a solution which meets all the boundary conditions by writing un €f& u1 Ú

…€f e 



Ú

i gn 2

 y$

 % ! y $ y ! y

á Ú

N 1 2 n

 %

á Ú

N 1 2 n

k%

UN  n ! 1  2  N

2¨

¨ % € 2  (120) This solution equals u1 Ú 2 for n  1  2, and equals UN 2¨æ>a¨ 2 % € 2  for n  N % 1  2. The reason to introduce ¨ is that for n  N % 1  2, u  m n t ß UN . The Fourier transform of

2

i

N

N

2

this boundary condition is a delta function, and hard to work with formally. To resolve uncertainties, it is better to use instead the boundary condition uN

á Ú  t ,

UN e



u1 Ú

1 2

¦ññ t

(121)

and send ¨ to zero the end of the calculation. In what follows, frequent use will be made of the fact that ¨ is small. The most interesting variable is not u1 Ú 2 , but the distance between the bonds which will actually snap. For this reason define U  t +

u1 Ú

2

 t !

u 2

Ú  t 

1 2

2

 t %

u1 Ú 2  t % 1  2
4 2

Rewrite Eq. (112b) as u1å Ú

2

 t ,

1 2

% %

!

% u3Ú 2  t ! 1I
4 u3 Ú 2  t  u1 Ú 2  t % 1 Ê
BY! 4 u1 Ú 2  t  % u1 Ú 2  t ! 1 I
4 ! 2U  t ! 1 2
4r„> 1> 2
BY! t  2U  t r„BJ! t  36



 !

(122)

 bu1 Ú

ä

2



(123)

 Fourier transforming this expression using Eq. (120) and defining

&

U f…€f, now gives u1 Ú

2

ã

€f F €f!à 1 %



d € ei t U  t r„B

   U

ei>Ú

yN

³

1

N 1

N

N

Next, use Eq. (122) in the form

1%  U €f,

to obtain Writing

2

2

finally gives

U áf…€f Q €f

with

%

U



Q …€f+

…€f,



u1 Ú

U

UN Q0 Ñ

F …€fY!

UN N €

2

2¨

2

% ¨ 2



%

1

€f

…€f&"! %

U …€f+ U áf…€f

(124)

´

i

U €f F …€fY! 2  cos2 €f 4
« U

t

2z cos €f 2
4

 >Ú  

e



€f,"!

2

$  % ! y $  % ! y ! y

with F …€f&





(125) (126)

(127)

UN 2¨ N € 2%

¨ 2

…€f

(128) (129)

1 % 1 % ¨ i€ ¸ ¨ ! i€ Ô

F …€f 2 cos2 €f 4
0

(130) (131)

To obtain the right hand side of Eq. (130) one uses the facts that F  0 &"! 1  N, and that ¨ is very small, so that the right hand side of Eq. (130) is a delta function. The Wiener-Hopf technique [46] directs one to write Q Qá

Q €f&





…€f …€f

(132)

where Q is free of poles and zeroes in the lower complex € plane and Q á is free of poles and zeroes in the upper complex plane. One can carry out this decomposition with the explicit formula

&

i' £

Q ‘€f> exp lim

0

d€ j ln Q € j  ! i€ 2 ¤ i€

)(+*

ã

j

k

(133)

Now separate Eq. (130) into two pieces, one of which has poles only in the lower half plane, and one of which has poles only in the upper half plane: U á‘…€f Q0UN 1 ! Q áh€f Q  0  J! i€ % ¨Ð Q0UN 1 U …€f ! % ¨æ Q €f  Q  0   i€









37



(134)



Because the right and left hand sides of this equation have poles in opposite sections of the complex plane, they must separately equal a constant, . The constant must vanish, or U and U á will behave as a delta function near t  0. So

,

U



Q0 Q UN Q  0 ;a¨



€f,

and

Q0 Q ᑐ€f  0ºm¨¸! i€f

Q

U ᑐ€f, UN

€f % i€f

(135a)

(135b)



One now has an explicit solution for U …€f . Numerical evaluation of Eq. (133), and U  t  from Eq. (135) is fairly straightforward, using fast Fourier transforms. However, in carrying out the numerical transforms, it is important to analyze the behavior of the functions for large values of € . In cases where functions to be transformed decay as 1  i€ , this behavior is best subtracted off before the numerical transform is performed, with the appropriate step function added back analytically afterwards. Conversely, in cases where functions to be transformed have a step function discontinuity, it is best to subtract off the appropriate multiple of e t „> t  before the transform, adding on the appropriate multiple of 1 > 1 ! i€f afterwards. A solution of Eq. (135) constructed in this manner appears in Figure 15.



4.7 Relation between -

and .

Recall that making the transition from the nonlinear problem originally posed in Eq. (99a) to the linear problem in Eq. (112) relies on supposing that bonds along the crack line snap at time intervals of 1  2
. Because of the symmetries in Eq. (109), it is sufficient to guarantee that u  t , u f at t  0 (136)



All displacements are simply proportional to the boundary displacement UN , so Eq. (136) fixes a unique value of UN , and its dimensionless counterpart, . Once one assumes that the crack moves in steady state at a velocity
, there is a unique to make it possible. To obtain Eq. (136), one needs to require that t

£



lim 0

ã

d€ e 2¤

  U i t

€f,

uf



This integral can be evaluated by inspection. One knows that for positive t

(137)

ê

0,

 (138) i ! i€ t k U €f, 0 and that any function whose behavior for large € is 1  i€ has a step function discontinuity ã

d € exp

at the origin. Therefore, Eq. (137) and Eq. (135a) become uf



Since from Eq. (131) it follows that Q xp

UN Q0

&



Q Q

mp    0 

1, one sees from Eq. (133) that 38

(139)

Q



Qá

xp &

xp &

1

As a result, one has from Eq. (139) and the definition of Q



\



 0

Q0

(140)



given in Eq. (108) that (141)



To make this result more explicit, use Eq. (133) and the fact that Q R!@€f0 write Q

l

Q



¾

 0>

exp



exp

 0> ’

d € 1 ln Q …€  % ln Q J!h€ j j % i€ 2¤ 2 i Ó! i€ j d€ 1 Q €  % ln ³ j j 2¤ i ! 2i€ Q† € ´

*

¾ã ã

¾

Q0 exp

! ã

*

j

j

exp

*

2

j

d€ 1 Q € ln ³ j 2¤ 2i€ Q† €

 j ´ k 

j

i! ã

j

(142)

*

% €

Q …€ d€ 1 ln ³ j 2¤ 2i€ Q† …€

Placing Eq. (143) into Eq. (141) gives



j

 j kÀ

Q† €f to

j

2

ln Q  0 Ïk À

À j ´ 

(143)

(144)

j

In order to record a final expression that is correct not only for the Mode III model considered here, but for more general cases, rewrite Eq. (144) as



C exp

d€ j 1 2¤ 2i€

i! ã

Ý j

ln Q …€

j !

ln Q €

 k j Þ

(145)

where C is a constant of order unity that is determined by the geometry of the lattice, equaling 1 for the triangular lattice loaded in Mode III, but 2  \ 3 for a triangular lattice loaded in Mode I [38]. When written in this form, Eq. (145) is suitable for numerical evaluation, since there is no uncertainty relating to the phase of the logarithm. When b becomes sufficiently small, Q is real for real € except in the small neighborhood of isolated roots and poles that sit near the real € axis. Let riá be the roots of Q with negative imaginary part (since they belong with Q á ), ri the roots of Q with positive imaginary part, and similarly pi the poles of Q. Then one can rewrite Eq. (145) as



&



1C /00 2 r  pá rá p i

i

i

i

for b n

0



(146)

One may derive Eq. (146) as follows: away from a root or pole of Q, the integrand of Eq. (145) vanishes. Consider the neighborhood of a root r á of Q which falls to the real axis from the negative side as b n 0. For the sake of argument, take the imaginary part of this root to be ! ib. In the neighborhood of this root, say within a distance \ b, the integral to compute for Eq. (145) is

!

1 2¤

 á43 3 r

ãr

b b

d€ 2i€

j j

¾

ln €

!à r á ! ib2 ! ln …€ !à r á % ibr À  j

39

j

(147)

Defining €

jj

€ ! r á j

, and integrating by parts gives

! 21¤ ã -!

3

3

1 ln r á 2

d€ j j ln i € b 2i

b

%àé \

b

% j j rá k €



2ib

jj

2

%

(148)

b2

(149)

Similar integrals over other roots and poles of Q finally produce Eq. (146). Together with Eq. (135), Eq. (145) and Eq. (146) constitute the formal solution of the model. Since Q is a function of the steady state velocity
, Eq. (144) relates the external driving force on the system, , to the velocity of the crack
. The results of a calculation appear in Figure 14. 0.8 Steady states go unstable 0.6 c

65

0.4

Slowest steady state

0.2

Velocity gap

0.0

Linearly stable lattice-trapped states

1.0

1.1

1.2

1.3

1.4

Figure 14: The velocity of a crack
4 c, scaled by the sound speed c  \ 3  2, is plotted as a function of the driving force . The calculation is carried out using Eq. (146) for N  9.

4.8 Phonon emission: Right at  1 just enough energy is stored to the right of the crack tip to break all bonds along the crack line. However, all steady states occur for ê 1, so not all energy stored to the right of the crack tip ends up devoted to snapping bonds. The fate of the remaining energy depends upon the amount of dissipation b, and the distance from the crack tip one inspects. In the limit of vanishing dissipation b, traveling waves leave the crack tip and carry energy off in its wake; the amount of energy they contain becomes independent of b. Such a state is depicted in Figure 15, which shows a solution of Eq. (135) for
â 0 5,  N  9, and b  0 01. For all nonzero b, these traveling waves will eventually decay, and the  extra energy will have been absorbed by dissipation, but the value of b determines whether one views the process as microscopic or macroscopic. The frequencies of the radiation emitted by the crack have a simple physical interpretation as Cherenkov radiation. Consider the motion of a particle through a lattice, in which 40

5.0 4.0

u(t)

3.0 2.0 1.0 0.0 -50

-25

0

25

50

t Figure 15: A plot of U  t  for
Î 5, N  9, and b  0 01, produced by direct evaluation of   Eq. (135). Note that mass points are nearly motionless until just before the crack arrives, and that they oscillate afterwards for a time on the order of 1  b.

798 :

3.0 2.5

798 :

2.0

2

1.5

k

k

@ ?= > <;

5

1.0 0.5 0.0

0

2

6

4

8

10

k Figure 16: Graphical solution of Eq. (153), showing that for low velocities, a large number of resonances may be excited by a moving crack. 41

the phonons are described by the dispersion

€¯¦¡ k 

(150)

If the particle moves with constant velocity
 , and interacts with the various ions according to some function , then to linear order the motions of the ions can be described by a matrix D which describes their interactions with each other as

A

B

muå l

"!

©C

DFE

l



D

B

C

C

©  R l ! R l   ul  %

l



A B  R  !¿
 t   l

(151)

Multiplying everywhere by eik R , summing over l, letting K be reciprocal lattice vectors, and letting be the volume of a unit cell gives

G

l

©C

B

muå  k Y

D

B

C

C

 k u  k %

HD E  B G E eE á A 1

©

i

k K t

K

 k % K  

(152)

Inspection shows that the lattice frequencies excited in this way are those which in the extended zone scheme[ [47]] obey

€Î k&²
‘y k 

(153)

Pretending that the crack is a particle, one can use Eq. (153) to predict the phonons that the crack emits. There is another version of this argument that is both simpler and more general. The only way to transport radiation far from a crack tip is in traveling waves. However, in steady state, the traveling waves must obey symmetry (109a). In a general crystal with lattice vectors R and primitive vectors a , applying this requirement to a traveling wave exp ik y R ! i€ t k gives

i

$ E D E á E    E I á Ú !J%  e $ E D E   E  % 

e ik

R na

i

k t na

ik R i

kt

(154)

Assuming a and
 are parallel, Eq. (153) results again. There are two phonon dispersion relations to consider. One gives the conditions for propagating radiation far behind the crack tip, and the other gives the conditions for propagating radiation far ahead of the crack tip. Far behind the crack tip, all the bonds are broken. Finding waves that can travel in this case is the same as repeating the calculation that led to Eq. (128), but with U set to zero, since all bonds are broken, and with UN  0, since phonons can propagate without any driving term. Examining Eq. (128), one sees that the condition for surface phonons to propagate far behind the crack tip is F €f, 0. Similarly, far ahead of the crack tip no bonds are broken U should be set equal to U , and the condition for phonons is F …€f>! 2 cos2 €f 4
c 0. According to Eq. (131), the roots and poles of Q …€f are therefore the phonon frequencies behind and ahead of the crack, and these are the quantities appearing in Eq. (146). I do not know if Eq. (146) is more than approximately correct for particle interactions more general than ideally brittle bonds.





42

5.0

Displacement U

4.0 3.0 2.0 1.0 0.0 -50

-25

0

25

50

Time t

Figure 17: Behavior of U  t  for
Î 2, b  0 05, and N  9. Notice that at t  0, indicated   by the dashed line, u is decreasing, and that it had already reached height 1 earlier. This state is not physical.

4.9 Forbidden velocities After making sure that bonds along the crack line break when they are supposed to, it is necessary to verify that they have not been stretched enough to break earlier. That is, not only must the bond between u0 and u0 reach length 2u f at t  0, but this must be the á first time at which that bond stretches to a length greater than 2u f . For 0 Ü×
ÒÜ 0 3 (the  r2 precise value of the upper limit varies with b and N) that condition is violated. The states have the unphysical character shown in Figure 17. Masses rise above height u f for t less than 0, the bond connecting them to the lower line of masses remaining however intact, and then they descend, whereupon the bond snaps. Since the solution of Eq. (112) is unique, but does not in this case solve Eq. (99a), no solutions of Eq. (99a) exist at all at these velocities. Once the crack velocity has dropped below a lower critical value, all steady states one tries to compute have this character. This argument shows that no steady state in the sense of Eq. (109a) can exist. It is also possible to look analytically for solutions that are periodic, but travel two lattice spacings before repeating. No solutions of this type have yet been found. Numerically, one can verify that if a crack is allowed to propagate with right above the critical threshold, and is then very slowly lowered through the threshold, the crack stops propagating. It does not slow down noticeably; suddenly the moving crack emits a burst of radiation that carries off its kinetic energy, and stops in the space of an atom. That is why Figure 14 shows a velocity gap.



4.10 Nonlinear Instabilities It was assumed in the calculations predicting steady states that the only bonds which break are those which lie on the crack path. From the numerical solutions of Eq. (135), one can test this assumption; it fails above a critical value of . The sound speed c equals 43

\ 3 2, and velocities will be scaled by this value. For N  9, at a velocity of
c  c  666 , c  1 158 , the bond between u  0 1  2  and u  1 1  2  reaches a distance of 2u 2r  2r

f some short time after the bond between u  0 1  2  and u  0 ! 1  2  snaps. The steady

state solutions strained with larger values of are inconsistent; only dynamical solutions

more complicated than steady states, involving the breaking of bonds off the crack path, are possible. To investigate these states, one must return to Eq. (99) and numerically solve the model directly. These simulations have been carried out [37, 48] and some results are contained in Figure 18. These theoretical calculations were done in response to experiments [49–51, 51] that showed cracks The geometry of the instability is the same as the geometry of instabilities found experimentally in Plexiglas, shown in Figure 19. However, a detailed connection between the theory, which describes the onset of branching in a crystal, and the experiment, which observes mesoscopic branching in an amorphous polymer, has not fully been worked out.

K L 1 N 147 M OQP c L 0 N 645 K L 1 N 165 M OQP c L 0 N 630 K L 1 N 376 M OQP c L 0 N 624 K L 1 N 835 M OQP c L 0 N 775 Figure 18: Pictures of broken bonds left behind the crack tip at four different values of . The top figure shows the simple pattern of bonds broken by a steady-state crack. At a value of slightly above the critical one where horizontal bonds occasionally snap, the pattern is periodic. All velocities are measured relative to the sound speed c  \ 3  2. Notice that the average velocity can decrease relative to the steady state, although the external strain has increased. As the strain increases further, other periodic states can be found, and finally states with complicated spatial structure. The simulations are carried out in a strip with half-width N  9, of length 200 and b  0 01. The front and back ends of the strip have  short energy-absorbing regions to damp traveling waves. The simulation was performed by adding unbroken material to the front and lopping it off the back as the crack advanced. The diagram shows patterns of broken bonds left behind the crack tip. Just above the 44

Figure 19: Images of local crack branches in the x ! y plane in PMMA. The arrow, of length 250e m, indicates the direction of propagation. All figures are to scale with the path of the main crack in white. (top)
ÒÜà
c (center)
 1 18
c (bottom)
 1 45
c . (from [52])





threshold at which horizontal bonds begin to break, one expects the distance between these extra broken bonds to diverge. The reason is that breaking a horizontal bond takes energy from the crack and slows it below the critical value. The crack then tries once more to reach the steady state, and only in the last stages of the approach does another horizontal bond snap, beginning the process again. This scenario for instability is similar to that known as intermittency in the general framework of nonlinear dynamics [53] ; the system spends most of its time trying to reach a fixed point which the motion of a control parameter has caused to disappear. Heizler, Kessler, and Levine [43] have shown that when the snapping bond is replaced by a rounded potential, the steady crack states are linearly unstable in a more conventional way, the the instability is a Hopf bifurcation.

4.11 The connection to the Yoffe instability The basic reason for the branching instability seen above is the crystal analog of the Yoffe instability, working itself out on small scales. The Mode III calculation finds that the critical velocity for the instability to frustrated branching events is indeed close to the value of 0 6cR predicted by Yoffe in the continuum. The critical velocity seen experimentally in  amorphous materials is 1/3 of the wave speed, not 2/3. This discrepancy could be due to some combination of three factors. 1. The force law between atoms is actually much more complicated than ideal snapping bonds. Gao [54], has pointed out that the Rayleigh wave speed in the vicinity of 45

a crack tip may be significantly lower than the value of cR far away from the tip because material is being stretched beyond the range of validity of linear elasticity. 2. The experiments are at room temperature, while the calculations are at zero. 3. The experiments are in amorphous materials, while the calculations are in crystals.

4.12 Generalizing to Mode I: Let u i describe the displacement of a mass point from its equilibrium location. Assume that the energy of the system is a sum of terms depending upon two particles at a time, and linearize the energy to lowest order in particle displacements. Translational invariance demands that the forces between particles 0 and 1 depend only upon u 1 !àu 0 , which will be defined to be  1 . However, the force can be a general linear functional of  1 . A way to write such a general linear functional is to decompose the force between particles 0 and 1 into a component along dè 1 ï , and a component that is along dè 1 ð . The first component corresponds to central forces between atoms, while the second component is a non-central force. Non-central forces between atoms are the rule rather than the exception in real materials, as first appreciated by Born. Their quantum-mechanical origin is of no concern here, only the fact that they are not zero. Suppose that the restoring force parallel to the direction of equilibrium bonds is proportional to ì ï , while that perpendicular to this direction is proportional to ì ð . The force due to the displacement of the particle along  1 "u i 1 ¥ j 1 !¢u i ¥ j is

R

R



á

R

ì ï dè ï 1  R  1 y dè ï 1  % ì ð dè ð 1  R  1 y dè ð 1 

²ì ï  ! 21 \ 23 º  ! 21 R x1 \ 23 R y1  % ì ð  \ 23 12 ; \ 23 R x1 12 R y1 

(155) (156)

Adding up contributions from other particles in this way we get for the force due to neighbors

©6 © è R  j  m ny dè q j  ì (157) q dq j 

 j ª 1 q ª ïJ¥ ð By varying the constants ì ï and ì ð , one can obtain any desired values of shear and F  m n ,

longitudinal wave speeds, which are given by c2l  ct2 

3a2 aì 8m 3a2  3ì 8m

ð % 3ì ï  ð % ì ï

(158a) (158b)

where m is the mass of each particle. In addition to forces between neighbors, it is possible to add complicated dissipative functions depending upon particle velocities. In Ref. [38], a term was added to the equations so as to reproduce the experimentally measured frequency dependence of sound attenuation in Plexiglas. There is a slight technical restriction in the calculations of Ref. [38]; right on the crack line, forces are required to be central. 46

Microcracking instability

0.6

Minimum velocity steady state

Height=80 Height=160

ct

65

0.8

0.4 Velocity gap 0.2

Lattice trapping

0.0 1.0

1.2

1.4

1.6

1.8

2.0

Figure 20: Crack velocity versus loading parameter  KI  KIc for strips 80 and 160 atoms high, in limit of vanishing dissipation, computed for Mode I crack by methods of following section. will be defined precisely in Eq. (108). The spring constants have values ì ï  2 6 and ì ð  0 26, so that the non-central forces are relatively small. The   fact that the results have become independent of the height of the strip for such small numbers of atoms in the vertical direction suggests that relatively small molecular dynamics simulations can be used to obtain results appropriate for the macroscopic limit. E. Gerde has found a way to overcome this technical limitation, and the results do not change noticeably. There is no universal curve describing Mode I fracture. Many details in the relation between loading and crack velocity depend upon ratios of the sound speeds, and upon the frequency dependence of dissipation. In the limit of central forces, ì ð  0, it turns out to be difficult but not impossible to have cracks in steady state. This means that the range of loads for which cracks can run in a stable fashion is small, and depends in detail upon the amount of dissipation. Crack motion is greatly stabilized by having nonzero ì ð . Figure 20 displays a representative result with some implications for the design of molecular dynamics simulations. In the limit of vanishing dissipation, the relation between crack speed and dimensionless loading  KI  KIc becomes nearly independent of the number of vertical rows of atoms for strips around 40 rows high. The definition of is that it describes the vertical displacement imposed upon the top and bottom of the strip, but scaled so that  1 when these boundaries have been stretched just far enough apart so that the energy stored per length to the right of the crack tip just equals the energy cost per length of extending the crack. It is defined precisely by Eq. (108) Analytical work has continued on these classes of models, and there are many results I will not refer to much. Kessler and Levine have conducted the most careful studies of lattices with nonlinear force laws and dissipation acting in concert [43, 44]. Their findings emphasize the point that crack dynamics are terribly sensitive to features of the tip that do appear in any continuum form of the problem. With Eric Gerde, I applied some of these ideas to interface fracture and to friction [45]. 47

5 Scaling ideas and molecular dynamics Because fracture involves both large and small scales, because it is of practical importance, and because it became fashionable for reasons that are hard to explain precisely, there is an ongoing effort of many groups to study fracture through numerical methods that reach down to the atoms. These studies fall into roughly two groups: 1. Calculations employing classical molecular dynamics. The largest studies employ a billion atoms or more [55–59]. 2. Calculations employing continuum mechanics, but with special spatial regions where the analysis descends to molecular dynamics, and perhaps even, in small regions, quantum mechanics [60–65] My own prejudice is toward carrying out molecular dynamics in the spirit of fracture mechanics. That is, instead of trying to simulate a macroscopic sample in its entirety in the midst of one computer program, I prefer to use computer calculations to tabulate quantities that linear elastic fracture mechanics needs in order to make predictions. If one assumes that cracks move along a straight line, then what one then has to compute is fracture energy as function of crack velocity. In what follows, I will show how this task can be accomplished with numerical systems of relatively small size. The essential idea is to abstract from the analytical solutions in lattices two scaling laws, one for space and one for time. The time scaling law is easy. Simply look for steady states. Once one finds numerically a steady state obeying the symmetry (109a), one knows what it does forever, and extrapolating up to time scales of milliseconds is no problem. The spatial scaling is slightly trickier. In order to relate samples of different size to one another, let Gc be the Griffith energy density [12,66]; that is, twice the crack surface energy density, a lower bound on the energy per unit area required for a perfectly efficient crack to propagate along a certain plane. One can then define a dimensionless measure of loading to be  ’ G  Gc (159)

where G is the fracture energy density already introduced, i.e. the elastic strain energy stored per unit area (in the fracture plane) ahead of the crack. Analytical solutions for the ideal brittle solid show that the relationship between and crack velocity becomes independent of the height of the strip (number of planes stacked vertically) for surprisingly small strips [67]; a strip 80 atoms high has for all practical purposes reached the infinite limit, Fig 20. Guidance for conducting computationally expensive, moderately realistic, molecular dynamics simulations can be had from this simple result: The very rapid convergence of the main quantity of physical interest allows one to obtain physically meaningful results from simulations that are considerably smaller than many being carried out these days . This approach can also be viewed as an alternative to methods that join together atoms and continua [60–65]. If computational resources were infinite, one would still need to choose properly the sample geometry in order to compare with and inform experiment and theory. With limited computational resources, choice of sample geometry and scale are of paramount importance: when performing molecular dynamics simulations of fracture, one should not ask 48

how large a system one can work with, but rather what is the smallest system that will give results scalable to the macroscopic level. A good scaling argument can transform a ‘grand challenge’ problem into something that’s better focused and less of a challenge. Size does matter: a smaller spatial scale allows one to follow the evolution of systems for longer times; and with a strip geometry one can then analyze fully steady state behavior and connect with theory and experiment.

5.1 Molecular dynamics The molecular dynamics method applies to any system of particles with some prescribed inter-particle potential. It is based upon the following equation: F



ma (

(160)

that is it consists of integrating Newton’s equations of motion for all particles in lock step over a series of time steps, the size of the step being chosen small enough to give converged dynamics. Time step integration is often done using the Verlet algorithm [68–70]. Lattice theory makes predictions that are hard to observe in experiment, and one of the reasons for doing computer simulations is to relate the two. The simulations I will describe here in detail were carried out in silicon. Silicon is extremely brittle, and highquality macroscopic single crystal wafers are cheap. It is therefore an excellent candidate for laboratory fracture experiments. Silicon is also technologically of great importance, and as a result is one of the most studied materials. One measure of this is that there are over thirty effective (although rather ineffective) interatomic potentials for silicon in the literature [71, 72]. Furthermore, transmission electron micrographs of cracks in silicon wafers reveal atomically sharp crack tips [13,73]. Silicon therefore is an obvious candidate for molecular dynamics investigations of dynamic fracture, and an appropriate setting for testing atomic–level understanding. If there is any brittle element for which one ought to have the knowledge base making it possible to calculate fracture properties in detail, and the motivation to carry the work through, it would be silicon.

5.2 Interatomic potentials The equation underlying materials physics is not in doubt. It is Schr¨odinger’s equation. This equation can be solved analytically for the hydrogen atom, numerically for the helium atom, and with reliable approximate methods of quantum chemistry for small molecules. For any solid of interest in the study of materials, all hope of controlled approximations must be abandoned, and the Schr¨odinger equation is brutally reduced to tractable form in a way that is refined by comparison with experiment. Quantitative methods that employ such approximations from the beginning are called ab initio. Despite their origin, the ab inito methods are the most reliable technique available for numerical treatment of materials. However, they are restricted to perfect crystals, and become unmanageable when the unit cell contains more than around 1000 atoms. The main motivation, then, for constructing empirical or effective ‘classical’ interatomic potentials is speed of computation and the ability to work with relatively large numbers of particles. 49

With effective potentials, 107 atoms followed for a few tens of nanoseconds is achiev able. This difference in computational scales becomes important when modeling processes which require a minimum of 105 atoms to capture just some of the complex underlying  physics — processes involving fracture, dislocation loops, grain boundaries, or amorphousto-crystal transitions, for example. Another motivation for constructing effective interatomic potentials is that they make the complex physics of what are fundamentally quantum mechanical phenomena more physically intuitive, so that one may interpret the results of atomistic simulations in terms of simple principles of chemical bonding [72].

5.3 Realistic potentials for silicon? Solid silicon is covalent and has the open diamond crystal structure. If only two-body potentials operated among the atoms, one would expect the crystal to collapse in on itself to form a close-packed structure, thereby reducing its energy. Covalent systems, however, are characterized by restoring forces between pairs of contiguous interatomic bonds. That is, pairs of bonds with an atom in common want to maintain a preferred angle between them. These extra forces are what stabilize the open diamond structure in silicon, carbon, grey-tin, and germanium, for example. The lowest-order way of capturing this property with effective interatomic potentials is to go beyond binary bonding and include a threebody term in the system’s Hamiltonian. The Stillinger-Weber (two- and three-body) potential [74] has proved to be very popular and durable in the literature. It gives excellent elastic properties, and captures well the nonlinear physics involved in heating and melting. It therefore is a reasonable starting point for conducting molecular dynamics fracture simulations in silicon. Unfortunately,

WYSkU

l!T

s urt ommq pn VkU

ijT S"U S!T

T"U

V"T WXU!U Z[Y\^]"_a`cbehdf g

WYS!T

Figure 21: Crack tip blunting in Stillinger-Weber silicon: two dislocations open up at the tip, preventing it from advancing. The crack is pointing in the direction 1 1† 0k in the plane i (110), and the system is loaded with a strain parameter  1 6.



50

the potential will not yield fracture along the experimentally preferred fracture planes  111  and  110  : At low or moderate strains a crack will not move at all. What happens is that two dislocations open up at the crack tip, blunting it and preventing it from advancing, Fig 21. One can play around with giving a transverse opening velocity to a select few atoms around the tip. But to no avail. The crack simply will not crack. At very high strains, the crack tip region melts. Stillinger-Weber does give a type of fracture along the  100  plane which is quite rough on the atomic scale, Fig 22. Abraham et al call this brittle fracture [58]. There is as yet little consensus on a precise definition of brittleness. The experimental results of Lawn and Hockey [13,73] for fracture along  111  in silicon show, however, that it is possible to have atomically sharp fracture, i.e. where the newly created fracture surfaces are atomically flat. The experimental evidence for fracture along  100  , on the other hand, is scant and inconclusive [75], and that Stillinger-Weber yields fracture along this plane, albeit in a rough manner, might even be yet another indication of the potential’s shortcomings. It is possible to get cracks going, in an ideal brittle manner, along  111  and  110  with the Stillinger-Weber potential by increasing the restoring forces between pairs of bonds, i.e. by increasing the stability of the tetrahedra in the diamond lattice, and thus making the crystal more brittle. This can be done by scaling a parameter, ± , the coupling constant in the three-body term. Once ± has been changed in this way, the potential acquires a new name, and is known as the Inadvertently Modified Stillinger-Weber potential (IMSW). Originally, Stillinger and Weber set ±U 21 (dimensionless). However, by doubling this, as in the pioneering inadvertent modification, one can obtain fast brittle fracture. With a fast crack running, if one quasistatically decreases ± , the crack arrests well before one reaches ±z 21, Fig 23. Although ±z 42 gives fracture phenomenology in reasonable accord

ŒŒ  ’

‰ ‹ˆŠ a‘  ƒ…‡ † ƒ „ a Ž Œ  ŒIaŽ

ŒIa

ŒI‘a ŒI’ v w xXyaz|{~}k€‚ 

Œ Ž“Œ 

Œ ŽŽ

Figure 22: A two-dimensional projection of rough cracking along the (100) plane in Stillinger-Weber silicon. The system is loaded with a strain parameter  1  7, corresponding to a fracture energy density of almost three times the Griffith energy density. The average crack speed is 1.9 km/s. 51

– F§ ¨Y©4§ © •k•” ÌÊ¹Ë ¥ ÈÀ·É ÅÆ Ç ¢ ýĠÀ¯½¿Á • ½¾ ” ¢Y£ ¢X¤ ¥k” ¥k¢ ¥ – ¦ – §«ª^©4§ ¢j•Y• © ÌÊ¹Ë ¥ ÈÀ·É Å Æ ½Ç ¢ ÃÄ À¯½¿Á • ½¾ ” ¢Y¢ ¢ – ¢X£ ¢X¤ ¥k” ¥"¢ ¦

—™—™˜›˜Ÿž š|œœ ¡ ¥k£ ¥k¤ – ” – ¢

—™—™˜›˜›š|š|œœ ¶ ¬·¬ ­¯³¹®±¸?° º ²»³ ¬·­¯®±°¼² ³ —™˜Ÿž œ ž ´ µ ¬­¯®±° ² ³ ¥ – ¥Y£ ¥k¤ – ” – ¢

Figure 23: Crack velocity profile (a) along (110) 11† 0k , and (b) along (211) 11† 1† k in modii i fied Stillinger-Weber silicon: ± scales the strength of angular forces between pairs of bonds that stabilize the diamond lattice structure. Stillinger and Weber’s value for ± is 21. However, leaving all other parameters unchanged, cracks will not propagate unless one uses a larger value of ± . specifies the loading. with experiment, it has the adverse effect of raising the melt temperature above 3500 K, whereas experimentally the melt happens at 1685 K. The Young moduli also get shifted. Moduli results from tensile tests on small samples, for three lattice directions, are given in Table 2 for both the original and modified Stillinger-Weber (SW) potentials, along with the corresponding experimental values. The story of searches for potentials after learning about the Stillinger-Weber potential and its inadvertent modifications is long and alternately depressing and amusing. Dominic Holland and I devoted a fair amount of attention to a potential of Kaxiras and Bazant [72,76–78], and then one of Chelikowsky and Phillips [79]. We learned that these potentials and many others have a defect that can easily be understood in a qualitative way. Imagine taking a uniform crystal and simply expanding the lattice constant. This is a 52

computation that can be done easily by density functional theory, and one can compare the results with the same calculation repeated with a classical potential. Then repeat the computation putting the crystal in uniform shear. Most classical potentials have too large a peak when the crystal is expanded, and resist shearing too little. The height of the peak resisting expansion can be as much as 5 times larger than the corresponding peak predicted by density functional theory. Getting these features of atomic force laws wrong is terrible if one wants to predict fracture properties at all accurately. If shearing is easier than pulling apart in tension, a crystal will emit dislocations by shearing at the tip rather than breaking in tension. In searching for a potential that would correctly reproduce the qualitative features obtained from density functional theory, I eventually settled upon the Modified Embedded Atom Method of Baskes [80, 81]. This potential is based upon a number of physical ideas that I find quite persuasve, particularly the idea of screening. The density functional results for uniform expansion are built into the potential from the beginning, Without any tuning or adjustment, it does an acceptable job of describing the resistance to shear, and in computer simulations it does permit the motion of brittle cracks. Still, there are some arbitrary features of the potential that are a bit unnerving, and from a quantitative point of view it does not turn out in the end to be any better than the IMSW with which our studies began.

5.4 Results of zero Kelvin calculations in silicon I will now describe the results of performing simulations on 100-200,000 atoms in silicon. Most of the calculations are with the modified Stillinger-Weber potential, but toward the end I also quote some results from MEAM. Questions to be answered: (1) Are there loads where cracks are attracted to steady states? (2) Do cracks emit phonons at the predicted frequencies? (3) Do cracks refuse to travel below a minimum velocity
1 ê 0?, and (4) Do they go unstable above an upper load c ? The answer to all questions is yes. Fig 25, however, which shows atomic motions after the crack has been traveling for over 0.24 ns, as anticipated by the theory of ideal brittle fracture the crack has reached steady state with velocity
=3460 m/s, which means that the vertical displacement z R of an atom originally at crystal location R is related to the vertical displacement zR na of an atom á n lattice spacings a axè to the right by

E E

Í

E á E t%

zR

E100 E110 E111

na

E

na I
4, zR  t 



E

(161)

Silicon Original SW IMSW GPa 130 114 172 GPa 169 139 189 GPa 188 151 201

Table 2: Elastic constants of silicon, comparing Stillinger-Weber, modified StillingerWeber potentials, and experiment. 53

6 5

Bazant-Kaxiras

˚ Force (eV/A)

4 3

Stillinger-Weber

2

Density functional

1 0 1.0

(A)

10

Energy per atom (eV)

8

Î

a a0

1.4

1.6

Density Functional Theory Modified Stillinger Weber Original Stillinger Weber Modified Embedded Atom Method

6 4 2 0 1

(B)

1.2

1.5

2 Extension Ratio

2.5

3

Figure 24: (A) Comparison of density functional theory, EDIP, and the Stillinger-Weber potential for the restoring force during uniform expansion. MEAM is identical to the density functional result by construction. (B) Same as (A), but now repeated for uniform shear. For a range of loads , Eq. 161 applies for any pair of atoms, whatever their separation along the crack surface. In order to obtain the perfect periodicity shown in Fig. 25, the crack was allowed to run first for 60,000 time steps so as to come into equilibrium with the waves it sends towards top and bottom boundaries. The longer and/or higher the system the longer it will take to reach a steady state. The 240 ps required to reach the steady state depicted in Fig 25 is about an order of magnitude longer than the duration of most large-scale molecular dynamics simulations of fracture. Using a reduced spatial system size, as validated by the scaling argument, is what makes this possible. This becomes crucially important when proceeding through a sequence of steady states, which requires quasistatic changes in the loading, or G. To obtain a full set of results, like Fig 26, a crack actually will travel tens of microns, or for times on the order of a tenth of a microsecond. To achieve this, one needs not only efficient code and a high performance computer, but also a physically motivated smallest computational cell — a minimum thin strip on a conveyor belt.

54

Steady State after 0.24 ns at ∆=1.6 vertical profile of two atoms 160 Angstroms apart

5

vertical displacement (Angsrtoms)

4

3

2

1

0

1

2

3 time (ps)

4

5

Figure 25: Crack that has been traveling for over 0.24 ns. That the two overlapped  curves are almost completely indistinguishable shows the crack has reached steady state, according to Eq. 161, and is emitting phonons in accord with Eq. 153.

5.5 Along Ï 110 Ð

The relation between velocity 
and load for cracks exposing  110  and traveling along k is shown in Fig. 26. The crack velocity smoothly decreases as decreases, until at i
@11† 02256 m/s and  1 258, the crack abruptly comes to a halt. Raising again, the crack  does not begin to move until  1 366, a value that is sensitive to residual vibrations in  the crystal, but the rising curve then perfectly overlaps the descending one. Crack speed continues to rise smoothly until
Î 3586 m/s,  2 2, at which point steady state motion  becomes unstable. When a crack goes unstable, complicated phenomena such as formation of small branches, emission of dislocations, changes in the plane of propagation can occur, and intermittency where the crack makes repeated attempts at branching.

5.6 Along Ï 111 Ð : Crackons

The relation between velocity
 and load 011† k is shown in Fig. 27. For 1 44 Ü

i



Ü

for cracks along  111  and traveling along 2 2, the crack has stable steady states, and



55

ê

ÖF×|رټÚFÛÜ Ú¹ØIÖFÚÝÖßÞIÚFÖßà

ûúø ù é öî ÷ Õ ñîóô ëðòõï ëëìí Ñ Ó Ñ±Ò Ó

á

Ñ±Ò Ô

ØIâ·ã¹äæå ç Ú

Õ ÒÓ è

Õ ÒÔ

Figure 26: Relation between crack speed
and loading for crack along  110  at 0 K. As descends, velocity drops abruptly to zero at a lower critical value, and as ascends resumption of crack motion is hysteretic. For convergence check, system size was doubled along x and z, and
  measured.  





üký"þæÿ






  #

" " !





   



 

 



 



  

 

 





  



 





Figure 27: Relation between crack speed
and loading for cracking along  111  01 1† k at i 0 K. Dotted lines indicate forbidden velocities. The lower figure shows crackons for 1 18 Ü  Ü 1  44: the crack is able to expose many%$different states lying along many hysteresis loops. Ideal steady states are unstable above 2 2.



56

for ê 2  2 it goes unstable in a similar manner to cracks along  110 . However, for 1 175 Ü 1 44 the dynamics of the crack exhibits a number of interesting features that have notÜ been  seen previously, and for which there is not yet a complete theoretical description. There is a variety of different dynamical states available for each value of , where the crack travels at different speeds. Each of these states corresponds to a plateau in
  ; can change by as much as one fifth of the amount needed to go from arrest to instability and the crack velocity does not alter within numerical resolution. When the crack finally decides to accelerate out of the plateau, it may jump by over 1 km/s and reach an upper plateau to within a few m/s. On cyclical loading the same plateaus are always reached. All of these transitions are hysteretic, as depicted in Fig. 27. The different states emit noticeably different phonons; on a given plateau, the phonon frequencies appear fixed and their amplitude changes, while between plateaus the frequencies change in accord with Eq. 153. This is crackon behavior. All these phenomena are easily disguised if strain rates are too high. Resolving all the fine structure visible in Fig. 27 required Ü 8 e s 1 , or ä hz  c Ü 10 5 .

*ä

*





crack speed (km/s)

4

unstable

3

2

1

0

1.2

1.4

1.6 &

1.8

2.0

2.2

Figure 28: Crack velocity profile for Si '$ (111) 21† 1† k at 0 K. Dashed lines indicate forbidden i velocities. The crack is unstable above 1 9.



5.7 Crack behavior at nonzero temperatures in silicon Temperature implies energy fluctuations in time, so that if the crack gets trapped due to an energy fluctuation that reduces the fracture energy, further kinetic fluctuations may subsequently enable the crack to move on, giving rise to the possibility of creep. This suggests that the strength of lattice trapping should be a function of temperature. A series of simulations at 50 Kelvin intervals were carried out to investigate the effect of temperature on lattice trapping along  111  011† k . The results are shown in Fig 29. i Note also in Fig 30 that as is decreased below 1 the crack actually heals and travels backwards. This is perfectly reasonable: there hasn’t been any nonuniform surface damage, and no oxide layer has formed.

57

/ )

476

4 )

82496

4 L

IKJ G

>H

)

).404:6

)

)8+4:6

)

/2404:6

)

/+8+4:6

4 F

C D ;E A,@B 4

>? ; = ;< 4

4 3+404:6 ) 4

)+*,).-

)+* /+/

)0* /21

)+* 3+4

)+* 325

(

Figure 29: Velocity gap vanishes near 200 K. Dotted lines indicate forbidden velocities. All low-lying velocity states exist at room temperature. Cracking is along (111) 01 1† k .

i

5.8 Comparisons with experiment I turn to a number of questions that allow a comparison of theory and experiment. First, is there a velocity gap in brittle amorphous materials at room temperature? The answer seems pretty clearly to be “no;” data to this effect appear in Figure 31. There is no detailed dynamical theory for the motion of cracks in very brittle amorphous materials, so the experiments are the best guide to what happens. However, the point of performing numerical work in crystalline silicon was the possibility of carrying out experiments to compare directly. The first were done by Hauch et al. [82], with additional experiments by Cramer, Wanner, and Gumbsch [83], and more recently by Deegan, who has been working with me to obtain data at cryogenic temperatures. Room-temperature runs were performed along  111  1† 1† 2k in order to allow direct comi parison with experiment [82] (the zero Kelvin runs for this lattice direction are in Fig 28). ½ ½ ˚ 3, The results are in Fig. 32, and were obtained for a thin strip of size 532 15 154 A periodic along the thin axis. As before, new material was added ahead of the crack tip and ˚ of the old material lopped off at the tail every time the crack advanced to within 200 A forward end of the strip. In this fashion, the crack traveled 7 e m during the course of the simulations as G was varied between 5 and 14 J/m2 . The room-temperature experiments for fracture in silicon wafers along (111) 1† 1† 2k are described in [82]. These experiments i are difficult to perform because of the high Young’s modulus and brittleness of silicon. 58

S

`a.b [2cd+d2^Ufe^] ag c Whci^ZV o

T2U2VXWZY\[2]^

‚

€ }

si~

Q

P.N m

| y z p{

wixv

‚

sut p

P.N M

€ P }

pqr

si~ |

M2N m

M M2N M

b ^jW b ^YkWhcU+l

n

P

M+N O

PN Q

P.N Q\S

P.N Q O _

P.N R _

P.N oQ

QN M

Q2N S

Figure 30: Crack velocity profile with respect to quasistatic loading silicon at 300 K. All low-lying velocity states exist. 300

along  111  01 1† k in

i

PMMA

200 100 0

(m/sec)

300

„

200

Velocity

0

100

100

200

300

Homalite-100

0 0

100

200

300

400

500

75 Glass

50 25 0 0

400

800

1200

1600

Time t ( ƒ sec)

Figure 31: Velocity records for experiments in the double cantilever beam configuration. In PMMA and Homalite-100 crack arrest was achieved. In both materials the crack appears to decelerate smoothly but with increasing deceleration until it arrests. In glass crack arrest was not achieved in this configuration. However there is no sign of an initial velocity jump. This can be attributed to the sharp seed cracks that can be generated in glass. 59

4

4

3 2

0

0

2

4

6

8

10

12

14

16

(km/s)

Length [cm]

6

18

Velocity [km/s]

(A)

2

…

4 3

Experiment Molecular Dynamics

1

2 1

0

0 0

2

4

6

8

10

12

Time[µs]

14

16

18

0

(B)

5 10 15 Energy flow to crack tip G [J/m2]

20

5

Velocity [km/s]

4

3

2

Experiment (111, Hauch) Experiment (110, Cramer) Experiment (111, Deegan, ramp)

1

0 0

20

40

60

80

2

G [J/m ]

(C) Figure 32: Experimental and numerical determination of the crack velocity
as a function of fracture energy G along  111  1† 1† 2k at room temperature. The fuzz indicates the spread i in simulation velocity data. Experimental data from Hauch et al. [82], Cramer et al. [83] and R. Deegan.

60

The highest experimental and numerical crack velocities shown in Fig. 32 are reasonably close, but the minimum fracture energies at which a crack propagates differ: 2.3 J  m 2 in the experiments and 5.2 J  m2 in the simulations. Since the scale of crack velocities in a material is bounded by sound speeds [12], which are given correctly by the StillingerWeber potential, it is not surprising that the experimental and computational crack velocity scales agree. Furthermore, the potential gives the correct cohesive energy of silicon (but an inaccurate cohesive energy curve — see Sec 5.3), leading to the agreement in numerical and experimental energy scales. However, the nonlinear parts of the potential involved in stretching and rupturing bonds play an important role in determining the actual fracture energies and crack velocities, in particular where the crack arrests and what its highest velocity will be. The quantitative disagreements shown here point to a shortcoming of the nonlinear parts of the potential, which have not received much attention. The modified Stillinger-Weber potential cannot be expected to reproduce correctly experimental data. But it highlights the control over brittleness in the three-body term, and the inadequacy of the potential tails in the two-body term. These are complicated matters, and yet perhaps only hinting at greater difficulties on the road to a better potential. The lowest fracture energy density at which a crack propagated in the experiments was close to the Griffith energy density for a (111) plane, 2.2 J/m2 [84]. Since a crack cannot travel with less energy, there must be a narrow range of fracture$ energy over which the crack velocity rises rapidly from zero to the lowest value measured, 2 km/s. This phenomenon is also seen in glass and polymers [38]. Because of the extreme precision required at the boundaries, experiments with silicon are not yet capable of settling the matter of whether this sharp velocity rise signals a velocity gap. However, as shown in Fig 29, in numerical silicon, the velocity gap is temperature dependent, vanishing above 200 K. At 300 K there do exist ‘steady states’ at all velocities between 0 and 3 km/s. The silicon experiments covered a range of fracture energy densities, 2 ! 16 J  m 2 , in which the cracks produced very smooth surfaces. Thus, cracks in silicon can dissipate large amounts of energy, more than seven times the amount needed to create a clean cleavage through the whole crystal, without leaving behind any large scale damage on the fracture surfaces. Investigation by atomic force microscopy shows that for low fracture energies the fracture surfaces are flat on the nanometer scale, while at higher energies the surfaces have pronounced features. These features, however, are smooth on the micron scale, and account for height variations on the order of 30 nm over an area of 16 e m 2 . The roughness gives an area increase of only 0.1% above that of a flat cleaved surface. This extra surface  cannot account for the sevenfold increase in dissipated energy. The simulations, however, indicate that most of the energy can be carried off in lattice vibrations.

5.9 Final thoughts I begin by posing three questions:

 



How can one calculate the energy needed for a crack to propagate? When is the motion stable or unstable? How does a crack choose its path? 61

Crack velocity [km/s]

Theory and experiment in Silicon

3 2 1 0

0

2

4 8 6 2 Fracture energy [J/m ]

MEAM, 0 kelvin MEAM, 77 kelvin Hauch et al, (111) IMSW, 300 kelvin Cramer et al, (110) Deegan, (111) 10 12

Figure 33: Molecular dynamics runs at 300 K with Stillinger-Weber potentials and 0K and 77 K runs compared with experiments of Hauch et al, Cramer, Wanner, and Gumbsch [83], and Deegan.

Figure 34: Atomic force microscopy pictures of (111) fracture surfaces in silicon at three different energies: (a) G  4 3 J/m2 , (b) G  7 2 J/m2 , and (c) G  14 4 J/m2 . The crack    has traveled from left to right in a 1† 1† 2k direction. The diagonal surface markings in (b) i and (c) are approximately in a 11† 0k direction, but it is unclear what causes them. Pictures i by Rachel Mahaffy; experiments and figure composition by Jens Hauch. I will close by commenting on how well they are yet answered. For ideal brittle crystals, such as silicon, the energy needed to propagate can now be calculated. The scaling problems needed to find the effect of microscopic features on macroscopic behavior have been solved by learning from analytically solvable models. However, quantitative comparison of theory and experiment is not better than around a factor of two. Furthermore, in amorphous materials such as glass there is no theory, and in more complex materials with intermediate microstructures it is not clear what predictive value atomic–scale considerations bring. The stability of the ideal brittle crystals against microscopic branching can be determined. However, the instabilities observed in real materials are on much larger than atomic scale, and they have not been followed far or explained. 62

Finally concerning crack paths, there are phenomenological rules that appear to have much experimental support for cracks that travel slowly in isotropic materials. For cracks that move quickly, it is quite uncertain what the rules ought to be, and neither theory nor experiment have gone very far. The interest and difficulty of these problems should be sufficient to ensure that fracture continues to interest physicists for some time to come.

References [1] B. L. Averbach, Fracture, an Advanced Treatise, vol. I, ch. Some physical aspects of fracture, pp. 441–471. New York: Academic Press, 1968. [2] I. S. Grigoriev and E. Z. Meilkhov, eds., Handbook of Physical Quantities. Boca Raton: CRC Press, 1997. [3] J. W. Dally, “Dynamic photoelastic studies of fracture,” Experimental Mechanics, vol. 19, pp. 349–361, 1979. [4] N. Mott, “Brittle fracture in mild steel plates,” Engineering, vol. 165, pp. 16–18, 1947. [5] E. Dulaney and W. Brace, “Velocity behavior of a growing crack,” Journal of Applied Physics, vol. 31, pp. 2233–2266, 1960. [6] A. Stroh, “A theory of the fracture of metals,” Philosophical Magazine, vol. 6, pp. 418–465, 1957. Supplement: Advances in Physics. [7] E. H. Yoffe, “The moving Griffith crack,” Philosophical Magazine, vol. 42, pp. 739– 750, 1951. [8] T. L. Anderson, Fracture Mechanics: Fundamentals and Applications. Boca Raton: CRC Press, 1991. [9] L. I. Slepyan, Models and Phenomena in Fracture Mechanics. Berlin: Springer, 2002. [10] K. B. Broberg, Cracks and Fracture. San Diego: Academic Press, 1999. [11] M. F. Kanninen and C. Popelar, Advanced Fracture Mechanics. New York: Oxford, 1985. [12] L. B. Freund, Dynamic Fracture Mechanics. Cambridge: Cambridge University Press, 1990. [13] B. Lawn, Fracture in Brittle Solids. Cambridge: Cambridge University Press, second ed., 1993. [14] D. Broek, The Practical Use of Fracture Mechanics. Dordrecht: Kluwer Academic Press, 3rd ed., 1994.

63

[15] G. R. Irwin, “Fracture mechanics,” in Structural Mechanics (J. N. Goodier and N. J. Hoff, eds.), pp. 557–591, Elmsford, NY: Pergamon Press, 1960. [16] E. Orowan, “Energy criteria of fracture,” Weld. Res. Supp., vol. 34, p. 157, 1955. [17] B. Cotterell and J. R. Rice, “Slightly curved or kinked cracks,” International Journal of Fracture, vol. 14, p. 155, 1980. [18] G. R. Irwin, “Fracture,” in Handbuch der Physik, vol. 6, p. 551, Berlin: SpringerVerlag, 1958. [19] G. R. Irwin, “Analysis of stresses and strains near the end of a crack traversing a plate,” Journal of Applied Mechanics, vol. 24, pp. 361–364, 1957. [20] N. I. Muskhelishvili, Some Basic Problems of the Mathematical Theory of Elasticity. New York: P. Noordoff and Company, 1933. (English Translation). [21] H. Gao, Y. Huang, and F. F. Abraham, “Continuum and atomistic studies of intersonic crack propagation,” JMPS, vol. 49, pp. 2113–32, 2001. [22] A. J. Rosakis, “Intersonic shear cracks and fault ruptures,” Advances in Physics, vol. 51, pp. 1189–257, 2002. [23] R. V. Gol’dstein and R. Salganik, “Brittle fracture of solids with arbitrary cracks,” International Journal of Fracture, vol. 10, pp. 507–523, 1974. [24] J. A. Hodgdon and J. P. Sethna, “Derivation of a general three-dimensional crackpropagation law: a generalization of the principle of local symmetry.,” Physical Review B, vol. 47, pp. 4831–40, 1993. [25] A. Yuse and M. Sano, “Transition between crack patterns in quenched glass plates,” Nature, vol. 362, pp. 329–31, 1993. [26] O. Ronsin, F. Heslot, and B. Perrin, “Experimental study of quasistatic brittle crack propagation,” Physical Review Letters, vol. 75, pp. 2252–5, 1995. [27] M. Marder, “Instability of a crack in a heated strip,” Physical Review E, vol. R51, pp. 49–52, 1994. [28] H. A. Bahr, A. Gerbatsch, U. Bahr, and H. J. Weiss, “Oscillatory instability in thermal cracking: a first-order phase-transition phenomenon,” Physical Review E, vol. 52, p. 240, 1995. [29] M. Adda-Bedia and Y. Pomeau, “Crack instabilities in a heated glass strip,” Physical Review E, vol. 52, pp. 4105–4113, 1995. [30] R. C. Ball and H. Larralde, “Three-dimensional stability analysis of planar straight cracks propagating quasistatically under type i loading,” International Journal of Fracture, vol. 71, pp. 365–77, 1995. 64

[31] H. Larralde and R. C. Ball, “The shape of slowly growing cracks,” Europhysics Letters, vol. 30, pp. 87–92, 1995. [32] J. R. Rice, Fracture, vol. 2, ch. Chapter 3, pp. 191–311. Academic Press, 1968. [33] J. S. Langer and A. E. Lobkovsky, “Critical examination of cohesive-zone models in the theory of dynamic fracture,” Journal of the Mechanics and Physics of Solids, vol. 46, no. 9, pp. 1521–56, 1998. [34] A. Karma, D. Kessler, and H. Levine, “Phase-field model of Mode III dynamic fracture,” Physical Review Letters, vol. 87, p. 045501, 2001. [35] I. S. Aranson, V. A. Kalatsky, and V. M. Vinokur, “Continuum field theory of crack propagation,” Physical Review Letters, no. 85, pp. 118–121, 2000. [36] L. Slepyan, “Dynamics of a crack in a lattice,” Soviet Physics Doklady, vol. 26, pp. 538–540, 1981. [37] M. Marder and X. Liu, “Instability in lattice fracture,” Physical Review Letters, vol. 71, pp. 2417–2420, 1993. [38] M. Marder and S. Gross, “Origin of crack tip instabilities,” Journal of the Mechanics and Physics of Solids, vol. 43, pp. 1–48, 1995. [39] R. Thomson, C. Hsieh, and V. Rana, “Lattice trapping of fracture cracks,” Journal of Applied Physics, vol. 42, no. 8, pp. 3154–3160, 1971. [40] R. Thomson, “The physics of fracture,” Solid State Physics, vol. 39, pp. 1–129, 1986. [41] D. Holland and M. Marder, “Ideal brittle fracture of silicon studied with molecular dynamics,” Physical Review Letters, vol. 80, pp. 746–9, 1998. [42] S. A. Kulakhmetova, V. A. Saraikin, and L. I. Slepyan, “Plane problem of a crack in a lattice,” Mechanics of Solids, vol. 19, pp. 102–108, 1984. [43] S. I. Heizler, D. A. Kessler, and H. Levine, “Mode–I fracture in a nonlinear lattice with viscoelastic forces,” PRE, vol. 66, pp. 016126/1–10, 2002. [44] D. A. Kessler and H. Levine, “Nonlinear lattice model of viscoelastic mode III fracture,” PRE, vol. 63, pp. 016118/1–9, 2001. [45] E. Gerde and M. Marder, “Friction and fracture,” Nature, vol. 413, pp. 285–288, 2001. [46] B. Noble, Methods Based on the Wiener-Hopf Technique for the Solution of Partial Differential Equations. New York: Pergamon, 1958. [47] N. Ashcroft and N. D. Mermin, Solid State Physics. Saunders, 1976. [48] X. Liu, Dynamics of fracture propagation. PhD thesis, University of Texas at Austin, 1993. 65

[49] K. Ravi-Chandar and W. G. Knauss, “Dynamic crack-tip stress under stress wave loading-a comparison of theory and experiment,” International Journal of Fracture, vol. 20, pp. 209–222, 1982. [50] K. Ravi-Chandar and W. Knauss, “An experimental investigation into dynamic fracture: I. Crack initiation and arrest,” International Journal of Fracture, vol. 25, pp. 247–262, 1984. [51] K. Ravi-Chandar and W. G. Knauss, “An experimental investigation into dynamic fracture: II. Microstructural aspects,” International Journal of Fracture, vol. 26, pp. 65–80, 1984. [52] E. Sharon, S. P. Gross, and J. Fineberg, “Local crack branching as a mechanism for instability in dynamic fracture,” Physical Review Letters, vol. 74, pp. 5146–5154, 1995. [53] P. Manneville, Dissipative Structures and Weak Turbulence. Boston: Academic Press, 1990. [54] H. Gao, “Surface roughening and branching instabilities in dynamic fracture,” Journal of the Mechanics and Physics of Solids, vol. 41, pp. 457–486, 1993. [55] P. Gumbsch, S. J. Zhou, and B. L. Holian, “Molecular dynamics investigation of dynamic crack stability,” Physical Review B, vol. 55, pp. 3445–3455, 1997. [56] S. J. Zhou, D. M. Beazley, P. S. Lomdahl, and B. L. Holian, “Large-scale molecular dynamics simulations of three- dimensional ductile failure,” Physical Review Letters, vol. 78, pp. 479–482, 1997. [57] F. F. Abraham, “Portrait of a crack: rapid fracture mechanics using parallel molecular dynamics,” IEEE Computational Science and Engineering, vol. 4, pp. 66–77, 1997. [58] F. F. Abraham, “Simulating materials failure using parallel molecular dynamics,” in Some New Directions In Science On Computers (G. Bhanot, S. Chen, and P. Seiden, eds.), pp. 91–113, World Scientific, 1997. [59] R. K. Kalia, A. Nakano, K. Tsuruta, and P. Vashishta, “Morphology of pores and interfaces and mechanical behavior of nanocluster-assembled silicon nitride ceramic,” Physical Review Letters, vol. 78, pp. 689–692, 1997. [60] J. Q. B. F. F. Abraham, N. Bernstein and D. Hess, “Dynamic fracture of silicon: concurrent simulation of quantum electrons, classical atoms, and the continuum solid.,” Mater. REs. Soc Bulletin, vol. 25, no. 5, pp. 27–32, 2000. [61] V. B. Shenoy, R. Miller, E. B. Tadmor, and R. Phillips, “Quasicontinuum models of interfacial structure and deformation,” Physical Review Letters, vol. 80, pp. 742–5, 1998.

66

[62] H. Rafii-Tabar, L. Hua, and M. Cross, “A multi-scale atomistic-continuum modelling of crack propagation in a two-dimensional macroscopic plate,” vol. 10, pp. 2375–, 1998. [63] E. B. Tadmor, M. Ortiz, and R. Phillips, “Quasicontinuum analysis of defects in solids,” Philosophical Magazine A, vol. 73, pp. 1529–63, 1996. [64] T. Honglai and Y. Wei, “Atomistic/continuum simulation of interfacial fracture. ii. atomistic/dislocation/continuum simulation,” vol. 10, pp. 237–49, 1994. [65] T. Honglai and Y. Wei, “Atomistic/continuum simulation of interfacial fracture. i. atomistic simulation,” vol. 10, pp. 150–61, 1994. [66] A. Griffith, “The phenomena of rupture and flow in solids,” Mechanical Engineering, vol. A221, pp. 163–198, 1920. [67] M. Marder, “Small-scale simulations of fracture,” in Fracture – Instability Dynamics, Scaling, and Ductile/Brittle Behavior (R. L. B. Selinger, J. J. Mecholsky, A. E. Carlsson, and E. R. Fuller, eds.), (Pittsburgh), pp. 289–296, Materials Research Society, 1996. [68] L. Verlet, “Computer ”experiments” on classical fluids. i. thermodynamical properties of lennard-jones molecules,” Physical Review, vol. 159, p. 98, 1967. [69] D. Frenkel and B. Smit, Understanding Molecular Simulation. New York: Academic Press, 1996. [70] M. Allen and D. Tildesley, Computer Simulation of Liquids. Clarendon, Oxford, 1987. [71] H. Balamane, T. Halicioglu, and W. A. Tiller, “Comparative study of silicon empirical interatomic potentials,” Physical Review B, vol. 46, no. 4, pp. 2250–2279, 1992. [72] M. Z. Bazant, E. Kaxiras, and J. F. Justo, “Environment-dependent interatomic potential for bulk silicon,” Physical Review B, vol. 56, pp. 8542–??, 1997. [73] B. Lawn, B. Hockey, and S. Wiederhorn J. Mater. Sci., vol. 15, p. 1207, 1980. [74] F. H. Stillinger and T. A. Weber, “Computer simulation of local order in condensed phases of silicon,” Physical Review B, vol. 31, p. 5262, 1985. [75] C. P. Chen and M. H. Leipold, “Fracture toughness of silicon,” American Ceramics Society Bulletin, vol. 59, pp. 469–472, 1980. [76] M. Z. Bazant and E. Kaxiras, “Modeling of covalent bonding in solids by inversion of cohesive energy curves,” Physical Review Letters, vol. 77, no. 21, p. 4370, 1996. [77] J. Justo, M. Z. Bazant, E. Kaxiras, V. Bulatov, and S. Yip Physical Review B, vol. 58, p. 2539, 1998. 67

[78] E. K. M. Z. Bazant and J. F. Justo, “The environment-dependent interatomic potential applied to silicon disordered structures and phase transitions,” Mat. Res. Soc. Proc., vol. 491, p. 339, 1997. [79] J. Chelikowsky and J. Phillips Physical Review B, vol. 41, p. 5735, 1990. [80] R. Ravelo and M. I. Baskes, “Free energy calculations of Cu-Sn interfaces,” in Thermodynamics and Kinetics of Phase Transformations. Symposium (J. S. Im, B. Park, A. L. Greer, and G. B. Stephenson, eds.), (Pittsburgh), pp. 287–93, Materials Research Society, 1996. [81] M. I. Baskes, “Modified embedded-atom potentials for cubic materials and impurities,” PRB, vol. 46, pp. 2727–2742, 1992. [82] J. Hauch, D. Holland, M. Marder, and H. L. Swinney, “Dynamic fracture in singlecrystal silicon,” Physical Review Letters, vol. 82, pp. 3823–3826, 1999. [83] T. Cramer, A. Wanner, and P. Gumbsch, “Energy dissipation and path instabilities in dynamic fracture of silicon single crystals,” PRL, vol. 85, pp. 788–91, 2000. [84] J. C. H. Spence, Y. M. Huang, and O. Sankey, “Lattice trapping and surface reconstruction for silicon cleavage on (111). ab initio quantum molecular dynamics calculations,” Acta Metallurgica, vol. 41, pp. 2815–2824, 1993.

68