Abstract

This paper presents numerical simulations of high-pressure biaxial tests on breakable granular soils with the discrete element method. The 2D setting is more economic in terms of computational cost, which allows simulation with a larger number of particles with a wider size distribution. The results of breakable and unbreakable agglomerates show that particle breakage has a significant influence on the macro- and micromechanical behaviors of the assembly. Higher confining pressure and larger axial strain result in the variation of particle grading and agglomerate numbers. The evolution of bond breakage during shearing makes it possible to trace the failure process and breakage mechanism at the microlevel. The breakage energy is found to account for a small fraction of total energy input compared with friction energy. A hyperbolic correlation between relative particle breakage and total energy input per unit volume was established regardless of the influence of confining pressure.

1. Introduction

Particle breakage frequently occurs in granular materials as the exterior energy acting on granular particle exceeds its strength. This phenomenon prevalently encountered in civil, petroleum, and mining engineering, such as crushing in end-bearing piles, rock-fill dams, oil wells, and mine construction [15]. Even some natural soils composed of high-strength mineral materials, such as deep underground granular soils, may undergo serious breakage under high pressure [6]. Up until now, more than 500 vertical shafts exceeding 500 meters were constructed in northern China. It is inevitable that the disturbance in the process of construction dramatically changes the stress field in soils, including serious particle breakage. As a cohesionless granular soil, the property of sand particle is relatively simpler than clay or silt, such as its coarser particle size and poor plasticity. Thus, this paper makes a deep and systematical research on particle breakage under high-pressure shear test, taking sand particle as a research object.

Grading change due to particle breakage has a great influence on the macroscopic behavior of a granular soil. This is true because particle breakage may create a more drastic change in the internal structure than can be achieved by rearrangement [7]. For conventional triaxial tests, crushing at sliding contacts or breakage of particles affects stress-strain behavior and decreases the rate of dilation [813]. Bolton [14] showed that strength and dilatancy of sand are influenced by relative density and stress level, relating to particle breakage. Coop et al. found that a constant grading can be reached at very large strains on carbonate sands, even quartz sands at low stress levels, were subjected to small amounts of particle breakage [15, 16]. Wu et al. [17] and Hyodo et al. [18] found that the large amount of particle breakage resulting from the high-level strain was induced by particle transformation and rotation during shearing. Yu [12] conducted a great deal of triaxial tests under various influence factors (e.g., confining pressure, initial void ratio, and drainage condition) and discussed the effects of these factors on particle breakage and whole volume change. From what have been investigated above, particle breakage is principally influenced by particle strength and effective stress state, especially at high pressures. Despite these considerable studies, the micromechanics of particle breakage in granular soils is not well understood because it is difficult to observe the crack propagation, movement, and internal interaction at the particle level. Alikarami et al. [11] successfully applied X-ray microtomography to the crushing process at the scale of particle and described the evolution of shear bands. Hall et al. [19] quantified the onset and evolution of localized deformation processes in sand with particle-scale resolution and confirmed the importance of particle rotations associated to strain localization with X-ray microtomography imaging. However, the grading change during shearing was not considered in these studies in which limited loads were permitted in this environment.

Discrete element method (DEM) has proved to be a powerful numerical methodology developed and applied for simulating granular materials. For the crushable soils, different methods have been used in DEM: (1) forming agglomerates by cementing elementary balls (rigid and unbreakable). The bonds between elementary balls are deformable and breakable [7, 2024], and (2) replacing large particles by a group of small ones when a certain strength is reached [2528]. In addition, combined DEM and FEM method was adopted as an effective method for simulation of polygon-shaped particles [29, 30]. Raisianzadeh et al. [31] proposed a combined DEM and XFEM approach to simulate particle interaction by DEM and the breakage analysis by XFEM of angular particles. The crushing mechanism of two-dimensional circular particles has been studied with the first method. Moreno et al. [32] investigated the influence of impact angle on the breakage characteristics of circular agglomerates. Kun and Herrmann [33] adopted PFC2D to simulate the crushing process and obtained size distribution of fragments for an agglomerate under impact. Ueda et al. [23] investigated the effects of particle shape on the breakage behavior of granular materials. The researches mentioned above have greatly improved our understanding of the crushing mechanics in granular materials, but a limited number of particles with a poor size distribution were allowed.

The main purpose of this study is to simulate the process of particle breakage in high-pressure biaxial tests on sands with large amounts of particles, using PFC2D with the first method. The following study provides both macro- and micromechanical behaviors for an assembly of DEM agglomerates in the sand samples. First, the macromechanical behaviors are discussed to prove the correctness and feasibility of this method, in which particle breakage is verified as an indispensable phenomenon in high-pressure shear tests. Then, the micromechanical behaviors, such as anisotropy of force distribution, crack evolution, and energy dissipation due to particle breakage, are investigated at great length. Unbreakable agglomerates are also examined in a parallel study to investigate the effects of particle breakage on the strength, dilatancy, and anisotropy.

2. Numerical Simulation Methodology

2.1. DEM Method for a Sand Sample

The contact models used in this study are the same with those in Xu et al. [24]. The linear parallel bond model provides the behavior of two interfaces: (1) an infinitesimal, linear elastic (no-tension), and frictional interface that carries a force and a finite-size, linear elastic, and (2) bonded interface that carries a force and moment. The first interface is equivalent to the linear model: it does not resist relative rotation, and slip is accommodated by imposing a Coulomb limit on the shear force. The second interface is called a parallel bond because when bonded, it acts in parallel with the first interface. When the second interface is bonded, it resists relative rotation, and its behavior is linear elastic until the strength limit is exceeded and the bond breaks, making it unbonded [34]. In the linear model, the force-displacement law for the normal and tangential components is given bywhere kn and ks are the constant normal and shear stiffnesses, is the surface gap, and Δδs is the adjusted relative increment of shear displacement. The linear stiffnesses, kn and ks, can be inherited from the contacting particles; thus, the linear stiffnesses are converted into the following equations:where the superscripts (1) and (2) denote the properties of particle 1 and 2, respectively. Moreover, the linear stiffnesses can be described by the elastic constants of Young’s modulus (E) and Poisson’s ratio (ν). E and ν are emergent properties that can be related to the effective modulus (E∗) and the normal-to-shear stiffness ratio (κ∗) at the contact. E∗ and κ∗, as the arguments of the linear model method, are given bywhere L is the distance between the centers of contacting particles or the distance between the center of the particle and wall and A is the projection area of the smaller diameter of contacting particles in PFC2D, as shown in Figure 1.

A parallel bond provides the mechanical behavior of a finite-sized piece of cement-like material deposited between contacting particles. This makes parallel bonds that can transmit both force and moment between particles. The parallel bond is defined by bond effective modulus (), bond normal-to-shear stiffness (), tensile strength (σc), and cohesion (c); if the tensile strength or shear-strength limit is exceeded, then the bond breaks.

2.2. Modelling of a Sand Agglomerate

In two-dimensional DEM simulation of biaxial test, thousands of agglomerates need to bond together to form a numerical sample. An agglomerate comprising several elementary balls will increase the computation time significantly due to the additional searches for bonds. For this reason, agglomerates with less number of elementary balls were preferred. Not only does this have the advantage of less computation time but also can get reasonable results [23, 35]. As shown in Figure 2, the agglomerate consists of a circular assembly of 7 elementary balls, and on this basis, the simulation results would be reasonable and the computation time would still be shortened. For the naturally granular soils, the strength of granular soils may vary widely due to the complex internal flaws. It has been shown that strengths of granular soils follow a Weibull distribution [36, 37]. Robertson and Bolton [20] firstly performed crushing tests on single agglomerates with random bonding flaws. Cheng et al. [22] randomly removed 20% of elementary balls in an agglomerate to provide a statistical variability to the strength of the agglomerates. Laufer [38] investigated high-pressure compression tests on brittle grains and found that the crushing strengths of the grains are directly proportional to the bond strengths, including normal distributed bond strengths, degraded bonds, and uniformly increased/decreased bond strengths.

In this study, increased bond strengths of tensile strength (σc) and cohesion (c) were set to get the random strength of the agglomerates. The mean value of σc and c was 2.3 × 107 N/m2 and 7 × 106 N/m2, respectively. Agglomerate details of a typical sand grain are given in Table 1. Figure 3 shows a typical simulation result from the compression test; it is compared with the crushing test of a silica sand investigated by Nakata et al. [22, 39]. Their variable trends have no great difference, both of them appear as lower peaks (at a and b in Figure 3) before reaching their final peak stresses. From the result of numerical simulation, evolution of boundary energy and energy dissipated by breakage are also obtained at the microlevel. Bolton et al. [7] considered the changes in four components of work and energy during simulated tests on DEM agglomerates: δW, work input at the boundaries; δS, energy stored as elastic strain energy; δB, energy dissipated by breakage; and δF, energy dissipated by friction and damping due to the rearrangement of balls. Applying the first law of thermodynamics to these components we obtain:

Similarly, the linear parallel bond model in PFC provides the following energy partitions: δW, work input at the test boundaries; δS, energy stored as elastic strain energy; δf, energy dissipated by frictional slip; and δd, energy dissipated by damping so that energy dissipated by breakage δB, can be obtained:

Figure 4 shows that energy dissipated by friction and breakage is a small fraction compared with the total boundary energy in this specific agglomerate. Boundary energy variation with axial strain agrees well with a square-law increase. Energy dissipated by friction increases slowly followed by a sharp increase at the maximum peak stress. The stepped increase of energy dissipated by breakage indicates that a broken bond appears at an inflection point and only bond breakage can increase breakage energy. The evolution of energy can be divided into two typical stages, prebreakage stage and breakage stage (at stages 1 and 2, respectively). At stage 1, almost all of the boundary energy is instantaneously stored as strain energy without bond breakage. At stage 2, with the breakage of the agglomerate, a small fraction of boundary energy is dissipated by breakage, and more friction is added between the elementary balls. However, most of the boundary energy is still stored as strain energy in this specific agglomerate.

The survival probability of a batch of 30 agglomerates is calculated using the mean rank position and following the procedures in Cheng et al. [22]:where i is the rank position of a grain when the grains are sorted into ascending order of peak stress and N is the number of tests (N = 30). The Weibull distribution can be used to describe the variability in tensile strength of apparently identical test-pieces of a brittle material [40]:where σ0 is the characteristic stress at which 1/e or 37% of the grains survive and m is the Weibull modulus. To estimate the Weibull modulus, (5) can be rewritten as

The tensile strength of the agglomerates is assessed with (7):where F represents the greatest force and d is the distance between the compression platens at the beginning of the crushing test.

Figure 5 shows that a similar shape of the survival probability curve is obtained for both numerical simulation and Weibull distribution when the crushing peak stress is normalized with σ0. The Weibull modulus is approximately 4.83, which sensibly lies in the region of 2.5 to 5.0 for sand grains reported by Robertson and Bolton [20].

2.3. Sample Preparation

To simulate dense sand in biaxial tests, an assembly of 1758 unbreakable agglomerates, which is called clump in PFC, was firstly generated in four walls, and its original particle size distribution is shown in Figure 6. The initial dimension of the specimen was 39.1 mm × 80 mm. They were placed randomly and cycled to equilibrium so as to exclude overlaps. During this stage, effective modulus was set to a low value, and friction coefficient between clumps was reduced to zero. The next step was to replace all pebbles with balls (a clump is a rigid collection of n rigid spherical pebbles), which have the same diameters and coordinates while replacing the pebbles. The final effective modulus of the balls of the agglomerates was set to final value (4 × 108 N/m2). Then, the unbreakable agglomerates translated into breakable agglomerates, which is called cluster in PFC.

After the specimen was generated, it was isotropically compressed until a complete equilibrium condition was reached. The top and bottom walls moved towards each other with a constant speed of 0.1 mm/s, and the lateral walls were controlled by a servo algorithm to keep the confining stresses constant during the biaxial tests. The confining stresses used in this study were 2 MPa, 3 MPa, 4 MPa, 5 MPa, and 6 MPa.

2.4. Biaxial Compression Test

Subsequently, the sample was sheared in a strain-controlled mode by moving the top and bottom walls towards each other at a specific velocity of 0.5 mm/s. Therefore, the servo-controlled mode on the top and bottom walls was canceled, and the state of lateral walls remained unchanged. The simulation results of macro- and micromechanical behaviors will be discussed, respectively, in the following sections.

3. Macromechanical Results

3.1. Results of Breakable Agglomerates

As a vitally important parameter, confining pressure has a significant effect on granular soil behavior. A series of numerical biaxial tests were conducted after isotropic consolidation under constant confining pressures (σc = 2 MPa, 3 MPa, 4 MPa, 5 MPa, and 6 MPa).

As shown in Figure 7(a), the sample stress-strain responses show a linear initial trend followed by a curvature to the peak stress. As anticipated, higher confining pressures led to higher mobilised deviator stresses. The relationship between volumetric strain and axial strain is shown in Figure 7(b). For all samples, the volumes contract evidently at first, and the volumetric strains reduce with increasing confining pressures. With the development of shearing, all samples begin to expand until finally. At relatively low confining pressure levels (σc = 2 MPa, 3 MPa), the sample exhibits evident dilatancy after contraction. Under higher confining pressures (σc = 5 MPa, 6 MPa), more obvious volumetric contraction occurs as the axial stains reach nearly 10%. This is due to the closer contact and more dramatic particle breakage, which makes particles break into smaller ones. Therefore, smaller particles fill the original voids of sample, which consequently has a more compression behavior.

The peak friction angle φp is the maximum value of mobilised friction angle φm:where σa and σc are the major and minor effective stresses, respectively. As shown in Figure 8, the peak friction angles are 27.59°, 27.45°, 25.66°, 25.04°, and 23.65° with the increase of confining pressures, and the reduction trend observed in this study is consistent with experimental tests by Corfdir and Sulem [41].

3.2. Results of Unbreakable Agglomerates

In order to understand more deeply the importance of particle breakage, shearing tests with unbreakable agglomerates were simulated. Corresponding results for unbreakable agglomerates are given in Figure 9. It can be seen that the assemblies with unbreakable agglomerates have a more dilative behavior. For both deviator stress and volumetric responses, the difference between breakable and unbreakable agglomerates becomes more evident as the confining pressure becomes larger. This is due to the tendency of more breakage under higher confining pressures. This behavior agrees with experimental results [42] and DEM simulation results [7, 26, 43].

Figure 8 shows the comparison of peak friction angles between breakable and unbreakable agglomerates under different confining pressures. In contrast, an almost constant friction angle for an unbreakable agglomerate is obtained in a wide range of confining pressure. Thus, the reduction trend of friction angle for breakable agglomerates can be attributed to breakage of sand particles during shearing.

4. Micromechanical Results

4.1. Evolution of Anisotropy

During shearing, granular soils present continuous change in the evolution of interparticle forces and internal geometry. Furthermore, changes in the number of contacts and distribution of orientations are happening synchronously. Rothenburg and Bathurst [44] showed that the distributions of normal and tangential contact forces can be approximated by the following analytical expressions:where the coefficients an and at define the magnitude of the directional variation of normal and tangential forces, θn is the major principal direction of normal force, θt is the direction in which tangential force is zero on average, and is the average normal force.

Figure 10 shows the rose diagrams of normal force in the breakable assembly with σc = 4 MPa. At the beginning (axial strain εa = 0%), the assembly is geometrically isotropic after isotropic compression. With the growth of axial strain, the shape shows a polar distribution of the anisotropy, which grows to be deformed such as a peanut. The coefficient of normal force anisotropy (an) and maximum normal force () is found to increase during the initial state up to 4% of axial strain whereas the macroscopic softening begins. After this stage, an and exhibit a reduction trend after the maximum value; this is due to loss of contacts and the loss of capacity of particle chains to sustain high forces [43].

The corresponding tangential components of contact forces are shown in Figure 11. At the initial state (axial strain = 0%), the tangential force approaches zero with σc = 4 MPa. This reflects the fact that the assembly is geometrically isotropic, and zero tangential force does not imply the absence of tangential force. As simulation continues, the angular distribution of tangential forces shows symmetrical peak values oriented at about 45° to principal stress direction. Similarly, the coefficient of tangential force anisotropy (at) and maximum tangential force () is found to increase to a peak value followed by a slow decline in magnitude. The initial increase in at is due to the translational movement of agglomerates. As the number of contacts approaches the ultimate value, the agglomerates are more likely to mobile and rotate with loss of contacts, and this loss results in a slowly release in tangential force.

Figure 12 presents the comparison of anisotropy in breakable and unbreakable samples at the same axial strain of 4%. Because the comparison of anisotropy in breakable and unbreakable samples with the same confining pressure is very similar, the paper only presents the comparison of anisotropy under the confining pressures of 4 MPa, 5 MPa, and 6 MPa. Through the comparison of Figures 12(a) and 12(b), the anisotropy coefficient of normal force (an) in unbreakable samples is larger than that in breakable samples with the same confining pressure. Also in Figures 12(c) and 12(d), the anisotropy coefficient of tangential force (at) shows the same trend by contrasting the difference in breakable and unbreakable samples. The results follow a very similar trend in Hosseininia and Mirghasemi [43].

4.2. Particle Breakage Behavior

In real triaxial tests, it is difficult to obtain the evolution of particle size distribution during shearing. Yu [12] carried out a great deal of triaxial tests on silica sand with the purpose of investigating the characteristic of particle breakage. To obtain the particle size distributions, all tests were stopped at certain axial strains (i.e., 0%, 5%, 10%, 15%, and 20%). As an effective tool, DEM simulation plays an important role in this study. As shown in Figure 13, the red and black bond stripes between balls represent linear parallel bond and linear bond, respectively, and the dotted line represents bond breakage or a microcrack in an agglomerate. When the microcracks appear, as shown in Figure 13(b), the parallel bond breaks and the bond material is removed from the model along with its accompanying force, moment, and stiffness. However, particle A still belongs to the original agglomerate, and there is no change in the particle grading. Figure 13(c) shows a penetrating crack in an agglomerate, from which particle A and particle B are thoroughly separated. Naturally, a new agglomerate consisting of particle A and particle B is generated.

The final particle grading for all tests are plotted in Figure 14. It is clearly observed that the particle grading transforms to a relatively more well-graded distribution along with the increase of confining pressures. This is similar to the evolution of real particle size distributions. Meanwhile, it can also be found that the size of large-diameter agglomerates remains nearly constant, and the small-diameter ones undergo significant changes. This similar phenomenon has also been studied in experimental tests [12, 45] and numerical simulations [46].

In this study, the Hardin’s relative breakage index, Br, was adopted to quantify the particle breakage. As shown in Figure 15, particle breakage increases with an increasing axial strain, and it is more obvious particle that breakage occurs under higher confining pressure. This trend agrees with the observation of experimental tests presented by Luzzani and Coop [15] and Lade et al. [47].

In the process of shearing, the original agglomerates simultaneously break up into several new agglomerates, which may crush again after rearrangement. Figure 16 shows the evolution of agglomerate numbers with axial strain under different confining pressures. It can be seen that the number of agglomerates almost linearly increases with axial strain after axial strain reaches nearly 4%. When the confining pressure is 2 MPa, the agglomerate numbers increase from 1758 to 1981. At the end of the test, only 223 new agglomerates separate from the original agglomerates. With the further increase of confining pressures, there are 983, 1682, 2617, and 3361 newly separated agglomerates generated at the end of tests, respectively. From this point of view, the extent of breakage has a strong relationship with confining pressure.

Form the numerical simulation, we can also obtain information about the evolution of microcracks, which represent bond breakage and internal fracture in agglomerates. As the number of broken bond increases, an original agglomerate gradually splits into several fragments. Figure 17 shows the increment of broken bonds in every 1000 steps and evolution of broken bonds with axial strain under different confining pressures. In Figure 17, there is no bond breakage initially until the axial strains reach nearly 1.5%–2.5%. It is noticeable that as the confining pressure increases, the time of the first bond breakage occurrence is earlier. This is due to the closer contact and greater contact force between agglomerates at the initially condition with an increase in confining pressures. After this, the evolution of accumulated broken bonds increases at the same axial strain, the trend agrees well with a hyperbolic increase. It is also shown that the increment of broken bonds is the same as acoustic emission (AE) count variation with time, which can explain the failure process and breakage mechanism of sand. Generally, the entire processes were featured by substantial bond breakage due to progressive crack growth in agglomerates. Similar observations were also reported by Chen et al. [48] and Qin et al. [49] in triaxial tests of granular soils.

4.3. Evolution of Energy

As mentioned by Lade et al. [47], the total amount of energy input per unit volume during a shear test is the sum of energy input during isotropic compression and shearing stages:where ET is the total energy input per unit volume, EC is the energy input during isotropic compression, and ES is the energy input during shearing. Miura and O-hara [50] and Liu et al. [51] proposed to replace total energy input with plastic work because the elastic energy in samples with particle breakage is neglected compared to the magnitude of plastic work. Herein, boundary energy (EB = ET × volume) indicating work done by all walls is adopted to monitor the total energy input, which can be obtained from DEM simulation conveniently.

Figure 18 presents evolutions of total energy input and energy dissipated by breakage (Eb) during shearing. There exists a roughly linear relationship between energy input and axial strain. In the initial stage of shearing, energy input after isotropic compression is larger with larger confining pressures, and energies dissipated by breakage from all samples are around zero until the axial strain reaches approximately 4%. Then, energies dissipated by breakage increase slowly with the increase of axial strain. It can also be seen that energy input and energies dissipated by breakage increase with an increasing confining pressure. This trend agrees well with the experimental observations presented by Miura and O-hara [50].

Figure 19 indicates the change in the ratio of energies dissipated by breakage and energy input (Eb/EB) with developing axial strain. It can be seen that there are three stages in these curves. In the first stage, at small axial strain, energies dissipated by breakage are zero because particle breakage has not yet happened or particle breakage is negligible; after the axial strain comes to 4%, the ratio begins to increase rapidly with an increase in the axial strain; and in the final stage, the rate of the ratio increases gradually slow. It is more noteworthy that the energy dissipated by breakage is a small fraction within the total energy input with σc = 2 MPa. In contrast, at higher confining pressures, larger fraction of breakage energy within the total input energy happens during shearing. However, although particle breakage is necessary for sand under high pressures, most of the plastic work is dissipated by friction [52]. As Bolton et al. [7] mentioned, bond breakage opens up new degrees of freedom and allows more energy to be dissipated during consequent rearrangements.

5. Correlation between Particle Breakage and Energy Input per Unit Volume

A large number of experiments have investigated the correlation between particle breakage and total energy input per unit volume [12, 47, 50]. Herein, the correlations between relative particle breakage (Br) and total energy input per unit volume (ET) under different confining pressures are presented in Figure 20. It is notable that the relative particle breakage for all tests increases with the increase of energy input. In the initial stage, particle breakage is very small and it occurs at certain energy, which increases from the beginning of isotropic compression. The result is consistent with the experimental observation by Yu [12], Hyodo et al. [18], and DEM result simulated by Liu et al. [53]. The data distribution may be fitted well by a hyperbolic equation of the following form:where ET is the total energy input including isotropic compression and shearing stages and a and b are the hyperbolic curve-fit parameters. It can be seen that the extent of particle breakage is approximately related to the total energy input, regardless of initial confining pressure.

As described above, the total energy input includes effects of friction and rearrangement between agglomerates, which governs the behavior at lower confining pressures where very little particle breakage exists, as shown by the low slope at the early stage of breakage in Figure 20. Nevertheless, the energy dissipated by breakage begins to play a significant role in the behavior with the increase of confining pressure and shear strain. This is because the agglomerates contact with each other tightly at higher stresses, and the effect of particle breakage on shearing deformation becomes more significant.

6. Conclusion

This paper discussed the results of DEM simulation of macro- and micromechanical behaviors of breakable granular soils under biaxial tests. Granular soils with random Gaussian distribution are modeled as breakable agglomerates bonded by elementary balls. Agglomerates with increased bond strengths were tested to ensure the feasibility and validity of the simulation method, which satisfies Weibull’s statistical distribution well. The presented model has the advantage of less computation time attributed to the simple equation and search procedure for bonds and contacts. Therefore, simulation with an assembly of more than 1700 breakable agglomerates was successfully conducted within the limited domain.

A comparison between the results for assemblies of breakable and unbreakable agglomerates was studied. The difference between the comparison results can be attributed to particle breakage during shearing. The evolution of anisotropy from isotropic compression to the end of shearing was examined, and shear deformation of the assembly produced an evident anisotropy in the distribution and magnitude of normal and tangential forces. Higher confining pressure and larger axial strain were found to have a significant influence on the variation of particle grading, agglomerate numbers, and broken bond increments. The ratio of energy dissipated by breakage (Eb) and total energy input (EB) indicates that energy dissipated by breakage is a small fraction of energy input compared with energy dissipated by friction. A hyperbolic correlation between relative particle breakage and total energy input per unit volume was established regardless of the influence of confining pressure. Overall, the consistency of DEM simulation and previous experimental findings indicate that the simulation method can help us to improve the quantitative knowledge of particle breakage behaviors of granular soils in 2D discrete element analyses.

Data Availability

The data used to support the findings of this study are available from the corresponding author upon request.

Conflicts of Interest

The authors declare that there are no conflicts of interest regarding the publication of this paper.

Acknowledgments

The authors gratefully acknowledge the support received from the National Natural Science Foundation of China (Grant 51804300), the National Key Research and Development Program of China (Grant 2016YFC0600904), and the China Postdoctoral Science Foundation (Grant 2018M632418).