Propagation Mechanism of Brittle Cracks
We employ a recently developed model that allows the study of two-dimensional brittle crack propagation under fixed grip boundary conditions. The crack development highlights the importance of voids which appear ahead of the crack as observed in experiments on the nano-scale. The appearance of these voids is responsible for roughening the crack path on small scales, in agreement with theoretical expectations. With increasing speed of propagation one observes the branching instabilities that were reported in experiments. The simulations allow understanding the phenomena by analyzing the elastic stress field that accompanies the crack dynamics.
In laboratory experiments one studies crack propagation in elastic media by holding a stretched slab of material with a grip, and initiating the crack by making a small notch at one end of the sample 99HHMS (). When the notch exceeds the Griffith’s length 21Gri (), the crack propagates with increasing speed until there is a balance between the release of elastic energy and the creation of surface energy 98Fre (); 99MF (). It turned out that molecular dynamics simulations failed systematically to mimic this set up, forcing simulators to pull continuously on the boundaries to achieve a propagating crack 11HKL (). With fixed grip boundary conditions crack tended to slow down and stop propagating. The inability to simulate a brittle crack by molecular dynamics gave rise to claims that it would be impossible to advance a brittle crack without continuous stretching. The fundamental reason for this long-standing problem was understood recently. In Ref. 11DKPZ () it was shown that the solution of this conundrum lies in the range of the inter-particle potential. By changing the potential range one can go from a brittle to a ductile crack, and the latter is stopped by plastic dissipation. Only a brittle crack can support itself with a fixed grip. Brittleness was shown to be guaranteed by choosing a potential range that is of the order of the inter-particle distance. In this paper we present results of molecular dynamics simulations of brittle cracks under a fixed grip that are based on this understanding.
One of the interesting issues that can be studied by such simulations is the proposed existence of voids that open up before the crack, acting as pointers for further propagation. This issue had been quite controversial. While some experiments, and in particular those pertaining to ultra-slow crack propagation in glass in wet atmosphere, strongly indicated the appearance of voids ahead of the crack Bou (), other experiments failed to observe such voids, and claimed that they are irrelevant to the propagation of brittle cracks. In a theoretical study, Ref. 04BMP () constructed a model of void-dominated crack propagation and attempted to explain the observed roughness of brittle cracks by the randomness associated with the positioning of the voids ahead of the crack. The numerical simulations shown below appear to vindicate this approach almost verbatim. Voids do appear, and their appearance is random as stipulated in Ref. 04BMP (), and see Sect. 5 for more details.
Another issue of interest is the instabilities of the propagating crack when its velocity increases. The simulations show very clearly branching instabilities with subsequent competition between the two branches that typically result in the death of one branch in favor of the other.
In Sect. II we present the details of the preparation of the samples for the molecular dynamics simulations. In Sect. III we show that our cracks are brittle by demonstrating that they start running precisely when the Griffith’s criterion met. In Sect. IV we discuss the appearance of voids ahead of the crack and the resulting roughening of the crack path. Finally, in Sect. V we present a brief discussion of the observed instabilities. The last section presents a short summary and indicates the road ahead.
Ii Molecular Dynamics
For the numerical experiments we employ a generic glass former in 2-dimensions in the form of a 50-50 binary mixture of ‘small’ and ‘large’ particles, chosen to avoid any crystallization. In fact the particles interact by inter-particle potentials as shown in Fig. 1, with the analytic form
Here is the length where the potential attain it’s minimum, and is the cut-off length for which the potential vanishes. The coefficients and are chosen such that the repulsive and attractive parts of the potential are continuous with two derivatives at the potential minimum and the potential goes to zero continuously at with two continuous derivatives as well. The interaction length-scale between any two particles and is , and for two ‘small’ particles, one ‘large’ and one ‘small’ particle and two ‘large’ particle respectively. The unit of length is set to be the interaction length scale of two small particles, is the unit of energy and .
We used a modified Berendsen thermostat which couples a constant number of particles to the bath, regardless of the system size 10KLPZ (). The preparation protocol starts with equilibrating the sample at high temperature and pressure, using periodic boundary conditions. The systems are then cooled to temperature while keeping the pressure high, in order to avoid the creation of any holes in the material. Then the pressure is reduced to zero, , such that the periodic boundary conditions could be removed; particles forming the right and left walls were frozen, but the upper and lower boundaries were rendered free. The brittle crack experiment starts by loading the system uniaxially (with a constant velocity such that ) until a desired stress is reached, and then the side walls are held fixed. A cut is then implemented by the cancelation of forces crossing an imaginary line of desired length which starts at the lower boundary. The evolution of the crack is simulated by molecular dynamics coupled to a heat bath at and requires no further loading of the system. The results presented below were created using two different geometries, one of width 1000 and length 3000, and the other being a square sample of 1600. The density is the zero temperature and zero pressure density .
Iii The Griffith Length
Upon making the initial cut, the system always has a microscopic plastic response which however stops if the length of the cut is smaller than the Griffith length. To demonstrate that our running cracks are indeed brittle we test the Griffith’s criterion in our simulation.
The length of the initial cut that evolves spontaneously into a crack was determined by Griffith 21Gri (), comparing the energy release from the stress field to the energy consumed by the creation of two new fresh surfaces by the advancing crack. With the stress exerted on the slab being , the critical length is determined by the bulk modulus and the surface energy per unit length via the relation
In our simulations we are able to measure the critical length for a given loaded system. Since all the other material parameters appearing in Eq. 2 are measured independently, we can show that our model crack propagation agrees with the expected physics of brittle cracks.
The stress field is simply determined by our grip boundary conditions. The bulk modulus is read from the linear part of a uniaxial straining experiment shown in Fig. 3.
To estimate the surface energy in our system we computed the total energy before and after making the initial cut of length . Taking the difference and dividing by we get the estimate
At this point we employed 50 independent samples that were equilibrated at high temperature and quenched to the glassy state. The glassy samples were then loaded quasi-statically to strain values in the range . In order to probe the Griffith’s Length we tracked the behavior of our sample during the first 30 time units of our simulation to see whether a crack has started to run as a result of the introduced cut. For each sample we repeated this procedure for various values of . The results, (averaged over 50 samples) are presented in the bottom panel of Fig. 4. We show the length of the evolving crack, denoted as , as a function the initial cut size . Where the length begins to run determines for the conditions at hand.
In the upper panel of the same figure we plot vs. . The resulting linear plot is in agreement with Eq. (2), with the slope being close to the theoretical expectation of .
Iv Roughening resulting from voids ahead of the crack
The present numerical simulations allows the verification of the role of voids that form ahead of the crack, determining, at least on short length scales, the path chosen by the crack as it develops in the amorphous solid. We note that an amorphous solids is not an ideal elastic material, in which mathematically straight cracks are possible. The randomness of the material must show up in the geometry of the crack in this way or another. The voids ahead of the crack serve as pointers for the forthcoming propagation of the crack, and as is shown below, their random positions determines the roughening of the crack path.
In Fig. 5 we show a typical crack tip with the voids that opened ahead.
The physical reason why voids may prefer to open ahead of the crack and not on the crack tip was explained in Ref. 04BMP (). The argument is based on the fact that voids open due to plastic yield, and they do this where the pressure is maximal. Near the crack tip there is a process zone where the pressure is increasing going outward, until one hits the maximal pressure curve which connects with the outer elastic solution, see Fig. 6. In this figure we show the average pressure as measured in the numerical simulation in a fixed window tracking in time the crack tip.
In particular one should pay attention to the lower right panel in Fig. 6 which shows that the maximal pressure lies ahead of the tip, at a distance from the tip. This figure should be compared with Fig.2 of Ref. 04BMP () which it directly vindicates.
Obviously, the void that opens ahead of the tip can have a broadly distributed position around the forward direction (where the probability to form a void is maximal). An actual measurement of this distribution for the present simulation is shown in Fig. 7.
Shown are the actual positions of the voids as they appear ahead of the crack during the time evolution, and in the lower panels the probability distribution for the void to appear at an angle with respect to the forward direction, and at a distance along the forward direction. These data should be compared with Fig. 3 of Ref. 04BMP ().
Two points are worth stressing: first, the angular distribution of void positions will be the source of roughening of the crack - the probability to fall along the forward direction is not high enough compared to positive or negative angles with the respect to the forward direction. Second, once a void appears on, say, a positive angle, the next void will have an even higher probability for a positive angle, meaning that the rough crack is expected to show a roughening exponent higher than 1/2. Indeed, such a persistent random walk is always expected to show exponents higher than 1/2, whereas anti-persistent random walks are characterized by a roughening exponent smaller than 1/2.
The roughening exponent was measure here in the usual way, i.e. considering the crack as a graph and determining according to
For a self-affine graph the scaling exponent is defined via the scaling relation
We show this quantity in a log-log plot in Fig. 8. As expected the function exhibits a persistent scaling exponent of for scales , and then a cross over to a random graph without persistence or anti-persistence for higher scales.
Crack dynamics are governed by the balance of energy at the crack tip. The influx of elastic energy through the crack-tip is used to create surface energy behind the advancing tip. With increasing crack velocity there is not sufficient surface to store the energy that is released. Therefore the system resorts to instabilities in the form of branching and oscillations in order to increase the amount of surface created per unit length on the axis of propagation. In our simulations we find that when the velocity of the crack tip reaches about of the Rayleigh speed, one begins to observe crack branching. We note that our system never develops two independently propagating cracks, but rather reduced the velocity of propagation through attempted branchings. Oscillations were not observed. A typical branching event is shown in Fig. 9, where we see the two crack tips as they still grow simultaneously. The competition between them always results in the demise of one of the growing cracks in favor of the other, until the next branching event.
Vi summary and discussion
In summary, we have demonstrated that molecular dynamics can be usefully employed to study brittle crack propagation by selecting an appropriate interparticle potential. The resulting cracks are growing by nucleating voids ahead of the crack tip in much the same way that was anticipated by Ref. Bou () and theorized in Ref. 04BMP (). This mechanism of growth is responsible for roughening of the crack path on small scales, but as long as side branching does not commence, the crack is globally flat on macroscopic scales, as expected theoretically 06BBP ().
In future work molecular dynamics simulations will be employed in three dimensions where the issue of micro-branching 99MF () can be studied and compared with experimental results.
Acknowledgements.We thank Laurent Boué for useful discussion. This work had been supported in part by the German Israeli Foundation, the Israel Science Foundation and the European Research Council under an ”ideas” grant STANPAS.
- (1) See for example: J. A. Hauch and M. P. Marder, Int. J. Fract. 90, 133 (1998); J.A. Hauch, D. Holland, M.P. Marder and H. Swinney, Phys Rev, Lett, 92, 3823 (1999).
- (2) A.A. Griffith, Phil. Trans. Roy. Soc.London, 221, 163 (1921).
- (3) L. B. Freund, Dynamic Fracture Mechanics (Cambridge, London, 1998).
- (4) M. Marder and J. Fineberg, Phys. Rep. 313, 1 (1999).
- (5) S.I. Heizler, D.A. Kessler and H. Levine, Phys Rev E 84, 026102 (2011).
- (6) O. Dauchot, S. Karmakar, I. Procaccia and J. Zylberg, Phys. Rev. E 84, 046105 (2011)
- (7) E. Bouchaud, J. Phys. Condens. Matter 9, 4319 (1997); F. Ce´larie´, S. Prades, D. Bonamy, L. Ferrero, E. Bouchaud, C. Guillot, and C. Marlie‘re, Phys. Rev. Lett. 90, 075504 (2003).
- (8) E. Bouchbinder, J. Mathiesen and I. Procaccia, Phys. Rev. Lett., 92, 245505 (2004).
- (9) S. Karmakar, E. Lerner, I. Procaccia and J. Zylberg,Phys. Rev. E 82, 031301 (2010).
- (10) I. Ben-Dayan, E. Bouchbinder and I. Procaccia, Phys. Rev. E 74, 046102 (2006).