Analysis of One-Dimensional Compression under a Wide Range of Stress with Densely Arrayed BPM

Tao Zhang; Ke Xu; Wenxiong Huang

doi:10.32604/cmes.2021.015406

1College of Mechanics and Materials, Hohai University, Nanjing, 211100, China
2School of Civil Engineering, Xuchang University, Xuchang, 461000, China
*Corresponding Author: Tao Zhang. Email: zht_hhu@163.com
Received: 16 December 2020; Accepted: 19 February 2021

Abstract: In this paper, the densely arrayed bonded particle model is proposed for simulation of granular materials with discrete element method (DEM) considering particle crushing. This model can solve the problem of pore calculation after the grains are crushed, and reduce the producing time of specimen. In this work, several one-dimensional compressing simulations are carried out to investigate the effect of particle crushing on mechanical properties of granular materials under a wide range of stress. The results show that the crushing process of granular materials can be divided into four different stages according to er-logσy curves. At the end of the second stage, there exists a yield point, after which the physical and mechanical properties of specimens will change significantly. Under extremely high stress, particle crushing will wipe some initial information of specimens, and specimens with different initial gradings and void ratios present some similar characteristics. Particle crushing has great influence on grading, lateral pressure coefficient and compressibility of granular materials, and introduce extra irreversible volume deformation, which is necessary to be considered in modelling of granular materials in wide stress range.

Keywords: Densely arrayed BPM; particle crushing; gradation curve; lateral pressure coefficient; compressibility; extra volume deformation

Particle crushing plays an important role in the variation of granular material properties, such as gradation, compressibility and hydraulic conductivity [1]. In order to improve engineering safety and durability of granular materials, numerous experimental and theoretical studies on particle crushing have been conducted in the past decades.

Laboratory experiment is the most direct way to study the mechanical response of granular materials. Since Terzaghi et al. [2] and Athy [3] found the phenomenon of particle crushing in 1920s, various experiments, such as uniaxial compression [4], triaxial compression [5] and impact test [6], had been conducted to investigate the particle crushing laws, as well as the relationship between particle crushing and material properties from both micro and macro perspectives [7,8]. McDowell et al. [9] and Nakata et al. [10,11] pointed out that, for single particles, the reference strength and particle size are negatively correlated, and particle strength distribution can be described by Weibull distribution [12]. Tsoungui et al. [13] found that smaller particles are more easily to be crushed under compression tests of granular specimens. Valdes et al. [14] proposed a concept of yield stress for granular materials, and studied the gradation evolution of particles during crushing process. Though X-ray and CT technology have been introduced into particle crushing study [15–18], the change of some mechanical properties remains difficult to be acquired directly, and for the limitation of equipment and budget, complex experimental conditions are usually unavailable in laboratories.

Numerical experiment is an important supplemental method for laboratory tests. DEM, proposed by Cundall et al. [19], has been widely used to investigate the mechanical responses of granular materials [20,21]. By using DEM, the mechanical response of specimens becomes visible, and complicated and extreme loading conditions are available. After decades of development, DEM has been well-established, and some DEM software, such as PFC, Yade and MatDEM [22], have been widely used in the modeling of granular materials.

In the simulation of particle crushing with DEM, the key is to build crushable particle models reasonably. Two types of crushable particle models have been proposed [23]: the replacement method [24] and the agglomerate method [25,26]. In the former method, whether a particle is crushed is judged by a pre-defined criterion, and the crushed particles will be replaced by a per-set cluster of fragments, as for instance, Tsoungui et al. [13], Lobo-Guerrero et al. [27,28], McDowell et al. [29], Ciantia et al. [30] and Cil et al. [16,31]. The advantage of replacement method is its high computation efficiency, and the disadvantage is that, the crushing criterion and constitution of the replacing cluster are subjective and controversial. Besides, the replacement method can only be used in the study of granular specimens, and is inapplicable to single particle crushing issues. In agglomerate method, a macro crushable particle is formed by a cluster of bonded micro particles, and the failure of bonds will result in the breakage of macro particles. The agglomerate method can be used in single particle crushing tests, moreover, complicated particle shapes, not limited to sphere and circle, are also available. Unlike the replacement method, subjective crushing criterion and fragment formation are not required in agglomerate method. One of the most popular agglomerate models is bonded particle model (BPM), presented by Potyondy et al. [26]. Ge et al. [32,33] validated BPM with polymer crushing experiments with 3D printing technology, and revealed the failure laws of inter-particle bonds. Thornton et al. [34], Cheng et al. [35,36], Bolton et al. [37], and Alonso et al. [38] studied the problems of single particle crushing, complicated stress paths loading and other granular material crushing issues with BPM. In BPM simulation, the damage evolution and crack growth of granular materials are visible, while the huge computing scale and long model generating time bring about the consequence of lower computational efficiency, especially for 3D modeling.

The problem of void ratio, which is often ignored in the present replacement and agglomerate methods, is of great importance in particle crushing simulation with DEM. In replacement method, if a crushed particle is replaced by circular fragments within its outer boundary, the conservation of mass will fail, and specimen void ratio will change abruptly. Inversely, if the law of mass conservation is kept, large unbalanced contact force will be introduced due to the inevitable overlaps between replacing fragments [39]. In agglomerate method, the crushable clusters are porous, and the internal porosities of clusters are always greater than 0.2. When one cluster is crushed, the internal pore may be exposed and becomes external pore which can be filled by other fine fragments. While in fact, the internal pores in most geomaterial particles are too tiny to be filled by any fragment. It is worth noting that, both of the volumes of solid parts and internal pores should be considered in the calculation of fragment equivalent diameter, otherwise, the equivalent grain size will be smaller than the actual fragment size. Therefore, it is necessary to introduce corrections to current DEM models considering particle crushing.

In this paper, a special BPM, densely arrayed BPM, is proposed to simulate the crushing of particles. This method can calculate the releasing volume of inner voids more accurately and accelerate the generation of granular samples. Basic features of this model, especially the correction of fragment equivalent diameter and internal pore, are firstly introduced in detail. Then one-dimensional compression tests are simulated with densely arrayed BPM to characterize the development of particle crushing and its effects on mechanical properties granular materials in a wide stress range. The simulation results show that, the compression tests can be divided into four stages with different characteristics, and specimens with different initial states present similar mechanical properties under high stress level conditions. The influence of particle crushing on granular materials mainly focuses on fragment gradation, lateral pressure coefficient, compressibility and extra deformation.

Densely arrayed BPM is formed by the same sized micro-elements arrayed in the densest way. The unified and specific geometrical relationship between those micro-elements makes it possible to estimate the volume of internal pores. How to stack balls in the densest way is known as “Kepler’s Conjecture,” which has been strictly proved that triangular lattice is the densest stack of equal-sized disks [40], and the face-centered-cubic (FCC), as well as the hexagonal-close-packed (HCP) are the densest stacks of equal-sized balls [41], as shown in Fig. 1. In densely arrayed BPM, the geometrical relationship between micro-elements is simple and clear, and it is easy to get coordinate of each micro-element with iterative algorithm. Besides, overlaps and contact force between micro-elements are not exist in the initial state, and clumps are in stress equilibrium state. Therefore, it costs little time for cluster to reach a balance state.

Parallel bond is widely utilized in BPM, and the behavior of bonds can be defined by five parameters [25,26]: normal and shear bond stiffness, kn and ks (stress/displacement); normal and shear strength, σcb and τcb (stress); and bond radius r¯. Details of the bond behavior are introduced in Potyondy’s work [42]. Once the maximum tensile stress σmaxb exceeds normal strength σcb, or the maximum shear stress τmaxb exceeds shear strength τcb, bonds will break, and the acting force, moment and stiffness of bonds will vanish. The selection of parallel bond parameters is based on the mechanical response of granular materials in experiments [32,33] and the fundamental principle of DEM [1]. If τcb is set as infinity, shear failure will be excluded, and similarly, if σcb is set as infinity, tensile failure will be excluded. Therefore, both tensile and shear failure will be possible if σcb and τcb are set reasonably.

In order to solve the problem of calculation error of equivalent grain size, correction of equivalent particle diameter is introduced in this work. In the present numerical study, the statistics of fragment size is based on calculation of diameters of the equivalent volume or circle area (DECA). The equivalent diameter Deq of a fragment can be calculated as:

where Veq is the equivalent volume of a fragment, and Aeq is the equivalent area. It is worth noting that the equivalent volume Veq and the total solid volume of micro-elements VS of a fragment are not equal. As shown in Fig. 2, the volume of a unit cell includes both solid part and internal voids. If only take VS into consideration, the calculated equivalent diameter will be smaller than actual fragment size. The relationship between solid volume VS and equivalent volume Veq can be written as following:

In densely arrayed BPM, formula (2) is universally applicable, even for one single free micro-element. Otherwise, conservation of mass and volume will be violated after particle crushing. In other words, the space taken by every micro-element is assumed a little larger than its actual size. In the following simulations, DECAs of particles and fragments are all calculated by formulas (1) and (2).

Figure 2: The difference between Veq and VS, Aeq and AS: (a) a unit cell of FCC takes the space of a regular hexahedron; (b) HCP takes the space of a hexagonal prism; (c) in 2D models, a unit cell takes the area of a regular hexagon

Densely arrayed BPM can avoid the problem of internal pore releasing. Due to the dense arrangement of elements, the released inner voids will not be filled by any fragments after crushing, as Fig. 3 shows. As the internal pore (the red part) is still a forbidden zone for micro elements, it can be treated as a special solid material rather than a pore. Based on the same theory, the volume of solid parts can also be corrected by formula (2). Densely arrayed BPM clusters include no internal pores in some sense, and if necessary, such as modeling a porous particle, internal pores can be obtained by deleting certain micro elements.

Figure 3: Inner pores will be released during crushing, but no fragment can enter the red region

It should be noted that, for the absence of theoretical dense arrangements of polydisperse spherical granular systems and non-spherical granular systems, the analyses about DECA correction and void ratio, as well as the problem of void releasing in this work, are only appropriate for densely arrayed BPM. The correction of polydisperse spherical and non-spherical granular systems still demands further research.

For densely arrayed BPM, the average stress of a particle can be defined in Ehlers’ way [43]:

where σ¯ij is a non-symmetric tensor that represents the average stress of a particle [44], Nc is the number of particle contact points, xi(c) and Fj(c) are the coordinate components and contact force components acting on contact points, V is the volume or area of a particle. In plane problems, with the help of Mohr’s circle, σmax and σmin can be obtained by following [45]:

where σc=σ11+σ222, and τc=σ12+σ212. The max norm of shear stress τmax can be calculated by:

τmax= max(∥τc+(σ11-σc)2+(σ12-τc)2∥,∥τc-(σ11-σc)2+(σ12-τc)2∥).(5)

To investigate the crushing rules of granular materials, DEM simulations of silica one-dimension compression tests were carried out. The contact and bond parameters were calibrated by a single particle crushing test, and for the purpose of saving computing resource, all the numerical tests were modeled in 2D condition.

In this section, a particle was placed and crushed between two rigid parallel plates, as shown in Fig. 4a. The diameter of the particle is 1.9 mm, and it is formed by 388 micro-elements with radius of 0.05 mm. Bond parameters used in this test are listed in Tab. 1. In order to reduce overlaps, a non-linear contact model, Hertz-Mindlin model, is selected as the ball-ball contact model, and the contact stiffness parameter is relatively high. The shear modulus of micro-elements is 2,800 GPa, and the Poisson’s ratio is 0.3. The density of micro-element is 2,630 kg/m3, and the friction factor is 0.8. The particle was crushed into two large and several small fragments as shown in Fig. 4. In Fig. 4b, the red, yellow and black short lines represent activated bonds, tensile failure and shear failure, respectively. It can be seen that bond failure penetrates the particle and mainly occurs within a rectangular region in the middle area.

Figure 4: Single particle crushing test: (a) Initial state and crushed state; (b) Inside bond failure

The reaction force vs. wall displacement curve in the test is presented in Fig. 5. The numerical result matches well with the experimental data of quartz obtained by Nakata et al. [10]. The peak reaction force in Nakata’s experiment is about 165 N, and the value in DEM simulation is 162 N, with a relative error of about 2%. Though there is a certain deviation between numerical and experimental results in the second half, the reaction forces in numerical and experimental tests both drop rapidly after their peak points, and the residual strength of both tests keep in low level, which is in accordance with the characteristic of brittle fracture.

Figure 5: Comparison between densely arrayed BPM simulation and Nakata’s experiment [10]

In this section, two numerical tests were designed to study the evolution of particle crushing in granular specimens under one-dimensional compression, as shown in Figs. 6a and 6b. Both rectangular crushing chambers were formed by four rigid plates. The diameters of crushable particles range from 0.25 to 2.0 mm in specimen 1, and 1.4 to 1.7 mm in specimen 2. All the material properties are the same with those used in the single particle crushing test above. Both specimens were generated in natural stacking method. The marginal parts of both specimens are very unconsolidated, as shown in Fig. 6a. Initial void ratios of specimens are about 0.6. The loading was applied by moving the top plate downwards at a speed of 0.1 mm/min, and the final state of both samples are shown in Figs. 6c and 6d. Due to the nonlinear ball-ball contact model and high contact stiffness, overlaps between micro-elements are not obvious, and present indistinctive effect on deformation of specimens throughout the tests.

Figure 6: One-dimension compression tests with densely arrayed BPM. (a) Specimen 1: d=0.25~2.0 mm. (b) Specimen 2: d=1.4~1.7 mm. (c) Final state of specimen 1. (d) Final state of specimen 2

On account of the different definition of void ratio in 2D and 3D models, and for the convenience of study, the concept of relative void ratio er is proposed as:

where e is the current void ratio, e0 is the initial void ratio. The change of er can be regarded as an index of volume change in both 2D and 3D specimens. The er-log σy curves of both specimens are drawn in a semilogarithmic coordinate system and compared with Nakata’s experimental data [46], as shown in Fig. 7. It can be seen that the simulation results are close to experimental data. Although the quantitative relationship between 2D and 3D models is ambiguous, it is still worth drawing some qualitative conclusions through the simplified numerical experiments.

Figure 7: The comparison between er-logσy curves of DEM simulations and Nakata’s experiment [46]

The compression process can be divided into four stages with distinctive features, as shown in Fig. 8. In Stage I, from point A to B, the stress grows slowly and keeps at a low level, and no bond failure occurs. The deformation of specimen in this stage contains two parts: restorable part (elastic deformation of bonds), and unrestorable part (the sliding of clusters). The framework structure of clumps is still precarious. In Stage II, from point B to C, the mean vertical stress σy and the number of bond failure increase slowly. The structures of two specimens are rebuilt in response to particle cracking and sliding. In Stage III, from point C to D, a large number of bond-failures occur, which occupy 80% of the whole tests. The specimens start compressive hardening, which makes the vertical stress increase monotonically. The slope of er-log σy curve increases markedly after point C. At the end of this stage, all small particles are crushed, with only several larger particles remaining relatively intact. The gaps between large particles are filled with tiny debris, and the specimens become much denser than initial state. This phenomenon coincides with the experiment observations of Tsoungui [13].

Figure 8: Void ratio and bond failure curves of one-dimension compression tests. (a) Specimen 1 (b) Specimen 2

The upper limit of loading stress in laboratory tests is always no more than 100 MPa [10,16], but in numerical simulations, higher loading stress is available, which makes it possible to look into the mechanical response of granular materials in a wide stress range. In Stage IV, specimens show some interesting phenomenon. From point D to E, the vertical stress of specimens increases by more than twenty times, while the amount of bond failures increases only by about twenty percent. Specimen stiffness improves significantly after point D, and er-log σy curves become much flatter than those in Stage III. At the end of compression tests, the relative void ratios of the two specimens are about 0.05, and the specimens are hard to be compressed even under extremely high vertical stress.

The form of particle crushing in each stage differs from each other, and can be mainly classified into three categories [47]: (1) fracture, a particle breaks into several parts in similar size; (2) attrition, a particle breaks into a large fragment and several small ones; (3) abrasion, only a little of powder is worn off from the particle. Fig. 9 images the state of bond failure during compression, and the meaning of colors in Fig. 9 is the same with Fig. 4b. It can be seen that, in Stage II, bond failures only occur in individual particles, and the main crushing type in this period is fracture and attrition, as Figs. 9c and 9h show. In Stage III, the shape and size of particle fragments become random and no longer smooth, then attrition and abrasion become the main types of crushing, as Figs. 9d and 9i show. Medium-sized particles, especially those under bad contact conditions or concentrated stress, are crushed, and small particles that remain intact in the gaps of big ones break into several parts after the specimen structural rearrangement. In the last stage, the specimen structure and force chain become stable, and abrasion becomes the main crushing type, as Figs. 9e and 9j show.

Figure 9: The crushing process of two specimens. (a) Point A1 (b) point B1 (c) point C1 (d) point D1 (e) point E1 (f) point A2 (g) point B2 (h) point C2 (i) point D2 (j) point E2

The change of specimen gradation is one of the most significant influence of particle crushing. For the convenience of analysis, a normalized dimensionless parameter β, proposed by Valdes et al. [14], is introduced to represent the level of stress as following:

in which the yield stress σyield corresponds with the maximum curvature point in er-log σy curve in semi-logarithmic coordinate system. For example, in Fig. 8, point C can be regarded as the yield point, and the yield stresses of two specimens are 8.17 and 14.27 MPa, respectively.

Fig. 10 images the gradation evolution with different β. The evolution characteristics of these curves correspond with the crushing types of each stage mentioned above. Gradation curves of both specimens change significantly in Stage III, and the cumulated weight of medium and small-sized fragments presents a remarkable increase, especially in specimen 2. In Stage IV, although the increase of vertical stress is tremendous, the change of gradation curve is inapparent, and only the number of tiny fragments shows a tiny increase.

Figure 10: Fragment cumulated weight distributions of two specimens with different β

In spite of the difference of initial gradings, the gradations of two specimens evolve in similar rules, and the maximum DECAs of two specimens keep almost unchanged. This phenomenon also can be found in Figs. 6 and 9. Several large particles with large coordination numbers remain almost intact at the end, and only a few micro-elements are worn off. According to some conclusions of single particle tests [9–11], the strength of larger particles tends to be lower than smaller ones, but in actual compression tests, larger particles are more likely to survive at last. It is for the reason that particle bearing capacity not only depends on its material and shape, but also its coordinate conditions. Stress distribution in large particles with larger coordination number is always more uniform than that in small ones. As non-uniform stress distribution will raise the level of shear stress, it is a disadvantage to the carrying capacity of grains. Besides, it is worth nothing that, with the increasing vertical stress, the two series of curves with different initial gradations evolve convergently. This phenomenon indicates that, after being fully crushed, the size distributions of granular materials with different initial grading conditions share some common characteristics. In other words, under high stress level, there might be a theoretical ultimate gradation curve for granular materials with different initial state.

Lateral pressure coefficient K0 is an important mechanical parameter of geotechnical materials. In this work, K0 was calculated as the ratio of mean normal stresses acting on lateral and axial walls. Fig. 11 records the K0 curves of two specimens. It can be seen that both curves oscillate randomly, and present inconspicuous evolution rule in the first two loading stages. Evolution of lateral pressure coefficient in Stage III is much obvious and can be divided into two parts. Firstly, particle crushing lowers the stress level acting on sidewalls. When relative void ratio is greater than about 0.3, K0 keeps between 0.2 and 0.4. Secondly, in the latter half of Stage III, the lateral pressure coefficients of two specimens show a general increasing trend and remain at about 0.7.

Due to particle crushing and specimen densification, K0 curves keep fluctuating in a small range throughout the test. In spite of the difference under low stress level, K0 curves of two specimens still show some common characteristics. Evolution laws of K0 under high stress level prove again that properties of granular materials with different initial states become similar after being fully crushed.

The change of lateral pressure coefficient results from three aspects. Firstly, high pressure and particle crushing make specimen denser than initial state, and the improvement of relative compaction can influence specimen lateral pressure coefficient. Secondly, particle crushing changes the size distribution of grains, which is closely associated with specimen lateral pressure coefficient. Thirdly, particle crushing and structure rebuilding will induce a change in specimen frictional angle, which is one of the determinants of specimen lateral pressure coefficient.

Particle crushing can impact the compressibility and resiliency of granular materials, and induce extra deformation, which may bring safety threat to rockfill. The issue of compressibility can be investigated by comparison between unloading and reloading experiments, and the problem of extra deformation can be investigated by comparison between crushable and uncrushable compression simulations.

In this section, an unloading and reloading test was carried out based on specimen 1to study the impact of particle crushing on compressibility and resiliency of granular materials. The relative void ratio and bond failure curves under cycle loading are marked in green in Fig. 12. Simulation result shows that, the unloading springback is very tiny, and for the emergence of extra irreversible deformation, there exists a wide gap between the unloading and initial loading curves. Deformation of granular materials is almost unrecoverable after being crushed. In the early stage of reloading, the reloading curve retraces the unloading curve, and the slopes of both er-log σy curves are generally the same. In the second half of reloading, the gap between two curves is widened gradually. The divarication appears before the occurrence of new bond failure, and then bond failure enlarges the difference between unloading and reloading curves. Two factors contribute to the gap growth: densification of granular materials and particle crushing. Densification widely exists in cycle loading tests of geotechnical materials, even for uncrushable materials [48,49], and particle crushing expands the gap.

To study the extra deformation caused by particle crushing, a compression test of uncrushable specimen was carried out based on specimen 1, and the bond strength parameters of uncrushable specimen were set to be extremely high values to avoid any bond failures. The er-log σy curve of uncrushable specimen is drawn in Fig. 13a, with a comparison with the crushing result of a crushable specimen. The two curves are coincident before point B. The void ratio decrement of uncrushable specimen is inconspicuous when vertical stress is less than about 1,000 MPa, and in this period, it is mainly caused by particle sliding, rotating and rearrangement. While under high stress level, the void ratio shows a significant drop, as elastic deformation of particles plays a leading role in this period. Besides, the yield stress of uncrushable specimen is about 1,000 MPa. It is interesting that the slopes of two er-log σy curves are similar after yield points. As most granular materials are fragile, the distortion of uncrushable particle after the yield point is not a typical deformation and insignificant in the analysis of extra crushing deformation.

K0 curves of uncrushable and crushable specimens are drawn in Fig. 13b. Similar to the er-log σy curves, K0 curves overlap with each other in the early stage, and exhibit different increasing features after the appearance of bond failures. The K0 curve of uncrushable specimen keeps an increasing trend in general, and doesn’t present obverse changing stages throughout the tests, while the curve of crushable specimen shows a moderate growth, especially in the third stage. This phenomenon indicates that particle crushing slows the growth of stress acting on side walls.

Figure 13: Simulation results of uncrushable and crushable specimens: (a) er-logσy curves; (b) K0 curves

The gap between crushable and uncrushable curves represents the extra deformation brought by particle crushing. The extra decrement of relative void ratio with vertical stress less than 1,000 MPa is shown in Fig. 14, segmented by the four stages in Fig. 8a. The extra deformation changes indistinctively in the first two stages, but presents a rapid increase in the third stage, and reaches its peak at the end of Stage III. The curve of extra void ratio decrement flattens gradually in the latter half of Stage III. In Stage IV, the deformation of crushable specimen becomes stable, and uncrushable specimen deformation remains low before the yield point. Therefore, the extra void ratio presents an indistinctive decrement in this stage. It is can be concluded that, the general tendency of extra deformation growth rate firstly ascends and descends at last. As unrecoverable extra deformation is similar to plastic deformation to some degree, the strain of crushable granular materials can be calculated by summing the extra strain to the strain of uncrushable model.

This paper mainly focuses on DEM simulations of particle crushing in a wide stress range. As a customized bonded-particle model, the densely arrayed BPM can solve the issue of inner-pore releasing and equivalent diameter calculating by introducing simple and reasonable correction. Evolution of granular material properties under compression is analyzed, and the following conclusions can be drawn:

1. The compressing process of granular materials can be divided into four different stages. Particle crushing mainly occurs after the yield point in the third stage. Hence, the change of material properties in this stage is more significant. Particle crushing and its effects become insignificant when β is larger than about 50.

2. Under high stress level, properties of granular materials become stabilized and convergent, and the information of initial void ratio and gradation will be wiped by extreme high-stress. The simulation results indicate a possible existence of theoretically ultimate state for fully crashed granular materials.

3. Particle crushing has great influence on lateral pressure coefficient of granular materials, and reduces the level of stress acting on side walls. K0 curves of two specimens share similar characteristics when extensive crushing occurs. K0 curves show an upward trend generally after the yield point, and reach a stable state at the last stage when particles are fully crushed.

4. Particle crushing brings considerable extra and unrecoverable deformation to granular materials. The value of extra deformation firstly increases and then decreases with the crushing process, and reaches its peak at the terminus of the third stage. The evolution law of extra crushing deformation may help the constitutive modeling of granular materials considering particle crushing.

In brief, with the modification of granular specimen generating, as well as the correction of void ratio and DECA statistics, densely arrayed BPM is an appropriate method for numerical analyses of particle crushing. But the computational cost will be higher as more micro-elements are employed to form a macro crushable specimen. However, the simulation results and conclusions will bring inspiration to the constitutive modelling of crushable granular materials.

Funding Statement: The authors wish to thank the National Natural Science Foundation of China (No. 11772117), the Fundamental Research Funds for the Central Universities (No. 2015B37414), Henan Scientific and Technical Project under Grant (No. 192102310480) and Key Scientific Research Project of Colleges and Universities in Henan Province (CN) (21B560015) for financial support.

Conflicts of Interest: The authors declare that they have no conflicts of interest to report regarding the present study.

References

1. Lisjak, A., Grasselli, G. (2014). A review of discrete modeling techniques for fracturing processes in discontinuous rock masses. Journal of Rock Mechanics and Geotechnical Engineering, 6(4), 301–314. DOI 10.1016/j.jrmge.2013.12.007. [Google Scholar] [CrossRef]

2. Terzaghi, K., Peck, R. B., Mesri, G. (1995). Soil mechanics in engineering practice. 3rded. USA: John Wiley & Sons, Inc. [Google Scholar]

3. Athy, L. F. (1930). Density, porosity and compaction of sedimentary rocks. Aapg Bulletin, 14(1), 1–24. DOI 10.1306/3D93289E-16B1-11D7-8645000102C1856D. [Google Scholar] [CrossRef]

4. Hendron, A. J. (1963). The behavior of sand in one-dimensional compression (Ph.D. Thesis). University of Illinois, USA. [Google Scholar]

5. Lee, K. L., Farhoomand, I. (1967). Compressibility and crushing of granular soil in anisotropic triaxial compression. Canadian Geotechnical Journal, 4(1), 68–86. DOI 10.1139/t67-012. [Google Scholar] [CrossRef]

6. Lv, Y., Wang, Y., Zuo, D. (2019). Effects of particle size on dynamic constitutive relation and energy absorption of calcareous sand. Powder Technology, 356(11), 21–30. DOI 10.1016/j.powtec.2019.07.088. [Google Scholar] [CrossRef]

7. Karimpour, H., Lade, P. V. (2010). Time effects relate to crushing in sand. Journal of Geotechnical and Geoenvironmental Engineering, 136(9), 1209–1219. DOI 10.1061/(ASCE)GT.1943-5606.0000335. [Google Scholar] [CrossRef]

8. Lee, K. L., Seed, H. B. (1967). Drained strength characteristics of sands. ASCE Soil Mechanics and Foundation Division Journal, 117–141. DOI 10.1061/JSFEAQ.0001048. [Google Scholar] [CrossRef]

9. McDowell, G. R., Bolton, M. D. (1998). On the micromechanics of crushable aggregates. Géotechnique, 48(5), 667–679. DOI 10.1680/geot.1998.48.5.667. [Google Scholar] [CrossRef]

10. Nakata, Y., Hyde, A. F., Hyodo, M., Murata, H. (2001). A probabilistic approach to sand particle crushing in the triaxial test. Géotechnique, 49(5), 567–583. DOI 10.1680/geot.1999.49.5.567. [Google Scholar] [CrossRef]

11. Nakata, Y., Kato, Y., Hyodo, M., Hyde, A. F. L., Murata, H. (2001). One-dimensional compression behaviour of uniformly graded sand related to single particle crushing strength. Soils and Foundations, 41(2), 39–51. DOI 10.3208/sandf.41.2_39. [Google Scholar] [CrossRef]

12. Weibull, W. A. (1950). A statistical distribution function of wide applicability. Journal of Applied Mechanics, 18, 293–297. [Google Scholar]

13. Tsoungui, O., Vallet, D., Charmet, J. C. (1999). Numerical model of crushing of grains inside two-dimensional granular materials. Powder Technology, 105(1–3), 190–198. DOI 10.1016/S0032-5910(99)00137-0. [Google Scholar] [CrossRef]

14. Valdes, J. R., Koprulu, E. (2008). Internal stability of crushed sands: Experimental study. Géotechnique, 58(8), 615–622. [Google Scholar]

15. Cil, M. B., Alshibli, K. A. (2012). 3D assessment of fracture of sand particles using discrete element method. Géotechnique Letters, 2(3), 161–166. DOI 10.1680/geolett.12.00024. [Google Scholar] [CrossRef]

16. Cil, M. B., Alshibli, K. A. (2014). 3D evolution of sand fracture under 1D compression. Géotechnique, 64(5), 351–364. [Google Scholar]

17. Karatza, Z., Andò, E., Papanicolopulos, S. A., Viggiani, G., Ooi, J. Y. (2019). Effect of particle morphology and contacts on particle breakage in a granular assembly studied using X-ray tomography. Granular Matter, 21(3), 44. DOI 10.1007/s10035-019-0898-2. [Google Scholar] [CrossRef]

18. Zhai, C., Pagan, D. C., Hurley, R. C. (2019). In situ X-ray tomography and 3D X-ray diffraction measurements of cemented granular materials. Journal of the Minerals, Metals & Materials Society, 72(1), 18–27. DOI 10.1007/s11837-019-03774-4. [Google Scholar] [CrossRef]

19. Cundall, P. A., Strack, O. D. L. (1979). A discrete numerical model for granular assemblies. Géotechnique, 29(1), 47–65. DOI 10.1680/geot.1979.29.1.47. [Google Scholar] [CrossRef]

20. Zheng, J., An, X., Huang, M. (2012). GPU-based parallel algorithm for particle contact detection and its application in self-compacting concrete flow simulations. Computers & Structures, 112–113(4), 193–204. DOI 10.1016/j.compstruc.2012.08.003. [Google Scholar] [CrossRef]

21. Zheng, Z., Zang, M., Chen, S., Zeng, H. (2018). A GPU-based DEM-FEM computational framework for tire-sand interaction simulations. Computers & Structures, 209(4), 74–92. DOI 10.1016/j.compstruc.2018.08.011. [Google Scholar] [CrossRef]

22. Liu, C., Xu, Q., Shi, B., Deng, S., Zhu, H. (2017). Mechanical properties and energy conversion of 3D close-packed lattice model for brittle rocks. Computers & Geosciences, 103, 12–20. DOI 10.1016/j.cageo.2017.03.003. [Google Scholar] [CrossRef]

23. Xiao, Y., Desai, C. S., Daouadji, A., Stuedlein, A. W., Liu, H. et al. (2020). Grain crushing in geoscience materials-key issues on crushing response, measurement and modeling: Review and preface. Geoence Frontiers, 11(2), 363–374. DOI 10.3389/fgene.2020.00363. [Google Scholar] [CrossRef]

24. Åström, J. A., Herrmann, H. J. (1998). Fragmentation of grains in a two-dimensional packing. Physics of Condensed Matter, 570(83), 551–554. DOI 10.1007/s100510050476. [Google Scholar] [CrossRef]

25. Potyondy, D., Autio, J. (2001). Bonded-particle simulations of the in-situ failure test at Olkiluoto. Rock Mechanics in the National Interest, 2, 1553–1560. [Google Scholar]

26. Potyondy, D. O., Cundall, P. A. (2004). A bonded-particle model for rock. International Journal of Rock Mechanics and Mining Sciences, 41(8), 1329–1364. DOI 10.1016/j.ijrmms.2004.09.011. [Google Scholar] [CrossRef]

27. Lobo-Guerrero, S., Vallejo, L. E. (2005). Discrete element method evaluation of granular crushing under direct shear test conditions. Journal of Geotechnical and Geoenvironmental Engineering, 131(10), 1295–1300. DOI 10.1061/(ASCE)1090-0241(2005)131:10(1295). [Google Scholar] [CrossRef]

28. Lobo-Guerrero, S., Vallejo, L. E., Vesga, L. F. (2006). Visualization of crushing evolution in granular materials under compression using DEM. International Journal of Geomechanics, 6(3), 195–200. DOI 10.1061/(ASCE)1532-3641(2006)6:3(195). [Google Scholar] [CrossRef]

29. McDowell, G. R., de Bono, J. P. (2013). On the micro mechanics of one-dimensional normal compression. Géotechnique, 63(11), 895–908. DOI 10.1680/geot.12.P.041. [Google Scholar] [CrossRef]

30. Ciantia, M. O., Arroyo, M., Calvetti, F., Gens, A. (2015). An approach to enhance efficiency of DEM modelling of soils with crushable grains. Géotechnique, 65(2), 91–110. DOI 10.1680/geot.13.P.218. [Google Scholar] [CrossRef]

31. Cil, M. B., Buscarnera, G. (2016). DEM assessment of scaling laws capturing the grain size dependence of yielding in granular soils. Granular Matter, 18(3). DOI 10.1007/s10035-016-0638-9. [Google Scholar] [CrossRef]

32. Ge, R., Wang, L., Zhou, Z. (2019). DEM analysis of compression breakage of 3D printed agglomerates with different structures. Powder Technology, 356, 1045–1058. DOI 10.1016/j.powtec.2019.08.113. [Google Scholar] [CrossRef]

33. Ge, R., Ghadiri, M., Bonakdar, T., Zheng, Q., Zhou, Z. et al. (2020). Deformation of 3D printed agglomerates: Multiscale experimental tests and DEM simulation. Chemical Engineering Science, 217(1558), 1–15. DOI 10.1016/j.ces.2020.115526. [Google Scholar] [CrossRef]

34. Thornton, C., Yin, K. K., Adams, M. J. (1996). Numerical simulation of the impact fracture and fragmentation of agglomerates. Journal of Physics D Applied Physics, 29(2), 424–435. DOI 10.1088/0022-3727/29/2/021. [Google Scholar] [CrossRef]

35. Cheng, Y. P., Nakata, Y., Bolton, M. D. (2003). Discrete element simulation of crushable soil. Géotechnique, 53(7), 633–641. DOI 10.1680/geot.2003.53.7.633. [Google Scholar] [CrossRef]

36. Cheng, Y. P., Bolton, M. D., Nakata, Y. (2004). Crushing and plastic deformation of soils simulated using DEM. Géotechnique, 54(2), 131–141. DOI 10.1680/geot.2004.54.2.131. [Google Scholar] [CrossRef]

37. Bolton, M. D., Nakata, Y., Cheng, Y. P. (2008). Micro-and macro-mechanical behaviour of DEM crushable materials. Géotechnique, 58(6), 471–480. DOI 10.1680/geot.2008.58.6.471. [Google Scholar] [CrossRef]

38. Alonso, E. E., Tapias, M., Gili, J. (2012). Scale effects in rockfill behaviour. Géotechnique Letters, 2(3), 155–160. DOI 10.1680/geolett.12.00025. [Google Scholar] [CrossRef]

39. Ben-Nun, O., Einav, I. (2010). The role of self-organization during confined comminution of granular materials. Philosophical Transactions of the Royal Society a Mathematical Physical and Engineering Sciences, 368(1910), 231–247. DOI 10.1098/rsta.2009.0205. [Google Scholar] [CrossRef]

40. Rogers, C. A. (1964). Packing and covering. England: Cambridge University Press. [Google Scholar]

41. Hales, T. (2005). A proof of the Kepler conjecture. Annals of Mathematics, 162(3), 1065–1185. DOI 10.4007/annals.2005.162.1065. [Google Scholar] [CrossRef]

42. Potyondy, D. O. (2015). The bonded-particle model as a tool for rock mechanics research and application: Current trends and future directions. Geosystem Engineering, 18(1), 1–28. DOI DOI 10.1080/12269328.2014.998346. [Google Scholar] [CrossRef]

43. Ehlers, W., Ramm, E., Diebels, S., D’Addetta, G. A. (2003). From particle ensembles to Cosserat continua: Homogenization of contact forces towards stresses and couple stresses. International Journal of Solids and Structures, 40(24), 6681–6702. DOI 10.1016/S0020-7683(03)00418-9. [Google Scholar] [CrossRef]

44. Huang, W., Huang, L., Sheng, D., Sloan, S. W. (2014). DEM modelling of shear localization in a plane Couette shear test of granular materials. Acta Geotechnica, 10(3), 389–397. DOI 10.1007/s11440-014-0348-6. [Google Scholar] [CrossRef]

45. Iordache, M. M., Willam, K. (1998). Localized failure analysis in elastoplastic Cosserat continua. Computer Methods in Applied Mechanics and Engineering, 151(3–4), 559–586. DOI 10.1016/S0045-7825(97)00166-7. [Google Scholar] [CrossRef]

46. Nakata, Y., Hyodo, M., Hyde, A. F. L., Kato, Y., Murata, H. (2001). Microscopic particle crushing of sand subjected to high pressure one-dimensional compression. Soils and Foundations, 41(1), 69–82. DOI 10.3208/sandf.41.69. [Google Scholar] [CrossRef]

47. Guyon, E., Troadec, J. P. (1994). Du sac de billes au tas de sable. France: Odile Jacob Science. [Google Scholar]

48. Liu, H., Xiao, J., Wang, P., Liu, G., Gao, M. et al. (2018). Experimental investigation of the characteristics of a granular ballast bed under cyclic longitudinal loading. Construction and Building Materials, 163(8), 214–224. DOI 10.1016/j.conbuildmat.2017.12.037. [Google Scholar] [CrossRef]

49. Wu, Z. X., Yin, Z. Y., Dano, C., Hicher, P. Y. (2020). Cyclic volumetric strain accumulation for sand under drained simple shear condition. Applied Ocean Research, 101, 102200. DOI 10.1016/j.apor.2020.102200. [Google Scholar] [CrossRef]