# Numerical Calculation of Fracture Seepage in Rough Rock and Analysis of Local Pressure Drop

Authors: Gang Chen; Ling Ma; Hongsheng Gong; Fengqiang Luo
Journal: Geofluids
 (2021)

---

## Abstract

The seepage performance of a rock mass mainly depends on the rock fractures developed in it. Numerical calculation method is a common method to study the permeability properties of fractures. Seepage in rock fractures is affected by various factors such as fracture aperture, roughness, and filling, among which aperture and roughness are the two most widely influenced factors. The Navier-Stokes (NS) equation can be solved directly for the seepage flow in rock fractures with good accuracy, but there are problems of large computational volume and slow solution speed. In this paper, the fracture aperture space data is substituted into the local cubic law as an aperture function to form a numerical calculation method for seepage in rough rock fractures, namely, the aperture function method (AFM). Comparing with the physical seepage experiments of rock fractures, the calculation results of AFM will produce a small amount of error under the low Reynolds number condition, but it can greatly improve the calculation efficiency. The high efficiency of calculation makes it possible to apply AFM to the calculation of large-scale 3D rough fracture network models. The pressure drop of fluid in the fracture has viscous pressure drop (VPD) and local pressure drop (LPD). VPD can be calculated using the AFM. After analyzing the results of solving the NS equation for fracture seepage, it is concluded that the LPD includes the pressure drop caused by area crowding in the recirculation zone (RZ), kinetic energy loss in the RZ, kinetic energy loss in the vortices, and other reasons.


---

## Body

## 1. Introduction

Rock fractures are widely distributed in subsurface rock masses, and for fluid transport, the permeability of rock fractures is significantly greater than that of the rock matrix [1]. The hydraulic properties of rock fractures play an important role in assessing the performance of subsurface engineering, such as geothermal energy development, enhanced recovery, nuclear waste disposal, and groundwater pollution control [2–6]. Fluid flow in rock fractures can be calculated using the Navier-Stokes (NS) equation. The NS equation is used to solve the fracture flow with high accuracy, but the fracture aperture dimension is much smaller compared to the dimension in the extension direction, which makes the modeling difficult and the calculation workload is large. Zhou et al. [7] used the NS equation to solve a real rock fracture seepage model with a size of 150mm×120mm, and the number of model elements exceeded 106 to ensure the computational accuracy.Due to the complexity of the NS equation solution, a concise fracture seepage equation, the cubic law, can be derived by neglecting the nonlinear terms in fluid flow. The same conclusion was obtained for the seepage experiments using two smooth flat plates [8]. The natural rock fracture surface is rough, which restricts the flow of fluid in the fracture, and using the cubic law to calculate the fracture flow would overestimate the permeability of the fracture [9, 10]. Through theoretical derivation and experimental research, some scholars proposed that under the condition of low Reynolds number (Re<<1), the effect of inertial flow of fluid can be neglected, and the local cubic law (Reynolds equation) can be used to calculate the fracture seepage flow, and its calculation results do not deviate much [11–14]. Zimmerman and Bodvarsson [12] found by single fracture seepage experiment that when Re is in the range of 1 ~ 10, there is a weak inertia effect zone within the fluid in the fracture, and the fracture seepage flow gradually deviates from the linear relationship with the hydraulic gradient but becomes a nonlinear relationship. By correcting the geometric mean aperture, arithmetic mean aperture, surface roughness factor, and curvature factor, the local cubic law can still be used to evaluate the fluid flow in the fracture [11, 14, 15]. Zhu et al. [16] concluded by theoretical derivation and experimental study that the cubic law still has good applicability in the laminar flow range (Re≤2300) when the flow of the fracture is long enough.Javadi et al. [17] and Liu et al. [18] found that there is a deviation between the calculated results of linear and nonlinear flow in the fracture, but the difference is not too big at low Reynolds number conditions. Javadi et al. [17] considered that there are two types of pressure drop during fluid flow in the fracture: one is viscous pressure drop (VPD) due to fluid viscosity; the other is the local pressure drop (LPD) due to the sudden change of aperture. For laminar flow in a smooth flat plate, the VPD can be calculated using the cubic law. For rough fractures, corrections for fracture wall roughness and fracture curvature are required.In this paper, the aperture function method (AFM) based on the local cubic law is introduced, which utilizes the fracture aperture function as the parameter of the numerical model, where one fracture wall surface is used as the geometric model, and the rough fracture seepage is calculated based on the local cubic law. The results of the AFM calculations are compared with the physical experiments and the results of the numerical model of rough fracture seepage based on the solution of the NS equation. Multiple pressure boundary conditions are applied to the inlet of the fracture models, and the rough fracture seepage flow is solved based on the NS equation, and the formation of the recirculation zone (RZ) in the fracture flow and the reasons for the influence on the seepage flow are analyzed in detail. The rotational motion of fluid flow is called vortex, and the presence of vortices has a significant effect on the motion of the fluid due to the vortices. This paper presents the mathematical definition and identification of vortices in fluids and identifies the vortices and their distribution in rough fracture flows. A comparative analysis of vortices in fracture flow and RZ is carried out, and the effect of vortices on fracture seepage flow is studied.


## 2. Methodology

### 2.1. Aperture Function Method (AFM)

#### 2.1.1. Aperture Interpolation Method

Rough rock fracture surface morphology is complex and can be obtained by using high precision laser scanning technology [19–21], and the obtained data can well reflect the real fracture surface morphology. For the obtained fracture surface data, let the fracture surface undulation height be z and use x and y to represent the fracture plane coordinates; then, the fracture surface geometry can be represented by the function z=gx,y in the Cartesian coordinate system. The 3D scanned fracture surface is processed to form a list of values that can be expressed by the function z=gx,y. For subsequent calculations, the x-direction is set to move in the shear direction. The upper fracture surface function is gUx,y, and the lower fracture surface function is gLx,y. To calculate the aperture between the two fracture walls, the following equation can be used [22–24].
(1)ex,y=gUx+us,y−gLx,y+uvgUx+us,y>gLx,y0gUx+us,y≤gLx,y,where ex,y is the upper and lower fracture wall aperture functions, us is the shear direction displacement, and uv is the normal displacement.Whenex,y does not correspond exactly to the spatial coordinates x,y of gUx,y (or gLx,y), it is necessary to use the interpolation method to find the aperture value at the desired spatial location. In Figure 1(a), the gi,j point on the geometric slit surface corresponds to the ei,j′ point in the aperture, but there is no aperture value at this point, so the interpolation method is used to obtain the aperture value at the ei,j′ point.Figure 1
Schematic diagram of fracture aperture interpolation method: (a) schematic diagram of the relationship between the grid points of the fracture wall and the grid points of the aperture; (b) solve the value of the aperture at a certain position (ex,y).The fracture aperture is a function of the spatial coordinatesx,y, so each of the four points can be expressed as x1,y1,ex1,y1, x1,y2,ex2,y1, x2,y1,ex2,y1, and x2,y2,ex2,y2 (Figure 1(b)). If we want to get the aperture value ex,y at a certain position in the middle of the four points, we can use linear interpolation to estimate it:First, interpolatex to find ex,y1 and ex,y2:
(2)ex,y1=ex1,y1⋅x2−xx2−x1+ex2,y1⋅x−x1x2−x1,ex,y2=ex1,y2⋅x2−xx2−x1+ex2,y2⋅x−x1x2−x1.Then, use the same idea to findex,y:
(3)ex,y=ex1,y1⋅y2−yy2−y1+ex,y2⋅y−y1y2−y1.The above equation gives the value of the aperture at any position within the spatial distribution of the aperture. Each pointgi,jx,y on the fracture wall grid can also be found corresponding to the desired aperture value.When discretizing the fracture numerical model, each fracture wall grid point corresponds to a tiny geometric grid cell. Each grid cell can find the corresponding fracture aperture data, and within this cell, the aperture value is a given value; then, the geometric grid cell and the aperture value can form a tiny virtual flat flow model. For each tiny virtual flat plate model, the fracture flow calculation can be performed using the cubic law, and the overall view is an application of the local cubic law. Since the fracture surface and fracture aperture values of each cell are a list of values that vary with space, a rough fracture model with two walls is formed. A triangular grid Figure2(a) or a rectangular grid Figure 2(b) can be used when the fracture data are discrete.Figure 2
Discretizing the rock fracture to the volumetric elements: (a) triangular discrete grid; (b) quadrilateral discrete grid.
(a)(b)


#### 2.1.2. From the NS to the Local Cubic Law Equations

The fluid flow in the fracture of the rock is very slow, and neglecting the effect of inertial forces, the Navier-Stokes (NS) equation can be derived from the [12, 25, 26] equation for the fracture flow in the laminar flow state, i.e., the Reynolds equation, also called the local cubic law.The two-dimensional space is expressed as follows:(4)∂∂xρge312μ∂h∂x+∂∂yρge312μ∂h∂y=0.Substituting the aperture functionex,y into Equation (4) yields a clearer expression of the equation [27].
(5)∂∂xρge3x,y12μ∂h∂x+∂∂yρge3x,y12μ∂h∂y=0,where h is the hydraulic head, μ is the dynamic viscosity, ρ is the fluid density, and g is the gravitational acceleration, and ex,y is the aperture function. Figure 2 shows the meaning of Equation (5) graphically by dividing the rough fracture surface into a series of connected small parallel plates.In this paper, the calculation method of using the aperture function to characterize the aperture value and solving the rough fracture flow based on the local cubic law is referred to as the aperture function method (AFM).


### 2.2. Recirculation Zone in Fracture Seepage Flow

Lee et al. [25], Zhou et al. [28], and Zhou et al. [7] found through physical and numerical experiments that in the steady laminar flow state, the fluid in the fracture will generate a recirculation zone (RZ) when flowing through the abrupt change of aperture. The fluid in the RZ also consumes part of the flow energy, and because the fluid in the RZ occurs back on its own, it does not produce effective flow at the outlet of the fracture. The size of the RZ is nonlinearly related to the fluid flow rate and the fracture roughness procedure, and its magnitude cannot be calculated using a simple equation. In addition, Lee et al. [29] found experimentally that under specific conditions, fluid can also produce slip phenomenon at the fracture wall, which also affects the overall seepage flow rate of the fracture.In order to further verify whether a fluid RZ is generated in the fracture model in this paper and to analyze whether the RZ effect is a linear or nonlinear problem, with other geometric and aperture conditions unchanged, a variety of pressure boundary conditions are set, and two methods, NS solution and AFM solution, are used to calculate and compare, respectively, (Table1). The JRC1 fracture model inlets were applied with pressures ranging from 10 Pa to 6000 Pa (−dh/dl=0.01~5.1) for a total of 52 submodels. JRC3 fracture model inlets were applied with pressures ranging from 10 Pa to 890 Pa (−dh/dl=0.01~0.76) for a total of 36 submodels. JRC6 fracture model inlets were applied 25 Pa~650 Pa pressure (−dh/dl=0.02~0.55), respectively, for a total of 26 sub-models. JRC9 fracture model inlets were applied with 10 Pa~1000 Pa pressure (-dh/dl =0.01~0.85), respectively, for a total of 34 submodels. The outlet pressure of all fracture models is 0 Pa. The fluid density in the model is 998.2 kg/m3, and the dynamic viscosity coefficient is 0.0010093 Pa.s. The calculation software was carried out using COMSOL Multiphysics finite element method.Table 1
Basic information of JRC1, JRC3, JRC6, and JRC9 numerical models.
SolutionNumber of elementsElement number ratioMesh size/mHydraulic gradients (-dh/dl)Number of submodelsJRC1NS9017918.763.97×10−8~1.33×10−50.01~5.152AFM1029303.36×10−5~3.11×10−434JRC3NS7433774.93.58×10−7~3.1×10−50.01~0.7636AFM1517901.82×10−5~2.78×10−436JRC6NS2657723.419.83×10−7~8.52×10−50.02~0.5526AFM778842.54×10−5~3.9×10−426JRC9NS1789003.31.80×10−7~6.04×10−50.01~0.8534AFM542027.44×10−5~6.88×10−434


### 2.3. Identification of Vortices in Seepage Flow

The phenomenon of vortices is prevalent in the fluid flow process. Intuitively, vortices represent the rotational motion of the fluid; that is, where there is fluid rotation, there are vortices, and conversely, where there are vortices, there is fluid rotation [30]. In order to give a quantitative mathematical definition of vortex, a parameter Ω is introduced, representing the ratio of the size of the rotating part of the vortex to the size of the total vortex, calculated as [31]
(6)Ω=BF2AF2+BF2+ε,where
(7)AF2=∂u∂x2+∂v∂y2+∂w∂z2+12∂u∂y2+12∂u∂z2+12∂v∂x2+12∂v∂z2+12∂w∂x2+12∂w∂y2+∂u∂y∂v∂x+∂u∂z∂w∂x+∂v∂z∂w∂y,BF2=12∂u∂y2+12∂u∂z2+12∂v∂x2+12∂v∂z2++12∂w∂x2+12∂w∂y2−∂u∂y∂v∂x−∂u∂z∂w∂x+∂v∂z∂w∂y,where F represents the Frobenius parametrization of the matrix, A is the symmetric part of the fluid velocity gradient tensor, and B is the antisymmetric tensor (vorticity).According to its physical meaning, it is obvious thatΩ takes values in the range 0≤Ω≤1, which can be interpreted as the concentration of vorticity, and more specifically, Ω represents the rigidity of fluid motion, and when Ω=1, it represents the fluid doing rigid body rotation. When Ω>0.5 represents that the antisymmetric tensor B dominates over the symmetric tensor A. Therefore, Ω slightly larger than 0.5 can be used as the criterion for vortex identification. Liu et al. [31] proposes to use Ω=0.52 to determine the boundary of vortex.


## 2.1. Aperture Function Method (AFM)

### 2.1.1. Aperture Interpolation Method

Rough rock fracture surface morphology is complex and can be obtained by using high precision laser scanning technology [19–21], and the obtained data can well reflect the real fracture surface morphology. For the obtained fracture surface data, let the fracture surface undulation height be z and use x and y to represent the fracture plane coordinates; then, the fracture surface geometry can be represented by the function z=gx,y in the Cartesian coordinate system. The 3D scanned fracture surface is processed to form a list of values that can be expressed by the function z=gx,y. For subsequent calculations, the x-direction is set to move in the shear direction. The upper fracture surface function is gUx,y, and the lower fracture surface function is gLx,y. To calculate the aperture between the two fracture walls, the following equation can be used [22–24].
(1)ex,y=gUx+us,y−gLx,y+uvgUx+us,y>gLx,y0gUx+us,y≤gLx,y,where ex,y is the upper and lower fracture wall aperture functions, us is the shear direction displacement, and uv is the normal displacement.Whenex,y does not correspond exactly to the spatial coordinates x,y of gUx,y (or gLx,y), it is necessary to use the interpolation method to find the aperture value at the desired spatial location. In Figure 1(a), the gi,j point on the geometric slit surface corresponds to the ei,j′ point in the aperture, but there is no aperture value at this point, so the interpolation method is used to obtain the aperture value at the ei,j′ point.Figure 1
Schematic diagram of fracture aperture interpolation method: (a) schematic diagram of the relationship between the grid points of the fracture wall and the grid points of the aperture; (b) solve the value of the aperture at a certain position (ex,y).The fracture aperture is a function of the spatial coordinatesx,y, so each of the four points can be expressed as x1,y1,ex1,y1, x1,y2,ex2,y1, x2,y1,ex2,y1, and x2,y2,ex2,y2 (Figure 1(b)). If we want to get the aperture value ex,y at a certain position in the middle of the four points, we can use linear interpolation to estimate it:First, interpolatex to find ex,y1 and ex,y2:
(2)ex,y1=ex1,y1⋅x2−xx2−x1+ex2,y1⋅x−x1x2−x1,ex,y2=ex1,y2⋅x2−xx2−x1+ex2,y2⋅x−x1x2−x1.Then, use the same idea to findex,y:
(3)ex,y=ex1,y1⋅y2−yy2−y1+ex,y2⋅y−y1y2−y1.The above equation gives the value of the aperture at any position within the spatial distribution of the aperture. Each pointgi,jx,y on the fracture wall grid can also be found corresponding to the desired aperture value.When discretizing the fracture numerical model, each fracture wall grid point corresponds to a tiny geometric grid cell. Each grid cell can find the corresponding fracture aperture data, and within this cell, the aperture value is a given value; then, the geometric grid cell and the aperture value can form a tiny virtual flat flow model. For each tiny virtual flat plate model, the fracture flow calculation can be performed using the cubic law, and the overall view is an application of the local cubic law. Since the fracture surface and fracture aperture values of each cell are a list of values that vary with space, a rough fracture model with two walls is formed. A triangular grid Figure2(a) or a rectangular grid Figure 2(b) can be used when the fracture data are discrete.Figure 2
Discretizing the rock fracture to the volumetric elements: (a) triangular discrete grid; (b) quadrilateral discrete grid.
(a)(b)


### 2.1.2. From the NS to the Local Cubic Law Equations

The fluid flow in the fracture of the rock is very slow, and neglecting the effect of inertial forces, the Navier-Stokes (NS) equation can be derived from the [12, 25, 26] equation for the fracture flow in the laminar flow state, i.e., the Reynolds equation, also called the local cubic law.The two-dimensional space is expressed as follows:(4)∂∂xρge312μ∂h∂x+∂∂yρge312μ∂h∂y=0.Substituting the aperture functionex,y into Equation (4) yields a clearer expression of the equation [27].
(5)∂∂xρge3x,y12μ∂h∂x+∂∂yρge3x,y12μ∂h∂y=0,where h is the hydraulic head, μ is the dynamic viscosity, ρ is the fluid density, and g is the gravitational acceleration, and ex,y is the aperture function. Figure 2 shows the meaning of Equation (5) graphically by dividing the rough fracture surface into a series of connected small parallel plates.In this paper, the calculation method of using the aperture function to characterize the aperture value and solving the rough fracture flow based on the local cubic law is referred to as the aperture function method (AFM).


## 2.1.1. Aperture Interpolation Method

Rough rock fracture surface morphology is complex and can be obtained by using high precision laser scanning technology [19–21], and the obtained data can well reflect the real fracture surface morphology. For the obtained fracture surface data, let the fracture surface undulation height be z and use x and y to represent the fracture plane coordinates; then, the fracture surface geometry can be represented by the function z=gx,y in the Cartesian coordinate system. The 3D scanned fracture surface is processed to form a list of values that can be expressed by the function z=gx,y. For subsequent calculations, the x-direction is set to move in the shear direction. The upper fracture surface function is gUx,y, and the lower fracture surface function is gLx,y. To calculate the aperture between the two fracture walls, the following equation can be used [22–24].
(1)ex,y=gUx+us,y−gLx,y+uvgUx+us,y>gLx,y0gUx+us,y≤gLx,y,where ex,y is the upper and lower fracture wall aperture functions, us is the shear direction displacement, and uv is the normal displacement.Whenex,y does not correspond exactly to the spatial coordinates x,y of gUx,y (or gLx,y), it is necessary to use the interpolation method to find the aperture value at the desired spatial location. In Figure 1(a), the gi,j point on the geometric slit surface corresponds to the ei,j′ point in the aperture, but there is no aperture value at this point, so the interpolation method is used to obtain the aperture value at the ei,j′ point.Figure 1
Schematic diagram of fracture aperture interpolation method: (a) schematic diagram of the relationship between the grid points of the fracture wall and the grid points of the aperture; (b) solve the value of the aperture at a certain position (ex,y).The fracture aperture is a function of the spatial coordinatesx,y, so each of the four points can be expressed as x1,y1,ex1,y1, x1,y2,ex2,y1, x2,y1,ex2,y1, and x2,y2,ex2,y2 (Figure 1(b)). If we want to get the aperture value ex,y at a certain position in the middle of the four points, we can use linear interpolation to estimate it:First, interpolatex to find ex,y1 and ex,y2:
(2)ex,y1=ex1,y1⋅x2−xx2−x1+ex2,y1⋅x−x1x2−x1,ex,y2=ex1,y2⋅x2−xx2−x1+ex2,y2⋅x−x1x2−x1.Then, use the same idea to findex,y:
(3)ex,y=ex1,y1⋅y2−yy2−y1+ex,y2⋅y−y1y2−y1.The above equation gives the value of the aperture at any position within the spatial distribution of the aperture. Each pointgi,jx,y on the fracture wall grid can also be found corresponding to the desired aperture value.When discretizing the fracture numerical model, each fracture wall grid point corresponds to a tiny geometric grid cell. Each grid cell can find the corresponding fracture aperture data, and within this cell, the aperture value is a given value; then, the geometric grid cell and the aperture value can form a tiny virtual flat flow model. For each tiny virtual flat plate model, the fracture flow calculation can be performed using the cubic law, and the overall view is an application of the local cubic law. Since the fracture surface and fracture aperture values of each cell are a list of values that vary with space, a rough fracture model with two walls is formed. A triangular grid Figure2(a) or a rectangular grid Figure 2(b) can be used when the fracture data are discrete.Figure 2
Discretizing the rock fracture to the volumetric elements: (a) triangular discrete grid; (b) quadrilateral discrete grid.
(a)(b)


## 2.1.2. From the NS to the Local Cubic Law Equations

The fluid flow in the fracture of the rock is very slow, and neglecting the effect of inertial forces, the Navier-Stokes (NS) equation can be derived from the [12, 25, 26] equation for the fracture flow in the laminar flow state, i.e., the Reynolds equation, also called the local cubic law.The two-dimensional space is expressed as follows:(4)∂∂xρge312μ∂h∂x+∂∂yρge312μ∂h∂y=0.Substituting the aperture functionex,y into Equation (4) yields a clearer expression of the equation [27].
(5)∂∂xρge3x,y12μ∂h∂x+∂∂yρge3x,y12μ∂h∂y=0,where h is the hydraulic head, μ is the dynamic viscosity, ρ is the fluid density, and g is the gravitational acceleration, and ex,y is the aperture function. Figure 2 shows the meaning of Equation (5) graphically by dividing the rough fracture surface into a series of connected small parallel plates.In this paper, the calculation method of using the aperture function to characterize the aperture value and solving the rough fracture flow based on the local cubic law is referred to as the aperture function method (AFM).


## 2.2. Recirculation Zone in Fracture Seepage Flow

Lee et al. [25], Zhou et al. [28], and Zhou et al. [7] found through physical and numerical experiments that in the steady laminar flow state, the fluid in the fracture will generate a recirculation zone (RZ) when flowing through the abrupt change of aperture. The fluid in the RZ also consumes part of the flow energy, and because the fluid in the RZ occurs back on its own, it does not produce effective flow at the outlet of the fracture. The size of the RZ is nonlinearly related to the fluid flow rate and the fracture roughness procedure, and its magnitude cannot be calculated using a simple equation. In addition, Lee et al. [29] found experimentally that under specific conditions, fluid can also produce slip phenomenon at the fracture wall, which also affects the overall seepage flow rate of the fracture.In order to further verify whether a fluid RZ is generated in the fracture model in this paper and to analyze whether the RZ effect is a linear or nonlinear problem, with other geometric and aperture conditions unchanged, a variety of pressure boundary conditions are set, and two methods, NS solution and AFM solution, are used to calculate and compare, respectively, (Table1). The JRC1 fracture model inlets were applied with pressures ranging from 10 Pa to 6000 Pa (−dh/dl=0.01~5.1) for a total of 52 submodels. JRC3 fracture model inlets were applied with pressures ranging from 10 Pa to 890 Pa (−dh/dl=0.01~0.76) for a total of 36 submodels. JRC6 fracture model inlets were applied 25 Pa~650 Pa pressure (−dh/dl=0.02~0.55), respectively, for a total of 26 sub-models. JRC9 fracture model inlets were applied with 10 Pa~1000 Pa pressure (-dh/dl =0.01~0.85), respectively, for a total of 34 submodels. The outlet pressure of all fracture models is 0 Pa. The fluid density in the model is 998.2 kg/m3, and the dynamic viscosity coefficient is 0.0010093 Pa.s. The calculation software was carried out using COMSOL Multiphysics finite element method.Table 1
Basic information of JRC1, JRC3, JRC6, and JRC9 numerical models.
SolutionNumber of elementsElement number ratioMesh size/mHydraulic gradients (-dh/dl)Number of submodelsJRC1NS9017918.763.97×10−8~1.33×10−50.01~5.152AFM1029303.36×10−5~3.11×10−434JRC3NS7433774.93.58×10−7~3.1×10−50.01~0.7636AFM1517901.82×10−5~2.78×10−436JRC6NS2657723.419.83×10−7~8.52×10−50.02~0.5526AFM778842.54×10−5~3.9×10−426JRC9NS1789003.31.80×10−7~6.04×10−50.01~0.8534AFM542027.44×10−5~6.88×10−434


## 2.3. Identification of Vortices in Seepage Flow

The phenomenon of vortices is prevalent in the fluid flow process. Intuitively, vortices represent the rotational motion of the fluid; that is, where there is fluid rotation, there are vortices, and conversely, where there are vortices, there is fluid rotation [30]. In order to give a quantitative mathematical definition of vortex, a parameter Ω is introduced, representing the ratio of the size of the rotating part of the vortex to the size of the total vortex, calculated as [31]
(6)Ω=BF2AF2+BF2+ε,where
(7)AF2=∂u∂x2+∂v∂y2+∂w∂z2+12∂u∂y2+12∂u∂z2+12∂v∂x2+12∂v∂z2+12∂w∂x2+12∂w∂y2+∂u∂y∂v∂x+∂u∂z∂w∂x+∂v∂z∂w∂y,BF2=12∂u∂y2+12∂u∂z2+12∂v∂x2+12∂v∂z2++12∂w∂x2+12∂w∂y2−∂u∂y∂v∂x−∂u∂z∂w∂x+∂v∂z∂w∂y,where F represents the Frobenius parametrization of the matrix, A is the symmetric part of the fluid velocity gradient tensor, and B is the antisymmetric tensor (vorticity).According to its physical meaning, it is obvious thatΩ takes values in the range 0≤Ω≤1, which can be interpreted as the concentration of vorticity, and more specifically, Ω represents the rigidity of fluid motion, and when Ω=1, it represents the fluid doing rigid body rotation. When Ω>0.5 represents that the antisymmetric tensor B dominates over the symmetric tensor A. Therefore, Ω slightly larger than 0.5 can be used as the criterion for vortex identification. Liu et al. [31] proposes to use Ω=0.52 to determine the boundary of vortex.


## 3. Results

### 3.1. Analysis of AFM Calculation Results

In order to verify the rationality of the AFM, this paper compares the physical experimental results of rough fracture fluid flow completed by Zhu et al. [16, 32]. The typical definition of rock fracture roughness curve (JRC) in rock mechanics was adopted for the experiment, and the 1st, 3rd, 6th, and 9th of these curves were selected and combined with a smooth flat plate to form four fracture models with different roughness, each with a fracture length of 100 mm and a minimum aperture value of 0.51 mm. The four models were numbered as JRC1, JRC3, JRC6, and JRC9, respectively. Multiple experiments were conducted for each artificial fracture sample, and the pressure gradient and flow rate were averaged.Based on the AFM, the numerical model of the four fractures was completed and calculated using the finite element numerical calculation software COMSOL Multiphysics. For better display, all the 4 fracture curves are stretched along the direction perpendicular to the hydraulic gradient in the numerical model to form 4 fracture space surfaces. Figure3 shows the numerical model, the aperture spatial distribution, and the fracture seepage field. It can be observed from the figure that the seepage capacity of the fracture is very obviously affected by the minimum aperture.Figure 3
(a, d, g, j) 4 fracture curves (stretched properly for display). (b, e, h, k) Spatial distribution of fracture aperture. (c, f, i, l) Each numerical model fracture surface velocity distribution.
(a)(b)(c)(d)(e)(f)(g)(h)(i)(j)(k)(l)The error values in Table2 are calculated as
(8)ε=qs−qcqs×100%,where ε represents the flow error, qs represents the experimentally measured flow rate per unit width of the fracture, and qc represents qp calculated according to the average aperture and the cube law or qn calculated using the AFM.Table 2
Comparative analysis of experiment results according to average aperture width and aperture interpolation function.
Specimen no.Ps/Paqs/10−3m2/sep (mm)qp (m2/s)εqp/%qn/10−3m2/sen (mm)εqn/%JRC1486.500.1650.910.268620.1850.7731298.200.0350.910.054550.0370.7737JRC3147.500.0971.350.2651740.1060.956998.200.0341.350.0881600.0350.9564JRC6290.200.1873.307.62339770.2170.96716149.000.1063.303.91435930.1120.9685JRC949.800.1204.373.03824320.1391.4991529.900.0814.371.82421460.0831.4992The subscripts represents the actual measurement result of the physical experiment; the subscript p represents the calculation result according to the average aperture and the cube law; the subscript n represents the calculation results of AFM. εqp and εqn represent the error calculated by the flow rate.Table2 shows the comparison of different calculation methods with experimental results. It can be seen that the results calculated by the average aperture and the overall cubic law deviate very much from the experimental measured results. When the fracture roughness increases, the error of the overall cubic law calculation results is even larger. The results calculated by the numerical model using the AFM are in good agreement.


### 3.2. Comparison of the Results of AFM and NS Method

The results of the flow calculated by the AFM are compared with the results of the NS equation solution. We plotted the flow rate variation curves with hydraulic gradient derived from the two calculation methods and plotted the error curves of the two calculation methods (Figure4).Figure 4
Flow per unit width (q) ~ hydraulic gradient (−dh/dl) curf. The solid circle indicates the NS calculation result curve, the solid square indicates the AFM calculation result curve, and the hollow triangle indicates the error in the calculation result of NS vs. AFM.The error curve shows that the AFM calculation result overestimates the permeability of the fracture, causing its calculated flow rate to be large, which is consistent with the results of Yeo et al. [9] and Bauget and Fourar’s study [10]. The errors of NS and AFM calculations are small at low hydraulic gradient and become larger as the hydraulic gradient gradually increases, which is consistent with the results of Oron and Berkowitz [11] and Zimmerman et al. [13]. The increase of fracture roughness also increases the flow error (the overall error of JRC9 model is the largest). The flow error variation also shows a nonlinear variation pattern, rather than a simple linear increase.


### 3.3. Recirculation Zone Analysis

#### 3.3.1. Evolution of the Recirculation Zone

The comparison of the flow results in Table2 shows that the results calculated by the AFM under the low hydraulic gradient conditions are in good agreement, but there is still a certain degree of error. The theory of fracture seepage calculation used in the AFM is the local cubic law, which makes a linear relationship between flow and pressure drop (corresponding to the hydraulic gradient) because it ignores the effect of fluid inertia forces. In contrast, Javadi et al. [17] considers the pressure drop during the flow of fluid in a fracture as VPD and LPD. VPD is due to fluid viscosity and can be calculated using the cubic law.LPD is the pressure drop (or energy loss) that occurs when the fluid flows through a part of the fracture where the aperture changes abruptly. When using the local cubic law-based AFM for seepage calculations, the full LPD cannot be calculated, which is the main reason for the deviation of the AFM calculation results. Zhu et al. [16] and Liu et al. [18] also found the phenomenon of LPD in their studies. The RZ in the fracture flow is one of the important reasons for the LPD. When using the NS solution method for fracture seepage calculations, the extent of the RZ can be determined in the fracture model using streamlines. In this paper, we take JRC3 and JRC6 models as examples to discuss the evolution of RZ and the influence of RZs on fracture seepage.Figures5(a) and 6(a) are discrete mesh of the JRC3 and JRC6 partial models, respectively. Figures 5(b) and 6(b) is a schematic outline of the two models. Figures 5(c)–5(h) and 6(c)–6(h) show the streamlines in the fracture under the typical hydraulic gradient conditions of JRC3 and JRC6, respectively, which can visualize whether there is a recirculation zone or not. Under the low hydraulic gradient (JRC3, −dh/dl=0.13; JRC6, −dh/dl=0.02) conditions, the fluid flow velocity in the fracture seepage field is slow, and the inertial force of the fluid has little influence, and no obvious RZ formation is observed (Figures 5(c) and 6(c)). With the increase of the applied hydraulic gradient, the influence of fluid inertial force gradually increases, and small RZs appear in the region of large changes in the fracture wall (Figures 5(f) and 6(f)), and the range of RZs is small and has little influence on the mainstream zone. The further increase of hydraulic gradient, the influence of fluid inertia force is further revealed, not only to expand the range of RZ but also the emergence of a number of RZs of different sizes, resulting in the narrowing of the mainstream zone, bending (Figures 5(g), 5(h), and 6(f)–6(h)). The RZ not only appears on the side of the rough fracture wall but also produces obvious RZ on the side of the smooth wall, which reflects the strong influence of fluid inertia force. Figures 5(i), 5(j), 6(i), and 6(j) are local streamlines enlargement diagrams; the arrows on the streamline show that the fluid in the RZ is obviously refluxed and does not enter the main flow zone. The fluid in the RZ consumes part of the energy of the fracture flow and generates LPD, but it cannot generate effective flow at the fracture outlet.Figure 5
(a) Discrete mesh (partial). (b) JRC3 model. (c–h) The streamline with different hydraulic gradients (-dh/dl). (i, j) The streamlines in the local area are enlarged, and the recirculation zone and main flow zone are marked.Figure 6
(a) Discrete mesh (partial). (b) JRC6 model. (c–h) The streamline with different hydraulic gradients (-dh/dl). (i, j) The streamlines in the local area are enlarged, and the recirculation zone and main flow zone are marked.Figures7 and 8 shows the variation of the fluid velocity profile with hydraulic gradients on a cut line in the JRC3 and JRC6 fractures. When the hydraulic gradient is small, the velocity distribution on the intercept line has a parabolic shape; when the hydraulic gradient continues to increase, the velocity profile has a bimodal shape due to the creation of the RZ, reflecting the formation of the RZ.Figure 7
JRC3 model. The profiles of the velocity magnitude (U) on the cut line varying with different hydraulic gradients (-dh/dl), and the filled color represents the value of U.Figure 8
JRC6 model. The profiles of the velocity magnitude (U) on the cut line varying with different hydraulic gradients (-dh/dl), and the filled color represents the value of U.The most classical mathematical description of nonlinear flow of a fluid in a fracture is Forchheimer’s law [13, 33]
(9)−∇p=μkv+ρβv2,where ∇p is the hydraulic gradient ∇p=ρgdh/dl, v is the average flow velocity of the fluid under hydraulic gradient conditions v=q/eh, and k is the intrinsic permeability of the fracture k=eh2/12. βm−1 is the inertial resistance coefficient of the fluid, eh is the equivalent hydraulic aperture of the fracture, μ is the viscosity coefficient of the fluid, and ρ is the fluid density. It should be noted that k is the permeability when the hydraulic gradient (−dh/dl) is tiny and the effect of fluid inertial forces can be neglected. When the hydraulic gradient increases, eh changes when there is the creation of a RZ. The keh and the coefficient β can be obtained from the applied hydraulic gradient and the flow fitting (Equation (9)). The fitted equations of the Forchheimer equation for the JRC3 and JRC6 models are labeled in Figure 9, respectively, and the fitted curves are perfect, reflecting the phenomenon that the nonlinear flow in the fracture becomes progressively more pronounced with the increase of the hydraulic gradient.Figure 9
(a) JRC3 model. (b) JRC6 model. Solid lines refer to the fitting curves using the Forchheimer equation, while the dashed lines refer to that using Darcy’s law.
(a)(b)Reynolds number is a quantified parameter to compare the proportional relationship between inertial and viscous forces in the fluid, and for the fluid passing through the fracture, the Reynolds number is calculated as [13]:
(10)Re=ρvDμ=ρQμw,where v is the average flow velocity of the fluid in the fracture, μ is the viscous coefficient of the fluid, ρ is the fluid density, and D is the characteristic length in the flow system; here, the average aperture of the fracture is taken. Q is the fracture flow rate, and w is the width of the fracutre sample perpendicular to the direction of the hydraulic gradient.We calculated the Reynolds number for each hydraulic gradient condition in the JRC3 and JRC6 fracture NS seepage model, respectively, and the variation of the Reynolds number is represented as a curve (blue solid line in Figure9). The obtained Reynolds number can only indicate the trend of increasing fluid inertia force with the hydraulic gradient, but it does not tell the kind of conditions in which significant nonlinear flow occurs.In order to quantify the nonlinear flow, Javadi et al. [34] and Zhou et al. [33] proposed to use the critical Reynolds number (Rec) to make a judgment. Rec indicates the beginning of the transition from fluid flow to non-Darcy flow and can be defined as the Reynolds number when the percentage of nonlinear pressure drop (βρv2) to total pressure drop (μv/k+βρv2) reaches the critical point α. Rec can be obtained from the following equation.
(11)Rec=αeh1−αβk,where α=5% and eh, β, and k are obtained from the Forchheimer’s law fitting equations for the two models JRC3 and JRC6. Using Equation (11), we obtained that Rec=109.27 (Figure 9(a)) for JRC3, corresponding to a hydraulic gradient −dh/dl=0.30 (Table 3) and Rec=74.03 (Figure 9(b)) for JRC6, corresponding to a hydraulic gradient −dh/dl=0.16 (Table 3). Referring to the meaning of the critical Reynolds number, it can be considered that the hydraulic gradient corresponding to the critical Reynolds number, theoretically the fluid inertia force for the overall pressure drop is equal to 5%; combined with the Equation (10) to calculate the variation of Reynolds number, the smaller the hydraulic gradient, the smaller the Reynolds number, the impact of the inertia force is also smaller.Table 3
Main calculation results of JRC3 and JRC6 numerical models.
SolutionCalculation duration/sCalculate duration ratioεqn/%Critical Reynolds numberRecHydraulic gradients corresponding toRec/(-dh/dl)JRC1NS1694158916.58-0.26~ 14.34143.060.39AFM19——JRC3NS1767126311.140.15~ 29.12109.270.3AFM28——JRC6NS586657333.130.73~ 32.4774.030.16AFM8——JRC9NS147721641.33-1.78~ 36.2139.60.0845AFM9——The critical Reynolds number (Rec) indicates the critical point between linear and nonlinear flow in fracture flow. The formation of the RZ is again an important factor for the effect of nonlinearity in rough fracture flow. Therefore, when the Reynolds number of fracture flow is less than the critical Reynolds number (Rec), it is not easy to form the RZ or only a small RZ is formed, and the linear flow is dominant at this time. And when greater than the critical Reynolds number (Rec), the formation of obvious RZ, and with the increase of hydraulic gradient and increase, the nonlinear flow in the fracture flow is more significant.


#### 3.3.2. The Effect of Recirculation Zone

In order to further study the influence of RZ on fracture seepage, we analyze the area occupied by the RZ and the kinetic energy of fluid in the RZ, respectively. The extent of the RZ is difficult to be directly defined, but there must be a velocity vector in the RZ that is opposite to the main seepage direction. The main seepage direction of all models in this paper is thex-direction, i.e., the range of velocity u<0 in the x-direction is part of the RZ. According to the interrelationship between the velocity vector and the flow line in Figure 10, it is found that the range of u<0 can be approximated to half of the area of the RZ. The area of u<0 under a certain hydraulic gradient is calculated and taken as two times its size as the area of the RZ, although there is a certain deviation from the real RZ area, but it can reflect the RZ area change law.Figure 10
The magenta line represents the streamline; the blue arrow represents the velocity vector in thex-direction of the model; the gray shaded part represents the area where the velocity in the x-direction is less than 0.The fluid in the RZ for circular flow, its internal fluid kinetic energy is calculated using a similar method. Calculate the kinetic energy of the fluid with velocityu<0 in the x-direction and then multiply it by 2 as the total kinetic energy of the fluid in the RZ. Then, analyze the law of kinetic energy change in the RZ. The kinetic energy of the fluid is calculated as follows:
(12)EΩ=∫ΩρΔVU22dV,where E is kinetic energy, Ω is the area of RZs, ρ is fluid density, U is fluid flow velocity (scalar), and ΔV is element volume.The main purpose of this paper is to analyze the proportion of the RZ area to the total fracture area and the change of the proportion of the kinetic energy of the RZs to the total kinetic energy of the fluid. The influence of the RZs is analyzed according to the proportional change curve. The calculated results of solving the NS equation for four models, including JRC1, JRC3, JRC6, and JRC9, are counted to show the area of the RZs, the ratio of kinetic energy to total area, and total kinetic energy under each hydraulic gradient, and the change curves are plotted.The JRC3 model is used as an example for illustration (Figure11(a)). The hydraulic gradient of JRC3 fracture model ranges from 0.009 to 2.55, and the ratio of the RZ area to the total area is 0.06% to 12.02%. With the increase of the hydraulic gradient, the percentage change of the RZs area is very obvious. The proportion of the kinetic energy of the RZ area to the total kinetic energy is 0.0002% ~ 0.21%. The kinetic energy percentage is much lower compared to the area, and under the condition of small hydraulic gradient, the effect of kinetic energy in the RZs can even be ignored. The remaining three models JRC1, JRC6, and JRC9 have similar characteristics (Figures 11(b)–11(d)).Figure 11
(a) JRC1 model; (b) JRC3 model; (a) JRC6 model; (b) JRC9 model. The solid red circle line represents the proportion of kinetic energy in the RZ to the total kinetic energy (%). The solid blue square line represents the ratio of the area of the RZ to the total area (%).
(a)(b)(c)(d)Therefore, it can be considered that the main effect of the RZs on the rough fracture seepage is that the fluid occurs back in the RZ by itself and does not produce an effective flow at the fracture outlet. The fluid flow in the RZs is slow, and the loss of kinetic energy of the fluid in the fracture is not significant. However, the RZs squeeze the seepage channel of the fracture, resulting in the reduction of the actual seepage capacity of the fracture. The area occupied by the RZ increases with the hydraulic gradient and shows a nonlinear variation.In addition, the area of RZs and the proportion of kinetic energy do not increase infinitely with the increase of hydraulic gradient, but gradually reach a peak, and then stabilize or slightly decrease. And the rougher the fracture is, the easier it is to reach the peak. In terms of the peak area share of the RZ, the four models JRC1, JRC3, JRC6, and JRC9 reach the peak at the hydraulic gradient of 5.1, 1.68, 0.76 (presumed), and 0.57, respectively.


### 3.4. Analysis of the Effect of Vortices

The RZ is circled according to the shape of the streamline, and its shape is influenced by the number of streamlines and the calculation method of the flow line. The RZ is intuitively seen as a backflow or rotation of the fluid, which is similar to the physical concept of a vortex. However, the vortex determined according to a strict mathematical definition differs significantly from the RZ (Figure12). Vortices are generated near the fracture wall, inside the RZ, at the transition between the reflux and main flow zones, or even within the main flow zone. That is, when there is a velocity difference within the fluid, vortices may be generated. Figure 12 shows that the range of vortices also includes part of the mainstream zone; that is, there is also local fluid rotation within the mainstream zone. This part of the vortex is different from the RZ; it will lose part of the fluid energy but at the same time can produce effective flow.Figure 12
Comparison of the RZ and vortex for JRC1, JRC3, JRC6, and JRC9 models with large hydraulic gradients.To further analyze the effect of vortices on the seepage in the fracture, the value of the vortex parameterΩ in the fluid within the fracture is calculated using Equation (6), and the kinetic energy of the fluid in the range Ω≥0.52 is calculated using Equation (12), along with its percentage of the total kinetic energy. Then, the variation of total kinetic energy, intravortex kinetic energy, and the proportion of intravortex kinetic energy of the fluid is shown as curves using the hydraulic gradient as the horizontal coordinate. Figure 13 shows the curves for the four fracture models. It can be seen that (a) the intravortex kinetic energy increases nonlinearly with the hydraulic gradient. (b) The variation of intravortex kinetic energy ratio is much more complicated. When the fluid Reynolds number is smaller than the critical Reynolds number (Rec), the intravortex kinetic energy ratio shows a decreasing trend. (c) The proportion of kinetic energy in the vortex first decreases, then rises, and then the curve changes to a smooth, indicating that the kinetic energy in the vortex is not completely lost energy, and the fluid in the vortex both rotational motion, but also with the mainstream of the flat motion. (d) Compared with the proportion of kinetic energy in the RZ, the proportion of kinetic energy in the vortex is significantly larger. In the JRC1 model, for example, the highest is 22.16%, and the lowest is 14.85%, so it is obvious that the proportion of kinetic energy in the vortex is much larger and has a greater impact.Figure 13
(a) JRC1 model; (b) JRC3 model; (a) JRC6 model; (b) JRC9 model. The solid blue squares represent the total kinetic energy in the fracture flow; the solid red circles represent the kinetic energy inside the vortex; the solid light blue hollow squares represent the percentage (%) of the kinetic energy inside the vortex to the total kinetic energy.
(a)(b)(c)(d)


## 3.1. Analysis of AFM Calculation Results

In order to verify the rationality of the AFM, this paper compares the physical experimental results of rough fracture fluid flow completed by Zhu et al. [16, 32]. The typical definition of rock fracture roughness curve (JRC) in rock mechanics was adopted for the experiment, and the 1st, 3rd, 6th, and 9th of these curves were selected and combined with a smooth flat plate to form four fracture models with different roughness, each with a fracture length of 100 mm and a minimum aperture value of 0.51 mm. The four models were numbered as JRC1, JRC3, JRC6, and JRC9, respectively. Multiple experiments were conducted for each artificial fracture sample, and the pressure gradient and flow rate were averaged.Based on the AFM, the numerical model of the four fractures was completed and calculated using the finite element numerical calculation software COMSOL Multiphysics. For better display, all the 4 fracture curves are stretched along the direction perpendicular to the hydraulic gradient in the numerical model to form 4 fracture space surfaces. Figure3 shows the numerical model, the aperture spatial distribution, and the fracture seepage field. It can be observed from the figure that the seepage capacity of the fracture is very obviously affected by the minimum aperture.Figure 3
(a, d, g, j) 4 fracture curves (stretched properly for display). (b, e, h, k) Spatial distribution of fracture aperture. (c, f, i, l) Each numerical model fracture surface velocity distribution.
(a)(b)(c)(d)(e)(f)(g)(h)(i)(j)(k)(l)The error values in Table2 are calculated as
(8)ε=qs−qcqs×100%,where ε represents the flow error, qs represents the experimentally measured flow rate per unit width of the fracture, and qc represents qp calculated according to the average aperture and the cube law or qn calculated using the AFM.Table 2
Comparative analysis of experiment results according to average aperture width and aperture interpolation function.
Specimen no.Ps/Paqs/10−3m2/sep (mm)qp (m2/s)εqp/%qn/10−3m2/sen (mm)εqn/%JRC1486.500.1650.910.268620.1850.7731298.200.0350.910.054550.0370.7737JRC3147.500.0971.350.2651740.1060.956998.200.0341.350.0881600.0350.9564JRC6290.200.1873.307.62339770.2170.96716149.000.1063.303.91435930.1120.9685JRC949.800.1204.373.03824320.1391.4991529.900.0814.371.82421460.0831.4992The subscripts represents the actual measurement result of the physical experiment; the subscript p represents the calculation result according to the average aperture and the cube law; the subscript n represents the calculation results of AFM. εqp and εqn represent the error calculated by the flow rate.Table2 shows the comparison of different calculation methods with experimental results. It can be seen that the results calculated by the average aperture and the overall cubic law deviate very much from the experimental measured results. When the fracture roughness increases, the error of the overall cubic law calculation results is even larger. The results calculated by the numerical model using the AFM are in good agreement.


## 3.2. Comparison of the Results of AFM and NS Method

The results of the flow calculated by the AFM are compared with the results of the NS equation solution. We plotted the flow rate variation curves with hydraulic gradient derived from the two calculation methods and plotted the error curves of the two calculation methods (Figure4).Figure 4
Flow per unit width (q) ~ hydraulic gradient (−dh/dl) curf. The solid circle indicates the NS calculation result curve, the solid square indicates the AFM calculation result curve, and the hollow triangle indicates the error in the calculation result of NS vs. AFM.The error curve shows that the AFM calculation result overestimates the permeability of the fracture, causing its calculated flow rate to be large, which is consistent with the results of Yeo et al. [9] and Bauget and Fourar’s study [10]. The errors of NS and AFM calculations are small at low hydraulic gradient and become larger as the hydraulic gradient gradually increases, which is consistent with the results of Oron and Berkowitz [11] and Zimmerman et al. [13]. The increase of fracture roughness also increases the flow error (the overall error of JRC9 model is the largest). The flow error variation also shows a nonlinear variation pattern, rather than a simple linear increase.


## 3.3. Recirculation Zone Analysis

### 3.3.1. Evolution of the Recirculation Zone

The comparison of the flow results in Table2 shows that the results calculated by the AFM under the low hydraulic gradient conditions are in good agreement, but there is still a certain degree of error. The theory of fracture seepage calculation used in the AFM is the local cubic law, which makes a linear relationship between flow and pressure drop (corresponding to the hydraulic gradient) because it ignores the effect of fluid inertia forces. In contrast, Javadi et al. [17] considers the pressure drop during the flow of fluid in a fracture as VPD and LPD. VPD is due to fluid viscosity and can be calculated using the cubic law.LPD is the pressure drop (or energy loss) that occurs when the fluid flows through a part of the fracture where the aperture changes abruptly. When using the local cubic law-based AFM for seepage calculations, the full LPD cannot be calculated, which is the main reason for the deviation of the AFM calculation results. Zhu et al. [16] and Liu et al. [18] also found the phenomenon of LPD in their studies. The RZ in the fracture flow is one of the important reasons for the LPD. When using the NS solution method for fracture seepage calculations, the extent of the RZ can be determined in the fracture model using streamlines. In this paper, we take JRC3 and JRC6 models as examples to discuss the evolution of RZ and the influence of RZs on fracture seepage.Figures5(a) and 6(a) are discrete mesh of the JRC3 and JRC6 partial models, respectively. Figures 5(b) and 6(b) is a schematic outline of the two models. Figures 5(c)–5(h) and 6(c)–6(h) show the streamlines in the fracture under the typical hydraulic gradient conditions of JRC3 and JRC6, respectively, which can visualize whether there is a recirculation zone or not. Under the low hydraulic gradient (JRC3, −dh/dl=0.13; JRC6, −dh/dl=0.02) conditions, the fluid flow velocity in the fracture seepage field is slow, and the inertial force of the fluid has little influence, and no obvious RZ formation is observed (Figures 5(c) and 6(c)). With the increase of the applied hydraulic gradient, the influence of fluid inertial force gradually increases, and small RZs appear in the region of large changes in the fracture wall (Figures 5(f) and 6(f)), and the range of RZs is small and has little influence on the mainstream zone. The further increase of hydraulic gradient, the influence of fluid inertia force is further revealed, not only to expand the range of RZ but also the emergence of a number of RZs of different sizes, resulting in the narrowing of the mainstream zone, bending (Figures 5(g), 5(h), and 6(f)–6(h)). The RZ not only appears on the side of the rough fracture wall but also produces obvious RZ on the side of the smooth wall, which reflects the strong influence of fluid inertia force. Figures 5(i), 5(j), 6(i), and 6(j) are local streamlines enlargement diagrams; the arrows on the streamline show that the fluid in the RZ is obviously refluxed and does not enter the main flow zone. The fluid in the RZ consumes part of the energy of the fracture flow and generates LPD, but it cannot generate effective flow at the fracture outlet.Figure 5
(a) Discrete mesh (partial). (b) JRC3 model. (c–h) The streamline with different hydraulic gradients (-dh/dl). (i, j) The streamlines in the local area are enlarged, and the recirculation zone and main flow zone are marked.Figure 6
(a) Discrete mesh (partial). (b) JRC6 model. (c–h) The streamline with different hydraulic gradients (-dh/dl). (i, j) The streamlines in the local area are enlarged, and the recirculation zone and main flow zone are marked.Figures7 and 8 shows the variation of the fluid velocity profile with hydraulic gradients on a cut line in the JRC3 and JRC6 fractures. When the hydraulic gradient is small, the velocity distribution on the intercept line has a parabolic shape; when the hydraulic gradient continues to increase, the velocity profile has a bimodal shape due to the creation of the RZ, reflecting the formation of the RZ.Figure 7
JRC3 model. The profiles of the velocity magnitude (U) on the cut line varying with different hydraulic gradients (-dh/dl), and the filled color represents the value of U.Figure 8
JRC6 model. The profiles of the velocity magnitude (U) on the cut line varying with different hydraulic gradients (-dh/dl), and the filled color represents the value of U.The most classical mathematical description of nonlinear flow of a fluid in a fracture is Forchheimer’s law [13, 33]
(9)−∇p=μkv+ρβv2,where ∇p is the hydraulic gradient ∇p=ρgdh/dl, v is the average flow velocity of the fluid under hydraulic gradient conditions v=q/eh, and k is the intrinsic permeability of the fracture k=eh2/12. βm−1 is the inertial resistance coefficient of the fluid, eh is the equivalent hydraulic aperture of the fracture, μ is the viscosity coefficient of the fluid, and ρ is the fluid density. It should be noted that k is the permeability when the hydraulic gradient (−dh/dl) is tiny and the effect of fluid inertial forces can be neglected. When the hydraulic gradient increases, eh changes when there is the creation of a RZ. The keh and the coefficient β can be obtained from the applied hydraulic gradient and the flow fitting (Equation (9)). The fitted equations of the Forchheimer equation for the JRC3 and JRC6 models are labeled in Figure 9, respectively, and the fitted curves are perfect, reflecting the phenomenon that the nonlinear flow in the fracture becomes progressively more pronounced with the increase of the hydraulic gradient.Figure 9
(a) JRC3 model. (b) JRC6 model. Solid lines refer to the fitting curves using the Forchheimer equation, while the dashed lines refer to that using Darcy’s law.
(a)(b)Reynolds number is a quantified parameter to compare the proportional relationship between inertial and viscous forces in the fluid, and for the fluid passing through the fracture, the Reynolds number is calculated as [13]:
(10)Re=ρvDμ=ρQμw,where v is the average flow velocity of the fluid in the fracture, μ is the viscous coefficient of the fluid, ρ is the fluid density, and D is the characteristic length in the flow system; here, the average aperture of the fracture is taken. Q is the fracture flow rate, and w is the width of the fracutre sample perpendicular to the direction of the hydraulic gradient.We calculated the Reynolds number for each hydraulic gradient condition in the JRC3 and JRC6 fracture NS seepage model, respectively, and the variation of the Reynolds number is represented as a curve (blue solid line in Figure9). The obtained Reynolds number can only indicate the trend of increasing fluid inertia force with the hydraulic gradient, but it does not tell the kind of conditions in which significant nonlinear flow occurs.In order to quantify the nonlinear flow, Javadi et al. [34] and Zhou et al. [33] proposed to use the critical Reynolds number (Rec) to make a judgment. Rec indicates the beginning of the transition from fluid flow to non-Darcy flow and can be defined as the Reynolds number when the percentage of nonlinear pressure drop (βρv2) to total pressure drop (μv/k+βρv2) reaches the critical point α. Rec can be obtained from the following equation.
(11)Rec=αeh1−αβk,where α=5% and eh, β, and k are obtained from the Forchheimer’s law fitting equations for the two models JRC3 and JRC6. Using Equation (11), we obtained that Rec=109.27 (Figure 9(a)) for JRC3, corresponding to a hydraulic gradient −dh/dl=0.30 (Table 3) and Rec=74.03 (Figure 9(b)) for JRC6, corresponding to a hydraulic gradient −dh/dl=0.16 (Table 3). Referring to the meaning of the critical Reynolds number, it can be considered that the hydraulic gradient corresponding to the critical Reynolds number, theoretically the fluid inertia force for the overall pressure drop is equal to 5%; combined with the Equation (10) to calculate the variation of Reynolds number, the smaller the hydraulic gradient, the smaller the Reynolds number, the impact of the inertia force is also smaller.Table 3
Main calculation results of JRC3 and JRC6 numerical models.
SolutionCalculation duration/sCalculate duration ratioεqn/%Critical Reynolds numberRecHydraulic gradients corresponding toRec/(-dh/dl)JRC1NS1694158916.58-0.26~ 14.34143.060.39AFM19——JRC3NS1767126311.140.15~ 29.12109.270.3AFM28——JRC6NS586657333.130.73~ 32.4774.030.16AFM8——JRC9NS147721641.33-1.78~ 36.2139.60.0845AFM9——The critical Reynolds number (Rec) indicates the critical point between linear and nonlinear flow in fracture flow. The formation of the RZ is again an important factor for the effect of nonlinearity in rough fracture flow. Therefore, when the Reynolds number of fracture flow is less than the critical Reynolds number (Rec), it is not easy to form the RZ or only a small RZ is formed, and the linear flow is dominant at this time. And when greater than the critical Reynolds number (Rec), the formation of obvious RZ, and with the increase of hydraulic gradient and increase, the nonlinear flow in the fracture flow is more significant.


### 3.3.2. The Effect of Recirculation Zone

In order to further study the influence of RZ on fracture seepage, we analyze the area occupied by the RZ and the kinetic energy of fluid in the RZ, respectively. The extent of the RZ is difficult to be directly defined, but there must be a velocity vector in the RZ that is opposite to the main seepage direction. The main seepage direction of all models in this paper is thex-direction, i.e., the range of velocity u<0 in the x-direction is part of the RZ. According to the interrelationship between the velocity vector and the flow line in Figure 10, it is found that the range of u<0 can be approximated to half of the area of the RZ. The area of u<0 under a certain hydraulic gradient is calculated and taken as two times its size as the area of the RZ, although there is a certain deviation from the real RZ area, but it can reflect the RZ area change law.Figure 10
The magenta line represents the streamline; the blue arrow represents the velocity vector in thex-direction of the model; the gray shaded part represents the area where the velocity in the x-direction is less than 0.The fluid in the RZ for circular flow, its internal fluid kinetic energy is calculated using a similar method. Calculate the kinetic energy of the fluid with velocityu<0 in the x-direction and then multiply it by 2 as the total kinetic energy of the fluid in the RZ. Then, analyze the law of kinetic energy change in the RZ. The kinetic energy of the fluid is calculated as follows:
(12)EΩ=∫ΩρΔVU22dV,where E is kinetic energy, Ω is the area of RZs, ρ is fluid density, U is fluid flow velocity (scalar), and ΔV is element volume.The main purpose of this paper is to analyze the proportion of the RZ area to the total fracture area and the change of the proportion of the kinetic energy of the RZs to the total kinetic energy of the fluid. The influence of the RZs is analyzed according to the proportional change curve. The calculated results of solving the NS equation for four models, including JRC1, JRC3, JRC6, and JRC9, are counted to show the area of the RZs, the ratio of kinetic energy to total area, and total kinetic energy under each hydraulic gradient, and the change curves are plotted.The JRC3 model is used as an example for illustration (Figure11(a)). The hydraulic gradient of JRC3 fracture model ranges from 0.009 to 2.55, and the ratio of the RZ area to the total area is 0.06% to 12.02%. With the increase of the hydraulic gradient, the percentage change of the RZs area is very obvious. The proportion of the kinetic energy of the RZ area to the total kinetic energy is 0.0002% ~ 0.21%. The kinetic energy percentage is much lower compared to the area, and under the condition of small hydraulic gradient, the effect of kinetic energy in the RZs can even be ignored. The remaining three models JRC1, JRC6, and JRC9 have similar characteristics (Figures 11(b)–11(d)).Figure 11
(a) JRC1 model; (b) JRC3 model; (a) JRC6 model; (b) JRC9 model. The solid red circle line represents the proportion of kinetic energy in the RZ to the total kinetic energy (%). The solid blue square line represents the ratio of the area of the RZ to the total area (%).
(a)(b)(c)(d)Therefore, it can be considered that the main effect of the RZs on the rough fracture seepage is that the fluid occurs back in the RZ by itself and does not produce an effective flow at the fracture outlet. The fluid flow in the RZs is slow, and the loss of kinetic energy of the fluid in the fracture is not significant. However, the RZs squeeze the seepage channel of the fracture, resulting in the reduction of the actual seepage capacity of the fracture. The area occupied by the RZ increases with the hydraulic gradient and shows a nonlinear variation.In addition, the area of RZs and the proportion of kinetic energy do not increase infinitely with the increase of hydraulic gradient, but gradually reach a peak, and then stabilize or slightly decrease. And the rougher the fracture is, the easier it is to reach the peak. In terms of the peak area share of the RZ, the four models JRC1, JRC3, JRC6, and JRC9 reach the peak at the hydraulic gradient of 5.1, 1.68, 0.76 (presumed), and 0.57, respectively.


## 3.3.1. Evolution of the Recirculation Zone

The comparison of the flow results in Table2 shows that the results calculated by the AFM under the low hydraulic gradient conditions are in good agreement, but there is still a certain degree of error. The theory of fracture seepage calculation used in the AFM is the local cubic law, which makes a linear relationship between flow and pressure drop (corresponding to the hydraulic gradient) because it ignores the effect of fluid inertia forces. In contrast, Javadi et al. [17] considers the pressure drop during the flow of fluid in a fracture as VPD and LPD. VPD is due to fluid viscosity and can be calculated using the cubic law.LPD is the pressure drop (or energy loss) that occurs when the fluid flows through a part of the fracture where the aperture changes abruptly. When using the local cubic law-based AFM for seepage calculations, the full LPD cannot be calculated, which is the main reason for the deviation of the AFM calculation results. Zhu et al. [16] and Liu et al. [18] also found the phenomenon of LPD in their studies. The RZ in the fracture flow is one of the important reasons for the LPD. When using the NS solution method for fracture seepage calculations, the extent of the RZ can be determined in the fracture model using streamlines. In this paper, we take JRC3 and JRC6 models as examples to discuss the evolution of RZ and the influence of RZs on fracture seepage.Figures5(a) and 6(a) are discrete mesh of the JRC3 and JRC6 partial models, respectively. Figures 5(b) and 6(b) is a schematic outline of the two models. Figures 5(c)–5(h) and 6(c)–6(h) show the streamlines in the fracture under the typical hydraulic gradient conditions of JRC3 and JRC6, respectively, which can visualize whether there is a recirculation zone or not. Under the low hydraulic gradient (JRC3, −dh/dl=0.13; JRC6, −dh/dl=0.02) conditions, the fluid flow velocity in the fracture seepage field is slow, and the inertial force of the fluid has little influence, and no obvious RZ formation is observed (Figures 5(c) and 6(c)). With the increase of the applied hydraulic gradient, the influence of fluid inertial force gradually increases, and small RZs appear in the region of large changes in the fracture wall (Figures 5(f) and 6(f)), and the range of RZs is small and has little influence on the mainstream zone. The further increase of hydraulic gradient, the influence of fluid inertia force is further revealed, not only to expand the range of RZ but also the emergence of a number of RZs of different sizes, resulting in the narrowing of the mainstream zone, bending (Figures 5(g), 5(h), and 6(f)–6(h)). The RZ not only appears on the side of the rough fracture wall but also produces obvious RZ on the side of the smooth wall, which reflects the strong influence of fluid inertia force. Figures 5(i), 5(j), 6(i), and 6(j) are local streamlines enlargement diagrams; the arrows on the streamline show that the fluid in the RZ is obviously refluxed and does not enter the main flow zone. The fluid in the RZ consumes part of the energy of the fracture flow and generates LPD, but it cannot generate effective flow at the fracture outlet.Figure 5
(a) Discrete mesh (partial). (b) JRC3 model. (c–h) The streamline with different hydraulic gradients (-dh/dl). (i, j) The streamlines in the local area are enlarged, and the recirculation zone and main flow zone are marked.Figure 6
(a) Discrete mesh (partial). (b) JRC6 model. (c–h) The streamline with different hydraulic gradients (-dh/dl). (i, j) The streamlines in the local area are enlarged, and the recirculation zone and main flow zone are marked.Figures7 and 8 shows the variation of the fluid velocity profile with hydraulic gradients on a cut line in the JRC3 and JRC6 fractures. When the hydraulic gradient is small, the velocity distribution on the intercept line has a parabolic shape; when the hydraulic gradient continues to increase, the velocity profile has a bimodal shape due to the creation of the RZ, reflecting the formation of the RZ.Figure 7
JRC3 model. The profiles of the velocity magnitude (U) on the cut line varying with different hydraulic gradients (-dh/dl), and the filled color represents the value of U.Figure 8
JRC6 model. The profiles of the velocity magnitude (U) on the cut line varying with different hydraulic gradients (-dh/dl), and the filled color represents the value of U.The most classical mathematical description of nonlinear flow of a fluid in a fracture is Forchheimer’s law [13, 33]
(9)−∇p=μkv+ρβv2,where ∇p is the hydraulic gradient ∇p=ρgdh/dl, v is the average flow velocity of the fluid under hydraulic gradient conditions v=q/eh, and k is the intrinsic permeability of the fracture k=eh2/12. βm−1 is the inertial resistance coefficient of the fluid, eh is the equivalent hydraulic aperture of the fracture, μ is the viscosity coefficient of the fluid, and ρ is the fluid density. It should be noted that k is the permeability when the hydraulic gradient (−dh/dl) is tiny and the effect of fluid inertial forces can be neglected. When the hydraulic gradient increases, eh changes when there is the creation of a RZ. The keh and the coefficient β can be obtained from the applied hydraulic gradient and the flow fitting (Equation (9)). The fitted equations of the Forchheimer equation for the JRC3 and JRC6 models are labeled in Figure 9, respectively, and the fitted curves are perfect, reflecting the phenomenon that the nonlinear flow in the fracture becomes progressively more pronounced with the increase of the hydraulic gradient.Figure 9
(a) JRC3 model. (b) JRC6 model. Solid lines refer to the fitting curves using the Forchheimer equation, while the dashed lines refer to that using Darcy’s law.
(a)(b)Reynolds number is a quantified parameter to compare the proportional relationship between inertial and viscous forces in the fluid, and for the fluid passing through the fracture, the Reynolds number is calculated as [13]:
(10)Re=ρvDμ=ρQμw,where v is the average flow velocity of the fluid in the fracture, μ is the viscous coefficient of the fluid, ρ is the fluid density, and D is the characteristic length in the flow system; here, the average aperture of the fracture is taken. Q is the fracture flow rate, and w is the width of the fracutre sample perpendicular to the direction of the hydraulic gradient.We calculated the Reynolds number for each hydraulic gradient condition in the JRC3 and JRC6 fracture NS seepage model, respectively, and the variation of the Reynolds number is represented as a curve (blue solid line in Figure9). The obtained Reynolds number can only indicate the trend of increasing fluid inertia force with the hydraulic gradient, but it does not tell the kind of conditions in which significant nonlinear flow occurs.In order to quantify the nonlinear flow, Javadi et al. [34] and Zhou et al. [33] proposed to use the critical Reynolds number (Rec) to make a judgment. Rec indicates the beginning of the transition from fluid flow to non-Darcy flow and can be defined as the Reynolds number when the percentage of nonlinear pressure drop (βρv2) to total pressure drop (μv/k+βρv2) reaches the critical point α. Rec can be obtained from the following equation.
(11)Rec=αeh1−αβk,where α=5% and eh, β, and k are obtained from the Forchheimer’s law fitting equations for the two models JRC3 and JRC6. Using Equation (11), we obtained that Rec=109.27 (Figure 9(a)) for JRC3, corresponding to a hydraulic gradient −dh/dl=0.30 (Table 3) and Rec=74.03 (Figure 9(b)) for JRC6, corresponding to a hydraulic gradient −dh/dl=0.16 (Table 3). Referring to the meaning of the critical Reynolds number, it can be considered that the hydraulic gradient corresponding to the critical Reynolds number, theoretically the fluid inertia force for the overall pressure drop is equal to 5%; combined with the Equation (10) to calculate the variation of Reynolds number, the smaller the hydraulic gradient, the smaller the Reynolds number, the impact of the inertia force is also smaller.Table 3
Main calculation results of JRC3 and JRC6 numerical models.
SolutionCalculation duration/sCalculate duration ratioεqn/%Critical Reynolds numberRecHydraulic gradients corresponding toRec/(-dh/dl)JRC1NS1694158916.58-0.26~ 14.34143.060.39AFM19——JRC3NS1767126311.140.15~ 29.12109.270.3AFM28——JRC6NS586657333.130.73~ 32.4774.030.16AFM8——JRC9NS147721641.33-1.78~ 36.2139.60.0845AFM9——The critical Reynolds number (Rec) indicates the critical point between linear and nonlinear flow in fracture flow. The formation of the RZ is again an important factor for the effect of nonlinearity in rough fracture flow. Therefore, when the Reynolds number of fracture flow is less than the critical Reynolds number (Rec), it is not easy to form the RZ or only a small RZ is formed, and the linear flow is dominant at this time. And when greater than the critical Reynolds number (Rec), the formation of obvious RZ, and with the increase of hydraulic gradient and increase, the nonlinear flow in the fracture flow is more significant.


## 3.3.2. The Effect of Recirculation Zone

In order to further study the influence of RZ on fracture seepage, we analyze the area occupied by the RZ and the kinetic energy of fluid in the RZ, respectively. The extent of the RZ is difficult to be directly defined, but there must be a velocity vector in the RZ that is opposite to the main seepage direction. The main seepage direction of all models in this paper is thex-direction, i.e., the range of velocity u<0 in the x-direction is part of the RZ. According to the interrelationship between the velocity vector and the flow line in Figure 10, it is found that the range of u<0 can be approximated to half of the area of the RZ. The area of u<0 under a certain hydraulic gradient is calculated and taken as two times its size as the area of the RZ, although there is a certain deviation from the real RZ area, but it can reflect the RZ area change law.Figure 10
The magenta line represents the streamline; the blue arrow represents the velocity vector in thex-direction of the model; the gray shaded part represents the area where the velocity in the x-direction is less than 0.The fluid in the RZ for circular flow, its internal fluid kinetic energy is calculated using a similar method. Calculate the kinetic energy of the fluid with velocityu<0 in the x-direction and then multiply it by 2 as the total kinetic energy of the fluid in the RZ. Then, analyze the law of kinetic energy change in the RZ. The kinetic energy of the fluid is calculated as follows:
(12)EΩ=∫ΩρΔVU22dV,where E is kinetic energy, Ω is the area of RZs, ρ is fluid density, U is fluid flow velocity (scalar), and ΔV is element volume.The main purpose of this paper is to analyze the proportion of the RZ area to the total fracture area and the change of the proportion of the kinetic energy of the RZs to the total kinetic energy of the fluid. The influence of the RZs is analyzed according to the proportional change curve. The calculated results of solving the NS equation for four models, including JRC1, JRC3, JRC6, and JRC9, are counted to show the area of the RZs, the ratio of kinetic energy to total area, and total kinetic energy under each hydraulic gradient, and the change curves are plotted.The JRC3 model is used as an example for illustration (Figure11(a)). The hydraulic gradient of JRC3 fracture model ranges from 0.009 to 2.55, and the ratio of the RZ area to the total area is 0.06% to 12.02%. With the increase of the hydraulic gradient, the percentage change of the RZs area is very obvious. The proportion of the kinetic energy of the RZ area to the total kinetic energy is 0.0002% ~ 0.21%. The kinetic energy percentage is much lower compared to the area, and under the condition of small hydraulic gradient, the effect of kinetic energy in the RZs can even be ignored. The remaining three models JRC1, JRC6, and JRC9 have similar characteristics (Figures 11(b)–11(d)).Figure 11
(a) JRC1 model; (b) JRC3 model; (a) JRC6 model; (b) JRC9 model. The solid red circle line represents the proportion of kinetic energy in the RZ to the total kinetic energy (%). The solid blue square line represents the ratio of the area of the RZ to the total area (%).
(a)(b)(c)(d)Therefore, it can be considered that the main effect of the RZs on the rough fracture seepage is that the fluid occurs back in the RZ by itself and does not produce an effective flow at the fracture outlet. The fluid flow in the RZs is slow, and the loss of kinetic energy of the fluid in the fracture is not significant. However, the RZs squeeze the seepage channel of the fracture, resulting in the reduction of the actual seepage capacity of the fracture. The area occupied by the RZ increases with the hydraulic gradient and shows a nonlinear variation.In addition, the area of RZs and the proportion of kinetic energy do not increase infinitely with the increase of hydraulic gradient, but gradually reach a peak, and then stabilize or slightly decrease. And the rougher the fracture is, the easier it is to reach the peak. In terms of the peak area share of the RZ, the four models JRC1, JRC3, JRC6, and JRC9 reach the peak at the hydraulic gradient of 5.1, 1.68, 0.76 (presumed), and 0.57, respectively.


## 3.4. Analysis of the Effect of Vortices

The RZ is circled according to the shape of the streamline, and its shape is influenced by the number of streamlines and the calculation method of the flow line. The RZ is intuitively seen as a backflow or rotation of the fluid, which is similar to the physical concept of a vortex. However, the vortex determined according to a strict mathematical definition differs significantly from the RZ (Figure12). Vortices are generated near the fracture wall, inside the RZ, at the transition between the reflux and main flow zones, or even within the main flow zone. That is, when there is a velocity difference within the fluid, vortices may be generated. Figure 12 shows that the range of vortices also includes part of the mainstream zone; that is, there is also local fluid rotation within the mainstream zone. This part of the vortex is different from the RZ; it will lose part of the fluid energy but at the same time can produce effective flow.Figure 12
Comparison of the RZ and vortex for JRC1, JRC3, JRC6, and JRC9 models with large hydraulic gradients.To further analyze the effect of vortices on the seepage in the fracture, the value of the vortex parameterΩ in the fluid within the fracture is calculated using Equation (6), and the kinetic energy of the fluid in the range Ω≥0.52 is calculated using Equation (12), along with its percentage of the total kinetic energy. Then, the variation of total kinetic energy, intravortex kinetic energy, and the proportion of intravortex kinetic energy of the fluid is shown as curves using the hydraulic gradient as the horizontal coordinate. Figure 13 shows the curves for the four fracture models. It can be seen that (a) the intravortex kinetic energy increases nonlinearly with the hydraulic gradient. (b) The variation of intravortex kinetic energy ratio is much more complicated. When the fluid Reynolds number is smaller than the critical Reynolds number (Rec), the intravortex kinetic energy ratio shows a decreasing trend. (c) The proportion of kinetic energy in the vortex first decreases, then rises, and then the curve changes to a smooth, indicating that the kinetic energy in the vortex is not completely lost energy, and the fluid in the vortex both rotational motion, but also with the mainstream of the flat motion. (d) Compared with the proportion of kinetic energy in the RZ, the proportion of kinetic energy in the vortex is significantly larger. In the JRC1 model, for example, the highest is 22.16%, and the lowest is 14.85%, so it is obvious that the proportion of kinetic energy in the vortex is much larger and has a greater impact.Figure 13
(a) JRC1 model; (b) JRC3 model; (a) JRC6 model; (b) JRC9 model. The solid blue squares represent the total kinetic energy in the fracture flow; the solid red circles represent the kinetic energy inside the vortex; the solid light blue hollow squares represent the percentage (%) of the kinetic energy inside the vortex to the total kinetic energy.
(a)(b)(c)(d)


## 4. Discussion

For the rough fracture seepage model, it is obvious that solving the NS equation directly can reflect the nonlinear flow state in the fracture more accurately; for example, the fit of Forchheimer equation is very high, and the corresponding calculation results are more accurate. However, its disadvantage is also very obvious; that is, the computational effort is very large. In comparison, the computational efficiency of AFM is much higher than the former.Comparing the four models, the number of discrete elements of the model is much larger than that of AFM when solving the NS equation directly (Table1). This is because AFM only needs to discrete one fracture wall surface. The former has 8.76 times, 4.90 times, 3.41 times, and 3.3 times more than the latter (Table 1).Numerical calculations were performed using the same computer (CPU: Intel(R) Xeon(R) CPU E5-1620 v3, RAM: 16G); the difference in computational efficiency between the two is huge. The JRC1, JRC3, JRC6, and JRC9 models took 169415 s, 176712 s, 58665 s, and 14772 s, respectively, to solve the NS equations (Table3), while the computation time using AFM was 19 s, 28 s, 8 s, and 9 s, respectively (Table 3).The huge computational efficiency advantage of AFM can complete more model computation work in a short time. Because the calculation time-consuming of a single rough fracture model is very short, the seepage calculation of a complex three-dimensional spatial rough fracture network (e.g., a fracture network composed of thousands of rough fractures) can be completed using only an ordinary computer, which is one of our future works. In addition, the geometric model of AFM requires only one fracture wall surface, and the modeling process of numerical model is simpler, while the aperture data can be produced more quickly with the help of third-party software or programming.From this point of view, the huge numerical calculation efficiency of AFM can offset the impact of its calculation error, and the error is not significant at low Reynolds number conditions. This will help researchers to carry out seepage calculations and research work on large scale rock fracture networks.


## 5. Conclusions

From the above analysis and description, we can get several conclusions as follows.(1)
Under the low Reynolds number condition, the numerical calculation work of seepage flow in rough fractures can be done by using the AFM, but the calculation results will produce some small errors(2)
Due to the abrupt change of the fracture wall to produce recirculation zone (RZ), the main effect of RZ is that it crowds the main flow range, making the effective seepage channel of the fracture narrower. The area of the RZs is related to the hydraulic gradient and varies nonlinearly. In addition, the fluid motion inside the RZ is slow, and the loss of kinetic energy of the fluid is not large, accounting for a small proportion of the overall kinetic energy(3)
In addition to the RZ in the rough fracture, vortices are formed, i.e., local fluid rotation. Vortex and RZ area are very different, the fluid in the vortex loss of kinetic energy, but can produce effective flow, and its changes are more complex(4)
Using the local cubic law to solve the rough fracture seepage flow, the main source of error in the calculation results is the LPD. The LPD in fracture seepage includes the area crowding in the RZ, kinetic energy loss in the RZ, kinetic energy loss in the vortex, and pressure drop caused by other reasons(5)
The main advantage of AFM is the efficiency of its numerical calculation. Based on this advantage, AFM can quickly complete the numerical calculation of a huge number of 3D rough fracture networksStrictly speaking, the fracture model in this paper is a single-wall rough fracture model, which has some shortcomings of its own; for example, the nonuniformity and anisotropy of the fracture are not sufficiently reflected. This is mainly due to the need to be consistent with the fracture physical model for comparative analysis. These will be improved gradually in our subsequent studies.



---

*Source: 1828175-2021-10-25.xml*