## Abstract

During confined comminution of granular materials a power-law grain size distribution (gsd) frequently evolves. We consider this power law as a hint for fractal topology if self-similar patterns appear across the scales. We demonstrate that this ultimate topology is mostly affected by the rules that define the self-organization of the fragment subunits, which agrees well with observations from simplistic models of cellular automata. There is, however, a major difference that highlights the novelty of the current work: here the conclusion is based on a comprehensive study using two-dimensional ‘crushable’ discrete-element simulations that do not neglect physical conservation laws. Motivated by the paradigm of self-organized criticality, we further demonstrate that in uniaxial compression the emerging ultimate fractal topology, as given by the fractal dimension, is generally insensitive to alteration of global index properties of initial porosity and initial gsd. Finally, we show that the fractal dimension in the confined crushing systems is approached irrespective of alteration of the criteria that define when particles crush.

## 1. Introduction

Patterns in many complex systems arise primarily from the self-organization of their subunits (Ashby 1962; Nicolis & Prigogine 1977; Cale *et al*. 1989; Solé & Manrubia 1995). Quite often these patterns develop dynamically towards an ultimate topology that can be described by a power-law set, or a fractal set if self-similarity appears across scales (Mandelbrot 1967; Bak 1996; Scanlon *et al*. 2007). A classical example of a power-law set is the ultimate grain size distribution (gsd), frequently observed in brittle granular materials under the extreme conditions of high-pressure shear within fault gouges (Turcotte 1986). Evidence that this power-law gsd is connected to a random self-similar topology was recorded by means of image analysis of the gouge material using increasing resolutions (Sammis *et al*. 1987), suggesting that the ultimate gsd is fractal (Biegel *et al*. 1989; Einav 2007*a*). The current paper systematically explores the leading role of the self-organization of the fragments in determining the fractal nature of the ultimate gsd in comminuting brittle granular materials.

Under confinement, granular materials transmit a fascinating set of force chains through the contact network (Jaeger 1966; Drescher & de Josselin de Jong 1972; Radjai *et al*. 1996; Makse *et al*. 2000). When several forces acting on a given particle satisfy critical conditions, this particle can crush and be replaced by a set of fragments. The fact that a single crushed particle is replaced by several smaller fragments suggests that—for a given macroscopic volume—the system can choose between increasing number of plausible topological configurations. This freedom in rearrangement enables the smaller fragments to fit tightly within smaller pores. As the system compacts and the macroscopic volume decreases, newer force chains are reborn and more grains fracture until an asymptotic state is reached where the frequency of fragmentation effectively stops. This asymptotic behaviour is typified by the emergence of an ultimate gsd. In this paper, we investigate the sensitivity of the ultimate gsd to the rules that define the self-organization of the fragment subunits during increasing confinement.

In comminuting poly-dispersed granular assemblies, smaller particles normally have fewer contacts; this makes them weaker in terms of their spatial position because they need to carry larger internal shear stress concentrations. At the same time, the smaller particles could be seen as stronger in terms of their material strength: smaller particles have smaller volume and thus tend to contain fewer internal defects (Weibull 1951). This competition leads to the ‘cushioning effect’ via topological self-organization that brings bigger particles to be surrounded by smaller particles (Sammis *et al*. 1987; Tsoungui *et al*. 1999). The cushioning effect determines the asymptote of the macroscopic constitutive behaviour. Another observation at the assembly level is that the stored elastic energy in the particles statistically scales proportional to the surface area of the particles (Einav 2007*b*). So, rather than the statistics of individual crushing events, it seems that, to a first order, the constitutive behaviour of crushable granular materials is determined by the collective statistics of reorganization and energy scaling.

Towards exploring the role of self-organization, we develop a model based on the discrete-element method (DEM) (Cundall & Strack 1979), at this stage in terms of a two-dimensional idealization. The extension of the DEM to crushable media has been explored in the past using various methods (Astrom & Herrmann 1998; Tsoungui *et al*. 1999; Robertson 2000; Jensen *et al*. 2001; Lang 2002; McDowell & Harireche 2002; Cheng *et al*. 2003, 2004; Lobo-Guerrero & Vallejo 2005). The DEM development in the current paper is quite unique in that our aim is not to come up with a so-called ‘generally best’ model. We do not believe that a single model can reflect a wide range of materials and loading scenarios. For that purpose, our philosophy is aimed at identifying the physical conclusions that are beyond the details of the models. We explore various alternative model developments, and question whether under specific conditions their use leads to similar results. If we can recognize such commonalities then significant conclusions can be drawn. Obviously, an unreasonable model can lead to absurd results and defeat the purpose, which is why we carefully formulate models that follow an attractive logic.

Motivated by the paradigm of self-organized criticality (Bak *et al*. 1988), we examine whether a fractal gsd emerges at ultimate conditions and question its sensitivity to the following changes.

*Detail A1*: the external/global index properties that characterize how the distinct elements are initially packed in the sample; here explored by varying the initial global porosity and the initial gsd.*Detail A2*: the criterion that governs when an individual particle crushes; here explored using two criteria for fracture modes I and II.

Simplistic models of cellular automata suggest that the nature of global patterns depends on the organization rules of the subunits (Wolfram 1984; Wooten 2001). Motivated by this conclusion, we investigate the sensitivity of the ultimate fractal gsd to yet another detail:

*Detail B*: the rules of self-organization, here defined by replacing a crushed particle by various sets of fragments.

As mentioned above the question of sensitivity (or insensitivity) of emerging patterns in complex systems to the changes in the above details was explored in the past using simplistic models of cellular automata. However, the following paper is novel in that our current examination is based on a comprehensive physical model.

## 2. Sample preparation and critical porosities

In order to validate the insensitivity of the ultimate topology to global index properties, i.e. the initial porosity and initial gsd (detail A1), we first need to establish a framework that enables us to generate initial samples that satisfy the required global properties. ‘Initial samples’ refer to granular samples that are at the verge of changing from fluid-like to solid-like materials. This ‘jamming transition’ (Liu & Nagel 1998; Corwin *et al*. 2005) corresponds to a particular critical porosity *η*_{c} (Jaeger 1966; Aharonov & Sparks 1999; Makse *et al*. 2000; Zhang & Makse 2005). It was found that the value of *η*_{c} is affected by both the friction coefficient *μ* and the gsd. In the context of the crushing simulations, we target the full solid-like response of the granular material during compaction (i.e. diminishing porosity from the point of transition). For that purpose, our crushing simulations will always start with an initial porosity *η*_{o} that equals *η*_{c}. It is principally possible to find a range of critical porosities, for any given gsd, bounded by maximum and minimum values, and , simply by adjusting the interparticle friction coefficient *μ* during the sample preparation. We employ the two-dimensional DEM and inject about 3000 particles of prescribed gsd randomly into a space defined by a box with known dimensions; in that way, we establish a prescribed initial porosity, which may or may not represent the critical porosity. The system is set running to a point where the overall kinetic energy is essentially negligible. At the end of this equalization phase, the average number of contacts per particle, i.e. the average coordination number *C*, is calculated. This is made for various designated global porosities and uniform gsd by numbers, in terms of the different ratios of *d*_{M}/*d*_{m} between the maximum and minimum grain sizes (figure 1). As expected, the average coordination number increases with decreasing porosity, with quite an abrupt jump along a narrow range of porosities. Following Makse *et al*. (2000), we employ the value of *C*_{c}=*D*+1 to determine the critical porosity at the jamming transition *η*_{c}, where *C*_{c} is the critical average coordination number at the transition and *D* is the packing dimension (i.e. we use *C*_{c}=3 for two-dimensional simulations).

Sample preparations are executed with the friction coefficients in the range *μ*=0–0.6. The crushing simulations, which will be presented later in this paper, will always correspond to a friction coefficient of *μ*=0.6. For that purpose, after the sample preparation phase has finished, we switch off the old preparation value of *μ* and turn on the *μ*=0.6. This simple preparation technique enables us to interpret the densest initial state using the minimum preparation coefficient *μ*=0, and the loosest initial state using the maximum preparation coefficient *μ*=0.6.

To obtain an effective representation of the full range of transitional critical porosities as a function of the various gsd values, it is essential to define an index property that can reflect the effects not only of the nominal grain size but also of the shape of the distribution. For that purpose we adopt the grading index, *ϑ*, which weighs the proximity of the initial gsd to an ultimate gsd, based on energy considerations and statistical homogenization (Einav 2007*b*). For three-dimensional assemblies, this is given in terms of the second-order moments of the distributions. For two-dimensional assemblies, we use
2.1
where *J*_{10} and *J*_{1u} are the first-order moments (average grain size) of the initial and ultimate gsd, respectively. The full range of transitional critical porosities as a function of the various gsd values, via *ϑ*, is plotted in figure 2. The ultimate gsd is assumed to be a power law with a coefficient *α*=1.3. The theoretical porosity of densest random mono-sized packing (*d*_{M}/*d*_{m}=1) by Tory & Jodrey (1983) is added to the plot. The porosity of *η*=0.215 is also added, corresponding to a square packing. The maximum grain diameter ratio used in our simulations is *d*_{M}/*d*_{m}=100, for which case we obtain *ϑ*=0.428 as marked by the minimal value in figure 2.

## 3. Developing two ‘competing’ grain-fracture criteria

Based on Irwin’s (1957) classical approach to fracture mechanics, the three modes of fracture are: the *opening* mode, the *in-plane shear* mode and the *out-of-plane shear* mode. For simulating two-dimensional granular materials, only the opening and the in-plane shear modes are accessible. In the following we present two criteria, namely criterion A and criterion B. Criterion A is driven by the *opening* fracture mode, by refining the work of Tsoungui *et al*. (1999). Criterion B is motivated by an in-plane shear fracture mode, by refining the ideas of Sukumaran *et al*. (2006). A simple schematic representation of these two criteria is portrayed in figure 3. It is important to note that these two criteria are not necessarily analytically complete and cannot capture the full complexity of the single-grain crushing problem. However, our interest is in the evolving physics and patterns that are beyond the details of this condition. We question the conditions with which the predictions using these two competing models converge and link this to the leading role of self-organization (detail A2). In using competing, we highlight that these two models both follow a distinctively different logic (namely, by capturing the effect of the corresponding fracture mode).

### (a) Criterion A: diametrical splitting via fracture mode I

We develop a model that is based on a refined version of the fracture criterion by Tsoungui *et al*. (1999), which in essence represents a mode I failure of a single particle subjected to an arbitrary set of contact forces. During our DEM simulations, we crush a particle if the following ‘Brazilian’ criterion is met:
3.1
where *S*=(*s*_{1}−*s*_{2})/2≥0 represents a nominal shearing force; *N*=(*s*_{1}+*s*_{2})/2≥0 denotes the nominal normal force; while *s*_{1} and *s*_{2} are the major and minor principal forces calculated via
3.2
and
3.3
This calculation requires us to construct the effective force tensor per particle using , with *N*_{c} being the grain’s coordination number, the superscript (*c*) denoting the *c*th contact and the subscripts (*i*) and (*j*) representing the *i*th and the *j*th components of the unit normal vector and contact force vector , respectively. The representation of this force transformation is given schematically in figure 4, where *γ* designates the angle between the major principal force *s*_{1} and the horizontal axis.

Finally, to complete equation (4.1), the critical conditions of a particle of diameter *d* are bounded by the critical force magnitude *F*_{crit}, using
3.4
Following Jaeger (1966), *σ*_{fM} represents the critical tensile stress for failure of the biggest particle with a size *d*_{M}, depending on its intrinsic material elastic moduli and strength. Since grains are rarely homogeneous but rather contain internal defects we account for the variability of the particle strengths based on the number of flaws within them, using Weibull statistics given by the reduction factor *f*_{W}. The cumulative probability of a two-dimensional particle with diameter *d* to survive a stress *σ*_{f}, attained by a two-force diametrical loading test, is given by McDowell *et al*. (1996)
3.5
where *m* is the Weibull modulus. Extracting the inverse cumulative distribution function from the above equation and using *σ*_{f}=*σ*_{fM}*f*_{W}, we obtain
3.6
Note that two physical limits are satisfied: (i) none of the particles can sustain an infinite tensile stress and (ii) all the particles will survive zero stress. However, when a grain is loaded symmetrically, for example as described by the scenario in figure 3*b*, the major and minor principal forces are the same. Therefore, the LHS of equation (3.1) is negative and the corresponding condition will never be satisfied. However, experience suggests that under such loading conditions particles may in fact crush. A competing model can therefore be proposed to enable this situation to be explored by introducing a second fracture criterion that is driven by the mode II of fracture.

### (b) Criterion B: in-plane shear via fracture mode II

To devise the second fracture criterion, which will enable a particle to crush under isotropic loading (i.e. when *s*_{1}=*s*_{2}), we compare the average of the normal component of the contact forces against a new critical force magnitude ,
3.7
A grain crushes if , given by
3.8
is greater than the . This criterion is illustrated in figure 5, which shows how an initial arbitrary set of forces is replaced, for calculation purposes, by a representative set of normal forces in isotropic configuration (equally spaced). The magnitude of each of the force components in the representative set is calculated using equation (3.8). As mentioned before, the previous criterion B is unable to predict failure in isotropic configurations, but this could develop through a mode II fracture (figure 3*b*). According to the extensive numerical study by Sukumaran *et al*. (2006), particles under isotropic configuration will crush when the following critical force magnitude is met (in relation to equation (3.7)):
3.9
Note that is given as *F*_{crit} (of our previous criterion) multiplied by two additional partial factors *f*_{CN} and *f*_{D}, which represent the dependency of the isotropic loading mode on the coordination number and the curvature of the loading particles. The precise definition of these factors is given in the following.

Assuming that particles split according to the maximum tensile stress and using the elastic solution by Timoshenko & Goodier (1970) of a particle squeezed diametrically between two loading grains, Sukumaran *et al*. (2006) proposed capturing the curvature effect by the following partial factor:
3.10
where *d* is the diameter of the grain under question and *D* is the diameter of the two loading grains. Since more than two particles will normally contact a grain under question, in our simulations, we take *D* as the average of the contacting grains. If the grain under question contacts a wall we use , i.e. in which case this factor is omitted since *f*_{D}→1.

The second factor in equation (3.9), *f*_{CN}, was introduced to account for the dependency of the isotropic mode of failure on the coordination number, in accordance with the empirical relation of Sukumaran *et al*. (2006) (for any *C* greater than 1)
3.11
where *C* is the coordination number of the crushing grain. This expression suggests that, under increasing coordination number, the factor *f*_{CN} increases such that it is harder to crush the particle. Also, it satisfies the logical boundary of *f*_{CN}=1 when *C*=2, corresponding to two-force diametrical loading.

## 4. Rules of self-organization via fragment configurations

In the previous two sections, we developed the framework to explore the latter insensitivity of global patterns to changes in global properties. The current section aims to set up a framework that enables study of the sensitivity of the ultimate fractal gsd to changes in the rules that control self-organization of the subunits (i.e. detail B). Changing these rules is governed here by altering the number of fragments that replace a grain once it has been crushed. We explore three different configurations of fragment replacements (see figure 6 for configurations I–III).

An attempt to prescribe circular fragments within the circumference of the parent crushed particle, without introducing overlaps between the fragments, will always fail to satisfy mass conservation. Otherwise, if overlaps are introduced unwanted numerical localization of continuous internal crushing from artificially large forces may occur. In order to resolve this situation Astrom & Herrmann (1998) proposed fitting the rest of the mass (i.e. the mass difference between the parent particle and the non-overlapping replacing fragments) in the nearby surrounding pockets. Our strategy is different. We conserve the mass by involving two phases. During the first phase the fragments are prescribed and randomly rotated, without conserving the mass, as illustrated in figure 6. However, a rapid linear expansion is then introduced in the second phase to gain back the overall solid mass.

Our mass conservation strategy is supported by the identification of two distinct dynamic timescales during crushing: (i) *the local timescale*, through the rapid rearrangement of the fragments that occurs immediately after a particle splits and (ii) *the global timescale*, from the dynamics of the boundary conditions (here, the moving walls). The artificial expansion time was carefully selected to be smaller than the local timescale, and in any case was less than 0.002 per cent of the total simulation time. Therefore, the artificial correction of the mass could be regarded as adequate for all practical purposes.

## 5. Results

By adopting the framework as established in §§2–4, we study the leading role of the rules of self-organization (detail B) in the crushing system and its general insensitivity to minor alterations of global index properties (detail A1) and fracture model (detail A2). The crushing system is defined by subjecting the model material to uniaxial compression. The following grain properties were selected, which compare well with experiments on brittle sands: normal and shear stiffness of 1×10^{9} N m^{−1}, particle density of 2000 kg m^{−3}, the Weibull coefficient *m*=4 while *σ*_{fM}=18 MPa; the largest grain is always *d*_{M}=6 mm, while the smallest initial grain size is set to obtain the required *d*_{M}/*d*_{m} of the uniform gsd by number. The following figures show results that combine data from two similar simulations that are different only with respect to their initial number of particles. The larger sample starts with more than 10 000 particles to enable smooth results to be obtained at the start of the loading. Once the number of fragments exceeds the computational limit, we then use the data obtained from a smaller sample simulation (initially with about 300 particles), which at this stage acquires enough particles to continue retrieving the desired smooth results. The horizontal walls are fixed with normal and shear stiffness of 1×10^{11} N m^{−1} that constrains the lateral expansion of the entire aggregate. Axial load is increased gradually by moving the upper and lower lateral walls with a relative inwards velocity that grows asymptotically to 0.3 m s^{−1}. This slow loading velocity was chosen to prevent dynamic effects due to non-dissipating wave propagation. Pressure is given by the average of the vertical and lateral stresses, both of which are calculated by summing the exerted forces from the particles on the corresponding walls divided by the length of these walls. We note that in our simulations the fragmentation begins at around 5–10 MPa, from which point the gsd evolves rapidly. Then, the frequency of the fragmentation gradually decays, showing an asymptotic behaviour in terms of the rate at which the gsd evolves.

### (a) Typical topology

Our simulations reveal fascinating patterns that evolve dynamically during the simulation, until crushing asymptotically ceases. A typical example of the evolving patterns is presented in figure 7, starting from the initial assembly shown in figure 7*a*. This example is based on fracture criterion A and fragment configuration (I); similar patterns are observed irrespective of the chosen fracture criteria and the global initial index properties, but with a varying magnitude depending on the chosen fragment configuration. Soon after the external loading begins, particles start to split (figure 7*b*) with a propagation of crushing across the sample. This is referred to here as the initial yielding point. We reveal a new element of behaviour, which we call the ‘non-local effect of crushing’ (see the localized area of the first generation of fragments in figure 7*b*). A fractured particle removes support from its adjacent particles, increasing their chances of being crushed. In criterion A, this occurs due to the increasing ratio of *S*/*N*. In criterion B, this occurs owing to reduced coordination number, and thus decreases in *f*_{CN}. Each generation of fragments is distinguished by a different colour, i.e. if a fragment re-crushes its colour changes. This process continues without limiting the number of the new fragments or the number of generations per single particle (figure 7*c*). Finally, figure 7*d* highlights the random self-similarity that emerges at the end of simulation across at least three scales. This is significantly higher than what could be observed in previous studies.

During the simulations, the gsd by mass and the gsd by number were extracted at pressure increments of 5 MPa. We examine the evolution of these distributions towards the theoretical power-law set. Since we have established in figure 7*d* that (random) self-similar topology does in fact appear across the scales, we refer to this ultimate power law as fractal. The assumption that fractal geometry evolves during confined granular crushing is well received (e.g. Sammis *et al*. 1987; McDowell *et al*. 1996; Einav 2007*a*). It has been proposed by Turcotte (1986) that the number of particles with a diameter Δ that are greater than *d* can be expressed by
5.1
where *α* is the fractal dimension and *k* is a normalization constant with respect to the total number of particles in the assembly. In two-dimensional systems, the expression in equation (6.1) converts to a cumulative gsd by mass of *F*(*d*)=(*d*/*d*_{M})^{2−α}, if the smallest possible grain size is not limited from below, where *d*_{M} is the maximum grain size. Based on energy consideration of fracture growths Kendall (1978) suggested that it should be impossible to crush a particle below a certain threshold. Taking this physical constraint into account, the cumulative fractal gsd by mass is given as (Einav 2007*a*)
5.2

### (b) Leading role of self-organization

Figure 8 presents the evolution of the gsd in three simulations, all involving the use of the same criterion A, but adopting the various configurations I–III. Figure 8*a* presents the gsd by number on a log–log scale. An additional regression line is added to best fit the final numerical distribution. The slope of this line corresponds to the fractal coefficient *α* (namely the fractal dimension) according to equation (5.1).

All three of the simulations do end up with a power-law scaling, but with varying fractal dimensions: *α*=1.3, 1.4 and 1.16, for configurations I, II and III, respectively. The conclusion is that the chosen rule of self-organization of the subunits—as it is being controlled by the fragment allocation—does affect the ultimate topology only through a change in the power-law coefficient. The same conclusion and fractal dimensions are recovered by applying equation (5.2), this time expressing linearity between *F*(*d*) and (*d*/*d*_{M})^{2−α}, as shown in figure 8*b*.

### (c) Insensitivity to mode of fracture

To investigate the insensitivity of the ultimate gsd to the fracture criterion of a single particle, we examine two simulations based on the distinct criteria A and B. Both simulations commence with a similar initial porosity and initial uniform gsd, where the largest particle is twice the initial smallest. Both of the simulations are based on configuration A of post-crushing allocation. Figure 9 demonstrates the evolution of the gsd during these two simulations. It is intriguing that the fractal dimension of the ultimate gsd is identical and equals *α*=1.3, irrespective of the fact that the simulations are based on distinguishably different modes of fracture.

### (d) Insensitivity to global index properties

Finally, we investigate the sensitivity of the ultimate fractal gsd to alterations in the initial global index properties of gsd and porosity. For that purpose, it is useful to look at the porosity–pressure response in figure 10 in the various simulations. Figure 10*a* presents these curves for three samples with similar initial uniform gsd (*d*_{M}/*d*_{m}=2), but corresponding to different initial porosities: , *η*_{o}=*η*_{c}(*μ*=0.2)=0.195 and . Figure 10*b* presents the porosity–pressure response of two samples with similar initial porosity but different initial uniform gsd: *d*_{M}/*d*_{m}=2 and 10. All of the porosity–pressure curves, in figure 10*a,b*, evolve towards a single unique line as the pressure increases. The corresponding evolutions of the gsd by mass, during all of these simulations, are displayed in figure 11, while adding the fractal line using equation (5.2) and the fractal dimension of *α*=1.3. While we observe slight variations from this fractal line, they are minimal, which suggests that the ultimate fractal gsd is generally insensitive to changes in the initial index properties.

## 6. Conclusions

This paper presents the development of various two-dimensional crushable discrete-element models. The models follow an established strategy that involves the replacement of a pre-crushed particle by post-crushing fragments, once a certain crushing criterion is met. Simulations of confined uniaxial compression were executed using two fracture criteria (motivated by the opening fracture mode I and the in-plane shear fracture mode II) and three alternative allocations of fragment configuration. In all of our simulations—irrespective of which fracture criterion and fragment configuration were chosen, or which global index properties of initial porosity and initial gsd were assigned—the gsd approached asymptotically towards an ultimate power-law distribution. Since the topological pattern of the particles in these simulations exhibited random self-similarity, we conclude that the ultimate gsd is indeed fractal.

In particular, we showed that the chosen rule of self-organization of the subunits—as it is being controlled by the fragment allocation—does affect the ultimate topology only through a change in the power-law coefficient. On the other hand, and motivated by the paradigm of self-organized criticality (Bak *et al*. 1988), we observed that during the compression of brittle granular materials the fractal dimension was generally insensitive to minor alterations of the global index properties that characterize the way the distinct particle elements are initially packed in the sample, i.e. initial porosity and initial gsd. Finally, the fractal dimension was also shown to be insensitive to the choice between the distinguishably different modes of fracture.

## Acknowledgements

O.B.-N. gratefully acknowledges the support of a J.W. and I.C.M. Roderick Scholarship. I.E. acknowledges the Australian Research Council through Discovery grant DP0986876. Both authors would like to thank Dr Pierre Rognon and Prof. Antoinette Tordesillas for many fruitful discussions.

## Footnotes

One contribution of 17 to a Theme Issue ‘Patterns in our planet: applications of multi-scale non-equilibrium thermodynamics to Earth-system science’.

- © 2010 The Royal Society