Study on the Fracture Distribution Law and the Influence of Discrete Fractures on the Stability of Roadway Surrounding Rock in the Sanshandao Coastal Gold Mine, China

Abstract

Cracks are critical for the deformation and failure of rock masses, but the effects of real cracks are rarely considered when evaluating the stability and safety of practical engineering. This paper presents a study on the application of fractures in the Sanshandao Gold Mine. Field investigation and statistical analysis methods were adopted to obtain the distribution laws of the cracks. Laboratory tests, MATLAB programming, and simulation using the software, GDEM (Gdem Technology, Beijing, China, Co., Ltd.), were employed to study the mechanical behaviors of rock masses with real fractures after excavation. The main results are as follows: (1) Three sets of highly discrete cracks were developed in the study area. Their inclination and dip can be approximately considered to follow a Gaussian distribution or uniform distribution. They had close ties to the three faults developed in the mining area. (2) Compared with the model that did not consider cracks and the model processed by the equivalence idea, the surrounding rock deformation caused by excavation of the model that considered real cracks was larger than that of the former and smaller than that of the latter. However, its influence range was larger than that of the other two models. The results show that it is reasonable to use three sets of discrete cracks to characterize the fracture distribution of the surrounding rock. In the evaluation of roadway stability, it is not advisable to use the equivalence method to deal with all the cracks. Considering a part of the cracks that are compatible with the size of the calculation model, a relatively accurate evaluation can be obtained in terms of the deformation, failure, and permeability changes of the surrounding rock.

1. Introduction

Rock mass is a natural geological body, composed of rock and structural planes. As a common structural plane, cracks play an important role in controlling the deformation and failure of rock mass, which seriously affects the stability of the surrounding rock in engineering. Additionally, cracks can provide moving space for fluids, posing safety hazards to many water-related human engineering activities. Therefore, finding the distribution laws of cracks and studying the development and evolution of them in human engineering are of great significance for safe production in mine engineering and the sustainable use of water conservation and hydropower engineering.

As cracks play a key role in ergonomic activities, many scientists have conducted a lot of research on the quantitative descriptions of the geometric elements of cracks and their statistical distribution laws [1,2,3,4]. Li et al. [5] used nuclear magnetic resonance (NMR) and X-ray computed tomography (X-CT) to evaluate the characteristics of the fractures of coal. Xie et al. [6] studied the quantitative description of the discrete fracture network (DFN), based on the statistical reconstruction method. Kulatilake and Wu [7] proposed a window statistical method to estimate the mean trace length on finite exposures. Priest and Hudson [8] pointed out that the fracture spacings obey the negative exponential distribution. Through a statistical analysis of real fractures, Wu [3] summarized that the fracture occurrences often obey a truncated normal distribution or log-normal distribution. On the other hand, scientists have carried out application research on cracks in engineering rock mass [9,10,11,12]. Song et al. [13] studied the influence of the spacing and azimuth of cracks on the stability of the roadway by laboratory tests. Elmo et al. [14] simulated the progressive failure process of natural fractured pillars by the hybrid finite–discrete element method, code: ELFEN. According to the statistical laws of the fractures of marble rocks, Huang et al. [15] used statistical methods to solve the representative elementary volume (REV) of discrete fractured rock masses. Combining DFN with the discontinuous deformation analysis (DDA) method, Liu et al. [16] analyzed the slope stability in water conservation and hydropower engineering.

In general, the engineering application research on jointed rock mass has the following characteristics: (1) Accurately obtaining a large number of geometry data of cracks is to some extent difficult. (2) Laboratory experiments have technical difficulties in portraying the DFN. Thus, the study of the discrete fractures mainly focuses on mathematical models and numerical calculations. (3) The DFN models are mostly used for seepage calculations, and the simulations of models for mechanical or hydro-mechanical coupling calculation usually employ the equivalent continuum approach that almost does not consider real cracks. (4) Conventional numerical analysis methods have many difficulties in dealing with cracks. Thus, new numerical methods, developed in recent years, such as the DDA and hybrid finite–discrete element methods/continuum–discontinuum element methods (FDEM, CDEM), should be further used to simulate the failure of jointed rock mass.

Numerous fractures, including faults and cracks, are developed in the Sanshandao Gold Mine. Meanwhile, there is abundant sea water above the mining area. Therefore, the study of the Sanshandao Gold Mine is of great significance for the safety of its production. Firstly, 1120 cracks were obtained by field investigation of the surrounding rock of roadways in the study area for many years. Then, the crack data were grouped and analyzed by statistical methods, and a DFN model was generated using the Monte Carlo method, which was implemented in MATLAB. A more reasonable meshing method was adopted when dealing with multi-fracture geometric models in the meshing software, Gmsh. Physical and mechanical parameters were obtained from laboratory tests. Finally, three models, including an intact rock model, DFN model, and equivalent continuum model were built to investigate the deformation, failure, and influence range of the surrounding rock of the roadway using the software, GDEM, which can consider initial cracks and the propagation and failure evolution of new fractures. The flow chart of this research can be described as shown in Figure 1.

2. Crack Investigation and Distribution Law of the Mining Area

2.1. Geological Background

The Sanshandao Gold Mine is located in Sanshandao Village, Laizhou City, Shandong Province, China. It is surrounded by the Bohai Sea on three sides and is connected to the land only in the southeast. The flat regional terrain is high in the east and south, low in the west and north, and dips to the northwest. The elevation of the ground is generally 2–60 m. The ground surface is mainly covered by Quaternary deposits, underlying the Neoarchean Jiaodong rock group including granite, granodiorite, gneiss, and monzonitic granite [17].

The geological structure of the mining area is dominated by fractures, and three faults, named F1, F2, and F3, are developed (Figure 2). Fault F1 is a compression–torsion fault, which controls the distribution of ore and is located in the northeastern section of the Sanshandao–Cangshang fault zone. It is covered by Quaternary deposits, tending to the southeast, with a dip of 33°–67°, which is mostly between 40°–50°. There is a black fault gouge, with a thickness of 10–40 cm, which developed in the fault zone, and a cataclastic rock belt, with a thickness of 1–10 m developed on both sides of the main section of the fault. Fault F2 is a left-handed torsional fault, with a length of about 600 m. It tends to the northwest, with a dip of 85°, and extends into the Bohai Sea in the north. It also crosses the pyrite quartz vein at a distance of 20 m. The torsional fault F3 is 1500 m long in the mining area and belongs to the northwestern end of the Sanyuan–Chenjia fault, extending westward to the Laizhou Bay. The strike of fault F3 is 300°, and it is almost upright. The faultage fracture zone is about 10–25 m wide and consists of several strips of diabase veins, cataclastic rocks, and breccia.

2.2. Crack Distribution

The occurrence of cracks in the surrounding rock of the main roadways has been measured by a geological compass for many years, and the inclination and dip data of 1120 cracks were obtained (the data are not listed one by one due to the limited length of the article). A pole isopycnal map of the cracks was drawn, based on the inclination and dip of all the fractures (Figure 3). As shown in Figure 3, in terms of inclination, more fractures developed in the second and fourth quadrants than in the first and third quadrants. The dip range of the former was wider than in the latter. Combining the occurrence of the three main faults in the mining area, the fractures can be roughly divided into three sets. The inclination of Set 1 of the fractures was distributed in the fourth quadrant, which was related to fault F1. The inclination of Set 2 of the fractures was distributed in the second quadrant, which was related to fault F2, and the inclination of Set 3 of the fractures was distributed in the first and third quadrants and had a large dip, which was related to the nearly erect fault F3. As shown in Figure 4, two frequency histograms of the inclination of the cracks were plotted, based on the inclination data. One contained all the data (Figure 4a), and the other only contained part of the data (Figure 4b). Figure 4a shows that the relative frequency had two maximum values at the positions of 120° and 300°, corresponding to Set 1 and Set 2, respectively. However, the relative frequency value was almost constant at the positions of 0°–70°, 170°–250°, and 350°–360°. This may be affected by Set 3. Thus, the cracks with a dip greater than 70° in this range are temporarily removed, and the remaining cracks were used to draw a frequency histogram of the inclination (Figure 4b). The regularity of fracture distribution in Figure 4b was much stronger than that in Figure 4a. Therefore, it can be concluded that the 216 excluded data belonged to the third group of cracks, while the remaining 904 cracks originated from the first and second groups. Therefore, it can be considered that the 216 removed cracks were from Set 3, and the remaining 904 cracks were from Set 1 and Set 2.

The frequency distribution histograms of the inclination of Set 1, Set 2, and Set 3 are plotted in Figure 5. As shown in Figure 5, the inclination of the first two sets of fractures can be approximately considered to follow a truncated normal distribution. The expectation and standard deviation of Set 1 were 123° and 28.7°, respectively (Figure 5a), and the inclination was generated in the range of 30°–200°. The expectation and standard deviation of Set 2 were 297° and 38.2°, respectively (Figure 5b), and the inclination was generated in the range of 0°–30° and 200°–360°. The inclination of Set 3 had a uniform distribution (Figure 5c), and the inclination was generated in the range of –10°–70° (merged by 350°–360° and 0°–70°) and 170°–250°.

The frequency distribution histograms of the dip of Set 1, Set 2, and Set 3 are plotted in Figure 6. As shown in Figure 6, the dip of the two sets of fractures can be approximately considered to follow a normal distribution in the interval of 0°–90°. The expectation and standard deviation of Set 1 were 68° and 25.2°, respectively (Figure 6a). The expectation and standard deviation of Set 2 were 62° and 18.6°, respectively (Figure 6b). The dip of Set 3 had a uniform distribution (Figure 6c), and the dip was generated in the range of 70°–90°.

Most importantly, the distribution laws of the three sets of fractures are summarized in

Table 1. Distribution laws of three sets of cracks.

Set of Cracks		Set 1	Set 2	Set 3
Number of Cracks		451	453	216
Inclination	Distribution law	Gaussian normal	Gaussian normal	Uniform
	Range (°)	(30, 200)	(0, 30), (200, 360)	(−10, 70), (170, 250)
	Expectation (°)	123	297	-
	Standard deviation (°)	28.7	38.2	-
Dip	Distribution law	Gaussian normal	Gaussian normal	Uniform
	Range (°)	(0, 90)	(0, 90)	(70, 90)
	Expectation (°)	68	62	-
	Standard deviation (°)	25.2	18.6	-

.

3. Physical and Mechanical Properties of the Surrounding Rocks

3.1. Laboratory Tests

The main lithology of the surrounding rock in the roadway is the sericite-quartz granitic cataclastic rock. Laboratory-scale samples were obtained by drilling sampling and rock sample processing. The size of these samples satisfied the following: (1) the length of the cylinder was at least 10 cm, and the ratio of the height–to–diameter was between 2 and 2.5; (2) the ratio of the height–to–diameter of the Brazilian disc was between 0.5 and 1.0. Density measurement, uniaxial compression, Brazilian splitting, and triaxial compression tests were carried out. The density was obtained from the average of the two samples. Figure 7 shows the testing machine used in the uniaxial compression, Brazilian splitting, and triaxial compression tests, as well as the rock samples before and after the tests. The uniaxial compression test was carried out on the WES-2000 digital display hydraulic universal testing machine (Figure 7a). Three cylindrical samples were tested to obtain the uniaxial tensile strength, elastic modulus, and Poisson’s ratio of the rock. The Brazilian splitting test was carried out on the WEP-600 universal testing machine (Figure 7b). Three disc samples were tested to obtain the tensile strength of the rock. The triaxial compression test was carried out on a TYS-500 rock triaxial test machine (Figure 7c). Five cylindrical samples were tested, with a confining pressure of 0, 5, 10, 15, and 20 MPa. The cohesion and internal friction angle of the rock were obtained by plotting 5 Mohr circles and fitting them with the strength envelope. The measured data are shown in

Table 2. Physical and mechanical properties of the surrounding rock.

Lithology	Density (g·cm–3)	Brazilian Split Test	Uniaxial Compression Test			Triaxial Compression Test
		Tensile Strength (MPa)	Uniaxial Compression Strength (MPa)	Young’s Modulus (GPa)	Poisson’s Ratio	Cohesion (MPa)	Friction Angle (°)
Sericite-Quartz Granitic Cataclastic Rock	2.627	3.43	27.73	7.13	0.120	10.98	46.25
	2.620	5.69	41.97	23.38	0.197
	-	6.62	63.24	22.08	0.348
Mean Value of Properties	2.624	5.25	44.31	17.53	0.222	-	-

.

3.2. Parameter Calibration in the Numerical Model

GDEM, developed on the basis of the continuum–discontinuum element method (CDEM), is a parallel computing numerical simulation software. Relying on a high-performance graphics processing unit (GPU), the software can increase the computational speed by two orders of magnitude and increase the computational scale by one order of magnitude. GDEM can simulate the progressive failure process of geological bodies and evaluate the local or global stability of materials in geotechnical engineering. It is widely used in mining, tunneling, water conservation engineering, and other fields. By coupling the finite element with the discrete element, the GDEM performs a finite element calculation inside the block and discrete element calculation at the block boundary. Moreover, it realizes the fracture and slippage of the block through a rupture of the normal and the tangential springs at the boundary [18] (Figure 8).

In addition to the density, elastic modulus, Poisson’s ratio, tensile strength, cohesion, and internal friction angle contained in Table 2, the GDEM also requires the parameters of the springs placed between the block elements, including the normal stiffness, tangential stiffness, friction, cohesion, tensile strength, and fracture energy. These parameters cannot be obtained directly from tests or are difficult to obtain. Therefore, they need to be calibrated according to the results of the laboratory test. That was a trial and error process. Taking the uniaxial compression test as an example in this paper (Figure 9), the stress–strain curves of numerical simulation and laboratory tests are shown in Figure 10. As shown in Figure 9a, a rock model with a width of 50 mm and a height of 100 mm was built for the uniaxial compression test. The bottom of the model was fixed, and the top of the model was loaded down at a constant speed of 2e–8 m/s. The final failure trajectories were formed as shown in Figure 9b, and the penetrating rupture surface at an angle in relation to the loading direction was substantially consistent with the results of the laboratory tests. As shown in Figure 10, the trend of the stress–strain curves of the simulation and laboratory tests were roughly the same, and the Young’s modulus (slope of the curve) of the simulation and uniaxial compression strength was basically equal to the corresponding average values of the laboratory tests. Thus, the parameters shown in Table 2 and

Table 3. Parameters of the springs placed between the block elements.

Parameter	Normal Stiffness (Pa/m)	Shear Stiffness (Pa/m)	Friction Angle (°)	Cohesion (MPa)	Tensile Strength (MPa)	Tensile Fracture Energy (Pa·m)	Shear Fracture Energy (Pa·m)
Value	2 × 1013	2 × 1013	40	9.8	4.8	20	140

can be considered as reasonable values of the rock samples at the laboratory scale.

4. DFN Model

4.1. Geometry Model

The inclination and dip distribution law of the three sets of cracks were obtained by the statistics of the fracture data. In determining the number of fractures generated in each set, the ratio of the number of each set to the number of all fractures and the calculation efficiency of the GDEM was considered. It was finally determined that 200, 200, and 100 cracks of the three sets were generated, respectively, in the calculation model with a size of 80 m × 80 m. It was assumed that the length of the trace obeyed a negative exponential distribution [3], and the trace length of cracks was generated in the range of 3 m to 11 m. Based on the Monte Carlo method, the generation and display programs of random cracks were implemented in MATLAB. The generated cracks were shown in Figure 11. The two-dimensional (2D) surface in Figure 11 was perpendicular to the strike of the horizontal roadway, and the dip of the cracks was converted from the true dip to the apparent dip on the surface.

4.2. Meshing of the Model

The calculation module of GDEM can handle the model with many discrete fractures very well, and the calculation efficiency is consistent with the model without these structural planes. However, the pre-processing function of GDEM is relatively weak, making it difficult to generate a computational model with complex cracks. Therefore, GDEM allows the grids generated by other software, which have a powerful pre-processing function (e.g., ANSYS) or other dedicated meshing software (e.g., GID and Gmsh), to be imported into it. Among these types of software, Gmsh can consider non-through cracks when meshing, and as Gmsh is a command operation software, it is convenient to use the point coordinates generated by MATLAB. However, when dealing with interacting cracks, GDEM usually does not assign mesh nodes at the intersections of cracks. This will cause several cracks to be broken into some disconnected fragments, as shown in Figure 12a,b, and the model will not be meshed when several cracks intersect with the boundary lines. Therefore, it was necessary to find the intersections between the cracks and boundary lines and the intersections between the cracks to make the meshing work successfully and achieve high-quality grids. Then, the boundaries and cracks were divided at the intersections (if they exist), and all the points and lines were renumbered. Finally, they were imported into Gmsh for meshing, as shown in Figure 12c,d. For a few cracks, the intersections can be calculated and renumbered manually. However, for a large number of cracks, computer programming should be used to improve the work efficiency. In this article, the corresponding calculation procedure was developed in MATLAB to facilitate the meshing of the geometry model.

The generated geometric model was imported into the software Gmsh for meshing. Then, a numerical calculation model with a size of 80 m × 80 m was obtained, which contained a roadway with a diameter of 4 m in the center of the model, as shown in Figure 13.

5. Influence of the DFN Model on Roadway Stability

5.1. Parameters and Scheme Designs for the Simulation

To study the influence of DFN on the stability of the main permeable midsection with a depth of 510 m, three different numerical simulation schemes were designed. Scheme 1 was a numerical calculation model without discrete fractures; Scheme 2 was a numerical calculation model, which had an added DFN, based on Scheme 1; Scheme 3 was a numerical calculation model, which weakened the rock mass parameters, based on Scheme 1. The parameters for the three schemes are shown in

Table 4. Physical and mechanical parameters for the three simulation schemes.

		Parameters for Scheme 1 and 2	Parameters for Scheme 3	Parameters of DFN for Scheme 2
				Set 1 and 2	Set 3
Block Elements	Young’s modulus (GPa)	1.75	0.88	-	-
	Poisson’s ratio (-)	0.20	0.20	-	-
	Density (g·cm−3)	2.60	2.60	-	-
	Cohesion (MPa)	1.01	1.00	-	-
	Friction angle (°)	35	25	-	-
	Tensile strength (MPa)	0.52	0.50	-	-
Springs	Cohesion (MPa)	1.10	1.00	1.40 × 10−2	7.00 × 10−3
	Friction angle (°)	30	22	10	10
	Tensile strength (MPa)	0.52	0.50	7.00 × 10−3	3.50 × 10−3
	Normal stiffness (Pa/m)	2 × 1012	1 × 1012	5 × 1010	5 × 1010
	Shear stiffness (Pa/m)	2 × 1012	1 × 1012	5 × 1010	5 × 1010
	Tensile fracture energy (Pa·m)	10.0	8.0	0.2	0.1
	Shear fracture energy (Pa·m)	60.0	48.0	1.2	0.6

. Compared with the rock parameters obtained in the laboratory tests, the parameters used in Schemes 1 and 2 were somewhat weakened. The reduction of the parameters was mainly due to two factors. The first was the size effect of the rock [20,21]. The sample size for the laboratory test was much smaller than the rock size for the practical engineering. Thus, the parameters obtained from the laboratory test should be adjusted when applied to engineering. The second factor was the existence of small cracks or defects in the engineering rock mass. While the discrete fractures were set in the model of Scheme 2, the size of these cracks was limited, which meant that only the cracks matching the size of the model were considered. The parameters of the model that contained smaller cracks needed to be reduced by the equivalence method. The parameters of the discrete fractures were drawn from a previous paper [22]. Similarly, the parameters were also obtained by repeated trial and error, which should ensure that the magnitude of the simulated deformation of the roadway was consistent with the actual deformation. As fault F3 is a tensile fracture, the parameters of Set 3 were relatively weaker than those of the other two sets in the numerical simulation. The degree of reduction of the rock mass parameters was determined according to a large number of trials, and by comparing the roof’s displacement of the roadway in the numerical simulation and practical monitoring. Compared with Scheme 1, the rock parameters in Scheme 3 were further weakened. It can be considered that, when the numerical model with cracks was equivalently replaced with a homogeneous and isotropic rock mass, only small fractures were considered in Scheme 1, while both small and large fractures were considered in Scheme 3.

The three schemes had the same boundary conditions. The normal displacement of each boundary of the model was limited. Horizontal and vertical geostresses, with gradients in the vertical direction, were applied to the model [23]. The horizontal geostress gradually increased from the top to the bottom of the model. Its magnitude at the top was 21.91 MPa, and the gradient was 0.0449 MPa. The gradient of the vertical geostress was mainly due to gravity. Its magnitude at the top was 12.26 MPa.

5.2. Modeling Results and Analyses

The similarity between Schemes 1, 2, and 3 was that the cracks that were smaller than the discrete fractures were treated equivalently. The difference is that the description of larger cracks was discrepant. Therefore, it can be found that, after considering the size effect and the equivalence of smaller cracks, Scheme 1 describes the intact rock without fractures, Scheme 2 describes the rock mass with real fractures, and Scheme 3 describes the rock mass considering the equivalence of the real fractures.

5.2.1. Deformation Features of the Roadway

Whether the surrounding rock of the roadway is stable or not will be directly reflected in the deformation of the roadway. Since the degree of necking of the roadway in the vertical direction was larger than that in the horizontal direction in the calculation process, the vertical displacement cloud diagrams of the three schemes were drawn, as shown in Figure 14.

Comparing Scheme 1 (Figure 14a) and Scheme 2 (Figure 14b), the deformations at the bottom of the roadway were not very different. However, the deformation at the top of the roadway in Scheme 2 was larger than that in Scheme 1. This was because there were weak faces around the roadway before excavation, and the faces were first broken after excavation. These broken faces acted as small faults, which had the properties of a displacement barrier effect and local amplification effect. Thus, the deformation was concentrated in the area between the weak faces and the roadway.

Comparing Scheme 1 (Figure 14a) and Scheme 3 (Figure 14c), the deformation of the roadway in Scheme 3 was much larger than that in Scheme 1. From a shallow perspective, this was due to the decrease of parameters such as the elastic modulus. From a deeper perspective, Scheme 3 considered a greater range of crack sizes and treated them equivalently. It was proved that the equivalent treatment of cracks was reasonable. However, the correctness of the method that transformed the model with large fractures equivalently into a homogeneous and isotropic rock mass, as well as the degree of reduction of the rock mass parameters, still needs further study.

Comparing Scheme 2 (Figure 14b) and Scheme 3 (Figure 14c), the deformation of the roadway in Scheme 3 was much larger than that in Scheme 2. This was because the size of cracks considered in Scheme 2 was limited, while Scheme 3 considered a greater range of crack sizes, and the maximum length of cracks in Scheme 3 may have already exceeded the scope of the model. If such a transfixion plane exists, the deformation along the plane is bound to be very large. However, the equivalence idea makes the directional deformation along the crack become a uniform necking deformation. Meanwhile, the equivalence idea will increase the overall deformation of the rock evenly, but delay the failure of the rock mass. If larger cracks were still treated with the equivalence idea, the anisotropy of the rock mass deformation and its control action on the local deformation and failure cannot be reflected well. Therefore, it is necessary to consider larger cracks as real fractures rather than equivalent treatments.

The displacements of key positions of the roadway deformation were monitored during the simulation. The vertical displacements of the top and bottom (positions A and B), and the horizontal displacements of the left and right (positions C and D) in the three schemes were recorded separately, as shown in Figure 15. The displacements were taken in the absolute value for ease of comparison in Figure 15. Both Scheme 1 (Figure 15a) and Scheme 3 (Figure 15c) did not consider real cracks and the deformation law, exhibited by their four monitoring points, was consistent. The order of deformation values, from large to small, was A, B, C, and D. The reason why the displacement of A was larger than that of B is that damage occurred around the roadway. Under the action of gravity, the top of the roadway moved toward the free surface, resulting in a larger displacement of position A. The reason why the displacements of A and B were generally larger than that of C and D is that the damage range in the vertical direction was larger than that in the horizontal direction (Figure 15). The difference in deformation between C and D may be related to the size and shape of the mesh. In a case where the stress conditions and geometric conditions are symmetrical, the displacements of C and D should be equal. This can be illustrated by the fact that the curves of C and D gradually became closer together in Figure 15c. However, the deformation law shown by Scheme 2 (Figure 15b) was different. The order was A, D, B, and C. Due to the presence of weak faces near the upper right side of the roadway, the displacements of A and D were amplified, which fully illustrated the importance of discrete cracks.

5.2.2. Influence Scope of Excavation

As shown in Figure 16, the failure modes of the fractures of the three schemes were plotted. There were three failure modes in the three schemes: tensile failure, shear failure, and combined tensile–shear failure. Under the same stress and boundary conditions, the influence scopes caused by excavation of the three schemes were different. Due to the existence of cracks, the influence range of Scheme 2 (Figure 16b) was larger than that of Scheme 1 (Figure 16a). While the mechanical properties of Scheme 3 were weaker than that of Scheme 2, the influence range of the former was much smaller than that of the latter, considering the actual fracture. That is to say, weakening the parameters of the rock only increased the influence range around the roadway. It did not reflect the differential damage caused by the discrete cracks to the rock mass. The deformation of the cracks in Scheme 2 was not very large, but the failure should not be neglected.

Figure 17 shows a nephogram of the tensile damage factor, minimum principal stress, shear damage factor, and shear stress of Scheme 2. The length of each subplot is six times as long as the diameter of the roadway. The discrete fractures set in the model and the new cracks, generated after excavation, had different degrees of tensile damage (Figure 17a) and shear damage (Figure 17c). As shown in Figure 17b,d, the minimum principal stress was concentrated at the intact rock between the fractures and the shear stress was concentrated at the crack tip. The phenomenon of stress concentration indicates that the rock bridge between the fractures will continue to be destroyed under the long-term effect of mining disturbance or geostress. Permeable channels may appear in the jointed rock masses, resulting in an increase in the permeability of the surrounding rock of the roadway. Therefore, it is crucial to take the real cracks into account for the stability of the surrounding rock of the mine that is affected by a large amount of seawater and bedrock fissure water.

6. Summary and Conclusions

The attitudes of 1120 cracks were obtained according to a large number of actual investigations on cracks in the surrounding rock of the Sanshandao coastal Gold Mine. Moreover, the distribution law of these cracks was obtained by statistical methods. A 2D geometric model of the roadway surrounding rock with DFN was generated by MATLAB programming, and the model was imported into the software, Gmsh, for meshing. The deformation and failure of the surrounding rock, caused by excavation, were simulated by the software GDEM. The main conclusions obtained are as follows:

(1) There are three sets of cracks in the surrounding rock of the main roadway in the study area, corresponding to three major faults, F1, F2, and F3, developed in the mining area. The inclination of the first two sets is southeast and northwest, respectively, and both of them have a large dip. Their inclination and dip can be approximately considered to follow a truncated Gaussian distribution. The inclination of the third set, which is almost erect, is southwest and northeast, and their inclination and dip can be approximately considered to follow a uniform distribution. While the dispersion of the three sets of fractures is high, the distribution law provides the possibility for the quantification of fractures and numerical calculation of the stability of surrounding rock considering real fractures.

(2) Fractures are crucial in the stability evaluation of the surrounding rock of the roadway. The deformation of the surrounding rock with the real cracks is larger than that without cracks under the same conditions. These fractures have, to some extent, a barrier effect on the propagation of deformation, and the presence of them will amplify the deformation of the rock between the cracks and roadway, thereby accelerating the destruction of the surrounding rock. Meanwhile, the failure of the real cracks directly increases the influence scope caused by excavation. It is not only the change in displacement, but more importantly, the change in the rock mass failure direction and rock mass permeability in that direction. The penetrating permeable channels are fatal for the safety of the production in mines, especially the Sanshandao Gold Mine, which has abundant seawater over the mine. The distribution of cracks should be considered when designing a water plugging scheme.

(3) The dual media model should be used to evaluate the stability of the engineering rock mass. Small cracks are treated with an equivalent continuum model and large cracks are treated with a discrete fracture network model. The directionality of deformation, failure, and permeability could not be reflected using only the equivalent continuum method. Due to the limited computing power of the current commercial software, and the inability to obtain cracks of various sizes by field investigation, it is unachievable to evaluate the stability of the surrounding rock using only the DFN method. In case the computing power is satisfied, more multi-scale cracks should be considered as much as possible. Even if the statistical results of the cracks have a certain degree of deviation from the distribution of the crack in the actual rock mass, the evaluation results of the model with discrete fractures are more reasonable than that of the weakened homogeneous and isotropic model.