On the packing and crushing of granular materials

This paper is a study of the dependence of the volume of voids in a granular material on the particle size distribution. It has previously been proposed that the volume of voids is proportional to the volume of the smallest particles. In a particle size distribution which is progressively becoming wider (e.g. as occurs due to crushing during the compression of sand), the smallest size of particle decreases, yet there are only ever a few of these particles out of many thousands or millions. This paper attempts to identify which particles govern the overall density of a granular material, and a new definition of the ‘smallest particles’ is proposed. These particles are shown to govern the void space in a range of simulations of spherical and non-spherical crushable particles. The theory also applies to idealised Apollonian sphere


Introduction
The packing characteristics of granular materials have been studied for over a century. It is generally accepted that the comparative density of any granular material is a function of particle shape, size range and size distribution ( Burmister, 1938 ), however predictions for the density or porosity for a granular assembly are usually empirical, and are only semi-analytical at best (e.g. Latham et al., 2002 ).
The packing of granular material has a direct influence over its engineering properties. For sands, the density determines the dilatancy and peak strength when sheared ( Bolton, 1986 ); the size and interconnectivity of the voids govern the permeability ( Hazen, 1892 ); whilst the distribution of interparticle contacts affects the induced particle stresses which govern particle breakage ( McDowell and Bolton, 1998 ). In addition to geotechnical engineers, achieving or predicting maximum packing densities is of relevance to the mining and construction industries.
Given the dependence of packing characteristics on the particle size distribution, particle crushing, which changes the distribution therefore also has a major influence on the constitutive behaviour. The importance of crushing in sands for example has long been acknowledged; when a sand is subjected to monotonic increasing stress (normal compression), it is well known that after a high enough stress is reached, particle crushing begins and causes a permanent decrease in volume (e.g. Nakata et al., 2001 ). Particle crushing also has a direct influence on the dilatancy of sand when sheared, and is associated with either a significant volumetric contraction or an increase in pore pressure if the pore fluid cannot escape.
Recent Discrete Element Method (DEM) simulations incorporating particle crushing have enabled all the main features of sand behaviour to be reproduced. This has led to new micro mechanical insights into the normal compression and shearing of sands, critical states, and yielding; and has further highlighted the influence that particle crushing and the particle size distribution (PSD) has on the observed macroscopic behaviour for sand (e.g. de Bono and McDowell, 2018a ).
The work presented here focuses on how the PSD controls the overall volume of a soil sample, with a particular focus on how crushing, and therefore a change in the size distribution, governs changes in overall volume.
A background to the theory linking the particle size distribution (in particular the smallest particles) to the current voids ratio or porosity is given in the next section. Section 3 provides a description of the numerical simulations, followed by results and a quantitative analysis in Section 4 . Section 5 provides some discussion followed by concluding remarks.

Background
When sand is subjected to increasing stress (e.g. isotropic or one-dimensional), after yielding the sample follows a unique normal compression line (NCL) in volume-stress space, known to be caused by particle crushing. Following Pestana andWhittle (1995), McDowell (2005) compression line for a sand should be linear on log e -log σ axes, where e is the voids ratio, defined as the ratio between the volume of voids ( V V ) and volume of solids ( V S ): e = V V / V S , and σ is the applied stress. The above theory was based on the assumption of a fractal particle size distribution and the kinematics of particle fracture. This was subsequently investigated using DEM ( McDowell and de Bono, 2013 ) by performing simulations of crushable sands, and led to a new compression law in which the normal compression line can be expressed as: where σ is the applied stress and b describes the size-effect on strength for the particles: where σ av is the average crushing strength for particles of diameter d . The strength here corresponds to bulk failure, i.e. catastrophic splitting into smaller fragments, and can be measured from crushing single particles between platens. The parameter b describes the rate at which the average strength increases with decreasing particle size (a common observation for most brittle materials), and is related to the distribution of imperfections/flaws in the material. The slope of the NCL is therefore a function of the particle strength size-effect. In deriving Eq. (1) , two major steps were necessary: firstly, to relate the applied stress to the current size distribution of particles; and secondly to relate the current particle size distribution to the current void space. The first step was achieved by assuming that following yielding, the strength of the smallest (strongest) particles, σ sm is proportional to the applied stress: which accounts for the hardening of the soil. If a fractal distribution develops via the continuous fracturing of the smallest grains, these grains will be in self-similar configurations as stress increases, loaded in the same manner, so the strengths must be proportional to the applied stress. Substituting Eq. (1) into Eq. (2) then links the applied stress to size of the smallest particles d sm : The second step-relating the particle size distribution to the current void space will now be revisited, and is the subject of a more detailed analysis. A fractal particle size distribution is usually defined as one in which (e.g. Turcotte, 1986 ): where N ( l > d ) is the number of particles with diameter greater than a 'cut-off' particle size d , and D is the fractal dimension. For a discrete, hierarchical size distribution, in which there is a constant ratio between consecutive particle sizes, the number of particles of a given size d can be expressed as ( McDowell, 2005;Palmer and Sanderson, 1991 ): The collective volume of particles of size d can then be expressed as: and the cumulative volume of the smallest particles (size d sm ) can therefore be written as: A fundamental assumption used by McDowell (2005) was that the volume of the smallest particles, V sm is proportional to the available void space ( V v ): and so: For the confined compression of a granular material, the volume of solids ( V S ) remains constant, meaning the volume of voids is proportional to the voids ratio: V V ∝ e , and therefore: Substituting Eq. (10) into Eq. (3) , and taking the typical fractal dimension of 2.5 (e.g. Turcotte, 1986 ) therefore gives the compression law in Eq. (1) , which is able to correctly predict the NCLs for a range of real sands ( de Bono and McDowell, 2018b ). The assumption in Eq. (8) however has been difficult to quantitatively confirm in any simulations due to the nature of crushing. For crushed sand, there is invariably a numerically insignificant number of the smallest particles, irrespective of their actual size. For example, there are typically only 2 fragments of the smallest physical size in samples of > 30,0 0 0 particles. The cumulative volume of these 2 particles clearly provides no information about the overall pore space or its interconnectivity.
The rationale behind the assumption in Eq. (8) was that as the smallest particles become smaller and fill the available voids, the available voids become smaller, and the typical void space must be more closely related to the volume of one of those particles as opposed to the larger particle(s). Incidentally, it has been shown previously that for space-filling sphere packings (such as Apollonian packing), the porosity of an assembly of spheres can be expressed as ( Anishchik and Medvedev, 1995;Manna and Herrmann, 1991 ): where n ( l > d ) is the porosity of the packing calculated only considering particles greater than size d, derived assuming that the volume of voids tends to zero as the smallest particle size tends to zero: The porosity n of a material is defined as the ratio between the volume of voids and the total volume, i.e. n = V V / V T . Putting the minimum particle size into Eq. (11) therefore gives the overall porosity of the packing: Apollonian sphere packing is the process of iteratively filling the voids between existing spheres with the largest possible mutuallytangent sphere (e.g. Anishchik and Medvedev, 1995 ). The insertion of new spheres means the total enclosing volume V T remains constant whilst the cumulative volume of particles V S increases. Therefore the porosity is proportional to the volume of voids: n ∝ V V , which substituting into Eq. (14) gives Eq. (10) , providing an analytical basis to the simple assumption invoked in the compression law.
In a previous attempt to quantify the volume of the smallest particles, from simulations of crushable spheres it was found that the most suitable method was using the coordination number, defining the so-called smallest particles as those with minimal contacts (as opposed to using some finite physical measure). This approach seemed logical as it was capable of tracking a continuously and gradually evolving PSD, whilst also taking into account any local variations in the PSD through the sample. Thus, defining the 'smallest particles' as those (spheres) with 5 or fewer contacts, there emerged a direct proportionality between the cumulative volume of these smallest particles, V sm and the overall voids ratio e , for all states on the normal compression line, and was interpreted as confirming the assumption that the volume of the smallest particles is proportional to the void space. It was also found that this relation held true for all states on the critical state line, which justified why the critical state line is parallel to the NCL ( de Bono and McDowell, 2018a ).
In this previous analysis, which used spheres, it was speculated that the value of 5 contacts was unimportant and assumed to be a function of particle shape. Although ostensibly arbitrary, the value of 5 was taken based on the fact that higher coordination numbers are only possible if a particle is surrounded by smaller particles. For a dense random packing of spheres, the minimum number of contacts a load-carrying particle may have is 4. The only feasible way for a sphere to be in contact with a much greater number of surrounding particles is if they are much smaller. Although theoretically the greatest number of contacts possible with equallysized neighbouring spheres is 12 (strictly for hexagonal close packing), for real, random packings this number will be much lower.
This work seeks to define a more rigorous definition of the 'smallest particles', avoiding the use of an arbitrary value of contacts by investigating non-spherical particles, and further explore what governs the overall density of a granular material. New isotropic compression simulations are presented for crushable sands using a range of particle shapes, followed by a brief analysis of sphere-packing.

DEM model
Isotropic compression simulations were performed on crushable samples with varying particle shapes, for the purpose of achieving a decrease in overall volume caused by a changing particle size distribution. The principles and modelling procedure behind all simulations are identical to previous work by the authors, the only difference being the use of non-spherical particle shapes. Modelling details such as the contact model, simulation procedure, breakage mechanism and criterion, and all other mechanical properties are given in several prior publications (e.g. de Bono and McDowell, 2018c ). These details have no effect on the present analysis, which is only concerned with volumetrics, so will be omitted for brevity; only relevant details and the new particle shapes will be discussed.
All simulations were performed using PFC ( Itasca, 2014 ). Particles were modelled using 'clumps'-comprised of overlapping 'subspheres'. All particles are crushable, and isotropic compression simulations were performed on each sample in order to generate a normal compression line and fractal particle size distribution.
A total of 5 different shapes were used including spheres. Three simple clump shapes consisting of 2 or 3 sub-spheres are illustrated in Fig. 1 . Particles were allowed to undergo bulk failure, whereby breaking particles were replaced by 2 smaller fragments of the same shape, obeying conservation of mass (and temporarily overlapping). The final clump simulation uses a combination of more complex particle shapes (the motivation behind this case was to minimise fragment overlap when placing new fragments, however this is outside the current scope and the details are irrelevant to the present analysis). This case involves a cycle of 3 different particle shapes, whereby a particle splits into two fragments of the next shape in the cycle (meaning each shape has a different aspect ratio), which repeats itself, shown in Fig. 2 .
The compression simulations were performed on cylindrical triaxial samples, illustrated in Fig. 3 . These were enclosed vertically by rigid platens and laterally by a flexible, faceted cylindrical membrane, which were used to apply the axial and radial stresses. The majority of the initial samples (before isotropic compression) consisted of approximately 800 randomly packed particles, with spherical equivalent diameter d e = 2 mm (diameter of a sphere of equal volume), shown in Fig. 3 (a). Although this initial number of particles is small, this was to enable simulations to be performed in an acceptable timeframe. The crushing that occurs in all simulations results in a great quantity of (mainly fine) particles, e.g. Fig. 3 (b). An important feature of all simulations is that there was no comminution limit-i.e. no lower limit to particle size. This allows the unadulterated evolution of realistic fractal particle size distributions. A consequence of this however is that the numerical timestep used in the simulations (influenced by d min ) becomes very small ( Itasca, 2014 ). Thus, to enable the simulations to be performed whilst not imposing artificial conditions, the initial number of particles was small, helping minimise the calculation time once extensive crushing has occurred.

Normal compression
The stress-strain results which are used to perform volumetric analysis are presented in Fig. 4 . Fig. 4 (a) shows the normal compression results for all clump simulations. The simulations demonstrate the familiar sand behaviour: an initial elastic response, until the onset of crushing causes yielding, followed by a more dramatic decrease in volume. All simulations demonstrate tend towards the same slope, which is to be expected as all materials have the same size-effect on average particle strengths. Fig. 4 (b) shows the normal compression results obtained previously for spheres, as well as all critical states that were obtained from shear tests on samples at various stresses ( de Bono and McDowell, 2018c ). Also shown is the normal compression results for a new simulation of spheres using a larger sample (60 0 0 initial spheres rising to > 150,0 0 0).
All initial samples have identical, uniform initial PSDs. The final PSDs from all normal compression simulations are shown in Fig. 5 in the conventional manner, showing the varying degrees of breakage that have occurred, which correlate with the decrease in volume (taking into account the different initial volumes). Similar PSDs exist at the critical states, with more crushing associated with higher stresses and lower voids ratios ( de Bono and McDowell, 2018c ). All particle size distributions tend towards a fractal distribution with a fractal dimension of 2.5, plotted in Fig. 6 . The fractal dimension is most clear in Fig. 6 (b), and once D = 2.5, remains constant as the span of fractal particle sizes broadens.

Void space as a function of the smallest particles
A similar approach to defining the 'smallest particles' based on the coordination number is used here. The average coordination number for the various clump simulations varies greatly from that for spheres. Before crushing, the average coordination number for spheres is around 6, whilst that for clumps varies between 6 and 12 depending on the shape (with greater variation once crushing begins). The prism-shaped clumps ( D ) invariably exhibit the highest coordination numbers, which should be expected due to the increased capacity for multiple contacts between two particle surfaces.
It therefore seemed logical to use some relative measure of coordination number, and it was found that a more robust definition of the smallest particles is achieved by comparing the number of contacts a particle has with smaller and bigger neighbouring particles. If Z sm is the coordination number with smaller particles, and Z bg is the coordination number with bigger particles then: if Z sm ≤ Z bg smallest particle , representative of the local void space if Z sm > Z bg many contacts , not representative of the void space (15)  The use of particle contact information to determine whether or not a particle is effectively the smallest means that these 'smallest particles' will include particles with different physical sizes, and has the advantage of taking into account local differences in the sample grading.
Calculating the cumulative volume of all such smallest particles as V sm , demonstrates a unique direct proportionality with the volume of voids ( Eq. (8) ), and therefore the voids ratio. For the simulations here, V sm is of the order 10 -6 , so the plotted values have been normalised by the volume of solids: V sm * = V sm / V s . Fig. 7 (a) plots e as a function of V sm * for the clump simulations. As soon as breakage initiates, all clump simulations follow a linear path directed towards the origin, as the volume of the smallest particles and overall void space decrease. Before any breakage the samples consist of a single size of particle, all being simultaneously the smallest and largest, and V sm * = 1, and the relationship does clearly not apply. The prism clumps ( D ) appear to have a different coefficient of proportionality, which is likely to be a function of the particle shape, although the other clump simulations appear approximately coincident. Fig. 7 (b) shows the equivalent data for the simulations using spheres; the use of V sm * facilitates a comparison between the two simulations with different sample size. All simulations appear to have the same constant of proportionality. The large samples of spheres displays almost perfect linearity, suggesting that the noise in the data calculated from other simulations could be due to the relatively small sample size.
One observation is that the relation V sm ∝ e holds for a wide range of states, including states immediately following yielding, before the PSD becomes fractal in character; V sm ∝ e appears to apply to all states except the initial states when the samples consists of a single particle size. The same constant of proportionality between the compression simulations and the critical states suggests a unique relation between V sm and e , independent of the stress path or applied stress ratio. To investigate this further, the same analysis was applied to results from a previous simulation which was used to explore the elastic-plastic stress space, also plotted in Fig. 7 (b). The stress-path of this simulation is summarised in the inset, full details showing the location of the yield surface can be found in de Bono and McDowell (2018d) . For this case, V sm ∝ e remains true even for elastic states when particle crushing is not occurring.

Application to sphere packing
To further investigate the proposed theory-that particles with mainly larger neighbours are representative of the voids-the analysis will now be applied to an Apollonian sphere packing. Apollonian sphere packing is the process of iteratively filling the voids between existing spheres with the largest possible mutuallytangent sphere (e.g. Anishchik and Medvedev, 1995 ). For this process, it has been shown that the emergent PSD is also fractal with D = 2.5, the same as observed from crushing.
An Apollonian packing was created starting with 4 initial spheres of unit radius, see Fig. 8 (a). Following 9 iterations, the final arrangement consisted of ≈50 0,0 0 0 spheres, shown in Fig.  8 (b) and (c). In this case the volume of the smallest particles, V sm , should be proportional to the porosity n rather than the voids ratio, as the volume of solids is increased with each iteration.  Despite the apparently regular process of sphere insertion, the resulting particle size distribution is far more continuous than that obtained from the compression simulations, in which particles undergo hierarchical splitting. Of the 473,039 particles, there are over 10,0 0 0 different sizes; 4 spheres with the maximum size d = 2, and 24 spheres with the minimum size d ≈ 0.001824. The ultimate PSD is shown in Fig. 9 (a), showing a clear similarity with those obtained from the crushing simulations (as well as experimentally), exhibiting a fractal dimension of 2.5 across a broad range of sizes. The calculated values of V sm are plotted against porosity in Fig. 9 (b), which is plotted from the 1st iteration onwards and annotated with the total number of spheres. The graph shows a very good correlation using this definition of the smallest particles with seemingly perfect proportionality emerging as the number of spheres and range of sizes increase. For the first 2 iterations there is clearly an inadequate number of spheres and/or range of sizes, and the relationship between n and V sm does not apply.

Discussion
The initial motivation for trying to relate the PSD to the current voids ratio was in the derivation of a compression law, able to predict the compressibility of crushable soil. The follow-up analysis here appears to reveal a more fundamental relation between the smallest particles, which are representative of the surrounding voids, and the overall packing density or voids.
The proposed definition is that these 'smallest particles' are defined as those with predominantly larger neighbouring particles. The cumulative volume of these effectively-smallest particles appears to correlate directly with the overall pore space, and appears to be valid irrespective of stress path or whether the PSD is fractal, as long as there is some range of physical particle sizes. The fact that crushing leads to a fractal distribution, allows a prediction to be made of the quantity and volume of these particles, and provides the link needed for the compression law. In a fractal distribution, at the fine end of the PSD, the geometrical configuration is always the same-for example in an Apollonian packing the smallest particles are always mutually tangent to exactly 4 surrounding larger spheres. In other words, the distribution of contacts is self-similar. However this self-similarity will be different for different stress paths-for example, samples compressed under onedimensional conditions, or critical states, are loaded anisotropically and will have a different distribution of contacts, and the 'smallest  particles' will be different in each case. For spheres, in Fig. 7 (b) the overall void space appears to be a unique function of these smallest particles; the constant of proportionality for V sm ∝ e appears the same for all loading configurations, independent of the contact distribution. This insensitivity to stress path or stress ratio is most supported by the stress-path simulation, which was subjected to positive and negative shear stresses causing plastic strains, and then isotropic compression, and during all of this V sm remains proportional to e .
Going forward, it would seem sensible to attribute a new term to these particles, such as the 'mechanically-smallest' or 'effective' particles, or some such expression, avoiding repeated and potentially confusing use of the term smallest particles. In any case, it is hoped that the proposed definition of the effectively-smallest particles, which govern the overall porosity will form the basis of further analytical solutions or predictions of granular densities. It has already been shown that if it can be predicted (or it is known) how the PSD changes, then the change in volume of these smallest particles can be inferred, which under the same loading conditions can be taken as a measure of the change in voids ratio. Knowing which particles control the overall density or porosity could have possible applications in issues of internal erosion (in which fine particles are gradually removed, changing the soils engineering properties), as well as in predicting changes in permeability due to changes in grading (e.g. caused by crushing). The permeability of any granular or soil sample depends on more factors than simply the void space, however empirical solutions typically use values such as d 10 2 ( Hazen, 1892 ) or similar-a more suitable linear dimension could be d sm .

Conclusion
This work aimed to clarify what governs the void space in granular soil. It has previously been speculated that the volume of voids is linked to the volume of the smallest particles. A new definition of 'smallest particles' has been proposed, where the smallest particles are those that have relatively fewer contacts with smaller particles compared to larger particles. The volume of such particles is proportional to the void space, and this is true for all particle shapes investigated, as well as idealised Apollonian packings of spheres.
For soil, the strength of the smallest particles gives the hardening law as stress increases. It appears that the normal compression of soil is the breakage of particles in order to achieve an efficient space-filling packing, similar to Apollonian packing and with the same fractal dimension.