Abstract

In this paper, two basic assumptions are introduced: (1) The number and length distribution of fractures in fractured rock mass are in accordance with the fractal law. (2) Fluid seepage in the fractures satisfies the cubic law. Based on these two assumptions, the fractal model of parallel seepage and radial seepage in fractured rock mass is established, and the seepage tensor of fracture network which reflects the geometric characteristics and fractal characteristics of fracture network under two kinds of seepage is derived. The influence of fracture geometry and fractal characteristics on permeability is analyzed, and the validity and accuracy of the model are verified by comparing the calculated results of the theoretical model and physical model test. The results show that the permeability coefficient of fracture network is a function of the geometric (maximum crack length , fractured horizontal projection length , diameter calculation section porosity , fracture strike , and fracture angle ) and fractal characteristics (fracture network fractal dimension and seepage flow fractal dimension ). With the increase of fractal dimension , the permeability coefficient increases. With the increase of , the permeability coefficient decreases rapidly. And the larger the (), the greater the change of permeability coefficient with .

1. Introduction

A large number of fractures are distributed in the natural rock mass and artificially disturbed rock mass, which provide channels for fluid seepage. To determine the permeability characteristic of fractured rock is important in the field of geology, geotechnical engineering, petrochemical resource exploitation, groundwater resource development and protection, nuclear waste disposal, and storage of carbon dioxide and others [19]. Generally, it is assumed that the fluid flows only in the interconnected fractures, and the permeability characteristics of equivalent fracture network is determined by analyzing the permeability characteristics of the different characteristic fractures [10, 11]. This method is called the discrete fracture network (DFN) method, which has been widely used and developed in recent decades [1216].

The permeability characteristics of the fracture network mainly depend on the fracture characteristics (spatial distribution, density, connectivity, etc.) at the macroscopic level and the fracture characteristics (length, gap width, directivity, roughness, etc.) at the microscopic level. Liu et al. [8] reviewed the current research on the influence of geometrical characteristics of fractured rock mass on the permeability of two-dimensional discrete fracture network and summed up nine parameters which have great influence on the permeability of fracture network. It is the length and distribution of fractures, the width and distribution, fractured surface roughness, dead fracture, fracture cross point, hydraulic gradient, stress condition, anisotropy, and size effect and listed the analytical expressions of the relationship between the relevant fracture parameters and permeability. de Dreuzy and Philippe [17] studied the effect of fracture length and gap width distribution on the infiltration characteristics of two-dimensional random fracture networks. Based on the fractal network statistics and fractal characteristics, a practical method to determine the permeability of fracture network was proposed by Jafari and Babadagli [18, 19]. Sensitivity analysis of the permeability of each parameter was carried out. The results showed that the fracture density and length have the greatest influence on the permeability of fracture network. The permeability of fracture network with exponential distribution was studied by Rossen and Gu [20]. It is found that the single large scale fracture plays a major role in the permeability of the fracture network. The permeability of the study area increases with the increase of the study size. Under normal circumstances, it is very difficult to quantitatively describe the fracture characteristics of the fracture network due to its complexity [21, 22]. Studies have shown that natural fractures have fractal properties that have the same effect on the permeability of fracture networks as well as fracture geometries [18, 2326].

A series of studies have been made on the influence of geometric characteristics and fractal features of fractures on permeability. A two-dimensional random fracture network model considering the characteristics of fracture network (density, length, gap width, directionality, and connectivity) was established. The fracture network statistics and the relationship between fractal characteristics and fracture network permeability were obtained through the multiple regression analysis and artificial neural network processing. An estimation method was proposed to estimate the permeability of multilayer complex fracture network, in order to obtain a good estimate result, which was recommended for comprehensive utilization of drill core data (1D), bedrock exposed data (2D), and drilling test data (3D). Another key factor that affects the permeability of fracture network is the connectivity between fractures. Based on the fractal geometry and percolation theory, the two-dimensional fracture network connectivity was defined, and the relationship between fracture network permeability and fractal dimension and dimensionless fracture density and percolation threshold was given and analyzed. It is pointed out that the estimated value obtained by fractal dimension of fractured network was more accurate than that of fractal dimension of fracture intersections, fracture connectivity rate, and intersection fractal dimension in the X-direction and Y-direction.

In recent years, many scholars have established the fractal model of pore and fractured media based on the fractal characteristics of porosity and fractures and have made great achievements in studying fluid seepage, solute transport, and heat conduction [2732]. Zheng and Yu [33] deduced the seepage characteristics of the gas in the porous media model composed of porous rock and fractal tree fractures. It is shown that the pore fractal dimension, pore bending degree, porosity, the ratio of maximum pore diameter to fractal tree fracture network length, diameter ratio, bifurcation angle, and bifurcation level are the key factors that have significant effects on gas permeability. Yun et al. [34] analyzed the plane radial seepage and parallel seepage of Newtonian fluid in porous media and deduced the permeability coefficient and flow rate and velocity expression in the two seepage cases. The fractal theory and the Monte Carlo method were used to establish the probability model of radial percolation in porous media by Xu et al. [35]. The results showed that the effective radial permeability coefficient decreased drastically with the radial distance increase and the porosity and pore fractal dimension of the radial seepage interface had a significant effect on the effective permeability coefficient. Based on the fractal geometry theory and the laminar cube law of fluid in fractures, a fractal model of seepage in fractured rock mass was deduced by Miao et al. [36]. The theoretical model showed that the permeability coefficient of fractured rock mass was a function of fracture fractal dimension, porosity, fracture density, maximum fracture length, gap width, fracture direction, and inclination angle. Miao et al. [37] assumed that the fractures in the fractured rock mass were randomly distributed, the fluid flow in the fracture was in accordance with the law of cube, the pores in the rock mass connected into tortuous channels, the fracture length and pore diameter distribution had fractal characteristics, and a fractal model of seepage flow in a two-hole media model was established. The Monte Carlo method was used to generate the fractured network, and a fractal flow model was established to reflect the geometrical characteristics of fractured rock mass. It is revealed that the bending degree, the gap width and the random number that reflects fractal regularity of fracture length distribution, had a significant influence on the permeability of fracture network [38]. The seepage fractal model of discrete fracture network was established by using the fractal dimension of fracture geometrical distribution and the fractal dimension of flow line which reflects the surface roughness of fractures. By simulating and calculating the equivalent permeability coefficient of fluid flowing through different geometric characteristics of discrete fractures, it was found that when the fractal dimension of the fractured network is less than 1.5, the seepage of the fracture network is mainly controlled by the small fracture with the length less than the width of the fracture network. With the increase of fractal dimension, the effect of long fractures on fracture network seepage is increasing [7, 8].

The permeability of fractured rock mass is influenced by the heterogeneity of fracture direction, inclination, and length distribution, which lead to the fact that the permeability of the fractured rock mass is directional, and many scholars have carried on fruitful research in this area [9, 10, 12, 13, 39, 40]. Based on the previous research, this paper is aimed at establishing a parallel seepage and radial seepage fractal model of fractured rock mass considering the effect of fracture surface roughness in Section 2 and deducing the seepage tensor of the fracture network which reflects the geometrical and fractal characteristics of fracture network in Section 3. The influence of fracture geometry and fractal characteristics on the permeability is analyzed in Section 4. The validity and accuracy of the model are verified by comparing the theoretical model and the physical model water injection test in Section 5.

2. Fractal Characteristics of Fractured Network of Rock Mass

Studies have shown that the cumulative surface area distribution of the Earth’s surface islands is subject to power distribution [36, 41], i.e., where is the total number of islands with area greater than constant and is the fractal dimension representing the area distribution of the island. On this basis, is used to represent the area of the largest island [42]. Equation (2) is given:

The relationship between the equivalent gap width and the crack length can be expressed [36, 4345]: where is the equivalent gap width; is the proportional coefficient, which is related to the mechanical properties of the surrounding rock and in the range of 0.001-0.1 [37, 44]; is the fracture length; and is a constant that reflects the fracture characteristics, which ranges from 0.5 to 2.0 [37, 44, 45]. When , it indicates that the equivalent gap width and the length of the fracture are linearly distributed. The fractured network has self-similarity and fractal characteristics [37, 38, 43, 45]. Equation (3) can be written as follows:

Many studies have shown that the fracture length distribution satisfies the fractal law [15, 36, 4549]. Therefore, equation (2), which describes the distribution law of island area, is used to describe the distribution of fractures area in fractured media: where and are the maximum gap width and maximum length, respectively, and and are the fracture width and length, respectively. Equation (4) is substituted into equation (5), and equation (6) is given as follows: where is the fractal dimension of the fracture length, for the two-dimensional problem and for the three-dimensional problem . By substituting in equation (6) with , the total number of fractures in the fracture network can be expressed as follows:

In general, the number of fractures in the fracture network is large, so equation (6) can be approximated as a continuous differential equation and equation (8) can be obtained by solving differential in equation (6):

Equation (8) indicates the number of fractures in [], and the negative sign indicates that the number of fractures decreases with increasing fracture length. Equation (8) divided by equation (7) is equal to the following equation: where is the probability density function, which is satisfied by equation (10) by the definition of the probability function:

Therefore, equation (11) is given as follows:

In general, in equation (11) is a necessary condition for the fractal network to show fractal characteristics. Many researchers [22, 38] used as a threshold for the fractal model considering the two-dimensional fracture network seepage. This criterion is also used in this paper.

The relationship between porosity and fractal dimension can be expressed [37, 50, 51]: where is the fractal dimension of the fracture length; is the porosity of the fracture network; and are the minimum and maximum values of the fracture length; is the European dimension, in which for the two-dimensional problem and for the three-dimensional problem .

According to the definition of porosity, the following equation is available: where is the porosity of the fracture network, is the area of the cross section, and is the pore area.

By solving simultaneously equations (12) and (14), the following equation becomes available:

By substituting equation (13) into equation (15), the following equation becomes available:

3. Fractal Model of Permeability of Fractured Rock Mass

3.1. Formula Derivation

In a three-dimensional space, the direction of the fracture is determined by the direction and inclination, as shown in Figure 1. The -axis and the -axis are the two coordinate directions of the horizontal plane, which is in the same direction as the geodetic coordinate; that is, the -axis is and the -axis is ; is the angle between the trend line and the -axis, θ is the angle between the fracture plane and the horizontal plane (inclination), is the fracture length, is the fracture width, is the seepage length of the fluid in the fracture considering the fracture roughness, is the linear length in the direction of fracture seepage, is the tangential component of the hydraulic gradient in the fracture, and is the normal component of the hydraulic gradient in the fracture. is the linear length of projected to the -axis direction.

Usually, there are a large number of cracks in the fractured rock mass, and it is impossible to determine the direction and inclination of the fractures one by one. Studies have shown that the directionality of many fractures in a given area is not exactly the same but can usually show a tendency [36, 52]. For example, the results of 1878 fractures in the study showed an average inclination of 70° and an average orientation of N-S [53]. Furthermore, the fracture network studied at present is only caused by mining. Compared with the natural fracture networks affected by complex geological factors, the fracture network studied here is more regular. Therefore, we take the trend and dip angle in the calculation model for a certain range of statistical average.

It is assumed that the flow of fluid in the fracture can be described by the cubic law [39, 44, 54]: where is a single fracture flow, is the dynamic viscosity coefficient, is the fracture width, is the fracture length, is the length of the seepage of the fluid in the fracture considering the fracture roughness, and ΔP is the pressure difference across the fracture.

Because of the rough surface of the fracture, the fluid flow path in the fracture is a curve, which leads to the extension of the flow path and the decrease of the effective flow capacity, as shown in Figure 1. The relationship between and can be expressed as follows: where is the fractal dimension of the flow line of the seepage flow. reflects the nonlinearity of the streamline, and the streamline is a straight line when , .

By substituting equations (4) and (18) into equation (17), the following equation becomes available:

In a three-dimensional space, the hydraulic gradient is divided into normal and tangential components along the fracture, as shown in Figure 1. Only the tangential hydraulic gradient produces the seepage flow, that is,

By substituting into equation (20), the following equation becomes available:

It is assumed that the angle between the normal direction of the fracture surface and the coordinate axes is , , and , respectively. The normal direction vector of the fracture surface can be expressed as

At the same time, the hydraulic gradient vector is decomposed along the fracture plane normal and tangential components [55]:

Taking equation (24) into consideration,

The hydraulic gradient is decomposed in the direction of the coordinate axis; equation (25) is available.

According to theorem of vector operation, equation (26) is available.

Equations (22), (23), (24), (25), and (26) are simultaneously solved; equation (27) is available.

Equation (27) is written in a matrix form:

The total seepage flow can be obtained by integrating seepage flow of single fracture along the fracture length on the calculated cross section.

By substituting equations (8), (21), and (28) into equation (29), equation (30) becomes available:

Considering , , , , equation (30) is simplified and rewritten as a matrix:

Figure 1 shows that the relationship between the direction cosine of fracture normal plane and the direction and inclination of the geodetic coordinate is as follows:

By substituting equation (32) into equation (31), equation (33) can be obtained as follows:

The permeability coefficient can be expressed as

By substituting equations (16) and (33) into equation (34), equation (35) is obtained as follows:

3.2. Radial Seepage

In a three-dimensional space, the radial flow diagram is shown in Figure 2. Map coordinates and parameter definitions are consistent with Figure 1. The - and -axes are two coordinate directions, which are consistent with the direction of geodetic coordinates. It is assumed the -axis is , the -axis is . is the angle between the strike line and the -axis, and is the angle between the fracture plane and the horizontal plane (inclination). is the fracture length, is the fracture width, is the seepage length of the fluid in the fracture considering the fracture roughness, is taken as the straight line length along fracture seepage direction, is radius of the hole, is the tangential component of the hydraulic gradient in the fracture, and is the normal component of the hydraulic gradient in the fracture.

It is assumed that the flow of fluid in the fracture can be described by the cubic law [39, 44, 54]. Equation (17) can be rewritten as follows under the condition of radial flow: where is a single fracture flow, is the dynamic viscous coefficient, is the fracture width, is the fracture length, and is the water pressure. Due to the rugged surface of the fracture, the fluid seepage path in the fracture is a curve, resulting in a decrease in the flow path and an effective overcurrent capacity, as shown in Figure 2.

The relation between and can be rewritten according to equation (18) as follows:

Equation (38) can be obtained by substituting equations (4) and (37) into equation (36) as follows:

Similarly, the hydraulic gradient vector will be decomposed along the fracture plane in normal and tangential components [55]. Only the hydraulic gradient along the tangential direction of the fracture will produce seepage flow [12], i.e.,

On the inner side of the radial hole wall, the seepage flow of single fracture is integrated along the fracture length, and the total seepage flow can be obtained.

Equation (16) is substituted into equation (40) and is rewritten in a matrix form:

Combining equations (41) and (34) and being replaced by , equation (42) is obtained as follows:

The permeability coefficient can be expressed as follows at the inner side of the radial hole wall:

From equations (42) and (43), the dimensionless permeability coefficient is obtained as follows:

Equation (44) shows that the dimensionless permeability coefficient is closely related to and . It decreases with the increase of and decreases with the increase of , as shown in Figure 3.

3.3. Parameter Impact Analysis

This section mainly analyzes the influence of the parameters of equations (35) and (42) on the permeability coefficient , and the values of each parameter are shown in Table 1.

Figure 4 shows the relationship between the fractal dimension and the permeability coefficient in the case of different , and the other parameters are , , and . As shown in Figure 4, the permeability coefficient increases with the increase of fractal dimension . And with the increase of , the permeability coefficient decreases. The results show that the increase of the fracture surface roughness leads to the increase of the flow path and the decrease of flow velocity.

Figure 5 shows the relationship between and permeability coefficient in the case of different , and the other parameters are , , and . As shown in Figure 5, with the increase of , the permeability coefficient decreases rapidly. The overall permeability coefficient increases with the increase of , and the larger the (), the greater the change rate of with , which means that when , the weakening effect of the fracture surface roughness on the permeability of the fracture network increases rapidly.

As can be seen from Figures 4 and 5, and have a significant effect on the permeability of the fracture network, but they have different physical meanings and influences. The larger the value of , the more complex the fractured network, the higher the fracture density [18], and the greater the permeability coefficient. reflects the flexural degree of the fluid seepage path in the fracture due to fracture surface roughness. With the increase of , the seepage path is lengthened and the effective overcurrent capacity of the fracture is reduced, which leads to the decrease of the permeability coefficient of the fracture network.

Figure 6 shows the relationship between porosity and permeability coefficient in the case of different , and the other parameters are , , and . As shown in the figure, the permeability coefficient increases with the increase of porosity . In general, it can be seen from Equation (12) and Figure 7 that, assuming , is positively correlated with , and thus the increase of porosity and increase have identity, and all reflect the increase of seepage cross section. It also shows in Figure 6 that the overall permeability coefficient decreases with the increase of , and the larger the (), the smaller the change rate of with , which indicates that when , the roughness of the fracture surface will greatly offset the effect of increased permeability of the fracture network due to the increase of the seepage area.

Figure 8 shows the relationship between fractal dimension and permeability coefficient of the fracture network under the condition of the change of the fracture dip angle , and the other parameters are , , and . As shown in the figure, the permeability coefficient increases with the increase of porosity , which is the same as that of Figure 4. At the same time, as the fracture dip angle increases, the permeability coefficient decreases.

Figure 9 shows the relationship between the fracture dip angle and the permeability coefficient , and the other parameters are , , , and . As shown in Figure 9, the permeability coefficient decreases with the increasing of fracture dip angle from 0 to . Figure 10 shows the relationship between the fracture strike and the permeability coefficient , and the other parameters are , , , and . As shown in Figure 10, the permeability coefficient increases as the dip angle of fracture increases from to . The variation of the permeability coefficient component with the fracture strike and dip angle reflects the nonuniformity and seepage directivity of fractured rock mass.

Figure 11 shows the relationship between the maximum fracture length and permeability coefficient under condition of different . The other parameters are , , , and . As shown in Figure 11, the permeability coefficient increases with the increase of . At the same time, with the increase of , the overall permeability coefficient will increase, and with the increase of (), the change rate of with get sharply greater, indicating that when , with the increase of fractal dimension, the influence of long fractures on fracture network seepage is increasing, which is consistent with the results of Li et al. [22].

Figure 12 shows the relationship between and the permeability coefficient , and the other parameters are , , , and . As shown in Figure 12, the permeability coefficient increases with the increase of . Considering equation (4), this also reflects the relationship between the fracture width and the permeability coefficient to a certain extent.

Figure 13 shows the relationship between and permeability coefficient , and the other parameters are , , , , and . As shown in Figure 13, the permeability coefficient decreases gradually with the increase of . This is because the increase in is the growth of the seepage path, which means that the seepage resistance increases and the hydraulic gradient decreases.

4. Experimental Verification

4.1. Physical Model Test

Figure 14 shows a large physical model (partial) profile. According to the generalized stratigraphic parameters shown in Table 2, the similar materials arranged according to a similarity ratio are poured up and down, and the bottom pouring height of no. 2 coal seam and no. 5 coal seam was 2.25 m and 0.5 m, respectively. The thickness of the layer is 0.05 meters, the net length of working face is 2.1 meters, and the reserved pillar width is 0.3 meters. The displacement and stress changes of the overlying strata and the expansion and evolution of the fractures are studied under the conditions of different burial depth and mining speed. Figure 15 shows the distribution of the fractures of the overlying strata after the mining of no. 5 coal seam.

4.2. Water Injection Test in the Physical Model

In order to obtain the permeability parameters of the rock mass in different fracture zones and verify the fractal model of the fractured rock mass, the water injection test in the physical model was carried out after the physical model excavation was completed. In order to compare with the evolution of fracture expansion, the working face is divided into seven zones along the excavation direction, and the width of each zone is 30 cm. In order to reduce the interaction of the boundary effect and the water injection holes, the water injection holes were arranged in plum blossom in the middle of the partition, as shown in Figure 16.

The location of water injection test sites is carried out at the top of the model in accordance with Figure 16 when constructing, and then the hydraulic drill is set up to conduct drilling construction, as shown in Figures 17 and 18. The water injection test is carried out by the top-down stratified water injection method to obtain the permeability parameters of different plane positions and different elevation rocks in the mining face. Meanwhile, the corresponding fracture distribution images are obtained by drilling imaging an analyzer, and the fracture images are processed digitally to obtain the fracture expansion images, as shown in Figure 19. The fracture expansion images are analyzed to extract the fracture geometric characteristic parameters, and thus the fractal characteristic parameters, and , are calculated by the box counting method [51, 5659].

5. Analysis and Discussion

According to the field test conditions, the continuous head water injection test is used to determine the seepage per unit time in the case of constant water head and then calculate the permeability coefficient of the simulated strata according to equation (44). where is the permeability coefficient of the soil layer (cm/s), is the stable water flow (cm3/s), is the test water head (cm), is the length of the water injection section (cm), and is the water injection test hole radius (cm).

Fracture development images are analyzed in the previous section; then the fractal characteristic parameters are extracted, and the relevant parameters are substituted into equation (42). The calculated results are compared with the values calculated by equation (43), as shown in Table 3. And the fractal characteristics are calculated, as shown in Figure 20. The relevant parameters are substituted into equation (41); then the calculated results are compared with the measured values for the water injection test, as shown in Figure 21. The calculated values of and are compared with the measured values of water injection, as shown in Figure 22. From Figures 2022, it shows that the theoretical calculation is in good agreement with the experimental results, which proves the validity of the theoretical model.

6. Conclusions

In this paper, a fractal model for characterizing hydraulic properties of fractured rock mass under mining influence considering fracture surface roughness is established, and the infiltration tensor of fracture network is deduced under two kinds of seepage conditions. The effects of fractal geometry and fractal characteristics on the permeability are analyzed. The validity and accuracy of the model are verified by comparing the theoretical model with the physical model water injection test. The results show the following: (1)The permeability coefficient of the fracture network is a function of the geometric of fracture network. In general, the maximum fracture length , calculation section porosity and fracture network fractal dimension reflect the related indicators of the calculation section, such as the fracture length, discharge area, and density, which play a role in promoting permeability. The horizontal projection length and the fractal dimension of the seepage flow line reflect the obstacle of the surface roughness and the increase of the seepage path to the seepage fluid, which inhibit the permeability(2)With the increase of fractal dimension , the permeability coefficient increases. While with the increase of , the permeability coefficient decreases rapidly, and the larger the (), the greater the change of permeability coefficient with , which indicates that when , the weakening effect of fracture surface roughness on the permeability of fracture network increases rapidly(3) and have a greater impact on the permeability of the fracture network, but their physical meanings and influences are different. The larger the value, the more complex the fractured network in the cross section, and the higher the fracture density, thus the higher the permeability coefficient. reflects the rupture degree of the fluid seepage path in the fracture due to the rough surface of the fracture. With the increase of , the seepage path is lengthened and the effective overcurrent capacity of the fracture is reduced, which leads to the decrease of the permeability coefficient of the fracture network(4)With the increase of porosity , the permeability coefficient increases. While with the increase of , the permeability coefficient decreases. The higher the (), the smaller the change of permeability coefficient with the change of , which shows that when , the fracture surface roughness will be greatly offset due to the increase of water seepage area leading to increasing effect on permeability coefficient(5)With the increase of , the permeability coefficient increases. Simultaneously, with the increase of , the overall permeability coefficient will increase, and with the increase of (), the rate of permeability coefficient with change sharply increases, which shows that when , with the increase of fractal dimension, the influence of long fractures on fracture network seepage is increasing

Data Availability

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

Conflicts of Interest

The authors declare that they have no conflicts of interest.

Acknowledgments

The authors would like to acknowledge the support from the National Key Research and Development Plan (Grant No. 2016YFC0501104), the National Natural Science Foundation of China (Grant Nos. 51522903, 51479094, and 41772246), and the Key Research and Development Plan of Ningxia Hui Autonomous Region (Grant No. 2018BCG01003).