Effect of Particle Shape on Mechanical Behaviors of Rocks: A Numerical Study Using Clumped Particle Model

Since rocks are aggregates of mineral particles, the effect of mineral microstructure on macroscopic mechanical behaviors of rocks is inneglectable. Rock samples of four different particle shapes are established in this study based on clumped particle model, and a sphericity index is used to quantify particle shape. Model parameters for simulation in PFC are obtained by triaxial compression test of quartz sandstone, and simulation of triaxial compression test is then conducted on four rock samples with different particle shapes. It is seen from the results that stress thresholds of rock samples such as crack initiation stress, crack damage stress, and peak stress decrease with the increasing of the sphericity index. The increase of sphericity leads to a drop of elastic modulus and a rise in Poisson ratio, while the decreasing sphericity usually results in the increase of cohesion and internal friction angle. Based on volume change of rock samples during simulation of triaxial compression test, variation of dilation angle with plastic strain is also studied.


Introduction
Rocks are made from minerals, which are various in chemical composition and crystal morphology. Most Rocks in nature are composed of irregular mineral particles strongly bonded together [1]. The shapes of mineral particles are crucial for mechanical behaviors of rocks. Particle shape has been recognized to affect the mechanics characters of granular material, which has been revealed in several publications [2,3]. The internal friction angle of different shaped particles was investigated through triaxial compression test by Shinohara et al. [4]. Dodds [5] expounded the influences of particle shape and stiffness effects on soil behavior. Liu et al. [6] quantified particle shapes by digital image method and expounded the correlation of mechanical performance and shape factors defined in their paper. Since various particle shapes in a material, it is very difficult to examine the influence of a specific particle shape. Johanson [7] used plastic pellets of different shapes (round, heart, and stars) coated with soft Tacky Wax to make samples, respectively. Although this method can obtain consistent samples of distinct shape, it is not proper to produce rock or sand samples because of different performance between plastic pellets and these materials.
Many scholars studied mechanical behaviors of granular matter using numerical simulation method [8][9][10]. Compared with other continuum approaches, the particle mechanics approach based on discrete element method can reproduce the processes of fracture initiation and growth at quasimicro-or macroscopic levels and treat the problems relating to discontinuous large deformation with some simple assumptions and constitutive models [11,12]. Particle flow codes (PFC2D and PFC3D) are the most widely used particle mechanics codes [1,11]. Numerical tests by Kock showed that composition and texture of sediments were relevant to frictional strength and development of shear zone [13]. Furthermore, with the help of PFC2D, peak strength, internal friction angle, and thickness of shear zones were also proven to be connected with particle shape [14][15][16].
However, researches referring to this issue are still not wide and enough; a systematic study to this problem is 2 The Scientific World Journal of urgent demand for understanding mechanism of deformation and failure from quasi-microlevels. And past study results concentrated on materials such as sand or soil, but a few papers discussed how particle shapes in rock affect the mechanical behaviors of rocks [17]. Unlike loose particle materials sand and soil, rock particles are strongly bonded together. The cements in rocks cause different interaction mechanisms between rock particles. Consequently, Potyondy and Cundall [1] proposed a bonded-particle model for rock, where rock is represented by a dense packing of spherical particles that are bonded together. For different particle shapes, there is no a unified quantitative evaluation method due to their complexity. Some results described particle shape mainly based on the qualitative approach, lacking quantitative analysis of particle shapes.
Hence, the main purpose of this paper is to examine the influences of particle shape on the mechanical behaviors of rocks. From the result of quartz sandstone triaxial compression test and mineral particle shape in quartz sandstone, four representative particle shapes were created. We use sphericity index as the particle shape factor to characterize four representative particles. Then mechanical behaviors of four samples formed by four representative particles were studied, respectively.

Basic Theory of Particle Flow Method
A general particle flow model simulates the mechanical behavior of materials based on discrete element method. The Newton's laws of motion and the force-displacement law provide the fundamental relationship among force, displacement and particle motion.
Mechanical behaviors of materials are simulated in terms of the movement of each particle and the interparticle forces acting at each contact point. At each contact point, contact behaviors consist of stiffness, slip, and bond [18].
In the contact normal direction, the stiffness behavior provides the relation among the contact normal force component , the total normal displacements , and the contact normal stiffness (unit: Pa/m) . In the contact tangential direction, the stiffness behavior relates the increment of shear force Δ , and the increment of shear displacement Δ , as follows: where and are the contact normal stiffness and tangent stiffness (unit: Pa/m), respectively. In this paper, the linear contact models are adopted, where and are independent of displacement. The slip behavior relates slip condition, given by where max is maximum allowable shear contact force and is the friction coefficient. At every contact point, if | | > max , then slip is allowed to occur. And during the next calculation cycle, the magnitude of is set to max . The bond behavior acts as a kind of glue joining the two particles. PFC3D can simulate not only the loose particle materials such as sand and soil but also the rock-like materials where particles are strongly bonded together [1]. A contact-bond model and a parallel-bonded model are two widely used bond behaviors. Compared with the contactbonded model, the parallel-bonded model is more suitable for rock, as the parallel-bonded model can transmit both forces and moments at contact point between particles, and this better meets the interaction mechanism of particles. And the stiffness in the parallel-bonded model consists of the contact stiffness and bonding stiffness. In the parallelbonded model, the bond breakage will cause bonding stiffness inactive, which results in stiffness reduction. In this respect, the performance of stiffness reduction is consistent with experiments where the rock-like materials may failure in either tension or shearing with an associated reduction in stiffness [12,19]. So the parallel-bonded model is used in this paper.
In a 3D model, the bonding areas of two particles with parallel bond are equivalent to a disk (Figure 1(a)). , , , and are tensile force, shear force, bending moment, and twisting moment acted on the bonding cross-section, respectively. The maximum normal stresses max and shear stresses max carried by the bonding material can be written as [18] where and are the area and moment of inertia, respectively, of the parallel-bond cross-section. And is the polar moment of inertia of the parallel-bond cross section. is the radius of the disk (Figure 1(b)). The parallel-bond breaks when either normal stresses or shear stresses exceed corresponding maximum stresses.
The Scientific World Journal 3

Generation of Representative
Particles and Samples. The mineral particles formed by different morphological structures of crystals in rocks are irregular and often have no fixed shapes. In a numerical simulation of granular material, some representative particles substitute for the real particles [20]. According to the characteristic of quartz sandstone particles, four kinds of representative particles were built to simulate complex particle shapes in quartz sandstone, as shown in Figure 2(b). A common feature of these representative particles adopted in this paper is that they are all made up of a big spherical particle and several small spherical particles. The big spherical particle acts as the main body of representative particles, and the several small spherical particles act as the rugged edges on the surface of rock particles. Figure 2(b) shows the relationship of size between the big spherical particle and the small ones. The samples in Figure 2(c) were formed by four representative particles respectively, and they were used to research mechanical behaviors of different particle shapes. The sample 5 ( Figure 2(a)) was created by combinations of the representative particles from 1#-4# to simulate real quartz sandstone. The representative particle 1# is a ball, which can be created directly in PFC3D. The representative particles 2#-4# can be generated using clumped particles. The clumped particles model is widely used model to generate nonspherical particles [15,16]. The clumped particles are as a whole consisting of two or more balls in PFC3D. The clumped particles can create a group of cement particles that behave as a single rigid body, and these clumped particles may overlap as a deformable body that will not break apart regardless of the forces acting upon it [18].
Sample 1 was first created; then the balls in sample 1 were replaced with representative particles 2#-4#, which produced samples 2-4, respectively. The replacements obey the following three principles [18]. (1) Replace ball with clumped particle, each of which has the same volume as the ball that it replaces. The porosity of sample remains the same after replacement. (2) Each clumped particle is oriented randomly by rotating them about the axes by a random angle.
(3) Volume-based centroid of the clumped particle coincides with centroid of the replaced ball.
Considering the generation approach of samples 2-4, the number of balls in sample 4 is four times greater than in sample 1. Although balls in a clumped particle is treated as a single rigid body in calculation cycle, the replacement process mentioned previously will cost too much time if too many balls are generated in the sample 1. To improve calculation speed, we suggest not using too many balls in sample 1.

Quantitative Description of Particle
Shapes. Shape factors to quantify particle shapes are discussed in detail within a number of literatures [6,21]. Roundness and aspect ratios are very commonly used shape factors [3,9]. To two-dimensional particles, Kong and Peng [14] suggested to determine shape  factors by (4) after analyzing the mechanical behaviors from mesoscopic levels: where is shape parameter, 1 and 2 are two parameters related to roundness of particle and roughness of surface, respectively, with and being corresponding weighting coefficients. The Fourier analysis technique also an important method to characterize particle shapes. The digitized particle outline is described by Fourier analysis which is originally based on the -, -coordinate detection method for contour curves [4,22].
In this paper, we discuss the mechanical response of three-dimensional particles. Combined with previous shape factors method, considering the feasibility and utility of method, we define sphericity as a shape factors for particles where s is surface area of a sphere whose volume is the same with the particle. is the surface area of particle. is 4 The Scientific World Journal  sphericity of particle. Equation (5) can be expressed by the volume and surface area of particle, where is the particle volume In order to obtain the volume and surface area of representative particles, four representative particles are created in AutoCAD, and their volumes and surface areas are easy to get by AutoCAD. Then the sphericities of four representative particles are counted by (4). The result is shown in Table 1. As is shown in Table 1, actually, sphericity expresses the degree of similarity between the sphere and particle. The closer the shape of a particle comes to a sphere, the nearer the sphericity approximates to 1.
All the samples are cylinders with height of 80 mm and diameter of 40 mm. The minimal radius of the balls is 2.5 mm, and the ratio of the maximum radius to the minimal radius is 1.5. Figure 3 shows the location of load planes and a sample. The upper plane and lower plane are two load planes,  Figure 5). In laboratory experiment, confining pressure is 8 MPa. The quartz sandstone specimens are loaded until destruction, a rupture plane can be seen from Figure 4. Meanwhile, the numerical sample 5 ( Figure 6) is used to simulate triaxial compression test in PFC. Figure 6(a) is sample 5 before the numerical test, and Figure 6(b) shows the distribution of microcracks at the peak stress. The microcracks are showed in white. When a bond breaks either for tension or shearing failure, a new microcrack forms in PFC. From Figure 6(b), we can see the microcracks are concentrated in the diagonal, top, and bottom of the sample. In theory, the microcracks should mainly distribute in the maximum shear stress plane. But in this numerical test, the loads from load planes transmit to the sample only by the contact point of particles and load planes which results in stress concentration easily at the interfaces of sample and load planes. And the bond strength is not homogeneous the standard deviations of bond strengths is showed in Table 2.
In the process of numerical calibration, when the microscopic parameters of sample 5 correspond with the parameters of quartz sandstone, the stress-strain curves both experiment and numerical simulation should be similar ( Figure 6). Figure 7 shows the comparison of experiment and numerical simulation. As can be seen in Figure 7, two curves come close before stress get to peak strength. There are slight differences in postpeak region. The main reason is small stiffness of test machine and slow servo response speed of test machine, which cause the large intervals and uneven The Scientific World Journal  characteristics of measuring point. So it is reasonable to exist slight differences. The agreement between numerical simulation and laboratory results is acceptable (Figure 7), which shows microscopic parameters adopted in the numerical simulation are suitable for the quartz sandstone. The corresponding model parameters are listed in Table 2.

The Influences of Particle Shapes on Strength.
To study the relation of particle shapes and mechanical behaviors, numerical triaxial compression test of samples 1-4 ( Figure 2) was made under different confining pressures with model parameters listed in Table 2.
Stress-strain curves with confining pressures 2-20 MPa are shown in Figures 8(a)-8(d). The sphericity of particles is noted in parentheses. As can be seen in Figure 8, particle shapes affect stress-strain curve considerably. This impact manifests mainly in the relation between peak strength and particle shapes. Specifically, peak strength decreases with the increasing of sphericity of particles under the same confining pressure. The modulus of elasticity of samples varies with particle shapes, which can be seen from the slope of curves before peak point. The sample with smaller particles sphericity has a higher modulus of elasticity.
Comparing stress-strain curves of different confining pressure (Figures 8(a)-8(d)), we found the influences of particle shapes on residual strength depend on confining pressure. With a high confining pressure and a small sphericity of particles, stress drop at peak point is inconspicuous. Under such condition, the residual strength is close to peak strength. As can be seen from Figure 8 At mesoscale level, the results can be explained as follow. Before stress go to peak point, particles bond and interlock together. The smaller the particles sphericity, the greater the degree of interlocking is, correspondingly, the greater the overall peak strength. In the post-peak stage, bonds of particles crack. Particles with big sphericity approaching the spheroid, slip and roll easily, which leads to rapid stress drop. While for small sphericity particles, the effect of interlocking, and friction is stronger, which makes samples remain a higher residual strength in the postpeak stage.

The Influences of Particle Shapes on Crack Initiation
Stress and Crack Damage Stress. In the failure process of rock, crack initiation stress and crack damage stress are two important indicators. The crack initiation stress ci marks crack initiation and stable propagation. When stress exceeds crack damage stress cd , the crack growth is unstable, and cd is also the beginning of rock dilation.  Figure 9: Schematic diagram of stress-strain curves of rocks.
According to the research results of progressive failure process by Martin [23], Eberhardt [24], Eberhardt et al. [25], and Diederichs et al. [26], there are three major ways to ascertain ci , as is shown in Figure 9. (1) Acoustic Emission Test (AE). Eberhardt put out that ci is the stresses when new AE counts first rise above background. (2) The stressvolumetric strain curves. When stresses get to ci , the stressvolumetric strain curves deviate from the elastic line. (3) Crack volumetric strains. The crack volumetric strains deviate from zero with ci coming. The crack volumetric strains are defined as (7) by Martin [23]: where cv and V are the crack volumetric strains and volumetric strains, respectively, ] is poisson ratio, and is elastic modulus. The crack damage stresses cd can be ascertained by AE or the stress-volumetric strain curves. When stress reaches cd , AE counts curves show transition and AE counts increase rapidly. In addition, from the stress-volumetric strain curves we can observe volumetric strains rate is near zero and a transition point arises with the crack damage stresses, cd coming, which are shown in Figure 9.
This paper determined the crack initiation stress and crack damage stress of rock samples by crack volumetric strains and the stress-volumetric strain curves. Relevant curves of four samples under confining pressure 15 MPa are shown in Figure 10. As mentioned before, the stress that crack volumetric strains deviate from zero is ci , marked with green circle in Figure 10. And the cd can be obtained from volumetric strain reversal, marked with red circle in Figure 8.
Further, Figure 11 shows how the crack initiation stress and crack damage stress vary with particles sphericity. As we can see from Figure 11, the crack initiation stress and crack damage stress reduce with particles sphericity increasing. In other words, particle shapes affect crack initiation and unstable crack growth in the process of rock failure. With a smaller particles sphericity, rocks can support a higher load before appearing crack initiation and unstable crack growth.

The Influences of Particle Shapes on Cohesion and Internal Friction Angle. Cohesion and internal friction angle
are two important material parameters of the Mohr-Coulomb strength theory. The Mohr-Coulomb yield criterion is expressed by (8) using principal stresses. The Mohr-Coulomb yield criterion assumes the cohesion and internal friction angle of rock are both constant, but this assumption is limitation in practical application. Since some scholars [27][28][29] found cohesion and internal friction angle of rock were not constant in triaxial compression test, where 1 and 3 are maximum principal stress and minimum principal stress, respectively. Supposing rock samples yield at peak strength, and strength and confining pressure are substituted into (8), we get (9): Equation (9) shows the relation between the peak strength and confining pressure under triaxial condition.
Based on numerical triaxial compression tests with the confining pressure of 2 MPa, 8 MPa, 15 MP, and 20 MPa, − curve of sample 1 is displayed ( Figure 12). As is shown by Figure 12, there is a good linear relation between and . Then − curve can be obtained by beeline fitting, which is given in Figure 12. Comparing the linear equation in Figure 12 and (9), cohesion and internal friction angle can be calculated by the gradient and the constant of linear equation (see (10)): = arcsin ( − 1 + 1 ) , Along the same ways, cohesion and internal friction angle of other samples are acquired. Figure 13 shows how cohesion and internal friction angle vary with particles sphericity. As can be found in Figure 13, basically, cohesion and internal friction angle of samples reduce with particles sphericity rising. When particles sphericity is between 0.92 and 0.96, the variation of internal friction angle is gentle.  -direction. In numerical triaxial compression tests, elastic modulus and Poisson ratio ] can be calculated as follows:

The Influences of Particle Shapes on Elastic Modulus and
where Δ , Δ , and Δ are strain increments in -,and -directions, respectively. And Δ is volumetric strain increment, Δ is stress increment in direction. Ii is obviously that and ] change with load process. But, for simplicity sake, elastic modulus and Poisson ratio ] are subject to the point of half peak stress. Figure 14 shows how elastic modulus and Poisson ratio vary with particles sphericity under confining pressure 15 MPa. As particles sphericity increase, elastic modulus of samples drops and Poisson ratio rises. As can be seen from the magnitude of Poisson ratio, there are only small differences between Poisson ratio in spite of sphericity influence.

Dilation Effect of Samples.
In the process of numerical triaxial compression test, the samples are compressed first, and then they are dilated. In continuum mechanics, dilation angle is the widely used parameter to discribe the dilation effect. Dilation angle is not a constant during the deformation process of rock [30]. But, considering the rock materials may not obey Drucker's stability postulate [31], which recites that the work of the external agency on the displacement produced must where is a plastic multiplier and is plastic potential. One of the commonly used plastic potential assumptions states that [32] = 1 − ( , ) 3 , where is dilation angle and 1 , 3 are the maximum and minimum principal stress, respectively. is the stress tensor and is the plastic parameter. Here, can be expressed as shear plastic strain, as follows: The dilation angle of rock can be obtained by (15) [33], wherėV anḋ1 are volumetric and axial plastic strain increments, respectively, Volumetric and axial plastic strain can be obtained by loading and unloading cycles [32,34]. Since the dilation effect in prepeak stage is not obvious, we focus on evolvement laws of dilation angle in postpeak stage. Figure 15 shows the stress-strain curve from loading and unloading cycles under confining pressure 15 MPa. The numerical test process is as follows.
The confining pressure and axial stress of 15 MPa are applied on samples first, in which case samples work in elastic behavior. That is plastic strains still do not occur; this state is called original state. Then axial stresses increase to load samples and axial strain is recorded from the original state. When axial strains reach to 0.008, 0.009, 0.0100, 0.011, 0.012, 0.013, 0.014, 0.015, 0.016, and 0.017, unloading, respectively until samples return to the original state. The strains in original state are axial plastic strains 1 , marked with red circle in Figure 15, and then load to the next unloading point, and so forth. In the process of loading and unloading cycle, volumetric strains are recorded, as it is shown in Figure 16.
The volumetric strains V are marked with red circle in Figure 16. Equation (15) is conveniently written as incremental form, as follows: The Scientific World Journal 11 (a) (b) (c) Figure 18: Simulating complex particle shape.
Given the deformation conditions of triaxial compression test, the following formula is deduced.
Dilation angles can be calculated according to (16). The volumetric strains V and axial plastic strain 1 can be obtained from loading and unloading cycles in Figures 15 and  16. Shear plastic strain can be calculated using (14) and (17). Figure 17 shows the evolvement process of dilation angle with shear plastic strain. As can be seen from Figure 17, dilation angles increase rapidly in the beginning. Dilation angles reach a maximum around shear plastic strain 0.003. Then, there was no significant change on dilation angles. According to the research of Alejano and Alonso [32] and Zhao and Cai [34], dilation angles will start to fall when shear plastic strain is larger. As shown in Figure 17, although dilation angles of four samples have a little difference, there are no essential distinctions on trends of curves among four samples. That is, particle shapes only impact size of dilation angles rather than trends of curves.

Discussions
The real shape of mineral particles in rocks is quite complex, while the shape of representative particles used in this paper is simple. Theoretically, the method forming representative particles by clumped particles is also effective for complex particles ( Figure 18). The red lines represent the planar contour for a complex particle, and the black lines represent the clumped particles in Figure 18. As is shown in Figure 18(c), the clumped particles can be a good approximation of the complex particle. That is to say, the real shapes of mineral particles in rocks can be simulated using clumped particles in PFC. Deserved to be mentioned, when the overlap parts in a clumped particle involve over two particles (the green part in Figure 18(c)), the clump volume automatically calculated in PFC is inaccurate [18]. In this case, tools such as AutoCAD can be use to compute the clump volume.
The sphericity is not the only shape factor to characterize particle shapes. As a matter of fact, the number and type of representative particles used in this paper are limited. To be more precise, for this kind of representative particles, which are formed by a big spherical particle and several small spherical particles, the sphericity is a relatively efficient shape factor. However, for other kinds of particles, such as elongated particles, aspect ratios and other shape factors may be more effective. It is verified in our study that there are good correlations between the main mechanical parameters of samples and the sphericity of particles. But, the sphericity does not significantly influence the residual strength and dilation effect. And, when particles shapes get more complex, the relation between the sphericity and the internal friction angle becomes unclear. The main reason for these problems may lie in the fact that the sphericity only expresses the degree of similarity between a sphere and particle. In other words, although the sphericity expresses much of morphological characters especially for simple particle shapes, some morphological characters of particle cannot be entirely reflected by the sphericity. So we need to take into account more parameters to further express particle shapes.

Conclusions
Our numerical experiments show that the mechanical behaviors of rock are influenced by their particle shapes; the primary conclusions are as follows.
(1) The sphericity index is an applicable shape factor to measure particles shapes. The sphericity describes the proximity of a particle to a sphere, which directly influences the interlocking ability of particles.
(2) The crack initiation stress, crack damage stress, and peak stress of rock are affected by particle shapes. Specifically, three stress indices decrease with increasing of sphericity. The increasing sphericity also leads to smaller elastic modulus and larger poisson ratio. And the decreasing sphericity causes particle interlocking of different degrees which restrains slip and rotation, consequently, cohesion and internal friction angle rise.
(3) To samples of different particle shapes, the evolvement process of dilation angle with shear plastic strain was researched. The result shows the trends of dilation angle changing with shear plastic strain are similar, but values of dilation angle have a little difference.