DEM Analysis of Single-Particle Crushing Considering the Inhomogeneity of Material Properties

Crushing characteristics of single particles are the basis of granular material simulation with discrete element method (DEM). To improve the universality and precision of crushable DEM model, inhomogeneous stiffness and strength properties are introduced into the bonded particle method, with which the Weibull distribution and size effect of particle strength can be reproduced without deleting elementary balls. The issues of particle strength and carrying capacity under complex contact conditions are investigated in this work by symmetric loading tests, asymmetric loading tests, and ball–ball loading tests. Results of numerical experiments indicate that particle carrying capacity is significantly influenced by coordination numbers, the symmetry of contact points, as well as the relative size of its neighbors. Contact conditions also show impact on single-particle crushing categories and the origin position of inner particle cracks. The existing stress indexes and assumptions of particle crushing criterion are proved to be inappropriate for general loading cases. Both the inherent inhomogeneity and contact conditions of particles should be taken into consideration in the simulation of granular materials.


Introduction
Particle crushing is a key issue in the research of granular materials in civil engineering, mineral industry, and chemical industry. Large amounts of experiments have been carried out to investigate the crushing laws of granular materials since last century [1][2][3][4][5][6][7][8]. In recent decades, with the development of experimental science, some new technologies, such as X-ray technology and tomography technology, have been introduced into particle crushing studies [9][10][11]. With the aid of these technologies, the inside situation of granular specimens becomes visible in laboratory tests, while the evolution of inner flaws remains difficult to be traced clearly for the limitation of resolution ratio of X-ray and CT graphics. Besides, considering the limitation of equipment and research expenditure, some complex loading conditions are still unable to be achieved in laboratory. Hence, numerical experiments, in which the crushing process can be easily observed in multiscale views, become an important and effective supplemental means of particle crushing research.
In the agglomerate method, a cluster of elementary balls are stuck together to form a macrobreakable grain. By reasonable designing of fragment size and shape, the simulation of particle crushing process can be very close to laboratory experiments [18,28,29]. Besides, the agglomerate method is also appropriate for the simulation of single-particle crushing. A famed agglomerate model, proposed by Potyondy [30], is known as the bonded particle method (BPM). Zhang et al. [31] proposed densely arranged BPM by introducing densely arrayed elementary balls into the BPM to simplify the calculation process and avoid the problem of inner void releasing. The replacement method was firstly put forward in MD by Tsoungui [24] and then introduced into the DEM. In this method, the crushed grains will be replaced by a group of smaller grains [24,25,32]. Compared with the agglomerate method, the number of particles needed in the replacement method is much smaller, which accelerates the computation speed of the replacement model [33]. It is worth noting that the grain failure criterion and replacing scheme are the most important issues in the replacement method [34,35]. Several crushing criteria have been purposed and introduced into the replacement method [4,22,36]. For both agglomerate and replacement methods, the accuracy of DEM modeling is mainly based on the rationality and reliability of laws and characteristics of single-particle crushing strength and carrying capacity.
The Brazilian disc test is the most convenient and widely accepted way to investigate the crushing strength of single particles [1][2][3][4]. Plenty of Brazilian disc tests have been conducted in the past decades, and statistical data show that the distribution of grain crushing strength can be described by the Weibull distribution [1,4,6,37]. Experimental results indicate that the Weibull characteristic strength of large-sized particles is lower than that of the smaller ones. The explanation for this phenomenon is that large-sized particles contain more inner defects, such as tiny flaws and weak planes, which may grow into macro-cracks under loadings [5,38]. In order to simulate the size effect of particles, Lim et al. [39] proposed a method, in which the original flaws can be produced by deleting 0% to 25% elementary balls randomly. Wang et al. [40] and Kuang et al. [41,42] reproduced the distribution of single-particle strength based on Lim's method and investigated the relationship between particle strength and coordination number, while in fact, the inner porosity of most granular materials is less than 5%, and the introduction of original flaws by removing elementary balls will result in a statistical deviation of both void ratio and equivalent fragment diameter in DEM simulations. Therefore, the rationality of deleting elementary balls deserves further discussing.
Strength is an indigenous mechanical property of particles, and carrying capacity is the external expression of particle strength. Particle carrying capacity depends on both material strength and contact or loading conditions. Experimental investigations on the relationship between particle carrying capacity and its contact conditions are difficult to be conducted for the lack of specialized experimental equipment. Numerical experiments mainly focus on the effect of particle coordination number on carrying capacity. Some numerical investigations indicate that particle carrying capacity depends linearly on its coordination numbers, no matter whether the distribution of contact points is symmetric or not [40].
In this work, inhomogeneous bond stiffness and strength are introduced into the BPM, and the Weibull distribution, as well as the size effect of single-particle strength, is reproduced with this method. Several numerical experiments under complex contact conditions, including symmetric loading tests, asymmetric loading tests, and ball-ball Brazilian tests, are carried out to investigate the relationship between single-particle carrying capacity and contact conditions. The results show some interesting rules that differ from the conventional assumptions, and suggestions on the crushing criterion of particle crushing simulation with DEM are concluded at the end of this work.

Stress and Carrying Capacity
Granular material is a kind of non-continuum medium; hence, the stress of particles needs to be defined specially. In the DEM theory, the stress tensor of particles can be defined as Cundall and Strack [43] suggested:σ whereσ ij is the asymmetric average stress tensor of an individual ball; N c is the number of particle contact points, or the coordination number of a particle; x (c) i and F (c) j are the coordinate component and contact force component acting on the contact point, respectively; and V is the volume or area of particles.
Simplified particle stress formula is preferred in laboratory experiments, and tensile stress is assumed as the primary cause of particle crushing in the Brazilian tests. Jaeger [44] proposed a widely used formula to define the tensile strength of single particles: where F f is the maximum reaction force acting on two loading points, and d is the diameter of a particle. It should be noted that Eq. (2) is only appropriate for the Brazilian tests, and the application of Eq.
(2) under other contact conditions may result in intrinsic mistakes as particle crushing strength and criterion under complex loading conditions remain controversial. In some opinions, particle crushing is closely related to the inner tensile stress. Wang et al. [45] came up with a formula to measure the maximum tensile stress σ cen at the grain center: in which σ max and σ min are the maximum and minimum principal stress, respectively, while in other opinions, shear stress is regarded as the determining factor of particle crushing. McDowell [46] took the peak octahedral stress q as the stress index of particle strength: in which σ 1 , σ 2 , and σ 3 are the three principal stresses. The principal stresses in Eqs. (3) and (4) can be derived from Eq. (1). It should be noticed that, in both Wang's and McDowell's theories, the criterion stress of a particle is assumed to be a constant even under different loading paths. Particle carrying capacity can be measured by the maximum loading force [40][41][42]. In this work, the average loading force F ca is selected to represent particle carrying capacity: where N c is the coordination number of a particle and F c is the contact force vector acting on a contact point.

Single-Particle Crushing Strength
As a kind of heterogeneous materials, an individual particle always contains several kinds of minerals, binders, and internal flaws, and these defects always induce a dispersion of mechanical properties of grains. For the purpose of describing crushing strength of particles, a particle model based on densely arrayed BPM considering the inhomogeneity of material properties is introduced in this section, and the distribution and size effect of particle crushing strength are reproduced by simulations of the Brazilian test.

DEM Modeling with Inhomogeneous Properties
As the deviation of both void ratio and fragment equivalent diameter will be introduced, it is inappropriate to simulate the inhomogeneous clump by deleting elementary balls. In this section, the inhomogeneity of materials is achieved by introducing inhomogeneous material parameters into densely arrayed BPM, with which the inherent nature of particle crushing strength can be reproduced without introducing undesirable statistical deviation. As granular material is a kind of friction materials, the cohesive force between particles and contacting objects can be ignored, and the contacting force can be described by friction coefficients and contacting stiffness. The parallel-bond model is selected as the clump bond model in this work for its advantage in the DEM simulation for geotechnical materials [20,47]. The parallel bond can be defined by six parameters [48]: bond normal stiffnessk n , bond shear stiffnessk s , tensile strengthσ c , cohesion c, friction angel φ, and bond radiusr. The bond shear strengthτ c can be calculated from cohesion, friction angel and normal stress. If the maximum tensile stressσ max exceeds bond normal strengthσ c , or the maximum shear stressτ max exceeds bond shear strengthτ c , parallel bonds will break, and the accompanying force, moment and stiffness of those bonds will vanish. The selection of parameters in this work takes both DEM theory and experimental results into consideration, and follows the suggestions of Lisjak [17] and Potyondy [20] at the same time. After a number of attempts, parameters of elementary balls and bonds are finally determined, as shown in Table 1. The inhomogeneously distributed stiffness and strength of each bond obey Gaussian distribution, with values ranging from 0.4 to 1.6 times of the average values in Table 1, as Fig. 1a shows. Bonds with low strength or low stiffness represent weak planes inside particles, and the initial structural defects can be introduced by those inhomogeneous bond parameters. Seventy Brazilian tests are carried out to reproduce the distribution of particle crushing strength, as Fig. 1b shows. The diameters of particles range from 1.4 to 1.7 mm, and the diameter of uniformly sized elementary balls is 0.1 mm. The loading speed of top wall is 1.0 mm/min toward the particle center.

Strength Distribution and Size Effect
Numerous experimental studies indicate that single-particle crushing strength obeys the Weibull distribution law. Nakata [1] proposed a formula to describe the strength distribution of uniformly sized particles: Quartz tests [3] DEM simulations Probability of surviving (%) where P s is the survival probability of a particle under stress σ, m is the Weibull modulus, and σ 0 is the characteristic stress under which 37% of the tested particles survive.
The normalized crushing strength of seventy particles is drawn in Fig. 2, in contrast to Nakata's experimental results [3]. The features of two curves are similar. The Weibull modulus and characteristic strength are 3.51 and 32.92 MPa in the DEM simulation and 3.04 and 30.96 MPa in Nakata's experiment, respectively. The numerical result indicates that, in a statistical sense, particle strength distribution can be reproduced by densely arranged BPM by considering the inhomogeneity of bond parameters.
The characteristic strength of particles is size-dependent. In order to investigate the relationship between particle diameter and characteristic strength, more Brazilian test simulations are conducted with different-sized particles. The relationship between particle size and characteristic strength, as well as Weibull modulus are plotted in double logarithmic coordinates as Fig. 3 shows, in comparison with the quartz tests of Nakata [3]. It can be seen from Fig. 3a that particle characteristic tensile strength is exponentially dependent on particle size, and similar relationship is also reflected between Weibull modulus and particle size.
Therefore, in summary, by introducing inhomogeneous bond properties to the DEM modeling, particle strength distribution and size effect can be reproduced well with the BPM. The DEM simulation results present similar characteristics and changing tendency with laboratory tests, indicating that the densely arrayed BPM considering inhomogeneous stiffness and strength properties is reliable and robust in single-particle crushing tests. Unless particularly stated, the models and material properties in the following numerical experiments are the same with those in this section.

Wall-Ball Loading
The amount and distribution of coordination points present significant influence on particle carrying capacity. In this section, a number of wall-ball loading tests, including symmetric loading tests and asymmetric loading tests, are conducted to investigate the relationship between particle carrying capacity and its contact conditions.

Symmetric Loading
As the simplest contact condition, symmetric contact condition is chosen to investigate the relationship between particle carrying capacity and coordination number. Five groups of spherical particles are crushed in symmetric loading path with regular polyhedron boxes, and the shapes of boxes are regular tetrahedron, hexahedron, octahedron, dodecahedron and icosahedron, respectively. In these simulations, the peak average reaction force F ca , the maximum octahedral stress q, and the maximum   grain center tensile stress σ cen are taken into statistics. Similar to particle tensile strength, particle carrying capacity can also be described by the Weibull distribution. Figure 4 records the characteristic values of F ca0 in five groups of simulations. It can be noticed that F ca0 monotonously increases with N c , and particle carrying capacity and its coordination number are almost linearly related. Figure 5 exhibits the relationship between characteristic octahedral stress q 0 and characteristic grain center tensile stress σ cen0 in symmetric loading tests. Two stress indexes present the same linearly changing trends with the increase of coordination number. It should be pointed out that the grain center tensile stress σ cen in all symmetric loading tests are negative, which indicates that no tensile stress exists in particle center when crushing occurs. Therefore, the assumption that particle crushing is caused by particle central tensile stress seems to be incorrect. Besides, the approximately linear relationship reflected in Fig. 5 indicates that stress indexes are inappropriate for particle crushing judgment. Because when particles are about to break under different loading conditions, stress indexes cannot keep constant, which runs counter to the definition and assumption of crushing criterion.
A particle, generated with random seed 1031, is selected to investigate the crushing behavior of single particles. The crushing states of five particles when F ca reaches their peaks are shown in Fig. 6, in which the small disks inside clumps are discrete fracture network (DFN), and the red and black disks  represent tensile and shear bond failures, respectively. Coordination number also shows influence on the category of particle crushing. The inner particle crushing can be divided into three categories [49]: fracture, attrition and abrasion. In Fig. 6, the fractures in tetrahedron, hexahedron, and octahedron mainly concentrate in certain areas near the contact points, and the main category of failure is attrition. While the distribution of fractures in dodecahedron and icosahedron is more homogeneous, and with the increase of coordination number, the type of particle crushing changes from fracture and attrition to abrasion gradually.
The peak values of F ca and the corresponding amounts of DFN are listed in Table 2. It can be noticed that the selected particle shows some features that differ from statistical laws. Unlike the characteristic reaction force F ca0 , the peak reaction force shows a wavelike rising tendency with the increase of coordination number. The peak reaction forces of regular tetrahedron, hexahedron, and octahedron keep at relatively low level, as well as the number of fractures, while the peak average reaction force in dodecahedron and icosahedron tests exhibits a leapfrogging increase, and the particle carrying capacities are much higher than those in other three tests.
To sum up, the numerical results of symmetric loading experiments indicate that coordination number has significant influences on the failure category and carrying capacity of particles. Particles with a larger coordination number are harder to be crushed integrally and can bear more cracks before losing their strengths. Particle carrying capacity, as well as stress indexes, shows a positive and almost linear relation with the coordination number in symmetric tests, which indicates that the stress indexes cannot be regarded as a particle crushing criterion. Besides, for one single particle, this monotone and positive correlation may not be correct all the time. Size can only decide the characteristic strength of particles, rather than the crushing strength and carrying capacity for one certain particle.

Asymmetric Loading
The particle crushing issues discussed above are all under absolutely symmetric contact conditions. Although the simplified crushing criterions are very popular in the replacement method, the symmetric loading condition is infrequent in actual grain crushing cases. In this part, asymmetric loading tests, including asymmetric multipoint loading tests and non-parallel Brazilian tests, are carried out to investigate the influence of loading asymmetry on particle crushing.
The asymmetric multipoint loading tests are realized by deleting one of the loading walls in symmetric multipoint tests randomly. As the tetrahedron box cannot hold a sphere clump any more after losing a plane, it is not considered in the asymmetric tests. The characteristic reaction force in asymmetric loading tests is recorded in Fig. 7, compared with symmetric loading tests. The average force F ca and particle coordination number are positively related in the asymmetric tests, similar to those in Fig. 4. The values of F ca in asymmetric tests are much lower than those in symmetric tests. This is in accordance with the common sense, but differs from the results of former research [40].
Asymmetric contact conditions always induce asymmetric stress distribution. The stress indexes in asymmetric tests are recorded in Fig. 8. The characteristic octahedral stress q 0 presents an increase with the increase of coordination number, while the grain center tensile stress σ cen shows totally different changing tendency. In the four series of simulations, there exists positive σ cen in Fig. 8b, which means that tensile stress appears in some particles. However, as the sign of σ cen keeps changing, it is impossible to describe σ cen with the Weibull distribution. Besides, the values of stress indexes are much lower than those in the symmetric tests, and also do not keep constant.
DFNs of the particle with random seed 1031 under asymmetric multipoint loading conditions are imaged in Fig. 9. The reaction force, bond failure events, and stress indexes of this particle are listed in Table 3. Compared with Fig. 6 and Table 2, particles under asymmetric contact conditions lose much carrying capacity, but hold more fractures before being crushed. Similar to the symmetric tests, the category of particle crushing in asymmetric tests changes from fracture and attrition to abrasion with the increase of coordination number.
In some replacement methods, particles with the coordination number of two are assumed to obtain the same carrying capacity ignoring the symmetry of contact points. In this work, the assumption is   tested by non-parallel Brazilian experiments. Several groups of spherical particles are crushed by nonparallel loading planes, and the asymmetry is characterized by the dip angle α of the top plane, which ranges from 2 • to 12 • , as Fig. 10 shows. DFN connects two loading points and the center of particles in non-parallel Brazilian tests. The characteristic reaction force in non-parallel Brazilian tests is recorded in Fig. 11. It can be seen that F ca0 of non-parallel Brazilian tests is lower than that of the conventional Brazilian tests and reduces monotonically with the increase of dip angle. Moreover, the angle between two planes shows remarkable influence on particle carrying capacity, and the assumption that particles with two contact points have the same carrying capacity is inconsiderable. The results of asymmetric multipoint loading tests and non-parallel Brazilian tests indicate that asymmetric contact conditions will lower the carrying capacity of particles. Stress indexes in asymmetric loading tests are also not constants, which is in accordance with the results derived from Sect. 4.1.
The assumption that particles with two coordination points have constant carrying capacity is also invalid in this section. In other words, criterion stress may only be applicable for some specific occasions, but no longer appropriate for general particle crushing cases.

Ball-Ball Loading
For granular materials, the proportion of wall-ball contact is relatively small, while the ball-ball contact accounts for the majority of contact. The crushing strength and carrying capacity of particles are associated with the relative size of their neighbors, and the relative size can be quantified by radius ratio R r as follows: R r = r n r (7) in which r and r n are the radii of an individual particle and its neighboring particle, respectively. The ball-ball Brazilian tests are conducted as Fig. 12 shows, in which a crushable clump is clamped and crushed by two rigid balls. Eleven groups of particles with radius ratio ranging from 0.1 to 10 are crushed in the ball-ball Brazilian tests, and their characteristic reaction forces are recorded in Fig. 13. The characteristic values of F ca0 monotonously increase with the radius ratio R r , which is consistent with the understanding that particles are more easily to be crushed by sharper objects. The relationship between F ca0 and R r is approximately linear in semi-logarithmic coordinate. DFNs of the particle with random seed 1031 in the ball-ball tests are imaged in Fig. 14. Most bond failures occur within a narrow columnar space, and the crushing category of particles in all ball-ball Brazilian tests is fracture absolutely. The evolution of DFN distribution indicates that radius ratio shows an impact on the occurrence and evolution of particle inner cracks.
The position where internal fracture initiates in the Brazilian test remains controversial. Hudson [50] found that the initiation of crack always locates near the loading points in the Brazilian disc tests, and this point of view is supported by the work of Klaus [51] and Swab [52]. Another viewpoint is that the position of the first microcrack locates at the maximum shear stress point, which is about 0.2 times of particle radius away from the loading points, and this viewpoint receives theoretical and experimental support from Markides [53], Li [54], and Lin [55]. Both opinions above are manifested by the ball-ball crushing tests in this work. It can be noticed in Fig. 8 that cracks grow from both ends to the center of sphere clumps if R r is less than 2.0, while bond failures mainly concentrate in the central area of particles if R r is greater than 2.0. Moreover, the fracture region moves from two loading points to the particle center with the increase of R r , which indicates that the crack growth paths are not unique, but significantly affected by contact conditions.

Conclusion
This paper mainly focuses on investigating the problem of single-particle crushing strength and the influence of particle contact conditions on carrying capacity. From the numerical experiments and analysis, the following conclusions can be drawn.
1. Densely arrayed BPM considering the inhomogeneous bond parameters shows its advantage in reproducing the Weibull distribution and size effect of single-particle crushing strength. This method does not need making flaws by deleting elementary balls or introducing extra inner voids. 2. Particle carrying capacity increases linearly with coordination number in both symmetric and asymmetric multi-position loading tests. The asymmetry of contact conditions reduces particle carrying capacity significantly. Both results of symmetric and asymmetric tests show that the relationship between particle carrying capacity and coordination number cannot be simply described as linear.
The octahedral stress and grain center tensile stress are inaccurate as a particle crushing criterion. For a single particle, equal coordination number does not mean the same carrying capacity. The oversimplified assumptions and crushing stress criteria are proved to be inconsistent with facts. 3. Contact conditions present influence on particle crushing categories. For particles with lower coordination number, fracture and attrition are the main crushing categories, while with the increase of coordination number, abrasion becomes the main crushing category, and particles can bear more cracks before losing their carrying capacity in such cases. 4. The particle carrying capacity is exponentially related with the relative size of its neighbors. Besides, radius ratio presents obvious effect on the origin and evolution of particle inner cracks. The crack growth paths in ball-ball loading tests are not unique, and the origin position of initial cracks moves gradually from loading points to particle center with the increase of radius ratio.
The DEM numerical studies considering inhomogeneous bond parameters exhibit notable laws and features of single-particle crushing in this work. Particle crushing strength and carrying capacity are not only related to particle inherent nature, but also influenced by contact conditions. This work may help establish more reasonable particle crushing criteria, and the connections between crushing issues of single particles and granular materials still demand further study.