Modeling the particle breakage of rockfill materials with the cohesive crack model

Computers and Geotechnics 61 (2014) 132–143

Contents lists available at ScienceDirect

Computers and Geotechnics journal homepage: www.elsevier.com/locate/compgeo

Modeling the particle breakage of rockfill materials with the cohesive crack model Gang Ma, Wei Zhou ⇑, Xiao-Lin Chang State Key Laboratory of Water Resources and Hydropower Engineering Science, Wuhan University, Wuhan 430072, China

a r t i c l e

i n f o

Article history: Received 19 December 2013 Received in revised form 7 May 2014 Accepted 10 May 2014

Keywords: Rockfill materials Particle breakage Combined FDEM Cohesive crack model Force chain

a b s t r a c t A practical combined finite–discrete element method was developed to simulate the breakage of irregularly shaped particles in granular geomaterials, e.g., rockfill. Using this method, each particle is discretized into a finite element mesh. The potential fracture paths are represented by pre-inserted cohesive interface elements (CIEs) with a progressive damage model. The Mohr–Coulomb model with a tension cut-off is employed as the damage initiation criterion to rupture the predominant failure mode occurs at the particle scale. Two series of biaxial tests were simulated for both the breakable and unbreakable particle assemblies. The two assemblies have identical configurations, with the exception that the former is inserted with CIEs and is breakable. The simulated stress–strain–dilation responses obtained for both assemblies are in agreement with experimental observations. We present a comprehensive study of the role of particle breakage on the mechanical behavior of rockfill materials at both the macroscopic and microscopic scales. The underlying mechanism of particle breakage can be explained by the force chain in the assemblies. Ó 2014 Elsevier Ltd. All rights reserved.

1. Introduction Rockfill material is increasingly used in many geotechnical engineering applications, particularly in high rockfill dams. In contrast to sand and gravel, rockfill grains generally possess a lower crushing strength and are subjected to a higher contact force due to their larger size and coarse gradation. The role of particle breakage in the mechanical behavior of granular geomaterials has been experimentally investigated by many researchers [1–4]. In parallel with experimental studies, another contribution to this field comes from the numerical modeling of granular materials using the discrete element method (DEM). Previously published works regarding DEM and particle breakage were mainly concerned with the breakage of subparticles that are joined by bonding or cohesive forces [5–11]. Another approach is to replace the particle fulfilling a predefined failure criterion with an equivalent group of smaller particles [12–15]. Both techniques for a particle breakage simulation are mainly employed, with either discs in 2D or spheres in 3D. For particles with general shapes, ⇑ Corresponding author. Address: State Key Laboratory of Water Resources and Hydropower Engineering Science, Wuhan University, 8 Donghu South Road, Wuhan 430072, China. Tel.: +86 027 68773778. E-mail addresses: [email protected] (W. Zhou), [email protected] (X.-L. Chang). http://dx.doi.org/10.1016/j.compgeo.2014.05.006 0266-352X/Ó 2014 Elsevier Ltd. All rights reserved.

the realization of particle breakage follows the idea proposed by Potapov and Cambell [16]. In their simulations, a breakable particle is created by bonding unbreakable and nondeformable subparticles via cohesive bonds; if the bond between these subparticles breaks, breakage occurs [17,18]. A natural evolution of the DEM is to accurately replicate the deformability and crushability of real granular materials and to provide a robust model for contact interaction that is insensitive to the complexity of the considered particle shape. For these reasons, Munjiza et al. proposed a method that merges the capabilities of the finite element method (FEM) with discrete element implementations [19]. In the combined FDEM, each particle is discretized into a finite element mesh, and the deformation is solved computationally according to the rheological behavior of the material, where the particle movements and interactions are handled as they would be in the DEM [20]. One of the remarkable benefits of this combination is its ability to handle the contact interactions of particles with arbitrary shapes, independent of the contact configuration. This combined FDEM has been applied to granular materials by several researchers [21–25]. In addition to the combined FDEM, the multi-particle finite element method (MPFEM) and meshed discrete element method (MDEM) also allow for finite element simulation at the particle scale [26,27]. Basically, the combined FDEM, MDEM, and MPFEM have similar theoretical backgrounds.

G. Ma et al. / Computers and Geotechnics 61 (2014) 132–143

The main purpose of this paper is to develop a practical combined FE/DE method for modeling particle breakage in granular geomaterials, e.g., rockfill materials. Our method uses the cohesive interface element (CIE) in the general-purpose finite element software ABAQUS to model potential particle fracture or fragmentation. First, densely packed polygonal particles are prepared to represent the rockfill sample. The CIEs are inserted into the initial FE mesh of each particle using a specially designed and efficient algorithm. Then, the cohesive crack model and a progressive damage model are introduced. Once the initial damage criterion is satisfied, the CIE starts to degrade and gradually loses its carrying capability. Finally, two series of biaxial tests are simulated with both the breakable and unbreakable particle assemblies. The input parameters of these simulations are chosen on a purely empirical basis and without calibration to any specific type of parent rock. The simulation results are analyzed on both the macro- and micro-scales to investigate the role of particle breakage in the mechanical behavior of rockfill materials.

2. Combined finite–discrete element method In the combined approach, each particle is discretized into an FE mesh, as shown in Fig. 1. These meshes define the shapes of the discrete elements. The contact between the interacting particles is defined commonly in DEM. Once the contacts have been detected, a linear contact model is employed to calculate both the normal and tangential contact forces between two particles. The normal component is linear-elastic with a no-tension limit, and the shear component is linear-elastic with a Coulomb friction limit. In this paper, the combined FDEM-based simulations of granular materials are performed using the explicit module of ABAQUS [28]. The explicit integration scheme and general contact capability of this module make it appropriate for a large number of actual and potential contacts for particles undergoing large deformation. In addition, the ABAQUS software provides the user with an extensive array of secondary development interfaces that allow them to adapt ABAQUS to their particular analysis requirements. A user-defined CIE with the theoretical formulation described below was developed and implemented into ABAQUS using the user defined element subroutine. Several python scripts were also developed to facilitate the pre- and post-processing processes.

Fig. 1. Typical particle assembly used in combined FDEM simulations.

133

3. Particle breakage with the cohesive crack model 3.1. Modeling procedures The proposed method includes the following steps:  Generating a realistic description of the particle assembly for discrete modeling. A Voronoi tessellation-based approach is adopted to generate the 2D assembly of polygonal particles [29]. Because the simulated particle breakage patterns are dependent upon the initial FE mesh of particles, the selection of the element type is relevant to the accuracy of the simulation. Compared with the use of quadrilateral elements, the simulated breakage patterns of an FE mesh with triangular elements have less mesh sensitivity. However, the locking problem arising from linear triangular elements reduces the simulation accuracy. In this paper, quadrilateral elements are preferred in a fine FE mesh so that the particle breakage can be simulated with adequate accuracy and relatively less mesh sensitivity.  Compressing the initially loose particle assembly until the desired void ratio is reached. This process is accomplished using a combined FDEM simulation.  Inserting CIEs into the initial FE mesh of the dense particle assembly using a python-based computer program. This paper uses four-node CIEs with zero in-plane thickness. The FE mesh of a single particle is taken as an example to clarify the method used to realize this procedure (Fig. 2). The details of the cohesive element insertion algorithm can be found in the study by Yang et al. [30].  Numerical simulations of biaxial tests of the polygon assemblies are performed using the ABAQUS/Explicit solver. The detailed modeling procedures will be discussed later. 3.2. Cohesive crack model The cohesive crack model assumes the existence of a fracture process zone (FPZ) in front of the crack tip, in which energy dissipation occurs during fracture, as illustrated in Fig. 3. In the FPZ for a 2D case, tractions exist in the normal direction tn and shear direction ts across the crack surface, and the corresponding relative displacements are dn and ds. The nonlinear relationship between relative displacement and traction inside the FPZ is characterized by the cohesive crack model. Fig. 4 shows a typical bilinear traction–separation response. A linear elastic traction–separation law is assumed to model the initially undamaged material prior to damage. Once damage is initiated, the failure of the elements is characterized by the progressive degradation of the material stiffness, which is driven by a damage process. The initial CIE stiffness is defined by the Young’s modulus and elastic shear modulus of the intact material. The initial normal 0 0 stiffness kn and shear stiffness ks should be large enough to

Fig. 2. Inserting CIEs in the initial mesh (which is exaggerated for clarity).

134

G. Ma et al. / Computers and Geotechnics 61 (2014) 132–143

Fig. 3. Characterization of the cohesive crack model.

The fracture energy under mixed-mode loading is GT ¼ GcI þ Gshear , where Gshear ¼ GcII þ GcIII is the fracture energy for shear loading proposed by Li [34]. The Benzeggagh–Kenane fracture criterion is practically useful when the critical fracture energies during deformation along the first and second shear directions are the same, i.e., GcII ¼ GcIII . Taking the linear damage evolution (in Fig. 2) as an example, the damage variable, D, is given by

D¼ dfm Fig. 4. Constitutive relations of the cohesive interface elements.

produce a negligible relative displacement before the shear slip or tensile debonding occurs but should not be high enough to cause a numerical ill-condition [31]. In the case of a zero-thickness CIE, the initial stiffness is given as follows after several trials: 0

kn ¼

10E 0 10G 10E ; ks ¼ ¼ le le 2le ð1 þ uÞ

ð1Þ

where le is the average element size for the whole mesh and E, G, and u are the Young’s modulus, shear modulus, and Poisson’s ratio of the surrounding solid elements, respectively. The CIEs presented are based upon the cohesive crack model. The resilient feature of the CIEs is based upon the damage mechanics framework, in which the stiffness values, kn and ks, upon unloading and reloading are degraded as a result of the progressive irreversible damage. The damage is characterized by a scalar variable, D, which represents the overall damage of the crack under mixed-mode loading. The evolution of the damage under a combination of normal and shear deformations across the interface is defined by the effective displacement dm and effective traction teff:

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi hdn i2 þ d2s qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ¼ htn i2 þ t2s 

hdn i ¼

ð2Þ

dn ; dn  0 ðtensionÞ 0;

dn < 0 ðcompressionÞ

G ¼

GcI

þ

ðGcII



GcI Þ

c

¼ 2G

ð5Þ

=2Gc t 0eff

where dmax refers to the maximum value of the effective relative m displacement attained during the loading history. d0m and dfm are the effective relative displacements at damage initiation and complete failure, respectively. t 0eff is the effective traction at damage initiation. The damage variable monotonically increases from 0 to 1 upon further loading, after the initiation of damage. The damage initiation and evolution will degrade the unloading and reloading stiffness values, kn and ks: 0

kn ¼ ð1  DÞkn

ð6Þ

0

ks ¼ ð1  DÞks

The corresponding tractions are

( tn ¼

0

ð1  DÞkn dn

dn  0 ðtensionÞ

0 kn

dn < 0 ðno damage to compressive stiffnessÞ 0

t s ¼ ð1  DÞks ds ð7Þ

 g Gshear GT

Shearing damage Zone/Mode II

ð3Þ

where h i is the Macaulay bracket. This definition is based upon the fact that compression is not a main source of mechanical damage. The most widely used criterion to predict damage propagation under mixed-mode loading is the expression proposed by Benzeggagh and Kenane [32,33]. The Benzeggagh–Kenane fracture criterion for the critical energy Gc is c

f 0 dmax m ðdm  dm Þ

ts

dm ¼ t eff

  0 dfm dmax m  dm

ð4Þ

Region of admissible stresses Tensile damage zone/Mode I

c

ft

tn

Fig. 5. Damage initiation criteria based on the Mohr–Coulomb model with a tension cut-off.

G. Ma et al. / Computers and Geotechnics 61 (2014) 132–143

The damage initiation criterion used here is the Mohr–Coulomb criterion, with a tension cut-off (Fig. 5). For a given CIE, the maximum normal traction is limited by the tensile strength, ft. The maximum shear traction is characterized by the normal traction, tn, cohesion, c, and internal friction angle, u. Damage is assumed to initiate when the following expression is satisfied.

F t ¼ t n  ft

ð8Þ

F s ¼ jt s j þ tn tan u  c

The input parameters for the CIEs include the friction, cohesion, normal stiffness, shear stiffness, and tensile strength. Once the overall damage variable reaches its upper bound at all of its material points and none of its material points are in compression, the fully damaged CIE is removed from the FE mesh.

Fig. 6. Schematic diagram of the biaxial test and numerical sample.

135

4. Procedures and input parameters of numerical modeling for a biaxial test A 600 mm  1200 mm numerical sample is enclosed by four rigid, frictionless walls and is composed of a polydisperse assemblage of 3676 irregular polygons (Fig. 6). The particle diameter varies between dmin ¼ 10mm and dmax ¼ 30mm. The total area of particles with a set diameter, d, is uniformly distributed within the range of dmin  dmax (see Fig. 7a). The size distribution obtained for the particle assembly is shown in Fig. 7b. The void ratio of the numerical sample is 0.175. Initially, the particles are randomly generated in a rectangle that is large enough to have a solid fraction of 0.5. This is followed by isotropic compression until the desired void ratio is reached. During sample preparation, both the gravitational constant and inter-particle friction coefficient are set to zero, and the particle breakage is also disabled. At the end of sample preparation, both the contact orientation and normal contact force are isotropically distributed, as confirmed by the observations of corresponding rose diagrams. A numerical biaxial test is performed using the following protocol. Following sample generation, the numerical sample is subjected to hydrostatic compaction under gravity-free conditions. Then, the sample is shared by moving the top and bottom boundaries toward each other under strain-controlled conditions, with constant confining pressure acting on the lateral walls. The loading speed is slow enough to ensure that the test is conducted under quasi-static conditions. Mass scaling is used for computational efficiency in quasi-static simulation of the biaxial tests. The amount of mass scaling is well selected to ensure that changes in the mass and consequent increases in the inertial forces do not alter the solution significantly. The combined FDEM simulations require various parameters, e.g., the normal stiffness, shear stiffness, inter-particle friction coefficient, particle strength, elastic modulus, and so on. The selection of these parameters is very important to ensure the accuracy of the simulation. In this paper, however, the values for these parameters are chosen on an empirical basis because the paper is intended to investigate the mechanical behavior qualitatively rather than quantitatively. The micro-parameters used in this paper are summarized in Table 1. Under biaxial conditions with vertical compression, r1 P r2, where r1 and r2 are the principal stress values. The confining pressure is measured on the horizontal walls, and the axial stress is calculated from the static equilibrium of the loading boundaries by computing the resultant of the vertical forces acting on them.

Fig. 7. Particle size distribution.

136

G. Ma et al. / Computers and Geotechnics 61 (2014) 132–143

were conducted with four levels of confining pressures: 0.5, 1.0, 1.5, and 2.0 MPa.

Table 1 Set of parameters used in simulations. Parameter

Value

Basic parameters Mass density (kg/m3) Elastic modulus (GPa) Poisson’ s ratio

2600 20 0.2

Parameters of the cohesive crack model Normal stiffness of CIE (N/m3) Shear stiffness of CIE (N/m3) Tensile strength (MPa) Cohesion intercept (MPa) Friction angle of cohesive element Mode-I fracture energy (N/m) Mode-II fracture energy (N/m)

4.0  1012 1.667  1012 15 34.97 40° 200 1000

Parameters of the contact model Inter-particle friction coefficient Normal contact stiffness Kn (N/m3) Stiffness ratio (Kt/Kn) Particle/wall friction coefficient

0.5 2.0  109 0.25 0

5.1. Macroscopic results

The mean and deviator stresses are expressed as p ¼ ðr1 þ r2 Þ=2 and q = r1  r2, respectively. The cumulative vertical, horizontal, and shear strains, e1, e2, and eq, respectively, are defined by

e1 ¼

Dh h0

ð9Þ

e2 ¼

Dw w0

ð10Þ

eq ¼ e1  e2

ð11Þ

where h0 is the initial height and Dh = h0  h is the total downward displacement. w0 is the initial width and Dw = w0  w is the total change of the box width. The cumulative volumetric strain ev is then defined and given by

ev ¼

DV V0

ð12Þ

where V0 = w0h0 is the initial volume and DV = V0  V is the cumulative change of the solid fraction. 5. Numerical simulation results To investigate the role of particle breakage on the mechanical behavior of rockfill materials, two series of numerical biaxial tests

The stress–strain curves derived from the numerical biaxial tests are characterized by several macroscopic parameters. The parameter Ei corresponds to the initial tangential slope of the stress–strain curve. The peak angles of shearing resistance, umax, and dilation, wmax, are defined and given by

ðr1 =r2 Þmax  1 ðr1 =r2 Þmax þ 1 ðde1 =de2 Þmax þ 1 sinðwmax Þ ¼  ðde1 =de2 Þmax  1 sinðumax Þ ¼

ð13Þ

Figs. 8 and 9 illustrate the simulation results in terms of the deviator stress and volumetric strain versus the axial strain. In both tests, increasing the confining pressure resulted in a greater deviator stress and initial tangential modulus, in line with experimental observations (Figs. 8a and 9a). For low confining pressures, failure occurs at a relatively low strain level, and there is a monotonic evolution in the response as the confining pressure increases. The volumetric strain versus axial strain clearly shows that during the initial shearing phase, the assembly has a contractive behavior and a limited volume compression occurs. As the shearing loading proceeds, the dilation gradually occupies a dominant position, which leads to a volume expansion (Figs. 8b and 9b). Under higher confining pressures, the assembly becomes more compressible and indicates that it has a dilating tendency at higher strains. However, under lower confining pressures, the assembly has a dominant dilative behavior. Fig. 10 shows a direct comparison of the simulation results for both the breakable and unbreakable systems. The figure shows qualitatively that the simulated stress–strain–dilation responses obtained for both assemblies are typical of those observed in laboratory tests. The initial tangential modulus and shear strength are much higher for the unbreakable assembly, which exhibits an obvious peak at an axial strain of approximately 4%, followed by a strain softening behavior extending to the residual state. Furthermore, large dilation is observed by an initial slight compression followed by a significant volume expansion. This appreciable increase in the peak shear resistance is primarily caused by the strong interlocking of unbreakable particles and, consequently, an increase in the rate of dilation. However, after the peak occurs, the interlocking effect disappears and the shear resistance drops considerably. In contrast, the breakable particles cannot bear the force imposed

Fig. 8. Assembly of unbreakable particles: the relationships between deviator stress (a) and volumetric strain (b) and axial strain.

G. Ma et al. / Computers and Geotechnics 61 (2014) 132–143

137

Fig. 9. Assembly of breakable particles: the relationships between deviator stress (a) and volumetric strain (b) and axial strain.

Fig. 10. Simulation results for two assemblies under a confining pressure of 2.0 MPa: (a) deviator stress versus axial strain and (b) volumetric strain versus shear strain.

on them and subsequently break. This causes the assembly to have a more contractive behavior and mild post-peak strain softening. At large strain levels, both assemblies deform at a constant volume, which is accompanied by a relatively constant deviator stress, independent of the crushability of the particle assemblies. A conventional explanation is that the shear behavior of a dense granular packing is essentially controlled by the competition between two opposing mechanisms: shear-induced dilation and particle breakage [10]. Because of the extra energy required for the particles to dilate against the confining pressure, a dilative assembly has greater strength than a contractive assembly. Furthermore, the dilation of the assembly is accompanied by a volume expansion, and both the shear stress and dilation rate decrease until they reach critical values. In contrast, particle breakage has an opposing effect on the dilative mechanism and suppresses the mobilization of the assembly dilation. Fig. 11 shows the variation in the peak friction angle as well as the maximum angle of dilation, with the corresponding confining pressures for both series of simulations, where uu is the interparticle friction angle. The peak friction angle and dilation angle are calculated using Eq. (13). It is important to observe that almost the same peak friction angle is obtained for the unbreakable assembly with different confining pressures. In tests with breakable assemblies, the frictional angles tend to decrease with an increase in the confining pressure. This observation leads to the important conclusion that the stress dependence of the peak friction angle of granular geomaterials is primarily caused by particle crushability. Other features are captured by the curves:

Fig. 11. Variation in the peak friction angle (umax), maximum dilatancy angle (wmax), and interparticle friction angle (uu) with changing confining pressure.

 Regarding the dilation characteristics, it is clear that the maximum angle of dilation decreases with an increase in the confining pressure for a breakable assembly. In the case of an unbreakable assembly, however, it dilates at a relatively constant rate under different confining pressures.

138

G. Ma et al. / Computers and Geotechnics 61 (2014) 132–143

 Both the peak friction angle and maximum dilation angle converge to the same value at an extremely low confining pressure, where particle breakage is largely absent.

F ij ¼

5.2. Micromechanical analysis

PðhÞ ¼

Nc ðhÞ Nc

ð14Þ

where Nc is the total number of contacts. In a similar way, by decomposing the contact force into the normal and tangential components in a local contact frame (n, t), where t is an orthonormal unit vector along the tangential direction, we can compute the angular averages, f n ðhÞ and f t ðhÞ, as follows:

f n ðhÞ ¼

1 X c f Nc ðhÞ c2SðhÞ n

f t ðhÞ ¼

1 X c f Nc ðhÞ c2SðhÞ t

ð15Þ

where fnc and ftc are, respectively, the normal and tangential components of the contact force at contact c. These directional functions can be approximated by their second-order Fourier expansions [35,36].

1 ½1 þ ac cos 2ðh  hc Þ 2p f n ðhÞ ¼ f 0 ½1 þ an cos 2ðh  hn Þ f t ðhÞ ¼ f 0 at sin 2ðh  ht Þ P ðhÞ ¼

ð16Þ

where f 0 is the average contact force; ac, an, and at represent the anisotropies of the contact orientations and the normal and tangential contact forces, respectively; and hc, hn, and ht are the corresponding privileged directions. Although the anisotropies ac, an, and at can be calculated by fitting the measured values of PðhÞ, f n ðhÞ, and f t ðhÞ, respectively, to the approximations presented in Eq. (16), it is more convenient to calculate these parameters using fabric tensor and force tensor, given by [37,38].

Fig. 12. Geometry of the contact between two polygons.

PðhÞni ðhÞnj ðhÞ dh ¼

h

Hnij ¼

In this section, we analyze the particle-scale information of two packings of breakable and unbreakable particles in terms of the fabric and force anisotropies and contact force distributions. For an assembly of polydisperse polygonal particles, the fabric and contact network can be described using the directional distributions of the contact orientations, branch vectors, and contact normal and tangential forces. This will allow us to quantify the effect of particle breakage and the connection between particle breakage and macroscopic shear strength. For the 2D case, the unit vector, n, is described using a single angle, h (see Fig. 12). Let S(h) be the set of contacts with direction h 2 ½h  dh=2; h þ dh=2 for the angle increment dh and for which N c ðhÞ is the number of contacts in S(h). The probability distribution function, P(h), of the contact orientations is given by

Z

Htij ¼

Z

1 X ni nj Nc c2Nc

h

X fn ni nj f n ðhÞn ðhÞn ðhÞ dh ¼ 1 i j Nc c2Nc 1 þ ackl nk nl

h

X ft t i nj f ðhÞn ðhÞt ðhÞ dh ¼ 1 t i j Nc c2Nc 1 þ ackl nk nl

Z

ð17Þ

where the subscripts i and j represent the Cartesian components and ni denotes the component of the unit contact normal in direction i for contact c. acij ¼ 4F 0ij , F 0ij is the deviatoric part of fabric tensor Fij. Using the above equations, the anisotropies ac, an, and at can be obtained using the following relations:

ac ¼ 2

F1  F2 ; trðFÞ

an ¼ 2

Hn1  Hn2 ; trðHn Þ

at ¼ 2

Ht1  Ht2 trðHt Þ

ð18Þ

where the subscripts 1 and 2 refer to the principal values of each tensor, tr(F) = 1, and trðHn Þ ¼ trðHt Þ ¼ f 0 . Fig. 13 shows the evolution of ac, an, at, and ac + an + at with axial strain for the two packings. Similar trends are observed for all anisotropies, which initially increase with increasing axial strain and then relax to residual values after reaching a peak. The contact normal and normal force distributions evolve to the largest degree of anisotropy with the occurrence of the peak mobilized shear stress. However, the tangential force anisotropy, at, reaches a maximum value before the loading strain corresponding to the peak shear stress. In both assemblies, the contact orientation anisotropy, ac, and contact normal force anisotropy, an, represent the major contribution to the overall anisotropy, ac + an + at. In contrast, the tangential force anisotropy, at, contributes no more than 10% of the overall anisotropy. The overall anisotropy, ac + an + at (exclusive of the branch vector anisotropy), correlates well with the simulated normalized shear stress, q/p, with respect to both trend and quantity. The anisotropies in the unbreakable packing always remain above those in the breakable packing. This indicates that the unbreakable assembly is inclined to develop larger fabric and force anisotropies than the breakable one, which subsequently results in a larger macroscopic shear strength. With the recent development of the force chain theory, it has become well accepted that the mechanical behavior of granular materials is inherently governed by the continual collapse of old force chains and the formation of new force chains, where force chains are defined as columnar particle structures carrying larger than average loads that are laterally confined by weak network particles [39,40]. The common failure mode of force chains is considered to be bulking, which occurs suddenly when the laterally weak network particles cannot provide enough confinement. Fig. 14 shows the probability density function (PDF) and polar representation of the normal contact forces for two packings of unbreakable and breakable particles. It was observed that in both assemblies, the majority of the particles carry less than the average load, and the number of particles carrying more than the average load decreases exponentially with an increasing contact force. We see that the homogeneity of the normal forces becomes higher when the particles are breakable, which means that the proportion of weak contacts represent the majority of all contact force interactions. The unbreakable assembly provides a more inherently stable structure, which results in an almost instantaneous bulking of the force chains. The drop in deviator stress is most often the result of the instantaneous collapse of force chains by bulking (as shown in the inset of Fig. 10a). In contrast, for the breakable assembly, the particles along the force chains can break, which causes these columnar structures to yield, even though they have not reached their bulking limit. Hence, the force chain collapse is more gradual. This is the most obvious explanation of why the assembly of

G. Ma et al. / Computers and Geotechnics 61 (2014) 132–143

139

Fig. 13. Evolution of anisotropies, with the axial strain for two assemblies: (a) ac, (b) an, (c) at, (d) ac + an + at.

breakable particles does not show a significant drop in deviator stress. The polar representations illustrate the distribution of the normal contact force during progressive loading. At any stage of shearing, the normal contact force of the breakable packing has a smaller magnitude and a capsule-shaped distribution. In addition, the distributions of normal contact forces in an unbreakable packing are peanut shape distributions, which indicate a higher degree of shear induced anisotropy. Following the force chain theory, the fabric anisotropy comes from the fact that, in biaxial tests, the

strong force chains carrying the majority of the load are aligned approximately in the vertical direction, while the magnitude of the normal contact forces in the horizontal direction remain relatively small. One fundamental aspect of the behavior of granular materials is their inter-particle Here, we consider the friction mobiliza friction.  tion index, If ¼ ftc =lfnc , at each contact, where l is the interparticle friction coefficient and ftc and fnc are the actual values of the normal and tangential forces, respectively, at contact c.

Fig. 14. Evolution of the probability density function and polar representations of normal contact forces for two assemblies: (a) 4% axial strain and (b) 8% axial strain.

140

G. Ma et al. / Computers and Geotechnics 61 (2014) 132–143

Fig. 15. Evolutions of the average friction mobilization in the whole contact network, strong subnetwork, and weak subnetwork for two assemblies: (a) unbreakable assembly and (b) breakable assembly.

Fig. 16. Evolution of the accumulated fraction of broken CIEs with the axial strain and spatial distribution of AEs at the end of loading under different confining pressures.

Fig. 18. Evolution of broken CIEs with axial strain under a confining pressure of 2.0 MPa.

Fig. 19. Particle magnification of the breakable assembly at the end of loading.

Fig. 17. Evolution of the damage dissipation energy in the course of shearing under different confining pressures.

Fig. 15 displays the evolution of the average friction mobilization index for the whole contact network, strong subnetwork, and weak subnetwork. For the strong and weak contact, we use the criterion described by Azéma et al. [36], i.e., we consider a particle to belong

to the strong contact network if the normal contact force is larger than the average value from the assembly of all particles; otherwise, it belongs to the weak network. As shown in Fig. 15, the frictional force is not uniformly mobilized at all the contacts. The strong subnetwork is composed of weakly mobilized contacts, while the majority of the strongly mobilized contacts belong to the weak subnetwork. This means that contacts with higher

G. Ma et al. / Computers and Geotechnics 61 (2014) 132–143 Table 2 Strength parameters for three modeled materials. Strength parameters

BT100

BT150

BT200

Compressive strength (MPa) Tensile strength (MPa) Cohesion intercept (MPa) Friction angle of cohesive element

100 10 23.31 40

150 15 34.97 40

200 20 46.63 40

normal forces are not necessarily associated with larger mobilized frictional forces. This was also evidenced in Azéma and Radjaï’s work, which showed that not only the average friction mobilization index but also the number of highly mobilized contacts are larger in the weak force network [41]. In both assemblies, the strong and weak subnetworks approach the friction mobilization to a similar extent, which demonstrates the tendency of granular packing to develop uniformly mobilized contacts in the critical state. As observed, the unbreakable assembly has the ability to mobilize larger frictional forces compared to the breakable assembly. 5.3. Particle breakage behavior This section comprehensively considers the breakage behavior of the modeled material. Fig. 16 shows the evolution of the accumulated fraction of broken CIEs for simulations with different confining pressures. This figure demonstrates that higher confining pressures produce a larger amount of particle breakage. The particle breakage has a negligible effect on the mechanical behav-

141

ior of granular materials at a small stress level. However, in the case of a high stress level, this effect is significant and cannot be ignored. The energy dissipated by the damage and failure of CIEs confirms this conclusion (see Fig. 17). Fig. 16 also displays the spatial distribution of the broken CIEs at the end of loading. Each broken CIE is shown by a dot. In all simulations, the breakage occurs mainly in the region of the shear band, a tendency which becomes more evident when the stress level increases. Fig. 18 shows the evolution of both the number and accumulated fraction of the broken CIEs in the course of loading with a confining pressure of 2.0 MPa. Each red bar is the total number of broken CIEs happened in every 0.002 of axial strain. It is interesting to observe that the particle breakage process is synchronized with the development of the macroscopic stress. The number of broken CIEs rapidly increases at a small strain level and experiences a reduction after the peak value. Fig. 19 displays a partial magnification of the breakable assembly at the end of loading. 5.4. Effect of particle strength on mechanical behavior of rockfill materials Three different biaxial tests are carried out on identical particle assemblies with different particle strength parameters to investigate the effect of particle strength on the mechanical behavior of rockfill materials. The strength parameters are shown in Table 2. The friction angle of the cohesive element and the ratio of uniaxial tensile strength to uniaxial compressive strength are held constant in this comparison. In the following, to more easily reference the simulations, we adopt a code with the format ‘‘BTa’’ for a particular

Fig. 20. Evolution of the deviator stress and volumetric strain with the axial strain for different strengths of particles.

Fig. 21. Evolution of the accumulated fraction of broken CIEs and damage dissipation energy with the axial strain for different strengths of particles.

142

G. Ma et al. / Computers and Geotechnics 61 (2014) 132–143

simulation, in which ‘‘BT’’ represents biaxial test and ‘‘a’’ is the value of the uniaxial compressive strength, e.g., BT100. All of these simulations are performed under the same conditions, with a confining pressure of 2.0 MPa. As illustrated in Figs. 20 and 21, the macroscopic responses of an assembly with higher strength are characterized by strong strain-softening and dilation. The amount of particle breakage is reduced with increasing particle strength, which results in a higher shear resistance and rate of dilation.

the assembly. In general, the results herein show that the proposed method is capable of reproducing the complex mechanical behavior of granular geomaterials and can accurately and conveniently simulate the particle breakage phenomenon.

Acknowledgment This work was supported by the National Natural Science Foundation of China (Grant No. 51379161, 51322905).

6. Conclusions In this paper, a practical combined FE/DE method has been developed to simulate particle breakage in granular geomaterials, e.g., rockfill materials. In this combined approach, the calculations of the particle deformation, stress, and contact force are formulated on the basis of the FE discretization of each particle. To simulate particle breakage, zero-thickness CIEs with a well-established progressive damage model are inserted between the adjacent solid elements. When the overall damage variable of a CIE reaches the upper bound, the degraded CIE is removed from the FE mesh configuration, and a new contact is then established between the previously joined solid elements. There are two prominent benefits of using this method. First, an arbitrary particle shape can easily be introduced without any modification of the contact configuration because both the contact detection and interaction algorithms are formulated on the basis of the FE discretization instead of the particle geometry. Second, a vast range of alternative (e.g., nonlinear constitutive or internally fracturing) properties can be introduced for individual particles. There are additional advantages to this method, e.g., both rigid and flexible boundaries are available. The simulated stress–strain–dilation responses for both breakable and unbreakable assemblies qualitatively agreed with experimental observations, which demonstrated that the proposed method is capable of reproducing the fundamental features of modeled granular geomaterials. The macroscopic response of the breakable assembly was characterized by a contractive and mild strain softening behavior. The unbreakable assembly exhibited higher peak shear strength, followed by strong strain softening. Additionally, the results indicated that particle breakage reduced the macroscopic friction angle but increased the deformability of modeled materials. The anisotropies in the unbreakable packing always remained above those in the breakable packing. This indicated that the unbreakable assembly was inclined to develop larger fabric and force anisotropies than the breakable one, which further resulted in a larger macroscopic shear strength. The probability distribution of the normal contact force revealed that the distribution inhomogeneity was more obvious in the unbreakable assembly. For the breakable assembly, the strong force chains tended to collapse before reaching their buckling limits, which resulted in a polar distribution of the contact force, with a smaller magnitude and less anisotropy. By exploring the interparticle friction characteristics, we found that the unbreakable assembly had the ability to develop a larger mobilized frictional force at contact. The numerical results showed that higher confining pressures resulted in a larger amount of particle breakage. Additionally, the particles in the shear band were more likely to break. An interesting observation is that the particle breakage process was synchronized with the development of the macroscopic stress. This tendency may be attributed to the initiation of shear and tensile failures during the shearing process. The effect of the parent rock strength on the particle breakage and mechanical behavior of granular geomaterials was also discussed. It was found that the amount of breakage increased in weak rocks, which caused a decrease in the peak shear strength and an increase in the compressibility of

References [1] Hardin BO. Crushing of soil particles. J Geotech Eng 1985;111(10):1177–92. [2] McDowell GR, Bolton MD, Robertson D. The fractal crushing of granular materials. J Mech Phys Solids 1996;44(12):2079–101. [3] Lade PV, Yamamuro JA, Bopp PA. Significance of particle crushing in granular materials. J Geotech Eng 1996;122(4):309–16. [4] Varadarajan A, Sharma KG, Venkatachalam K, et al. Testing and modeling two rockfill materials. J Geotech Geoenviron 2003;129(3):206–18. [5] Robertson D, Bolton MD. DEM simulations of crushable grains and soils. Proc. Powders and Grains, Sendai; 2001. p. 623–6. [6] Jensen RP, Plesha ME, Edil TB, Bosscher PJ, Kahla NB. DEM simulation of particle damage in granular media-structure interfaces. Int J Geomech 2001;1(1):21–39. [7] Cheng YP, Nakata Y, Bolton MD. Discrete element simulation of crushable soil. Geotechnique 2003;53(7):633–41. [8] Harireche O, McDowell GR. Discrete element modelling of cyclic loading of crushable aggreates. Granul Matter 2003;5(3):147–51. [9] Bolton MD, Nakata Y, Cheng YP. Micro-and macro-mechanical behaviour of DEM crushable materials. Geotechnique 2008;58(6):471–80. [10] Wang J, Yan H. On the role of particle breakage in the shear failure behavior of granular soils by DEM. Int J Numer Anal Meth Geomech 2011;37:832–54. [11] Alaei E, Mahboubi A. A discrete model for simulating shear strength and deformation behaviour of rockfill material, considering the particle breakage phenomenon. Granul Matter 2012;14(6):707–17. [12] Tsoungui O, Vallet D, Charmet JC. Numerical model of crushing of grains inside two-dimensional granular materials. Powder Technol 1999;105(1):190–8. [13] Lobo-Guerrero S, Vallejo LE. Discrete element method evaluation of granular crushing under direct shear test conditions. J Geotech Geoenviron Eng 2005;131(10):1295–300. [14] Brosh T, Kalman H, Levy A. Fragments spawning and interaction models for DEM breakage simulation. Granul Matter 2011;13(6):765–76. [15] Bruchmüller J, van Wachem BGM, Gu S, Luo KH. Modelling discrete fragmentation of brittle particles. Powder Technol 2011;208(3):731–9. [16] Potapov AV, Campbell CS. Computer simulation of shear-induced particle attrition. Powder Technol 1997;94(2):109–22. [17] Seyedi Hosseininia E, Mirghasemi AA. Numerical simulation of breakage of two-dimensional polygon-shaped particles using discrete element method. Powder Technol 2006;166(2):100–12. [18] Galindo-Torres SA, Pedroso DM, Williams DJ, Li L. Breaking processes in threedimensional bonded granular materials with general shapes. Comput Phys Commun 2012;183(2):266–77. [19] Munjiza A, Owen DRJ, Bicanic N. A combined finite–discrete element method in transient dynamics of fracturing solids. Eng Comput 1995;12(2):145–74. [20] Munjiza A. The combined finite–discrete element method. New York: Wiley; 2004. [21] Ransing RS, Gethin DT, Khoei AR, Mosbah P, Lewis RW. Powder compaction modelling via the discrete and finite element method. Mater Des 2000;21(4):263–9. [22] Gethin DT, Ransing RS, Lewis RW, Dutko M, Crook AJL. Numerical comparison of a deformable discrete element model and an equivalent continuum analysis for the compaction of ductile porous material. Comput Struct 2001;79(13):1287–94. [23] Lewis RW, Gethin DT, Yang XS, Rowe RC. A combined finite–discrete element method for simulating pharmaceutical powder tableting. Int J Numer Meth Eng 2005;62(7):853–69. [24] Gethin DT, Yang XS, Lewis RW. A two dimensional combined discrete and finite element scheme for simulating the flow and compaction of systems comprising irregular particulates. Comput Meth Appl Mech Eng 2006;195(41):5552–65. [25] Frenning G. An efficient finite/discrete element procedure for simulating compression of 3D particle assemblies. Comput Meth Appl Mech Eng 2008;197(49):4266–72. [26] Procopio AT, Zavaliangos A. Simulation of multi-axial compaction of granular media from loose to high relative densities. J Mech Phys Solids 2005;53(7):1523–51. [27] Harthong B, Jérier JF, Dorémus P, Imbault D, Donzé FV. Modeling of highdensity compaction of granular materials by the discrete element method. Int J Solids Struct 2009;46(18):3357–64. [28] Systèmes D. Abaqus 6.10: Analysis user’s manual. Providence, RI: Dassault Systèmes Simulia Corp; 2010.

G. Ma et al. / Computers and Geotechnics 61 (2014) 132–143 [29] Mollon G, Zhao J. Fourier–Voronoi-based generation of realistic samples for discrete modelling of granular materials. Granul Matter 2012;14(5):621–38. [30] Yang ZJ, Su XT, Chen JF, Liu GH. Monte Carlo simulation of complex cohesive fracture in random heterogeneous quasi-brittle materials. Int J Solids Struct 2009;46(17):3222–34. [31] Day RA, Potts DM. Zero thickness interface elements—numerical stability and application. Int J Numer Anal Methods Geomech 1994;18(10):689–708. [32] Benzeggagh ML, Kenane M. Measurement of mixed-mode delamination fracture toughness of unidirectional glass/epoxy composites with mixedmode bending apparatus. Compos Sci Technol 1996;56(4):439–49. [33] Camanho PP, Davila CG, De Moura MF. Numerical simulation of mixed-mode progressive delamination in composite materials. J Compos Mater 2003;37(16):1415–38. [34] Li J. Three-dimensional effect in the prediction of flange delamination in composite skin-stringer pull-off specimens. J Compos Tech Res 2002; 24(3):180–7.

143

[35] Kruyt NP. Contact forces in anisotropic frictional granular materials. Int J Solids Struct 2003;40(13):3537–56. [36] Azéma E, Radjai F, Peyroux R, Saussine G. Force transmission in a packing of pentagonal particles. Phys Rev E 2007;76(1):011301. [37] Ouadfel H, Rothenburg L. Stress–force–fabric’ relationship for assemblies of ellipsoids. Mech Mater 2001;33(4):201–21. [38] Azéma E, Radjaï F. Stress–strain behavior and geometrical properties of packings of elongated particles. Phys Rev E 2010;81(5):051304. [39] Muthuswamy M, Tordesillas A. How do interparticle contact friction, packing density and degree of polydispersity affect force propagation in particulate assemblies? J Stat Mech-Theory E 2006;2006(09):P09003. [40] Tordesillas A, Zhang J, Behringer R. Buckling force chains in dense granular assemblies: physical and numerical experiments. Geomech Geoeng: Int J 2009;4(1):3–16. [41] Azéma E, Radjaï F. Force chains and contact network topology in sheared packings of elongated particles. Phys Rev E 2012;85(3):031303.