Particle breakage criteria in discrete-element modelling Open Access

Particle crushing is usually modelled using the discrete-element method (DEM) and employing one of two methods: agglomerates or replaceable particles. Using agglomerates involves representing individual soil grains by groups of smaller sub-particles that are bonded together and can fragment as/when the bonds are broken. The replacement method, favoured by the present authors, involves modelling grains with single particles and replacing them with smaller fragments once some characteristic stress within the original particle is deemed to have overcome the particle strength.

The most obvious advantage of using the replacement method is computational efficiency; the number of particles in such a model is equal to the number of soil grains modelled, with no need for large quantities of smaller sub-particles, and there is no arbitrary comminution limit imposed by the existence of elementary particles. The replacement method also avoids the problems inherent in determining the current voids ratio for an aggregate of agglomerates, which are porous and have internal voids. Although useful in a qualitative sense, agglomerates are limited in their ability to correctly and quantitatively model the evolution of voids ratio.

A prerequisite for the replacement method is a suitable breakage criterion, that is, a measure of some characteristic particle stress which can be related to experimentally obtained particle strengths. Such strengths are typically measured by crushing single particles diametrically between flat platens (McDowell & Amon, 2000; Nakata et al., 2001; McDowell, 2002). Jaeger (1967) suggested that the tensile strength of particles could be measured as

where F is the diametric force at failure and d is the particle size. Measured in this way, the average particle strength, σ_av, is usually found to be a function of size, with smaller particles exhibiting higher average strengths and therefore being statistically stronger. This is usually expressed in the form

where the parameter b is a constant and is a function of the statistical distribution of flaws in the material.

The key task then is relating available particle strength data, for the case of diametric compression, to some characteristic measure of particle stress, which may result from any number of complex loading configurations. Any suitable measure of characteristic particle stress in DEM must be easily linked to the stresses measured experimentally; in other words, for diametric compression, the characteristic stress should be proportional to F/d². Furthermore, the ideal measure of particle stress must be physically reasonable and give the correct results with regard to experimental data; that is, it should lead to the emergence of a fractal particle size distribution (PSD) during compression, and the evolution of a normal compression line (NCL) when the (logarithm of) voids ratio is plotted against the logarithm of effective stress. Following McDowell & Bolton (1998), the emergence of a fractal PSD implies that any suitable breakage regime must take into consideration the coordination number, whereby smaller particles (which have higher strengths but fewer contacts) suffer higher stresses than comparatively larger particles (lower strengths but more contacts) – otherwise, if it were simply the weakest particles that are most likely to crush, then the result would be a uniform matrix of fine particles, behaviour which is not evident in the geotechnical literature.

In their previous work, the present authors' used the octahedral shear stress, q, as the characteristic particle stress (and therefore to determine whether a particle should break), defined as

and calculated from the average principal stresses. Although the internal stresses within a loaded sphere are not uniform, and vary with position, it is generally accepted that the maximum tensile stress is proportional to F/d² (e.g. Hiramatsu & Oka, 1966; Jaeger, 1967; Chau et al., 2000). In the DEM software used, it was found that q is proportional to F/d², for a particle subject to diametric compression by forces F, so the octahedral shear stress could easily be related to the strengths according to equation (1). The software PFC3D (Itasca, 2015) returns the average stress tensor, $σij$⁠, for a particle according to

where V is the particle volume, N_c is the total number of contacts on the particle, x^(c) and x^(p) are locations of the contact and particle respectively, and F_j^(c,p) is the force acting on the particle at contact (c). Hence, for two equal and opposite loads F acting on a particle, the major principal stress is

while the other two principal stresses are zero. Therefore the octahedral shear stress q is approximately 0·9 F/d² (McDowell et al., 2013).

The octahedral shear stress also provides a convenient way of taking into account multiple contacts; it seems logical that a particle is more likely to break if it has few contacts and a large shear stress, whereas if a particle has many contacts, and is uniformly loaded with a low shear stress (but potentially large hydrostatic stress), the particle will be less likely to break. This reasoning agrees with the widely accepted notion that coordination number plays a pivotal role in the emergence of fractal PSDs from the crushing of granular materials (Turcotte, 1986; Palmer & Sanderson, 1991; Steacy & Sammis, 1991; McDowell & Daniell, 2001). The use of equation (3) means that a particle with three orthogonal pairs of equal opposing forces would not break, as q = 0. In reality, such a particle might fail, although one would expect larger forces would be required when compared to the case of two (diametric) contact forces (e.g. Tsoungui et al., 1999; Ben-Nun & Einav, 2010); however, if six alternative orthogonal and equal forces are superposed and superposed again then the particle, under a large hydrostatic stress but zero octahedral shear stress, would be unlikely to break – so the desired effect is for all intents and purposes achieved by using q.

By using the octahedral shear stress within particles to govern breakage, and an accurate size effect on strength, the authors reproduced the correct behaviour when simulating the normal compression of silica sand, with the correct NCL and realistic PSDs (McDowell & de Bono, 2013). However, there remains some uncertainty over the use of the octahedral shear stress, as it is not practical to obtain measures of particle strengths for an unlimited combination of contact forces experimentally in order to validate use of the q. Additionally, there appears to be no clear consensus in the literature as to what measure of stress should be used (a point often raised by the referees of the present authors' previously published work, so this paper is essentially a detailed study provoked by the thoughts of previous referees); similar works by other researchers employ a variety of breakage criteria. Some of these criteria will now be summarised, although the list is by no means exhaustive.

One of the earliest attempts to model particle crushing in DEM using the replacement method was by Åström & Herrmann (1998). In two dimensions, they investigated two breakage criteria, one based on the total pressure from all compressive contact forces on a particle, the other based on the largest contact force acting on a particle. The first regime, using the total pressure, led to unstable breakage that was concentrated in a single location (which was mitigated somewhat by the inclusion of gravity). However, it is difficult to gauge whether this unstable breakage was a function solely of the breakage criteria, or the lack of a size-hardening effect or their replacement mechanism. Their latter regime, using the largest contact force on a particle to govern breakage, resulted in more stable breakage, and their results suggested that an increasing number of contacts reduces the magnitude of associate forces. This breakage criterion was also later used by Couroyer et al. (2000) in three dimensions. However, a difficulty in assessing these criteria lies in the fact that no size effect on strength was present, and Couroyer et al. (2000) did not replace broken particles.

In a different approach, Tsoungui et al. (1999) calculated the principal stresses for each two-dimensional (2D) particle, meaning that arbitrary sets of contact forces could be represented simply by two pairs of opposing stresses or forces, which they assumed analogous to biaxial loading. They used finite-element analysis to compare the maximum tensile stress in a particle under diametric loading, to the maximum tensile stress from biaxial loading. They found that the presence of minor principal forces reduced the maximum tensile stress, and therefore their characteristic stress was a function of both the major and minor principal forces. Employing a hardening law of the form given in equation (2), their simulations (and therefore their breakage criterion) resulted in realistic PSDs; however, with increasing stress their material reached a state of reduced breakage, most probably due to the comminution limit they imposed on the fragments.

Later, Lobo-Guerrero & Vallejo (2005) and Lobo-Guerrero et al. (2006) elected to calculate particle stress using the maximum contact force, as F_max/d (in two dimensions), and employing a size effect on particle strength. Use of the maximum contact forces does not consider any additional contacts on a particle; however, they took particle coordination number into account by only allowing particles with three or fewer contacts to break. This imposition is somewhat artificial, however, and may prevent a true fractal distribution from emerging – if only the smallest particles may break (larger particles will have more contacts), eventually larger particles (although only very few of them) will need to fragment in order to maintain fractal proportions. Hence this obfuscates the suitability of their breakage criteria, as does the fact that their model did not obey conservation of mass when replacing broken particles, meaning it would not be capable of modelling the evolution of voids ratio. This breakage criterion was also used by Marketos & Bolton (2009) and Elghezal et al. (2013) but in three dimensions, measuring particle stress as F_max/d², with similar difficulties associated with judging the suitability of this criterion.

In another approach, Ben-Nun & Einav (2008) used the average normal contact force, F_ave, as the basis of their breakage criterion. They measured the characteristic stress in any particle as F_ave/d (in two dimensions) but then used a number of empirical multiplicative factors to account for the effects that the particle curvature, imperfections and coordination number have on the induced stress and likelihood of breakage, with an appropriate hardening law for particle strengths. While the use of F_ave does not give attention to the number of contact forces, their position or magnitude, their coordination number multiplier meant that the probability of breakage decreased with an increasing number of contacts, giving preference to smaller particles, but not preventing larger particles from breaking. Their compression simulations resulted in encouraging fractal PSDs, and notably their simulations obeyed conservation of mass and did not have an arbitrary comminution limit. Subsequently, Ben-Nun & Einav (2010) compared their breakage rule to that used previously by Tsoungui et al. (1999) (based on a measure of average shear stress), and found that both breakage rules resulted in similar fractal PSDs.

Esnault & Roux (2013) presented a DEM study which compared the use of a Von Mises criterion (using the average particle stresses) with the stress calculated from the largest contact force on a particle (F_max/d²). However, like many others, they allowed a 52% volume loss with each breakage, which, as they mentioned, neglects the mechanical role of smaller particles, something which is of greater interest herein.

More recently, several other researchers have also opted to use the maximum contact force to calculate the particle stress using F_max/d². Among them, Hanley et al. (2015) presented large-scale DEM results with a focus on the triaxial behaviour of crushable sand. The ultimate PSDs after shearing appeared realistic, suggesting that their choice of criterion could be appropriate for use in DEM, although they started with an already well-graded distribution. Ciantia et al. (2015) also used a criterion of the form F/d², while also considering varying contact area. This was following work by Russell & Muir Wood (2009) and Russell et al. (2009), who, using a modified Von Mises type criterion (Christensen, 2000) demonstrated that the maximum mobilised stress (not the maximum tensile stress) inside a sphere is a function solely of the largest compressive contact force, occurring almost directly below it, independent of all other contacts. However, in the context of leading to fractal PSDs, it is difficult to gauge whether this was an appropriate breakage rule, as, like some of the aforementioned work, their simulations involved a comminution limit and did not obey conservation of mass – which inevitably influences the evolving/resultant PSD (as well as the mechanical response).

Tapias et al. (2015) used yet another set of criteria for particle breakage, employing what can be described as a combination of agglomerates and the replacement method. As mentioned above, there are several problems associated with agglomerates when modelling large-scale problems; nonetheless, they calculated particle stress from the largest contact force, and allowed their agglomerates to fragment when one of two conditions was met. The first of these was largely similar to those outlined above, when the stress intensity factor (∝ F_max/d²) reached the fracture toughness for the particle. The second condition was when an initial particle flaw (which were distributed randomly, and limited by particle size) had grown to the size of the particle. Initial flaws were allowed to grow under any stress intensity resulting from contacts, meaning that given enough time, any loaded particle could eventually fail, although larger initial flaws and greater intensities would result in faster crack growth. However, the limitations of agglomerates outlined previously hinder any assessment of such breakage criteria.

Hence, there is no agreement in the literature as to what the characteristic stress that governs breakage should be. The focus of this paper is therefore to examine alternative criteria for crushing to ascertain whether each gives rise to (a) a realistic NCL and (b) a fractal distribution of particle sizes with the correct fractal dimension. In this way, it is hoped that further micro mechanics of normal compression can be exposed and that a suitable breakage criterion can be established, based on comparison with existing experimental data.

The simulations presented here all use a mono-disperse, cylindrical sample, initially 20 mm × 20 mm, subjected to one-dimensional normal compression by applying vertical stress increments with a resulting vertical strain. When a particle is deemed to have broken, it is replaced by two smaller spherical fragments, equal in size, which together have the same volume as the original sphere (details of the breakage criteria investigated will be given shortly). Particles split into two fragments, with the size ratio between any new fragment and its ‘parent’ particle constant. Newly placed fragments overlap, aligned in the direction of the minor principal stress axis of the breaking particle, and within the original particle boundaries. Breakage is implemented by checking all particles at once, and any particle in which the stress is equal to or greater than the strength is replaced by new fragments. The fragments are then allowed to move apart, by completing a number of time steps during which no breakages may occur. The approach is fully described by McDowell & de Bono (2013).

The first step in the sequential modelling procedure is to apply a macroscopic stress increment to the upper wall of the sample. Particles are then checked and allowed to break. Any particles that break are replaced by fragments, which are then allowed to move apart, releasing the energy induced by the artificial overlap. This continues until no further breakages occur, after which the macroscopic stress is re-applied. After achieving a macroscopic stress with no subsequent breakage, the next stress increment is applied and the simulation continues. This process continues until the simulation reaches a point where the quantity and size of the smallest particles renders the time step too small to be computationally economical. The macroscopic stress increment used is 125 kPa, and maximum velocity of the upper boundary is limited to 0·1 m/s. Model specifics are given in Table 1. For further details of the modelling procedure, including discussion of its limitations, the breakage mechanism, readers are directed to earlier publications (de Bono, 2013; McDowell & de Bono, 2013; McDowell et al., 2013).

The results consist of four simulations, each using a different breakage regime, that is, a different method of measuring particle stress and the associated strengths. The first simulation uses the octahedral shear stress, q, as the characteristic measure of stress, and is the same as that used previously by the authors (McDowell & de Bono, 2013; de Bono & McDowell, 2014, 2015). A further two simulations use the mean particle stress, p, and the major principal particle stress, σ₁, to determine particle breakage. The remaining simulation will use a stress criterion calculated from the maximum contact force, with the particle stress calculated as σ_Fm = F_max/d², referred to herein as the F_max stress. The F_max stress is one of the more commonly used criteria in the literature, although it does not relate to any characteristic stress calculated from the stress tensor for the sphere (readily retrieved in the DEM software), but rather to a critical stress within the sphere (e.g. Chau et al., 2000; Russell & Muir Wood, 2009). The F_max stress will be of particular interest due to a combination of its relatively common use and the uncertainty of the role of coordination number in such a regime. It has been suggested that using the F_max stress to govern breakage indirectly accounts for the number of contacts on a particle, whereby (larger) particles with many contacts will, on average, have smaller contact forces. However, it remains to be seen if the number of contacts on a particle has the correct effect on the maximum contact force, especially since much of the published work using the F_max stress as a criterion imposes artificial conditions relating to the coordination number.

The particle strengths used here are taken from McDowell (2002), who reported the strengths of various sizes of silica sand particles. The average crushing force and average crushing stress calculated according to equation (1) are reproduced in Table 2. McDowell found that the strengths of the silica sand particles fitted Weibull distributions. A Weibull distribution is described by a characteristic value, in this case referred to as the Weibull strength, σ₀ (a value whereby 37% of particles are stronger); and the Weibull modulus, m, which defines the variability. The Weibull strength is similar to, and proportional to the mean, whereas the modulus is directly related to the coefficient of variation. The Weibull strengths for each size of particle are also given in Table 2.

For the silica sand, the Weibull modulus m, a material property, was found to be approximately 3·3, and from Weibull's survival probability for a block of material under tension, it is possible to express the size effect on strength as

assuming that ‘bulk fracture’ dominates (if surface-initiated fracture dominates then this equation would use d^−2/m). This equation is the hardening law used in all simulations, and is used when attributing random strengths to new fragments; in other words, the Weibull strength of any particle size is determined by equation (6), and together with the Weibull modulus, defines the distribution from which random strengths are attributed. Details of the strengths and stresses for each simulation will now be outlined, and are summarised in Table 3.

This simulation uses the octahedral shear stress as the characteristic particle stress, defined by equation (3). For a sphere loaded by two flat walls in PFC3D, the value of induced octahedral shear stress is approximately equal to 0·9F/d², hence the initial 2 mm particles are given a Weibull strength of 37·5 MPa (in terms of octahedral shear stress). These 2 mm particles are given random strengths drawn from a Weibull distribution defined by characteristic value of 37·5 MPa and a modulus of 3·3. From equation (6), it can be seen that particles of the next smaller size, d = 1·59 mm, will have a characteristic strength of 46·2 MPa.

This simulation uses the mean particle stress, p, to govern breakage. This breakage criterion, like one of those investigated by Ben-Nun & Einav (2010), enables particles to crush under hydrostatic stress. The mean particle stress is calculated from the principal stresses: (σ₁ + σ₂ + σ₃)/3, which are obtained from the (average) stress tensor for the particle. From equations (4), (5) and (3), it can be deduced that for diametric compression, p is approximately equal to 0·64F/d². From the strengths in Table 2, this means that the initial particles (2 mm) are given a characteristic Weibull strength (in terms of p) of 26·5 MPa.

This simulation uses the major principal particle stress, σ₁, to govern breakage, obtained from the average particle stress tensor. From equation (5), for diametric loading, the major principal stress, σ₁, is equal to approximately 1·91F/d². Therefore the initial 2 mm particles have a characteristic Weibull strength of 79·6 MPa (in terms of σ₁).

In this simulation, the measure of stress used to govern breakage, σ_Fm, is calculated from the largest contact force, as F_max/d², where F_max is the largest normal contact force acting on the particle. Hence the Weibull strength for the 2 mm particles can be taken directly from McDowell (2002), and is 41·7 MPa.

The compression results for the four simulations are presented in Fig. 1, alongside the experimental results for the silica sand from which the strengths are taken. As discussed in the authors' previous work (McDowell & de Bono, 2013), the NCL can be described by the following law

where the slope is given by (1/2b). The parameter b describes the hardening law for the particles, and is the exponent in equation (6) (in these simulations, b = −0·91). It has been shown that this compression law correctly describes the slope of the NCL for a range of hardening laws (McDowell & de Bono, 2013). Using the strength data used here, this law predicts NCLs with slopes of approximately 0·5. An ideal line with this slope is included in Fig. 1, and it can be seen that the q- and F-simulations (as well as the experimental results) demonstrate the best agreement with this compression law. The σ₁-simulation appears to demonstrate some agreement, but at high stresses the slope of the NCL beings to change, whereas the p-simulation does not display a linear NCL. These observations are reflected in the R² values, obtained individually for each NCL and best-fit trend lines: the q-simulation reveals an R² value of 0·92 (the trend line shown), the F-simulation reveals 0·93, the σ₁-simulation 0·85 and the p-simulation gives 0·11. The simulations are terminated when the computational time step becomes too small, which although it occurs at different stresses in the four simulations, is invariably when there is a large number of particles covering a very wide range of sizes.

The most prominent observation from Fig. 1 is that the four simulations exhibit different yield stresses. The p-simulation yields first, at around a vertical stress of 5 MPa, whereas the q-simulation yields last at approximately σ_v = 11 MPa. The yield stresses for the four simulations do not correlate with the initial particle strengths, owing to the different measures of stress increasing at different rates. For each simulation, the average (mean) characteristic stress for all particles is calculated and plotted against applied vertical stress in Fig. 2, for when the stage before yielding occurs (when the four samples are physically similar, and consist of the same quantity of solely 2 mm particles).

Figure 3(a) displays the voids ratios plotted against the average particle stress normalised by initial characteristic strength (Table 2) for each simulation. In this case the yield points for all simulations coincide at a normalised particle stress of approximately 0·3. This echoes McDowell & Humphreys (2002), who observed yielding at applied stresses of approximately 25% of the characteristic particle strength for different materials, due to strong force chains forming through approximately one quarter (the proportionality depends on particle angularity (Nakata et al., 2001)). Yielding signifies the onset of particle breakage. Figure 3(b) shows the total number of particles in each simulation plotted against the normalised particle stress: major crushing begins at the same point in each simulation. The average coordination number at yield is approximately the same for all simulations (≈ 5), although there is little scope for variation due to the fact that all samples initially consist of mono-disperse spheres in the densest possible packing.

The compression law in equation (7) was based on the assumption that a fractal PSD with a fractal dimension of approximately 2·5 emerges during compression. This was based on numerous historical experimental observations (Turcotte, 1986; McDowell & Daniell, 2001). Progressive PSDs for the four materials are plotted in Fig. 4 at 1 MPa intervals, although some are coincident. For the q-simulation, the PSD broadens and results in a smooth PSD at the final stress of 45 MPa. Although not immediately clear from this style of plot, the grading curve does become fractal in character. The F-simulation displays similar behaviour, and exhibits realistic PSDs with increasing stress. The remaining p- and σ₁-simulations, however, display an increased degree of crushing of the largest particles, and the PSDs are not fractal in character.

To assess the fractal nature and obtain the fractal dimension, PSDs are often plotted on logarithmic axes (Turcotte, 1986). Because the PSDs are discrete, they can be plotted in terms of the number of particles of each particular size, given in Fig. 5. In this plot, the fractal dimension should emerge as the slope. These graphs show that the q- and F-simulations result in realistic crushed PSDs; for both simulations the final fractal dimension is approximately 2·5, and the PSDs appear linear across a range of sizes. This agrees with experimental data for crushed granular materials (Turcotte, 1986; Palmer & Sanderson, 1991), and therefore suggests that both these measures of stress are appropriate to govern breakage.

The PSD for the p-simulation appears skewed on the log-log scale in Fig. 5(b), with no distinguished linearity (although a best-fit trend line is given to indicate the gradient). The PSD in this simulation can be categorised as uniform, with the majority of particles lying within a narrow size range compared to the q- and F-simulations. For the σ₁-simulation, the grading curve is not quite linear, and in any case the slope is too steep, with a gradient in excess of 3, beyond the usual limit for a fractal dimension (Turcotte, 1986; Sammis et al., 1987; Palmer & Sanderson, 1991; McDowell & Daniell, 2001).

The evolution of the fractal dimensions obtained from the q- and F-simulations is plotted in Fig. 6. The values of the fractal dimension, D, appear to approach a constant value, despite extensive crushing continuing to occur once a value of approximately between 2·0 and 2·5 has been reached.

In response to one of the referee's interest in work and new surface created, there now follows an analysis of this issue. From the definition of a fractal, it can be stated that the number of particles N (of size L), equal or greater than a size d, is proportional to d^−D. For the q- and F-simulations, the fractal dimension D = 2·5, as shown earlier in Fig. 5. Considering the surface area of any particle is proportional to d², then the total surface area of particles equal to or larger than size d can be expressed as

Hence, the total surface area of all particles can be found by using d_s in this equation

Considering that the current vertical stress is proportional to the strength of the smallest particles (McDowell & Bolton, 1998; McDowell & de Bono, 2013), and recalling the hardening law in equation (6), it is possible to relate the stress to smallest particle size, d_s as σ ∝ d_s^−3/3·3. If this is taken here as approximately σ ∝ d_s⁻¹ (with b = 1), combining this with the surface area then leads to a total particle surface area

The new surface area resulting from crushing is often linked to the plastic work performed on the system during compression, for example, McDowell & Bolton (1998). Russell (2011) also related the surface area to the energy put into the system, which led to alternative compression laws to that given in equation (7).

The plastic work done on the sample (neglecting elastic strains) is

For the q- and F-simulations, which obey the compression law

according to equation (7) with b = 1, differentiating equation (12) gives

Separation of the variables and multiplication by σ gives

The total work done, W, during compression is then given by

Upon integration, substitution of the surface area, S_T, using equation (10) into the integral σ^0·5 leads to

that is, the plastic work is proportional to the increase in surface area. This relation is plotted in Fig. 7 for the two simulations. The work plotted (actually work done per unit volume of solids) is cumulative, calculated for each increment as σde (using the increment average stress). The new surface areas created in both simulations appear to agree quite well with equation (16).

The stress–strain behaviour and resulting PSDs show that both the octahedral shear stress, q, and F_max stress, σ_Fm, appear suitable for governing particle breakage; whereas the mean particle stress, p, and the major principal particle stress, σ₁, appear unsuitable. These latter two regimes result in compression lines that do not conform well with the compression law and unrealistic PSDs.

During normal compression, the decrease in voids ratio is facilitated by the continual crushing of the smallest particles, with new, smaller fragments filling the voids. Simultaneously, for a fractal distribution to emerge and remain consistent, the largest particles must remain mostly intact; only a very small number of them must occasionally break to maintain the fractal dimension (e.g. assuming a constant D = 2·5, one particle of size d = 2 mm should break for approximately every 3000 particles of d = 0·25 mm that break).

Smaller particles are statistically stronger than larger particles, hence they require larger stresses to break. In all cases, the calculated particle stresses involve the term d⁻², so it seems reasonable that a typical contact force will result in larger stresses in smaller particles, which is conducive to breakage of the smallest particles (although clearly the quantity and magnitude of contact forces will also exert influence on the particle stress). However, at high values of σ_v, there must be some factor which protects the larger particles, mitigating the induced particle stress and preventing substantial quantities of the largest (and weakest) particles from breaking.

It was anticipated that, for the p-simulation, the mean particle stress would not be appropriate to govern breakage, as the many contact forces acting on large particles will contribute to large principal stresses (equation (4)), hence large mean particle stress, and therefore it would consistently be the largest particles breaking, resulting in a persistently uniform PSD and hindering any decrease in voids ratio. This indeed appears to be the case, with nearly all of the original 2 mm particles breaking immediately following yield, followed by nearly all of the next biggest size (1·59 mm) breaking. At the final stress of 17 MPa, more than 90% of the mass is finer than 1·26 mm. As shown in Fig. 1, this simulation displays reduced compressibility. In other words, if particles break under hydrostatic stress, large particles continue to break, giving a comminution limit type effect and because the PSD tends to a uniform mush of decreasing sized fines, the voids ratio will tend towards being constant; this is exactly evident in Fig. 1.

The behaviour of the σ₁-simulation, interestingly, appears to be between that of the q- and p-simulations. Like the p-simulation, the largest particles, with many contact forces acting on them, will be subject to large σ₁. Hence the use of σ₁ to govern breakage means the largest particles will be very susceptible to breakage. However, in comparison with the p-simulation, this effect is not as pronounced, owing to the lack of consideration of σ₂ and σ₃; smaller particles, with few contacts, will also have large σ₁, even though the other principal stresses may be low. This fits with the observed intermediate behaviour, with evident crushing of the largest particles (but not as substantial), yet a slightly greater portion of smaller particles when compared to the p-simulation.

The F-simulation, by contrast, demonstrates realistic behaviour and results that are similar to the q-simulation. This at first seems surprising, as the use of the single largest contact force (σ_Fm = F_max/d²) is indiscriminate to the number of contacts acting on a particle, and would not appear to minimise the crushing of the largest, weakest particles. This behaviour, and the similarity with the q-simulation, can be understood by considering how comparatively large and small particles are loaded.

Small particles have few contacts and will generally not be loaded isotropically, so if F_max is large, then both q and F_max/d² will be large, owing to the absence of additional (lateral) contacts. In both simulations the smallest particles will have high stresses and continually break. The factor preventing excessive breakage of large particles in the F-simulation can be attributed to the contact forces required to break particles. In this case the particle stress is a function of the single largest contact force only, so the size effect on strength can be rewritten in terms of force as F_m,0 ∝ d^1·1. It is evident that larger particles require greater critical forces to break. Hence, statistically, a large particle will break if it comes into contact with an equal- or larger-sized particle which can sustain the contact force. Statistically, if a large particle is in contact with a much smaller one, the mutual maximum contact force will cause the smaller particle to break. This implies that once there are few enough of the largest and dispersed particles remaining, they would not be susceptible to breakage. The PSD for this simulation in Fig. 5(d) shows limited evidence of this: there is a larger proportion of 2 mm particles than expected, and after further crushing of the smaller sizes this would dislocate the fractal distribution.

To illustrate these phenomena, one can consider an idealised loaded particle, such as in Fig. 8. Fig. 8(a) shows a single particle, size d = 1 m, subjected to unit forces. For diametric loading, the different measures of stress within this sphere are: q = 0·9, σ₁ = 1·91, p = 0·64 and σ_Fm = 1 (units of Pa), as labelled. If this sphere is now subjected to additional unit forces, but in an isotropic manner (i.e. three orthogonal pairs), as demonstrated in Fig. 8(b), then the shear stress q decreases, while there is no change in the major principal stress, σ₁, or the F_max stress, σ_Fm. The mean stress, p, however, increases. If this particle is subjected to an additional set of isotropic, orthogonal forces, offset by some rotation, shown in Fig. 8(c), then both the octahedral shear stress, q, and the F_max stress, σ_Fm, do not change, while the major principal stress, σ₁, and mean stress, p, both double in magnitude. This is similar to what happens in the simulations – where the largest particles accumulate contacts as crushing progresses. However, as noted, in such cases the contact forces will generally be smaller in magnitude; such a case is shown in Fig. 8(d), which shows an identical configuration to Fig. 8(c), but with smaller forces. Again the octahedral shear stress will be zero. The F_max stress is also smaller, but σ₁ and p are still relatively large compared to the case in Fig. 8(a), despite smaller contact forces. Although highly idealised, this serves to show that, for a given size of particle, those with fewest contacts will, in general, be under the largest shear stress, q, and the largest F_max stress, σ_Fm, whereas those with the most contacts, even if the contact forces are smaller, will exhibit the largest major principal and mean stresses. Meanwhile, the smallest particles will always be loaded with the fewest contacts, such as in Fig. 8(e), which shows a smaller particle loaded diametrically. In this case, all measures of stress are larger than an equivalently loaded larger sphere in Fig. 8(a), despite smaller forces. However in comparison to Fig. 8(d), which represents the state of larger particles, the small particle in Fig. 8(e) has much larger q, σ_Fm and σ₁, but a smaller value of p. Finally, it is also worth noting that in the q- and F-simulations, for a large particle to break, there would need to be either anisotropic or relatively large contact forces acting on it, respectively. These scenarios are illustrated in Fig. 8(f), which is similar to Fig. 8(d) apart from one pair of opposing contact forces being larger in magnitude. Compared to Fig. 8(d), in this case, q and σ_Fm (and σ₁) have increased significantly, whereas p has increased by a much lesser extent. Hence, there is a clear similarity in the q- and F-simulations, whereby the smallest particles (regardless of their actual size), are consistently subjected to the fewest contacts but largest stresses.

Figure 9 shows how the average coordination number increases with particle size, the data corresponding to points where the total number of particles is approximately 25 000 in all simulations. The same trend is seen for all simulations. The smallest particles have average coordination numbers less than 1; this is because a portion of them are not in contact with other particles, while those that are carrying load have between three and five contacts. The four sets of data appear almost identical, despite quite different gradings; the q- and F-simulations at this point consist of fractal PSDs, in contrast to the p-simulation, which comprises a narrower distribution of particle sizes (so although the average coordination number with respect to particle size is similar in all simulations, the frequency distributions of coordination numbers are different). The observed trends are the same at any point during the simulations, that is, regardless of the range of sizes of particles or their quantity, the smallest particles always have the same number of contacts, whereas the largest particles gain more as crushing progresses.

The average characteristic particle stress as a function of particle size is plotted in Fig. 10; the data correspond to the final vertical stress (which is different for each simulation). The average particle stresses in this instance are calculated only considering load-carrying particles (particles with 0–1 contacts are neglected, if they were not, then the average stresses for the smallest sizes of particles would be lower). The data are plotted on logarithmic axes and, for comparison, the average particle strengths are also shown. These are the actual particle strengths, which differ slightly from the theoretical strengths according to equation (6) (which would dictate a slope of − 0·91), owing to the fact that particles of a particular size have a distribution of strengths, and the weaker particles crush first. Nonetheless, although different in magnitudes, the average particle strengths follow the same relationship with size in all simulations. For the q- and F-simulations, the average particle stress increases steadily with reducing particle size, and the average stress curve is approximately parallel with the average strengths. The data for p-simulation, however, display no linear trend between average stress and particle size, and although the very smallest particles exhibit slightly higher average stresses, the values for most sizes are similar. The difference between the average stress and strength in this (p) simulation is narrowest for the largest 2 mm particles, and in general the average particle stress increases with particle size relative to its strength.

The σ₁-simulation displays intermediate behaviour; smaller particles are under larger stresses than the larger particles, with average stress increasing with reducing size, but the rate of increase is less than the q- and F-simulations. As such, a disproportionate quantity of the large particles break as the overall stress increases.

The average stresses for individual particle sizes are plotted against the applied stress in Fig. 11. Only the first six sizes of particle are shown for clarity. In all simulations, the particle stresses increase with applied stress; however, there are markedly different patterns of behaviour between the q- and F-simulations and the p- and σ₁-simulations. In the former, the average stress in the largest, 2 mm particles increases approximately linearly, until smaller particles come into existence, at which point the rate of increase reduces, with little change beyond a vertical stress of 20 MPa. The average stress of the next-smallest size of particle, d = 1·59 mm, at first increases, to a value higher than that of the 2 mm particles, but then the rate of increase also reduces; this pattern continues with subsequently smaller sizes bearing increasingly larger stresses. It appears that the average stress for any size of particle stops increasing rapidly once there are many smaller particles in existence, at which point these particles will become surrounded and have additional contacts. For the q-simulation, this will be due to the additional contacts mitigating the induced octahedral shear stress. For the F-simulation, this will be due to the scarcity of contacts with similar or larger-sized particles, which are able to bear large contact forces.

In contrast, in the p- and σ₁-simulations, it appears that the average stress for all particle sizes increases linearly with increasing vertical stress. In the p-simulation, the six particle sizes shown exhibit similar values of average stress, which will be of most significance for the largest and weakest particles, and this is consistent with the PSDs shown in Fig. 4(b), which show very few of the largest particles remaining after compression. The σ₁-simulation shows similar behaviour, with the average stress in the largest particles increasing rapidly with increasing vertical stress; however, the smallest particles are under moderately greater stresses. This is consistent with the notion of intermediate behaviour between the q- and p-simulations.

These graphs, as well as those in the previous two figures, demonstrate the similarity between the q- and F-simulations, where the average stress in the smallest particles (whatever their actual size) increases approximately linearly with increasing σ_v, whereas the average stresses in the larger particles are mitigated and tend to approximately steady values. So in both cases, the smallest particles, always with the fewest number of contacts, are consistently subjected to the largest stresses, causing many of them to break, creating new ‘smallest’ particles, which then start ‘protecting’ the previously smallest particles. This process is continuous, and is the same regardless of scale.

This work has sought to clarify the suitability of various breakage criteria when modelling particle breakage in DEM, which was achieved by investigating the use of four different common measures of stress to determine whether or not a particle should break. In each case, the strengths used were from experimental single particle crushing tests on a silica sand.

The simulation using the octahedral shear stress resulted in the correct macroscopic stress–strain behaviour, as well as a realistic fractal PSD, which was consistent with all previous work by the authors, which used the same breakage regime. Notably, the simulation using the F_max stress: F_max/d², also produced the correct macroscopic behaviour. Although this is one of the more common breakage criteria used in DEM by other researchers, it has often been chosen arbitrarily with no justification, or owing to the theory that the maximum critical stress within a solid sphere is a function only of the maximum contact force (e.g. Chau et al., 2000; Russell et al., 2009; Ciantia et al., 2015). However, it was unclear whether this criterion would result in the correct behaviour in DEM. Both of these two breakage criteria were shown to result in continuous breakage of the smallest particles with increasing applied stress, while minimising breakage of the largest and weakest particles, leading to the correct NCL and fractal PSDs.

For the simulation that used the octahedral shear stress, q, the additional contacts acting on comparatively larger particles serve to mitigate the induced stress, thus enabling these weaker particles to survive when the applied vertical stress is high. This effect means that it is the particles with fewest contacts that will be most likely to be subjected to high shear stresses. This in turn means that particles in contact with those much larger will be subjected to high shear stresses, as large neighbouring particles physically prevent additional particles from forming contacts.

For the simulation that considered the particle stress calculated from the maximum contact force (F_max/d²), a similar effect was evident. Once comparatively large particles were in contact with many surrounding smaller particles, the only way in which these particles would break is if they were in contact with equal- or larger-sized particles. Hence, the use of both of these alternative measures of stress to determine breakage results in a regime whereby contacts between same- or larger-sized particles is minimised, and both lead to the emergence of a fractal distribution with a dimension of 2·5, consistent with experimental results.

Of the other two criteria investigated: the mean particle stress, p, and the major principal stress, σ₁, both resulted in poor modelling of the macroscopic behaviour in comparison to experimental results. The use of mean particle stress, p, meant that the largest particles, which have many contacts (all of which contribute to the mean stress), will always be the most likely to break, and as such, it is primarily the largest particles that are continuously breaking; this resulted in a uniform PSD and the evolution of a constant voids ratio. The use of σ₁ in governing breakage resulted in somewhat intermediate behaviour. In this case the largest particles, with many contact forces contributing to large principal stresses, will be under large values of σ₁, but so too will many of the smallest particles, which have fewer contacts and therefore low values of σ₂ and σ₃.

This paper has therefore served as a detailed investigation into the factors affecting particle breakage during normal compression and has provided insight into the evolution of fractal distributions of particle sizes which accompany a NCL.

The authors are grateful to the Engineering and Physical Sciences Research Council for their financial support through Research Grant EP/L019779/1.

b

size hardening parameter

D

fractal dimension

d

particle size (diameter)

d_s

smallest particle size

e

voids ratio

F

force

F_max

largest contact force

G

particle shear modulus

L

length

m

Weibull modulus

N

number of particles

N_c

total number of contacts on particle

p

mean stress

q

octahedral shear stress

S

surface area

S_T

total particle surface area

V

particle volume

W

work

x

position vector

ν

Poisson ratio

σ

stress

σ_av

average particle strength

σ_Fm = F_max/d²

maximum contact force stress

$σij$

particle stress tensor

σ_v

vertical stress

σ₀

characteristic Weibull strength

σ₁, σ₂, σ₃

principal stresses