Simulation Method and Application of Three-Dimensional DFN for Rock Mass Based on Monte-Carlo Technique

Abstract

In this study, the authors simulate a polygonal discrete fracture network (DFN) in rock masses. The probability models of the relevant geological parameters, including the orientation, trace length, volume density, and coordinates of the centroid, are firstly developed as fractures are in the shape of rectangles. In the process, the probability distribution of rectangular fractures with side lengths as random variables is introduced and described in terms of mean trace lengths on the basis of the probability model of disk-shaped fracture with the diameter as the random variable. The relationship between the volume density and the linear density of rectangular fractures is given for a negative exponential distribution. Following this, the coordinates of the vertices of fractures are derived based on spatial algebraic geometry, and the data for the three-dimensional DFN model are generated using the Monte-Carlo technique. The resulting three-dimensional DFN is visualized by calling the Open GL graphics database in the environment of Visual C, and the process of implementation of the DFN simulation is given. Finally, the validity of the simulation is verified by applying it to engineering practice.

1. Introduction

Discontinuity is a significant mechanical break or fracture of negligibly small tensile strength in a rock, and it has a lower shear strength and higher conductivity of fluids than the rock itself. The properties of discontinuities are closely related to the mechanical properties of the rock mass, and are often used as a controlling factor to determine the stability and safety of rock engineering [1,2,3]. Therefore, the development of a reasonable and effective DFN model of rock mass has attracted research interest because it can provide a clear distribution of fractures in the context of stability analysis and calculations in rock engineering.

In general, the DFN model consists of three basic parts, an orientation distribution model of the discontinuities, a dimension model of the discontinuities and a spatial distribution model of the discontinuities [4,5]. For the dimension model, although many scholars have attempted to describe the relationship between the trace of the outcrop of discontinuities and the dimension of discontinuities, our understanding of the shape of the fracture is not entirely clear due to the complexity of fracture formation and the limitations associated with field observations [6,7]. Two main types of models have been used thus far to formulate the DFN: the Baecher disk model and the polygon model. The former was introduced by Baecher et al. [8], and is characterized by the assumption that the fracture is shaped as a disk or an ellipsoidal disk while ignoring its thickness. The polygon model was derived by Veneziano [9] in the process of generating a fracture model through a system of Poisson lines. Dershowitz and Einstein [10] reached a conclusion consistent with that of Veneziano. They assumed that the fracture is a polygon with three to twenty sides, the side lengths of these polygons can be equal or extensible, and the lengths of these extensions can be defined by the aspect ratio. Because the modeling is simple and can satisfy the needs of the project, the Baecher disk model has been widely used in the analysis of rock mechanics [11,12,13,14,15,16,17].

However, the rock mass is an inhomogeneous, discontinuous, and anisotropic medium. The shape of the fracture is related to the fracture mechanics and energy at the tip of the crack. It usually extends to form polygons rather than disks. Moreover, field investigations have indicated that a polygon is a more appropriate shape to represent fractures. Zhang and Xu [18] assumed that the shape of a single fracture tends to be circular in homogeneous rock masses, but when the rock mass contains multiple fractures and its anisotropy is prominent, the fractures should be assumed to be polygonal due to the influence of mutual cutting among them. Li et al. [19] coupled a network of polygonal fractures with a large-scale geological model to establish a refined model of the rock structure to identify and analyze random blocks. Fu et al. [20] introduced a statistical method for generating 3D polygonal fractures in order to identify a Complex Block System. Pan et al. [21] used an improved control circle algorithm to describe the spatial dimensions and shapes of fractures, and used this to establish a fracture network to provide a base model for the simulation-based analysis of the diffusion of grout. Graciela et al. [22] proposed a new method to generate realistic 3D Stochastic DFN by getting fracture plane points from Virtual Outcrop Models and defining clusters using a K-means algorithm. Current field observations suggest that polygonal fractures often do not have more than six sides. The number of sides of a fracture in polygon models proposed by researchers is also generally set at between four and six [5,19,23].

In this paper, considering the operability of the simulation method of the fracture network, a DFN model of rock mass with rectangle as the fracture shape is proposed by integrating the Monte-Carlo stochastic method, spatial algebraic geometry theory and 3D graphic visualization technology based on Open GL. It is structured as follows: in Section 2, the simulation method of three-dimensional DFN based on the Monte-Carlo technique is illustrated. Specifically, probability models of geological parameters on fractures are established, including the model of fracture orientation, dimension, density, and center point. The simulation process of three-dimensional DFN based on Monte-Carlo technology is displayed. In Section 3, the implementation process of the simulation method of three-dimensional DFN is given, and the above DFN model is applied on engineering practice. In Section 4, the discussion and conclusions are presented.

2. Simulation Method of Three-Dimensional DFN Based on Monte-Carlo Technique

The spatial geometric characteristics of rock fractures are very complex, and in order to facilitate the generation of the DFN based on the probability models of geological parameters, the following assumptions are adopted: (1) the upper and lower boundaries of the rectangular fracture are parallel to the horizontal plane, that is, the two endpoints of the upper boundary have equal Z-coordinate values, and the two endpoints of the lower boundary also have equal Z-coordinate values; (2) the variables of two side lengths of the rectangular fracture are independent, and obey the same probability distribution.

2.1. Construction of Probability Models for Geological Parameters of Rock Fractures

2.1.1. Model of Fracture Orientation

Einstein and Baecher [24] summarized five forms of the distribution of fracture orientation: uniform distribution, Fisher distribution, elliptical distribution, Bingham distribution, and normal distribution. Wang et al. [25] examined statistics on the fracture orientation of several large hydraulic projects in China, and concluded that the angle and direction of the dips of fractures are independent of each other, and degrade the bivariate normal distribution of the orientation of the fracture to an ordinary 1D normal distribution. This treatment not only simplifies the difficulty of the analysis, but also gives satisfactory results. In this paper, we use this treatment. The angle ${\alpha }_i$ and direction ${\beta }_i$ of the dips of fractures conforming to a random normal distribution are as follows,

$$\{\begin{array}{l}{\alpha }_i={\sigma }_{\alpha }{N}_i+\overline{\alpha }\\ {\beta }_i={\sigma }_{\beta }{N}_i+\overline{\beta }\\ N_i=\sqrt{\frac{12}{n}}\left({\displaystyle \sum _{i=1}^n{R}_i-\frac{n}{2}}\right)\end{array}$$

(1)

where $\overline{\alpha }$ and $\overline{\beta }$ are the means of samples, respectively, ${\sigma }_{\alpha }$ and ${\sigma }_{\beta }$ are their mean squared deviations, respectively, and $R_i$ is a random variable that follows a uniform distribution on the interval [0, 1].

2.1.2. Model of Fracture Dimension

The mean trace length is an important parameter reflecting the dimensions of fractures. It is related to the statistical distribution of the diameter of the fracture in the Baecher disk model. Log-normal and negative exponential distributions are typically used to model the distribution of the trace length of a fracture. Given that field observations have shown that the number of fractures increases significantly with the reduction in the trace length, we conclude that the negative exponential distribution is more consistent with the actual distribution of the trace length of fractures. Warburton [6] established the relationship between the distribution of trace lengths of fractures and the distribution of disk diameters in the Baecher disk model for both measured line sampling and statistical window sampling, respectively. Villaesusa and Brown [26] further proposed an equation for estimating the statistical parameters of the distribution of the diameter of fractures by using the measured trace lengths in the case of a negative exponential distribution,

$$\overline{D}=\frac{\pi }{8}\overline{l}$$

(2)

where $\overline{D}$ is the mean diameter of fractures, and $\overline{l}$ is their mean trace length.

We assume that the expected areas of the rectangular and disk fractures are equal to establish a relationship between the lengths of the edges of the rectangular fracture and the diameter of the disk fracture. Specifically, the two side lengths of the rectangular fracture are set to be independent, and both satisfy a negative exponential distribution with parameter $\overline{E}$. Then, according to the previous assumptions, we can obtain the following:

$$\overline{E}=\frac{\sqrt{\pi }}{2}\overline{D}$$

(3)

By substituting the above equation into the negative exponential distribution, the two side lengths $a_i$ and $b_i$ of the rectangular fracture can be acquired as follows:

$$\{\begin{array}{l}{a}_i=-\frac{\sqrt{\pi }}{2}\overline{D}\ln (1-R_i)\\ b_i=-\frac{\sqrt{\pi }}{2}\overline{D}\ln (1-R_i)\end{array}$$

(4)

where $R_i$ is a random variable that follows a uniform distribution in the interval [0, 1].

2.1.3. Model of Fracture Density

Fracture density is an important index to measure the degree of development of a fracture, and is generally divided into three types, linear density, surface density and bulk density. Oda [27] derived the relationship between bulk density and linear density using a fracture tensor under the assumption that the distribution of orientation of the disk fracture and its diameter are independent of each other. Kulatilake and Wu [28] established a correlation between bulk density and surface density based on the relationship between traces within sample windows and the windows themselves. However, due to the complexity of the calculation process, this method is rarely used.

Wu et al. [29] found that the bulk density of a fracture is relatively constant in a relatively statistically homogeneous region, independent of its orientation, and is related to the linear density as follows:

$${\lambda }_v=\frac{{\lambda }_d}{2\pi {\displaystyle {\int }_0^{\infty }R}{\displaystyle {\int }_R^{\infty }f(r)drdR}}$$

(5)

where ${\lambda }_v$ is the bulk density of the fracture; ${\lambda }_d$ is their linear density, and $f(r)$ is the distribution function of the probability density of the radius of the disk.

Substituting the probability density function of the radius of the disk that satisfies the negative exponential distribution and Equation (2) into Equation (5) yields the following:

$${\lambda }_v=\frac{\overline{{\lambda }_d}}{2\pi {\lambda }^2}=\frac{32\overline{{\lambda }_d}}{{\pi }^3{\overline{l}}^2}$$

(6)

In addition, we assume that the bulk density of the fracture is independent of its shape. Thus, rectangular fractures can be assumed to have approximately the same bulk density as disk fractures in a relatively statistically homogeneous region. Therefore, the number of rectangular fractures in a given study area can be expressed as follows:

$$N=V{\lambda }_v$$

(7)

where $N$ is the number of rectangular fractures contained in the area, $V$ is the volume of the rock mass in the study area, and ${\lambda }_v$ is the bulk density of fractures.

The number of fractures in each group in the study area can be calculated by repeating the above process several times.

2.1.4. Model of Center Point of Fractures

The locations of fractures are frequently described by specifying representative points (e.g., center points) on each fracture, and the spatial distribution of these representative points can be assumed in a statistical process. Currently, the most widely used method for describing the locations of fractures is the Poisson process. It is generally accepted that any small non-overlapping block in rock masses has an equal probability of containing the center point of a fracture. In order to determine the center points of fractures, we divide the study area V into N non-overlapping sub-regions V_i (i = 1, 2, …, N) of equal volume. The center point of the ith fracture is assumed to be randomly and uniformly distributed within V_i. Then, if the range of coordinates of sub-region V_i is $x_i\le x\le x_{i+1}$, $y_i\le y\le y_{i+1}$ and $z_i\le z\le z_{i+1}$, and the coordinates of the center of the fracture are ($x_{ic},y_{ic},z_{ic}$), the following can be obtained:

$$\{\begin{array}{l}{x}_{ic}=x_i+(x_{i+1}-x_i)R_x\\ y_{ic}=y_i+(y_{i+1}-y_i)R_y\\ z_{ic}=z_i+(z_{i+1}-z_i)R_z\end{array}$$

(8)

where $R_x,R_y,R_z$ are random variables that follow a uniform distribution in the interval [0, 1].

2.2. Simulation of 3D DFN of Rock Mass

2.2.1. Monte-Carlo Technique

The Monte-Carlo technique is the inverse process of sampling based on the theory of probability. It involves conducting statistical experiments with random numbers, and using the generated random numbers that conform to some random distribution function as the basis for application. The application of the Monte-Carlo technique to the simulation of 3D DFN is usually based on probability models of the geological parameters of rock fracture to generate a specified number of random numbers to represent the direction of dip, the dip angle, the side length, and coordinates of the center that conform to the form of the corresponding probability distribution. The coordinates of the vertices of each rectangular fracture are then obtained to generate a 3D DFN that conforms to the probability models of the geological parameters of rock fractures.

2.2.2. Determining of the Dimensions of the Simulation

A cube with a certain volume is first set as the domain of examination of the 3D DFN model. The origin of the coordinate system is located at the center of the cube, the x-axis is geographically east, the y-axis is the geographic north, and the z-axis is the height. To avoid the reduction in the density of fractures near the boundary, the generated domain of the fracture network is extended by a factor of 1 in the average trace length of fractures along each of the three dimensions compared with the study domain. After generating the fracture network in the domain, the part of the fracture outside the domain is removed.

2.2.3. Geometric Representation of Rectangular Fractures

Under the assumption that the upper and lower sides of the rectangular fracture are parallel to the horizontal plane and the two side lengths are independently and identically distributed, the dip angle ${\alpha }_i$ and dip direction ${\beta }_i$ of the ith fracture are obtained from Equation (1), the two side lengths $a_i$ and $b_i$ of the ith fracture are obtained from Equation (4), and the coordinates ($x_{io},y_{io},z_{io}$) of its center are obtained from Equation (8). Then, the four vertices of the ith fracture can be obtained by the vector of geometric transformation.

As shown in Figure 1, we let the unit normal vector of the fracture be $\stackrel{\rightharpoonup }{e_i}$, the unit vector of the edge of $\stackrel{\rightharpoonup }{AB}$ be $\stackrel{\rightharpoonup }{e_{AB}}$, the unit vector of the edge $\stackrel{\rightharpoonup }{BC}$ be $\stackrel{\rightharpoonup }{e_{BC}}$, the unit normal vector of the z-axis be $\stackrel{\rightharpoonup }{e_k}\left(0,0,1\right)$, and the coordinates of point A be $(x_{ia},y_{ia},z_{ia})$. The perpendicular vector from the center O of the fracture toward $\stackrel{\rightharpoonup }{AB}$ intersects $\stackrel{\rightharpoonup }{AB}$ at point E, to form a closed vector $\mathsf{\Delta }OEA$. The following relationship obtains in the vector $\mathsf{\Delta }$:

$$\stackrel{\rightharpoonup }{OA}=\stackrel{\rightharpoonup }{OE}+\stackrel{\rightharpoonup }{EA}=(x_{ia}-x_{io},y_{ia}-y_{io},z_{ia}-z_{io})$$

(9)

where

$$\stackrel{\rightharpoonup }{OE}=-\frac{b_i}{2}\stackrel{\rightharpoonup }{e_{BC}}$$

(10)

$$\stackrel{\rightharpoonup }{EA}=-\frac{a_i}{2}\stackrel{\rightharpoonup }{e_{AB}}$$

(11)

It follows from the relationship between the normal vector and the angle and direction of the dip that $\stackrel{\rightharpoonup }{e_i}=(\sin {\alpha }_i\sin {\beta }_i,\sin {\alpha }_i\cos {\beta }_i,\cos {\alpha }_i)$.

From the perpendicular relationship between $\stackrel{\rightharpoonup }{AB}$, $\stackrel{\rightharpoonup }{e_i}$ and $\stackrel{\rightharpoonup }{e_k}$, it follows that

$$\stackrel{\rightharpoonup }{e_{AB}}=\frac{\stackrel{\rightharpoonup }{e_k}\times \stackrel{\rightharpoonup }{e_i}}{|\stackrel{\rightharpoonup }{e_k}\times \stackrel{\rightharpoonup }{e_i}|}=(-\cos {\beta }_i,\sin {\beta }_i,0)$$

(12)

Similarly, it follows from the perpendicular relationship between $\stackrel{\rightharpoonup }{BC}$, $\stackrel{\rightharpoonup }{e_i}$ and $\stackrel{\rightharpoonup }{e_{AB}}$ that

$$\stackrel{\rightharpoonup }{e_{BC}}=\frac{\stackrel{\rightharpoonup }{e_i}\times \stackrel{\rightharpoonup }{e_{AB}}}{|\stackrel{\rightharpoonup }{e_i}\times \stackrel{\rightharpoonup }{e_{AB}}|}=(-\sin {\beta }_i\cos {\alpha }_i,-\cos {\beta }_i\cos {\alpha }_i,\sin {\alpha }_i)$$

(13)

The coordinates of point A$(x_{ia},y_{ia},z_{ia})$ can be determined by substituting Equations (10)–(13) into Equation (9),

$$\{\begin{array}{l}{x}_{ia}=x_{io}+\frac{b_i\sin {\beta }_i\cos {\alpha }_i+a_i\cos {\beta }_i}{2}\\ y_{ia}=y_{io}+\frac{b_i\cos {\beta }_i\cos {\alpha }_i-a_i\sin {\beta }_i}{2}\\ z_{ia}=z_{io}-\frac{b_i\sin {\alpha }_i}{2}\end{array}$$

(14)

Likewise, the coordinates of points $B(x_{ib},{\mathrm{y}}_{ib},z_{ib})$, $C(x_{ic},y_{ic},z_{ic})$ and $D(x_{id},y_{id},z_{id})$ can be calculated, respectively, as follows:

$$\{\begin{array}{l}{x}_{ib}=x_{io}+\frac{b_i\sin {\beta }_i\cos {\alpha }_i-a_i\cos {\beta }_i}{2}\\ y_{ib}=y_{io}+\frac{b_i\cos {\beta }_i\cos {\alpha }_i+a_i\sin {\beta }_i}{2}\\ z_{ib}=z_{io}-\frac{b_i\sin {\alpha }_i}{2}\end{array}$$

(15)

$$\{\begin{array}{l}{x}_{ic}=x_{io}-\frac{b_i\sin {\beta }_i\cos {\alpha }_i+a_i\cos {\beta }_i}{2}\\ y_{ic}=y_{io}-\frac{b_i\cos {\beta }_i\cos {\alpha }_i-a_i\sin {\beta }_i}{2}\\ z_{ic}=z_{io}+\frac{b_i\sin {\alpha }_i}{2}\end{array}$$

(16)

$$\{\begin{array}{l}{x}_{id}=x_{io}-\frac{b_i\sin {\beta }_i\cos {\alpha }_i-a_i\cos {\beta }_i}{2}\\ y_{id}=y_{io}-\frac{b_i\cos {\beta }_i\cos {\alpha }_i+a_i\sin {\beta }_i}{2}\\ z_{id}=z_{io}+\frac{b_i\sin {\alpha }_i}{2}\end{array}$$

(17)

At this point, the probability models of the geological parameters of rock fractures and related parameters obtained from the field observation can be used to generate the corresponding 3D DFN model.

2.2.4. Visualization of the DFN

Although the morphological details of fractures can be acquired using algebraic geometry for the above-generated fracture network, using data as the form of expression seems very abstract and unintuitive, and is not conducive to verifying the validity of the results of the simulation. We use Open-GL to solve the problem of visualizing the fracture network of rock mass, and write a corresponding program to better visualize the 3D DFN. It is also possible to rotate, translate, and scale graphs of the 3D DFN at any angle, and to then save the graph as a specified path.

The specific flow of the program is shown in Figure 2, and it can be seen that the program includes three parts of basic operations. The first part is the initialization of window and data, that is, setting the size, color and depth information of the display window and reading in group data of fractures, and then storing the information containing the geometry and color of fractures in the display list for the preparation of the graphic batch processing. The second part is for the geometry, that is, geometric elements are assembled inside the program, and the anti-aliasing process is used to prevent the graphics from going out of shape and to measure the distance of fractures from the observation point by means of a depth test to determine the geometric relationship between the fractures. The third part is the graphics screen display and image saving operations, that is, the pixel information of the geometry is determined by rasterization, and the rendering of same-nature lighting is achieved by setting the type, color, position and material properties of the light source. On the basis of the above, the graphics are displayed on the screen by suitable model changes, projection methods and viewport settings, and then the graphics are saved in accordance with the requirements.

3. Implementation and Application of DFN Simulation of the Rock Mass

3.1. Implementation

We implemented the 3D DFN model based on the Monte-Carlo technique on the MATLAB platform and the Open GL graphics toolkit in the Visual C environment. The composition and functions of the simulation program were as follows.

(1) Data input module for fractures

The main function of this module is to import data, such as the domain of generation and the scale of the study domain of the DFN model, the state of grouping of fractures, and the parameters of the geological model, into the module for the generation of the DFN in the form of human–machine interaction based on data from field observations of rock masses and experiments on them.

(2) Generation module of the DFN

This module is the core part of the simulated fracture network, and its main function is to generate a certain number of random numbers to represent the dip angle, dip direction, side length, and coordinates of the center of rectangular fractures in the domain of generation according to the parameters of the geological model and information on the area of the simulation domain. The coordinates of the vertices of each rectangular fracture are then obtained according to the formula for the coordinates of the fracture vertex to generate each group of fractures and obtain the data for the fracture network. Finally, the theory of spatial algebraic geometry is used to determine the intersectional tangential relationship between the boundary of the study domain and nearby fractures. The parts of the fracture outside the boundary are removed. The entire process was performed in MATLAB.

(3) Output module of the DFN

This module is mainly concerned with the output of the data generated by the DFN and the visualization of the corresponding data files. The data structure of the DFN generated by MATLAB is converted into a data interface, and the converted data files are imported into a visualization program written in the Visual C environment with the help of an Open GL graphics database to display the 3D DFN. The flow of implementation of the simulation of the DFN, as well as the relationship between the modules, are shown in Figure 3.

3.2. Example of Application

We applied the proposed method to the high rock-slope engineering of the Baiyunebo eastern iron open pit in Inner Mongolia to illustrate its feasibility, and analyzed the results of the simulation.

According to the relevant literature [30], the fractures of this slope can be divided into four groups, and the characteristic statistical values of the parameters of the geological model of each group are shown in

Table 1. Characteristic statistical values of the geological parameters of fractures.

Group Number of Fractures	Dip Direction/(°)		Dip Angle/(°)		Trace Length/(m)	Linear Density/(1/m3)
	Mean	Standard Deviation	Mean	Standard Deviation	Mean
A	232.50	23.38	58.00	16.17	11.90	0.32
B	323.41	17.24	68.92	14.10	12.20	0.35
C	76.40	26.01	75.57	8.14	9.78	0.45
D	171.43	12.30	70.77	6.94	10.03	0.87

. The scales of the study domain in the x, y, and z axes were all taken to be 20 m.

Figure 4 shows the results of the 3D DFN model based on the proposed method using the data in Table 1. Figure 4a reflects the spatial cross-cutting relationship of each group of rectangular fractures, and Figure 4b reflects the spatial cross-cutting relationship between fractures and the boundary of the study domain.

Figure 5 shows the results of the DFN model in comparison with the measured fractures. Figure 5a reflects the distribution of the measured fractures on an outcropping surface of the high rock slope of the Baiyunebo eastern iron open pit in Inner Mongolia, and Figure 5b shows the distribution of fractures obtained by the DFN model on the boundary of the study domain. In general, the DFN simulation ought to identify statistical laws and construct a sample-based fracture network according to random information, where this is not equivalent to a real fracture network and is only statistically similar to it in the mathematical sense. From this perspective, the 3D DFN model based on the Monte-Carlo technique exhibited good statistical similarity i to the actual distribution of fractures at the engineering site on a 2D representative surface in terms of the orientation, scale, and density of fractures. This proves that the proposed method of simulation of the 3D DFN is reasonable and effective.

4. Conclusions

The Monte-Carlo stochastic method, spatial algebraic geometry theory and graphic visualization technique base on Open GL are used to model the three-dimensional DFN of the rock mass with rectangle as the fracture shape. The main discussion and conclusions are as follows:

(1) A probabilistic model of geological parameters based on rectangular fractures was developed here using statistical theory, and the form of the coordinates of vertices of the rectangular fractures was derived;

(2) The process of implementation of the simulation method of the rectangular DFN based on the Monte-Carlo technique was detailed using the MATALB platform with the Open GL graphics toolkit in the Visual C environment;

(3) The simulation of DFN can only achieve “probabilistic statistical simulation”, and it is difficult to achieve “entity simulation” at present. Therefore, the error in simulation of fracture distribution is inevitable, but it can be reduced by improving the accuracy of sampling, correcting sampling bias, and improving simulation algorithms;

(4) Since the DFN is constructed based on the Monte-Carlo stochastic simulation method of probability theory, it is recommended to construct a certain number of DFNs that conform to the corresponding distribution laws in practice to study the fluctuation range of the mechanical properties of rock mass, and to provide a reliable basis for the correct evaluation of the quality of engineering rock mass.

In summary, the results of the proposed model of the DFN were statistically similar to those obtained by on-site measurement in terms of the orientation, scale, and density of fractures in rock slope engineering on a representative 2D surface. This provides the foundation for subsequent research on the quantitative analysis of fractures in rock masses, such as the characteristics of their deformation and strength, the effects of the number of their dimensions, anisotropy, representative elementary volume, the identification of key blocks, the effects of excavation, and mechanisms of fracture of the rock mass. The DFN model in this manuscript is valid for close fractures, and open fractures will be considered in the next stage.