Design of agglomerates using Weibull distribution to simulate crushable particles in the discrete element method

Abstract This paper focuses on the use of the agglomerate technique to simulate crushable particles in the Discrete Element Method. A novel approach is proposed to design a Weibullian agglomerate by mimicking flaws within the crushable particle. The particle is designed with a constant number of sub-spheres in contrast to the approach widely used in the literature. However, the adhesion bonds between sub-spheres within the particle are selected randomly from a normal distribution. The normal distribution is designed to generate negative adhesion values, which are replaced by zero adhesion to mimic flaws within the particle. It is shown that the particle designed in this fashion exhibits a tensile strength that follows the Weibull probability function. This includes the effect of particle size that is remarkably captured quantitatively. Finally, a simple method is proposed to derive the parameters of the adhesion normal distribution from the Weibull parameters determined experimentally on single particle diametral compression tests.


Introduction
Particle crushing, or particle breakage, is a critical aspect of the mechanical response of granular materials.Silica sand undergoes crushing when subjected to high isotropic and deviatoric stresses as occurs, for example, around driven piles.Particle crushing significantly affects shaft and tip resistance of the pile in the short and long term (creep) (Leung et al., 1996;Klotz & Coop, 2001;Yang et al., 2010).Carbonate sand undergoes crushing even at low stresses and this affects its compression and shear response (Coop, 1990;Coop et al., 2004;Tarantino & Hyde, 2005).Particle crushing under low isotropic and deviatoric stresses occurs, for example, below shallow footings (Dijkstra et al., 2013).
Many aspects of the behavior of crushable granular materials can be adequately investigated by the Discrete Element Method (DEM).DEM allows the observation of their response at the scale of the particles and the investigation of the relation between small scale phenomena and macroscopic behavior.
There are two main traditional numerical techniques used to simulate particle crushing with the DEM.The first one consists in replacing the original particle by smaller ones when a stress limit is reached (Ciantia et al., 2016;Hanely et al., 2015;Clearly & Sinnott, 2015).The main advantage of the replacement technique is a relatively small computational cost, since the number of particles only increases when breakage occurs.But this comes at the cost of representativeness as the method considers breakage to occur always in the same manner and in one single stage, whereas variations are observed in nature.Other important disadvantage of this approach consists in its inability of ensuring a mass conservation after the breakage process.
The second technique consists in representing particles as agglomerates, i.e. a group of sub-spheres "glued" to each other by adhesion forces.When the force between these sub-spheres exceeds the adhesion force, the contact is irreversibly broken.Thus, crushing results in fragmentation due to the separation of sub-spheres.The main advantage of the agglomerate technique is the better representativeness of crushing processes.For example, it allows fragmentation of the particle to occur before the particle eventually crushes.In addition, particle shape can be better reproduced, which is not the case in the replacement technique even when polyhedrons or super quadrics particles are used.However, computational costs may be a restraining factor to the use of this method.The 'agglomerate' technique is explored further in this paper.
The design of the agglomerate is a critical aspect in DEM modelling of crushable materials and several approaches have been proposed in the literature.The agglomerate design includes the number of sub-spheres used to build the particle,

Abstract
This paper focuses on the use of the agglomerate technique to simulate crushable particles in the Discrete Element Method.A novel approach is proposed to design a Weibullian agglomerate by mimicking flaws within the crushable particle.The particle is designed with a constant number of sub-spheres in contrast to the approach widely used in the literature.However, the adhesion bonds between sub-spheres within the particle are selected randomly from a normal distribution.The normal distribution is designed to generate negative adhesion values, which are replaced by zero adhesion to mimic flaws within the particle.It is shown that the particle designed in this fashion exhibits a tensile strength that follows the Weibull probability function.This includes the effect of particle size that is remarkably captured quantitatively.Finally, a simple method is proposed to derive the parameters of the adhesion normal distribution from the Weibull parameters determined experimentally on single particle diametral compression tests.the adhesion forces used to bond the sub-spheres, the approach used to simulate flaws within the particle, and the criteria used to validate the agglomerate design and calibrate the 'bonding' parameters.Wang & Arson (2016) and Afshar et al. (2017) designed agglomerates by assuming that bonds between sub-spheres within a given agglomerate have the same strength and validated the agglomerate design and bonding parameters against single particle crushing tests in diametral compression.The limitation of this approach is that the variability of the tensile strength observed experimentally on particles of similar size sampled from the same batch is not taken into account.
Indeed, tensile strength of equally sized particles of same provenance shows variability and, as a result, tensile strength of individual grains is acknowledged to be a stochastic variable (Brzesowsky et al., 2011).The variability is attributed to flaws within the particles, which vary randomly from one particle to another.According to Weibull (1951), Weibull's distribution generally captures satisfactorily the variability of the grain tensile strength since it is developed on the assumption that the resistance of a structure is determined by its weakest link or component (Brzesowsky et al., 2011).
To simulate the variability of particle tensile strength, flaws need to be introduced randomly within the DEM agglomerate.A popular approach consists in generating flaws by omitting a percentage of sub-spheres during the particle construction as first adopted by Robertson (2000).This approach has been pursued by several authors since the random omission of subspheres within the agglomerate allows reproducing a Weibull distribution of survival probability (McDowell & Harierche, 2002;Lim & McDowell, 2007;Bolton et al., 2008;Wang & Yan, 2013;Cheng et al., 2003;McDowell & Bolton, 1998).
While the omission of sub-spheres is successful in generating a Weibullian agglomerate, its capability to capture particle size effects is more controversial.McDowell & Harierche (2002) had to remove an excessive number of sub-spheres to reproduce quantitatively size effects observed experimentally.However, this made it difficult to get diametral fast fracture on these agglomerates during single particle crushing tests because the agglomerates were too porous and crumbled.Later on, Lim & McDowell (2007) observed that the coordination number is the key to reproduce size effect correctly.It was found that a minimum of around 500 sub-spheres was necessary to give an approximately constant coordination number and the correct size effect.However, this poses problems when testing an assembly of agglomerates due to the large number of sub-spheres required to represent correctly the agglomerate.More in general, the design of the agglomerate based on subsphere removal to simulate flaws appears to be based on trials and errors and no rigorous approaches have been proposed so far to 'design' the percentage of sub-spheres to be removed.
This paper presents an alternative approach to design Weibullian agglomerates.These are designed with the same number of sub-spheres.However, the adhesion forces between sub-spheres within an agglomerate are not taken constant but are assumed to vary according to a normal distribution.Different agglomerates are then characterised by different randomly distributed adhesion forces that, however, are sampled from the same probabilistic distribution.Two key aspects are addressed in the paper to validate this new approach, i) the capability to generate a Weibulllian distribution of tensile strength including size effects and ii) the development of an objective criterion to infer the parameters of the normal distribution from experimental data on single particle crushing test.

Background
Tensile strength of single grains can be measured indirectly by the diametral compression test, where single particles are placed between two rigid plates and subjected to uniaxial compression until breakage (McDowell & Bolton, 1998).Tensile strength is generally found to vary according to the Weibull distribution.The probability of survival of a particle with volume 0 V is represented by the following equation: ( ) where 0 σ is the stress value at which 37% of the particles with volume 0 V survives and the exponent m is the Weibull modulus, which is a measure of the variability of the tensile strength (m decreases as tensile strength variability increases).
As the crushing of grains is initiated by the propagation of pre-existent flaws under tensile stress, particle tensile strength is dependent on the number of flaws, and ultimately on its size, as probability of finding a flaw increases with size (Brzesowsky et al., 2011).Hence, small particles are stronger than the larger ones, and this effect of particle dimensions can be understood as a direct consequence of the tensile strength statistical distribution (McDowell & Bolton, 1998).Weibull (1951) accounts for particle size through the relation between the particle diameter d and a reference diameter d 0 , as stated in the following equation: Equation 2 can be re-written as follows where * 0 σ is the size-dependent tensile strength associated with 37% survival probability.Equation 2 makes implicitly two strong assumptions.The Weibull modulus m is independent of the grain size and the size effect is only controlled by the Weibull modulus m via the exponent 3/m, which can be determined in principle by testing a single grain size.These two assumptions have relatively wide experimental evidence including quartz sand (Nakata et al., 1999), Quiou calcareous sand (McDowell & Amon, 2000), Leighton Buzzard silica sand (McDowell, 2002), and Chongqing sandstone gravel (Xiao et al., 2018).

Particle design
A 0.6 mm diameter virtual particle of Dogs Bay sand (Tarantino & Hyde, 2005) was designed with the agglomerate technique.One particle consists of spheres glued to each other by adhesion forces as illustrated in Figure 1.As the particle is compressed, forces between the spheres increase; when a limit is reached, adhesion is lost to represent fragmentation or total crushing of the particle.
To introduce flaws in the particles, adhesion between spheres was assumed to have a normal distribution with initial mean of 500 kPa and standard deviation of 700 kPa.These values allow the existence of negative adhesion as shown in Figure 2.Such values were replaced by zero since negative values do not hold any physical meaning.These zero values allow mimicking flaws within the particle.A methodology is latter explained to choose the most appropriate values for the average normal strength and its standard deviation.
Two adhesion values are attributed to each contact, a normal and a tangential component.These are attributed independently, i.e. one component can be randomly assigned a value equal to zero whereas the other one can be assigned randomly a value different from zero.The contacts fall apart when both the tangential and normal adhesion are exceeded.
Normal and tangential stiffness, k n and k t respectively, between sub-spheres were defined according to the standard values calculated by YADE, a Discrete Element Method opensource code (Šmilauer et al., 2015).The normal stiffness k n is related to fictitious Young's modulus of the particles' material, E i , and the distance between the centre of the particle and the contact point, proportional to r i , according to the expression: The value assumed to E i was 100 MPa.
The tangential stiffness is determined as a given fraction of computed k n (Šmilauer et al., 2015).

Virtual diametral compression test
Diametral compression tests were simulated with YADE.The particle was placed on a plane while a platen on the other side compressed the particle.The simulation scheme is illustrated in Figure 3.An example of typical forcedisplacement curves is shown in Figure 4 for particles of 0.6 mm diameter assembled with 500 sub-spheres of 0.030 mm diameter.All the curves in Figure 4 were obtained from the same agglomerate, i.e., same shape and same number of spheres.The difference between the results happens because each one of the agglomerates have a unique bond distribution  between its spheres resulting in different resistances.It can be observed that a peak is well defined, and this peak force was taken to characterize the particle tensile strength.
The contact between the platens and the particle subspheres in contact with the basal plane and the platen were assumed to be frictionless and the basal plane and the platens were assumed to be infinitely stiff.
Two different contact models were used to reach the simulation desired behavior.The first one is related with the interaction between the platens and the particle.The "FrictMat_FrictMat_FrictPhys" contact model adopting friction as zero simulating a frictionless and cohesionless interaction.The other contact model used is related to the interactions between the bonded spheres inside the agglomerate."CohFrictMat_CohFrictMat_CohFrictPhys" was used to simulate these cohesive interactions.

Representativeness of the tensile strength sample
Since particle tensile strength is defined in probabilistic terms, it is important to define an adequate size of the sample required to represent the tensile strength population.McDowell (2001) suggested that the number of tests should be between 10% and 15% of the real population, depending on the standard deviation of the tensile strength.
The population of analyses required to characterize Weibull's probability was defined with the Monte Carlo Method, which is adequate for problems with random variables.This method was selected because of the randomness generated by the normal distribution of adhesion forces, since a new configuration of adhesion forces (sampled from the same probabilistic distribution) is generated in each test.
The methodology consisted in simulating initially 100 compression tests.Then, additional sets of 50 tests were executed.For each cumulative set of tests, the probabilistic distribution of the survival probability P s (d) of the particle of size d was drawn in Figure 5.It was observed that the survival probability stabilized for a number of tests ≥ 300, which is then assumed to be the minimum number of virtual crushing tests to be carried out to return a statistical representative sample of tensile strength.

Effect of the number of sub-spheres
The number of sub-spheres forming the particle introduces potentially a scale effect that was also investigated.Particles formed with 500, 1000, 2000, and 3000 sub-spheres were tested in diametral compression as illustrated in Figure 6.The minimum number of sub-spheres was selected based on Lim & McDowell (2007), who established a minimum of 500 sub-spheres to capture size-effects according to Weibull distribution.A number of 600 compression tests were simulated  on each particle to characterize their statistical distribution of survival probability.As shown in Figure 7, the number of sub-spheres does not essentially affect the distribution of survival probability.Only the distribution associated with 2000 sub-spheres deviates slightly from the other three distributions but there is not a trend as the number of subspheres is increased.

Effect of platen velocity
To test the effect of platen velocity, the same particle 0.6 mm diameter was compressed at the velocities of 0.01, 0.1, and 1 mm/s.The results of these tests in terms of distribution of survival probability are shown in Figure 8 and show that the platen velocity does not influence the diametral compression test results.A velocity of 0.1 mm/s was adopted for all the tests presented on the following.

Survival probability of particle of given size
Figure 9 shows the particle survival distribution data for particles of same size (0.6 mm) and formed with a different   number of sub-spheres (500, 1000 and 2000 respectively) fitted by the Weibull function given by Equation 2. The data for each particle in Figure 9 were generated by 600 virtual diametral compression tests.The excellent fitting clearly shows that the approach put forward in this paper based on random adhesion between sub-spheres is successful in generating a Weibullian survival distribution, which is consistent with experimental data from single particle crushing tests.
Figure 10 shows the evolution of the Weibull modulus m and the 37% survival tensile strength σ 0 with the number of virtual diametral tests used to generate the tensile strength data.The Weibull parameters tend to stabilize for a relatively large number of tests as discussed in a previous section.It is also worth noticing that number of sub-spheres used to generate the 0.6 mm particle (500, 1000, and 2000) has a negligible effect on the Weibull parameters.It was concluded that any of number of sub-spheres in this range could be adopted successfully to represent a Weibullian particle.This finding is used in the next section where the effect of particle size is investigated.

Particle size effect
Weibull (1951) introduced the effect of particle size in his equations to account for the higher probability of crushing of larger particles.It should therefore be demonstrated that the agglomerate design approach put forward in this paper is capable of capturing qualitatively and quantitatively the influence of the size of the compressed particle.
Studies on size effects were performed by testing particles of 0.48, 0.60 and 0.72 mm diameter constructed by using 500, 1000, and 1700 sub-spheres respectively, as shown in Figure 11.The size of the sub-spheres was maintained the same in each particle and equal to 0.025 mm.
Each particle was compressed 600 times to derive a Weibull probability distribution.The survival probability distribution for these three particles is shown in Figure 12.It can be inferred that the particle size effect is captured correctly.The distribution of the largest particle is shifted to the right, indicating that the largest particle is more susceptible to crushing.
It is interesting to plot the evolution of the Weibull modulus m with the number of tests for the three particles as shown in Figure 13.When the number of tests becomes sufficiently large, the Weibull modulus m is found to be essentially independent of the particle size.This is a feature required to the agglomerate in accordance with experimental results as discussed in Section 2.
Figure 14 shows the 37% survival tensile strength σ 0 * plotted against the particle size d.The slope of the curve log d-log σ 0 * is equal to 1.40, which equals the ratio ∼ 3/m = 1.4 as derived from Figure 13.The approach based on the normal distribution of adhesion (Figure 2) therefore proves capable of capturing quantitatively the effect of particle size according to the Weibull theory (Equations 2 and 4).

Calibration of bonding parameters
After attesting the validity of the approach proposed to design agglomerates, it is useful to discuss how the parameters of the normal distribution of the adhesion forces, namely the Average Av the Standard Deviation SD can be calibrated to match the behaviour of real particles in diametral compression observed experimentally.
Figure 15 shows the effect of standard deviation to average ratio SD/Av on the Weibull modulus m and the 37% survival tensile strength σ 0 for a given average Av.It can be observed that the ratio SD/Av controls both the parameters m and σ 0 of the Weibull distribution.Figure 16 shows the effect of average Av on the Weibull modulus m and the 37% survival tensile strength σ 0 for a given ratio SD/Av.In this case, the average Av seems to have essentially no effect on the Weibull modulus m but only on the 37% survival tensile strength σ 0 .This can be used    Let us assume in this example that the diametral compression tests on crushable grains return m = 3.2 and σ 0 = 2900 kPa.A value of Av = 500 kPa is initially tentatively selected and the dependency of m and σ 0 upon the ratio SD/Av is explored numerically as shown in Figure 15. Figure 15a shows that m = 3.2 is associated with value of SD/Av = 1, which is value selected for the normal distribution of adhesion forces.For this ratio SD/Av = 1, the dependency of m and σ 0 upon the average is explored numerically in Figure 16. Figure 16b shows that σ 0 = 2900 kPa is associated with Av = 700 kPa, which is the value eventually selected for the normal distribution of adhesion forces.
Because m and σ 0 are not controlled independently by SD/Av and Av, this new value of Av = 700 kPa returns a value of m = 3.5 as shown in Figure 16b that is slightly different from the target one initially inputted in Figure 15a.However, this difference is not significant and can be deemed acceptable for most practical applications.If not, the iteration can be started again by considering Av = 700 kPa rather than Av = 500 kPa in Figure 15.

Conclusion
This paper has presented an approach to design crushable particles in the Discrete Element Method.The particle is designed with a given number of subspheres, and the adhesion bonds between sub-spheres are selected randomly from a normal distribution.The normal distribution is designed to generate negative adhesion values, which are then replaced by zero adhesion to mimic flaws within the particle.
The effects of number of tests, platen velocity, and number of sub-spheres composing the particle have been preliminary investigated to make sure that numerical results were not biased by the numerical setting.Numerical results have shown that the tensile strength of particle of given size follows the Weibull distribution and that the dependency of the 37% survival tensile strength on particle size is also consistent with the Weibull distribution.It is worthy noticing that all tests were carried out for only one agglomerate shape, positioned in one direction.
The approach proposed in the paper is therefore capable of producing a Weibullian agglomerate in a simple and robust fashion.Finally, a simple method is proposed to derive the parameters of the adhesion normal distribution from the Weibull parameters determined experimentally on single particle diametral compression tests.

Figure 2 .
Figure 2. Normal distribution of adhesion forces between spheres.

Figure 5 .
Figure 5.Effect of number of tests on the statistical distribution of survival probability P s (d).

Figure 9 .
Figure 9. Survival statistical distribution for particle 0.6 mm size fitted with Weibull probability function.

Figure 12 .
Figure 12.Statistical distribution of survival probability for particles having 0.48, 0.60, and 0.72 mm diameter fitted with the Weibull probability function.

Figure 13 .
Figure 13.Effect of particle size on Weibull modulus m.

Figure 15 .
Figure 15.Effect of standard deviation to average ratio SD/Av on (a) Weibull modulus m and (b) 37% survival tensile strength σ 0 .