Rock Mech Rock Eng DOI 10.1007/s00603-017-1200-8

ORIGINAL PAPER

A Digital Image-Based Discrete Fracture Network Model and Its Numerical Investigation of Direct Shear Tests Peitao Wang1,2 • Meifeng Cai1,2 • Fenhua Ren1,2 • Changhong Li1,2 Tianhong Yang3

•

Received: 6 June 2016 / Accepted: 10 March 2017  Springer-Verlag Wien 2017

Abstract This paper develops a numerical approach to determine the mechanical behavior of discrete fractures network (DFN) models based on digital image processing technique and particle flow code (PFC2D). A series of direct shear tests of jointed rocks were numerically performed to study the effect of normal stress, friction coefficient and joint bond strength on the mechanical behavior of joint rock and evaluate the influence of micro-parameters on the shear properties of jointed rocks using the proposed approach. The complete shear stress–displacement curve of the DFN model under direct shear tests was presented to evaluate the failure processes of jointed rock. The results show that the peak and residual strength are sensitive to normal stress. A higher normal stress has a greater effect on the initiation and propagation of cracks. Additionally, an increase in the bond strength ratio results in an increase in the number of both shear and normal cracks. The friction coefficient was also found to have a significant influence on the shear strength and shear cracks. Increasing in the friction coefficient resulted in the decreasing in the initiation of normal cracks. The unique contribution of this paper is the proposed modeling

& Peitao Wang [email protected] 1

Key Laboratory of High-Efficient Mining and Safety of Metal Mines (Ministry of Education of China), University of Science and Technology Beijing, Beijing 100083, China

2

School of Civil and Resource Engineering, University of Science and Technology Beijing, Beijing 100083, China

3

Key Laboratory of Ministry of Education on Safe Mining of Deep Metal Mines, Northeastern University, Shenyang 110819, China

technique to simulate the mechanical behavior of jointed rock mass based on particle mechanics approaches. Keywords Jointed rock mass  Discrete fracture network  Digital image process  DFN modeling

1 Introduction Fractured rock masses are discontinuous and contain structures of varying scales, e.g., faults, bedding planes, joints or fractures. The existence and geomechanical properties of discontinuity like joints or fractures influence the mechanical behavior of a rock mass (Baghbanan and Jing 2007; Xu et al. 2013; Wasantha et al. 2015; Li et al. 2016; Johansson 2016). In a natural rock mass, the geometry of joints, including orientation, density, size and trace length, is very complex. The geometry of rock causes serious considerations in numerous research areas, including mining engineering, underground excavation and oil storage engineering. Thus, the characterization of complex joints is extremely important especially for mechanical research in numerical modeling (Jiang et al. 2006; Liu et al. 2009; Sun and Zhao 2010; Maleki 2011). Bahaaddini et al. (2013b) performed numerical simulations to study the shear behavior of rough joints using the Particle Flow Code in 2 Dimensions (PFC2D). They recorded the damage evolution during the loading process in terms of tensile and shear fracturing of rock asperities. Ivars et al. (2011) described a synthetic rock mass modeling technique investigating the mechanical behavior of fractured rock masses. The in situ joint network was represented by the discrete fractures network (DFN) model (Dershowitz and Einstein 1988), and the corresponding mechanical properties were studied using PFC3D code. According to the

123

P. Wang et al.

research, it was pointed out that the DFN simulation could be considered the most logical choice to explicitly represent the structure of the in situ joint fabric. PFC, developed by Itasca (2004), has been shown to have advantages in simulations of the mechanical behavior of rocks (Hunt et al. 2003; Tran et al. 2009; Lambert et al. 2010; Zhang et al. 2011) and discontinuities in fractured rock masses (Wang et al. 2003; Cho et al. 2012; Bahaaddini et al. 2013a, 2014, 2016; Sarfarazi et al. 2014; Park and Min 2015; Potyondy 2015). There are two major modeling methods for generating rock assemblies, the expanding method and the exploding method (Itasca 2004). Micro-parameters, such as particle radius, particle number and specimen porosity, were the key influencers on which these models were established. Basic rock specimens could be modeled using these two intrinsic methods, while complex fractured rock masses cannot be modeled directly based on the existing methods. Some researchers proposed new methods for modeling rock materials. Jiang et al. (2003), for example, introduced an efficient technique called the multilayer with undercompaction method to generate homogeneous specimens using PFC2D. This method could provide an efficient way to generate rock specimens of differing densities, ranging from loose to dense states. Wang et al. (2013a, b) provided a new method to characterize jointed rock masses. The joint network was generated based on the results of an in situ survey, and the corresponding mechanical properties of contacts between joint elements were identified. The corresponding uniaxial compression strength and failure pattern were also examined in their study. Wang et al. (2016) established a DFN model based on the Monte Carlo method using a contactbonded model. They discussed the scale effects of elasticity and failure patterns based on DFN models using varied sample sizes. Vallejos et al. (2016) introduced a methodology for generating a deterministic DFN model in PFC that was based on the veins mapped at the surface of each core sample and discussed the impact of stiffness and strength of vein networks on the behavior of rock masses. The goal of modeling using PFC2D code is to make the numerical model of jointed rock mass reflect the actual geometrical distribution accurately and efficiently. The digital image processing (DIP) technique is a practical tool for determining the geometrical distribution of joints based on the measurements of digital images of the rock surface (Chen et al. 2016). These data can be efficiently processed by automatic algorithms. The modeling of jointed rock masses is robust because the practical distribution of structures can be imported directly into the numerical models. The differences in the grades of material colors could be used to distinguish material types. Many researchers have used digital imaging to characterize rock structures, map the structures into various numerical codes

123

and examine the mechanical behavior of rocks by explicitly taking into account the geometrical fabric of jointed rocks. Reid and Harrison (2000) proposed a methodology for semi-automated detection of discontinuity traces in grayscale digital images of a rock mass. They demonstrated the efficacy of this image processing method and produced a discontinuity trace map to support their methodology. Zhu et al. (2006) applied DIP to characterize the actual spatial distribution of different materials and components in the rock sample using rock failure process analysis (RFPA) code and clarified the efficiency of this approach in mechanical studies of fractured rocks. Xu et al. (2008) developed a DIP technique based on the finite element method to study soil–rock mixture (SRM). In their study, they built a mesostructured concept model of SRM that represented the heterogeneity of SRM and analyzed the mechanical properties. Ghabraie et al. (2015) presented an efficient DIP technique for monitoring strata movement during the mining process. With the progress of measuring technologies, computerized tomography (CT) becomes a commonly used method for studying the internal structures or fractures of rocks (Yun et al. 2013; Cnudde and Boone 2013; Cai et al. 2014). Moreover, CT enables the visualization of internal microstructures. Therefore, it is frequently applied in combination with DIP techniques to characterize rock structures (Dai 2011; Huang et al. 2015; Yu et al. 2016). Images obtained from in situ surveys or CT measuring improves the efficiency and accuracy of rock models. The combination of digital images and PFC modeling can be an efficient way to characterize rocks. Until recently, however, only a few studies have been reported in the field of concrete modeling. For jointed rocks, Asadi et al. (2012) proposed a modeling method to create a roughness joint in the PFC2D model. They used a defined function y = f(x) or coordinate points (x, y) to represent a chosen 2D joint profile with the method. They also showed the effects of profile roughness on shear behaviors like peak shear strength and shear failure patterns. Generally speaking, the distribution and geometry of joints in natural rock mass are random and complex. Therefore, a particularly defined function cannot efficiently characterize the DFN. In this paper, an image processing technique that uses MATLAB code to capture the joint elements of a fractured rock mass based on in situ surveys was employed. The geometrical fabric of the joints and matrix was identified from the image of discrete fractures (IDF) and integrated into a well-established code (PFC2D). Then, a numerical model was built that can take into account the actual shape and geometrical distribution of materials within the mass. The direct shear tests were carried out on jointed rocks and mechanical behaviors such as failure process, micro-crack

A Digital Image-Based Discrete Fracture Network Model and Its Numerical Investigation of…

initiation and propagation, and contact force distribution were simulated and analyzed.

2 The Proposed Approach 2.1 Capture of Image of Discrete Fractures The establishment of a 3D DFN geological model based on the Monte Carlo method was reported by Wang et al. (2013a, b) and Yang et al. (2015) However, a statistical distribution of the geometry of the discrete fractures cannot be used directly in PFC modeling. Generally, IDF from in situ surveys or geological models can provide useful information about fractures. In this section, the IDF was obtained from the 3D DFN model proposed by Yang et al. (2015), as shown in Fig. 1. An arbitrary profile from the 3D DFN model was selected, and the corresponding 2D image of the DFN model was obtained. The information from the 2D IDF can then be processed using MATLAB code and the PFC model based on IDF established. 2.2 Principles of the Technique As shown in Fig. 2, the intensity of each pixel in the digital image was used to identify joint and intact rock. Pixels in darker gray represent joint elements, and pixels in white represent rock materials. The edges between the joints and matrix are obvious, as is clear in the image. Using the IDF shown in Fig. 2 as an example, the process is described as follows. Firstly, the image was treated using MATLAB code, from which a matrix named Gray [Col, Row] would be established. The size of the image shown in Fig. 2 is 351 9 351 pixels, and the value of the gray gradient is from 0 to 255. The matrix with 351 columns and 351 rows will be given, and the coordinates of particles in PFC will be determined by the column and row. If a threshold value of 250 is used, then all pixels having a relative intensity lower than 250 were identified, as shown in Fig. 2. For example,

the gray value of a pixel in column 10, row 10 is 0, then the ball of x = 10 and y = 10 will be identified as a joint element. Different gray values represent different types of materials. Then, the coordinates of intact rocks and joint elements will be automatically exported as a TXT file which PFC code can easily call. There were a total of 123,201 elements in the model, with 111,058 rock elements and 12,143 joint elements, respectively. It should be noticed that the threshold was used as a filter so that only those pixels that met the threshold were retained. For each image, a threshold specific to each detector should be established in order to obtain maximum recognition of the joint trace. In order to establish a proper model in PFC2D, the particle radius should be calibrated first. The size of the DFN model shown in Fig. 2 is 10 m 9 10 m, and the corresponding resolution size is 351 9 351 pixels. The radius of particles could be defined as shown in Eq. (1): n R0 ¼ ; ð1Þ 2N where R0 is the particle radius, n is the actual size of the DFN mode and N is the number of column of the IDF. The row and column of the pixel matrix Gray [Col, Row] can be treated as the y-coordinate and x-coordinate, respectively. The corresponding PFC model established by the proposed method is shown in Fig. 3. Since a uniform particle distribution (Rmax/Rmin = 1) may cause a formation of distinct blocks of particles with a specific orientation or structure, which may influence the fracture patterns, rock should generally be presented as a dense packing of non-uniform-sized circular particles (Potyondy and Cundall 2004). The proposed approach can then generate specimens using either uniform distribution or Weibull distribution. Equation (2) shows the relation of uniform distribution. R ¼ ravg þ rvar  rand;

ð2Þ

where R is the particle radius; ravg is the average value of particle radius; rvar is the variation of particle radius; and

σ

O Investigation of in-situ structures

Establish of 3D discrete An arbitrary section of 3D DFN fractures network (DFN) model model

ε Mechanical analysis of 2D DFN model

Fig. 1 The sketch map of discrete fractures network

123

P. Wang et al.

251 0 251 251

251 251 0 251 251 251 251 251 0

251

0 251 251 251

251 251 251

0

251 251

0 251 251 251

0

251 251

251 251 251 0 251 251 251

0 251

251 251

0

251

0 251 251 251

251 0

0

251 251

0 251

0 251 251

0 251

Joint element Rock element

Fig. 2 The principle of processing digital images (rock face area of 10 m 9 10 m) Fig. 3 The DFN model based on the image processing method

rand is the uniformly distributed random numbers within the region of [0,1]. Using Weibull distribution (Weibull 1951) with the proposed technique is shown in Eqs. (3) and (4): m r m1 r /ðRÞ ¼ ð Þ exp½ð Þm  ð3Þ ravg ravg ravg 1

R ¼ ravg ½ln ð1  randÞm ;

ð4Þ

where R is the element radius; ravg is the mean value of particle radius; m is the homogeneity index of the rock; and rand is the uniformly distributed random numbers within the region of [0, 1]. According to the definition, the higher value of m, the more heterogeneous radius distribution in the material, and vice versa (Tang et al. 2000; Feng et al. 2006; Zhu et al. 2006). Fig. 4 Image of discrete fractures with a scale size of 5 m 9 5 m

123

A Digital Image-Based Discrete Fracture Network Model and Its Numerical Investigation of… Fig. 5 Three types of rock assemblies with different particle radius distributions and the corresponding contact bonds distribution

Figure 4 shows an IDF cut from the image in Fig. 2 with scale size of 5 m 9 5 m. The corresponding specimens obey different particle size distributions, and the corresponding enlarged images are shown in Fig. 5.

3 PFC Modeling and Numerical Simulations 3.1 The Models from IDF The DFN model based on the 5 m 9 5 m IDF in Fig. 4 is shown in Fig. 6. In this model, there are a total of 28,900 particles, with 3451 joint elements and 25,449 rock elements, respectively. Among the particles assembled in the DFN model, particles in blue are intact rock and those in red are joint elements. The particle sizes in the proposed

modeling technique depend mainly on the resolution of digital DFN images as well as the scale size of the survey area. Each pixel in the IDF represents an individual particle. The more the image pixels are, the higher the number of particles will become. It is sure that the smaller the particle size becomes, the higher computational accuracy will be obtained. Although assembles with much smaller particles in millimeter scale or smaller can be established using the proposed modeling method, the PFC2D method cannot simulate the rock model because of the restriction of current computing ability of computers. Thus, it should be pointed out that this deficiency may create errors in computing joints shear strength. In addition, factors that affect the rock mass behavior at the grain scale (e.g., grain size, porosity), particle interlocking issues and micro-roughness created by large size particles, are not addressed explicitly

123

P. Wang et al. 2

Fig. 6 Initiated and interconnected DFN model 3

40

30

4

1

in the present work. Thus, limitations exit when dealing with the mechanical behaviors of jointed rock mass because the unexpected issues may bring some significant unrealistic effects. With the increase in computing power, it will be feasible to carry out simulations with assemblies of much smaller particles to represent the realistic grain size distribution seen in rock minerals. 3.2 Numerical Direct Shear Test on DFN Models 3.2.1 Model Preparation In order to avoid packing, the particle radii were chosen to be uniformly distributed in the model. Generally, the mechanical parameters of joint elements are relatively weaker than those of intact rock elements. The weaker the joint element is, the more the joint will control the mechanical characteristics of a jointed rock mass. The quantitative ratios of joint elements to rock elements, however, have not been reported. Thus, the stiffness and bond strength of joint elements are generally determined by trial and error. In most cases, the ratios of mechanical properties like elastic modulus and strength between joints Table 1 Micro-parameters of the contact bond model and the joint elements model

123

and rocks have been set at 1–20% from previous studies (Jia and Tang 2008; Li et al. 2011; Xu et al. 2013; Yang et al. 2015). Asadi et al. (2012) discussed different failure modes during shearing of the models with a single fracture using PFC2D method. In those studies, zero contact bond strength was assigned to the joint particles, and different failure modes including asperity sliding, cutoff and asperity degradation were explicitly observed. The understanding of rock fracture shear behavior was also expanded through observation under microscopic conditions. Kim et al. (2016) investigated the micro-mechanical behaviors of transversely isotropic rocks and provided a set of microparameters for two types of rocks. The mechanical parameters of intact rock and joint elements were chosen according to the models by Kim et al. (2016), Asadi et al. (2012) and Wang et al. (2016). According to the reports of Wang et al. (2016), a reasonable comparison of elasticity, strength and failure patterns of jointed rocks was obtained. Thus, the ratio of the bond strength of joint elements to rock elements using contact-bonded model was set to be 5% according to Wang et al. (2016). The stiffness, bond strength and friction coefficients of rock and joint elements are listed in Table 1.

Micro-parameters

Intact rock elements

Joint elements

Minimum particle radius, Rmin (m)

0.147

0.147

Particle size ratio, Rmax/Rmin Number of particles

1.66 25,449

1.66 3451

Contact elastic modulus (GPa)

75

60

Contact stiffness ratio, kn/ks

1.5

1.5

Friction coefficient

0.5

0.3

Bond normal strength (MPa)

5.0

0.25

Bond shear strength (MPa)

5.0

0.25

A Digital Image-Based Discrete Fracture Network Model and Its Numerical Investigation of…

3.2.2 Shear Loading Set Up The failure behavior of jointed rock under shear loading has attracted much attention among geotechnical researchers. Sarfarazi et al. (2014) performed numeric shear testing on the rocks with various types of overlap joints and emphasized that the failure patterns were highly dependent on the type of joint overlap. Thus, we conducted a shear test on the proposed DFN model in this section. In a numerical test, the normal stress applied to the sample is 0.10 MPa. Shear loadings are then applied to the sample by moving the left and right wall in the positive and negative x directions, respectively, with a sufficiently low velocity (i.e., 0.020 m/s) to ensure quasi-static equilibrium (Ghazvinian et al. 2011). Shear displacement is measured by tracing the horizontal displacements of walls 3 and 4, as shown in Fig. 6. Shear force is recorded by averaging the reaction forces on walls 3 and 4. The normal stress is applied to wall 2 through accurate servo controlling. 3.3 Results of Shear Loading Tests 3.3.1 Fracture Pattern and Shear Loading Behavior

Shear stress / MPa

A reliable DFN model using PFC method should require a significantly large number of particles to simulate largescale rock mass models. Considering the restriction of current computing ability of computers, many researchers have discussed the effect of the ratio of model scale to particle size distribution. Koyama and Jing (2007) have pointed it out that the diameter of the specimen should be at least 20 times of the largest mineral grain size in the rock. In this condition, the uncertainties of using the model will be minimized to an acceptable level. Ding et al. (2013) also investigated the effect of model scale (L) and particle 60

3500

50

3000

size (d) distribution on the simulated macroscopic mechanical properties. They recommend a reasonable ratio of 25 (L/d) for the simulation of rock specimens. The ratio of sample size to mineral grain size in our manuscript is about 34 and should be reasonable to model the jointed rock mass according to Koyama and Jing (2007) and Ding et al. (2013). Figure 7 shows the fracture pattern of the DFN model and the corresponding shear stress–displacement relationship and accumulative micro-cracks. As shown in Fig. 7a, there is a strong relationship between the sharp increase in the number of micro-cracks and the increase in shear stress. The accumulative numbers of both shear and tensile cracks increased dramatically when the shear displacement was lower than 0.8 m, at which point the peak shear stress had been reached. The rates of curves of cracks against shear displacement all decrease after peak shear strength. Moreover, more shear cracks were observed than normal cracks. Unlike the shear behavior of jointed rocks with smooth planes (Zhang et al. 2006), the shear stress–displacement curve of the DFN model showed a fluctuating drop after a peak shear strength of 56.57 MPa. The reason for this stress drop is the jump phenomenon of shear fracture propagation in the DFN model. After a fracture propagates, a stress relaxation occurs and fracture propagation stops. As shear loading increases, the concentration of shear stress in the fracture tip increases and the fracture propagation starts again. The process continues until the shear stress reaches a residual shear stress related to the friction resulting from the roughness of the shearing surface, as shown in Fig. 7b. 3.3.2 Contact Force Chains Figure 8 shows the distribution of contact bond force in the models during the shear loading process. The blue and red

2500

40

2000 30

Roughness shear surface

1500 20 Stress-strain curve Shear cracks

10

Normal cracks

1000 500

Joint element Rock matrix

0 0 0.0 0.2 0.4 0.6 0.8 1.0 1.2 1.4 1.6 1.8 2.0

Shear displacement / m

(a)

(b)

Fig. 7 Mechanical behaviors of the DFN model under direct shear test

123

P. Wang et al. Fig. 8 Evolution of the contact force chains in the models during shear loading

lines represent the compressive and tensile forces in the model, respectively. The width of the force chain represents the magnitude of the contact force. As shown in Fig. 8, the contact forces in the specimens were relatively evenly distributed in the early stages of testing. With the increasing in shear loading, the contact force chains were concentrated when the sample approached the peak strength, as shown in Fig. 8d. Figure 9 shows the distribution of contact forces when the shear displacement was 1.5 m. According to this image, the compressive force was mainly distributed along chains between particles. However, the tensile force between particles was uniformly distributed in a single unit of the specimen.

However, when the stress reached peak shear strength, breakages of joints bonds were observed. It is also shown that the shear band along the shear direction is very complex. The main reason should be the influence of particles interlocking when across several crossing joint sets. It is shown that shear cracks occurred mainly near the left and right loading walls, and no shear crack coalescence was observed at this stage. When the shear displacement reached 0.5 m, coalescence of macro-fractures occurred along the shearing direction. When the shear displacement reached 1.5 m, the contact bonds along the shearing direction almost broke. The residual strength of the DFN model was controlled by the coefficient friction of particles along the roughness surface.

3.3.3 Contact Bonds Breakage 3.3.4 Particle Velocity Vectors Figure 10 shows the simulated breakage of contact bonds of particles at different loading stages. The black lines represent the contact bond between particles. In the early stages, the contact bonds were normally distributed. Fig. 9 Numerical results of the contact forces distribution. Blue lines indicate compression force, and red lines indicate tension force (color figure online)

123

When the shearing began, the porosity of the granular materials changed. The granules redistribution happened and formed a vortex area. Figure 11 depicts the velocity

A Digital Image-Based Discrete Fracture Network Model and Its Numerical Investigation of…

(a)

(b)

(c)

(d)

Fig. 10 Development of contact bonds distribution during the shear test. a Initial model. b The model at peak shear strength. c The model at shear displacement = 0.5 m. d The model at shear displacement = 1.5 m

2 3

40

30

4

3

3

40

30

2

2

2

4

3

40

30

4

40

30

4

1

1

1

1

(a)

(b)

(c)

(d)

Fig. 11 Distribution of velocity vectors in the samples during shear loading

vectors, which were consistently distributed near the top and bottom wall in the early stage. However, the velocity vectors along the shearing direction became obviously nonuniform during the shear loading test, and, as also shown in Fig. 11, a jump phenomenon of shear fracture propagation was observed during the shear process. Furthermore, a vortex phenomenon in the velocity vectors appeared along the shearing direction during the test, as shown in Fig. 12. The particles moved in the opposite direction during shear loading from the concentration of shear stress in the fracture tip. After the fractures near this zone propagated, stress relaxation occurred and

fracture propagation stopped. Then, a new concentration of shear stress in another fracture tip occurred, and a jump phenomenon of shear fracture propagation started again. Along with the development of the shear progress, the radius of shear vertex areas gradually shrank or even disappeared, eventually merged to form a continuous curve strip zone throughout the whole surface of the particle materials, which could be called a shear band. From this simulation, we gain a better understanding of the mechanics of material shear strength. This phenomenon may be one of the principal reasons for the fluctuating drop of the shear stress–displacement curve shown in Fig. 7a.

Fig. 12 The vortex phenomenon of velocity vectors during shear loading

123

P. Wang et al.

3.3.5 Cracks Initiation and Propagation

2

Particles are bonded together at contact points in the PFC model. The bond breaks when the shear force exceeds the bond strength. Each bond breakage is considered as a micro-crack in a rock sample modeled in PFC2D. Two types of cracks—shear cracks and tensile cracks—could be simulated in the PFC2D model. A shear crack occurs when a bond’s shear strength is exceeded. A tensile crack occurs when a bond’s normal strength is exceeded. Micro-cracks were monitored throughout the shear loading experiment. The locations of micro-cracks that occurred under different shear loading stages are shown in Fig. 13. Each circle represents a micro-crack. The circles in red and black denote normal cracks and shear cracks, respectively. The diameter of the circle represents the relative magnitude of the released energy. Note that micro-cracks were distributed throughout the specimen prior to macro-fracture nucleation, and it was difficult to predict where the cracks would initiate at the early stage (Fig. 13a–c). As testing proceeded, a few micro-cracks still were initiated throughout the specimen, while more events were clustered along the shear direction, as shown in Fig. 13d–f. A potential nucleation site in the failure process had formed. The micro-fractures began to be clustered at the stages shown in Fig. 13e and f and became clearly localized in the stages shown in Fig. 13g and h. Moreover, shear cracks appeared during the shear loading process along the joint directions. Finally, the interior shear micro-cracks formed a rotated Z-shape, as shown in Fig. 13. It is important to note that although the macro-fractures were shear failure modes, most of them were distributed in a highly stressed compressive zone. As for physical experimental studies, crack initiation and stress distribution could not be directly observed until recently. A better understanding of the mechanical

3

40

30

4

1

Fig. 14 The contact force distribution and cracks in the specimen after shear failure

behaviors under shear test can be obtained by contrasting stress and cracks in the specimen. The micro-cracks and contact force chains at a shear displacement of 1.50 m are shown in Fig. 14. The normal cracks, which are indicated by red, are distributed mainly in accordance with the trend of compressive forces. Moreover, the normal cracks occur mainly at contact bonds of joint elements. The shear cracks, indicated in black, are distributed mainly along the shearing direction, and more shear cracks are observed than normal cracks according to Fig. 14.

4 Further Discussions The mechanical properties of the DFN model under direct shear loading are presented. Generally, the shear behavior of a fractured rock mass depends on the normal stress acting perpendicular to the loading direction and the

(a)

(b)

(c)

(d)

(e)

(f)

(g)

(h)

(i)

(j)

(k)

(l)

Normal cracks

Shear cracks

Fig. 13 Cracks distribution during the shear test. Red and black balls represent tensile cracks and shear cracks, respectively (color figure online)

123

A Digital Image-Based Discrete Fracture Network Model and Its Numerical Investigation of…

stresses increased, and the rate of changes followed a similar pattern. It should be noted that these parameters showed more sensitive responses when the normal stress was lower. As shown in Fig. 16a, the peak shear strength increased dramatically when the normal stress was lower than 1.0 MPa. Shear displacements at peak strength decreased with increasing in normal stress. The shear modulus was higher after an increase in normal stress, as shown in Fig. 15a, indicating that the DFN model reached peak strength at a lower shear displacement. Figure 16b shows the number of shear cracks, normal cracks and total cracks occurred during the shear test. Approximately increasing and decreasing trends were observed for shear and normal micro-cracks, respectively, under various normal stresses. Furthermore, a higher number of normal micro-cracks than shear cracks occurred at the same normal stress. Therefore, it can be concluded that the higher normal stress, the greater impact on the initiation and propagation of cracks in DFN models.

mechanical properties of the rock materials. This section discusses the shear behaviors of the proposed DFN model under various conditions. 4.1 Effect of Normal Stresses The normal stresses applied to the sample were 0.05, 0.1, 0.5, 1.0, 2.0 and 5.0 MPa. Figure 15 shows the shear stress and the total number of micro-cracks occurred with shear displacement at the corresponding normal stresses. The results of shear strength agree with the general observations from laboratory tests for rock materials. Figure 16a shows the changes in the peak strength, residual strength and peak displacements under normal stresses. For comparison, the residual strength in this study was set as the stress at which the shear displacement was 1.0 m. The peak shear strength and residual shear strength varied from 48.8 to 70.3 MPa and from 19.1 to 44.3 MPa, respectively. Both of the factors increased as normal

9000

70 8000

Number of micro-cracks

Shear stress / MPa

60 50 40 30 Normal stress=0.05 MPa Normal stress=0.1 MPa Normal stress=0.5 MPa Normal stress=1.0 MPa Normal stress=2.0 MPa Normal stress=5.0 MPa

20 10

7000 6000 5000 4000 Normal stress=0.05 MPa

3000

0.2

0.4

Normal stress=0.5 MPa

2000

Normal stress=1.0 MPa Normal stress=2.0 MPa

1000

Normal stress=5.0 MPa

0 0.0

Normal stress=0.1 MPa

0.6

0.8

1.0

0

1.2

0.0

0.1

0.2

Shear displacement / m

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1.0

1.1

1.2

Shear displacement / m

(a)

(b)

Fig. 15 Influence of normal stress on the direct shear strength and micro-cracks of the DFN models 80

8000

0.25

7000

60 50

0.15

40 0.1

30 Peak strength

20

0

1

2

3

Normal stress / MPa

(a)

4

5000 4000 3000 2000 Shear cracks Normal cracks Total cracks

1000

Peak displacement

0

6000

0.05

Residual strength

10

Displacement / mm

0.2

Number of cracks

70

Shear strength / MPa

Fig. 16 Peak strength, residual strength and peak shear displacement versus normal stress

5

0

0

0

1

2

3

4

5

Normal stress / MPa

(b)

123

P. Wang et al.

4.2 Effect of Friction Coefficient Park and Song (2009) carried out an extensive series of direct shear tests using PFC3D. They examined the shear behavior of a rock sample with a rough joint, with specific consideration of different joint roughness coefficients (JRC). The friction coefficient has been found to be the most important factor governing the shear strength and dilation angle. In this section, the effect of the friction coefficient is examined by carrying out direct shear tests on DFN models with different friction coefficients ranging from 0.0 to 1.5 (f = 0.0, 0.1, 0.3, 0.5, 0.8, 1.0,1.2, 1.5). Figures 17 show the shearing stress and micro-cracks at the normal stress of 0.1 MPa but with various friction coefficients. The shear modulus increases as the friction coefficient increases, as shown in Fig. 17a. As shown in

70

Friction coefficient=0.0 Friction coefficient=0.3 Friction coefficient=0.8

Friction coefficient=0.1 Friction coefficient=0.5 Friction coefficient=1.0

Friction coefficient=1.2

Friction coefficient=1.5

Fig. 17b, the number of the micro-cracks indicates sensitivity in the relationship between cracking and the friction coefficient. Figure 18a shows the peak and residual shear strength as related to the friction coefficient, indicating that both the peak and residual strength increased with the friction coefficient. The values vary from 35.74 to 60.89 MPa and from 7.20 to 37.49 MPa for peak and residual strength, respectively. When the friction coefficient is lower than 0.5, the peak and residual shear strength would be influenced more by a significant increase in the friction coefficient. However, the increase in the friction coefficient has no measurable influence on the displacement at the peak shear strength. Additionally, the number of shear cracks increases linearly with the increasing in the friction coefficient, as shown in Fig. 18b. Moreover, the increasing in friction coefficient decreases the initiation of normal cracks.

9000

Friction coefficient=0.1 Friction coefficient=0.5 Friction coefficient=1.0 Friction coefficient=1.5

8000

Number of micro cracks

60 50

Shear stress / MPa

Friction coefficient=0.0 Friction coefficient=0.3 Friction coefficient=0.8 Friction coefficient=1.2

40 30 20

7000 6000 5000 4000 3000 2000

10

1000 0

0 0.0

0.2

0.4

0.6

0.8

1.0

0.0

1.2

0.1

0.2

0.3

0.4

0.5

0.6

0.7

Shear displacement / m

Shear displacement / m

(a)

(b)

0.8

0.9

1.0

1.1

0.35

60

0.3

50

0.25

40

0.2

30

0.15

20

Peak strength Residual strength

10 0

123

Peak displacement

0.1 0.05 0

0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 1.1 1.2 1.3 1.4 1.5

9000 8000

Number of cracks

70

Displacement / mm

Fig. 18 Influence of the material friction coefficients on the shear stress–strain relationship

Shear strength / MPa

Fig. 17 Influence of the material friction coefficient on the shear stress–displacement relationship

7000 6000 5000 4000 3000 2000

Shear cracks Normal cracks Total cracks

1000 0

0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 1.1 1.2 1.3 1.4 1.5

Friction coefficient

Friction coefficient

(a)

(b)

0.26

8000

60

0.24

7000

0.22

55

0.2

50

0.18 45 0.16 40 Peak strength

35

Residual strength

0.14 0.12

Number of cracks

65

Displacement / mm

Shear strength / MPa

A Digital Image-Based Discrete Fracture Network Model and Its Numerical Investigation of…

6000 5000 4000 3000 2000 Shear cracks Normal cracks Total cracks

1000

Peak displacement

30

0

20 40 60 80 100 120 140 160 180 200 220 240 260 280 300

0.1

0

Ratio of bond strength of joint to rock (%)

0

20

40

60

80 100 120 140 160 180 200 220 240 260 280 300

Ratio of bond strength of joint to rock (%)

(a)

(b)

Fig. 19 Shear behaviors and cracks versus ratios of bond strength of joint to rock

4.3 Effect of Joint Bond Strength The strength of the joint bond is the key factor on the mechanical properties of DFN models. Wang et al. (2016) have conducted a comparison of experimental and numerical tests on stratified rocks to verify the feasibility of the joint model. They also discussed the influence of the ratio of joint bond strength to intact rock on the uniaxial compressive strength and proposed proper intervals of ratios for testing joint bond strength. They defined K and S as the ratio of bond strength between the joint model and the intact rock model. The value of K and S are defined in Eq. (5) and (6) as follows: K¼

n bond 0  100% n bond

ð5Þ

S¼

s bond 0  100%; s bond

ð6Þ

where n_bond0 is the normal bond strength of the joint model; n_bond is the normal bond strength of the intact rock; s_bond0 is the shear bond strength of the joint model; and s_bond is the shear bond strength of the intact rock. Generally, the strength of joint elements is lower than that of intact rock materials. However, the mechanical parameters of the joint layer could be higher than those for intact rocks, as is the case in some kinds of sandstone (Kim et al. 2016). In this section, fifteen numerical tests were performed to investigate the influence of joint bond strength on shear behavior in DFN models. The values of K (S) were set as 0, 1, 3, 5, 10, 20, 30, 50, 80, 100, 110, 120, 150 and 200%. Figure 19a shows the ratio of peak and residual shear strength against the bond strength at a normal stress of 0.1 MPa. Both the peak shear strength and residual strength increase as the bond strength ratio increases. The bond strength ratio has no measurable effect on the

displacement at peak shear strength. As noted in Fig. 19b, the increase in the bond strength ratio resulted in an increase in the number of both shear and normal cracks. The contact bond breakages in DFN models with varied strength ratios of the joint element to intact rock are shown in Fig. 20. The macro-fracture patterns of models with K(S) of 0, 5 and 30% are relatively similar, and the bonds of joint elements break in these cases. When K(S) is 100%, the DFN model can be treated as an intact rock, with macro-fractures occurring mainly along the shearing direction. In the scenarios where K(S) is 150 and 300%, the bond strength of the joint elements is relatively higher than that of intact rocks. The joints in these models can be treated as cables, as shown in Li et al. (2015), and the peak shear strength of the DFN models is higher, as shown in Fig. 19a. It should be noted that the present work is a 2D analysis and this model is more suitable for plane problems, such as discontinuities with a large spatial distribution in tunnels (as shown in Fig. 21a, b). The 3D problem can be simplified as a plane problem as shown in Fig. 21c. Thus, PFC2D could be a logical choice for simulating the behavior of a 2D DFN model. It should be emphasized that boundary effects may exist if the spatial distribution of joints is not large enough.

5 Conclusions A new technique of establishing the DFN model based on digital IDF and PFC2D code was introduced. The images of fractures were obtained from 3D DFN models based on in situ observations. The failure progress of the DFN models was visually represented by the contact bond chains, contact force distribution and micro-cracks. The shear behaviors of a DFN model at different normal

123

P. Wang et al.

(a)

(b)

0%

(c)

5%

(d)

30%

(e)

100%

(f)

150%

300%

Fig. 20 Numerical results with different K(S) values after failure under direct shear testing

Joint distribuon

Rock tunnel

Fractures

Joint idenficaon

Joint traces

(a) Rock face in a tunnel

(b) Joint identification

DFN for 2D analysis

Height

Joint sets Width

(c) Spatial distribution of the joints Fig. 21 Sketch map of a large spatial distribution of the joints in a tunnel

123

A Digital Image-Based Discrete Fracture Network Model and Its Numerical Investigation of…

stresses, friction coefficients and bond strength were then numerically examined using the contact bond model. The proposed technique could efficiently establish a DFN model in PFC2D based on the image of discrete fractures. In shear testing on the DFN model, a vortex phenomenon was observed in velocity vectors along the shear direction during shear loading. The shear cracks were distributed mainly along the shear direction with a rotated Z-shape, and the normal cracks occurred mainly along the joint elements. The peak shear strength was controlled by the normal stress, friction coefficient and material strength, as expected. The shear modulus of the DFN model increased as the normal stresses increased. A lower shear displacement at the peak strength was found when the normal stress was higher, and the peak and residual strength showed more sensitive responses when the normal stress was smaller. Higher normal stress had a greater effect on the initiation and propagation of cracks of DFN models. The friction coefficient was found to have a significant influence on the shear strength, with the number of shear cracks increasing linearly with increases in the friction coefficient. Moreover, the increase in friction coefficient decreased the initiation of normal cracks. The peak shear strength and residual strength increased with increases in the bond strength ratio, although the bond strength ratio did not have a measurable effect on displacement at peak shear strength. The increase in the bond strength ratio also resulted in an increase in the number of both shear and normal cracks. The digital image-based DFN model has at present a number of unresolved issues that have been highlighted throughout the paper: grain size and porosity, particle interlocking issues and micro-roughness created by large size particles, and the effect of scale on rock block have not been considered. As a consequence, there is a need for further development and validation. Acknowledgements This research was supported by the State Key Research and Development Program of China (Grant Nos. 2016YFC0600703 and 2016YFC801602), the National Natural Science Foundation of China (Grant No. 51604017), China Postdoctoral Science Foundation (Grant No. 2016M591079) and the Fundamental Research Funds for the Central Universities (Grant No. FRF-TP-15109A1). All of these supports are gratefully acknowledged. We thank Prof. Xu for editing and correcting the English. We also thank the editor, Prof. Barla, and the anonymous reviewers for greatly improving the manuscript. Conflict of interest The authors declare that they have no conflict of interest.

References Asadi SM, Rasouli V, Barla G (2012) A bonded particle model simulation of shear strength and asperity degradation for rough rock fractures. Rock Mech Rock Eng 45:649–675

Baghbanan A, Jing LR (2007) Hydraulic properties of fractured rock masses with correlated fracture length and aperture. Int J Rock Mech Min Sci 44:704–719 Bahaaddini M, Sharrock G, Hebblewhite BK (2013a) Numerical investigation of the effect of joint geometrical parameters on the mechanical properties of a non-persistent jointed rock mass under uniaxial compression. Comput Geotech 49:206–225 Bahaaddini M, Sharrock G, Hebblewhite BK (2013b) Numerical direct shear tests to model the shear behaviour of rock joints. Comput Geotech 51:101–115 Bahaaddini M, Hagan PC, Mitra R, Hebblewhite BK (2014) Scale effect on the shear behaviour of rock joints based on a numerical study. Eng Geol 181:212–223 Bahaaddini M, Hagan PC, Mitra R, Khosravi MH (2016) Experimental and numerical study of asperity degradation in the direct shear test. Eng Geol 204:41–52 Cai YD, Liu DM, Mathews PJ, Pan ZJ, Elsworth D, Yao YB, Li JQ, Guo XQ (2014) Permeability evolution in fractured coalcombining triaxial confinement with X-ray computed tomography, acoustic emission and ultrasonic techniques. Int J Coal Geol 122:91–104 Chen SJ, Zhu WC, Yu QL, Liu XG (2016) Characterization of anisotropy of joint surface roughness and aperture by variogram approach based on digital image processing technique. Rock Mech Rock Eng 49:855–876 Cho JW, Kim H, Jeon S, Min KB (2012) Deformation and strength anisotropy of Asan gneiss, Boryeong shale, and Yeoncheon schist. Int J Rock Mech Min Sci 50:158–169 Cnudde V, Boone MN (2013) High-resolution X-ray computed tomography in geosciences: a review of the current technology and applications. Earth Sci Rev 123:1–17 Dai QL (2011) Two and three-dimensional micromechanical viscoelastic finite element modeling of stone-based materials with X-ray computed tomography images. Constr Build Mater 25:1102–1114 Dershowitz WS, Einstein HH (1988) Characterizing rock joint geometry with joint system models. Rock Mech Rock Eng 21:21–51 Ding X, Zhang L, Zhu H, Zhang Q (2013) Effect of model scale and particle size distribution on PFC3D simulation results. Rock Mech Rock Eng 6:1–18 Feng XT, Pan PZ, Zhou H (2006) Simulation of the rock microfracturing process under uniaxial compression using an elasto-plastic cellular automation. Int J Rock Mech Min Sci 43:1091–1108 Ghabraie B, Ren G, Smith J, Holden L (2015) Application of 3D laser scanner, optical transducers and digital image processing techniques in physical modelling of mining-related strata movement. Int J Rock Mech Min Sci 80:219–230 Ghazvinian A, Sarfarazi V, Schubert W, Blumel M (2011) A study of the failure mechanism of planar non-persistent open joints using PFC2D. Rock Mech Rock Eng 45:677–693 Huang YJ, Yang ZJ, Ren WY, Liu GH, Zhang CZ (2015) 3D mesoscale fracture modelling and validation of concrete based on in situ X-ray computed tomography images using damage plasticity model. Int J Solids Strut 67:340–352 Hunt SP, Meyers AG, Louchnikov V (2003) Modelling the Kaiser effect and deformation rate analysis in sandstone using the discrete element method. Comput Geotech 30(7):611–621 Itasca (2004) Particle Flow Code in 2-Dimensions: Problem Solving with PFC2D, version 3.1. Itasca Consulting Group, Inc., Minneapolis, MN Ivars DM, Pierce ME, Darcel C, Reyes-Montes J, Potyondy DO, Young RP, Cundall PA (2011) The synthetic rock mass approach for jointed rock mass modelling. Int J Rock Mech Min Sci 48:219–244

123

P. Wang et al. Jia P, Tang CA (2008) Numerical study on failure mechanism of tunnel in jointed rock mass. Tunn Undergr Space Technol 23:500–507 Jiang MJ, Konrad JM, Leroueil S (2003) An efficient technique for generating homogeneous specimens for DEM studies. Comput Geotech 30:579–597 Jiang YJ, Tanabashi Y, Li B, Xiao J (2006) Influence of geometrical distribution of rock joints on deformational behavior of underground opening. Tunn Undergr Space Technol 21:485–491 Johansson F (2016) Influence of scale and matedness on the peak shear strength of fresh, unweathered rock joints. Int J Rock Mech Min Sci 82:36–47 Kim KY, Zhuang L, Yang H, Kim H, Min KB (2016) Strength anisotropy of Berea sandstone: results of X-ray computed tomography, compression tests, and discrete modeling. Rock Mech Rock Eng 49:1201–1210 Koyama T, Jing L (2007) Effects of model scale and particle size on micro-mechanical properties and failure processes of rocks—a particle mechanics approach. Eng Anal Bound Elem 5:458–472 Lambert C, Buzzi O, Giacomini A (2010) Influence of calcium leaching on the mechanical behavior of a rock-mortar interface: a DEM analysis. Comput Geotech 37(3):258–266 Li LC, Yang TH, Liang ZZ, Zhu WC, Tang CA (2011) Numerical investigation of groundwater outbursts near faults in underground coal mines. Int J Coal Geol 85:276–288 Li XW, Nemcik J, Mirzaghorbanali A, Aziz N, Rasekh H (2015) Analytical model of shear behaviour of a fully grouted cable bolt subjected to shearing. Int J Rock Mech Min Sci 80:31–39 Li LC, Xia YJ, Huang B, Zhang LY, Li M, Li AS (2016) The behaviour of fracture growth in sedimentary rocks: a numerical study based on hydraulic fracturing processes. Energies 9(3):169–197. doi:10.3390/en9030169 Liu XY, Zhang CY, Liu QS, Birkholzer J (2009) Multiple-point statistical prediction on fracture networks at Yucca Mountain. Environ Geol 57:1361–1370 Maleki MR (2011) Study of the engineering geological problems of the Havasan Dam, with emphasis on clay-filled joints in the right abutment. Rock Mech Rock Eng 44:695–710 Park B, Min KB (2015) Bonded-particle discrete element modeling of mechanical behavior of transversely isotropic rock. Int J Rock Mech Min Sci 76:243–255 Park JW, Song JJ (2009) Numerical simulation of a direct shear test on a rock joint using a bonded-particle model. Int J Rock Mech Min Sci 46:1315–1328 Potyondy DO (2015) The bonded-particle model as a tool for rock mechanics research and application: current trends and future directions. Geosystem Eng 18:1–28 Potyondy DO, Cundall PA (2004) A bonded-particle model for rock. Int J Rock Mech Min Sci 41:1329–1364 Reid TR, Harrison JP (2000) A semi-automated methodology for discontinuity trace detection in digital images of rock mass exposures. Int J Rock Mech Min Sci 37:1073–1089 Sarfarazi V, Ghazvinian A, Schubert W, Blumel M, Nejati HR (2014) Numerical simulation of the process of fracture of echelon rock joints. Rock Mech Rock Eng 47:1355–1371 Sun JP, Zhao ZY (2010) Effects of anisotropic permeability of fractured rock masses on underground oil storage caverns. Tunn Undergr Space Technol 25:629–637

123

Tang CA, Liu H, Lee PKK, Tsui Y, Tham LG (2000) Numerical studies of the influence of microstructure on rock failure in uniaxial compression-part I: effect of heterogeneity. Int J Rock Mech Min Sci 37:555–569 Tran TH, Ve´nier R, Cambou B (2009) Discrete modelling of rockageing in rockfill dams. Comput Geotech 36(1–2):264–275 Vallejos AJ, Suzuki K, Brzovic A, Ivars DM (2016) Application of synthetic rock mass modeling to veined core-size samples. Int J Rock Mech Min Sci 81:47–61 Wang C, Tannant DD, Lilly PA (2003) Numerical analysis of the stability of heavily jointed rock slopes using PFC2D. Int J Rock Mech Min Sci 40:415–424 Wang PT, Yang TH, Yu QL, Liu HL, Zhang PH (2013a) Characterization on jointed rock masses based on PFC2D. Front Struct Civ Eng 7:32–38 Wang PT, Yang TH, Xu T, Yu QL, Liu HL (2013b) A model of anisotropic property of seepage and stress for jointed rock mass. J App Math 2013:1–19 Wang PT, Yang TH, Xu T, Cai MF, Li CH (2016) Numerical analysis on scale effect of elasticity, strength and failure patterns of jointed rock masses. Geosci J 20(4):539–549 Wasantha PLP, Ranjith PG, Zhang QB, Tao Xu (2015) Do joint geometrical properties influence the fracturing behavior of jointed rock? An investigation through joint orientation. J Geomech Geophys Geol Energy Geol Res 1:3–14 Weibull W (1951) A statistical distribution function of wide applicability. J App Mech 18:293–297 Xu WJ, Yue ZQ, Hu RL (2008) Study on the mesostructure and mesomechanical characteristics of the soil–rock mixture using digital image processing based finite element method. Int J Rock Mech Min Sci 45:749–762 Xu T, Ranjith PG, Wasantha PLP, Zhao J, Tang CA, Zhu WC (2013) Influence of the geometry of partially-spanning joints on mechanical properties of rock in uniaxial compression. Eng Geol 167:134–147 Yang TH, Wang PT, Xu T, Yu QL, Zhang PH, Shi WH, Hu GJ (2015) Anisotropic characteristics of fractured rock mass and a case study in Shirengou Metal Mine in China. Tunn Undergr Space Technol 48:129–139 Yu QL, Yang SQ, Ranjith PG, Zhu WC, Yang TH (2016) Numerical modeling of jointed rock under compressive loading using X-ray computerized tomography. Rock Mech Rock Eng 49:877–891 Yun TS, Jeong YJ, Kim KY, Min KB (2013) Evaluation of rock anisotropy using 3D X-ray computed tomography. Eng Geol 163:11–19 Zhang HQ, Zhao ZY, Tang CA, Song L (2006) Numerical study of shear behavior of intermittent rock joints with different geometrical parameters. Int J Rock Mech Min Sci 43:802–816 Zhang ZX, Hu XY, Scott KD (2011) A discrete numerical approach for modeling face stability in slurry shield tunnelling in soft soils. Comput Geotech 38(1):94–104 Zhu WC, Liu J, Yang TH, Sheng JC, Elsworth D (2006) Effects of local rock heterogeneities on the hydromechanics of fractured rocks using a digital-image-based technique. Int J Rock Mech Min Sci 43:1182–1199