Advances in Civil Engineering
Volume 2018 (2018), Article ID 4294501, 19 pages
https://doi.org/10.1155/2018/4294501

Research Article

Experimental and Numerical Study of the Failure Behavior of Intermittent Rock Joints Subjected to Direct Shear Load

Xiang Fan,1,2,3 Kaihui Li,⁴ Hongpeng Lai,¹ Qihua Zhao,² and Zhenhua Sun1,5

¹School of Highway, Chang’an University, Xi’an 710064, China
²State Key Laboratory of Geohazard Prevention and Geoenvironment Protection, Chengdu University of Technology, Chengdu 610059, China
³Key Laboratory of Bridge Engineering Safety Control by Hunan Province, Department of Education, Changsha University of Science & Technology, Changsha 410114, China
⁴Department of Civil and Environmental Engineering, The Hong Kong Polytechnic University, Hong Kong, China
⁵Shaoxing Transportation Investment Group Co. Ltd., Shaoxing, Zhejiang 312000, China

Correspondence should be addressed to Qihua Zhao; zhqh@163.com

Received 3 March 2018; Revised 15 May 2018; Accepted 23 May 2018; Published 10 July 2018

Academic Editor: Wei Wu

Copyright © 2018 Xiang Fan et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Abstract.

Two series of intermittent rock joints containing three joints arranged along the central shear axis were considered in this study. The failure behavior under direct shear loads was investigated by means of both physical tests and numerical simulations. The cracking behavior was found to be distinctly associated with the joint arrangement. Several types of main and secondary cracks were identified. The variation trends of the crack initiation stress ratio with inclination angle were analyzed and found to be partly different for the two series of intermittent joints. The whole fracturing process was characterized by three phases. Not all samples have to experience all three phases. The second phase is alternative and can be reflected by the shearing curve. Hence, two types of shearing curves, including single and double peaks, were identified. The double peak is due to the extrusion or sawteeth cutting in the second phase. Moreover, the numerical micromechanical analysis was performed to explain the shear behavior using the contact force and microcrack within the specimen. Based on the numerically measured local stresses, maximum and minimum principal stresses around the middle joint at crack initiation stress and peak shear stress were analyzed.

1. Introduction

Discontinuous joints are widely known to exist within natural rock masses. Types of joints play a dominant role in the mechanical behavior of jointed rock such as the shear or compressive behavior. In the past decades, many approaches including the laboratory test, numerical simulation, and theoretical analysis, have been conducted to get a better understanding of the failure behavior of jointed rocks subjected to various loadings [1–7]. Deep knowledge of rock mass failure characteristics benefits the engineering construction, for instance, mining roadway, tunnel, hydropower dam, and so on, so the relevant investigations still attract intense interest.

Recently, for the discontinuous joint (intermittent joint) in the solid rock under a compressive loading, the failure behavior has attracted much attention in the literature [8–14]. The laboratory test is known to be the most widely used method to investigate cracking behavior. In general, open or closed joints were precast in rock-like materials by inserting thin sheets [8–12]. In some studies, the hydraulic cutting technology was adopted to produce open joints in the real rock [15]. For a rock block containing intermittent joints, the geometrical parameters, joint orientation, and length are considered as factors directly influencing cracking. Abundant observations show that the crack often initiates at the joint tip, and the crack initiation, propagation, and coalescence are dramatically dependent on the joint arrangement [8, 9, 15–20]. Tensile and shear cracks are two basic crack initiation modes, but the propagation and coalescence trajectory is influenced by the properties of neighboring joints such as joint spacing and the persistent ratio. Many kinds of combinations of influential factors incorporating joint length, orientation, spacing, and number were further considered. Therefore, the majority of crack patterns and failure modes were derived from different cases. In addition to the summary of cracking types, the relevant mechanism was also theoretically or numerically revealed [21, 22].

Obviously, shear and compression are two quite different ways of loading. The shear plane mainly bears the shear load in shear, while the whole block bears the compressive load in compression. For the rock block with the same joint arrangement, the observed failure behavior under shear loading may be different from the failure behavior under compression. During shearing, the failure characteristic is dependent mainly on the joints close to the shear plane. Moreover, cracking features are related to the joint geometrical parameters. Using a specially designed shear device, Gehle and Kutter [23] conducted a series of direct shear tests on intermittent joints, and three phases, including crack initiation, friction, and dilation sliding, were identified. The shear failure behavior of the intermittent joint depends significantly on the joint orientation. Gerolymatou and Triantafyllidis [24] prepared gypsum specimens with different joint orientations and carried out a large-shear-displacement direct shear test, they found that the shear strength as well as the crack pattern is influenced by both joint spacing and orientation. As stated for the shear failure of the intermittent joint, the research focuses mainly on the cracking types and the relationship with joint orientation. Even though the visual cracks were identified, the internal stress distribution surrounding a joint was not well exhibited under shearing. Therefore, a more comprehensive understanding of the intermittent joint under shearing is necessary. But it is quite difficult to derive the internal stress for a physical test. However, a numerical approach covers the shortage.

Numerical approaches are often used to further explain the physical failure phenomenon. To date, continuum methods, discontinuum methods, and hybrid continuum/discontinuum methods have mostly been used to simulate the fracture behavior of rock masses. The differences between these methods have been discussed by researchers [4, 25, 26]. Continuum methods cannot present explicitly the crack initiation, propagation, and nucleation because mesh grids in the problem domain should remain in the same neighborhood during the whole deformation or transport process and thus block rotation so that opening of the mesh will not happen. In contrast, discontinuum methods are able to model zone detachment, large displacements, and block rotation, and therefore, the fracturing process in rock or rock-like materials can be presented explicitly by microcracks formed in models based on the discrete element method (DEM) [5, 27–32]. Additionally, based on the rock failure analysis process, the numerical direct shear test was performed to investigate the shear behavior of intermittent joints with different geometrical parameters. Meanwhile, the whole shear failure process was visually recorded, which is the greatest advantage relative to the physical test [33]. Among the developed numerical approaches, particle flow code (PFC) is useful. Using the particle flow code approach, the internal stress evolution could be well exhibited, especially the internal stress evolution surrounding an open joint. Moreover, the magnitude of the internal stress can also be derived and can be used to account for the cracking mechanism. According to the published reports regarding the shear failure of intermittent joints, particle flow code has not yet been used to compute the principal stresses around an open joint in shear. Additionally, macrocrack propagation is composed of microcracks. However, the microcracks are impossible to identify in a physical test.

As stated above, the shear failure behavior of the intermittent rock joint has not been investigated currently relative to the compressive failure behavior. To further explore the shear failure mechanism of intermittent rock joints, physical and numerical tests both were conducted in the study. Two series of intermittent joint were molded using gypsum material by changing the joint arrangement along the shear axis. With the aid of the bonded particle model, the increase in microcracks was related to macrocrack propagation, and the internal stress distribution around the joint was computed to explain the shear failure mechanism.

2. Physical Test and Numerical Simulation

2.1. Joint Arrangement

In the present work, two series of joint configurations were considered. The dimensions of the prismatic specimen were 150 mm × 150 mm × 30 mm (height × length × width). Three joints were arranged within each sample along the shear plane. To better describe samples hereafter, the two series were named A and B. In Series A illustrated in Figure 1(a), two joints were located at the two lateral edges and kept the inclination angle of 0° fixed while the middle one was located at the sample centroid. The middle joint inclination angle, βA, rotating in the clockwise direction, was set from 0° to 150° with an interval of 30°. In Series B illustrated in Figure 1(b), all three joints rotated along the joint center in the counterclockwise direction. The inclination angle, βB, was also set from 0° to 150° with an interval of 30°. Each joint length (L) was set to 20 mm. The shear directions for Series A and B are shown in Figures 1(a) and 1(b), respectively. The ligament connecting the two tips of each neighboring joint is called the rock bridge. The detailed joint configuration can be found in Table 1.

(a)

(b)

Figure 1: Schematic diagram of the joint configuration for (a) Series A and (b) Series B.

Table 1: Joint inclination angle and sample number.

Series A
Series B

βA
Sample number
βB
Sample number

0°
A0
0°
B0

30°
A30
30°
B30

60°
A60
60°
B60

90°
A90
90°
B90

120°
A120
120°
B120

150°
A150
150°
B150

2.2. Sample Preparation

Gypsum is a brittle material like rock material, so it is usually used as a rock-like material to model the behavior of rock specimens [9, 10, 23, 24]. First, gypsum powder and water were mixed at a ratio of 1 : 0.35 in weight and then stirred well. Next, the liquid gypsum was quickly cast into steel molds and vibrated. The liquid gypsum often consolidated in less than twenty minutes. For jointed specimens, the mica sheets should be inserted according to the designed location and inclination angle. Since gypsum is the air setting material, the samples do not need a water spray cure.

2.3. Physical Test

A direct shear test was performed to investigate the failure characteristics of the jointed rock. The test was conducted on the shear-rheology instrument RYL-600 shown in Figure 2. This instrument consists mainly of the control system, shear loading system, normal loading system, and input/output devices including a computer and printer. The jointed sample was placed in the shear box first, and then the normal stress was applied at 2.5 MPa with a speed of 50 N/s. Then, the shear loading was applied at a speed of 2 mm/min.

Figure 2: Testing instrument.

2.4. Numerical Simulation

Among the numerical approaches mentioned in Section 1, DEM simulations offer another merit of allowing microscale behavior such as contact forces, particle displacement vectors, and internal stresses to be obtained. Considering these inherent advantages of simulating the fracturing process in rock or rock-like materials, Particle Flow Code (PFC) was used to explore the rupture mechanism of rock-like specimens with joints in shearing. The parallel-bonded particle model (PBPM) was employed in this study, which has been extensively used and validated in simulating the fracturing process of rock or rock-like materials in a flood of literature reviews [5, 16, 22, 28, 29]. Since the microproperties of the real physical material are not determined directly in the laboratory, calibrating the PFC synthetic specimen requires an iterative trial-and-error procedure. The microparameters for the PFC model were determined by conducting uniaxial compression and Brazilian numerical simulations on synthetic intact specimens without any preexisting flaws or joints. As the specific calibration procedure is rather mature, it will not be introduced in this paper. Please refer to the PFC2D manual for more details concerning microparameter calibration [34]. The model porosity is not equivalent to the porosity of real intact rock attributed to the utilization of circle disks to represent the real grains [35]. Moreover, the mode of tensile or shear failure of a contact is dependent on the applied local stress exceeding either the tensile or shear strength of the contact. After an iterative procedure, all the microparameters were determined as listed in Table 2. The PFC synthetic intact specimen is composed of 43,298 balls and 97,872 parallel-bonded bonds. The macroproperties of the PFC intact specimen are comparable to those tested by the physical experiment as listed in Table 3.

Table 2: Microparameters in bonded particle model.

Microparameters
Values

Minimum particle radius, Rmin (mm)
0.3

Particle size ratio, Rmax/Rmin
1.66

Particle bulk density, ρ (kg/m³)
1560

Particle contact modulus, Ec (GPa)
0.6

Ratio of particle normal to shear stiffness, kn/ks
2.0

Particle friction coefficient, μ
0.5

Parallel-bond modulus, Ec (GPa)
0.6

Ratio of parallel-bond normal to shear stiffness,
2.0

Parallel-bond normal strength, mean σm (MPa)
9.0

Parallel-bond normal strength, standard deviation σsd (MPa)
2.0

Parallel-bond shear strength, mean τm (MPa)
5.4

Parallel-bond shear strength, standard deviation τsd (MPa)
1.0

Table 3: Comparison of experimental and numerical results.

Experimental
Numerical

Uniaxial compressive strength (MPa)
10.63
10.62

Young’s modulus (GPa)
0.89
0.89

Poisson’s ratio
0.25
0.246

Tensile strength (MPa)
1.58
1.77

The numerical direct shear setup (Figure 3(a)) is a simplification of the experimental direct shear setup in the laboratory. The gap between the upper and lower parts of the shear box is 2 mm. All of the contacts between the wall and particle are frictionless. First, an initial velocity was exerted on the top wall of the shear box (Wall 2) until the normal stress reached the desired value. The normal stress was obtained through averaging the reaction forces acting along both the top and bottom walls of the shear box (Wall 2 and Wall 1) divided by the length of the synthetic specimen. Then, the shear displacement was applied to the upper part of the shear box (Wall 3 and Wall 4) by a fixed shearing velocity when the normal stress remained constant. To ensure that the whole shearing process was in a quasistatic loading condition, a shearing velocity of 0.006 m/s was adopted. The loading rate in the simulation is different from the loading rate in the physical experiment. For instance, the loading rate of 0.006 m/s in this situation means approximately 1.5 × 10⁻⁹ m/step because the timestep is 2.47 × 10⁻⁷ s. To be consistent with the experiment, the lower part of the shear box (Walls 1, 5, and 6) was kept fixed throughout the shearing. The average shear stress was calculated by the reaction forces exerted along Wall 6 divided by the specimen length. Contact forces denoted by force chains exist between particles. Compressive and tensile contact forces are marked black and green, respectively, as shown in Figures 3(b) and 3(c). Before loading, the compressive contact force exists basically uniformly inside samples, but the tensile contact force exists mainly near the sample edges. For the jointed specimens, the open joints were created by deleting particles located in the corresponding regions shown in Figure 4(a). The arrangement of joints is the same in the physical flawed specimens. Even though open joints were created, compressive and tensile contact force distributions (Figures 4(b) and 4(c)) are still identical to the compressive and tensile contact force distributions for the intact model, respectively. No contact force concentration is observed around joints.

Figure 3: Intact numerical model: (a) schematic diagram of direct shear setup; (b) compressive contact force chain; (c) tensile contact force chain.

Figure 4: Jointed numerical model (A30): (a) intermittent joint model; (b) compressive contact force chain; (c) tensile contact force chain.

3. Analysis of Shear Failure Behavior

3.1. Category of Crack Patterns

From the published reports [1, 20, 23, 28, 30], we know that in the condition of compressive loading, the crack often initiates at the joint tips, and the crack initiation type is related to the joint orientation. The loading modes of shear and compression are quite different. Under shear loading, the shear plane is subjected mainly to the shear stress. For coplanar and non-coplanar joints on the shear plane, internal stresses around joints resulting from the shear loading are different, so that cracking behaviors are different. In the study, a main crack (MC) and secondary crack (SC) were observed. The main crack is the category of crack connecting the tips of two neighboring joints. The SC is the kind of crack initiating at the joint tips before or after the peak shear stress (τp). The SC seldom connects neighboring joints except in B120 and B150 cases. Both Samples A0 and B0 can be regarded as the coplanar intermittent joints. The joint orientation is consistent with the shear direction as well as the shear plane. Therefore, cracks are first initiated at the tips and then easily propagated along the shear plane. This MC is named MC-I, as shown in Table 4. In this case, since the direct shear failure leads to the rock bridge rupture, the generated cracks are shear cracks, and no SC was observed. From 30° to 150°, one (Series A) or three joints (Series B) are inclined to the shear direction, and the rock bridge is longer relative to the rock bridge of 0° samples. The cracking is quite different from the cracking in 0° cases. In Series A, since only the middle joint orientation varied, cracks initiated from both tips of the middle joint and one tip of the edge joints, and then propagated oppositely and finally coalesced. In the case of the joint arrangement, as shown in Figure 5, cracks always propagated through A (solid line), not B (dashed line), that is, the cracking trajectory between two neighboring joints followed the shortest way. This MC type is named MC-II. In Series B from 30° to 150°, MCs always also initiated at the joint tips, and cracks obliquely propagated towards another tip of adjacent joints. For another two joints close to the sample edge, cracks initiated at the joint tip approaching the sample edge and propagated towards the sample edge (see MC-III in Table 4). If the shear loading continued after MCs completely coalesced, in some cases of Series A, an SC may be generated at the lower middle joint tip (SC-I) such as that shown in Figure 6(a). For B30 and B60, SCs initiated at the upper and/or lower joint tips and propagated in the direction approximately perpendicular to the joint orientation (Figure 6(b)). They are also regarded as SC-I. SC-I generation is due to mutual squeezing between upper and lower blocks. For B90, B120, and B150, SCs (SC-II) also initiated at the upper and/or lower joint tips, but the propagating direction of SC-II is approximately consistent with the shear direction that is quite different from the shear direction of SC-I. SC-II generation is due to the tensile internal stress acting at the joint tips before the peak shear stress, which will be analyzed in Section 3.6. Once the MCs have completely coalesced, regular interfaces with one sawtooth in Series A or several in Series B were created. The inclination angle of the sawtooth varies with the joint orientation. A detailed description of the failure of the sawtooth interface will be given in the following context. The MC cannot be simply regarded as a pure tensile or shear crack. Actually, the MC is a mixed mode, as shown in Figure 6(a). The shear dust and the tensile crack could be observed on one MC trajectory. The MC mode is related to the joint orientation. Certainly, two types of SCs are both generated at the sawtooth roots of the upper or lower block under the tensile stress, so all SCs are regarded as tensile. The mechanisms of SC-I and SC-II are different, which will be numerically explained in Section 3.4.

Table 4: Category of crack patterns.

0°
30°, 60°
90°, 120°, 150°

MC
Series A

Series B

SC
Series A

Series B

Figure 5: Schematic diagram of the main cracking trajectory in Series A.

(a)

(b)

Figure 6: Typical cracks observed in the physical test (T stands for tensile crack; S stands for shear crack). (a) A60. (b) B60.

3.2. Curve Types of the Shear Stress-Shear Displacement

The shear stress-shear displacement curve (shearing curve) reveals the variation tendency of shear resistance. In addition, the cracking process can be indicated partly by the shearing curve. For instance, the SC propagation may create a second peak shear stress. It will be explained in detail in the following context. Among all curves of shear stress versus shear displacement, two typical curve types are identified. As shown in Figure 7(a), the shear stress monotonically increases, with the shear displacement increasing before the peak value. After one certain shear displacement, the shear stress reaches the peak value. Actually, the shear displacement at the peak shear stress is dependent on the joint configuration. After the peak value, the shear stress drops rapidly and then enters into a long residual stage. Due to the continuous friction and sliding behavior on the rupture interface, within the residual stage, the shear curve is not smooth but slightly fluctuating. This described curve type has one single peak only, so it is named a single-peak curve (Type I). The single-peak curve was mostly found in the shearing process of Samples A0, A90, A120, A150, and B0. The other curve type (Type II) is shown in Figure 7(b). Like the single-peak curve prior to the peak, the shear stress also purely increased with the shear displacement increasing. After the first peak, the shear stress dropped approximately vertically to some degree. With new cracks generated or the extrusion on the sawtooth interface after MC coalescence, the shear stress rose up again after experiencing a short shear stress valley. Due to the second increase in the shear resistance, a second shear stress peak was generated. The second peak shear stress is as a rule lower than the first peak shear stress. After the second peak, the shear stress gradually comes into the fluctuated residual stage. According to the shear curve features of Type II, it is named a double-peak curve. The type II curve can be found during the shearing process of Samples A30, A60, B30, B60, B90, B120, and B150. This curve type has two peaks, which is the biggest difference between the two curve types. The difference between the two curve types denotes the fracturing difference. According to the crack patterns described above, after the MC coalescence between neighboring joints, the SC propagation or extrusion between the upper and lower blocks on the sawtooth face finally results in the second peak. The curve shape indirectly reveals the fracturing process during the shear test, which will be described in detail in the following section. A shear displacement that is long enough can bring about a complete shear curve.

(a)

(b)

Figure 7: Typical shear curves for (a) single peak and (b) double peak.

3.3. Crack Initiation Stress Ratio

Fracturing may be revealed by the shear stress-shear displacement curve to some extent, as stated above. In the numerical simulation, crack initiation and coalescence can easily be recorded. According to published reports, the crack initiation stress is associated with the joint configuration. In particle flow code, the crack initiation stress is defined as the stress when the number of microcracks reaches 1% of the peak stress, so the crack initiation stress is associated with the fracturing process. The crack initiation stress ratio is defined as the ratio of the crack initiation stress to the peak shear stress. The real crack initiation stress is very difficult to be accurately obtained. Instead, it is easy to be recorded using the numerical approach. Figure 8 shows that the crack initiation stress ratio varies with the inclination angle. The crack initiation stress ratio for 30° samples in both series is the lowest ratio, indicating that the most further shear stress is needed from the first crack appearance to sample failure when the joint inclination equals 30°. The trend of the crack initiation stress ratio for Series A and Series B is similar from 0° to 90°, decreasing first and then increasing. The trend is distinguished from the trend obtained in compression tests, which increase from 0° to 75° [34]. Thus, it demonstrates further that the fracturing behavior of jointed rocks is related to the loading condition. However, the variation trend is different between Series A and Series B from 90° to 150°. The total length of sheared rock bridges from the crack initiation stress to the peak stress may be contributed by this difference. For instance, in Series A, the MC propagation and coalescense are along the shortest rock bridge (rotating compared with that of samples with lower joint inclinations), and the total length of the sheared rock bridges decreases with the joint inclination increasing from 90° to 150°. In contrast, the total length of the sheared rock bridges for Series B increases with the joint inclination increasing from 90° to 150°, as shown in Table 4. Therefore, the crack initiation stress ratio for Series A rises while the crack initiation stress ratio for Series B is reduced from 90° to 150°.

Figure 8: Crack initiation stress ratio of Series A and Series B samples.

3.4. Phase Division of the Cracking Process

On the basis of the function of cracking recording in PFC, the cracking process is described based on numerical results using four typical fracturing behaviors. To clearly state the shear failure process, according to the fracturing features at different stages, four typical fracturing moments, that is, visual crack initiation (Point A), one MC coalescence (Point B), complete crack coalescence (Point C), and the residual stage (Point D), are presented to stand for the whole fracturing process. As mentioned above, particles in the numerical model are cemented by the parallel bond, which may break in the shear or normal direction when the local stress acting on the parallel bond exceeds the predetermined shear or tensile strength of the bond. Each broken bond creates one microcrack that finally forms visual cracks. Each microcrack is marked in red for tensile failure or blue for shear failure. The number of microcracks can be recorded in the shearing process.

According to the definition of crack initiation stress, cracks are not obviously distinguishable at the initiation stress. The cracks are assemblies of a few microcracks, and the cracking trend is hard to identify. Therefore, to more distinctly evaluate the initial crack, visual cracks are presented and the corresponding shear stress is marked on the shearing curve. A combination of Figures 9 and 10 shows the fracturing process of A30 at different shear stresses. When the shear stress was approaching 0.9 τp, it was observed that cracks first initiated at the tips of the middle joint and the other two edge joints (Point A, Figure 10(a)). Once the cracks started initiating, cracks grew rapidly and first coalesced at the left rock bridge from Point A to Point B (Figure 10(b)). Meanwhile, in the process of left rock bridge coalescence, the shear stress reached the peak. From Point B to Point C, the shear stress sharply deceased. Simultaneously, cracks went through the right rock bridge (Figure 10(c)). Figures 10(c)–10(e) show that cracks coalesced following the shortest trajectory illustrated in Figure 5. Comparing Figures 10(c) and 10(d), the numerical crack coalescence is consistent with the experimental crack coalescence. On the shearing curve, after cracks completely coalesced, the curve experienced a long shear stress valley, and then the shear stress gradually rose again. Complete MCs break the sample into upper and lower blocks with a single sawtooth (Figure 10(e)). Under the shear load, the sawtooth may be sheared off under a high normal stress or the upper block climbs up along the contact face under a low normal stress. In the present work, the strength of the modeling material, gypsum, is quite low. For A30, under a high normal stress, despite climbing is limited, when cracks coalesced completely, the upper block climbed a little, so the normal displacement increased as shown in Figure 9. During climbing, the sawtooth face is extruded, which creates a second peak shear stress. Due to extrusion, the SC (Figure 10(f)) may be generated at the lower tip of a middle joint like A60 shown in Figure 6(a). Before Point A, the number of microcracks increased a little, but from Point A to Point C, the number of microcracks increased rapidly. After Point C, the number of microcracks increased slowly. Failure of A60 is analogous to the failure of A30.

Figure 9: Shear curves of Sample A30.

(a)

(b)

(c)

(d)

(e)

(f)

Figure 10: Fracturing process of Sample A30. (a) Point A. (b) Point B. (c) Point C. (d) Real failure. (e) Single sawtooth. (f) Point D.

As shown in Figures 11 and 12, cracks first initiated at approximately 0.8 τp (Point A, Figure 12(a)). The crack initiation locations are identical to the crack initiation locations of A30. The MC first coalesced at the left rock bridge from Point A to Point B (Figure 12(b)), within which the shear stress reached the peak value. Cracks completely coalesced from Point B to Point C, that is, the rapid shear stress declining stage (Figure 12(c)). Comparing Figures 12(c) and 12(d), the numerical crack coalescence is in good agreement with the experimental crack coalescence. As Figures 12(c) and 12(d) show, the crack propagation still followed the shortest trajectory. The MC coalescence separated the sample into upper and lower blocks and created an interface with one sawtooth (Figure 12(e)). Since the upper block moved away from the middle joint, the joint gap got wider and wider with the shear displacement increasing. Compared with Figure 10, since no climbing took place, the normal displacement did not keep increasing like Sample A30. The peak normal displacement may be resulted from that the upper block slide of the irregular interface. Namely, the interface between the upper and lower blocks is composed of particles, the interface is just a little fluctuate because of the arrangement of unequal radii particles. Therefore, no second peak and extrusion on the sawtooth face were observed. The variation trend of Sample A150 is like the variation trend of A30. Before the visual crack initiation, the crack number increased slowly. Cracks sharply increased from initiation to crack coalescence. During the residual stage, the microcrack number increased slightly. Failure of A120 and A90 is analogous to the failure of A150.

Figure 11: Shear curves of Sample A150.

(a)

(b)

(c)

(d)

(e)

Figure 12: Fracturing process of Sample A150. (a) Point A. (b) Point B. (c) Point C. (d) Real failure. (e) Single sawtooth.

A combination of Figures 13 and 14 shows the failure process of Sample B0. When the shear stress is close to the peak value (Point A, Figure 13), cracks initiated at the joint tips (Figure 14(a)). Once cracks initiated, cracks grew sharply, and completely coalesced along the shear plane (Figure 14(b)). Certainly, the peak shear stress was obtained from Point A to Point B, namely, during the crack coalescence. Comparing Figures 14(b) and 14(c), a good agreement between numerical and experimental failures can be found. For 0° cases, the shear plane was directly subjected to the shear stress, so the failure occurred along the shear direction on the shear plane. For the cases of A0 and B0, which are coplanar joints, the joint orientations are consistent with the shear plane. Once the cracks coalesced, A0 and B0 can be regarded as smooth joints. In this case, the upper block just slides along the ruptured interface. The normal displacement is also quite small. The reason is the same like the above-stated A150. If the shear load continued to apply, the friction on the shear plane contributed to a stable residual strength (Figure 14(d)). Since no SC was generated, the single-peak shear stress was observed. Though the joint arrangement in Series B is slightly different from the joint arrangement in Series A, the failure of A0 is like the failure of B0. The variation trend of Sample B0 is also like the variation trend above, but microcracks are basically observed on the shear plane.

Figure 13: Shear curves of Sample B0.

(a)

(b)

(c)

(d)

Figure 14: Fracturing process of Sample B0. (a) Point A. (b) Point B. (c) Real failure. (d) Point C.

For Sample B90, the joints are perpendicular to the shear plane as well as the shear direction. The shear plane traverses the midpoints of three joints. At a shear stress of approximately 0.9 τp (Figure 15), visual tensile cracks, namely, SCs, first initiated at most joint tips (Figure 16(a)). For one identical joint, two tensile cracks initiating at the upper and lower tips propagated towards opposite directions. As shown in Figure 16(b), even though MCs have not completely coalesced yet, the shear stress has reached the peak and then starts decreasing. During the shear stress decreasing stage, MCs completely coalesced obliquely (Figure 16(c)). Meanwhile, the shear stress was still dropping. Comparing Figures 16(c) and 16(d), the main and SCs are basically identical for numerical and experimental failures. Complete MCs separated the original sample into two blocks with a sawtooth interface (Figure 16(e)). The sawteeth mutually interlocked and prevented the upper block from climbing along the interface. With the continuous shear load, some sawteeth could be sheared off along the initial tensile cracks (Figure 16(f)), so that a second shear stress peak was created (Figure 15). The normal displacement trend with shear displacement for B90 is different from B0. The normal displacement increased with increasing shear displacement. It is because the upper block climbed slightly along the sawteeth interface. Due to cracking along the first initial SC, the second peak is lower than the first peak. The variation trend of the microcrack number is different from the aforementioned variation trends. The variation trend will be described in Section 3.7 in detail.

Figure 15: Shear curves of Sample B90.

(a)

(b)

(c)

(d)

(e)

(f)

Figure 16: Fracturing process of Sample B90. (a) Point A. (b) Point B. (c) Point C. (d) Real failure. (e) Multiple sawteeth. (f) Point D.

Figures 15 and 16 represent typical ruptures. Then, samples withJ different joint arrangements are known to experience different failure processes. The failure process is closely related to the joint arrangement. According to the whole fracturing process of all cases, three phases can be summarized as in the following sections.

3.4.1. First Phase: Cracking

The cracking phase includes the crack initiation, propagation, and coalescence as well as the SC initiation. All cases in Series A and Series B have to experience the cracking phase. Certainly, cracking types are distinguished. In Series A, the MC coalescence types are the same, but the SC initiation is influenced by the combination of shear direction and joint orientation. That is why the SC initiation can be observed for A30 and A60 but not for A90, A120, and A150. In Series B, MCs coalesced obliquely between two neighboring joints. SCs were generated in all cases except B0. However, the cracking mechanism of SCs for B30 and B60 is different from the cracking mechanism of B90, B120, and B150. Cracking is also influenced by the combination of shear direction and joint orientation.

3.4.2. Second Phase: Extrusion and Sawtooth Cutting

According to the abovementioned fracturing process, the second failure phase is obviously not necessary for all samples. As stated above, parts of samples such as A30 and A60, B30 and B60 experienced the extrusion at joint faces as shown in Figure 17(a) and Figure 17(b). A sawtooth interface was created due to MC coalescence. Sawteeth on the upper and lower blocks interlocked, under continuous shear loading. The extrusion effect took place on the contact face between the sawteeth. For B90, B120, and B150, after MC coalescences, some sawteeth were cut off along the first initial tensile cracks like that shown in Figure 17(c).

(a)

(b)

(c)

Figure 17: Schematic diagram of extrusion and sawtooth cutting.

3.4.3. Third Phase: Friction and Sliding

When the shear stress experienced a first and second peak, it came into the residual stage. On the interface, the friction and sliding effect played the main role. For A0 and B0, once cracks completely coalesced, blocks oppositely slid on the interface, so they did not experience the second phase and directly came into the third phase. For some cases, they came into the third phase after sawtooth cutting or extrusion. The sliding and friction effect on the interface contributes to the residual shear stress. For the smooth interface such as A0 and B0, the residual shear stress is low relative to that for the case in which the interface is composed of small fragments. Fragments result in undulate residual shear stress.

3.5. Contact Force Evolution around Joints

During the particle movement under the external loading, a great number of parallel bonds were subjected to compressive stress, but some parallel bonds may be subjected to tensile stress. Contact force including compression and tension is used to describe the mutual force action between particles. As shown in Section 3.4, different fracturing features can be observed at different shear stresses. Like the contact force illustrated in Figures 3 and 4, initial compressive contact forces are dispersively and uniformly distributed inside the whole sample while the tensile contact force mainly exists close to the sample edges before the normal stress was applied. As normal and shear stresses were applied, the external load forced particles to move, so the contact force distribution gradually changed. At the beginning of shear loading, compressive contact forces concentrated at joint tips (Figures 18(a1) and 18(b1)) while tensile contact forces concentrated on both sides of the joints (Figure 18(a2)). Then, compressive contact forces dramatically concentrated at rock bridges (Figure 18(c1)) so that cracks coalesced at rock bridges (Figure 10(c)). At this moment, tensile cracks first initiated at the middle joints due to tensile contact forces (Figure 18(b2)). During the MC coalescence, tensile contact forces also concentrated at rock bridges (Figure 18(c2)). The contact force concentration at rock bridges and joint tips indicates that the large internal stress is located there, which can well explain why MCs always coalesced along rock bridges. Throughout the contact force evolution, before MC coalescence, contact forces never took place on the middle joint face, and contact forces mainly changed near the shear plane. However, when rock bridges were through-going, the compressive contact force concentration occurred at the location of the middle joint (Figure 10(f)), which can account for the extrusion on the middle joint face (Figure 18(d1)). As shown in Figure 18(d2)), tensile contact forces concentrated at the middle joint, which resulted in a new SC (Figure 10(f)).

(a)

(b)

Figure 18: Contact force evolution for Sample A30 (shear direction: right to left): (a) compressive contact force; (b) tensile contact force.

Figure 19 shows the contact force evolution of Sample B90. At the beginning of loading, within the whole sample, compressive contact forces have not yet distinctly presented at joint tips, but tensile contact forces started intensively locating along joints (Figures 19(a1) and 19(a2)). Then, compressive contact forces concentrated at the joint tips. Meanwhile, tensile contact forces also presented at the joint tips. Cracks initiated at the joint tips, and tensile contact force concentration can account for the SC initiation (Figure 16(b)), which will be analyzed further using the principal stress around the middle joint in Section 3.6. Then, compressive contact forces mainly concentrated at rock bridges (Figures 19(c1) and 19(c2)) like the rock bridge shown in Figure 18(c1), which finally resulted in the rock bridge coalescence, that is, the MC coalescence. Tensile contact forces concentrated mainly at tensile crack tips (Figures 19(b2)–19(b4)), and then the tensile crack propagation followed the tensile contact force location. The tensile contact force concentration guided SCs to propagate.

(a)

(b)

Figure 19: Contact force evolution of B90 (shear direction: left to right): (a) compressive contact force; (b) tensile contact force.

The contact force evolution throughout the duration of shearing can be expected to explain the fracturing process. The location of the contact force concentration may be where cracks initiated, that is, where the local stress is high. The contact force concentration qualitatively represents the location of high internal stress and also indicates where cracks may initiate and propagate. Even though fracturing processes are various for samples with different joint orientations, the contact force evolution has to experience the following phases. The compressive contact force evolution should experience at least three stages, including the first uniform distribution, the concentration at joint tips and rock bridges, and the slight concentration on the failure plane. The tensile contact force evolution is partly different from the compressive contact force evolution, experiencing mainly a scattered force distribution close to the sample edges, concentration on joint sides and rock bridges, and concentration at the tensile crack tips. Therefore, the compressive contact force evolution is good to explain the MC propagation and coalescence, while the tensile contact force evolution is for the SC initiation and propagation.

3.6. Analysis of Principal Stresses around the Middle Joint

As shown in Figures 10, 12 and 16, cracks often initiated at the joint tips. To explain crack initiation at the joint tips and understand the stress evolution around joints, particularly the middle joint, the internal stress evolution of A30 and B90 was investigated as examples. Owing to the discreteness of the PFC model, stress, which is a continuum variable, actually does not exist anywhere in the particle assembly. To overcome this problem, an averaging procedure is adopted, which also bridges the gap between the microscale and a continuum. Therefore, the stress at one point is the average value in a measurement region of volume , computed as the following formula:

where is the number of contacts within the measurement region, and are the contact vector and the branch vector joining the centroids of the two contact ends, respectively, and represents the outer product. includes both normal and shear force components at contact but neglects the parallel-bond moment [34].

On the periphery of the middle joint, 16 measurement circles are arranged. The schematic diagrams of the measurement circle arrangement are shown in Figure 20. The diameter of each measurement circle is 3 mm. Each circle contains approximately 14 particles. Within one measurement circle, σxx, σyy, and τxy can be obtained, according to the following equation:

(a)

(b)

Figure 20: Schematic diagram of the measurement circle arrangement.

Then, the maximum principal stress () and minimum principal stress () for each measurement circle are computed. The positive and negative signs represent the tensile and compression stresses, respectively.

The principal stress distribution around the middle joint is illustrated using two typical samples, that is, A30 and B90. As seen in Figures 21(a) and 21(b), and are found to be varied at different locations in A30, and the highest stresses take place at circles 1 and 9, which are located at the joint tips. At circles 1 and 9, both the maximum and minimum principal stresses are always positive ones, which means that they are under the biaxial compression and tend to be shear failure. In other circles, the magnitude of computed principal stresses is much lower than the magnitude of computed principal stresses at circles 1 and 9. Furthermore, the magnitude of minimum principal stresses at circles 8 and 16 is the highest, exceeding the tensile strength of the synthetic sample. Combining the maximum principal stresses at these two locations (Figure 21(a)), circle 8 is under axial tension, and the stress state at circle 16 resembles that in the Brazilian tests. These well explain the appearance of tensile microcracks at circles 8 and 16 (Figure 10(a)). Circles 1 and 16 and circles 8 and 9 are located at the left-top and right-bottom of the middle joint, respectively. Because the shearing direction for the upper part is right to left, the location of circles 1 and 16 relative to the joint is along the shearing direction. This rule is also suitable for the location of circles 8 and 9. Taken together, all of these factors drive the crack initiation at the joint tips and guide the crack propagation direction.

(a)

(b)

Figure 21: Principal stresses around middle joint for A30 at (a) crack initiation stress and (b) peak shear stress.

Figure 22 shows the maximum and minimum principal stresses at different circle locations in B90 for the different shearing stages. Like sample A30, at circles 1 and 16 and 8 and 9 in sample B90, both the maximum and minimum principal stresses remain positive and their magnitude is much larger than the magnitude at other measurement circles. They are under biaxial compression and inclined to fail by shearing. Circles 1 and 16 and 8 and 9 are located at the right-top and left-bottom of the middle joint, respectively, which, relative to the joint, are along the shearing direction. It can be used to explain cracks formed at the joint tips and guide the crack propagation direction. In addition, the stress states at circles 2, 10, and 11 (Figure 22(a)) are analogous to the Mode I fracture, and the magnitude of minimum principal stresses at these locations exceeds the tensile strength of the synthetic sample. Thus, these two factors result in the formation of tensile microcracks at these locations (Figure 16(a)). Namely, it well accounts for the formation of SC-II cracks.

(a)

(b)

Figure 22: Principal stresses around middle joint for B90 at (a) crack initiation stress and (b) peak shear stress.

Overall, the stress concentration at the joint tips is associated with the crack initiation, and the location relative to the joint is along the shearing direction, which guides the crack propagation direction. Moreover, the formation of SC-II cracks is attributed to the tensile stress.

3.7. Variation Trend of the Number of Microcracks

In the numerical modeling process, the increase in the microcrack number is directly related to the cracking behavior. Cracking behaviors are various for different joint-arranged samples. The visual crack during the shear failure is an assembly of microcracks. Since the macrocrack propagation is associated with the shear displacement as shown in Figures 10, 12, 14, and 16, the increase in the microcrack number is also associated with the shear displacement as well as the fracturing behavior. Figure 23 shows the varying relationship between microcrack number and shear displacement. At different stages of shear displacement, growth rates for the microcrack number are different, which can be evaluated by the slope of the microcrack number versus the shear displacement curve. A big slope means a rapid growth rate. At the beginning of shearing, since the internal stress has not reached one value that can break parallel bonds, when the shear displacement increased from 0 to 1 mm, even though the shear stress dramatically increased with the shear displacement, the microcrack number obviously increased a little. This increase can be named a slow increase phase. This phase corresponds to the crack initiation. When the shear displacement ranged within 1 to 2 mm marked with two vertical red dashed lines, the microcrack number sharply increased approximately vertically, which corresponds to the MC propagation and coalescence like that shown in Figures 12 and 14. The approximately vertical phase is called the first rapid increase. For Series A, the first rapid increase in microcrack number was within 1 to 2 mm, but for Series B, the first rapid increase in microcrack number finished in less than 2 mm, such as the dashed curves of B30 and B60. The first rapid increase mainly corresponds to MC propagation and coalescence. In some cases of Series B such as the curves for B30 and B60, due to the SC propagation or extrusion, the second rapid increase could be observed but hardly took place in Series A. As the shear displacement continued after MC coalescence, due to the friction effect on the rupture interface, the microcrack number still continued to increase, but its growth rate became slow relative to the growth rate of the first rapid increase phase, which can be evaluated upon the curve trend. Generally, four phases of microcrack number increases, that is, slow increase, first rapid increase, second rapid increase, and continuous increase, are summarized. The second rapid increase can be observed in a few cases in Series B.

Figure 23: Microcrack number versus shear displacement.

4. Conclusions

In the present work, two series of intermittent joints with different joint orientations were used to carry out direct shear tests. Experimental and numerical approaches were used to explain the macro- and microfracture mechanisms. Some interesting conclusions are drawn as follows:

(1)Main and SCs are identified according to mechanism and crack pattern. They both initiate at the joint tips. MCs connect the tips of two adjacent joints. In Series A, MCs always propagate following the shortest trajectory. In Series B, they obliquely go through rock bridges. SCs may initiate prior to or after the peak shear stress and propagate in the direction that is approximately perpendicular to the joint orientation or parallel to the shear direction. SCs are all tensile cracks, while the MC is usually a mixed mode of tensile and shear cracks.(2)The fracturing process is significantly related to the joint orientation. The whole fracturing process is identified as three phases, which are characterized by cracking, extrusion or sawtooth cutting, and friction and sliding, respectively. But three phases are not undergone for all cases. The first and third phases are necessary for all cases, but the second phase is necessary only for A30, A60, B30, B60, B90, B120, and B150.(3)Shearing curves are characterized by the number of peaks, so two types, single-peak and double-peak curves, are identified. A single-peak curve is derived from A0, A90, A120, A150, and B0, while the double-peak curve is derived from A30, A60, B30, B60, B90, B120, and B150.(4)The variation trend of the crack initiation stress ratio is associated with the loading condition and the joint inclination. From the joint inclination of 0° to 90°, the crack initiation stress ratio under the shearing decreases first reaching the valley at 30° and then increases, which is different from the crack initiation stress ratio under the compression loading. From the joint inclination of 90° to 150°, the crack initiation stress ratio keeps increasing for Series A while decreasing for Series B. In essence, this difference originates from the variation of crack propagation and coalescense paths.(5)The compressive contact force evolution may undergo initial uniform distribution, concentration at the joint tips, concentration at rock bridges, and an extrusion surface. The tensile contact force evolution may experience initial scattered distribution close to the sample edges, concentration on both joint sides, with concentration at the rock bridges and SC tips. The compressive contact force evolution is good to explain MC propagation and coalescence, while the tensile contact force evolution is good for SC initiation and propagation.(6)Principal stresses around joints are associated with crack initiation and propagation. The highest maximum and minimum principal stresses take place at the joint tips and guide crack propagation, accounting well for crack initiation at the joint tips and propagation trajectory.(7)At the beginning of shearing, the number of microcracks increases slowly within a certain shear displacement mainly ranging from 0 to 1 mm. The increase in microcracks can be classified into four phases, namely, slow increase, first rapid increase, second rapid increase, and continuous increase. The four increase phases correspond to crack initiation, MC propagation, extrusion or SC propagation, and sliding on failure interface. The second rapid increase is observed only in the shearing of several samples that may undergo the second phase of shearing.

As analyzed above, the failure of discontinuous joints under a shear load mainly takes place on the shear plane. The failure characteristic is dependent significantly on the joint orientation and joint configuration. The presented results can be used in other rock mechanics applications, for instance, the sliding analysis of rock slope with discontinuous joints and the stability of an underground cavern with jointed rock mass.

Data Availability

The data used to support the findings of this study are available from the first author upon request.

Conflicts of Interest

The authors declare no conflicts of interest.

Acknowledgments

This study received its funding from project (2016M602742) supported by China Postdoctoral Science Foundation (Project nos. 51474249 and 51774322) supported by National Natural Science Foundation of China; Natural Science Basic Research Plan in Shaanxi Province of China (Program no. S2018-JC-QN-1500); Opening Fund of State Key Laboratory of Geohazard Prevention and Geoenvironment Protection (SKLGP2016K009); Open Fund of Key Laboratory of Bridge Engineering Safety Control by Hunan Province, Department of Education (Changsha University of Science & Technology) (15KB02); Research Fund of Department of Transportation in Zhejiang Province (2016019); and Traffic Construction Research Funds of Shanxi Province (2017-1-4). Kaihui Li would like to acknowledge the financial support from The Hong Kong Polytechnic University (RUF4). The authors wish to acknowledge these supports.

References

A. Bobet and H. H. Einstein, “Fracture coalescence in rock-type materials under uniaxial and biaxial compression,” International Journal of Rock Mechanics and Mining Sciences, vol. 35, no. 7, pp. 863–888, 1998.

E. Z. Lajtai, “Brittle fracture in compression,” International Journal of Fracture, vol. 10, no. 4, pp. 525–536, 1974.

H. S. Jang and B. A. Jang, “New method for shear strength determination of unfilled, unweathered rock joint,” Rock Mechanics and Rock Engineering, vol. 48, no. 4, pp. 1515–1534, 2015.

Q. S. Bai, S. H. Tu, and C. Zhang, “DEM investigation of the fracture mechanism of rock disc containing hole(s) and its influence on tensile strength,” Theoretical and Applied Fracture Mechanics, vol. 86, pp. 197–216, 2016.

X. Fan, P. H. S. W. Kulatilake, and X. Chen, “Mechanical behavior of rock-like jointed blocks with multi-non-persistent joints under uniaxial loading: a particle mechanics approach,” Engineering Geology, vol. 190, no. 5, pp. 17–32, 2015.

H. Lin, P. Cao, and Y. Wang, “Numerical simulation of a layered rock under triaxial compression,” International Journal of Rock Mechanics and Mining Sciences, vol. 60, no. 6, pp. 12–18, 2013.

X. Fan, P. H. S. W. Kulatilake, X. Chen, and P. Cao, “Crack initiation stress and strain of jointed rock containing multi-cracks under uniaxial compressive loading: a particle flow code approach,” Journal of Central South University, vol. 22, no. 2, pp. 638–645, 2015.

L. N. Y. Wong and H. Q. Li, “Numerical study on coalescence of two pre-existing coplanar flaws in rock,” International Journal of Solids and Structures, vol. 50, no. 22-23, pp. 3685–3706, 2013.

X. Chen, Z. Liao, and X. Peng, “Cracking process of rock mass models under uniaxial compression,” Journal of Central South University, vol. 20, no. 6, pp. 1661–1678, 2013.

M. Prudencio and M. V. S. Jan, “Strength and failure modes of rock mass models with non-persistent joints,” International Journal of Rock Mechanics and Mining Sciences, vol. 44, no. 6, pp. 890–902, 2007.

X. Fan, R. Chen, H. Lin, H. P. Lai, C. Y. Zhang, and Q. H. Zhao, “Cracking and failure in rock specimen containing combined flaw and hole under uniaxial compression,” Advances in Civil Engineering, vol. 2018, Article ID 9818250, 15 pages, 2018.

Y. L. Zhao, L. Y. Zhang, W. J. Wang, C. Z. Pu, W. Wan, and J. Z. Tang, “Cracking and stress–strain behavior of rock-like material containing two flaws under uniaxial compression,” Rock Mechanics and Rock Engineering, vol. 49, no. 7, pp. 2665–2687, 2016.

Y. L. Zhao, Y. X. Wang, W. J. Wang, W. Wan, and J. Z. Tang, “Modeling of non-linear rheological behavior of hard rock using triaxial rheological experiment,” International Journal of Rock Mechanics and Mining Sciences, vol. 93, pp. 66–75, 2017.

H. Lin, W. Xiong, and Q. Yan, “Three-dimensional effect of tensile strength in the standard Brazilian test considering contact length,” Geotechnical Testing Journal, vol. 39, no. 1, pp. 137–143, 2016.

S. Q. Yang, “Experimental study on deformation, peak strength and crack damage behavior of hollow sandstone under conventional triaxial compression,” Engineering Geology, vol. 213, pp. 11–24, 2016.

L. N. Y. Wong and X. P. Zhang, “Size effects on cracking behavior of flaw-containing specimens under compressive loading,” Rock Mechanics and Rock Engineering, vol. 47, no. 5, pp. 1921–1930, 2014.

X. P. Zhang and L. N. Y. Wong, “Cracking processes in rock-like material containing a single flaw under uniaxial compression: a numerical study based on parallel bonded-particle model approach,” Rock Mechanics and Rock Engineering, vol. 45, no. 5, pp. 711–737, 2012.

A. Bobet, “The initiation of secondary cracks in compression,” Engineering Fracture Mechanics, vol. 66, no. 2, pp. 187–219, 2000.

H. Lin, W. Xiong, and Q. Yan, “Modified formula for the tensile strength as obtained by the flattened Brazilian disk test,” Rock Mechanics and Rock Engineering, vol. 49, no. 4, pp. 1579–1586, 2016.

H. Lin, W. Xiong, W. Zhong, and C. Xia, “Location of the crack initiation points in the Brazilian disc test,” Geotechnical and Geological Engineering, vol. 32, no. 5, pp. 1339–1345, 2014.

H. Cheng, X. Zhou, J. Zhu, and Q. Qian, “The Effects of crack openings on crack initiation, propagation and coalescence behavior in rock-like materials under uniaxial compression,” Rock Mechanics and Rock Engineering, vol. 49, no. 9, pp. 3481–3494, 2016.

Y. H. Huang, S. Q. Yang, P. G. Ranjith, and J. Zhao, “Strength failure behavior and crack evolution mechanism of granite containing pre-existing non-coplanar holes: experimental study and particle flow modeling,” Computers and Geotechnics, vol. 88, pp. 182–198, 2017.

C. Gehle and H. K. Kutter, “Breakage and shear behaviour of intermittent rock joints,” International Journal of Rock Mechanics and Mining Sciences, vol. 40, no. 5, pp. 687–700, 2003.

E. Gerolymatou and T. Triantafyllidis, “Shearing of materials with intermittent joints,” Rock Mechanics and Rock Engineering, vol. 49, no. 7, pp. 2689–2700, 2016.

L. Jing, “A review of techniques, advances and outstanding issues in numerical modelling for rock mechanics and rock engineering,” International Journal of Rock Mechanics and Mining Sciences, vol. 40, no. 3, pp. 283–353, 2003.

E. Karampinos, J. Hadjigeorgiou, J. Hazzard, and P. Turcotte, “Discrete element modelling of the buckling phenomenon in deep hard rock mines,” International Journal of Rock Mechanics and Mining Sciences, vol. 80, pp. 346–356, 2015.

J. Jin, P. Cao, Y. Chen, C. Pu, D. Mao, and X. Fan, “Influence of single flaw on the failure process and energy mechanics of rock-like material,” Computers and Geotechnics, vol. 86, pp. 150–162, 2017.

D. O. Potyondy and P. A. Cundall, “A bonded-particle model for rock,” International Journal of Rock Mechanics and Mining Sciences, vol. 41, no. 8, pp. 1329–1364, 2004.

X. Fan, K. H. Li, H. P. Lai, Y. L. Xie, R. H. Cao, and J. Zheng, “Internal stress distribution and cracking around flaws and openings of rock block under uniaxial compression: a particle mechanics approach,” Computers and Geotechnics, vol. 102, pp. 28–38, 2018.

C. Cheng, X. Chen, and S. Zhang, “Multi-peak deformation behavior of jointed rock mass under uniaxial compression: insight from particle flow modeling,” Engineering Geology, vol. 213, pp. 25–45, 2016.

M. Bahaaddini, P. C. Hagan, R. Mitra, and B. K. Hebblewhite, “Parametric study of smooth joint parameters on the shear behaviour of rock joints,” Rock Mechanics and Rock Engineering, vol. 48, no. 3, pp. 923–940, 2015.

M. Bahaaddini, G. Sharrock, and B. K. Hebblewhite, “Numerical direct shear tests to model the shear behaviour of rock joints,” Computers and Geotechnics, vol. 51, pp. 101–115, 2013.

H. Q. Zhang, Z. Y. Zhao, C. A. Tang, and L. Song, “Numerical study of shear behavior of intermittent rock joints with different geometrical parameters,” International Journal of Rock Mechanics and Mining Sciences, vol. 43, no. 5, pp. 802–816, 2006.

Itasca Consulting Group, PFC3D User’s Manual, Itasca Consulting Group, Minneapolis, MN, USA, 2003.

X. Ding, L. Zhang, H. Zhu, and Q. Zhang, “Effect of model scale and particle size distribution on PFC3D simulation results,” Rock Mechanics and Rock Engineering, vol. 47, no. 6, pp. 2139–2156, 2014.