Discrete element simulation of crushable rockfill materials

Lei SHAO*, Shi-chun CHI Liang-jing ZHOU, Yu-zan WANG

1. Faculty of Infrastructure Engineering, Dalian University of Technology, Dalian 116024, P. R. China

2. Hydrochina Chengdu Engineering Corporation, Chengdu 610072, P. R. China

3. Changjiang Institute of Survey, Planning, Design and Research, Changjiang Water Resources Commission, Wuhan 430010, P. R. China

Discrete element simulation of crushable rockfill materials

Lei SHAO*1,2, Shi-chun CHI1, Liang-jing ZHOU3, Yu-zan WANG1

1. Faculty of Infrastructure Engineering, Dalian University of Technology, Dalian 116024, P. R. China

2. Hydrochina Chengdu Engineering Corporation, Chengdu 610072, P. R. China

3. Changjiang Institute of Survey, Planning, Design and Research, Changjiang Water Resources Commission, Wuhan 430010, P. R. China

A discrete element method was used to study the evolution of particle crushing in a rockfill sample subjected to triaxial shear. A simple procedure was developed to generate clusters with arbitrary shapes, which resembled real rockfill particles. A theoretical method was developed to define the failure criterion for an individual particle subjected to an arbitrary set of contact forces. Then, a series of numerical tests of large-scale drained triaxial tests were conducted to simulate the behaviors of the rockfill sample. Finally, we examined the development of micro-characteristics such as particle crushing, contact characteristics, porosity, deformation, movement, and energy dissipation. The simulation results were partially compared with the laboratory experiments, and good agreement was achieved, demonstrating that the particle crushing model proposed can be used to simulate the drained triaxial test of rockfill materials. Based on a comparison of macro behaviors of the rockfill sample and micro structures of the particles, the microscopic mechanism of the rockfill materials subjected to triaxial shear was determined qualitatively. It is shown that the crushing rate, rather than the number of crushed particles, can be used to reflect the relationship between macro- and micro-mechanical characteristics of rockfill materials. These research results further develop our understanding of the deformation mechanism of rockfill materials.

discrete element simulation; triaxial shear test; rockfill material; geometric construction of particles; particle flow code; microscopic mechanism

1 Introduction

Rockfill materials suffer high stress in dam engineering due to the increasing height of rockfill dams. More and more crushed rockfill particles under high compressive pressure result in stress adjustment and location rearrangement inside a dam. Therefore, particle crushing of rockfill materials becomes increasingly significant. Understanding the crushing mechanisms inside the granular material under compression may provide an efficient key to study the properties of rockfill particle crushing. The crushability of rockfill particles has been studied experimentally by some scholars (Guo et al. 1997; Zhang and Si 1982; Fu et al. 2009).In addition to laboratory tests, the discrete element method (DEM) can be used to study their crushability. The DEM combined with the bonded particle method is a scientific tool used to investigate the complex micromechanism that leads to macroscopic behaviors of the granular material. It can also be used for laboratory tests, providing a better visualization of the crushable behavior.

The macroscopic behaviors of the soil mass are affected by the particle shape. It is well known that a rockfill particle has certain characteristics, such as a larger diameter, irregular shape, and larger locked-in force. Therefore, the complex shape of rockfill materials can be considered another issue affecting rockfill behaviors, both in DEM simulations and laboratory experiments (Cho et al. 2006; Potyondy and Cundall 2004; Shao et al. 2009). Many methods have been provided to simulate the shape of particles from simple to complex in numerical simulations using the DEM. For example, a random distribution model for a rockfill sample proposed by Ma et al. (2011) was generated by means of a random simulation technique. Lu and McDowell (2007) used overlapping balls to form clumps by controlling the sphericity, angularity, and surface texture. The particle flow code method allows elementary particles to be bonded to form clusters with arbitrary shapes. Cheng et al. (2003) made the agglomerates from regular crystalline assemblies of balls using PFC3D. The methods proposed by these scholars had their own characteristics, but the algorithms were complex. In this study, a simple procedure was developed to generate clusters resembling real rockfill particles by means of the random simulation technique based on the Mathematica package.

When the numerical sample is generated, it can be used to study the crushing mechanisms of rockfill materials. However, the DEM does not allow particle crushing. Different solutions have been proposed to overcome the limitation of having no particle crushing (Lobo-Guerrero et al. 2006). One solution is to treat each particle as an agglomerate built by bonding uniform or non-uniform smaller particles. In order to establish a failure criterion, the bonding strength between the particles forming agglomerates can be specified (Lobo-Guerrero et al. 2006). This method is capable of simulating particle crushing and providing better interlocking. McDowell and Harireche (2002) and Cheng et al. (2003) applied this method to simulation of the crushing phenomenon of sand particles subjected to compression.

However, these methods do not reflect the effects of particle sizes on particle crushing. Rockfill is a granular material containing many initial cracks and defects. Its tensile strength is much lesser than its shear strength and compressive strength. It has been observed that the particle fracture is always characterized by the propagation of a transversal crack inside the particle. This propagation is mainly caused by the high stress concentration inside the particle, whereas its initiation depends on the nature of the material and the loading conditions. Tsoungui et al. (1999a) proposed a theoretical model to study the crushing mechanisms of particles inside the granular material under the oedometric compression. Lobo-Guerrero and Vallejo (2005) and Lobo-Guerrero et al. (2006) visualized the evolution of particle crushingsubjected to uniaxial compression under direct shear conditions using a simplified failure criterion. However, these studies were just based on a two-dimensional particle flow code of disc particles.

In this study, we used a simplified three-dimensional numerical crushing model for rockfill materials on the basis of a tensile failure criterion and used the sphere cluster as the basic unit to generate an aggregation of particles so as to simulate the large-scale drained triaxial shear test for rockfill. Based on a comparison of the obtained numerical results with those of the experimental tests, we visualized the particle crushing, contact characteristics, porosity, deformation, movement, and energy dissipation.

2 Modeling particle shape by DEM

In this study, a numerical sample of spheres was generated at first according to the gradation curve (Shao et al. 2009), and the information such as the ID number, radius, and central coordinate of each sphere was recorded. Each sphere center was taken as a starting point to generate a regular tetrahedron in its sub-space. Then, irregularly shaped clusters were randomly generated by a breadth-first algorithm based on the cores of these tetrahedrons, with a dense packing of uniform-sized spherical particles bonded together at their contact points. The equivalent radii of clusters were changed according to geometric topological relations and the principle of equal substitution. As these clusters were intended to represent rockfill particles, the numerical sample for rockfill particles was generated.

Locked-in forces can arise if overlaps exist between spheres when bonds are established. A point-contact state between spheres in the cluster must be maintained during its generation, which means that the ratio of the overlap to sphere radius must be less than 10–6. Otherwise, when the maximum tensile stress of a particle exceeds its ultimate tensile stress, a bond will be broken, and the strain energy will be released and transformed into a large amount of kinetic energy, which can lead to violent disturbance of adjacent particles. That will affect the effectiveness of numerical experiments. To avoid this effect, a regular assembly of spheres in the cluster must be ensured at point contacts, without initial overlap. The regular packing may help minimize the space between the spheres and cause the clusters to achieve a dense state (Lobo-Guerrero and Vallejo 2005). This method has the following advantages: random shape of clusters, minimum space between spheres, point contact between the internal spheres of clusters, and representing particle crushing with bond breaking. Fig. 1 shows some typical rockfill particles.

Fig. 1 Typical clusters generated

3 Failure model fulfilling a predefined tensile failure criterion

Usually, the DEM does not allow particle crushing. As mentioned above, a solution has been proposed to overcome the limitation of having no particle crushing when working with DEM codes. In this way, we need to build a numerical sample containing a number of clusters, and treat each cluster as a porous agglomerate generated by bonding uniform smaller particles, following the method described in section 2.

The solution adopted in this study uses a new simplified tensile failure criterion that can be easily implemented with DEM codes to study the evolution of particle crushing. This failure criterion is based on the magnitude of loads applied to particles, the particle coordination number, and the particle size. If particles fulfill a predefined tensile failure criterion, they will be deleted and replaced by a group of new particles following certain rules that agree with the experimental observations.

To implement the solution, we used a theoretical approach first to derive formulae of the maximum tensile stress on an individual sphere subjected to an arbitrary set of contact forces on the basis of the elastic half-space theory. Then, the particle crushing criterion was defined from the knowledge of particle fracture under diametral compression. In the particular case when the fragmentation is caused by diametral compression, it has been shown that the compressive force required to crush a spheroidal particle follows a power law, and thus the critical tensile stress can be obtained. Based on the maximum tensile stress and critical tensile stress, we propose a simple particle crushing model of rockfill materials.

3.1 Derivation of formula for maximum tensile stress

In the case of an individual particle subjected to an arbitrary number of contact forces, the calculation of the stress distribution inside the particle and the prediction of its fracture condition remain a difficult task (Lobo-Guerrero and Vallejo 2005). In this study, we considered a particle with a radius R, which was compressed diametrically by two opposite forces, as illustrated in Fig. 2. It is well known that, inside a particle under compression, the local stress distribution is inhomogeneous: in the inner zones, one or two principal stress components are tractions, while other regions support purely compressive stresses (Tsoungui et al. 1999a). The spatial distribution of these stresses varies with the geometry of the particle. Assuming that the particle was a brittle material, Kienzler and Schmitt (1990) studied the stress field inside an elastic sphere under diametral compression based on three-dimensional finite element calculations. They confirmed that the azimuthal tensile stress (Tsoungui et al. 1999b) was mainly responsible for the particle fracture. Furthermore, in this study we derived the theoretical formulae of the maximum tensile stress of the stress field based on the elastic half-space theory.

Fig. 2Analysis of equivalent process of arbitrary contact forces

Fig. 3Force analysis of hemisphere

According to Hertz (Wang and Wang 1998), under the assumption that the contact surface is formed by the loads on two spheres at their respective centers, the radiusaof the contact surface is related to the radiiR1andR2of the spheres. The maximum contact pressureq0on the contact surface is related to the mean contact pressureqaccording to the following relationship (Wang and Wang 1998):

We obtain a simple relationship ofq0andathrough Eqs. (3) and (4):

It is assumed that the radius R2tends to infinity (R2→∞), that is, k2→0. Eq. (5) can be simplified as follows:

Then, we find, importantly, that the ratio of the maximum contact pressure q0to the radius a of the contact surface have nothing to do with the mode of force action when a?R1and a?R2. The ratio depends only on the radius of the sphere, Poisson’s ratio, and the elastic modulus.

The value of a is obtained as follows:

Substituting Eq (7) into (4), we can obtain the relationship of the mean contact pressure q and the point load P as follows:

It is well known that formulae for solving the stress value in the elastic half-space have been derived by Boussinesq (Wang and Wang 1998). Based on cylindrical coordinates, the formulae for the stress at a depth of z below the center point of the contact surface are derived according to Boussinesq formulae through the application of the superposition principle:

where γ is Poisson’s ratio of the material. Substituting Eqs. (7) and (8) into Eqs. (9) and (10), as z=R1, we can obtain the stress values σz, σr, and σθat the center of the sphere. Because of the symmetry, the stress values under a couple of opposed loads can be written by transformingσz, σr, and σθto rectangular coordinates as follows:

3.2 Critical tensile stress

The particle crushing criterion was defined from the knowledge of the particle fracture under diametral compression. In the particular case when the fragmentation is caused by diametral compression, it has been shown that the compressive forceFcritrequired to crush a spheroidal particle with a radiusRis given by a power law:Fcrit=K0Rα, where the parameterK0can be defined as a pre-factor, which only depends on the material nature (Tsoungui et al. 1999b). Concerning the exponentα, it is dependent on the particle dimension space and Weibull’s modulus of particles. In the case of a brittle material such as ordinary chalk, limestone, cement, or technical ceramics, the Weibull’s moduli lie between 5 and 20, which implies that the value ofαvaries within the range of 1.4 to 1.85 in the threedimensional space (Tsoungui et al. 1999b).

4 DEM simulation

4.1 Laboratory experiment

A series of laboratory tests, including consolidated drained triaxial shear tests and particle analysis, were conducted on rockfill samples of an earth-rockfill dam in China (Gong 2005).

Rockfill materials were taken from a quarry near the dam site. They were composed of the slightly weathered granite. The designed porosity was 0.22, the specific gravity was 2.68, the dry density was 2.09 g/cm3, and the maximum particle size was 800 mm. Considering that the maximum particle size of the sample permissible was 60 mm according to theCode ofWater Injection Test for Water Resources and Hydropower Engineering(the value of the diameter of the experiment instrument divided by the permissible maximum particle size of the sample is greater than or equal to 5), the designed gradation should be scaled. The gradation of the test sample is shown in Table 1, wheredis the particle size.

Table 1Gradation of samples and laboratory test results

The large-scale triaxial apparatus (SJ-70) was used. The diameter of the test sample was 302 mm, and the height was 655 mm. The sample was loaded and tamped into seven levels. Water came into the equipment from the bottom until the sample was saturated under the confining pressure of 20 kPa. The extra water was first drained away from the top of the apparatus. Then, the sample was consolidated under the confining pressures of 0.1 MPa, 0.4 MPa, 0.8 MPa, 1.2 MPa, and 1.6 MPa. The results of laboratory tests after the drained triaxial test are shown in Table 1.

4.2 DEM simulation procedures

Fig. 4 shows the arrangement of the simulated triaxial test for spheres. The numerical sample had a diameter of 302 mm and a height of 655 mm. It was formed by creating 1 062 initial clusters, containing 8 617 small spheres, and 17 522 contact bonds.

Fig. 5 shows the initial gradation of the laboratory test and the numerical test, wherenis the percentage of particles by weight with sizes smaller thand. The samples were loaded by specifying the velocities of the top and bottom walls. The stresses were computed by taking the average forces on walls divided by appropriate areas. Throughout the loading process, the confining stress was kept constant by adjusting the velocity of the radial wall using a numerical servomechanism implemented throughout the DEM procedures. The simple contact bond served to limit the total normal and shear forces. The maximum tensile force that the bond can sustain and the maximum shear force it can withstand before crushing were specified when the bond was created. The bond would break if either of these values was exceeded, meaning that the rockfill particle would partially or completely fail.

Fig. 4Numerical sample

Fig. 5Particle size distribution analytical curves

4.3 Guidelines for choosing macro- and micro-property parameters

The macro-property parameters for rockfill materials include Young’s modulusE, Poisson’s ratioγ, the densityρ, pre-factorK0, and exponentα. The values of these macro-property parameters were determined by experiments on the rockfill property, which were 7.0 GPa forE, 0.25 forγ, 2.68 g/cm3forρ, 7.52 × 106forK0, and 1.4 forαin this study.

The micro-property parameters for particles include the normal contact stiffnesskn, the shear contact stiffnessks, and the friction coefficientμ. It is possible to choose micro-property parameters for particles leading to macroscopic behaviors which resemble the results of laboratory tests. Considering that the initial Young’s modulus of the material is linearly related to the value of the contact stiffness, the normal stiffness of numerical rockfill particles (kn) and those of the top and bottom walls (knw1) were set to the same values. The radial wall was given a stiffness (knw2) one-tenth of the particle’s stiffness (kn), in order to simulate a soft confinement. Considering that Poisson’s ratioγdepends on the ratio of the shear contactstiffness to normal contact stiffness, the ratio was set to 0.25. The peak strength of the material depends both on the friction coefficient (μ), which causes the material to exhibit a plastic or weakly softening behavior, and the dilatancy effect (Shao et al. 2009). Therefore, the friction coefficient was set to 0.38 during subsequent loading. Meanwhile, the bond strength affected the characteristics of particle crushing during loading.

During the numerical triaxial shear test, the micro-property parameters were chosen by comparing the experimental results with the numerical ones, which were determined to be 8.2 × 107N/m forknandknw1, 8.2 × 106N/m forknw2, and 2.05 × 107N/m forks. The shear contact stiffness of the top and bottom walls (ksw1) and that of the radial wall (ksw2) were determined to be zero.

4.4 Modeling results

The results of laboratory triaxial shear tests on the real rockfill samples and DEM simulations of clusters were compared. Fig. 6 shows the results of laboratory triaxial shear tests under different confining pressures (0.8, 1.2 and 1.6 MPa) and those obtained with the parameters above. It can be seen from the results of laboratory experiments and DEM simulations that the stress-strain behaviors of rockfill samples are highly nonlinear, and that the difference between the test results of the peak stress and those calculated with the DEM increases with the confining pressure. This means that the numerical sample has a confining pressure-dependent manner during loading. This is because the number of numerical particles in the test is far less than the number in reality, which leads to insufficient compactness and the occurrence of locked-in forces. On the basis of comprehensive comparison and analysis of the tests curves, we can conclude that the numerical tests can reflect the stress-strain behaviors of the rockfill materials.

Fig. 6Comparison of results from laboratory tests and numerical simulations

5 Evolution of microscopic characteristics

5.1 Evolution of particle crushing

Fig. 7Variations of particle crushing rate with axial strain under different confining pressures

Fig. 8Distribution of crushed particles with different axial strains under confining pressure of 1.2 MPa

The analysis above shows that the particle crushing mainly occurs on a relatively concentrated axial strain interval (about 5% to 8%) under different confining pressures. It isnoted that under a relatively higher confining pressure, where dilatancy is considerably higher (Fig. 6), the particle crushing starts at a smaller axial strain. According to the analysis on the evolution of the particle crushing rate, the development process of the rockfill particle crushing under triaxial conditions can be divided into four stages: in the first stage, rockfill particles are not crushed during the initial loading; in the second stage, a minimum amount of particles are crushed; in the third stage, the crushing shows an accelerating tendency; and in the fourth stage, a stable increase of the crushing rate can be observed. This phenomenon can verify the extrapolation of the particle crushing process in macroscopic laboratory experiments.

5.2 Development of contact characteristics

Figs. 9 and 10 show the development of porosity and contact forces under the confining pressure of 1.2 MPa. Through the comprehensive comparison of developments of stress-strain behaviors, particle crushing, contact forces, and porosity, we can obtain some new results.

Fig. 9Development of porosity under confining pressure of 1.2 MPa

Fig. 10Development of contact force under confining pressure of 1.2 MPa

It can be observed that during the first stage of loading, no particles are crushed. The maximum contact force rapidly increases to its peak, and the porosity decreases. Meanwhile, the volume contraction becomes prevalent. The rockfill behaviors in this stage are influenced by interspace filling, location adjustment, and rearrangement among particles. In the second stage, a small amount of particles are crushed. The maximum contact force first increases and then decreases. Meanwhile, the average contact force tends toward stability.

In the third stage, it can be seen that a large amount of particles are quickly crushed (Fig. 7). These behaviors are explained by the fact that the particles located in the upper and lower regions of the numerical sample have to bear larger loads, and their coordination numbers are low since they are closer to the walls.

In the fourth stage, the crushed particles maintain stable growth. Maximum contact force within the particles further decreases, and the average contact force slightly decreases, as shown in Fig. 10. Meanwhile, the porosity tends to increase in the middle region, and the sample is in volume expansion, as observed from volumetric strain-axial strain curves in Fig. 6. After the axial strain exceeds 12%, the porosity of the middle region keeps increasing, leading to a critical failure condition of the rockfill sample.

Correlation analysis indicates that an important factor affecting microscopic contact characteristics, porosity, and the macroscopic stress-strain curve is the crushing rate of particles. Its effect is much greater than that of the number of crushed particles.

5.3 Movement of particles

Fig. 11 shows that particles in the upper and lower regions of the sample (region 1) near the loading device move along the loading direction. The particles in the middle region of the sample (region 2) move outside in the direction of 30° to 45° from the horizontal plane. When the axial strain is about 2%, the movement shows apparent nonuniformity and the particles in region 2 move slowly. With the strain increasing, the density of the vector diagram of particle movement increases gradually, and the move tends to be faster, showing that the sample transforms gradually from volume contraction to volume expansion. After the strain increases to 10%, the particles in region 2 move quickly toward the outside, leading to structural changes in the sample.

Fig. 11Vector diagram of particle movement

5.4 Energy conversion

The transformation and dissipation of the elastic strain energy, frictional energy, and crushing energy within granular systems always occur during the triaxial shear test. Fig. 12 shows the development processes of the elastic strain and frictional energies against the axial strain, respectively. The trends of the curves agree with the results of Ma et al. (2011). Based on Fig. 12, we try to study the relationship between the energy dissipation and behavior of the rockfill sample during the four stages of particle crushing.

Fig. 12Energy dissipation of rockfill sample under confining pressure of 1.2 MPa

In the first stage, the elastic strain energy increases rapidly, with less dissipation of the frictional energy and without occurrence of particle crushing. In the second stage, the frictional and crushing energy dissipations increase, and the friction energy dissipation grows faster than the crushing energy dissipation. In this stage, the friction is the main factor influencing the behaviors of the particles. In the third stage, the frictional and crushing energy dissipations grow rapidly. Meanwhile, the effect of the elastic strain energy can be ignored in comparison with the friction and crushing energy dissipations. Pores in the granular sample are filled with a large number of smaller cracked particles, leading to an increase of the specific surface area of the sample. In the fourth stage, the developments of frictional and crushing energy dissipations are relatively stable.

5.5 Preliminary discussion on microscopic mechanism

Now, we present our findings on the relationship between the behaviors of the rockfill sample and the microscopic characteristics of friction, crushing, deformation, and movement of rockfill particles based on the analysis above during the four stages of the crushing process.

Throughout the loading process, the particles first compress each other, leading to a rapid increase of contact force and elastic strain energy. With the continuous loading, the rockfill sample tends to compact due to the voids among particles to be filled. After the whole sample is enhanced enough to reach a certain level of compactness, the growth of the strain energy slows down, and a small amount of particles are crushed. The location adjustment and rearrangement among particles develop sustainably, and the maximum contact force decreases. Meanwhile, the volume of the sample tends to dilatate. This analysis verifies the conclusion that the friction becomes a primary factor influencing the behaviors of the particles. When loaded into the third stage, the sample tends toward its densest state, and particle crushing shows a tendency of accelerated increasing. Meanwhile, the proportion of crushing energy dissipation gradually catches up with the frictional energy dissipation. After that, particle crushing increases at a steady rate. The main energy dissipations of rockfill during this stage are frictional and crushing energy dissipations. With the continuous loading, more and more particles are crushed at a stable high rate, and particles cross and roll over each other, causing structural changes in the sample and a decrease of the carrying capacity of rockfill, indicating that the sample reaches the critical failure status, and finally the rockfill sample tends to fail.

6 Conclusions

The study of stress and deformation characteristics of rockfill materials is necessary due to the frequent occurrences of geological disasters. Through microscopic visualization, a series of numerical tests based on the discrete element method were performed to study the mechanism of crushable rockfill materials. Some conclusions are as follows:

(1) The simulated results agree with the experimental results. It can be seen that theproposed tensile failure criterion is appropriate.

(2) The numerical tests reveal the evolution of microscopic characteristics of the rockfill sample such as the crushing rate, contact force, porosity, particle movement, location of crushed particles, and energy dissipation during the loading stages.

(3) The middle particles’ movement outward along the radial direction is a main factor leading to dilatancy. The development of the porosity in different regions confirms these observations.

(4) Correlation analysis indicates that an important factor affecting macroscopic and microscopic mechanisms is the crushing rate of particles. It is more important than the number of crushed particles.

(5) This work can offer some useful ideas: To study the properties of rockfill materials, from the view of numerical simulation, the irregular shape of rockfill particles and particle crushing characteristics must be considered. On the basis of the results of numerical experiments, the evolution of macroscopic mechanical behavior of rockfill in triaxial shear conditions is divided into four stages by analysis of the variations of the macrostructures and micro-parameters. It is worth further study to better understand the existence of the threshold for the four stages and the factors affecting the threshold in the form of the microscopic parameters and the external environment.

Cheng, Y. P., Nakata, Y., and Bolton, M. D. 2003. Discrete element simulation of crushable soil. Geotechnique, 53(7), 633-641. [doi:10.1680/geot.53.7.633.37389]

Cho, G. C., Dodds, J., and Santamarina, J. C. 2006. Particle shape effects on packing density, stiffness, and strength: Natural and crushed sands. Journal of Geotechnical and Geoenvironmental Engineering, 132(5), 591-602. [doi:10.1061/(ASCE)1090-0241(2006)132:5(591)]

Fu, H., Lin, H., and Cai, Z. Y. 2009. Influencing factors for particle breakage of coarse soil. Journal of Hohai University (Natural Sciences), 37(1), 75-79. (in Chinese)

Gong, W. Q. 2005. The Feasibility Study Stage of Sujiahekou Hydropower Station: The Test Study Report of Weathered Materials and Rock-fill in Location. Kunming: The Scientific Research Branch of Kunming Investigation and Design Institute, Hydrochina Kunming Engineering Corporation. (in Chinese).

Guo, X. L., Hu, H., and Bao, C. G. 1997. Experimental studies of the effects of grain breakage on the dilatancy and shear strength of rockfill. Chinese Journal of Geotechnical Engineering, 3(19), 83-88. (in Chinese)

Kienzler, P., and Schmitt, W. 1990. On single-particle comminution-numerical analysis of compressed spheres. Powder Technology, 61, 29-38. [doi:10.1016/0032-5910(90)80063-5]

Lobo-Guerrero, S., and Vallejo, L. E. 2005. Discrete element method evaluation of granular crushing under direct shear conditions. Journal of Geotechnical and Environmental Engineering, 131(10), 1295-1300. [doi:10.1061/(ASCE)1090-0241(2005)131:10(1295)]

Lobo-Guerrero, S., Vallejo, L. E., and 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)]

Lu, M., and McDowell, G. R. 2007. The importance of modelling ballast particle shape in the discrete element method. Granular Matter, (9), 69-80. [doi:10.1007/s10035-006-0021-3]

Ma, G., Zhou, W., Chang, X. L., and Zhou, C. B. 2011. 3D mesoscopic numerical simulation of triaxial sheartests for rockfill. Chinese Journal of Geotechnical Engineering, 33(5), 746-753. (in Chinese)

McDowell, G. R., and Harireche, O. 2002. Discrete element modeling of yielding and normal compression of sand. Geotechnique, 52(4), 299-304. [doi:10.1680/geot.52.4.299.41018]

Potyondy, D. O., and 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]

Shao, L., Chi, S. C., and Jia, Y. F. 2009. Meso-mechanical simulation of a large scale triaxial test of rockfill materials. Rock and Soil Mechanics, 30(s1), 239-243. (in Chinese)

Tsoungui, O., Vallet, D., and Charmet, J. C. 1999a. 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]

Tsoungui, O., Vallet, D., Charmet, J. C., and Roux, S. 1999b. Size effects in single grains fragmentation. Granular Matter, 2(1), 19-27. [doi:10.1007/s100350050030]

Wang, G. W., and Wang, Y. C. 1998. Elasticity. Shanghai: Shanghai Jiaotong University Press. (in Chinese)

Zhang, Q. Y., and Si, H. Y. 1982. Shear strength and stress-strain property of coarse grain soil determined by large scale triaxial compression tests. Journal of Hydraulic Engineering, (9), 22-31. (in Chinese)

(Edited by Ye SHI)

This work was supported by the National Key Basic Research Program of China (Grants No. 50879007 and 50979014), and the Specialized Research Fund for the Doctoral Program of Higher Education of China (Grant No. 20090041110016).

*Corresponding author (e-mail: sl-fly@163.com)

Received Feb. 14, 2012; accepted Jul. 26, 2012