- Email: [email protected]
A simple model for elastic fracture in thin films
Thin Solid Films, 151 (1987) 165-190 METALLURGICAL
A SIMPLE
AND
165
PROTECTIVE COATINGS
MODEL
FOR
ELASTIC
FRACTURE
IN THIN
FILMS
PAUL MEAKIN
Central Research and Development Department *, E. I. du Ponl de Nemours and Company, Wilmington, DE 19898 (U.S.A.) (Received
September
5, 1986; accepted
March 25, 1987)
A simple model for the cracking of thin films has been developed and investigated. The film is represented by a network of bonds and nodes which initially form a triangular lattice in which each node is at the junction of six bonds. The nodes are also attached to a rigid substrate by bonds which are not broken but have a relatively small force constant. Cracking is simulated by a sequence of thermally activated bond breaking events followed by mechanical relaxation. If the stretching force constant associated with the bonds in the surface layer is large, linear cracks propagate rapidly and, at a later stage, are connected by secondary cracks which become less and less linear as the strain in the surface layer is relaxed by the cracking process. Many of our simulations exhibit a distinct “initiation” period in which the bond breaking rate is low, followed by a period of rapid crack propagation in which a large fraction of the surface strain energy is released. This period is then followed by a second period of relatively slow bond breaking. During the first period isolated defects are formed, followed by linear cracks in the second period. During the third period non-linear cracks connect the linear cracks, and structures which resemble shear bands are formed. Results for the effects of finite relaxation rates for the elastic network (visco-elastic effects) are also presented.
1. INTRODUCTION
Thin films are now finding a broad variety of important applications including erosion and corrosion protection under extreme conditions, optical and electronic devices, wear reduction and a wide range of aesthetic applications. In many of these applications cracking of the surface layer is a serious form of failure which limits reliability and overall performance. One of the most familiar examples of this type of failure is the cracking and “checking” of paint films’-‘. This form of failure is now relatively uncommon (but often catastrophic when it does occur) in most organic coating systems. In other cases, cracking is a serious problem which inhibits the adoption of otherwise attractive new technologies. Cracking of the surface layer can * Contribution
number
004~6090/87/$3.50
4176.
0 Elsevier Sequoia/Printed
in The Netherlands
166
P. MEAKIN
lead to loss of required optical, electronic or mechanical properties in the surface layer itself or failure of the substrate as a result of corrosion, stress concentration etc. A key ingredient in the fracture of thin films which have been deposited onto more or less rigid substrates is the development of stress in the deposited film. This may be a result of solvent 10ss’-~, other changes in chemical compositioniU7, changes in temperature’, phase changes in the surface layer, deformation of the substrate etc. Another characteristic of thin film systems which makes them prone to failure by cracking is the (often very large) mismatch in properties between the surface layer and the substrate. Here a simple model is investigated which is intended to represent the case of a “hard” surface layer on a soft substrate or a hard surface layer bonded to a hard substrate by means of a layer of much softer material. Such systems include ceramic-coated metals, metal-coated polymers and paint systems in which the surface has become hardened as a result of oxidation and other (physical or chemical) changes. The growth of cracks in an elastic medium is similar to other non-equilibrium processes such as diffusion-limited aggregation (DLA). In fact, part of the motivation for this work was the similarity between crack growth and the DLA model of Witten and Sander’. This similarity becomes more apparent in the dielectric breakdown” and diffusion-limited growthl’ versions of the DLA model. In these models, the growth probability depends on the local value of a field 4 (electric potential in the dielectric breakdown model and concentration in the diffusion-limited growth model) near to the surface of the growing object. The field is obtained by solving Laplace’s equation (V’C$ = 0) subject to a fixed value of zero on the growing object and a constant value (say unity) on some boundary which encloses the growing object. In these stochastic growth models the growth probability is assumed to be some function of the local value of the field C$or its gradient. The DLA model and related non-equilibrium growth models have been studied extensively in recent years 9- I4 . Despite the apparent simplicity of these models and some new ideas”-’ 7, they are not yet well understood theoretically. The major difficulties seem to result from the fact that, whereas the growth probability at some point on the surface depends on the local value of the potential (or equivalently the potential gradient since 4 = 0 on the growing structure), the potential at each point depends on the entire structure which typically has a complex (often fractal18) geometry. The process of elastic fracture’ 9 can be thought of in similar terms. Here, we are concerned with a stress field 0 which cannot propagate across crack surfaces. As in the non-equilibrium growth and aggregation problem the determination of the local value of CTis a non-local problem made more difficult by the tensor nature of the stress field. As in the DLA and related models, we assume that the rate of crack propagation (and initiation) depends on the local stress field. Because of the theoretical difficulties associated with crack growth in systems with complex geometries, it seems probable that simple computer models (such as that described below) will prove to be valuable in guiding and stimulating theoretical work. In our model (described in detail below) the elastic medium is represented by a network of bonds and nodes. The bonds are broken randomly with probabilities which depend on the stress in each bond. After each bond breaking event, the system is relaxed to a new mechanical equilibrium and the new bond breaking probabilities
SIMPLE
MODEL
FOR ELASTIC
FRACTURE
IN THIN FILMS
167
are calculated. The model is kinetic in nature and assumes that bond breaking is a thermally activated process. In this respect, it is related to the early models of Zhurkov and coworkers20P22 and of Bueche and Halpin23P26. The model described below is also closely related to other models27P29 which have been used to explore the possible origins of fractal geometry” in fractal surfaces29-32 and failure in polymer fibers.33 2.
COMPUTER
MODELS
In all our simulations the surface layer is represented by an array of nodes and bonds. At the start of each simulation these nodes and bonds form a triangular lattice (with periodic boundary conditions) in which each node is connected to six nearest neighbors by bonds of length r = 1. The nodes are considered to be of zero size. At the start of a simulation the triangular lattice has a size of I x lfi/2 and the nodes can be mapped onto an 1x 1square lattice. There are l2 nodes and 312 bonds in the initial configuration. Each of the bonds is considered to be a hookean spring with an equilibrium length roe The force cr associated with each bond is given by 0 = k(r-r,)
(1)
where k is the force constant. For the simulations reported in this paper all the bonds are assumed to be equivalent and the equilibrium bond length is 0.90909 (an initial extension ratio of 1.1). This initial extension ratio is quite large but is not beyond the range found in real systems. To represent the interaction with the underlying substrate, each node is connected to its original position by a weak “bond” and the force associated with this bond for the ith node is given by 0 = k2(Si -&“)
(2)
here Si is the position of the ith node and 4” is its initial position. The force constant k, is much smaller than k (typically k2/k = 0.01). The initial configuration in a simulation is illustrated in Fig. 1. In this model the bonds in the surface layer are allowed to break. The “weak” bonds which connect the surface to the substrate are not broken. It would be easy to allow both types of bonds to break but this was not done in this initial investigation to reduce the number of parameters in the model to a minimum. After the original configuration has been set up bonds in the surface are picked at random and broken if a random number x uniformly distributed over the range 0 < x < 1 is smaller than Pi/P,,, where Pi is the bond breaking probability for the ith (selected) bond and P,,,,, is the maximum value for all the bond breaking probabilities. In this model the bond breaking probability is given by Pi = PO exp{k(;i;)2]
(3)
where k, is Boltzman’s constant and T is the temperature. This form for the bond breaking probabilities is motivated by the absolute reaction rate theory of Eyring34 and it is assumed that the activation energy for breaking a bond is reduced by the strain energy k(ri - ro)2/2. Other models for fracture3’ indicate that the reduction in
168
P. MEAKIN
Fig. 1. A schematic representation of the model. The nodes (large dots) are connected by strong bonds to form a triangular lattice. Each node is joined to the underlying substrate by a weak bond (-- -) at its original position at the start of the simulation. Throughout the simulation the distance from the nodes to the underlying substrate is constant (only horizontal motion is allowed). This figure shows the original configuration in which each node is associated with six bonds.
the activation energy is either linear in the stress or quadratic in the stress (as in eqn. (3)). For polymers36 a reduction in the activation energy which is linear in the stress has been used almost universally3’j. In this event eqn. (3) must be replaced by Pi=P,eXp i
4r - ro> k~ I P
In these simulations surface cracking is assumed to take place under isothermal conditions and the activation energies are expressed in units of k, T (reduced units). Some simulations have been carried out using eqn. (4) but most of our results have been obtained using eqn. (3) which is more consistent with the use of a harmonic potential in the model. For most real systems it is unlikely that either eqn. (3) or eqn. (4) would be realistic. It is assumed that, once a bond has been broken, it is incapable of reforming. After each bond has been broken, the entire network of bonds and nodes is returned to mechanical equilibrium. Standard relaxation methods are used37 and the relaxation procedure is continued until the maximum change in the position of any of the nodes is less than O.Ol(l - r,,) during a cycle in the relaxation of the network. In order to reduce the computer time, extra relaxation steps were used for those nodes within ten bonds of the last broken bond. Over-relaxation was used in all the simulations. Block relaxation could be used to improve the speed of the simulations but was not used to obtain the results reported here. The algorithm outlined above is not explicitly time dependent but can be used to obtain kinetic information if the “time” is incremented by l/NP,,,., after each attempt (successful or otherwise) to break a bond, where N is the number of bonds remaining and Pm,, is the maximum bond breaking probability (i.e. the rate constant for breaking the most strained bond). In the model described above it is implicitly
SIMPLE
MODEL
FOR ELASTIC
FRACTURE
IN THIN FILMS
169
assumed that the network relaxes to mechanical equilibrium at a much faster rate than the bond breaking process. Because of visco-elastic interactions within the surface layer, within the substrate and between the surface layer and the substrate this condition might not be met. Consequently, simulations have been carried out in which the displacement of the nodes in a time 6t is given by
where SSi is the displacement of the ith node,A is the force on the ith node and C is a constant. Simulations were carried out using eqn. (5) to update the coordinates of each node. The time increment 6t is made sufficiently small that the change SSi in the coordinate is only a small fraction (typically 0.01) of the displacement which would bring them into equilibrium with their neighbors. This crude integration of the equations of motion for the nodes is stable. More elaborate integration procedures are not justified because of the stochastic nature of the bond breaking process. 3.
RESULTS
Some results from a small-scale simulation carried out using the fast relaxation model described above are shown in Fig. 2. In this simulation 2500 nodes (and 7500 bonds) were used. The force constants k and k, were given values of 1600 and 16 respectively (unless indicated otherwise k,/k = 0.01 for all the simulations). A value of 1600 for k implies that an energy of 661.2ksT would be required to double the bond length. This might seem to be a high energy, but a typical C-C single bond has a force constant of about 5 mdyn A-’ (ref. 38) and a length of about 1.54 A 39. Such a bond would require about 1450ks T to double its length at 300 K provided that the bond remained harmonic at such large extension ratios. Figure 2 shows the network of bonds and nodes after 200,400,600 and 800 bonds have been broken. At first, only small defects (single broken bonds) are formed but soon a crack is formed. During the early stages of the simulation when the stress in the bonds at the crack tip is very high, the cracks propagate rapidly and the direction of propagation is determined by the bond which bears the largest stress. Under these conditions the crack grows very rapidly and linearly (Fig. 2(a)). As the first few cracks propagate, they relieve the stress in the surface layer. This causes the rate of crack growth to slow down and the new cracks become less linear (Fig. 2(b)). At this stage (after 400 bonds have been broken in this simulation), new features have appeared which resemble shear bands. The formation of these shear bands appears to be a precracking phenomenon in the sense that the shear bands evolve into cracks at later stages in the simulation (Fig. 2(c)). However, shear bands do not appear at the very early stages of the simulation. As the stress is relieved, new cracks are formed but they grow only relatively slowly (Fig. 2(d)), and at still later stages when most of the stress has been relieved by cracking the bond breaking process becomes much more random. Pictures such as those shown in Fig. 2(a) which shows “snapshots” of the whole network allow all the features of the system to be displayed (however, there are distortions owing to the graphics devices used). For large-scale simulations a very large number of coordinates would have to be stored and the individual bonds and nodes would not be visible. Instead, the positions of only the broken bonds in the
170
P. MEAKIN
initial (triangular lattice) network are displayed. If these bonds are displayed in the order in which they were broken, the whole sequence of crack initiation and growth events can be seen. This is a very efficient and effective way of displaying the simulation results (particularly in the form of a movie or video recording) but information concerning features such as the crack widths and distortions of the network are lost. Figure 3 shows some results obtained from large-scale simulations (40 000 and 90 000 nodes, 120000 and 270 000 bonds) at four different values for the force constant k (400,800, 1200 and 1600). The results shown here and in the rest of this paper are the original positions of the broken bonds.
(4
(b) Fig. 2 (continued).
SIMPLE MODEL FOR ELASTIC FRACTURE
IN THIN FILMS
171
(4 Fig. 2. Results obtained from a small-scale simulation (2500 nodes) with k = 1600 and k, = 16. The system is shown after (a) 200, (b) 400, (c) 600 and (d) 800 bonds have been broken. In (b) three of the features which resemble shear bands are indicated by the letter S.
Figure 4 shows some results obtained from a 40 000 node simulation using the parameters k = 800 and k, = 8. After 1000 bonds have been broken (Fig. 4(a)), two cracks have been initiated and at this stage they are quite small. The time required to break the first 1000 bonds is 2.36 x 1O-4 in the reduced time units of the simulation. The next 1000 bonds are broken in only 1.6 x 10e5 time units (Fig. 4(b)) and there are now three large cracks. When one crack approaches another crack which has already grown, its rate of propagation slows down because the stress in the region near to the second crack has been reduced. The next 2000 bonds are broken in 3.29 x lop4 time units and a network of cracks has been formed (Fig. 4(c)). At this
172
P. MEAKIN
IO 000
@I
1 I=200-
BONDS
k=gOO
BROKEN
I
ke=8
Fig. 3. The effects of varying the elastic constant k for the strong bonds on the simulation results. It should be noted that the configuration of the network of bonds and nodes at mechanical equilibrium depends on the ratio k/k, but does not depend on the magnitude of k and k,. The magnitude of k does, however, determine the bond breaking probabilities via eqn. (3). In these simulations k/k, was fixed at a
SIMPLE MODEL
5000
FOR ELASTIC
BONDS
FRACTURE
173
IN THIN FILMS
BROKEN
’
(4 II0000
BONDS
BROKEN
I
value of 100 and k was given values of 200, 400, 800 and 1600 in (a), (b), (c) and (d) respectively. The simulations shown in (a)-(c) were carried out using 40 000 nodes (120 000 bonds) and that in(d) (k = 1600) was obtained using 90000 nodes (270000 bonds).
174
P. MEAKIN
1000
/
BONDS BROKEN I ‘I-. -,t- _,
\-
I
I
c
:.
‘-
I~
.’
,’
,.__!-
-
;
.’ ,.
,-
_
_.. I
1
,’ .I I.
(4
t ri.36
(b)
Fig. 4 (continued).
.,’
x 10’‘4
,,_
I,
SIMPLE
MODEL
4000
FOR ELASTIC
BONDS
FRACTURE
IN THIN FILMS
175
BROKEN
Fig. 4. Four stages in a simulation carried out using the parameters k = 800, k, = 8 and r,Jr,, = 1.1. In this simulation 40 000 nodes were used (I = 200) and the bond breaking probabilities are given by eqn. (3). The broken bonds (mapped onto the original triangular lattice) after 1000,2000,4000 and 8000 bonds have been broken are shown in (a),(b),(c) and (d) respectively.
P. MEAKIN
stage, most of the strain energy in the surface layer has been released and the rate of bond breaking is slowing down. Figure 4(d) shows the system after 800 bonds have been broken. On the scale used in these figures, it is not possible to distinguish between cracks and shear bands. Figure 5 shows an expansion of the central region in Fig. 4(d) and it can be seen that at this stage a large number of shear bands have appeared.
Fig. 5. An expansion of the central 80 x 80 region in Fig. 4(d). This figure shows that at this stage in the simulation the stress field has become very complex (no longer purely dilational) and that many shear bands have been formed. On the scale of Fig. 4(d) it is not possible to distinguish cracks from “shear bands”.
The number of bonds broken are shown as a function of time in Fig. 6(a). It appears that a slow crack initiation stage is followed by a period of rapid crack propagation which releases a significant fraction of the total strain energy in the surface (Fig. 6(b)). After this stage has been reached, new cracks are initiated but they propagate more slowly because the stress is smaller. Figure 6(b) shows the time dependence of the strain energy associated with the surface area (the quantity actually plotted here is
2E =
5 k(ri-roj2
i=l
(6)
where ri is the length of the ith bond and r. is the equilibrium bond length). This quantity is twice the strain energy. The time dependence of the strain energy also exhibits the same distinct stages as does the number of bonds broken.
SIMPLE MODEL FOR ELASTIC FRACTURE
6000
IN THIN FILMS
177
-
(4 600 g
700-
W 600L E 500 5
400-
g
JOO200
t x 10
,
-
loot ‘0
,
,
.X05
0010
, .0015
,
,
,
,
0020
9025 t
0030
0035
@I
, .GiMO
, .0045
j DO50
Fig. 6. (a) Time dependence of the number of broken bonds and (b) the total strain energy in the surface. layer obtained from the simulation illustrated in Figs. 4 and 5 (k = 800). These figures indicate a slow initiation period in which little strain energy is released and only isolated defects are formed, followed by a period of rapid crack propagation in which about 30% of the strain energy is released. This rapid release of the stress in the surface layer slows down subsequent crack initiation and propagation events. As more cracking occurs, the stress becomes still smaller and further slowing down is observed.
Figure 7 shows the time dependence of the number ofbonds broken for a similar simulation carried out using a bond force constant of 1600 instead of 800. Now the crack propagation rate is very rapid and a series of steps can be seen, corresponding to the growth of each crack. At later stages in the simulation there are usually several cracks growing at the same time. In a much larger-scale simulation or in a macroscopic sample there would always be many cracks growing simultaneously and the steps in Fig. 7 would be smoothed out into a continuous curve. 5000,
1000
I
H
I
r
Fig. 7. Time dependence of the number of bonds broken in a simulation which is similar to that used to obtain Figs. 4-6 except that the elastic constants k and k, have been doubled to 1600 and 16 respectively. Now crack propagation is very fast in the early stages of the simulation and each crack propagation event is seen as a vertical step in the curve. As the simulation proceeds, the cracks and the corresponding step sizes become smaller and less frequent.
178
P. MEAKIN
Since a single simulation can give misleading results the calculations used to obtain Figs. 6(a) and 6(b) were repeated several times. Figure 8 shows results obtained from seven simulations using the parameters I = 200, k = 800, k, = 80, 1.1. Since these simulations required quite large amounts of computer time rOIreq = (about 4 h of CPU time on an IBM 3081 to break 4000 bonds), it was not practical to average more than seven simulations.
(b) Fig. 8. The results obtained dependence for the number
t from seven simulations similar to that used to generate Figs. 4-6: (a) the time of bonds broken from each simulation;(b) the averaged results.
An important parameter in this node is the force constant k, for the bonds which attach the surface layer to the substrate. Figure 9(a) shows some results of a simulation carried out using the parameters 1 = 200, k = 800 and k2 = 40. At the stage depicted in this figure 10 000 bonds have been broken and 0.00222 time units have elapsed. Figure 9(b) shows an expansion of the region of size 80 x 80 in the center of Fig. 9(a). Similar results are shown in Figs. 9(c) and 9(d) except that the force constant k, has been reduced to 1.6. A comparison of these figures indicates that a low “strength” for the bonds joining the surface layer to the substrate results in cracks which are quite linear whereas a large value for k, gives rise to cracks which are much more crooked. It should also be noted that the time required to break 10000 bonds is an order of magnitude smaller for k, = 40 than for k, = 1.6. The reason for this is that the rapidly growing linear cracks combined with the weak attachment to the substrate result in a more rapid reduction of the strain in the surface layer for k, = 1.6 than for k, = 40. Figures 10(a) and 10(b) show the same simulations at the stage where 2000 bonds have been broken. In Fig. 10(a) (k2 = 40) a few short cracks have developed but most of the broken bonds are associated with more or less isolated defects. The elapsed time is 4.42 x 10d4 time units. For a much
179
SIMPLE MODEL FOR ELASTIC FRACTURE IN THIN FILMS
-
-% -,
(b) Fig. 9 (continued).
_
\,.,
180
P. MEAKIN
-r ,-,_-I
(4
-_I
__ ‘I
I
Fig. 9. The effects of changing the force constant k, describing the interactions between the and the substrate. These figures were obtained using the parameters I = 200 (40 000 nodes), r-o/r,, = 1.1. In (a) and(b) k, = 40 whereas in(c) and (d) k, = 1.6; (b) and(d) show expansions 80 x 80 regions of (a) and (c) respectively. These figures show that weak attachment of the leads to the formation of (fast propagating) linear cracks. In all of these figures 10 000 bonds broken.
surface layer k = 800 and of the central surface layer (8.3%) have
SIMPLE MODEL FOR ELASTIC FRACTURE
04
2000
BONDS
BROKEN
2000
BONDS
BROKEN
k = 800
k2=1.6
IN THIN FILMS
181
t = 0.000266
Fig. 10. The simulations shown in Figs. 9(a) and 9(c) at the stage where 2000 bonds (1.67%) have broken. This figure shows that weak “adhesion” between the surface layer and the substrate leads to more rapid and more catastrophic failure.
182
P. MEAKIN
smaller value of k, (weaker bonds joining the surface layer to the substrate) several very long cracks have developed (Fig. 10(b)). Despite the fact that an order of magnitude more time is required to break 10 000 bonds using this model (compared with the model with k, = 40) less time is required to break 2000 bonds. The reason for this is that the weaker bonding to the substrate allows a larger stress concentration at the crack tips. Figure 11(a) shows the time dependence of the number of broken bonds obtained from the simulations used to generate Figs. 9 and 10. At first, the rate of bond breaking is equal in both the models. However, the stress concentration at the crack tips leads to rapid crack propagation for the case of weak bonding between the surface layer and the substrate. This generates long linear cracks which are very effective at reducing the stress in the surface layer and after these cracks have grown the rate of crack initiation and growth slows considerably. The formation of long cracks such as those shown in Fig. 10(b) is normally much more catastrophic for most applications than the formation of short cracks such as those shown in Fig. 10(a). The effects of strong and weak bonding to the substrate are illustrated further in Fig. 10(b) which shows the time dependence of the strain energy obtained from both simulations. The simulations described above required quite large amounts of CPU time (typically about 10 h on an IBM 3081) because of the need to relax (almost) completely the network of nodes and bonds after each bond breaking event. In the case of the “visco-elastic” model even more computer time is required and it is not practical to carry out simulations on the same scale. Figure 12 shows the results of three simulations carried out using this model with the parameters I = 100, k = 1600
100 OO (4
(b)
_i ,005
.OlO
,015 +
I .020
I ,025
I ,050
1
Fig. 11. Time dependence of the number of broken bonds and the strain obtained from the simulations used to generate Figs. 9 and 10.
energy
in the surface layer
SIMPLE
MODEL
FOR ELASTIC
FRACTURE
4000
BONDS BROKEN
3348
BONDS BROKEN
( L=jdd k=160b (W Fig. 12 (continued)
183
IN THIN FILMS
k@160 kc=i06 t=i.16
x1ik4
t
184
1
1906
P. MEAKIN
BONDS
BROKEN
Fig. 12. The effects of the rate of mechanical relaxation on cracking of the surface layer. These results were obtained from simulations carried out using 10000 nodes (1 = 100) with k 1600, k, = 160 and using relaxation rate constants C of lo’, lo6 and 10’ rolr,, = 1.1. (a), (b) and (c) were obtained respectively. The quantity C is defined in the text. This figure shows that fast relaxation accelerates failure of the surface layer in this model.
and k, = 160. The relaxation rate constants C used were given by kC = 105, lo6 and 10’ to obtain the results shown in Figs. 12(a), 12(b) and 12(c) respectively. At each time interval the node displacements are calculated using eqn. (5) and the time increments 6t were set equal to O.Ol/Ck so that each adjustment of the coordinates of a node moves it only a small fraction of its distance from local mechanical equilibrium. For the simulation used to obtain Fig. 12(c) O.Ol/kC = 10e9 so that the number of steps made by each node in the network during the simulation is 109t or 2.8 x lo4 (t is the time at the end of the simulation). Simulations have also been carried out in which the reduction in the activation energy which defines the bond breaking probabilities is linearly proportional to the magnitude of the strain (eqn. (3)). Figure 13 shows some results obtained using the model parameters 1 = 200 (40 000 nodes) and k,/k = 0.01 for three different values of the parameter u. In Fig. 13(a), which was obtained using a relatively small value for v, only very short cracks have been formed after 6000 bonds have been broken. To obtain Fig. 13(b) u was doubled to 100. Now an extensive network of cracks has formed after 6000 bonds have been broken. The first cracks to be formed are long and relatively straight but they are not completely linear. The parameter u was doubled again to obtain Fig. 13(c). During the earlier stages of this simulation isolated defects are formed. Then long almost completely linear cracks grow very
SIMPLE MODEL FOR ELASTIC FRACTURE
1 6000
BONDS
IN THIN FILMS
BROKEN
(4
L= 200
V= i.00
(b) Fig. 13 (continued).
t = 8.6 ~10’~
kp/i=~O.Ol
185
P. MEAKIN
Fig. in which in activation energy is linearly dependent on the bond I = 200, r,/r,, = 1.1 and k,/k = 0.01 with three different values for the acceleration parameter v defined in eon. (4). (a), (b) and (c) show results obtained with u = 50, u = 100 and u = 200 respectively.
rapidly. During the later stages of the simulation these cracks become connected by secondary cracks which become less and less linear as the strain in the surface layer is reduced. During the stage just after the formation of the first few cracks very long “shear bands” are formed which soon develop into cracks. It is apparent from Figs. 3 and 13 that the parameters k in eqn. (3) and u in eqn. (4) have quite similar qualitative effects. Figure 14 shows the distribution of forces (T in the remaining bonds in the network. These results were obtained from the set of simulations used to generate Fig. 8 (1 = 200, k = 800, k, = 80 and rO/req = 1.1). Here N(lng) is defined by N(ln a)6 In csbeing the number of bonds for which the natural logarithm of the force lies in the range In CJto In CJ+ 6 In (r. The logarithmic distribution of forces is shown after N, = 200,400, 800,160O and 3200 bonds have been broken. If N, < 500 (see Fig. 8), only isolated defects exist and the force distribution is sharply peaked about the initial force go (indicated by a vertical arrow in Fig. 14(a)). For larger values of NB extended defects (cracks) have formed and there is a very broad distribution of forces in the bonds forming the network. As the number N, of broken bonds increases, the peak in the distribution moves slowly to the left (toward lower forces) and the number of bonds with very small forces increases rapidly. At early times in the simulation all the bonds are in tension (cr is positive). However, after about 500 bonds have been broken (cracks have formed) a few of the bonds have lengths
SIMPLE MODEL FOR ELASTIC FRACTURE
IN THIN FILMS
15, 14 13 12 11 10
187
I,
y c9 z:s +'
6 5 4 2 -10
-9
-6
-7
-6
-5
-4
-3
-2
7,,,,,,,,,
-1
,
,
-3
-2
6?
5-
?4t 5 321 -12
-11
-10
-9
-6
-7
-6
-5
-4
Fig. 14. Distribution of forces c in a network (with k = 800 and k, = 8) after N, = 200,400,800, 1600 and 3200 bonds have been broken. At the start of the simulation rO/req = 1.1. In (a) the force distributions in the bonds with r-/r_ > 1 (stretched bonds) are shown and in(b) the distributions for those bonds with r/r.s < 1 (compressed bonds) are shown. N(ln e) 6 In c is the number of bonds for which the force lies in the rangeInutoIna+6lncr.
smaller than r. and the distribution Fig. 14(b). 4.
of forces among these bonds is shown in
DISCUSSION
A very simple model for the cracking of thin films has been presented. This model is too simple to represent the rich phenomenology found in most real systems. Consequently, no detailed comparisons with experimental results obtained for specific systems have been attempted. However, the model does contain some of the key ingredients found in almost all systems in which a surface layer has been deposited onto a substrate (a rate of failure which is sensitive to the local stress field combined with complex stress concentration effects due to the presence of many defects of different types). Structures which resemble isolated defects, shear bands and cracks are clearly seen in the simulations. In many cases, the isolated defects develop into shear bands which later provide pathways for crack propagation. In many of our simulations a slow crack initiation period is followed by a period of very fast crack propagation which in turn gives way to a period of new crack initiation and propagation on a longer time scale (because of the reduction in stress in the surface layer as a result of the cracking process). If the strain is large (and/or the parameters k in eqn. (3) or u in eqn. (4) are large), the first non-localized defects to be formed are straight cracks. Only after these cracks have formed and
188
P. MEAKIN
destroyed the symmetry of the initial strain do shear bands develop as the shear stress grows in local regions in the surface layer. Also, in the initial stages of the simulations in which the strain is high (and purely dilational), the cracks which are formed are linear. At later stages when the stress at the crack tips is smaller and the stress field is much more complex, the cracks become less and less linear. If the magnitude of the initial stress is small, the bond breaking process becomes more random (and much slower). Under these conditions crack propagation occurs at a later stage in the simulations and the cracks which do form are not linear. In the low strain limit we would have an essentially random bond breaking (percolation) process and the large cracks which would eventually form could be described in terms of a fractal dimensionality ‘a D of about 1.98. In most of our simulations we are in a crossover region between D = 1 (linear cracks) and D = 1.98 (random cracking). Also, the distribution of defect sizes is an interesting quantity which we have not examined quantitatively. For large surface strain we have a large number of isolated defects and a few very long cracks at the point of failure (the point at which the surface layer becomes divided into two pieces). For the random bond breaking process (percolation), it is well known 4o that at the point of fracture (the percolation threshold) the distribution of defect sizes can be described by a power law with an exponent r of about -2.05 in two dimensions. Before failure has occurred, the distribution ofdefects can be described as a power law (with the same value of r) with an exponential cut-off 41 : N,cc S’exp(-as)
(7)
where N, is the number of defects containing s broken bonds. In our intermediate case the situation is more complex. For example, the broken bonds in a shear band are not nearest neighbors but the shear band is clearly a well-defined extended defect. Depending on the “rules” chosen, a shear band could be considered as either a large defect containing s broken bonds or s defects each containing one bond. An investigation of these questions is deferred until later. Here we have considered only the case of an applied tensile stress in the surface film, There are several reasons for this. The most familiar examples of surface cracking do involve tensile stresses (the cracking of paint films, surface cracking in rapidly cooled ceramics, mud cracking etc.) Cracking as a result of compressive stresses is also common. However, failure mechanisms which are not well represented by our two-dimensional model are likely to occur in this case (buckling and separation from the underlying substrate etc.). Since our model for the stressstrain field is linear for small displacements, very similar results would be expected for compressive forces2’. There are many ways in which the model could be made richer and more realistic. However, in this initial study an attempt was made to keep the model as simple as possible and to minimize the number of parameters. In some simulations we observed that the propagation of a crack can be arrested by an encounter with an isolated defect. Clearly, an investigation of the effect of pre-existing defects on the model would be interesting. It is also possible that heterogeneities in the bond lengths, bond strengths and force constants could have a similar effect (or they may accelerate failure). It is also apparent that a more realistic representation of the
SIMPLE
MODEL
FOR ELASTIC
FRACTURE
IN THIN FILMS
189
bonding between the surface layer and the substrate would be desirable. The weak bonds joining the surface layer to the substrate could be made capable of failure (as well as the bonds in the surface layer itself). A multilayer model would make a more realistic representation of the surface, interlayer (if any) and the substrate possible. It should be noted that the failure mechanism has not been specified in the model. The bond breaking process could be due to purely physical and/or chemical processes. It would also be possible to add the transport of chemically active species to the model. Another simple modification of the model would be to use a different lattice in order to represent materials with different elastic coefficients (Lame coefficients42). Finally, the possibility of forming3” as well as breaking bonds could be’ included so that inelastic effects could be represented. This might be particularly important for the “weak” bonds joining the surface layer to the substrate. The relaxation algorithm used in connection with this work causes the stress and strain fields to approach equilibrium via a “diffusion” process. Consequently, the number of relaxation steps required for the influence of a newly created defect to propagate through the entire system is on the order of 1’ (104-105). Since the total number of relaxation steps in a typical simulation also lies in the range 104-105, it might seem surprising that these simulations succeed at all. In the early stages of the simulation only isolated defects are formed. In a two-dimensional system these defects are (weakly) localized (i.e. their influence on the stress field decays asymptotically as rm2 where r is the distance measured from the defect). Consequently, the stress field does not need to diffuse over very large distances to approach mechanical equilibrium closely. In fact, for isolated defects the extensive relaxation carried out in the vicinity of a broken bond brings the system close to equilibrium. The relaxation of long cracks (extended defects) is potentially much more dangerous. However, the “weak” bonding between the surface layer and the substrate strongly localized the effects of the crack on the stress field (the effect of the crack on the stress field decays exponentially with increasing distance from the crack). The decay length for the decay of the perturbation of an infinitely long linear crack is given by C(k/k2)“2 where C is a coefficient of order unity. In fact, the bonding between the substrate and the surface layer strongly (exponentially) localizes all defects. Without this localization and the use of over-relaxation, the simulations described above would not succeed. ACKNOWLEDGMENTS
I would like to acknowledge many helpful discussions with Yves Termonia during the course of this work. The development of a model of this type to simulate failure in paint films was suggested by J. Hochberg.
REFERENCES
1 2
Paints and Protective Coatings. Army TM 5-618, NA VFAC MO-IlO, Air Force AFM-85-3, (U.S. Department of Defence, Department of the Army, Navy and Air Force). G. A. Sodenberg, Finishing Materials and Methods, McKnight and McKnight, Bloomington, 1952.
1969 IL,
P. MEAKIN
4
6 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34
35
36 31 38 39 40 41 42
P. Nylen and E. Sunderland, Modern Surface Coatings: A Textbook of the Chemistry and Technology of Paints, Varnishes and Lacquers, Wiley-Interscience, New York, 1965. C. R. Martens, Technology of Paints, Varnishes and Lacquers, Rheinhold, New York, 1968. C. R. Bragdon, Film Formation, Film Properties and Film Deterioration, Wiley-Interscience, New York, 1958. W. Namgoony and J. S. Chan, Thin Solid Films, 120 (1984) 153. C. A. Snavely and D. A. Vaughan, J. Am. Chem. Sot., I7 (1949)3 13. A. J. Becker and J. H. Blanks, Thin Solid Films, 119 (1984) 241. T. A. Witten and L. M. Sander, Phys. Rev. Lett., 47(1981) 1400. L. Niemeyer, L. Pietronero and H. J. Wiesmann, Phys. Rev. Lett., 52 (1984) 1033. P. Meakin, J. Theor. Biol., 118 (1986) 101. F. Family and D. P. Landau (eds.), Kinetics of Aggregation and Gelation, Elsevier, Amsterdam, 1984. H. E. Stanley and N. Ostrowsky (eds.), On Growth and Form: Fractal and Non-Fractal Patterns in Physics, NATO Advanced Study Institute, Cargese, France, 1985, Nijhoff, Dordrecht, 1986. L. Pietronero and E. Tosatti (eds.), Proc. 6th Int. Symp. of the International Centre for Theoretical Physics, Fractals in Physics, North-Holland, Amsterdam, 1986. L. A. Turkevich and H. Scher, Phys. Rev. Lett., 55 (1985) 1026. R. C. Ball, R. M. Brady, G. Rossi and B. R. Thompson, Phys. Reo. Left., 55 (1985) 1406. T. C. Halsey, P. Meakin and I. Procaccia, Phys. Reo. Left., 56 (1986) 854. B. B. Mandelbrot, The Fractal Geometry of Nature, Freeman, San Francisco, CA, 1982. A. G. Atkins and Y. W. Mai, Elastic and Plastic Fracture: Metals, Polymers, Ceramics, Composites, Biological Materials, Ellis Harwood, Chichester, 1985. S. N. Zhurkov, Int. .I. Fract. Mech., 2 (1965) 210. S. N. Zhurkov and E. E. Tomashevsky, Proc. Conf on Physical Basis of Yield and Fracture, Institute of Physics, London, 1966. S. N. Zhurkov, V. S. Kuksenko, S. I. Veliev, M. A. Gezalov and M. P. Vershina, Proc. Conf on the Yield, Deformation and Fracture of Polymers, Cambridge, 1970. F. Beuche, J. Appl. Phys., 28 (1957) 784. F. Bueche, J. Appl. Phys., 29 (1958) 1231. F. Bueche and J. C. Halpin, J. Appl. Phys., 35 (1964) 36. J. C. Halpin, J. Appl. Phys., 35 (1964) 3133. E. Louis and F. Guinea, Europhys. Lett., 3 (1987) 871. Y. Termonia and P. Meakin, Nature (London), 320 (1986) 6061. Proc. Conference on Fragmentation Form and Flow in Fractured Media, in Ann. Israel Phys. SOC., 78 (1986). B. B. Mandelbrot, D. E. Passoja and A. J. Paullay, Nature (London), 308 (1984) 721. D. Avnir, D. Farin and P. Pfeiffer, Nature (London), 308 (1984) 261. D. Avnir, D. Farin and P. Pfeiffer, J. Colloid Interface Sci., IO3 (1985) 112. Y. Termonia, P. Meakin and P. Smith, Macromolecules, I8 (1985) 2246; 19 (1986) 1.54. S. Glasstone, K. J. Laidler and H. Eyring, The Theory ofRate Processes, McGraw-Hill, New York, 1941. T. Ree and H. Eyring, in F. R. Eirich (ed.), Rheology, Vol. 2, Academic Press, New York, 1958, p. 83. E. R. Fuller and R. M. Thomson, in R. C. Bradt, D. P. H. Hasselman and F. F. Lange (eds.), Fracture Mechanics of Ceramics, Vol. 4, Crack Growth and Microstructure, Plenum, New York, 1973. A. J. Kinloch and R. J. Young, Fracture Behauior of Polymers, Applied Science, London, 1983. D. M. de G. Allen, Reiaxation Methods, McGraw-Hill, New York, 1954. E. B. Wilson, Jr., J. C. Decius and P. C. Cross, Molecular Vibrations, McGraw-Hill, New York, 1955. Handbook of Chemistry and Physics, CRC Press, Cleveland, OH, 53rd edn., 1972, p. F179. D. Stauffer, Introduction to Percolation Theory, Taylor and Francis, London, 1985. D. Stauffer, Phys. Rev., 54 (1979) 1. D. P. Landau and E. M. Lifshitz, Theory ofElasficity, Pergamon, Oxford, 1975.
A SIMPLE
AND
165
PROTECTIVE COATINGS
MODEL
FOR
ELASTIC
FRACTURE
IN THIN
FILMS
PAUL MEAKIN
Central Research and Development Department *, E. I. du Ponl de Nemours and Company, Wilmington, DE 19898 (U.S.A.) (Received
September
5, 1986; accepted
March 25, 1987)
A simple model for the cracking of thin films has been developed and investigated. The film is represented by a network of bonds and nodes which initially form a triangular lattice in which each node is at the junction of six bonds. The nodes are also attached to a rigid substrate by bonds which are not broken but have a relatively small force constant. Cracking is simulated by a sequence of thermally activated bond breaking events followed by mechanical relaxation. If the stretching force constant associated with the bonds in the surface layer is large, linear cracks propagate rapidly and, at a later stage, are connected by secondary cracks which become less and less linear as the strain in the surface layer is relaxed by the cracking process. Many of our simulations exhibit a distinct “initiation” period in which the bond breaking rate is low, followed by a period of rapid crack propagation in which a large fraction of the surface strain energy is released. This period is then followed by a second period of relatively slow bond breaking. During the first period isolated defects are formed, followed by linear cracks in the second period. During the third period non-linear cracks connect the linear cracks, and structures which resemble shear bands are formed. Results for the effects of finite relaxation rates for the elastic network (visco-elastic effects) are also presented.
1. INTRODUCTION
Thin films are now finding a broad variety of important applications including erosion and corrosion protection under extreme conditions, optical and electronic devices, wear reduction and a wide range of aesthetic applications. In many of these applications cracking of the surface layer is a serious form of failure which limits reliability and overall performance. One of the most familiar examples of this type of failure is the cracking and “checking” of paint films’-‘. This form of failure is now relatively uncommon (but often catastrophic when it does occur) in most organic coating systems. In other cases, cracking is a serious problem which inhibits the adoption of otherwise attractive new technologies. Cracking of the surface layer can * Contribution
number
004~6090/87/$3.50
4176.
0 Elsevier Sequoia/Printed
in The Netherlands
166
P. MEAKIN
lead to loss of required optical, electronic or mechanical properties in the surface layer itself or failure of the substrate as a result of corrosion, stress concentration etc. A key ingredient in the fracture of thin films which have been deposited onto more or less rigid substrates is the development of stress in the deposited film. This may be a result of solvent 10ss’-~, other changes in chemical compositioniU7, changes in temperature’, phase changes in the surface layer, deformation of the substrate etc. Another characteristic of thin film systems which makes them prone to failure by cracking is the (often very large) mismatch in properties between the surface layer and the substrate. Here a simple model is investigated which is intended to represent the case of a “hard” surface layer on a soft substrate or a hard surface layer bonded to a hard substrate by means of a layer of much softer material. Such systems include ceramic-coated metals, metal-coated polymers and paint systems in which the surface has become hardened as a result of oxidation and other (physical or chemical) changes. The growth of cracks in an elastic medium is similar to other non-equilibrium processes such as diffusion-limited aggregation (DLA). In fact, part of the motivation for this work was the similarity between crack growth and the DLA model of Witten and Sander’. This similarity becomes more apparent in the dielectric breakdown” and diffusion-limited growthl’ versions of the DLA model. In these models, the growth probability depends on the local value of a field 4 (electric potential in the dielectric breakdown model and concentration in the diffusion-limited growth model) near to the surface of the growing object. The field is obtained by solving Laplace’s equation (V’C$ = 0) subject to a fixed value of zero on the growing object and a constant value (say unity) on some boundary which encloses the growing object. In these stochastic growth models the growth probability is assumed to be some function of the local value of the field C$or its gradient. The DLA model and related non-equilibrium growth models have been studied extensively in recent years 9- I4 . Despite the apparent simplicity of these models and some new ideas”-’ 7, they are not yet well understood theoretically. The major difficulties seem to result from the fact that, whereas the growth probability at some point on the surface depends on the local value of the potential (or equivalently the potential gradient since 4 = 0 on the growing structure), the potential at each point depends on the entire structure which typically has a complex (often fractal18) geometry. The process of elastic fracture’ 9 can be thought of in similar terms. Here, we are concerned with a stress field 0 which cannot propagate across crack surfaces. As in the non-equilibrium growth and aggregation problem the determination of the local value of CTis a non-local problem made more difficult by the tensor nature of the stress field. As in the DLA and related models, we assume that the rate of crack propagation (and initiation) depends on the local stress field. Because of the theoretical difficulties associated with crack growth in systems with complex geometries, it seems probable that simple computer models (such as that described below) will prove to be valuable in guiding and stimulating theoretical work. In our model (described in detail below) the elastic medium is represented by a network of bonds and nodes. The bonds are broken randomly with probabilities which depend on the stress in each bond. After each bond breaking event, the system is relaxed to a new mechanical equilibrium and the new bond breaking probabilities
SIMPLE
MODEL
FOR ELASTIC
FRACTURE
IN THIN FILMS
167
are calculated. The model is kinetic in nature and assumes that bond breaking is a thermally activated process. In this respect, it is related to the early models of Zhurkov and coworkers20P22 and of Bueche and Halpin23P26. The model described below is also closely related to other models27P29 which have been used to explore the possible origins of fractal geometry” in fractal surfaces29-32 and failure in polymer fibers.33 2.
COMPUTER
MODELS
In all our simulations the surface layer is represented by an array of nodes and bonds. At the start of each simulation these nodes and bonds form a triangular lattice (with periodic boundary conditions) in which each node is connected to six nearest neighbors by bonds of length r = 1. The nodes are considered to be of zero size. At the start of a simulation the triangular lattice has a size of I x lfi/2 and the nodes can be mapped onto an 1x 1square lattice. There are l2 nodes and 312 bonds in the initial configuration. Each of the bonds is considered to be a hookean spring with an equilibrium length roe The force cr associated with each bond is given by 0 = k(r-r,)
(1)
where k is the force constant. For the simulations reported in this paper all the bonds are assumed to be equivalent and the equilibrium bond length is 0.90909 (an initial extension ratio of 1.1). This initial extension ratio is quite large but is not beyond the range found in real systems. To represent the interaction with the underlying substrate, each node is connected to its original position by a weak “bond” and the force associated with this bond for the ith node is given by 0 = k2(Si -&“)
(2)
here Si is the position of the ith node and 4” is its initial position. The force constant k, is much smaller than k (typically k2/k = 0.01). The initial configuration in a simulation is illustrated in Fig. 1. In this model the bonds in the surface layer are allowed to break. The “weak” bonds which connect the surface to the substrate are not broken. It would be easy to allow both types of bonds to break but this was not done in this initial investigation to reduce the number of parameters in the model to a minimum. After the original configuration has been set up bonds in the surface are picked at random and broken if a random number x uniformly distributed over the range 0 < x < 1 is smaller than Pi/P,,, where Pi is the bond breaking probability for the ith (selected) bond and P,,,,, is the maximum value for all the bond breaking probabilities. In this model the bond breaking probability is given by Pi = PO exp{k(;i;)2]
(3)
where k, is Boltzman’s constant and T is the temperature. This form for the bond breaking probabilities is motivated by the absolute reaction rate theory of Eyring34 and it is assumed that the activation energy for breaking a bond is reduced by the strain energy k(ri - ro)2/2. Other models for fracture3’ indicate that the reduction in
168
P. MEAKIN
Fig. 1. A schematic representation of the model. The nodes (large dots) are connected by strong bonds to form a triangular lattice. Each node is joined to the underlying substrate by a weak bond (-- -) at its original position at the start of the simulation. Throughout the simulation the distance from the nodes to the underlying substrate is constant (only horizontal motion is allowed). This figure shows the original configuration in which each node is associated with six bonds.
the activation energy is either linear in the stress or quadratic in the stress (as in eqn. (3)). For polymers36 a reduction in the activation energy which is linear in the stress has been used almost universally3’j. In this event eqn. (3) must be replaced by Pi=P,eXp i
4r - ro> k~ I P
In these simulations surface cracking is assumed to take place under isothermal conditions and the activation energies are expressed in units of k, T (reduced units). Some simulations have been carried out using eqn. (4) but most of our results have been obtained using eqn. (3) which is more consistent with the use of a harmonic potential in the model. For most real systems it is unlikely that either eqn. (3) or eqn. (4) would be realistic. It is assumed that, once a bond has been broken, it is incapable of reforming. After each bond has been broken, the entire network of bonds and nodes is returned to mechanical equilibrium. Standard relaxation methods are used37 and the relaxation procedure is continued until the maximum change in the position of any of the nodes is less than O.Ol(l - r,,) during a cycle in the relaxation of the network. In order to reduce the computer time, extra relaxation steps were used for those nodes within ten bonds of the last broken bond. Over-relaxation was used in all the simulations. Block relaxation could be used to improve the speed of the simulations but was not used to obtain the results reported here. The algorithm outlined above is not explicitly time dependent but can be used to obtain kinetic information if the “time” is incremented by l/NP,,,., after each attempt (successful or otherwise) to break a bond, where N is the number of bonds remaining and Pm,, is the maximum bond breaking probability (i.e. the rate constant for breaking the most strained bond). In the model described above it is implicitly
SIMPLE
MODEL
FOR ELASTIC
FRACTURE
IN THIN FILMS
169
assumed that the network relaxes to mechanical equilibrium at a much faster rate than the bond breaking process. Because of visco-elastic interactions within the surface layer, within the substrate and between the surface layer and the substrate this condition might not be met. Consequently, simulations have been carried out in which the displacement of the nodes in a time 6t is given by
where SSi is the displacement of the ith node,A is the force on the ith node and C is a constant. Simulations were carried out using eqn. (5) to update the coordinates of each node. The time increment 6t is made sufficiently small that the change SSi in the coordinate is only a small fraction (typically 0.01) of the displacement which would bring them into equilibrium with their neighbors. This crude integration of the equations of motion for the nodes is stable. More elaborate integration procedures are not justified because of the stochastic nature of the bond breaking process. 3.
RESULTS
Some results from a small-scale simulation carried out using the fast relaxation model described above are shown in Fig. 2. In this simulation 2500 nodes (and 7500 bonds) were used. The force constants k and k, were given values of 1600 and 16 respectively (unless indicated otherwise k,/k = 0.01 for all the simulations). A value of 1600 for k implies that an energy of 661.2ksT would be required to double the bond length. This might seem to be a high energy, but a typical C-C single bond has a force constant of about 5 mdyn A-’ (ref. 38) and a length of about 1.54 A 39. Such a bond would require about 1450ks T to double its length at 300 K provided that the bond remained harmonic at such large extension ratios. Figure 2 shows the network of bonds and nodes after 200,400,600 and 800 bonds have been broken. At first, only small defects (single broken bonds) are formed but soon a crack is formed. During the early stages of the simulation when the stress in the bonds at the crack tip is very high, the cracks propagate rapidly and the direction of propagation is determined by the bond which bears the largest stress. Under these conditions the crack grows very rapidly and linearly (Fig. 2(a)). As the first few cracks propagate, they relieve the stress in the surface layer. This causes the rate of crack growth to slow down and the new cracks become less linear (Fig. 2(b)). At this stage (after 400 bonds have been broken in this simulation), new features have appeared which resemble shear bands. The formation of these shear bands appears to be a precracking phenomenon in the sense that the shear bands evolve into cracks at later stages in the simulation (Fig. 2(c)). However, shear bands do not appear at the very early stages of the simulation. As the stress is relieved, new cracks are formed but they grow only relatively slowly (Fig. 2(d)), and at still later stages when most of the stress has been relieved by cracking the bond breaking process becomes much more random. Pictures such as those shown in Fig. 2(a) which shows “snapshots” of the whole network allow all the features of the system to be displayed (however, there are distortions owing to the graphics devices used). For large-scale simulations a very large number of coordinates would have to be stored and the individual bonds and nodes would not be visible. Instead, the positions of only the broken bonds in the
170
P. MEAKIN
initial (triangular lattice) network are displayed. If these bonds are displayed in the order in which they were broken, the whole sequence of crack initiation and growth events can be seen. This is a very efficient and effective way of displaying the simulation results (particularly in the form of a movie or video recording) but information concerning features such as the crack widths and distortions of the network are lost. Figure 3 shows some results obtained from large-scale simulations (40 000 and 90 000 nodes, 120000 and 270 000 bonds) at four different values for the force constant k (400,800, 1200 and 1600). The results shown here and in the rest of this paper are the original positions of the broken bonds.
(4
(b) Fig. 2 (continued).
SIMPLE MODEL FOR ELASTIC FRACTURE
IN THIN FILMS
171
(4 Fig. 2. Results obtained from a small-scale simulation (2500 nodes) with k = 1600 and k, = 16. The system is shown after (a) 200, (b) 400, (c) 600 and (d) 800 bonds have been broken. In (b) three of the features which resemble shear bands are indicated by the letter S.
Figure 4 shows some results obtained from a 40 000 node simulation using the parameters k = 800 and k, = 8. After 1000 bonds have been broken (Fig. 4(a)), two cracks have been initiated and at this stage they are quite small. The time required to break the first 1000 bonds is 2.36 x 1O-4 in the reduced time units of the simulation. The next 1000 bonds are broken in only 1.6 x 10e5 time units (Fig. 4(b)) and there are now three large cracks. When one crack approaches another crack which has already grown, its rate of propagation slows down because the stress in the region near to the second crack has been reduced. The next 2000 bonds are broken in 3.29 x lop4 time units and a network of cracks has been formed (Fig. 4(c)). At this
172
P. MEAKIN
IO 000
@I
1 I=200-
BONDS
k=gOO
BROKEN
I
ke=8
Fig. 3. The effects of varying the elastic constant k for the strong bonds on the simulation results. It should be noted that the configuration of the network of bonds and nodes at mechanical equilibrium depends on the ratio k/k, but does not depend on the magnitude of k and k,. The magnitude of k does, however, determine the bond breaking probabilities via eqn. (3). In these simulations k/k, was fixed at a
SIMPLE MODEL
5000
FOR ELASTIC
BONDS
FRACTURE
173
IN THIN FILMS
BROKEN
’
(4 II0000
BONDS
BROKEN
I
value of 100 and k was given values of 200, 400, 800 and 1600 in (a), (b), (c) and (d) respectively. The simulations shown in (a)-(c) were carried out using 40 000 nodes (120 000 bonds) and that in(d) (k = 1600) was obtained using 90000 nodes (270000 bonds).
174
P. MEAKIN
1000
/
BONDS BROKEN I ‘I-. -,t- _,
\-
I
I
c
:.
‘-
I~
.’
,’
,.__!-
-
;
.’ ,.
,-
_
_.. I
1
,’ .I I.
(4
t ri.36
(b)
Fig. 4 (continued).
.,’
x 10’‘4
,,_
I,
SIMPLE
MODEL
4000
FOR ELASTIC
BONDS
FRACTURE
IN THIN FILMS
175
BROKEN
Fig. 4. Four stages in a simulation carried out using the parameters k = 800, k, = 8 and r,Jr,, = 1.1. In this simulation 40 000 nodes were used (I = 200) and the bond breaking probabilities are given by eqn. (3). The broken bonds (mapped onto the original triangular lattice) after 1000,2000,4000 and 8000 bonds have been broken are shown in (a),(b),(c) and (d) respectively.
P. MEAKIN
stage, most of the strain energy in the surface layer has been released and the rate of bond breaking is slowing down. Figure 4(d) shows the system after 800 bonds have been broken. On the scale used in these figures, it is not possible to distinguish between cracks and shear bands. Figure 5 shows an expansion of the central region in Fig. 4(d) and it can be seen that at this stage a large number of shear bands have appeared.
Fig. 5. An expansion of the central 80 x 80 region in Fig. 4(d). This figure shows that at this stage in the simulation the stress field has become very complex (no longer purely dilational) and that many shear bands have been formed. On the scale of Fig. 4(d) it is not possible to distinguish cracks from “shear bands”.
The number of bonds broken are shown as a function of time in Fig. 6(a). It appears that a slow crack initiation stage is followed by a period of rapid crack propagation which releases a significant fraction of the total strain energy in the surface (Fig. 6(b)). After this stage has been reached, new cracks are initiated but they propagate more slowly because the stress is smaller. Figure 6(b) shows the time dependence of the strain energy associated with the surface area (the quantity actually plotted here is
2E =
5 k(ri-roj2
i=l
(6)
where ri is the length of the ith bond and r. is the equilibrium bond length). This quantity is twice the strain energy. The time dependence of the strain energy also exhibits the same distinct stages as does the number of bonds broken.
SIMPLE MODEL FOR ELASTIC FRACTURE
6000
IN THIN FILMS
177
-
(4 600 g
700-
W 600L E 500 5
400-
g
JOO200
t x 10
,
-
loot ‘0
,
,
.X05
0010
, .0015
,
,
,
,
0020
9025 t
0030
0035
@I
, .GiMO
, .0045
j DO50
Fig. 6. (a) Time dependence of the number of broken bonds and (b) the total strain energy in the surface. layer obtained from the simulation illustrated in Figs. 4 and 5 (k = 800). These figures indicate a slow initiation period in which little strain energy is released and only isolated defects are formed, followed by a period of rapid crack propagation in which about 30% of the strain energy is released. This rapid release of the stress in the surface layer slows down subsequent crack initiation and propagation events. As more cracking occurs, the stress becomes still smaller and further slowing down is observed.
Figure 7 shows the time dependence of the number ofbonds broken for a similar simulation carried out using a bond force constant of 1600 instead of 800. Now the crack propagation rate is very rapid and a series of steps can be seen, corresponding to the growth of each crack. At later stages in the simulation there are usually several cracks growing at the same time. In a much larger-scale simulation or in a macroscopic sample there would always be many cracks growing simultaneously and the steps in Fig. 7 would be smoothed out into a continuous curve. 5000,
1000
I
H
I
r
Fig. 7. Time dependence of the number of bonds broken in a simulation which is similar to that used to obtain Figs. 4-6 except that the elastic constants k and k, have been doubled to 1600 and 16 respectively. Now crack propagation is very fast in the early stages of the simulation and each crack propagation event is seen as a vertical step in the curve. As the simulation proceeds, the cracks and the corresponding step sizes become smaller and less frequent.
178
P. MEAKIN
Since a single simulation can give misleading results the calculations used to obtain Figs. 6(a) and 6(b) were repeated several times. Figure 8 shows results obtained from seven simulations using the parameters I = 200, k = 800, k, = 80, 1.1. Since these simulations required quite large amounts of computer time rOIreq = (about 4 h of CPU time on an IBM 3081 to break 4000 bonds), it was not practical to average more than seven simulations.
(b) Fig. 8. The results obtained dependence for the number
t from seven simulations similar to that used to generate Figs. 4-6: (a) the time of bonds broken from each simulation;(b) the averaged results.
An important parameter in this node is the force constant k, for the bonds which attach the surface layer to the substrate. Figure 9(a) shows some results of a simulation carried out using the parameters 1 = 200, k = 800 and k2 = 40. At the stage depicted in this figure 10 000 bonds have been broken and 0.00222 time units have elapsed. Figure 9(b) shows an expansion of the region of size 80 x 80 in the center of Fig. 9(a). Similar results are shown in Figs. 9(c) and 9(d) except that the force constant k, has been reduced to 1.6. A comparison of these figures indicates that a low “strength” for the bonds joining the surface layer to the substrate results in cracks which are quite linear whereas a large value for k, gives rise to cracks which are much more crooked. It should also be noted that the time required to break 10000 bonds is an order of magnitude smaller for k, = 40 than for k, = 1.6. The reason for this is that the rapidly growing linear cracks combined with the weak attachment to the substrate result in a more rapid reduction of the strain in the surface layer for k, = 1.6 than for k, = 40. Figures 10(a) and 10(b) show the same simulations at the stage where 2000 bonds have been broken. In Fig. 10(a) (k2 = 40) a few short cracks have developed but most of the broken bonds are associated with more or less isolated defects. The elapsed time is 4.42 x 10d4 time units. For a much
179
SIMPLE MODEL FOR ELASTIC FRACTURE IN THIN FILMS
-
-% -,
(b) Fig. 9 (continued).
_
\,.,
180
P. MEAKIN
-r ,-,_-I
(4
-_I
__ ‘I
I
Fig. 9. The effects of changing the force constant k, describing the interactions between the and the substrate. These figures were obtained using the parameters I = 200 (40 000 nodes), r-o/r,, = 1.1. In (a) and(b) k, = 40 whereas in(c) and (d) k, = 1.6; (b) and(d) show expansions 80 x 80 regions of (a) and (c) respectively. These figures show that weak attachment of the leads to the formation of (fast propagating) linear cracks. In all of these figures 10 000 bonds broken.
surface layer k = 800 and of the central surface layer (8.3%) have
SIMPLE MODEL FOR ELASTIC FRACTURE
04
2000
BONDS
BROKEN
2000
BONDS
BROKEN
k = 800
k2=1.6
IN THIN FILMS
181
t = 0.000266
Fig. 10. The simulations shown in Figs. 9(a) and 9(c) at the stage where 2000 bonds (1.67%) have broken. This figure shows that weak “adhesion” between the surface layer and the substrate leads to more rapid and more catastrophic failure.
182
P. MEAKIN
smaller value of k, (weaker bonds joining the surface layer to the substrate) several very long cracks have developed (Fig. 10(b)). Despite the fact that an order of magnitude more time is required to break 10 000 bonds using this model (compared with the model with k, = 40) less time is required to break 2000 bonds. The reason for this is that the weaker bonding to the substrate allows a larger stress concentration at the crack tips. Figure 11(a) shows the time dependence of the number of broken bonds obtained from the simulations used to generate Figs. 9 and 10. At first, the rate of bond breaking is equal in both the models. However, the stress concentration at the crack tips leads to rapid crack propagation for the case of weak bonding between the surface layer and the substrate. This generates long linear cracks which are very effective at reducing the stress in the surface layer and after these cracks have grown the rate of crack initiation and growth slows considerably. The formation of long cracks such as those shown in Fig. 10(b) is normally much more catastrophic for most applications than the formation of short cracks such as those shown in Fig. 10(a). The effects of strong and weak bonding to the substrate are illustrated further in Fig. 10(b) which shows the time dependence of the strain energy obtained from both simulations. The simulations described above required quite large amounts of CPU time (typically about 10 h on an IBM 3081) because of the need to relax (almost) completely the network of nodes and bonds after each bond breaking event. In the case of the “visco-elastic” model even more computer time is required and it is not practical to carry out simulations on the same scale. Figure 12 shows the results of three simulations carried out using this model with the parameters I = 100, k = 1600
100 OO (4
(b)
_i ,005
.OlO
,015 +
I .020
I ,025
I ,050
1
Fig. 11. Time dependence of the number of broken bonds and the strain obtained from the simulations used to generate Figs. 9 and 10.
energy
in the surface layer
SIMPLE
MODEL
FOR ELASTIC
FRACTURE
4000
BONDS BROKEN
3348
BONDS BROKEN
( L=jdd k=160b (W Fig. 12 (continued)
183
IN THIN FILMS
k@160 kc=i06 t=i.16
x1ik4
t
184
1
1906
P. MEAKIN
BONDS
BROKEN
Fig. 12. The effects of the rate of mechanical relaxation on cracking of the surface layer. These results were obtained from simulations carried out using 10000 nodes (1 = 100) with k 1600, k, = 160 and using relaxation rate constants C of lo’, lo6 and 10’ rolr,, = 1.1. (a), (b) and (c) were obtained respectively. The quantity C is defined in the text. This figure shows that fast relaxation accelerates failure of the surface layer in this model.
and k, = 160. The relaxation rate constants C used were given by kC = 105, lo6 and 10’ to obtain the results shown in Figs. 12(a), 12(b) and 12(c) respectively. At each time interval the node displacements are calculated using eqn. (5) and the time increments 6t were set equal to O.Ol/Ck so that each adjustment of the coordinates of a node moves it only a small fraction of its distance from local mechanical equilibrium. For the simulation used to obtain Fig. 12(c) O.Ol/kC = 10e9 so that the number of steps made by each node in the network during the simulation is 109t or 2.8 x lo4 (t is the time at the end of the simulation). Simulations have also been carried out in which the reduction in the activation energy which defines the bond breaking probabilities is linearly proportional to the magnitude of the strain (eqn. (3)). Figure 13 shows some results obtained using the model parameters 1 = 200 (40 000 nodes) and k,/k = 0.01 for three different values of the parameter u. In Fig. 13(a), which was obtained using a relatively small value for v, only very short cracks have been formed after 6000 bonds have been broken. To obtain Fig. 13(b) u was doubled to 100. Now an extensive network of cracks has formed after 6000 bonds have been broken. The first cracks to be formed are long and relatively straight but they are not completely linear. The parameter u was doubled again to obtain Fig. 13(c). During the earlier stages of this simulation isolated defects are formed. Then long almost completely linear cracks grow very
SIMPLE MODEL FOR ELASTIC FRACTURE
1 6000
BONDS
IN THIN FILMS
BROKEN
(4
L= 200
V= i.00
(b) Fig. 13 (continued).
t = 8.6 ~10’~
kp/i=~O.Ol
185
P. MEAKIN
Fig. in which in activation energy is linearly dependent on the bond I = 200, r,/r,, = 1.1 and k,/k = 0.01 with three different values for the acceleration parameter v defined in eon. (4). (a), (b) and (c) show results obtained with u = 50, u = 100 and u = 200 respectively.
rapidly. During the later stages of the simulation these cracks become connected by secondary cracks which become less and less linear as the strain in the surface layer is reduced. During the stage just after the formation of the first few cracks very long “shear bands” are formed which soon develop into cracks. It is apparent from Figs. 3 and 13 that the parameters k in eqn. (3) and u in eqn. (4) have quite similar qualitative effects. Figure 14 shows the distribution of forces (T in the remaining bonds in the network. These results were obtained from the set of simulations used to generate Fig. 8 (1 = 200, k = 800, k, = 80 and rO/req = 1.1). Here N(lng) is defined by N(ln a)6 In csbeing the number of bonds for which the natural logarithm of the force lies in the range In CJto In CJ+ 6 In (r. The logarithmic distribution of forces is shown after N, = 200,400, 800,160O and 3200 bonds have been broken. If N, < 500 (see Fig. 8), only isolated defects exist and the force distribution is sharply peaked about the initial force go (indicated by a vertical arrow in Fig. 14(a)). For larger values of NB extended defects (cracks) have formed and there is a very broad distribution of forces in the bonds forming the network. As the number N, of broken bonds increases, the peak in the distribution moves slowly to the left (toward lower forces) and the number of bonds with very small forces increases rapidly. At early times in the simulation all the bonds are in tension (cr is positive). However, after about 500 bonds have been broken (cracks have formed) a few of the bonds have lengths
SIMPLE MODEL FOR ELASTIC FRACTURE
IN THIN FILMS
15, 14 13 12 11 10
187
I,
y c9 z:s +'
6 5 4 2 -10
-9
-6
-7
-6
-5
-4
-3
-2
7,,,,,,,,,
-1
,
,
-3
-2
6?
5-
?4t 5 321 -12
-11
-10
-9
-6
-7
-6
-5
-4
Fig. 14. Distribution of forces c in a network (with k = 800 and k, = 8) after N, = 200,400,800, 1600 and 3200 bonds have been broken. At the start of the simulation rO/req = 1.1. In (a) the force distributions in the bonds with r-/r_ > 1 (stretched bonds) are shown and in(b) the distributions for those bonds with r/r.s < 1 (compressed bonds) are shown. N(ln e) 6 In c is the number of bonds for which the force lies in the rangeInutoIna+6lncr.
smaller than r. and the distribution Fig. 14(b). 4.
of forces among these bonds is shown in
DISCUSSION
A very simple model for the cracking of thin films has been presented. This model is too simple to represent the rich phenomenology found in most real systems. Consequently, no detailed comparisons with experimental results obtained for specific systems have been attempted. However, the model does contain some of the key ingredients found in almost all systems in which a surface layer has been deposited onto a substrate (a rate of failure which is sensitive to the local stress field combined with complex stress concentration effects due to the presence of many defects of different types). Structures which resemble isolated defects, shear bands and cracks are clearly seen in the simulations. In many cases, the isolated defects develop into shear bands which later provide pathways for crack propagation. In many of our simulations a slow crack initiation period is followed by a period of very fast crack propagation which in turn gives way to a period of new crack initiation and propagation on a longer time scale (because of the reduction in stress in the surface layer as a result of the cracking process). If the strain is large (and/or the parameters k in eqn. (3) or u in eqn. (4) are large), the first non-localized defects to be formed are straight cracks. Only after these cracks have formed and
188
P. MEAKIN
destroyed the symmetry of the initial strain do shear bands develop as the shear stress grows in local regions in the surface layer. Also, in the initial stages of the simulations in which the strain is high (and purely dilational), the cracks which are formed are linear. At later stages when the stress at the crack tips is smaller and the stress field is much more complex, the cracks become less and less linear. If the magnitude of the initial stress is small, the bond breaking process becomes more random (and much slower). Under these conditions crack propagation occurs at a later stage in the simulations and the cracks which do form are not linear. In the low strain limit we would have an essentially random bond breaking (percolation) process and the large cracks which would eventually form could be described in terms of a fractal dimensionality ‘a D of about 1.98. In most of our simulations we are in a crossover region between D = 1 (linear cracks) and D = 1.98 (random cracking). Also, the distribution of defect sizes is an interesting quantity which we have not examined quantitatively. For large surface strain we have a large number of isolated defects and a few very long cracks at the point of failure (the point at which the surface layer becomes divided into two pieces). For the random bond breaking process (percolation), it is well known 4o that at the point of fracture (the percolation threshold) the distribution of defect sizes can be described by a power law with an exponent r of about -2.05 in two dimensions. Before failure has occurred, the distribution ofdefects can be described as a power law (with the same value of r) with an exponential cut-off 41 : N,cc S’exp(-as)
(7)
where N, is the number of defects containing s broken bonds. In our intermediate case the situation is more complex. For example, the broken bonds in a shear band are not nearest neighbors but the shear band is clearly a well-defined extended defect. Depending on the “rules” chosen, a shear band could be considered as either a large defect containing s broken bonds or s defects each containing one bond. An investigation of these questions is deferred until later. Here we have considered only the case of an applied tensile stress in the surface film, There are several reasons for this. The most familiar examples of surface cracking do involve tensile stresses (the cracking of paint films, surface cracking in rapidly cooled ceramics, mud cracking etc.) Cracking as a result of compressive stresses is also common. However, failure mechanisms which are not well represented by our two-dimensional model are likely to occur in this case (buckling and separation from the underlying substrate etc.). Since our model for the stressstrain field is linear for small displacements, very similar results would be expected for compressive forces2’. There are many ways in which the model could be made richer and more realistic. However, in this initial study an attempt was made to keep the model as simple as possible and to minimize the number of parameters. In some simulations we observed that the propagation of a crack can be arrested by an encounter with an isolated defect. Clearly, an investigation of the effect of pre-existing defects on the model would be interesting. It is also possible that heterogeneities in the bond lengths, bond strengths and force constants could have a similar effect (or they may accelerate failure). It is also apparent that a more realistic representation of the
SIMPLE
MODEL
FOR ELASTIC
FRACTURE
IN THIN FILMS
189
bonding between the surface layer and the substrate would be desirable. The weak bonds joining the surface layer to the substrate could be made capable of failure (as well as the bonds in the surface layer itself). A multilayer model would make a more realistic representation of the surface, interlayer (if any) and the substrate possible. It should be noted that the failure mechanism has not been specified in the model. The bond breaking process could be due to purely physical and/or chemical processes. It would also be possible to add the transport of chemically active species to the model. Another simple modification of the model would be to use a different lattice in order to represent materials with different elastic coefficients (Lame coefficients42). Finally, the possibility of forming3” as well as breaking bonds could be’ included so that inelastic effects could be represented. This might be particularly important for the “weak” bonds joining the surface layer to the substrate. The relaxation algorithm used in connection with this work causes the stress and strain fields to approach equilibrium via a “diffusion” process. Consequently, the number of relaxation steps required for the influence of a newly created defect to propagate through the entire system is on the order of 1’ (104-105). Since the total number of relaxation steps in a typical simulation also lies in the range 104-105, it might seem surprising that these simulations succeed at all. In the early stages of the simulation only isolated defects are formed. In a two-dimensional system these defects are (weakly) localized (i.e. their influence on the stress field decays asymptotically as rm2 where r is the distance measured from the defect). Consequently, the stress field does not need to diffuse over very large distances to approach mechanical equilibrium closely. In fact, for isolated defects the extensive relaxation carried out in the vicinity of a broken bond brings the system close to equilibrium. The relaxation of long cracks (extended defects) is potentially much more dangerous. However, the “weak” bonding between the surface layer and the substrate strongly localized the effects of the crack on the stress field (the effect of the crack on the stress field decays exponentially with increasing distance from the crack). The decay length for the decay of the perturbation of an infinitely long linear crack is given by C(k/k2)“2 where C is a coefficient of order unity. In fact, the bonding between the substrate and the surface layer strongly (exponentially) localizes all defects. Without this localization and the use of over-relaxation, the simulations described above would not succeed. ACKNOWLEDGMENTS
I would like to acknowledge many helpful discussions with Yves Termonia during the course of this work. The development of a model of this type to simulate failure in paint films was suggested by J. Hochberg.
REFERENCES
1 2
Paints and Protective Coatings. Army TM 5-618, NA VFAC MO-IlO, Air Force AFM-85-3, (U.S. Department of Defence, Department of the Army, Navy and Air Force). G. A. Sodenberg, Finishing Materials and Methods, McKnight and McKnight, Bloomington, 1952.
1969 IL,
P. MEAKIN
4
6 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34
35
36 31 38 39 40 41 42
P. Nylen and E. Sunderland, Modern Surface Coatings: A Textbook of the Chemistry and Technology of Paints, Varnishes and Lacquers, Wiley-Interscience, New York, 1965. C. R. Martens, Technology of Paints, Varnishes and Lacquers, Rheinhold, New York, 1968. C. R. Bragdon, Film Formation, Film Properties and Film Deterioration, Wiley-Interscience, New York, 1958. W. Namgoony and J. S. Chan, Thin Solid Films, 120 (1984) 153. C. A. Snavely and D. A. Vaughan, J. Am. Chem. Sot., I7 (1949)3 13. A. J. Becker and J. H. Blanks, Thin Solid Films, 119 (1984) 241. T. A. Witten and L. M. Sander, Phys. Rev. Lett., 47(1981) 1400. L. Niemeyer, L. Pietronero and H. J. Wiesmann, Phys. Rev. Lett., 52 (1984) 1033. P. Meakin, J. Theor. Biol., 118 (1986) 101. F. Family and D. P. Landau (eds.), Kinetics of Aggregation and Gelation, Elsevier, Amsterdam, 1984. H. E. Stanley and N. Ostrowsky (eds.), On Growth and Form: Fractal and Non-Fractal Patterns in Physics, NATO Advanced Study Institute, Cargese, France, 1985, Nijhoff, Dordrecht, 1986. L. Pietronero and E. Tosatti (eds.), Proc. 6th Int. Symp. of the International Centre for Theoretical Physics, Fractals in Physics, North-Holland, Amsterdam, 1986. L. A. Turkevich and H. Scher, Phys. Rev. Lett., 55 (1985) 1026. R. C. Ball, R. M. Brady, G. Rossi and B. R. Thompson, Phys. Reo. Left., 55 (1985) 1406. T. C. Halsey, P. Meakin and I. Procaccia, Phys. Reo. Left., 56 (1986) 854. B. B. Mandelbrot, The Fractal Geometry of Nature, Freeman, San Francisco, CA, 1982. A. G. Atkins and Y. W. Mai, Elastic and Plastic Fracture: Metals, Polymers, Ceramics, Composites, Biological Materials, Ellis Harwood, Chichester, 1985. S. N. Zhurkov, Int. .I. Fract. Mech., 2 (1965) 210. S. N. Zhurkov and E. E. Tomashevsky, Proc. Conf on Physical Basis of Yield and Fracture, Institute of Physics, London, 1966. S. N. Zhurkov, V. S. Kuksenko, S. I. Veliev, M. A. Gezalov and M. P. Vershina, Proc. Conf on the Yield, Deformation and Fracture of Polymers, Cambridge, 1970. F. Beuche, J. Appl. Phys., 28 (1957) 784. F. Bueche, J. Appl. Phys., 29 (1958) 1231. F. Bueche and J. C. Halpin, J. Appl. Phys., 35 (1964) 36. J. C. Halpin, J. Appl. Phys., 35 (1964) 3133. E. Louis and F. Guinea, Europhys. Lett., 3 (1987) 871. Y. Termonia and P. Meakin, Nature (London), 320 (1986) 6061. Proc. Conference on Fragmentation Form and Flow in Fractured Media, in Ann. Israel Phys. SOC., 78 (1986). B. B. Mandelbrot, D. E. Passoja and A. J. Paullay, Nature (London), 308 (1984) 721. D. Avnir, D. Farin and P. Pfeiffer, Nature (London), 308 (1984) 261. D. Avnir, D. Farin and P. Pfeiffer, J. Colloid Interface Sci., IO3 (1985) 112. Y. Termonia, P. Meakin and P. Smith, Macromolecules, I8 (1985) 2246; 19 (1986) 1.54. S. Glasstone, K. J. Laidler and H. Eyring, The Theory ofRate Processes, McGraw-Hill, New York, 1941. T. Ree and H. Eyring, in F. R. Eirich (ed.), Rheology, Vol. 2, Academic Press, New York, 1958, p. 83. E. R. Fuller and R. M. Thomson, in R. C. Bradt, D. P. H. Hasselman and F. F. Lange (eds.), Fracture Mechanics of Ceramics, Vol. 4, Crack Growth and Microstructure, Plenum, New York, 1973. A. J. Kinloch and R. J. Young, Fracture Behauior of Polymers, Applied Science, London, 1983. D. M. de G. Allen, Reiaxation Methods, McGraw-Hill, New York, 1954. E. B. Wilson, Jr., J. C. Decius and P. C. Cross, Molecular Vibrations, McGraw-Hill, New York, 1955. Handbook of Chemistry and Physics, CRC Press, Cleveland, OH, 53rd edn., 1972, p. F179. D. Stauffer, Introduction to Percolation Theory, Taylor and Francis, London, 1985. D. Stauffer, Phys. Rev., 54 (1979) 1. D. P. Landau and E. M. Lifshitz, Theory ofElasficity, Pergamon, Oxford, 1975.