DEM analyses of shear behaviour of rock joints by a novel bond contact model

The failure of rock joints is one of the potential causes for the local and general rock instability, which may trigger devastating geohazards such as landslide. In this paper, the Distinct Element Method (DEM) featured by a novel bond contact model was utilized to simulate shear behaviour of centre/non-coplanar rock joints. The DEM results show that the complete shear behaviour of jointed rock includes four stages: elastic shearing phase, crack propagation, the failure of rock bridges and the through-going discontinuity. The peak shear strength of centre joint increases as the joint connectivity rate decreases. For intermittent noncoplanar rock joints, as the inclination of the rock joints increases, its shear capacity decreases when the inclination angle is negative while increase when positive. Comparison with the experimental results proves the capability of this DEM model in capturing the mechanical properties of the jointed rocks.


Introduction
Shear failure of rock joints can induce the instability of rock mass, a common medium in the field of mining, oil exploitation, water conservancy and hydropower. Thus, the shear behavior of rock joints becomes a topic of the geological engineering. Lajtai [1][2] conducted direct shear tests to examine and analyze the strength of discontinuous rocks, and concluded that the total shear strength is composed of two parts, namely, the joint friction along the separated parts of the weakness plane and the fundamental shear strength and internal friction in solid bridges. Gehle and Kutter [3] investigated the shear behavior of intermittent rock joints, arguing that the failure pattern and the shear resistance of jointed rock are significantly affected by both the orientation of the joints and the normal stress. Cao et al. [4] further found that the rock specimen takes a turn from wing crack propagation failure to crack coalescence failure with the increase of the angle of the rock bridge.
Although theoretical and experimental studies have given some promising results of the shearing of rock joints, some very important microscopic information controlling the failure patterns still remain unknown. More and more researchers have used numerical methods to further investigate this problem [5][6][7][8][9][10][11]. Among these methods, the Distinct Element Method (DEM) proposed by Cundall and Strack 5 To whom any correspondence should be addressed.  [8] is a feasible method to study the shear behavior and failure mechanism of sand and rocks [9][10][11][12]. Park et al. [9] examined the effects of the geometrical features and the micro-properties of a joint on its shear behavior by using a bond model. Potyondy and Cundall [10] proposed a numerical model for rock and investigated the sensitivity of the results to micro-properties. Jiang et al. [11] investigated the shear behavior and strain localization in cemented sands by 2D DEM. More recently, Bewick et al. [12] investigated mechanism of rupture by comparing DEM direct shearing result with experimental one under different normal pressures.
In this study, a new DEM bond model proposed by one of the authors has been adopted to simulate mechanical response of jointed rock mass under direct shearing condition. The shear curves and shearing capacities of the DEM virtual rock joints are qualitatively compared with the experimental results. Also, the force chains inside the jointed rock sample are presented to provide a view in the shear behaiviour microscopically.

DEM bond contact model
The DEM bond contact model is proposed by Jiang et al. [13] for the purpose of simulating heavily cemented materials such as granite. The bond contact model was first derived theoretically [14][15][16] and then further verified by laboratory tests on bonded aluminum rods [17][18].
The total bond contact is assumed to be composed by two paralleled parts: inter-particle and interbond contacts. Contact forces are assumed to be transferred through not only the bond but also the inter-particle contact: (1) where F bn , F bs and M b are the normal force, the shear force and the moment transferred through the bond, respectively, while F pn , F ps and M p are the corresponding components transferred through the inter-particle contact. The inter-particle contact laws can be approached in reference [15] and the mechanical response of inter-bond contact can be expressed as: (4) where min[•] is the operator taking the minimum; u n , u s and θ are the overlap, the shear displacement and the relative rotation angle, respectively; k bn , k bs , k bm are the normal, shear and rolling contact stiffness of bonds, respectively; R bc , R bt , R bs and R br are the compressive, tensile, shear and rolling bond strengths, respectively.
The projections of the strength envelope in F s -F n , M-F n planes and critical normal force R cr can be expressed as: where f=(F bn +R bt )/(R bt +R bc ) is the stress ratio, which can be used to define the types of bond failures: when f = 0, compressive failure occurs; when 0 < f < 1, tensile-shear-torsional failure (F n < 0) or compressive-shear-torsional failure (F n > 0) occurs; when f = 1, tensile failure occurs; μ b and μ p are friction coefficients of bonds and particles; β b and β p are rolling resistance coefficients of bonds and particles; f s , f r , g s and g s are the envelope shape factors; is the common radius of two particles in contact and k pn is inter-particle normal contact stiffness. It should be noted that the strength in Eqs. (5) , , bn n bc n bn bn n bt n and (6) includes both of the components of bond and inter-particle contacts when F n > 0, which means that when F n < 0 the second term disappears as inter-particle contact cannot transfer tensile force. If forces reach strength envelope/critical values, they will fall to the residual state. The residual shear strength R sr and rolling strength R rr can be expressed by: which also contain two terms corresponding to bond and inter-particle components.

DEM modelling
The DEM rock sample containing 10,000 particles was generated by Under-compaction Method (UCM) proposed by Jiang et al. [18], with maximum, minimum and average diameters being 2.0, 0.5 and 1.3 mm respectively. The size of the specimen is 12.23 cm×7.82 cm with planar void ratio being 0.20. The synthetic intact rock material has already been calibrated by comparison with uniaxial compression and Brazilian test results of granite [19]. A set of joints in the given arrangement are then generated within the rectangular DEM sample, and the joints were identified by as a 2 mm thick layer of particles including no bond effect. The assigned parameters for both intact rock and joints are listed in Table 1.
During the direct shear test, shear load was applied to the DEM jointed rock sample by moving walls #2 and #3 towards each other at a quasi-static velocity while maintaining walls #1 and #4 under a constant normal pressure. For the centre joints, the influence of joint connectivity rate (e.g., 0.6, 0.8, 0.9) is investigated with normal stress being 1.0 MPa, while for intermittent non-coplanar joints with a given joint length l and spacing e, the effect of the inclination of joints i (varying from -75° to 90°) is explored.   Figure 3 illustrates the relationships between the shear stress/number of bond breakages and shear displacements for the specimen with g = 0.6, it shows that the deformation and failure process of DEM specimens can be divided into four stages: (1) elastic shearing phase: shear stress increases almost linearly with the increase of shear displacement, while the no bond breakages happen during this shearing phase; (2) crack propagation: as the shear displacement increases, shear stress fluctuates until reaching the peak stress, and a large number of bond breakages occur with tensile-shear-torsional failures dominating in number; (3) the failure of rock bridges: the shear stress quickly drops to a residual value after arriving the peak; (4) through-going discontinuity: the shear stress remains at a relatively low residual value in spite of some fluctuations, during this period of shearing, the number of the microscopic bond breakages consistently grows. The shear behaviour above is quite typical compared with that normally obtained from experimental results [20]. Figure 4 shows the relationships between the peak shear strength and joint connectivity rate, the shear strength in numerical and experimental results are normalized by their own peak stresses for the case with g = 0.6 separately. It can be observed that both of the peak shear stresses significantly decrease with the increase of normal stress, and the trends of the DEM and experimental results are similar. This phenomenon can be reasoned that growth of the proportion of joint weakens the shearing capacity within the shear axis, causing the decrease of the shear strength of the jointed samples.

Intermittent non-coplanar joints
Considering that all the cases with different inclinations share almost the same trend, only the shear curve of non-coplanar intermittent rock joints with inclination i = 45° is showed in Fig. 6. The shear curve of the non-coplanar joints is quite similar to that of the centre ones, except that no apparent yielding when approaching the peak stress. The reason for this might be that the number of the bond breakages rises exponentially when the jointed sample arrives yielding state, so the stress drops sharply simultaneously, but not the case for the centre joint. Normalized peak shear stress Joint inclination Experiment [3] DEM Figure 5 Relationship between shear stress/broken bond number and shear displacement by DEM simulation with joint inclination angle of 45° Figure 6 The relationship between the joint angle and the normalized peak value of shear stress Figure 6 illustrates the shear loading capacities of non-coplanar intermittent joint assemblages with different inclinations. The numerical and experimental results are both normalized by the maximum peak stresses obtained in experimental tests (i.e., -30° for experimental results). The results show that in the DEM simulation, the peak stress decrease as the inclination angle changes from -75° to 15°. It reaches its lower limit value at around 15° and then increase till 75°. The experimental results show almost the same trend for inclination angles varying from -30° to 75° except for the values near -45°. The difference can be partially explained that in the DEM simulation, the joints are filled with particles, whereas in the laboratory, the man-made joints are just open cracks with nothing inside.

Conclusions
The numerical simulation results indicate that the complete shear behaviour of jointed rock includes elastic shearing phase, crack propagation, the failure of rock bridges and the through-going discontinuity. For centre rock joints, the decrease of the joint connectivity rate will enhance the peak shear stress. As to intermittent non-coplanar rock joints, the relationship between shear stress and shear displacement presents the elastic-brittle-plastic behavior. The shear capacity drops when the joint inclination angle is negative but rises when positive with the increases of rock joints inclination. The failure process of intermittent non-coplanar rock joints begins with the transmission of compressive forces through the tips of rock joints and ends with the breakage of rock bridges. The comparison between DEM and experimental results proves that this bond contact model is capable of describing the shear behaviour of jointed rocks during direct shearing.