Investigating the Failure Mechanism of Jointed Rock Slopes Based on Discrete Element Method

This paper presents a comprehensive engineering method to investigate the failure mechanism of the jointed rock slopes. The ﬁeld geology survey is ﬁrst carried out to obtain the slope joint data. A joint network model considering the structural complexity of rock mass is generated in the PFC software. The synthetic rock mass (SRM) approach for simulating the mechanical behavior of jointed rock mass is employed, in which the ﬂat-jointed bonded-particle model (FJM) and smooth joint contact model (SJM) represent intact rock and joints, respectively. Subsequently, the eﬀect of microparameters on macromechanical properties of rock is investigated for parameter calibration. Moreover, the scale eﬀect is analyzed by multiscale numerical tests, and the representative elementary volume (REV) size in the selected research area is found as 16m × 16 m × 16m. The microparameters of the SRM model are calibrated to match the mechanical properties of the engineering rock mass. Finally, an engineering case from Shuichang open-pit mine is analyzed and the failure process of the slope during the excavation process from micro- to macroscale is obtained. It has been found that failure occurs at the bottom of the slope and gradually develops upwards. The overall failure of the slope is dominated by the shallow local tension fracture and wedge failure.


Introduction
e rock mass is a geological body with fractures complexly distributed. ese natural discontinuities such as joints and faults are ubiquitous in rock, which have an important impact on the mechanical behaviors of rock material and lead to the anisotropy, heterogeneity, and scale effect [1]. e mechanical properties and failure mechanism of the rock mass are always a major field of rock mechanics. e accurate description of geometrical characteristics of joint network is critical for the safe design, construction, and stability analysis in rock mass engineering. Moreover, investigating the mechanical properties and failure mechanism of jointed rock from micromechanical viewpoint can both enrich rock mass theory and guide the rock engineering practice.
Considerable research studies have been carried out on the joint network modelling and analysis methods of jointed rock mass [2][3][4][5][6][7][8][9][10][11][12][13][14][15][16][17][18][19]. Due to the complexity and uncertainty of joints system in the rock, the stochastic approach based on the theories of probability and statistics is widely used to describe joint geometries from limited measurements such as borehole logging and scanline [2], which is also the mathematical theoretical basis of numerical simulation of joint network. Beacher et al. [3,4] assumed joints to be planar discs and considered the geometrical properties (e.g., spacing, size, and aperture) as independent random variables obeying certain probability distributions. Hudson and Priest [5], Priest and Hudson [6,7], and Rives et al. [8] arrived at the conclusion that joint spacing usually followed a negative exponential, lognormal, or normal distribution depending on the degree of joint saturation in the network. Similarly, the properties of frequency [9], size [10,11], and aperture [12,13] were analyzed based on geological survey. Einstein and Baecher [14] processed measurement data using a rosette or stereogram to group the joints into different sets with their orientation feature by a uniform, normal, or Fisher distribution. Hammah et al. [15][16][17][18][19] presented separately different data analysis methods for optimized and determined dominant orientations, which promoted the progress of statistical method of joint. ree-dimensional modelling and analysis numerical methods have been developed with the computer-aided technology in the past two decades [20][21][22]. Xi et al. [23] investigated the initiation and propagation of precracks in rock. Einstein and Locsin [24] and Ivanova et al. [25] proposed 2D or 3D fracture simulation models, which promoted the application of modelling theory.
Particle-based models, which were originally introduced by Cundall and Strack [26] to simulate granular materials (e.g., soils and sands) and which gradually evolved to commercial code, i.e., particle flow code (PFC) [27], are now widely used for jointed rock mass. Conventional PFC modelling of intact rock considers the bonded-particle model (BPM) [28], which has two limitations, i.e., an unreasonable UCS/TS (compressive to tensile strength) ratio and low friction angle. Some options were suggested to solve these problems: changing the particle size distribution to reduce the porosity [29,30] or changing particle shape into cluster [28] and clumps [31]. Potyondy et al. [32] created a new model, flat-jointed bonded-particle model (FJM), to overcome these limitations and used it to reproduce a realistic hard rock. To suppress such an artificial roughness, smooth joint contact model (SJM), which was independent of the arrangement of local particles in contact, was proposed to simulate the geometric-mechanical characteristics of joints [33]. Based on the two well-established techniques, BPM and SJM, Ivars et al. [33] and Pierce et al. [34] developed a synthetic rock mass (SRM) to represent the mechanical behaviors of jointed rock mass and obtained the scale effect, anisotropy, and brittleness that cannot be obtained using empirical methods of property estimation. e uses of SRM are mostly related to simulating crack propagation and the slip of rock mass under selected loading conditions, such as direct shear test [35][36][37][38][39] and uniaxial and triaxial compression test [40][41][42]. Some researchers employed SRM to study the effect of sample size on rock mass strength and derive the equivalent rock mass properties [43][44][45]. SRM was also used to simulate the mechanical properties and failure mechanism of rock mass in engineering projects [46][47][48][49]. However, due to complicated geometry of joints, SRM has rarely been used for the slope stability analysis of open-pit mine and the failure mechanisms of jointed rock slope during the ore excavation have not been fully understood. e remainder of this paper is organized as follows. Section 2 reviews the main aspects and methodology of the components of SRM. In section 3, based on the statistic characteristics of joints, a numerical model of the jointed rock mass is built by PFC software. In section 4, the microparameters of FJM and SJM are calibrated to reproduce the macroparameters derived from laboratory tests. en, multiscale SRM models are established to investigate the scale effect and determine the representative elementary volume (REV) size by numerical tests. Next, the microparameters are recalibrated to represent the behaviors of rock mass according to the REV size and Hoek-Brown method. In Section 5, a jointed slope of Shuichang open-pit iron mine is selected as an engineering case, and the failure mechanism and evolution process of the slope are analyzed from a micromechanical viewpoint.

Synthetic Rock Mass (SRM) in PFC.
e synthetic rock mass (SRM) approach has been implemented in PFC software, where particles are rigid spherical bodies joined by deformable contacts, to solve problems by using the explicit formulation of the distinct element method (DEM).
is new technique can be used as a virtual laboratory to conduct numerical experiments to represent in a qualitative and quantitative mechanical behavior of jointed rock mass and has already made it possible to extend the volumes of rock at the scale of 10-100 m containing thousands of nonpersistent joints. In PFC software, the SRM brings together two wellestablished techniques: bonded-particle model for intact rock deformation and fracture and discrete fracture network (DFN) for joints. at is to say, the SRM model represents the jointed rock mass as an assembly of intact rock and an embedded DFN. e main components of SRM sample are shown in Figure 1. e following section presents the main aspects of each model.

2.2.
Representation of the Intact Rock. FJM (flat-jointed bonded-particle model) was proposed based on the distinct element method to simulate the mechanical behavior of intact rock in PFC software by Potyondy [32]. is model can overcome the limitations of the standard BPM model and well explain the mechanical behavior of various aspects of rock, especially for hard rock. In FJM, the intact rock is represented as assemblage of rigid circular (in 2D) or spherical (in 3D) particles and piece of finite-size interface between two notional surfaces. A typical flat joint contact ( Figure 2) can be deformed, partially damaged, and broken by the influence of external force. e particle can be considered as a skirted particle, which can provide interlocking and rotational resistance even after the interface breaking. erefore, a more reasonable unconfined compressive to tensile strength (UCS/TS) ratio for hard rock can be obtained through FJM.

Representation of the DFN.
e smooth joint contact model (SJM) simulates the behavior of a smooth interface, regardless of the local particle contact orientations along the interface [33]. is kind of contact is described as smooth because particle pairs joined by a smooth joint contact may overlap and "slide" past each other, instead of being forced to move around one another ( Figure 3).
A smooth joint can be envisioned as a set of elastic springs uniformly distributed over a circular cross section, centered at the contact point and oriented parallel with the 2 Advances in Civil Engineering joint plane. e area of the smooth-joint cross section is given by where R (1) and R (2) are the particle radii and λ is the radius multiplier.
When the unbonded SJM model is specified, the incremental force calculation is performed depending on the elastic portion of the displacement increment (ΔU e n and ΔU e s ). F n � F n 0 + k n AΔU e n , where k n and k s are the normal and shear stiffness, respectively, and (F n ) 0 and (F s ) 0 are the smooth joint normal force and shear force at the beginning of the timestep, respectively. If |F * s | ≤ μF n , then F s � F * s . Otherwise, F s � μF n , sliding is assumed to occur. While slipping, shear displacements produce an increase in normal force due to dilation: where μ is the friction coefficient and ψ is the dilation angle.

Geometrical Characteristics of Joints.
e natural joints often comprise complex network and dominate the mechanical behaviors of jointed rock mass. erefore, it is important to obtain the geometrical characteristic representation of complex two or three-dimensional joint network for modelling jointed rock mass. However, it is impossible to directly obtain the in situ joint data deep in the rock. e details of 3D structures of joint network are usually inferred from the observation information of lower dimensional limited exposures like borehole logging, outcrop scanline, or window mapping. A large number of field survey data show that the joints in the rock mass mostly appear in sets, and it is necessary to identify the dominant sets of joints which highly affect the mechanical behavior of rock.
In this paper, a field survey was carried out in the north of Shuichang open-pit iron mine by using borehole logging and scanline methods. A total of 311 sets of joint are measured, and the data show that joints in this area are well developed. In addition, some joints have potential adverse effects on slope stability because the strike and dip direction of joints are almost the same with those of the slope. Statistical analyses are conducted on the data by using rose diagram or stereogram so that joints can be grouped into different sets with their orientations (Figure 4).
By using the hierarchical cluster method, the measured data were grouped into four dominant orientations, namely, A, B, C, and D (Table 1). And the number of joints in these four orientations is 248, accounting for 79.7% of the total. e characteristic of distribution of joints can be described by geometric parameters such as dip direction, dip angle, and spacing. To provide necessary information for the Advances in Civil Engineering subsequent establishment of joint network model, probability distribution modes of those parameters are studied with orientation characterized by a uniform, normal, or negative exponential distribution, and the histogram and fitting curve are shown in Figure 5. e basic parameters and distributions for dominant orientations in Shuichang open-pit iron mine are obtained as shown in Table 2.

Generation and Testing of Joint Network Model.
e shape of the joint in three-dimensional model is regarded as planar discs of radius in PFC. In this case, the probability density function of the disc's radius is calculated as where λ is the mean value of negative exponential function. e mean radius of joint disc is given as where r is the mean radius of joint disc and l is the mean trace length. In a relatively homogeneous area, the volume density is independent of direction and relatively constant and has the following relationship with the average linear density [50]: Based on formulas (4) and (6), the volume density of the joint with the radius approximately complied with the negative exponential distribution can be obtained as where λ d is the average linear density of joints, λ V is the volume density of joints, d is the diameter, and r is the average radius of joints. According to formulas (4) to (7), the parameters of joint disc were calculated, as shown in Table 3.
e Monte Carlo method is a good tool for joint network modelling, which can approximately solve uncertain problems with a series of random numbers. In this paper, combined with the Monte Carlo method, the built-in language Fish provided by PFC is used to construct a threedimensional joint network model. e random numbers of joint disc parameters in accordance with corresponding probability distribution are generated in a 10 m × 10 m × 10 m cube space. Based on the OPEN GL technique, the joint network model in PFC is shown in Figure 6. en, the reliability of this model is tested by comparing the distribution of main joints in the effective area between the two-dimensional section of the model and the window of field survey. Figure 7 shows the joint distribution from the field survey and numerical model. e three-dimensional joint network model has a good statistical similarity with the real field distribution in attitude, scale, and density. e joints obtained from numerical models can well describe the real distribution characteristics of joints.

Calibration of the Microparameters
e input microparameters in PFC software cannot be measured directly through conventional laboratory tests, which also have highly nonlinear relationships with the parameters such as Young's modulus, uniaxial compression strength, and Poisson's ratio in continuum numerical simulation. erefore, parameter calibration is an essential part of simulation process in PFC to ensure that the microparameters can represent the macromechanical properties of rock well.
Before the calibration process, a series of conventional laboratory tests are carried out on Shuichang iron mine granite to obtain the mechanical parameters.
e mechanical test equipment is shown in Figure 8. ere are two sizes of experimental sample: 50 mm diameter and 100 mm height for compression tests and 50 mm diameter and 25 mm height for Brazilian tests, which are in good agreement with the sample sizes recommended by ISRM. Uniaxial compression tests and Brazilian tests are conducted at loading rates of 0.5-1.0 MPa/s and 0.1-0.3 MPa/s, respectively. Triaxial compression tests were stress path control with a constant loading rate of 0.5 MPa/s, and the confining pressures were 10, 15, 20, and 25 MPa. In addition, direct shear tests were conducted on samples to obtain the joint parameters. e sample size is 10 mm × 10 mm × 10 mm, with a constant horizontal shear rate of 0.01 m/s under different normal stresses of 300, 500, and 700 kPa. e mechanical parameters of the granite obtained from experimental tests are given in Table 4.

Calibration of the FJM.
As is known, calibration of the microparameters in PFC is a complicated process because it is impossible to describe the relationship between macroand microparameters with quantitative mathematical relation. e basic method to determine the microparameters is "trial-and-error." e specific method is as follows: by keeping other parameters constant, investigate the effect of one single parameter on simulation results, compare the macroparameters from laboratory tests results, and then repeat the above steps time and again until the parameter could be in good agreement with the observed macroscale response.   In this paper, the microparameters of the FJM in PFC, which affect the macromechanical properties of numerical sample, should be determined as shown in Table 5. ose microparameters are calibrated to match Young's modulus  Based on previous research, there are six microparameters that have significant effects on the numerical sample in the FJM, i.e., E c , k rat , μ, σ b , c b , and Φ b . erefore, a cylindrical sample model with a 50 mm diameter and a 100 mm height, which is in agreement with the rock sample size recommended by ISRM, is generated to simulate laboratory test (Figure 9). e sample model consisted of 14036 balls with radius varying from 1.00 to 1.66 mm.
In the FJM, the calibration process is time-consuming due to no one-to-one correspondence existing between macroparameters and microparameters. For example, the UCS, Poisson's ratio, and cohesion occur in varying degrees of change with the change of microparameter c b . To improve the efficiency and avoid adjusting parameters blindly, a L 25 (5 6 ) orthogonal design test and range analysis were used to investigate the response law of macroparameters and microparameters. e orthogonal test table and calculated macroparameters are given in Appendix (Tables 12 and 13). e results of range analysis are shown in Table 6, where R is the range value, to represent the influence of the microparameter on macroparameters, and a higher value indicates a higher ranking.
On the basis of aforementioned analytical results, a new improved trial-and-error method with high efficiency is used to calibrate the FJM. e final microparameters obtained are given in Table 7.
en, the microparameters were verified by using the numerical test. Figure 10 illustrates the stress-strain curves of uniaxial compressive test and failure modes from experimental result and numerical result. e stress-strain curves, both before and after the peak, are consistent with those of the experiment, and the failure patterns from the experiment and numerical test are quite similar. e comparison of basic mechanical parameters between the experimental and numerical results is shown in Table 8. e errors of Poisson's ratio, uniaxial compression strength, and tensile strength are 3.83%, 2.31%, and 3.16%, respectively. Young's modulus and internal friction angle are just 0.94% and 0.53% larger than the experimental results, respectively. e value of cohesion from simulation yields the largest errors, which is 9.64% higher than the experimental value, but that is still less than 10%. Obviously, the properties obtained from the numerical tests with calibrated parameters are in good agreement with those from experiment.

Calibration of the SJM.
e macroparameters of the SJM are usually determined by simulating the direct shear test sample of rock joint. In this paper, an unbonded SJM is selected, that is, the parameter "sj_state" is set to 0. As in the case of FJM calibration, a model with the same size as the experimental sample was generated in cube shear box, and a suitable set of microparameters of SJM was obtained based upon numerical direct shear tests matching the laboratory test response of the rock, as shown in Table 9.
e shear stress-shear displacement curves under different normal stresses from experimental and numerical direct shear tests are shown in Figure 11. e cohesion and internal friction angle of SJM can be calculated by using strength envelope according to those curves. e values of these two parameters are 30.58 kPa and 20.2°, which are 9.88% higher and 4.08% lower than the experimental value, respectively. Clearly, the calibration parameters can represent the macromechanical properties of rock mass well.

Scale Effect.
Due to the randomness in the distribution of joints, scale effect can substantially influence the rock strength, namely, the strength of a region decreases with increase in region size up to the point at which a representative size is reached. Compared with the engineering rock mass, the laboratory samples are usually small and do not contain systematic joints which affect the rock strength. For this reason, the microparameters of engineering rock mass in the SRM model should be considered for the potential impacts of joints and should be calibrated to match a strength that reflects the size of a typical rock mass, instead of a core sample. Hence, a concept named "representative elementary volume (REV)" is introduced, which is the minimum size with which the rock mechanical properties can be treated as equivalent continuous [51,52].
In order to investigate the scale effect and REV size of rock mass, seven SRM models with different sizes were established to conduct numerical uniaxial compression tests. ese cubical models with side lengths of 1, 2, 6, 10, 14, 16, Normal-to-shear stiffness ratio of both particle and bond μ Friction coefficient - Friction angle°A dvances in Civil Engineering and 20 m, respectively, are shown in Figure 12. Moreover, the geometric centers of these models were fixed to make calculated results comparable. e calibrated microparameters in Tables 7 and 9 were used to represent macromechanical properties of rock. Figure 13 shows the calculated stress-strain curve of multiscale models. As expected, the strength parameters decrease with the increase of model size, showing obvious scale effect. e results are presented in Figure 14.     Microparameter Subsequently, a series of sample models were established based on the REV size to conduct numerical tests. Making use of the method mentioned in Section 4.1, the microparameters of SRM were calibrated to match the mechanical parameters of engineering rock mass. e final microparameters of engineering rock mass in SRM are presented in Table 10.
e comparison of basic mechanical parameters between the Hoek-Brown method and numerical results is shown in Table 11. While the errors of the result are limited within 10%, it is considered to be a suitable set of microparameters to represent the macromechanical properties of engineering rock mass well.

Engineering Background.
e field joint survey was carried out on the north side of the open-pit mine slope, No. II mining area of Shuichang iron mine. Rock mass quality in this area is declined due to the well-developed discontinuous nature of rock such as fault and joint, causing local instability of slope. Records show multiple rock fall and landslips occurred in this area, especially at elevation 44m-116 m. A 3D numerical model is developed using the PFC software in order to investigate the micromechanism of failure of jointed rock slope during excavation process.

Construction of Jointed Rock Model.
e model takes the direction of perpendicular to the stope as the x-axis, extension direction of slope is y-axis, and vertical direction is zaxis. e lengths along the axis of x, y, and z are 140 m, 20 m, and 90 m, respectively. For obtaining reliable simulation results, the ratio of smallest characteristic length for the model to the median particle radius should be 50-100 [33,34,[53][54][55]. Considering the computational efficiency and accuracy, the minimum particle radius and particle size ratio are set as 0.3 m and 1.66, respectively. e four groups' dominant orientations of joints obtained in section 3.1 were added into the model successively. e model before excavation is composed of 81886 particles and 54219 cracks, as shown in Figure 17. e height and slope angle of bench in  Mechanical property  erefore, the number, mode, and propagation of crack can be used to analyze the failure evolution. In this paper, the model is excavated by four steps, one bench at a time. Subsequently, the failure mechanism during the excavation process was investigated from micromechanical viewpoint based on PFC software. Figure 18 illustrates the relationship between crack number and timestep. It should be noticed that generation of cracks is concentrated in a short time (about 4000 timestep) after excavation, which indicates that excavation may cause instantaneous failure of slope. Moreover, with the progress of excavation, the number of cracks is increasing to a total of   Hoek-Brown Mohr-Coulomb      14 Advances in Civil Engineering

Advances in Civil Engineering
Step 1 Step 2 Step 4 Step 3 Figure 19: Evolutionary process of microcrack. Advances in Civil Engineering 14800, and the number of tensile cracks is much more than that of shear cracks at any timestep. It can be inferred that the failure is essentially caused by the dominated tensile microcracks. Figure 19 shows the evolutionary process of microcracks during excavation, where the tensile crack and shear crack are shown in red and blue, respectively. In the first two steps, a small number of cracks appeared around the toe of the slope. From the third step, cracks begin to increase sharply in number and expand constantly inside the slope. At the end of the excavation process, the distribution of cracks reaches a depth of about 50 mm, and the tensile crack almost propagates through from the middle of the third bench to the toe of the fourth bench, forming a macro fracture zone.
Diagrams of velocity vector ( Figure 20) and displacement vector ( Figure 21) clearly show the distribution of potential unstable rock mass when the slope is unstable. e distribution is the same as that of microcrack. e velocity and displacement around the third and top of the fourth bench are larger than other regions. Figure 22(a) shows the horizontal displacement of particles at the end of the simulation. Particles with large displacement are mainly concentrated in the middle and lower part of the slope, some of which have been separated from parent rocks and form a rock fall. And due to the cross cuttings of joints, the main failure mode of slope is wedge failure, being relatively complete in the upper and seriously damaged in the lower of the whole slope, which is consistent with the actual situation, as shown in Figure 22(b).

Conclusions
In this paper, the SRM model was established in PFC software to represent the mechanical properties of rock mass from Shuichang open-pit iron mine. Based on the prepared model, the failure mechanism during excavation process was investigated from the micromechanical viewpoint. e conclusions of this study are as follows: (1) statistical analysis method can act as a useful tool to group the dominant orientations of joints and determine the geometric parameters such as dip angle, size, and trace length. e joint distributions of Shuichang slope were obtained by the field survey and quantitatively described in the numerical model. (2) e effects of different microparameters on macromechanical properties of rock are different. rough the orthogonal experiment and variance analysis, the calibration and optimization of the parameters were accomplished. e mechanical properties including stress-strain curve and failure pattern from numerical results were verified with those from experimental results.
(3) e existence of joints results in scale effect and anisotropic behavior of rock mass, and these properties tend to gradually weaken with the increasing rock block size. By carrying out numerical tests on multiscale SRM models, the representative elementary volume (REV) in the selected research area was obtained with the size 16 m × 16 m × 16 m. en, the microparameters of the SRM model were calibrated to match the mechanical parameters of the engineering rock mass. (4) SRM is an effective method to analyze the failure evolution of jointed rock slope. e failure of the slope was dominated by the tensile microcracks between bonded particles. e microcracks primarily occurred at the bottom of the slope and gradually developed upwards. In addition, microcracks were mainly distributed on the shallow part of the slope. After excavation, the wedge occurred in the middle and bottom part of the slope (Appendix) (Tables 12 and 13).  16 Advances in Civil Engineering