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

Authors: Xin Zhou; Jianlong Sheng; Ruili Lu; Zuyang Ye; Wang Luo
Journal: Advances in Civil Engineering
 (2021)

---

## Abstract

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. The 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. Thus, 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.


---

## Body

## 1. 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]. Fracture-dominated 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–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. The 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.The 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. The mathematical relationship between the nonlinear coefficient, fractal dimensions, and the standard deviation of the aperture is also quantified.


## 2. Geometrical Model of Rough Fractures

Previous studies [18, 25–27] have shown that the fracture roughness can be described by self-affinity, which can be simulated using fractional Brownian motion (fBm). The height of rough fracture surfaces is described by a continuous and single random function Z(x). The stationary increment [Z(x)−Z(x + Δ)] over the distance Δ follows a Gaussian distribution with mean zero and variance δ2. The self-affinity relating to fBm obeys the following expressions:(1)Zx−Zx+λΔ=0,δλΔ2=λ2HδΔ2⇒δλΔ=λHδΔ,δλΔ2=Zx−Zx+λΔ2,δΔ2=Zx−Zx+Δ2,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–30] is used to generate the fracture surfaces. The 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 distributionN0,δ02(2)
The 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 fromN0,δ12 should be added into all points, in which(2)δ12=δ0222H1−22H−2.(3)
Step (2) is repeated for each newly generated square, and the random values fromN0,δn2 should be added to all heights as well until 2n×2n squares are created at iteration(3)δn2=δn−1222H=δ0222Hn1−22H−2.(4)
No new square is inserted again, but we keep adding random values fromN0,δj2 to all points, in which(4)δj2=δj−1222H=δ0222Hj1−22H−2,where j = n + 1, n + 2, …, NM; NM is large enough, so δNM/δ0 is negligibleFigure 1
Illustration of the two-dimensional SRAM.
(a)(b)(c)The upper and lower surfaces of rough fracture are generated using the SRAM. The height of low surface isZ1(x,y), and the height of the upper surface can be calculated by the following formula:(5)Z2x,y=Z1x+Δx,y+Δy+u,where u is the mean aperture between the upper and lower surfaces. The distribution function of the fracture aperture can be expressed as(6)bx,y=Z2x,y−Z1x,y,ifZ2x,y>Z1x,y,0,else.To demonstrate the reliability of the algorithm used to characterize the fracture surface geometry, Figure2 shows the fracture aperture distribution diagram of four groups with different roughness. The side length of square L = 40 mm, δ = 0.15 mm, and fractal dimensions are 2.1, 2.2, 2.3, and 2.4, 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.Figure 2
Comparison of the morphology of fracture models forδ = 0.15 mm (red represents large aperture; blue represents small aperture; unit: mm). (a) D = 2.1. (b) D = 2.2. (c) D = 2.3. (d) D = 2.4.
(a)(b)(c)(d)


## 3. 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]:(7)Q=−wbh312μ∇p,where Q is the volumetric flow rate per unit time, bh is the hydraulic aperture of the fracture, w is the fracture width, μ is the dynamic viscosity, and ∇p=pin−pout/l is the pressure gradient between the inlet and outlet.


### 3.2. Forchheimer Equation

For larger flow rate, the flow behaviour of rough fractures is nonlinear due to the larger inertial effect. Thus, instead of equation (7), the Forchheimer equation [14, 15, 17, 19] is commonly applied to describe the nonlinear flow properties:(8)−∇P=AQ+BQ2,(9)A=μkAh=12μwbh3,B=βρw2bh2,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 (BQ2) is much less than the linear term (AQ), and equation (8) can be reduced to the cubic law.


### 3.3. Distinguishing the Flow Regime

In order to quantify the nonlinear flow behaviour of fractures, the Reynolds numberRe is calculated as follows:(10)Re=ρvbhμ=ρQμw,where v is the average velocity of the fracture inlet. The 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 factorE is defined according to the Forchheimer equation:(11)E=BQ2AQ+BQ2.The physical meaning ofE 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. The 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 Rec for the transition from linear to nonlinear flow of fracture flow yields(12)Rec=Aρ9Bμw.The smaller theRec is, the more significant the inertia effect is and the easier it is to be the nonlinear flow state.


## 3.1. Cubic Law

Based on the smooth parallel plate model, the flow behaviour in a single fracture should obey the cubic law [10]:(7)Q=−wbh312μ∇p,where Q is the volumetric flow rate per unit time, bh is the hydraulic aperture of the fracture, w is the fracture width, μ is the dynamic viscosity, and ∇p=pin−pout/l is the pressure gradient between the inlet and outlet.


## 3.2. Forchheimer Equation

For larger flow rate, the flow behaviour of rough fractures is nonlinear due to the larger inertial effect. Thus, instead of equation (7), the Forchheimer equation [14, 15, 17, 19] is commonly applied to describe the nonlinear flow properties:(8)−∇P=AQ+BQ2,(9)A=μkAh=12μwbh3,B=βρw2bh2,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 (BQ2) is much less than the linear term (AQ), and equation (8) can be reduced to the cubic law.


## 3.3. Distinguishing the Flow Regime

In order to quantify the nonlinear flow behaviour of fractures, the Reynolds numberRe is calculated as follows:(10)Re=ρvbhμ=ρQμw,where v is the average velocity of the fracture inlet. The 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 factorE is defined according to the Forchheimer equation:(11)E=BQ2AQ+BQ2.The physical meaning ofE 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. The 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 Rec for the transition from linear to nonlinear flow of fracture flow yields(12)Rec=Aρ9Bμw.The smaller theRec is, the more significant the inertia effect is and the easier it is to be the nonlinear flow state.


## 4. Numerical Simulation

### 4.1. 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:(13)ρu⋅∇u−μ∇2u=−∇P,(14)∇u=0,where u and ∇ are the velocity vector and Hamiltonian operator, respectively. In this study, the fluid density is 997.1 (kg/m3), 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].


### 4.2. Numerical Procedure

As shown in Figure3, 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.Figure 3
Modelling process of nonlinear flow through 3D rough fracture.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. The size of the fracture model is 40 mm × 40 mm, and the mean apertureu is 0.8 mm; the geometric parameters including the standard deviation of the aperture and fractal dimension are shown in Table 1. The left side of the fracture model is specified as inflow boundary, the range of Q is 3.586 ×10−8∼3.228 × 10−6·m3/s, the corresponding Re ranges from 1 to 90 according to equation (10). The 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.Table 1
The values for linear coefficientA and nonlinear coefficient B.
Standard deviation of aperture (mm)D = 2.1D = 2.2D = 2.3D = 2.4D = 2.5ABABABABABδ = 0.095.51111.18025.45361.27795.51181.36425.46481.41795.51471.5342δ = 0.125.68751.22895.64281.32945.61681.41945.69851.54785.66471.6495δ = 0.155.96251.29425.97631.41315.93161.55365.89041.64915.88621.8116δ = 0.186.27101.37496.16221.50756.25581.67586.09191.76666.13352.0007δ = 0.216.53531.48176.49331.59966.44951.77246.45311.97746.45462.2197Notes:A (108 kg·s−1·m−5); B (1014 kg·m−8).


## 4.1. 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:(13)ρu⋅∇u−μ∇2u=−∇P,(14)∇u=0,where u and ∇ are the velocity vector and Hamiltonian operator, respectively. In this study, the fluid density is 997.1 (kg/m3), 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].


## 4.2. Numerical Procedure

As shown in Figure3, 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.Figure 3
Modelling process of nonlinear flow through 3D rough fracture.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. The size of the fracture model is 40 mm × 40 mm, and the mean apertureu is 0.8 mm; the geometric parameters including the standard deviation of the aperture and fractal dimension are shown in Table 1. The left side of the fracture model is specified as inflow boundary, the range of Q is 3.586 ×10−8∼3.228 × 10−6·m3/s, the corresponding Re ranges from 1 to 90 according to equation (10). The 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.Table 1
The values for linear coefficientA and nonlinear coefficient B.
Standard deviation of aperture (mm)D = 2.1D = 2.2D = 2.3D = 2.4D = 2.5ABABABABABδ = 0.095.51111.18025.45361.27795.51181.36425.46481.41795.51471.5342δ = 0.125.68751.22895.64281.32945.61681.41945.69851.54785.66471.6495δ = 0.155.96251.29425.97631.41315.93161.55365.89041.64915.88621.8116δ = 0.186.27101.37496.16221.50756.25581.67586.09191.76666.13352.0007δ = 0.216.53531.48176.49331.59966.44951.77246.45311.97746.45462.2197Notes:A (108 kg·s−1·m−5); B (1014 kg·m−8).


## 5. Results and Discussion

### 5.1. Relationship between Flow Rate and Pressure Gradient

The relationships between flow rate and pressure gradient of five groups are shown in Figure4. The fitting linear coefficient A and nonlinear coefficient B of the Forchheimer equation are presented in Table 1, and the coefficients of determination R2 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.Figure 4
Relationship between flowrate and pressure gradient. (a)δ = 0.09 mm. (b) δ = 0.12 mm. (c) δ = 0.15 mm. (d) δ = 0.18 mm. (e) δ = 0.21 mm.
(a)(b)(c)(d)(e)In Figure4, with the increment of flowrate, the deviation between Forchheimer curves and cubic law increases drastically. This 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, Figure5 compared the five velocity sections along the x direction when D = 2.5 and Q = 5.3796 × 10−7·m3/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.Figure 5
Velocity distribution forD = 2.5, Q = 5.3796 × 10−7·m3/s. (a) δ = 0.09 mm. (b) δ = 0.21 mm.
(a)(b)


### 5.2. Analysis of Permeability of Rough Fractures

Based on equation (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 δ. The same phenomenon can also be found in Figure 6, the hydraulic aperture bh 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.Figure 6
Relationship between the hydraulic aperture and geometric characteristics.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 The section diagrams of streamline distribution are presented in Figure 7, for δ = 0.21 mm and Q = 2.152 × 10−6·m3/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. This 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. Therefore, 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].Figure 7
Diagram of streamline distribution forδ = 0.21 mm, Q = 2.152 × 10−6·m3/s (y = 20 mm, flow direction x = 25∼40 mm). (a) D = 2.1. (b) D = 2.3. (c) D = 2.5.
(a)(b)(c)


### 5.3. Analysis of Forchheimer-Flow Characteristics

TheB 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 bh, 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. Therefore, the two-parameter expression is proposed:(15)β=ηδD+mδ+nD,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 R2 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:(16)B=ρw2bh2517.5δD−967δ+49.39D.Figure 8
The relationship between the non-Darcy inertia coefficient and geometric characteristics.It should be noted that when the value ofbh decreases with the increment of δ, the growth of B is more significant. Therefore, the effect of δ on the nonlinear flow behavior of the rough fractures is greater than that of D.


### 5.4. Critical Reynolds Number

Rec can be used to evaluate the flow regime of fractures, representing the Reynolds number of the fluid flow going from the linear state to nonlinear state. When Re < Rec, the flow behavior satisfied the cubic law on account of the inertia effect being very small and the viscous effect having a dominant role of fluid flow in fractures; when Re > Rec, the flow behavior was controlled by the inertia effect, appearing Forchheimer-flow properties. As shown in Figure 9, for the same condition, Rec 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.Figure 9
Relationship between the critical Reynolds number and geometric characteristics.Previous studies [13, 14, 35–38] on Rec 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 Rec are between 9.01 and 14.47, which are almost consistent with the previous researches.Table 2
Literature summaries about critical Reynolds number of flow in rock fractures.
AuthorsMaterialsAperture (μm)RecStudy methodKonzuk and KueperTensile fractures of limestoneMean = 3812.8∼14.3Laboratory testJavadi et al.Tensile fractures of granite—0.001∼25Laboratory testZhang and NemcikTensile fractures of sandstone6.14∼18.953.5∼24.8Laboratory testZimmerman et al.3D fracturesMean = 1491∼10.20Numerical simulationChen et al.Sedimentary and intrusive rocks—25∼66Field testQian et al.Artificial parallel plate model1000∼2500245∼759Laboratory testThis paper3D rough fracturesMean = 8009.01∼14.47Numerical simulation


## 5.1. Relationship between Flow Rate and Pressure Gradient

The relationships between flow rate and pressure gradient of five groups are shown in Figure4. The fitting linear coefficient A and nonlinear coefficient B of the Forchheimer equation are presented in Table 1, and the coefficients of determination R2 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.Figure 4
Relationship between flowrate and pressure gradient. (a)δ = 0.09 mm. (b) δ = 0.12 mm. (c) δ = 0.15 mm. (d) δ = 0.18 mm. (e) δ = 0.21 mm.
(a)(b)(c)(d)(e)In Figure4, with the increment of flowrate, the deviation between Forchheimer curves and cubic law increases drastically. This 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, Figure5 compared the five velocity sections along the x direction when D = 2.5 and Q = 5.3796 × 10−7·m3/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.Figure 5
Velocity distribution forD = 2.5, Q = 5.3796 × 10−7·m3/s. (a) δ = 0.09 mm. (b) δ = 0.21 mm.
(a)(b)


## 5.2. Analysis of Permeability of Rough Fractures

Based on equation (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 δ. The same phenomenon can also be found in Figure 6, the hydraulic aperture bh 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.Figure 6
Relationship between the hydraulic aperture and geometric characteristics.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 The section diagrams of streamline distribution are presented in Figure 7, for δ = 0.21 mm and Q = 2.152 × 10−6·m3/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. This 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. Therefore, 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].Figure 7
Diagram of streamline distribution forδ = 0.21 mm, Q = 2.152 × 10−6·m3/s (y = 20 mm, flow direction x = 25∼40 mm). (a) D = 2.1. (b) D = 2.3. (c) D = 2.5.
(a)(b)(c)


## 5.3. Analysis of Forchheimer-Flow Characteristics

TheB 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 bh, 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. Therefore, the two-parameter expression is proposed:(15)β=ηδD+mδ+nD,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 R2 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:(16)B=ρw2bh2517.5δD−967δ+49.39D.Figure 8
The relationship between the non-Darcy inertia coefficient and geometric characteristics.It should be noted that when the value ofbh decreases with the increment of δ, the growth of B is more significant. Therefore, the effect of δ on the nonlinear flow behavior of the rough fractures is greater than that of D.


## 5.4. Critical Reynolds Number

Rec can be used to evaluate the flow regime of fractures, representing the Reynolds number of the fluid flow going from the linear state to nonlinear state. When Re < Rec, the flow behavior satisfied the cubic law on account of the inertia effect being very small and the viscous effect having a dominant role of fluid flow in fractures; when Re > Rec, the flow behavior was controlled by the inertia effect, appearing Forchheimer-flow properties. As shown in Figure 9, for the same condition, Rec 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.Figure 9
Relationship between the critical Reynolds number and geometric characteristics.Previous studies [13, 14, 35–38] on Rec 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 Rec are between 9.01 and 14.47, which are almost consistent with the previous researches.Table 2
Literature summaries about critical Reynolds number of flow in rock fractures.
AuthorsMaterialsAperture (μm)RecStudy methodKonzuk and KueperTensile fractures of limestoneMean = 3812.8∼14.3Laboratory testJavadi et al.Tensile fractures of granite—0.001∼25Laboratory testZhang and NemcikTensile fractures of sandstone6.14∼18.953.5∼24.8Laboratory testZimmerman et al.3D fracturesMean = 1491∼10.20Numerical simulationChen et al.Sedimentary and intrusive rocks—25∼66Field testQian et al.Artificial parallel plate model1000∼2500245∼759Laboratory testThis paper3D rough fracturesMean = 8009.01∼14.47Numerical simulation


## 6. 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. The numerical simulation of nonlinear flow analysis in self-affine 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. The main conclusions are as follows:(1)
The larger the fractal dimension, the greater the fracture surface roughness and the lower the correlation of the local aperture. The larger the standard deviation of the aperture, the greater the fluctuation of fracture surface and the more heterogeneous the aperture distribution.(2)
The nonlinear relationship between flow rate and pressure gradient in self-affine aperture-based fractures can be well described by the Forchheimer equation. The 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)
The 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. The standard deviation of the aperture has greater impact on the nonlinear flow behavior than the fractal dimension. The 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.



---

*Source: 6687878-2021-06-07.xml*