A Statistical Model of Aggregates Fragmentation

A statistical model of fragmentation of aggregates is proposed, based on the stochastic propagation of cracks through the body. The propagation rules are formulated on a lattice and mimic two important features of the process -- a crack moves against the stress gradient and its energy depletes as it grows. We perform numerical simulations of the model for two-dimensional lattice and reveal that the mass distribution for small and intermediate-size fragments obeys a power-law, F(m)\propto m^(-3/2), in agreement with experimental observations. We develop an analytical theory which explains the detected power-law and demonstrate that the overall fragment mass distribution in our model agrees qualitatively with that, observed in experiments.

Introduction. Fragmentation processes are ubiquitous in nature and play an important role in many industrial processes. Numerous examples range from manufacturing comminution to collision of cosmic bodies in space. Hence, much effort has been devoted to comprehend the nature of fragmentation in different branches of research -geophysics [1,2], astrophysics [3,4], engineering [5] and material- [6] or military science [7,8].
An intriguing common property of fragmentation is a power law mass-distribution of the fragments which seems to be independent of the spatial scale or nature of the parent bodies: When heavy ions collide with a target, or asteroids suffer a high-speed impact, both produce a power law size (or charge) distribution of debris. This law has been reported in numerous experimental studies, e.g. [9][10][11][12][13][14][15][16][17].
The theoretical description of fragmentation is developing along two different lines: one, based on continuum mechanics, another -on the general statistical methods, e.g. [9,18,19]. Finite element analysis (FEM) [20], numerical simulations of agglomerate collisions with smooth particle hydrodynamics (SMH) [21], or discrete element method (DEM) [14,17] -these are the current tools to treat continuum mechanical problem of colliding material bodies. Alternatively, statistical methods, as random walk [22] or stochastic simulation of crack-tip trajectories [23] were used to model fractal crack formation and propagation. Furthermore, there were attempts to tackle the problem analytically. In particular, a semi-empirical approach to the physics of catastrophic breakup of cosmic bodies has been developed [24]; similarly, a microstructural approach has been adopted [25] to model the fracture generation and propagation in concrete.
The observed universality of the fragmentation law for the objects, drastically different in material properties and dimension, most probably implies a common physical principle inherent to all fragmentation processes [26]. Moreover, it is reasonable to assume that this principle is of statistical nature. This motivates the analysis of a model, which being very simple on the microscopic level, reflects the most prominent features of the process and adequately represents its statistical properties. In the present Letter we propose a novel statistical model for fragmentation of aggregates -macroscopic bodies, comprised of a large amount of smaller macroscopic constituents. The problem of aggregate fragmentation arises in many areas of science and technology, in particular in planetary science, where the availability of an adequate fragmentation model is crucial for understanding of the formation and evolution of planetary rings, e.g. [27][28][29]. We show that our model, based on a few very simple physical rules reproduces the main aspects of the fragmentation process and gives the power-law distribution for the debris size. Moreover, it qualitatively agrees with the experimental data. We perform numerical simulations and develop a simple analytical theory, which explains the observed power-law distribution for the fragments.
Formulation of the model. We consider aggregates, which are not too large to neglect the gravitation forces and assume that particles are kept together by the adhesive bonds. These are significantly weaker than the forces attributed to chemical bonds, hence the formation of cracks occurs only along the contact lines between the particles. We assume that all particles are of the same size and material. Then to break any adhesive contact, the amount of energy, equal to E b , is required; to destroy simultaneously n bonds one needs n-times larger energy, nE b . The quantity E b depends on the particles size, their surface tension and material parameters, e.g. [30]. To simplify further the problem, we consider its two-dimensional version and assume that the spherical constituents form a square lattice, Fig. 1.
In a collision between aggregates, which causes their fragmentation, the energy of the relative motion E coll is transformed into the surface energy of broken bonds and the kinetic energy of debris. We will ignore the latter part of the energy and explicitly consider the former one. Therefore, the problem of fragmentation in our model is, essentially, reduced to the problem of distribution of the energy E coll among the broken adhesive bonds. We assume the bonds being destroyed due the crack's propagation. Namely, we adopt following rules for the crack growth: (i) A crack originates on the surface, on the impact site, where the maximum stress is expected. It propagates (in average) away from the surface against the stress gradient, that is, "downhill" the load, and never "uphill". (ii) Each time-step the crack elongates by a single neighboring bond, consuming the energy E b , while the crack tip performs a random walk with the direction randomly chosen among two or three possible ones, as explained in Fig. 1. (iii) The propagation of a certain crack terminates if its tip arrives at the surface, or meets another crack, mimicking a crack bifurcation. When a crack terminates and the rest-energy allows for further bond-breaking, a new crack is initiated randomly on the surface of the aggregate, or bifurcates from another crack. (iv) A mechanism limiting a total crack energy is introduced. This is because crack generation and propagation needs a sufficient load gradient which decreases with distance from the impact site. Thus, we assume that if the energy assigned to a crack, chosen randomly between zero and a maximum E cross , is exhausted, the crack stops and another one is created instead. We called this "crossing factor". (v) The breaking process continues as long as the remaining impact energy, E rest = E coll − nE b (n -current number of broken bonds) is sufficient to break further bonds. It terminates whenever the total impact energy E coll is dissipated in cracks. To break a bond between particles the energy E b is needed. A crack propagates either away from the loaded surface (vertically-down on the plots), against the stress gradient, or laterally, but never along the gradient. All possible directions have equal probability.
In the case (a) there exist three allowed directions, each with the probability p3 = 1/3: one away from the surface and two lateral directions. In the cases (b) and (c) only two directions with the probability p2 = 1/2 are possible, since crack can not move back or "uphill" the load.
In Fig. 2 a typical fragmentation pattern, obtained with the use of the above fragmentation rules, is shown. Next, we analyze the statistical properties of debris.
Fragment mass distribution. We define a fragment as a collection of constituent particles connected by adhe- sive bonds where a single particle consitutes a mass unit. Hence, for a fragment of m particles we assign the mass m. We have performed numerical experiments, controlling the size of the aggregate, the total energy spent in the process and the crossing factor. To improve the fragments' statistics we perform many runs with the same set of parameters. Fig. 3 shows the fragment mass-distribution. The most part of the distribution, except for very large fragments, may be accurately fitted by the power-law Interestingly, the exponent α ≃ 1.5 is the same as in 2D egg-shell crushing experiments [14]. The data have been collected for 10 000 runs. The line represents a power law mass-distribution with an exponent ≃ −3/2. The inset shows the distribution function P (x, t) ∝ t −3/2 , which gives the probability that two particles, initially separated by distance x, meet at time t; it illustrates the affinity of the two models (see the text). Points -numerical data for x = 20, line -theory.
Analytical model. The observed power-law fragment size distribution may be obtained analytically, if we notice that the proposed fragmentation model may be mapped onto the one-dimensional model of diffusing and annihilating particles, where the particles correspond to the tips of the cracks. Indeed, the direction normal to the surface, that is, the direction of the most rapid decay of the stress, say axis y (vertical downward in Fig.  1), may be mapped onto the time axis, since the reverse motion along this axis is forbidden. Then the location of the crack tips along the lateral direction, say along axis x (horizontal in Fig. 1) corresponds to the location of fictitious "particles" on the line. One step along y axis corresponds to one time step, while one step along x axis corresponds to a "particle" jump on the line (without loss of generality we can assume unit time and space steps).
It is easy to see that the fragmentation model with the rules illustrated in Fig. 1 is equivalent to the following model: (i) Each particle on a line at each time step can remain at the same site with the probability p 0 = 1/3, move one site left or right with the probability p 1 = 1/6, two sites left or right with p 2 = 1/12, etc., k sites left or right with the probability p k = 1/(3 · 2 k ); note that p 0 + 2 k p k = 1. (ii) When two particles meet, one of them annihilates while the other one continues to move with the same rules. (iii) The fragment sizes of the primary fragmentation model is equal to the area confined by the trajectories of the particles and the axis x, Fig. 2.
First we show that the "particles" perform onedimensional diffusion motion. Indeed, due to the lack of memory, each time step does not depend on the previous one (note that the original problem does have a memory with respect to previous steps). Therefore, the mean square displacement < [∆x(t)] 2 > is a sum of these quantities for each step, ∆ 2 , where that is, [∆x(t)] 2 = 4t = 2Dt. Hence, we have a diffusive motion with the diffusion coefficient D = 2. Now we compute the distribution of the fragments' areas. These correspond to the areas of the figures on the x−t plane confined by the x-axis and the diffusive trajectories of pairs of particles, which start to move at t = 0, and meet at t > 0; the initial separation of the particles is x. The area of the figures, that is, the fragment mass scales as m ∝ tx. Hence, the mass distribution reads: where the coefficient γ accounts for the surface density and a geometric factor (1/2 for triangles), which relates a fragment area and its dimensions x and t. The coefficient κ is the line density of particles on the x axis at the starting time t = 0 (i.e. the density of crack tips which start at y = 0); the random initial distribution of particles (crack tips) implies the Poisson distribution for the initial inter-particle distance, κ exp(−κx). P (x, t) gives the probability of two diffusing particles, initially separated by distance x, to meet at time t; where one of the particles annihilates. P (x, t) may be written as is the survival probability, i.e. that both particles have not annihilated before time t, provided they started to diffuse at t = 0 separated by the distance x. It is obtained from the solution of the diffusion problem with the adsorbing boundary condition for the relative motion of two particles, with the double diffusion coefficient 2D (see e.g. [31]). Indeed the integrand in Eq. (3) gives the probability distribution of the inter-particle distance y at time t, if the initial separation was x, provided this probability is identically zero for y = 0 (the annihilation condition). Integrating this probability over all distances y > 0 gives the survival probability; differentiation of P surv (x, t) with respect to time yields the probability that the particles meet exactly at time t. Substituting P (x, t) = x exp(−x 2 /8Dt)/ √ 32Dt 3 , obtained from Eq. (3) into Eq. (2), we arrive at the expansion, where A 0 = γ 1/2 Γ(5/2)/8κ 3/2 and generally, A l = (−1) l Γ(3l + 5/2)(γ/κ) 3l+1/2 /(8κ)(4γ) 2l l! is the coefficient at m −(3/2+l) . Hence, Eq. (4) explains the power-law distribution (1) obtained numerically. The above expansion holds, however, not for very small masses, to keep the continuum diffusion approach valid, and not for very large masses, to ignore the effects of energy depletion and finite size effects of the fragmenting pattern.
Comparison with experiments. So far we considered the size distribution of fragments which mass is much smaller than the initial mass of the parent body. In this limit the numerically detected and explained powerlaw distribution coincides with the experimental one [14]. Surprisingly, our model can qualitatively and sometimes even semi-quantitatively reproduce the whole size distribution of fragments in the impact experiments. This is illustrated in Fig. 4 for the experiments of exploding eggshells [14] and in Fig. 5 for the colliding basalt spheres [10]. For the former case our model demonstrates an overall qualitatively correct behavior for F (m), while for the latter case one can achieve a quantitative agreement with the most part of the experimental fragment size distribution. It is worth noting that the agreement is obtained without any fitting or scaling of our data. Moreover, the fact that our simple two-dimensional model can describe the fragmentation statistics for three-dimensional bodies indicates, that the main ingredients of our model -the diffusive propagation of cracks against the stress gradient and the depletion of energy with the cracks' growth -are indeed the basic features of a fragmentation process.  4: (Color online). Results of our model (red circles) compared with the experimental data (green squares) of Herrmann et al. [14]. In simulations we used a 600 × 600 matrix, the total energy of E coll = 40 000 E b , the crossing-factor of Ecross = 3 000 E b and the number of runs is 200. Note the lack of any fitting or scaling of our data. . Results of our model (red circles) compared with the experimental data (green circles) of Nakamura and Fujiwara [10]. In simulations we used a 80 × 80 matrix, the total energy of E coll = 1500 E b , the crossing-factor of Ecross = 100 E b ; the number of runs is 10.
In conclusion we suggest a simple fragmentation model, based on a lattice random walk of a crack with propagation rules which mimic a real fragmentation process: A crack moves against (or laterally to) the stress gradient and the energy of the crack exhausts along its growth. We performed numerical simulations of the model and observed that the mass distribution of small and intermediate fragments obeys a power-law. The exponent of the power law, equal to −3/2, agrees with the reported experimental value. We develop an analytical theory which explains the nature of the power-law in the fragmentation statistics and demonstrate that the proposed model gives a qualitative description for the overall fragment size distribution observed in experiments.