Fragmentation experiments

In order to understand the shape selection mechanism of fragments, we performed seven sets of experiments considering two cracking mechanisms which give rise to fragment formation: Weathering induced slow growth of cracks gradually removes pieces from the surface of rock walls under the action of gravity. Such surface fragments were populated either by collecting them from the foot of a cliff where they accumulated (experiment A), or by removing them from the mother rock by gently hammering the surface. The hammering experiments were repeated with limestone and dolomite at different geographical locations (experiments B, C, D,). Explosion experiments were performed in stone mines generating a rapid breakup of a rock wall into a large number of pieces (experiments E, F). Field measurements on geomaterials were complemented by laboratory experiments fragmenting gypsum cubes by manual hammering applied in a dynamic way (experiment G). Table 1 summarizes the most important parameters of the experiments, further details are presented in Methods. The number of fragments in the experiments was ranging from a few hundred to one thousand, providing reasonable statistics for the data evaluation. The maximal size L of the collected fragments was in the range 15 mm < L < 150 mm and we attempted to collect all available fragments in this range in a given area. For the quantitative characterization of the fragment ensembles we determined the mass of single pieces and characteristic quantities of the shape of fragments.

Table 1 Most important parameters of the seven experiments we performed. The table provides the type of fragmentation, the materials involved, the number of fragments (Num. Frag.) obtained and the size range of pieces in units of millimeter Full size table

The probability distribution of the mass of fragments p(m) is one of the most important characteristics of fragmentation phenomena1,2,5,12,14,15,16,17,18,19,20,21. Figure 2 demonstrates that for all sets of experiments p(m) has a power law functional form described by equation (1). The sampling method applied when collecting fragments implies cutoffs both for low and high mass values but in the intermediate mass regime of the distributions the power law form prevails in all cases. A very important feature of the results is that the distributions p(m) are described by a unique exponent τ = 1.7 ± 0.06 in spite of the strongly different cracking mechanisms governing fragment formation in dynamic fragmentation initiated by explosion and hammering and in weathering induced spallation. The value of the exponent τ falls close to the analytic prediction of the branching-merging scenario of dynamic cracks in three dimensions23,34,35: cracks propagating at a high speed undergo sequential branching and fragments are formed as material regions enclosed when sub-branches merge34,35. Assuming the self-similar nature of the developing crack-tree, power law distributed fragment masses are obtained with a universal exponent τ = (2D − 1)/D depending solely on the dimensionality D of the system. For three-dimensional bulk solids τ = 5/3 follows in a good agreement with our experimental findings. It is astonishing that the same universality class is recovered when the fragment population is generated by weathering (experiments A, B, C and D). The result can be explained by the fact that environmental aging typically gives rise to the growth of pre-existing cracks instead of creating new ones. Cracks slowly advance due to thermal expansion, hydration, frost induced expansion, etc, gradually removing pieces from the surface. Hence, the agreement with the prediction of the branching-merging scenario may imply that during their geological history the rocks considered have experienced dynamic loading.

Figure 2 Mass distribution of fragments p(m) obtained by experiments and by computer simulations. The mass m is normalized by the largest fragment mass m max of the data sets. The legend indicates how the fragment populations were generated, for more details of the experiments see Table 1 and Methods. The data sets were arbitrarily rescaled along the vertical axis to obtain collapse of the distributions. The bold line represents best fit with a power law followed by an exponential cutoff. Full size image

From the geometrical point of view, freshly created fragments are polyhedra with sharp edges and corners as it is illustrated in Fig. 1. As a first step we focus on the overall shape of rock pieces neglecting any surface texture created by the cracking mechanism. Later on the analysis will be refined by considering details of polyhedral shapes. In order to define an appropriate shape descriptor we construct the bounding box of each fragment with principal axis S < I < L (see Methods for details). The value of L provides a measure of the linear extension of fragments, while the dimensionless ratios S/L and I/L and their relations characterize the shape of pieces11,36,37,38,39,41,42,43.

Figure 3 presents the value of the shape parameters S/L and I/L as function of the fragment size L. The most remarkable feature of the results is that for all materials and fragmentation modes the same shape-size relation is obtained: for small values of L the shape parameters S/L and I/L rapidly approach one which implies that small sized fragments have an isotropic shape with S ≈ I ≈ L. As the fragment size L increases both S/L and I/L decrease and tend to finite constant values B S and B I , respectively. The size dependence of fragment shape in Fig. 3 can be well described by an exponential form

Here the cutoff length L c provides the fragment size where the isotropic shape is reached, while the scale parameter L 0 controls how fast the converges is to the unique shape of large fragments. Equations (2,3) also indicate that the characteristic length scales L 0 and L c have the same values for S and I. In Fig. 3 best fit was obtained with the parameter values A S = 0.45, A I = 0.4, B S = 0.43, B I = 0.67, L 0 = 7.2 mm and L c = 17.0 mm. The results imply that large enough fragments L ≫ L c + L 0 have a unique elongated shape characterized by the same value of the shape parameters S/L ≈ B S and I/L ≈ B I .

Figure 3 The dimensionless ratios S/L and I/L of the principal axis of fragments as a function of the value of the largest axis L. (The legend is the same as in Fig. 2.) For small fragment sizes both quantities approach one while for increasing L they decrease and converge to finite constant values indicated by the dashed horizontal lines. The bold lines represent fits with equations (2,3). Results of the model calculations are in a good agreement with the experimental findings. Full size image

Besides the shape-size relation of fragments the statistics of the occurrence of different shapes in a fragment ensemble is also of fundamental importance. The breakup models of NASA use the surface-to-volume ratio A/V to characterize the shape of objects generated by on-orbit fragmentation events30,31. In order not to have dependence on the extension of objects, following Ref. 32, we multiply this combination with the radius of gyration R g and define the shape parameter as S f = AR g /V. Assuming rectangular shape of fragments R g reads as and S f can be cast into the form . For cubic fragments S ≈ I ≈ L the shape parameter simplifies to S f = 3 independent of the size of fragments, while larger values S f > 3 characterize elongated shapes. To quantify the statistics of occurrence of different shapes in fragment populations, we determined the probability distribution p(S f ) of the shape parameter. In Fig. 4 again a robust behavior is obtained, i.e. in all experiments p(S f ) exhibits an exponential decay

The scale parameter proved to have a relatively low value which shows that the contribution of anisotropic shapes provides a considerable fraction to the fragment ensemble.

Figure 4 Probability distribution of the shape parameter S f on semilogarithmic plot. The straight line represents a fit with an exponential function. Both Model-Rect and Model-Poly reproduce the exponential functional form of p(S f ). Full size image

The dimensionless ratios S/L, I/L, as well as, the shape parameter S f carry rather limited information on the geometry of the fragment: for example, consider that at any fixed values for these parameters the fragment could be either ellipsoidal, tetrahedral, cuboid or even octahedral. To distinguish between these (and other) geometries with identical side ratios, we use the numbers n S and n U of stable and unstable static equilibrium points44 which can be counted in simple hand experiments and have been already used to assess geological data45. In the previous example, for ellipsoids, tetrahedra, cuboids and octahedra we get n S = 2, 4, 6, 8 and n U = 2, 4, 8, 6, respectively. The definition of stable and unstable equilibria is illustrated in Supplementary Fig. S1. Figure 5 presents the probability distributions of the number of stable p(n S ) and unstable p(n U ) points of fragments for all experiments. The remarkable result is that both n S and n U scatter over a broad range with a common minimum value of 2 and maxima reaching to n S = 9 and n U = 11, which demonstrate the richness of shapes not captured by other approaches. The data can be well fitted by a log-normal form with parameters µ = 1.54, σ = 0.22 and µ = 1.63, σ = 0.29, for the stable and unstable equilibrium points, respectively. (For details of the functional form of the distributions p(n S ) and p(n U ) see Supplementary Information.)

Figure 5 Probability distribution of the number of stable (a) and unstable (b) equilibrium points of fragments. The extended discrete model of fragmentation (Model-Poly) is able to reproduce the measured statistics of equilibria with a reasonable precision. Both the experimental and theoretical results are fitted with a lognormal distribution represented by the continuous and dashed lines, respectively. Full size image

Discrete stochastic model of fragmentation

In order to understand how the disorder of materials and the dynamics of disruption leads to the emergence of the observed scaling laws of fragment shapes, we constructed a generic stochastic model of fragmentation (Model-Rect). The model is based on a sequential binary breakup of fragments1,2,6,32,40 starting from a single body which has an initial cubic shape. At each step of the hierarchy fragments either break into two pieces of equal mass with probability 0 < p < 1 or they keep their current size until the end of the process with probability 1 − p. Iterating the stochastic dynamics for all active pieces a power law mass distribution of fragments is obtained. The model was calibrated by setting the value of the breaking probability to p = 0.8 which provides a very good agreement with the measured mass distribution of fragments in Fig. 2 (see also Methods and Supplementary Information).

To follow how the shape of pieces evolves, we assume that cracks always occur in the middle of edges cutting through the center of mass C of the body. Figure 6 illustrates that for any fragment (Fig. 6(a)) there are three possibilities for crack formation (Figs. 6(b, c, d)), which all result in two pieces of equal mass. (This breakup mechanism is in spirit the three-dimensional generalization of the fragmentation model of Refs. 41, 42.) To capture the effect that it is easier to break a body perpendicular to a longer extension, we choose among the three possibilities with a probability p c proportional to a power α of the length of the breaking edge l

This cracking mechanism retains the orthogonal form of pieces and gives rise to a complex dynamics of shape selection controlled by the exponent α, i.e. α = 0 implies completely random cracking while α > 0 favors the breaking of longer edges. Figures 3 and 4 demonstrate that setting α = 3 for the exponent of p c an excellent agreement is obtained with the experimental findings. The model reproduces the exponential convergence from the isotropic form of small fragments to a unique anisotropic shape of the large ones (Fig. 3) and it recovers the exponential distribution of the occurrence of different shapes (Fig. 4), as well.

Figure 6 Breaking mechanism of fragments in Model-Rect. At each iteration step fragments break into two pieces of equal mass such that the crack occurs perpendicular to one of the sides. After a fragment has been chosen for breaking, one of the three cases (b), (c), or (d) is selected with a probability depending on the length of the breaking edge. Full size image

In the case of orthogonal fragments (cuboids) considered above n S and n U have the values n S = 6 and n U = 8 which remain constant as the breakup proceeds so that this simple model cannot account for the geometric complexity encoded in the statistics of equilibria. To extend the model beyond orthogonal shapes we let the crack plane orientation vary continuously: the crack always crosses the center of mass of fragments, however, its normal vector n can take any direction according to a continuous probability distribution of the form of equation (5). In the model calculations the linear extension l of fragments is measured from the center of mass C and we create a list based on these distances from which our stochastic algorithm can uniformly sample. Each linear extension l i is associated with a unit vector n i . In case of a polyhedron P with F faces and V vertices we have several alternatives to create this list, out of which we mention three:

(a) the list l (f),i (i = 1, 2, … F) of the distances of planar faces, the associated unit vectors are the unit normals of the given face.

(b) the list l (v),i (i = 1, 2, … V) of the distances of the vertices, the corresponding unit vectors are collinear with these lines.

(c) the list l (r),i (i = 1, 2, … N) associated with N tangent planes, the orientations of which are uniformly random on the unit sphere.

The different definitions of the lists of length are illustrated in Supplementary Fig. S2. We remark that list (c) is closely associated with the support function of P47. It is relatively easy to see that while (a) underestimates the largest linear extension, (b) overestimates it and only (c) appears to be unbiased. For each type of list we numerically computed the α exponent providing the best fit to the experimental data and in perfect agreement with the above observation we found that in the three cases we get α ≈ 3, α ≈ 0.2 and α ≈ 1, respectively. Nevertheless, list (a) is a natural choice for the cuboid model (Model-Rect) so we used it in those computations. In case of the polyhedral model (Model-Poly) we used the (c)-list. This shows that the physically relevant value of the exponent is α = 1, while other values merely compensate for the bias in the definition of the linear extension. Hence, when fragment shapes are realistically captured the parameter α can be removed from the model.

Figure 7 presents that starting from an initially cubic body this cracking process gives rise to fragments with convex polyhedral shapes very similar to real fragments. It can be observed in Figs. 2, 3 and 4 that the polyhedral extension of the discrete stochastic fragmentation model (Model-Poly) provides also a very good description of the experimentally obtained fragment mass distributions, shape-size relations and of the statistics of different shapes and simultaneously it also reproduces the statistics of stable and unstable points of fragments (see Fig. 5). Figure 7 also illustrates that the counting of different types of equilibria of fragments can be reduced to the identification of those faces and vertices of the polyhedra which carry stable and unstable equilibrium points, respectively. Two fundamental aspects of equilibrium points show a remarkable match with mathematical predictions: the minimal value is 2 both for n S and n U , as predicted by the geometric robustness of such shapes46 with respect to truncations. The expected values for both quantities are below the corresponding values of a cuboid n S = 6, n U = 8, confirming the prediction based on stochastic geometry47.