Numerical Simulation of the Nonlinear Flow Properties in Self-Affine Aperture-Based Fractures

In order to study the effect of fracture geometry on the nonlinear flow properties in aperture-based fractures, a fractal model based on the self-affinity is proposed to characterize the three-dimensional geometry of rough-walled fractures. By solving the N–S (Navier–Stokes) equation directly, the relationships between the Forchheimer-flow characteristics, fractal dimension, and standard deviation of the aperture have been obtained. *e Forchheimer equation is validated to describe the nonlinear relationship between flow rate and pressure gradient. For lower flow rate, the influence of the fractal dimension almost can be ignored, but the linear coefficient increases and the hydraulic aperture decreases with increasing standard deviation of the aperture, respectively. For larger flow rate, the nonlinear coefficient increases with the growth of the standard deviation of the aperture and fractal dimension.*us, an empirical relationship between the nonlinear coefficient, fractal dimension, and standard deviation of aperture is proposed. In addition, the critical Reynolds number decreases with the increase of the standard deviation of the aperture and the fractal dimension, and the numerical results are generally consistent with the experimental data.


Introduction
Fractures are prevalent among the natural rock masses under geological action. Compared with the intact rock, the permeability of fractures is much larger, and fractures become the dominant channels for fluid flow [1,2]. Fracturedominated flow plays an important role in many practical situations such as seepage control for fractured media [3,4], hazardous waste isolation [5], slope engineering [6,7], underground tunneling [8,9], and many geological disasters.
Many efforts have been taken to investigate the hydraulic properties of fractures. Snow [10] has conceptualized the rough-walled fracture as the smooth parallel plate model, and the famous cubic law was derived. However, the surface roughness and aperture distribution of natural fractures are random and irregular [11,12]. Consequently, there is no absolute smooth fracture in natural rock masses and it is difficult to satisfy the establishment conditions of cubic law, which may lead to a large deviation in the prediction of fracture permeability. Based on the laboratory tests and numerical analysis, many researchers [13][14][15][16][17][18][19][20][21] found that the relationship between flow rate and pressure gradient in rough-walled fractures is nonlinear deviating from the cubic law.
To investigate the relationship between the nonlinear flow properties and the roughness of the fracture surface, Zhang and Tian [22] carried out numerical simulation of a single fracture with different roughness and tortuousness, and the results showed a certain deviation from the modified cubic law resulted by the flow tortuosity. Chen et al. [16] studied the relationship between the hydraulic aperture and Forchheimer equation's nonlinear coefficient by conducting flow tests under different confining pressures on twelve groups of granite cracks with different roughness. Yin et al. [19] analyzed the effects of the shearing process under a series of normal loads on nonlinear flow behavior in 3D rough-walled fractures with different roughness by the shear-flow test. e abovementioned investigations on the impact of roughness on the nonlinear flow properties all use the JRC (Joint Roughness Coefficient) to represent the fracture surface roughness. However, the value of the JRC is obtained based on observation experiences; thus, the roughness of fracture surfaces cannot be accurately and quantitatively represented.
Some other influence factors of rough fractures on the nonlinear flow properties were also considered. For example, Xia et al. [23] found that different contact conditions and the root-mean square height of fractures had an evident impact on apparent transmissivity based on laboratory observations. Tsang [24] explored the impact of tortuous of flow path on flow behavior through experiments and found that the smaller the aperture distribution, the greater the influence of tortuosity. Xiong et al. [17] designed a saturated seepage test of fracture under low flow rate and evaluated the role of roughness and aperture in affecting the nonlinear flow properties. It is not difficult to find that there are many factors affecting the nonlinear flow properties of rough fractures. However, the effect of the standard deviation of the fracture aperture on the nonlinear flow properties is rarely discussed.
e main objective of this study is to study the nonlinear flow behavior in self-affine aperture-based rock fractures based on the fractal theory and Navier-Stokes equations. Based on the fractal theory, the rough-walled fracture is reconstructed by fractional Brownian motion, and the surface roughness and aperture distribution can be characterized by fractal dimension. For rough-walled fractures with different surface roughness and the standard deviation of the aperture, the nonlinear relationship between flow rate and pressure gradient can be well described by the Forchheimer equation. e mathematical relationship between the nonlinear coefficient, fractal dimensions, and the standard deviation of the aperture is also quantified.

Geometrical Model of Rough Fractures
Previous studies [18,[25][26][27] have shown that the fracture roughness can be described by self-affinity, which can be simulated using fractional Brownian motion (fBm). e height of rough fracture surfaces is described by a continuous and single random function Z(x). e stationary increment [Z(x)−Z(x + Δ)] over the distance Δ follows a Gaussian distribution with mean zero and variance δ 2 . e self-affinity relating to fBm obeys the following expressions: where <·> is the mathematical expectation, x donates the coordinate component, and H represents the Hurst exponent varying from 0 to 1 and associated with the fractal dimension D by D � 3-H. λ is a constant, δ 2 λΔ and δ 2 Δ are the variance corresponding to the height variation of fracture surface with the distance of λΔ and Δ, respectively, and δ is the standard deviation of the aperture.
To construct the geometrical model of rough-walled fractures, in present study, the successive random addition method (SRAM) [28][29][30] is used to generate the fracture surfaces. e generated square area is shown in Figure 1, and the specific steps are as follows: (1) In the given single square region, the initial random values of the four corners, which are labeled as number 1, satisfy the Gaussian distribution N(0, δ 2 0 ) (2) e center point of the square and midpoints of each side are marked by number 2, and their height are the average value of the four corners initial height and the average value of two points of each side, respectively; at the same time, the random values from N(0, δ 2 1 ) should be added into all points, in which (2) (3) Step (2) is repeated for each newly generated square, and the random values from N(0, δ 2 n ) should be added to all heights as well until 2 n × 2 n squares are created at iteration (4) No new square is inserted again, but we keep adding random values from N(0, δ 2 j ) to all points, in which where j � n + 1, n + 2, . . ., NM; NM is large enough, so δ NM /δ 0 is negligible e upper and lower surfaces of rough fracture are generated using the SRAM. e height of low surface is Z 1 (x,y), and the height of the upper surface can be calculated by the following formula: where u is the mean aperture between the upper and lower surfaces. e distribution function of the fracture aperture can be expressed as To demonstrate the reliability of the algorithm used to characterize the fracture surface geometry, Figure 2 shows the fracture aperture distribution diagram of four groups with different roughness. e side length of square L � 40 mm, δ � 0. 15  Advances in Civil Engineering respectively. It can be seen that, with the increase of D, the maximum values of the aperture increase and the minimum values of the aperture decrease; the aperture distribution is more scattered, the variation of adjacent apertures is more fluctuated, and the degree of roughness is greater.

Theory of Flow Properties in Fractures
3.1. Cubic Law. Based on the smooth parallel plate model, the flow behaviour in a single fracture should obey the cubic law [10]: where Q is the volumetric flow rate per unit time, b h is the hydraulic aperture of the fracture, w is the fracture width, μ is the dynamic viscosity, and ∇p � (p in − p out )/l is the pressure gradient between the inlet and outlet.

Forchheimer Equation.
For larger flow rate, the flow behaviour of rough fractures is nonlinear due to the larger inertial effect. us, instead of equation (7), the Forchheimer equation [14,15,17,19] is commonly applied to describe the nonlinear flow properties: where A and B are the linear coefficient and nonlinear coefficient, respectively; k is the fracture permeability; β is the non-Darcy flow inertial coefficient which depends on the geometric characteristics of fracture surfaces [16,31]; and ρ is the fluid density. When the flow rate is small, the inertia force is much less than the viscosity force; in other words, the quadratic term (BQ 2 ) is much less than the linear term (AQ), and equation (8) can be reduced to the cubic law.

Distinguishing the Flow Regime.
In order to quantify the nonlinear flow behaviour of fractures, the Reynolds number Re is calculated as follows: where v is the average velocity of the fracture inlet. e Reynolds number Re represents the ratio of inertial force to viscous force in fracture flow, and the inertial effect is stronger with larger Re; thus, it is easier to enter the nonlinear flow regime.
Simultaneously, the non-Darcy effect factor E is defined according to the Forchheimer equation: e physical meaning of E is the proportion of pressure gradient caused by nonlinear flow in the total pressure gradient, and it is a dimensionless coefficient ranging from 0 to 1 representing the degree of nonlinear flow. e greater the E is, the stronger the nonlinear effect is. At present, E � 0.1 [15,23,32], and combined with equations (10) and (11), the critical Reynolds number Re c for the transition from linear to nonlinear flow of fracture flow yields e smaller the Re c is, the more significant the inertia effect is and the easier it is to be the nonlinear flow state.

Governing Equation.
For the incompressible Newtonian fluid with a constant viscosity coefficient flowing in the rough fractures, the fluid motion is governed by the N-S equation and the continuity equation: where u and ∇ are the velocity vector and Hamiltonian operator, respectively. In this study, the fluid density is 997.1 (kg/m 3 ), and the dynamic viscosity is 0.894 ×10 −3 (Pa s) for water at 25°C. Since the flow rate is small and generally does not exceed the boundary of laminar flow, the laminar interface in finite element software COMSOL Multiphysics is selected for the solution [33].

Numerical
Procedure. As shown in Figure 3, at first, three-dimensional rough fracture surfaces are generated based on the SRAM, and then, the solid fracture model is constructed by the Boolean operation. Next, the solid model is meshing to determine the microstructure of the computing model. Finally, the software COMSOL is employed to obtain a series of data of flow rate and pressure gradient which is fitted by the Forchheimer equation.
Considering the solution accuracy, calculation cost, and errors caused by different element sizes, the tetrahedral element size is set as 0.25 mm, and the number of grid elements ranges from 400 thousands to 600 thousands. e size of the fracture model is 40 mm × 40 mm, and the mean aperture u is 0.8 mm; the geometric parameters including the standard deviation of the aperture and fractal dimension are shown in Table 1. e left side of the fracture model is specified as inflow boundary, the range of Q is 3.586 ×10 −8 ∼3.228 × 10 −6 ·m 3 /s, the corresponding Re ranges from 1 to 90 according to equation (10). e right side of the fracture model is defined as outflow boundary and the pressure is set as zero, and the rest of the boundaries are impervious.

Relationship between Flow Rate and Pressure Gradient.
e relationships between flow rate and pressure gradient of five groups are shown in Figure 4. e fitting linear coefficient A and nonlinear coefficient B of the Forchheimer equation are presented in Table 1, and the coefficients of determination R 2 are both greater than 0.99, which indicates the nonlinear relationship between the flow rate and pressure gradient in self-affine aperture-based fractures can be well described by the Forchheimer equation.
In Figure 4, with the increment of flowrate, the deviation between Forchheimer curves and cubic law increases drastically.
is deviation is also strengthened with the increment of D and δ, which indicates that flow resistance will be greater for the rougher surface. Zeng and Grigg [32] suggested that the Forchheimer-flow properties were mainly caused by inertial force; in other words, the inertial term ρ(u · ∇u) in equation (13) was the main factor for the nonlinear behavior.
In order to analyze the primary factor of nonlinear flow behavior, Figure 5 compared the five velocity sections along the x direction when D � 2.5 and Q � 5.3796 × 10 −7 ·m 3 /s of two fracture models with different δ (0.09 mm and 0.21 mm). It can be seen that the velocity distribution is more scattered in the fracture with more heterogeneous aperture distribution, and the local velocity becomes relatively larger. As a result, the inertia effect of flow is enhanced, and the local energy dissipation is increased, which leads to the fluid entering nonlinear state. (9), the linear coefficient A is negatively correlated with the permeability of rough fractures. As shown in Table 1, A increases with the δ increase, indicating that the permeability will be lower when the distribution of the aperture is more discrete and irrelevant for a larger δ. e same phenomenon can also be found in Figure 6, the hydraulic aperture b h decreases apparently with the δ increases, and it is always less than the mean aperture 0.8 mm; therefore, the permeability of rough fractures is less than that of smooth fractures. Moreover, with the increment of δ, the flow paths become more tortuous and the permeability becomes smaller.

Analysis of Permeability of Rough Fractures. Based on equation
However, D has no obvious effect on A except for a slight fluctuation as shown in Table 1. In other words, when the flow is small, the fractal dimension has little effect on the flow capacity of rough cracks e section diagrams of streamline distribution are presented in Figure 7, for δ � 0.21 mm and Q � 2.152 × 10 −6 ·m 3 /s with D being 2.1, 2.3, and 2.5, y � 20 mm, and x ranging from 25 mm to 40 mm. For same δ, the streamline distribution is similar while the roughness of the fracture surfaces is greater for larger D, and there are some blank areas at the rough bulge. is is because the fluid flow in the fracture tends to the region with lower resistance and will bypass the high resistance area to form a dominant channel. erefore, the effective flow space is substantially equal and the flow capacity does not change significantly. When the flow rate continues to increase, eddy flow will appear in these blank areas, resulting in accelerated energy consumption and the permeability of fractures being reduced [34].

Analysis of Forchheimer-Flow Characteristics.
e B and β represent the degree of evolution of Forchheimer-flow properties, and the nonlinear flow is stronger with a greater value of B and β. In Table 1, the values of B increase with larger δ and D, which indicates that the degree of flow nonlinearity is more drastic with the rougher surfaces and more heterogeneous aperture distribution.
According to equation (9), B is proportional to β but inversely proportional to the square of b h , so it is difficult for B to fully reflect the impact of fracture morphology on flow properties of rough fractures. Instead, Figure 8 shows the relationship between δ, D, and β. It is observed that β increases linearly with δ and D. erefore, the two-parameter expression is proposed: where η, m, and n are empirical parameters. On the basis of nonlinear Levenberg-Marquardt algorithm, the fitting parameters η, m, and n are 517.5, −967, and 49.39, respectively, and the coefficient of determination R 2 is 0.9879, indicating that this expression could well explain the effect of morphology parameters on Forchheimer-flow properties. Combining with equation (12), the empirical model of nonlinear coefficient B is obtained: It should be noted that when the value of b h decreases with the increment of δ, the growth of B is more significant.    when Re > Re c , the flow behavior was controlled by the inertia effect, appearing Forchheimer-flow properties. As shown in Figure 9, for the same condition, Re c has a decreasing trend with the increment of δ and D, indicating that the fluid flow can more easily develop into nonlinear Forchheimer-flow resulting from the rougher surfaces and more heterogeneous aperture distribution.
Previous studies [13,14,[35][36][37][38] on Re c of flow in rock fractures are listed in the Table 2. In this study, for rough fractures of u � 0.8 mm, δ being 0.09∼0.21 mm and D ranging from 2.1 to 2.5, the calculated values of Re c are between 9.01 and 14.47, which are almost consistent with the previous researches.

Conclusions
In order to estimate the relationships between the nonlinear flow properties and self-affine fracture geometry, a systematic procedure with respect to the nonlinear flow analysis for the three-dimensional rough-walled fractures model is proposed based on the fractal characteristics, in which the roughness of fracture surfaces is characterized by the fractal dimension and the fractal model is generated by the SRAM. e numerical simulation of nonlinear flow analysis in selfaffine aperture-based fractures is presented based on the fractal theory and N-S equation, in which the spatial variation of fracture geometry including roughness and aperture distribution is characterized by the fractal dimension and N-S equation is solved by the software COMSOL. e main conclusions are as follows: (1) e larger the fractal dimension, the greater the fracture surface roughness and the lower the correlation of the local aperture. e larger the standard deviation of the aperture, the greater the fluctuation of fracture surface and the more heterogeneous the aperture distribution. (2) e nonlinear relationship between flow rate and pressure gradient in self-affine aperture-based fractures can be well described by the Forchheimer equation. e flow rate, standard deviation of the aperture, and fractal dimension can improve the deviation of pressure gradient from the cubic law. (3) With the increase of the standard deviation of the aperture, the linear coefficient of the Forchheimer equation is larger, and the hydraulic aperture becomes smaller. When the flow rate is small, fractal dimension has little effect on the permeability of fracture. (4) e nonlinear coefficient of the Forchheimer equation increases with the increment of the standard deviation of the aperture and fractal dimension, and the empirical expression between the nonlinear coefficient, the standard deviation of aperture, and the fractal dimension is proposed. e standard deviation of the aperture has greater impact on the nonlinear flow behavior than the fractal dimension. e critical Reynolds number decreases with the increase of standard deviation of the aperture and fractal dimension, and its measured range is 9.01∼14.47, which is almost consistent with the previous results.
Data Availability e data used in the present study can be made available upon request from the authors.

Conflicts of Interest
e authors declare that they have no conflicts of interest.