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 TheoriesofFracture: 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 Continuumfracturemechanics
2.1 2.2 2.3 2.4 2.5 2.6 2.7 2.8 2.9 2.10 2.11
2
. . . .
...... ...... ... . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .. .. .. .. .. .. ..
2 . . . . . . . . . .. .
3 3 4 6 9
Structure of fracture mechanics. . . . . . . . . . . . . . . . . . . . . . . . 9 Dissipation and the process zone . . . . . . . . . . . . . . . . . . . . . . . 9 Conventional fracture modes . . . . . . . . . . . . . . . . . . . . . . . . . 11 Stresses Around an Elliptical Hole . . . . . . . . . . . . . . . . . . . . . . 12 Stress Intensity Factor . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14 Universal singularities near the crack tip . . . . . . . . . . . . . . . . . . . 15 The relation of the stress intensity factor to energy flux . . . . . . . . . . . 16 Elasticity in Two Dimensions . . . . . . . . . . . . . . . . . . . . . . . . . 16 Steady States in Plane Stress . . . . . . . . . . . . . . . . . . . . . . . . . 19 Energy flux and limiting crack speed . . . . . . . . . . . . . . . . . . . . . 22 Crack paths . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23
3 IncorporatingPlasticity
25
4 AtomicTheoriesofFracture 26 4.1 Physical argument for velocity gap . . . . . . . . . . . . . . . . . . . . . . 28 4.2 Dynamic fracture of a lattice in anti-plane shear, Mode III . . . . . . . . . . 31 4.3 Definition and energetics of the model . . . . . . . . . . . . . . . . . . . . 31 4.4 Symmetries: . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 34 4.5 Eliminating the spatial index m: . . . . . . . . . . . . . . . . . . . . . . . 35
1
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: . . . . . . . 35 Relation between and .......................... 38 Phonon emission: . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 40 Forbidden velocities . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43 Nonlinear Instabilities . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43 The connection to the Yoffe instability . . . . . . . . . . . . . . . . . . . . 45 Generalizing to Mode I: . . . . . . . . . . . . . . . . . . . . . . . . . . . . 46
5 Scalingideasandmoleculardynamics 48 5.1 Molecular dynamics . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 49 5.2 Interatomic potentials . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 49 5.3 Realistic potentials for silicon? . . . . . . . . . . . . . . . . . . . . . . . . 50 5.4 Results of zero Kelvin calculations in silicon . . . . . . . . . . . . . . . . . 53 5.5 Along 110 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 55 5.6 Along 111 : Crackons . . . . . . . . . . . . . . . . . . . . . . . . . . . . 55 5.7 Crack behavior at nonzero temperatures in silicon . . . . . . . . . . . . . . 57 5.8 Comparisons with experiment . . . . . . . . . . . . . . . . . . . . . . . . 58 5.9 Final thoughts . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 61
1 1.1
Continuum Theories of Fracture: Scaling 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 exhaustiv ely, 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 perspecti ve 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 wit h a simple estimate of what the criti cal height shou ld be. The gravitational potential energy of the rock is Ah2 g 2 where is the densi ty. 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/m 2 . Taking the density to be is G 1 eV/ A 2000kg/m 3, the critical height at which it pays to divide the rock in two is 4G g
h
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 betwe en 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 betwee n different elemen ts and compounds, but the forces typically reach their maximum value when the distance between atoms increases by around 20% of their srcinal separation. The force needed to stretch a solid slightly is (Figure 1) F
EA
3
L
(2)
F
L d L
L
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 Theoretical Strength Practical Strength Practical/Theoretical GPa GPa GPa Iron 195-205 43-56 0.3 0.006 Copper 110-130 24-55 0.2 0.005 Titanium110 31 0.3 0.009 Silicon 110-160 45 0.7 0.01 Glass 70 37 0.4 0.01 Plexiglas3.6 3 0.05 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 breaking point is around F E (3) A 5 As shown in Table 1, this estimate is in error by orders of magnitud e. 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 greatl y 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 objec ts 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 atomicall y 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 1mm 2 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 contin uous 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 throu gh the cent er of the strip, and reli eves 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 1 2E 2 L
right
(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 G 2 L (5) per unit area crack adva nce. The reason for this assert ion 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 line ar mechan ics be applicable. Steady veloc ity 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 reflec ting from top and bottom bound aries. 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 relationsh ip 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 diame ter 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 PE CP l 2 (7) 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 tota l energy of the system containing a crack is given by CK l 2
2
l
CP l 2
qs
(8)
l
with qs
l
(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 surfa ce, which typically prevent it from happening. That the crack grow s at all is due to additional irreversible processes, sometimes chemical attack on the crack tip, sometimes vibration or other irregul ar 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 occ ur at the crack length l 0 where d
qs
0
dl
(10)
Using Eq. (9) I find that l0
2CP 7
(11)
s q
E
y g r e n e c it a t is s a u Q
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) an d Eq. (12) , with qs l0 gives t
CP CK
1
l0 l
max
1
l0 l
(13)
This equation predict s 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 dyna mics. 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 especia lly 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 Raylei gh wave speed. The reason that the dimensiona l 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 fract ure 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 reducti on of linear elasticity to two dimensions. Some basic concepts common to both static and dynamic cracks will then be discussed. These will inclu de 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 crack s. 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 secti on. The general strategy is to solve for the displac ement 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 itsel f. Thus, an equation of motion can be obtai ned 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 observ ed would appear to preclu de 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 elastic a real material support singu stresses, in this vicinity the assumption of linearly behaviorcannot must break down andlar 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
Figure 5: Structure of fracture mechanics. The crack tip is surrounded by a region in which the physics is unknow n. Outside this process zone is a region in which elast ic 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 forwar d, 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 are a of crack ext ension. The size of the proces s 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 elast icity. 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 beco mes. 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 linea r 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 require d to form a unit area of fracture surfa ce. 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 particula r 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 fract ure 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 plan e. 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
x z
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 materi al is strongly anisotrop ic. Both friction and earthquake s 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 exper imentally 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 tensi on, 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
0
u
(14)
The theory of complex variables can once again be brought to bear in order to find solutions. For the boundary proble m at hand, conformal mapping is the appropriate technique. Because u is a solution of Laplace’s equation, it can be represented by u
2 12
(15)
where is analytic, and x iy. The asymptotic behavior of is easy to deter mine. Far from the hole, the stress goes to a constant value , where is dimensionless. So u yz
y i
i 2 for
x
iy
i
x
yz
(16)
iy
(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 x t y t (18) travels around the boundary of the hole as t moves along the real axis. Then T
x
y
t
t
and
N
y
x
t
t
It is easy to verify that N T =0, and because T is tangent N must be normal.
(19)
are tangent vector and normal vectors along the edge of the hole, so requiring normal stress to vanish means that xz
N
yz
0
(20)
u
u
y
x
x
y
t
t
0
u x
u y
y t
x t
x ix t
ix
y
t
iy
iy
0
(21) (22)
Insert (15) for u.
t
This equality holds only on the boundary.
(23)
Because is arbitrary up to a constant anyway, there is no worry about dropping a constant of integration.
(24)
t
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
with
(25)
p is a number lying between 0 and 1.
lying on the unit circle, ei
and real. When p 0, the boundary is circular, and when p along the real axis. Considering as a function of , one has 1
(26) 1, the boundary is a cut
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 on the unit circle. Equation (27) is a relation betwee n two analytic functions that holds over the whole unit circle. The difference between and 1 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 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 for large . The relation between and is 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.
4p
2 Therefore, when
is large, large,
(28)
, and in accord with Eq. (17) one obtains for
i
(29)
Consulting Eq. (27) one concludes also that 1
for
i
(30)
which means that i i
for
0
for
0
(31) (32)
However, can have no other singularities within the unit circle, or else 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 i i
i
2
1
(33)
Add Eqs. (32) and (29).
1
Use Eq. (28) and 1
4p
2 1
i 1
2p 4p
1 2
1
4p
2
(34)
2 p.
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 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
Radius of curvature R
l
yz
y
x z
0
20
y
15
g n o l a ,
10
z y
s s e r t
5
S
0
3 2 Distance x
1
4
Figure 7: The stresses at the tips of an elliptical hole in a solid are much greater than those applied off at infinity. The plot shows a solution of yz derived from Eq. (35) for 1 and m 0 9. Note the 20-fold incre ase 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 u yz
x
y
x
2 x
2
x
2
as x
2
yz
becomes
(35)
The stress intensity factor K is a coefficient of the diverging stress, defined by K
In the case of Eq. (35), K
2.6
lim
r
0
2 r 2
yz
r
2 is the distance from the crack tip.
x
(36)
.
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 comp utations were rest ricted to the case of Mode III, antiplane shear, for which the algebr a 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 surroun ding tip develops square ity.in As first noted by Irwin [18],field the stress field at athe point r neara the crackroot tip,singular measured polar coordinates with the crack line corresponding to 0, takes the form i j
K
1 2 r 15
fi
j
(37)
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 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 quan tity 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 K fi j (38) i j 2 r 1
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 2 1 (39) G A K2 E 1 where is the Poisson ratio of the material and the three functions A 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 compl ete forma l development of fracture mechanics is very elaborate. For stati c 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 univer sal 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, 2 u 2 (40) u u t2 srcinally found by Navier. u is a field describing the displacement of each mass point from its srcinal location in an unstrained body, is the density, and the constants and are called the Lam´e constants, have dimensions of energy per volume, and are typically of order 10 10 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 redu ces 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 uz
xz
(41)
x
and
uz
yz
(42)
y
The equation of motion for uz is 1 c2
2 uz
t2
2
(43)
uz
where (44)
c
Therefore, the vertical displacement obeys the ordinary wave equation. A second case corre sponds 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
(45)
z f ux uy
and that ux and uy are independent of z. One can deduce the functio n f by noticing that must vanish on the face of the plate. This means that ux
uy
x
y
2 17
uz z
0
(46)
zz
at the surface of the plate, which implies that uz
f u x uy
2
z
So ux
uy
uz
x
y
z
2 2
ux
uy
x
y
ux
uy
x
y
(47)
(48)
and one can write u
u
u
x
x
x
where
(49)
2
(50)
2
and and now range only over x and y. Therefore, a thin plate obeys the equat ions 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, E
1
E
1
2
2
1
E 2 1
(51)
The following relation will be useful in discussing two-dimensional static problems . First note that u 2 (52) Second, taking the divergence of Eq. (40), one finds that 2 2
t2
2 Therfore,
u obeys the wave equation, with the longitudinal wave speed
2
cl
Similarly,
(53)
(54)
u also obeys the wave equation, but with the shear wave speed
(55)
c t
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 rotati on of the plane of polarizat ion must be of the form 1
K
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 stres ses 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 u 2 2 u u (57) x2 Divide u into transverse and longitudinal parts so that u
ut
19
ul
(58)
with ul
and ut
l
t
t
y
x
(59)
It follows immediately that 2
2
2
2
ul
x2
2
f
x2
for some function f which must be harmonic ( f x that 2
2
2
x2
y2
2
2
x2
y2
2
2
2
i f y is a function of x
(60)
ut
iy). We have then
2 l
0
(61)
2
0
(62)
t
where 2
2
1
(63a)
2 2
2
1
(63b)
Therefore the general form of the potentials is l
0 l
t
0 t
0 l
z
0 t
z
z
1 l
x
i y
1 l
x
i y
(64)
z
1 t
x
i y
1 t
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 simul1 1 taneously, and can be neglected. Defining z z and z z I have t z l z for u ux uy
z
z
i
i
z
z
(66a)
z
z
z
(66b)
i y
(67)
z
where z
x
i y
z
x
Equation Eq. (66) gives a general form for elastic problems which are in steady state moving at velocity . Define also z z and z z Then the stresses are given by xx xx
2
2
yy
2
yy
2
1
xy
2
2i
z 2
z z
1
z
(68a)
z 2
z
20
2i
z
z
(68b)
1
z
z
(68c)
we will also need the rotation, which means uy
ux
x
y
2
1
z
(69)
z
It is worth writing down the stresses directly as well:
xx
2
1
yy
1
2
2
z 2
2i
z z
z
2i
z
(70)
z z
(71)
z
The definitions of and in Eq. (63) have been used to simplify the expressi ons. To solve a general problem, one has to find the functions and which match boundary conditions. It is interesting to notice that when 0, the right hand side of Eq. (68a) goes to zero as well. Since one will be findin g the potentials from giv en stresses at the boundaries, must diverge as 1 , 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 assumpt ion is verified in all cases that can be worked explicitly.) Near the crack tip, I assume that z
Br
iBi z1
2
(72)
z
Dr
iDi z1
2
(73)
We first appeal to symmetry. Observe that ux
y
ux y
uy
y
(74)
u y
Placing Eq. (73) into Eq. (66) and using Eq. (74), we find immediately that B i Dr 0. Thus Br iDi z z (75) z1 2 z1 2 We also observe that the square roots in Eq. (73) must be interpreted as having their cuts along the negative x axis, corresponding to the crack. On the crack surface, we have two boundary conditions, which require that xy and yy vanish. Upon substituting Eq. (75) yy is satisfied identically for x into Eq. (71) we find that the condition upon However, substituting into Eq. (69) with y 0 we find that
xy
i 2 Br
2
21
1 Di
1 x
1 x
0 y
0.
(76)
Thus
Di
2
Br
2
(77)
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 K xx
2
1
2 D
1
2
2
1 z
2
1 z
1 z
4
1 z
(78a)
2
1
K 2 D
yy
K xy
2 D
2
2
2i
1
1 z
1 z
4
1 z
1 z
1 z
1 z
1 z
1 z
(78b) (78c)
with D
4
2 2
1
(78d)
The constant K is again called the stress intensity factor, and is given by lim x
K
0
2 x
(79)
yy
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 r cos
z
ri sin
r cos
z
ri sin
(80)
xy
(81)
and to evaluate the stress cos2
yy
sin2
xx
sin 2
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 dt
K
d
P
dt
dxdy
2
u u
1 u 2 x
(82)
The spatial integral is taken over an region which is static in the laboratory frame. So d dt
K
P
dxdy
22
u u
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 u
x
(84)
x
x
we have dxdy
u
u
x
dxdy x S
u
(85)
x
(86)
u
(87)
n
where the integral is now over the boundary of the system, and n is an outward unit normal. By using the asymptotic forms Eq. (78) for yy and the corresponding expressi ons 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 4
1 1
2 2
K2
(88)
2 2 The factor 4 1 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 geome-
try. 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 discu ss 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 conse quence 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 consequences for crack motion. Consider an initially straight crack, propagating along the x axis. The stress field component s xx and yy have the following form xx
yy
KI 2 r 1 KI 2 r 1
2
2
r1
T r1
2
(89a)
2
(89b)
The constant stress T is parallel to the crack at its tip. Cotterell and Rice show ed 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 stres ses induc ed by the nonun iform thermal expansion of the material, follows that of the glass plat e. At a critical pulli ng 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 princip le 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 calculation, first described by Rice [32]. Let us look at the ratio xx yy right on the crack line. From Eq. (78) it is xx yy
2
1 1 4 24
2
2 1
2 2 2
4
(90)
2
2 4
2
1
2
1
2 2
1
The Taylor expansion of Eq. (91) for low velocities 2
1
ct4
2cl2 ct2 cl
(91) is
cl4 ct cl
ct
(92)
and in fact Eq. (91) greater than unity for all . This result is surprising because it states 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 ana lyses. They do find some models in which a crack undergoes a Hopf bifurcation to an oscillation at a critical velocity. However, their “general conclusio n is that these cohesive- zone models are inherently unsatisfactory for use in dynamical studi es. They are extremel y 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 Levi ne [34] and Aranson, Kalatsky V inokur [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 resolutio n 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 interac tions 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): (93)
c
where
c
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 appr oach 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, Lennar d-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 proces s 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 introduci ng 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 res trictive? 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 gover ned by the structure of the stress fields in the near vicinity of the crack tip. These fields are solely deter mined 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 dynamic s 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 simul ations indicate that lattice trapping disappea rs 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 calcul ations seem to be remarkably robu st. 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 tempera ture. 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 qualitati ve 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 velo city 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 betw een two ato ms has just brok en. 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 M
2
(95)
b
where M is its mass. Assuming that the nex t bond, at distan ce a, breaks during some fraction 1 of half this oscillato ry 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 a2
Mu
2
u
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 wh en one of the springs connected to an atom is snapped, the effective spring constant describing its vibrations should decrease, so in Eq. (96). b . There is also an extr a factor of 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 comple tely 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 snap s backwards each time a bond conn ected 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 tempe rature, thermal fluctuations can bring extr a 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. d2
d1
d5
d 1 d2 d5 d6 d3 d4
d6
d4 d3 u2
u1 u0
u5 u3
u6
u4
Figure 11: The geometry of the lattice used for fracture calcula tions.
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 unstabl e. 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 preferabl e 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
xm
1
N
1
1
1
N xm
N 1
1
xm
1
1
1
1
xm
1
xm 1
1
xm
1
N
1
1 1
N
N
2N
1
Figure 12: In this simplest of all solv able atomic crack model s, 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 stron g springs break if they exceed a certain extension.
32
linear springs until the instant they snap at a separation of 2 u 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. The model is described by the equation
Figure 13: Dynamic fracture of a triangula r crystal in anti-pl ane shear. The crystal has 2 N 1 rows of atoms, with N 4, and 2 N 1 rows of bonds. Heights of spheres indicate displacements u m n of mass points out of the page once displace ments are imposed at the boundaries. The top line of spheres is displaced out of the page by amount u m N 1 2 2N 1, and the bottom line into the page by amount u m N 1 2 UN UN 2N 1. Lines connecting mass points indicate wheth er the displacement between them has exceeded the critical value of 2u f ; see Eq. (99b)).
u m n
bu
1 2
f u m n nearest neighbors m
u m n
(99a)
n
with f u
u 2u f
(99b)
u
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 2
(99c)
UN
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
(100)
1 2
and the height difference between mass points with adjacent values of n is Uright
UN N
(101)
1 2
Therefore, the energy stored per unit length in the 2 N 1 2 Upper Bonds Site 2 1 1 UN 2 2 2N 1 2 2 N 1 2 2Q0 UN
Rows
1 rows of bonds is
Spring Constant
2 Uright
(102) (103)
2
(104)
with Q0
1 2N 1
(105)
The energy required to snap two bonds each time the crack advances by a unit length is 1 2 1 2
2 Bonds Site 2
1 2
2u f
2
Spring Constant 2u2f
Separation at fracture
2
(106) (107)
Therefore, equating Eqs. (104) and (107) the proper dimensionless measure of the external driving is UN Q0 (108) uf 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 relate d 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 beca use the dista nce 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 srcinal 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 t . The closest one can come in a triangular crystal is by having the symmetry u m n t
1 nt
u m
1
(109a)
and also u m n t which implies in particular that u m 1 2 t
We have defined
4.5
n t
u m
1 2t
0 1 mod n
gn
1 2
u m
1 2
(109c)
if n 1 2 5 2 if n 3 2 7 2 in general
1 2 2
Eliminating the spatial index
(109b)
gn
(110)
:
m
Assuming that a crack is in steady state, we can therefore eliminate the variable m entirely from the equation of motion, by defining u 0 n t
un t
(111)
and write the equations of motion in steady state as un t
if n u1
1 2
un 1 t gn un t 1 un 1 t gn
1
1
un 1 t gn un t 1 un 1 t gn
6 un t 1
1
1
(112a)
bun 1
1 2, and
2
t
1 2
u3 u1 u
2 2
t t
1 2
1 t
4 u1 u1 2 t
2
t t
u3 u1 u
2 2
1 1
t t
1 2
t
1
u1
2
t
1 2
t
(112b) bu1
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 chosen 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
2
n
1 2. Fourier transforming Eq. (112a) in time gives 2
un
1 u 2 n 1 1 2 un 1 2 un 1
ib
ei ei ei
gn
1
1
6 gn
1
1
ei
gn
1
ei
gn
1
i
e
(113)
Let un
u1
ek
2
n 1 2
i gn
2
Substituting Eq. (114) into Eq. (113), and noticing that g n 2
u1
ib u1
2
1 2 u1 2 1 u 2 1 2 1 2 u1 2
2
2
k
i
e
2cosh k cos
ib
ei e ei
k
e
gn
2
gn
1
2
2
gn
1
6 gn
(114)
i
ei
2
2
gn
1
gn
2
gn
1
gn
2
e
i
e gn
1
1 gives
cos
3
(115) (116)
0
Defining 3
z
2
cos 2cos
ib
(117)
2
one has equivalently that y
z2
z
1
(118)
with ek
y
(119)
One can construct a solution which meets all the boundary conditions by writing N 1 2 n
un
u1
2
i gn 2
e
y
N 1 2 n
y y
yN
UN n
N
N
1 2
2
2
2
(120)
2 2 This solution equals u1 2 for n 1 2, and equals UN 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 this boundary cond ition is a delta function, and hard to work with formal ly. To resolve uncertainties, it is better to use instead the boundary condition
uN
1 2
t
t
UN e
(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 u1
U t
2
t
u
1 2
t
u1
2
t
2
u1
2
t
1 2
(122)
2
Rewrite Eq. (112b) as u1
2
t
1 2
u3 2 t u1 2 t 2U t
1
4 t
u1 2 t 2U t 1 2
36
u3 2 t 1 u1 2 t 1 1 2 t
bu1
2
(123)
Fourier transforming this expression using Eq. (120) and defining U
d ei t U t
t
2
UN
(124)
now gives u1
F
2
with
yN
F
1
ei
1
y
yN
y
U
2 2
N
2
(125)
N 1
2z cos
N
2
1
(126)
Next, use Eq. (122) in the form 1
U
i
e
2
u1
2
to obtain U
2 cos2
F
(127)
2
UN
4 U
N
2 2
2
(128)
Writing U
U
(129)
U
finally gives U
Q
U
1
UN Q0
with
1 i
i
F
Q
2cos 2
F
(131)
4
To obtain the right hand side of Eq. (130) one uses the facts that F 0 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
(130)
1 N , and that
(132)
Q
where Q is free of poles and zeroes in the lower complex plane and Q is free of poles and zeroe s in the upper complex plane. One can carry out this decom position with the explicit formula d
exp lim
Q
ln Q
2
0
i
i
(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
Q0UN
Q
Q
1
Q0UN Q
0
i
0
1 i U Q
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
UN
U
UN
Q0 Q
0
Q
(135a)
i
and Q0 Q
0
Q
(135b)
i
One now has an explicit solution for U . 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 afterwards. A solution of Eq. (135) constructe d in this manner appear s in Figure 15.
4.7
Relation between
and
Recall that making the transition from the nonlinear problem srcinally 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 f at t
u 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
d
lim
e
2
0
i t
U
This integral can be evaluated by inspection. One knows that for positive t d exp
(137)
uf
0
i tU
0, (138)
and that any function whose behavior for large is 1 i has a step function discontinu ity at the srcin. Therefore, Eq. (137) and Eq. (135a) become uf
Since from Eq. (131) it follows that Q
UN Q0
Q Q
0
1, one sees from Eq. (133) that 38
(139)
Q
1
Q
(140)
As a result, one has from Eq. (139) and the definition of Q
given in Eq. (108) that
0
(141)
Q0
To make this result more explicit, use Eq. (133) and the fact that write 1 ln Q i 2
d Q
0
exp
2 d
exp Q
0
1 2i
2
2
to
Q
ln Q (142)
i
ln
d
Q0 exp
Q
Q 2
Q
2
ln Q 0
Q 1 ln 2i Q
(143)
1 Q ln 2i Q
(144)
Placing Eq. (143) into Eq. (141) gives d
exp
2
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 d
C exp
2
1 2i
ln Q
ln Q
(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 r i 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 r i pi
C
for b
r i pi
0
(146)
One may derive Eq. (146) as follows: away from a root or pole of Q, the integrand of Eq. (145) vanish es. Consider the neighborho od of a root r of Q which falls to the real axis from the negative side as b 0. For the sake of argument, take the imaginary part of b, the integral this root to be ib. In the neighborhood of this root, say within a distance to compute for Eq. (145) is 1 2
r r
b b
d 2i
ln
r
ln 39
ib r
ib
(147)
Defining
r , and integrating by parts gives
1 2 1 ln r 2
b b
d
2i
ln
r
2ib 2
(148)
b2
(149)
b
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, appear in Figure 14.
, to the velocity of the crack
. The results of a calculation
0.8 Steady states go unstable 0.6 c
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 3 2, is plotted as c, scaled by the sound speed c 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 diss ipation 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
3.0 ) t ( u
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.
3.0 5 2.5
k
2
2.0 1.5
k
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 (150)
k
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 mul
Rl
D
Rl u l
Rl
l
t
(151)
l
l
Multiplying everywhere by eik R , summing over l , letting K be reciprocal lattice vectors, and letting be the volume of a unit cell gives mu k
k u k
D
1
ei
k
K t
k
K
(152)
K
Inspection shows that the lattice frequencies excited in this way are those which in the extended zone scheme[ [47]] obey k
(153)
k
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 argumen t that is both simpler and more general. The only way to transport radi ation 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 R i t gives e ik
R na
i
k t na
e ik R
i
k t
(154)
Assuming a and are parallel, Eq. (153) results again. There are two phonon dispersion relat ions to consider. One gives the conditi ons 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 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 2cos 2 4 0. According to Eq. (131), the roots and poles of Q 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 4.0 U
t n e m e c la sip D
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 2 u f at t 0, but this must be the first time at which that bond stretches to a length greater than 2 u f . For 0 03 (the precise value of the upper limit varies with b and N ) that condition is violated. The states u fintact, have 0, thethe unphysical character shown in lower Figureline 17.ofMasses above height for t less than bond connecting them to the massesrise remaining however and then they descend, whereupon the bond snaps. Since the solution of Eq. (112) is unique, but does not in this case sol ve Eq. (99a), no soluti ons of Eq. (99a) exist at all at th ese 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 repea ting. No soluti ons 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 whi ch 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 scale d 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 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 directl y. 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 connecti on 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. 1 147 c
0 645
1 165 c
0 630
1 376 c
0 624
1 835 c
0 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 struct ure. The simulations are carried out in a strip with half-width N 9, of length 200 and b 0 01. The fron t and back ends of the stri p 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 advance d. The diagram show s patterns of broken bonds left behind the crack tip. 44
Just above the
Figure 19: Images of local crack branches in the x y plane in PMMA. The arrow, of length 250 m, indicates the direction of propagation. All figures are to scale with the path of the main crack in white. (top) 1 18 c (bottom) 1 45 c . (from [52]) c (center) 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 framewo rk 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 instabi lity 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 veloc ity seen experimen tally 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 Rayle igh wave speed in the vici nity of
45
a crack tip may be significantly lower than the value of c R 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 calculati ons are at zero. 3. The experiments are in amorphous materi als, while the calcul ations are in crystals.
4.12
Generalizing to Mode I:
Let ui describe the displaceme nt of a mass point from its equilib rium 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 displace ments. Translational invariance demands that the forces between particles 0 and 1 depend only upon u 1 u0 , which will be defined to be 1 . However, the force can be a general linear funct ional 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 d1 , and a component that is along d1 . 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 appreciate d by Born. Their quantum-mechanical srcin 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 ui 1 j 1 ui j is 1 d
1
1
d
d
1
1
1
d
1
(155)
1 3 1 x 3 y 3 1 3 x 1 y (156) 2 2 2 1 2 1 2 2 2 1 2 1 Adding up contrib utions from other particles in this way we get for the force due to neighbors 6
F m n
q dq j
j
m n
dq j
(157)
j 1q
By varying the constants and , one can obtain any desired values of shear and longitudinal wave speeds, which are given by c2l ct2
3 a2 8m 3 a2 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 techn ical restriction in the calculations of Ref. [38]; right on the crack line, forces are required to be central. 46
0.8
Microcracking instability
0.6
Minimum velocity steady state
Height=80 Height=160
c t
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 para meter 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 preci sely in Eq. (108). The sprin g 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 dissipat ion. In the limit of central forces, 0, it turns out to be difficult but not impossib le 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 amoun t of dissipation. Crack motion is greatly stabilized by having nonzero . Figure 20 displays a represen tative result with some impli cations 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 Ilattices will not refer to much. Kessler anddissipatio Levine have condu the most careful studies of with nonlinear force laws and n acting incted concert [43, 44]. Their findings emphasize the point that crack dynamics are terribly sensitive to features of the tip that do appear in any conti nuum 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 studi es 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 enti rety in the midst of one computer program, I prefer to use computer calculations to tabulate quantities that linear elast ic 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 (159) G Gc where G is the fracture energ y density alread y 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 togeth er 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 dynami cs 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 scalabl e to the macroscopic level. A good scaling argumen t 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 reaso ns for doing comput er simul ations is to relat e the two. The simul ations I will describe here in detail were carried out in silicon. Silicon is extremely brit tle, and highquality macroscopic single crystal wafers are cheap. It is therefore an excellent candidate for laboratory fractur e 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 ther e are over thirty effective (although rather ineffective) interatomic potentials for silicon in the literature [71, 72]. Furthermore, transmissio n electron microgr aphs of cracks in silicon wafers reveal atomically sharp crack tips [13,73]. Silicon therefore is an obvious candidate for molecular dynamics investi gations 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 equat ion under lying 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 srcin, the ab inito methods are the most reliable techniq ue 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 achievable. 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 intera tomic 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 structu re 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 literatu re. It gives excellent elastic propert ies, and captures well the nonlinear physics involved in heating and melting. It therefore is a reasonable starting point for conducting molecular dynamic s fracture simulations in silicon. Unfortunately,
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 110 in the plane (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 fractu re, i.e. where the newly creat ed fracture surfaces are atomical ly 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 tetrahed ra in the diamond lattice, and thus making the crystal more brittle. This can be done by scaling a parameter, , the coupling consta nt in the thr ee-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 21 (dimensionless). However, by doubling this, as in the pioneering inadvertent modification, one can obtain fast brittl e fracture. With a fast crack running, if one quasistatically decreases , the crack arrests well before one reaches 21, Fi g 23. Although 42 gives fracture phenomenology in reasonable accord
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
Figure 23: Crack velocity profile (a) along (110) 110 , and (b) along (211) 111 in modified 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 srcinal 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 expandin g 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 calculat ion 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 Baske s [80, 81]. This potential is based upon a number of physical id eas 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 predict ed 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 srcinally at crystal location R is related to the vertical displacement zR na of an atom n lattice spacings a ax to the right by zR
na
t
na
zR t
(161)
Silicon Original SW IMSW E100 GPa 130 114 172 E110 GPa 169 139 189 151 201 E111 GPa 188 Table 2: Elastic constants of silicon, comparing Stillinger-Weber, modified StillingerWeber potentials, and experiment.
53
6 5 ) ˚A / V (e e c r o F
Bazant-Kaxiras
4 Stillinger-Weber
3
Density functional
2 1 0 1.0
1.4
1.2 a
1.6
a0
(A) 10 8
) V e ( m6 o t a r e p 4 y g r e n E2
0 1
(B)
Density Functional Theory Modified Stillinger Weber Original Stillinger Weber Modified Embedded Atom Method
1.5
2 2.5 Extension Ratio
3
Figure for 24: the (A)restoring Comparison ofduring densityuniform functional theory, EDIP, andisthe Stillinger-Weber potential force expansion. MEAM 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 obtai n the perfe ct period icity 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 cruci ally important when proceeding throug h 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 a tenth of acomputer, microsecond To achieve this, motivated one needs not only computational efficient code and a highofperformance but. also a physically smallest 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
4
s) m to rs g3 n A ( t n e m e c a l p is d2 l a c tri e v
1
0
12345 time (ps)
Figure 25: Crack that has been trav eling for over 0.24 ns. That the two ov erlapped curves are almost completely indistinguishable shows the crack has reached steady state, according to Eq. 161, and is emitting phonon s in accord with Eq. 153.
5.5
Along 110
The relation between velocity and load for cracks exposing 110 and traveling along 110 is shown in Fig. 26. The crack velocity smoothly decreases as decreases, until at 2256 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 overl aps 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 is shown in Fi g. 27. For 1 44
for cracks along 111 and traveling along 2 2, the crack has stable steady states, and 55
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.
Figure 27: Relation between crack speed and loading for cracking along 111 011 at 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 hystere sis 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 travel s 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 withi n a few m/s. On cyclical loading the same plate aus 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 behavio r. All these phenome na are easily disguised if strain rates are too high. Resolving all the fine struc ture visible in Fig. 27 requir ed 8 s 1 , or hz c 10 5 . 4
unstable
s)/ 3 m (k d e e p s 2 k c a r c
1
0
1.2
1.4
1.6
1.8
2.0
2.2
Figure 28: Crack velocity profile for Si (111) 211 at 0 K. Dashed lines indicate forbidden 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 strengt h of lattice trapping shoul d be a function of temperat ure. A series of simulations at 50 Kelvin intervals were carried out to investigate the effect of temperature on lattice trapping results are shown in crack Fig 29.actually heals and travels Note also inalong Fig 30 111 that as011 .isThe decreased below 1 the backwards. This is perfectly reasonable: there hasn’t been any nonuniform surface damage, and no oxide layer has formed.
57
Figure 29: Velocity gap vanishes near 200 K. Dotted lines indicate forbidden velociti es. All low-lying velocity states exist at room temperature. Cracking is along (111) 011 .
5.8
Comparisons with experiment
I turn to a number of questio ns that allow a comparison of theory and experime nt. First, is there a velocity gap in brittle amorpho us materials at room temperature? The answer seems pretty clearly to be “no;” data to this eff ect 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 silico n was the possibility of carrying out experiment s 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 112 in order to allow direct comparison 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 m during the course of the simulations as G was varied between 5 and 14 J/m 2 . The room-temperature experiments for fracture in silicon wafers along (111) 112 are described in [82]. These experiments are difficult to perform because of the high Young’s modulus and brittleness of silicon. 58
Figure 30: Crack velocity profile with respect to quasistatic loading silicon at 300 K. All low-lying velocity states exist.
along 111 011 in
300 PMMA 200 100 0 0 ) c 300 e s/ m 200 ( y it c o l e V
100
200
300
Homalite-100
100 0 0
100
200
300
400
500
75 Glass
50 25 0 0
400
800 Time t (
1200
1600
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
6
] m4 c [ th g n2 e L
3 s)/ m2 (k
0 0
2
4
6
8
10
12
14
16
18
4
] /s 3 m [k 2 y ti c o l 1 e V
Experiment Molecular Dynamics
1 0 0
0 0
2
4
6
8
10
12
14
16
18
5
10
15
20
Energy flow to crack tip G [J/m2]
Time[µs]
(A)
(B)
5
4
] /s m k [ ty i c lo e V
3
2
Experiment (111, Hauch) Experiment (110, Cramer) Experiment (111, Deegan, ramp)
1
0 0
20
40
60
80
G [J/m2]
(C) Figure 32: Experimental and numerical determination of the crack velocity as a function of fracture energy G along 111 112 at room temperature . The fuzz indicates the spread 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 giv es 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 receive d much attention. The modified Stillinger-Weber potential cannot be expected to reproduce correctly experimental data. But it highlights the control over brittlene ss in the three-body term, and the inadequacy of the potential tails in the two-body term. These are complica ted 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/m 2 [84]. Since a crack cannot travel with less energy, there must be a narrow range of fracture ene rgy 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 featur es. These features, however, are smooth on the micron scale, and account for height variations on the order of 30nm over an area of 16 m 2 . The roughness gives an area increase of only 0.1% above that of a flat cleaved surfa ce. 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 questio ns: 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
Theory and experiment in Silicon
s]/ 3 m [k ty i c2 lo e v k c a r1 C
0
MEAM, 0 kelvin
0
2
4
MEAM, 77 kelvin Hauch et al, (111) IMSW, 300 kelvin Cramer et al, (110) Deegan, (111) 10 12
8
6
2
Fracture energy [J/m ]
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/m 2 , (b) G 7 2 J/m 2 , and (c) G 14 4 J/m 2 . The crack has traveled from left to right in a 112 direction. The diagonal surface markings in (b) and (c) are approximately in a 110 direction, but it is unclear what causes them. Pictures 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. Meilkh ov, eds., Handbook of Physical Quantities . Boca Raton: CRC Press, 1997. [3] J. W. Dally, “Dynamic photoelasti c 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 fract ure 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. Popela r, 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, “Fractur e 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 curv ed or kinked cracks,” International Journal of Fracture, vol. 14, p. 155, 1980. [18] G. R. Irwin, “Fract ure,” in Handbuch der Physik , vol. 6, p. 551, Berli n: SpringerVerlag, 1958. [19] G. R. Irwin, “Analysis of stresses and strai ns 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 arbitrar y 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, “Experime ntal study of quasistatic brit tle 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 analy sis 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 instabili ties,” 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 fractur e,” Solid State Physics, vol. 39, pp. 1–129, 1986. [41] D. Holland and M. Marder, “Ideal brittle fracture of silicon studied with molecul ar dynamics,” Physical Review Letters, vol. 80, pp. 746–9, 1998. [42] S. A. Kulakhmetova, V. A. Saraikin, and L. I. Slepyan, “Plan e 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, “Nonl inear lattice model of viscoelas tic 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 dynami cs 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 molecula r 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. Tsurut a, 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. Bernst ein and D. Hess, “Dynamic frac ture 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 analy sis 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. Wiederho rn 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 cova lent 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 transition s,” 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 calculati ons 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-a tom 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 fractu re 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
