# Study on the REV Size of Fractured Rock in the Non-Darcy Flow Based on the Dual-Porosity Model

Authors: Yuan Wang; Yu-long Niu; Qiang Feng
Journal: Geofluids
 (2018)

---

## Abstract

For the problem of whether the representative elementary volume (REV) obtained in the Darcy flow is also applicable to the case of the non-Darcy flow, the study on the REV size within the non-Darcy flow is proposed tentatively. The concept of the REV in the non-Darcy flow is based on the definition of the REV. According to the determination of the REV in the Darcy flow, the intrinsic permeabilityk and non-Darcy coefficient β are used simultaneously for the determination of the REV in the non-Darcy flow. The pore pressure cohesive element (PPCE) is developed with the subroutine in ABAQUS. Then the simulation method of the Darcy and non-Darcy flow in the fractured rock mass is built using the PPCE. The proposed plan is examined through the comparison with existing research results. It is validated that this technic is efficient and accurate in simulating the Darcy and non-Darcy flow in the fractured rock mass. Combined with fracture networks generated by Monte Carlo Simulation technique, the PPCE is applied to the study on the REV size. Both conditions of the Darcy and non-Darcy flow are simulated for comparison. The simulation results of this model show that the REV of the non-Darcy flow is inconsistent with the REV of the Darcy flow, and the REV of the non-Darcy flow is more significant than the REV of the Darcy flow. The intrinsic permeability k tensors obtained in the Darcy flow and the non-Darcy flow are basically the same.


---

## Body

## 1. Introduction

Hydraulic property of the fractured rock mass plays a vital role in on-site engineering, like the underground basement, oil storage, CO2 storage, and giant dams. The numerical simulation method can provide essential references for such projects, among which the Finite Element Method (FEM) is the most sophisticated numerical method. In the simulation of flow in the porous media, the equal intrinsic permeability should be determined first. However, as for the fractured porous media, before determining the equivalent intrinsic permeability, there is another critical parameter, the representative elementary volume (REV) to be determined. The REV is defined as the minimum volume of the sampling domain, beyond which the intrinsic permeability of the sampling domain remains essentially constant [1]. Therefore, the determination of the REV size is extremely important and meaningful for understanding the hydraulic behavior of the fractured rock mass.In 1982, Long et al. [2] first studied the DFN (Discrete Fracture Network), in which the scale can be regarded as porous media, and two proposed qualifications: (i) the REV exists; (ii) the equivalent hydraulic property at the certain REV can be approximated by an intrinsic permeability tensor. Oda [3, 4] researched the fractured rock’s REV, anisotropy, and inherent permeability tensor using numerical simulation. Min et al. [5] investigated the REV by the Monte Carlo Simulation technique method for the generation of the DFN based on the results of a site characterization of the field and proposed the “CV” values (coefficient of variation), which is defined as the ratio of standard deviation over the mean value, and the prediction error can be adopted for the determination of the REV. Baghbanan and Jing [6] investigated the intrinsic permeability of fractured rock mass with the constant aperture when considering the correlation between the disturbed fracture aperture and the trace length. The influence of different aperture distributions on the REV is studied by comparing. Min et al. [7] and Baghbanan and Jing [8] examined the effects of stress on the intrinsic permeability and fluid flow in fractured rock mass. The stress ratio was mainly considered for the study, and the results showed that it became more difficult to establish an equivalent tensor and the REV with increasing stress ratios. It is also found that there is a significant difference between correlated and noncorrelated apertures and the fracture length distribution [8]. The studies on the REV in three dimensions have also been reported [9, 10]. The volumetric joint count calculation is applied to the determination of the REV in three dimensions [9]. Thus, in general, there have been many studies on the REV and the intrinsic permeability tensor of the fractured rock mass. However, almost all the reviews of the REV above rely on two fundamental assumptions: (i) the rock matrix is impermeable, and the fluid flow can only occur through fractures. Therefore, the DFN analysis is widely applied; (ii) the fluid flow through the fractures obeys the cubic law. It is also important to take the intrinsic permeability of the rock matrix into consideration [11]. Besides, naturally, the Darcy flow can occur at low flow rates and the fluid flow would enter the non-Darcy flow zone with the high-pressure gradient unaffectedly [10, 12–15]. Thus, considering the intrinsic permeability of the rock matrix and the non-Darcy flow in the fracture is necessary for the further study on the REV of the fractured rock mass.In 1960, Barenblatt et al. [16] first proposed the dual-porosity model for the simulation of the fluid flow in the fractured rock mass. In the dual-porosity model, the fluid can both flow in the rock matrix and fracture. In addition, there is a fluid exchange between them, which is more realistic. Many achievements have been made [11, 16–19]. Ren et al. [20] extended a unified pipe network method (UPNM), which is a dual-porosity model, and applied it to on-site engineering. The results show that the UPNM has excellent advantages in simulating free surface flows in fractured rock slopes and oil underground storage projects with the water curtain system. Excellent reflection of the pore pressure field has been found by adopting the dual-porosity model [11, 21]. Chen et al. [22, 23] proposed the composite element method (CEM) based on the dual-porosity model. The CEM is adopted to investigate the REV of the fractured rock mass and numerical simulation of the on-site project. However, at present, the study on the dual-porosity model is mostly based on self-development, and it is quite complicated to realize this point, which makes it difficult to become popular. Therefore, it becomes more necessary for the realization of the dual-porosity model in mature commercial software, and it can exert advantages by simplifying its application in on-site engineering and related researches.Moreover, many studies on the non-Darcy flow in the fracture have been carried out. Currently, it is considered that inertial force mainly causes the non-Darcy flow. Many experiments include it [10, 24, 25], and it is found that the fluid would present the “weak inertia regime” before totally entering the fully developed turbulence from the Darcy flow. Moreover, it is mentioned that the loss of initial inertial effect mainly results from fracture geometries. As for such inertial losses, the Forchheimer’s equation has been proved by many experiments, which describe such flow behavior very well. Besides, as for the study of flow behavior in the single fracture or intersecting fractures, many works have been done on the factors like equivalent hydraulic apertures, roughness, and pressure gradients [26]. Only a few works focused on the multifractures, especially the non-Darcy flow behavior in the fracture network and its REV size. Zimmerman et al. [27] measured and computed the intrinsic permeability of a natural sandstone with fractures. When the Reynolds number is larger than 20, both results of the experiment and numerical simulation exhibit a Forchheimer-type regime, in which the non-Darcy flow pressure drop is quadratic in the flow rate. Many non-Darcy flow experiments in the fracture network have also been conducted [28, 29]. It has been reported that when the hydraulic gradient is large enough, the inertial effect becomes nonnegligible and the nonlinear characteristic of the fluid flow should be taken into consideration necessarily [28]. What is more, the study on the Darcy and non-Darcy flow based on the dual-porosity model is rare. The 2D models, based on the dual-porosity model and mathematic analysis, are established by taking the fracture as interfaces within the Forchheimer flow while the flow in the surrounding matrix is such that Darcy’s law is adequate [30–32], but those studies only set few fractures in porous media. Generally, to take both the dual-porosity model and the non-Darcy flow behavior in the fracture network into consideration, the REV size should be further studied.Throughout the previous literature, on one hand, most studies on the REV of the 2D fractured rock rely on the discrete model, and it is assumed that the rock matrix is impermeable and the fluid flow only occurs in the fracture. Considering that the pressure in the matrix would affect the velocity of flow in the breach, the REV size of fractured rock based on the dual-porosity model is investigated. Besides, the simulation method of the Darcy and non-Darcy flow in the fractured rock mass is built, which is conducive to promotion. On the other hand, almost all the studies on the REV involve the fluid flow in the fracture obeying the Cubic law, and as for the nonlinear flow, most studies focus on the single fracture and the fracture network. The survey on the REV in the non-Darcy flow has not been conducted in previous work. Based on the review of the non-linear flow behavior, the REV size of fractured rock in the non-Darcy flow is tentatively investigated.In this paper, firstly, based on the dual-porosity model and the mathematics model [30], a simulation method is proposed with the pore pressure cohesive element (PPCE) in ABAQUS, which can simulate both the Darcy flow in the matrix and the non-Darcy flow in the rock fracture. According to the interpretation of the REV in the Darcy flow, the definition of the REV in the non-Darcy flow is given, when the concept of the REV in the non-Darcy flow is proposed based on the description of the REV. Then, focusing on the problem whether the representative elementary volume (REV) obtained in the Darcy flow is also suitable for the case of the non-Darcy flow, the REV size in the non-Darcy flow with the stochastic fracture network generated by the Monte Carlo method is investigated. In this paper, the REV of fracture rock in the Darcy flow is compared with that in the non-Darcy flow, when paying attention to the difference of REV sizes between the Darcy and non-Darcy flow.


## 2. The Simulation Theory and Method of the Darcy and Non-Darcy Flow Based on the Dual-Porosity Model

### 2.1. Darcy’s Law and Forchheimer’s Law

Considering that the velocity of the fluid in the porous rock matrix is much lower than in the fracture, the fluid flow in the matrix always obeys Darcy’s law. For the single-phase flow and incompressible flow in porous media, Darcy’s law can be expressed as (in the absence of gravity)(1)v+Kmμ∇P=0,where v[LT-1] is the flow velocity; Km[L2] is the intrinsic permeability of porous media; ∇P[ML-2T-2] is the pressure gradient in the porous media.When the flow velocity is low in the fracture, the relationship between the volumetric flow rate and the pressure gradient is linear, and it can be described as the Cubic law:(2)Q=-we312μ∇P=-KfAhμ∇P,where Q[L3T-1] is the volumetric flow rate or discharge; w[L] is the fracture aperture; e[L] is the aperture of the idealized parallel fracture, and it is the equivalent hydraulic aperture eh for the rough fracture; μ[ML-1T-1] is the dynamic viscosity of the fluid; Kf[L2] is the intrinsic permeability defined as Kf=e2/12; Ah is the cross-section area. When the slope of the flow rate and the pressure gradient no longer coincide, the flow would enter into the non-Darcy regime, and its constitutive relationship can be described as Forchheimer’s law:(3)-∇P=AQ+BQ2,where ∇P[ML-2T-2] is the pressure gradient along the flow direction; A[ML-5T-1] and B[ML-8] are the coefficients, which describe the energy losses of the flow caused by viscosity and inertia, respectively. A and B in (3) are commonly written as follows:(4)A=μkAh=12μweh3(5)B=βρAh2=βρw2eh2,where β[L-1] is called the non-Darcy coefficient, whose value is related to the aperture and roughness of the fracture. Reynolds number Re is defined as the ratio of inertial force to viscous force, while Re can be calculated as follows:(6)Re=ρQμw=ρvehμ,where ρ[ML-3] is the density of fluid.


### 2.2. The Pore Pressure Cohesive Element

The pore pressure cohesive element (PPCE, 6-node element) is developed from a cohesive element (4-node element), which can simulate the fluid flow in fractured rock [33]. Chen et al. [34] and Carrier and Granet [35] study hydraulic fractures with the PPCE, indicating the PPCE can simulate the fluid flow in fractured rock precisely. As shown in Figure 1, in the two-dimensional case, nodes 3, 4, 5, and 6 have not only a freedom degree of 8 (the default pore pressure of ABAQUS is 8), but also the freedom degree of 1 and 2 (the freedom degree of the displacement), whereas nodes 1 and 2 only have a freedom degree of 8. In the simulation, nodes 1, 3, and 5 (2, 4, and 6) have different pore pressure values, but six nodes have the same coordinates when the displacement is excluded, so that the entire unit looks like a line unit. The fluid flow in the PPCE includes the tangential flow within the fracture and the normal flow across the fracture, as shown in Figure 1, while Figure 2 shows the meshing examples for 2D intersecting fractures with the PPCE in ABAQUS.Figure 1
Fluid flow within the PPCE in 2D.Figure 2
The meshing example for 2D intersecting fractures with the PPCE.In ABAQUS, the default constitutive relationship of tangential flow is described as (2), and the constitutive relationship of normal flow is described as follows:(7)Q=cPi-Po,where c[M-1L2T] is the fluid leak-off coefficients at the top and bottom element surfaces; Pi and Po[ML-1T-2] are the midface pore pressure and the pore pressure on the sides, respectively.In the flow exchange of fractures and the matrix, the entire basin must follow the law of mass conservation. As shown in Figure3, m1 and m2 are porous media and f is the fracture. Vector n is a vector normal to f. In m1, m2, for single-phase incompressible fluids, it can be expressed as follows:(8)divv=0inm1,m2.In the flow exchange of fractures and the matrix, the law of mass conservation(9)divQf=v1·n1+v2·n2onf,where Qf=∫-d/2d/2vf,τdn, vf,τ being the tangential flow rate in the fracture; v1·n1, v2·n2 represent the flow rate into the fissure and out of the fracture, respectively.Figure 3
Fractured media.


### 2.3. The Analysis Method Based on the PPCE for Simulation of Non-Darcy in Fractured Rock

For the standard calculation in linear FEM, when taking the theory of flow for example, the relationship can be written as follows:(10)Q=KP,where Q is the volumetric flow rate of the node, K is hydraulic conductivity of the matrix, and P is pore pressure of the node.If the relationship between the flow rate and pressure is linear, then the intrinsic permeability matrix is a constant matrix. However, for the non-Darcy flow, the problem is complicated. The Forchheimer flow rateQ can be obtained by solving (4):(11)Q=2-∇PA+A2+4B-∇PThen,(12)K=2A+A2+4B-∇P.From (12), it is obvious that the hydraulic conductivity K is no longer a constant one as pore pressure changes. Figure 4 abstractly describes the nonlinear relationship between the flow rate and pressure and the “slope” of a certain point on the curve is the intrinsic permeability matrix K. There are many ways of solving this kind of nonlinear problem, such as the iterative method, the incremental method, and the incremental iteration method in FEM.Figure 4
Relationship between the flow rate and pore pressure.As for the non-Darcy flow, the constitutive default of tangential flow can be adapted with a subroutine (nodes 1 and 2 in Figure1), and then the case of the non-Darcy flow can be simulated in fractured media. For a line element with 2 nodes, (10) can be written as follows:(13)Q1Q2=1LK-K-KKP1P2,where Q1 and Q2 are the flow rate of the node,L is the length of the line element,K is hydraulic conductivity, and P1 and P2 are pore pressure of the node.In this study, the primary incremental method is used for the calculation of the non-Darcy flow. As shown in Figure5, specific steps are as follows. Each increment (ΔPl or ΔQl) is small enough in the calculation. At the first increment, the quadratic term in (3) is ignored, and then Pl-1 and Ql-1 can be obtained for the calculation of the intrinsic permeability matrix in the next increase. Repeat the above steps until the non-Darcy flow is steady. At the end of each increment, there is an error justification, Δε<ε. If the error is beyond the threshold, the increment should be smaller for precision.Figure 5
Relationship between the flow rate and pore pressure in the incremental method.


### 2.4. Validation

To validate the correctness and effectiveness of the PPCE and its subroutine, several cases with the Darcy flow and the Forchheimer flow are calculated to compare with previous studies.


#### 2.4.1. Fractured Media within the Darcy Flow in the Fracture

In this case, the steady state of fractured media is simulated. The dimensions and boundaries are shown in Figure6. The length of the fracture is L=8m, Kf=0.08 m2, and Km=8×10-6 m2.Figure 6
2D Fractured media with applied pressure at top and bottom.The discretization and pressure fields of this field obtained by ABAQUS-PPCE are compared with DFM (the discrete fracture model) and FM-Meshfree (fracture mapping with mesh-free), which are proposed by Lamb [36], as shown in Figures 7 and 8. To be specific, Figure 7 shows the approaches of four kinds of algorithms for the same model. The mesh nodes should be set on the fracture for the FM-Meshfree model while nodes and element edges need to be placed on the fracture for ABAQUS-PPCE and DFM. It can tell from Figure 8 that the pore pressure contours obtained by four algorithms are consistent, and Figure 9 shows the profile of pore pressure for a vertical section taken through the center of the fractured media for comparison. It can be found that ABAQUS-PPCE’s pore pressure curve and other curves are consistent with the result calculated in Lamb [36], which can validate that ABAQUS-PPCE has its effectiveness in fractured media within the Darcy flow in the fracture simulation.Figure 7
Domain discretization of each method.
(a)
ABAQUS-PPCE (b)
DFM (c)
FM-Meshfree (d)
FM-Meshfree (irregular nodal distribution)Figure 8
Pore pressure fields obtained by ABAQUS-PPCE and other methods.
(a)
ABAQUS-PPCE (b)
DFM (c)
FM-Meshfree (d)
FM-Meshfree (irregular nodal distribution)Figure 9
The pressure profile for a vertical section taken through the center of the fractured media.


#### 2.4.2. Free Surface Flow in the Fractured Rock Dam (The Darcy Flow)

In this case, ABAQUS-PPCE is used to simulate one of the most common problems in geotechnical engineering: unconfined flow analysis. The dam model’s dimensions and water elevations of the upstream and the downstream are shown in Figure10(a). To evaluate the effects of the fracture on free surface flow, a fracture with the aperture e=0.001m and the length L=2.5m is embedded in the dam as shown in Figure 10(b). The results obtained by ABAQUS-PPCE and Ren et al. [20] are compared in Figures 11 and 12. When the experiment result calculated by ABAQUS-PPCE is slightly higher than others (about 0.4 m), it can be found that the curves of free surface flow coincide with each other. Thus, it is reasonable to conclude that ABAQUS-PPCE is also correct and efficient in solving unconfined flow problems.Figure 10
Dimensions and boundaries of the dam model.
(a)
Porous media (b)
Fractured mediaFigure 11
Pore pressure fields obtained by ABAQUS-PPCE.
(a)
Porous media (b)
Fractured mediaFigure 12
Comparison of the results obtained by different methods.


#### 2.4.3. Fractured Media within Forchheimer Flow in the Fracture

Frih et al. [30, 32] validated the correctness of the interface model by being compared with the standard model, both in the case of the Darcy flow and in the case of the Forchheimer flow. The dimensions and boundaries are shown in Figure 13, where Kf=10-9m2, e=0.01 m (ignoring cubic law), β=10000m-1, and Km=10-12m2. The steady flow results of pore pressure and velocity simulated by ABAQUS-PPCE are shown in Figure 14. The comparison of tangential velocities between ABAQUS-PPCE and the interface model is given in Figure 15. According to the velocity along the fracture in Figure 15 and the recalculation equation (6), it can be seen that Re is about 3000. Combined with previous studies, the flow in this fracture is turbulent. It is clear that both the Darcy flow and Forchheimer flow of ABAQUS-PPCE results are close to the interface model. Besides, Figure 15 shows that the flow velocity in the fracture is affected by the pore pressure in porous media, which indicates it has the role of a barrier. Obviously, this situation cannot be achieved by the discrete media model, so the dual-media model has its advantages.Figure 13
Dimensions and boundaries of the model with one fracture.Figure 14
(a) The 3D pore pressure field and (b) velocity for the Forchheimer flow case.
(a) (b)Figure 15
Darcy and Forchheimer flow velocity in the fracture given by the interface model and ABAQUS-PPCE.Above all, by comparing ABAQUS-PPCE’s results with other existing examples, it can be found that ABAQUS-PPCE is effective and accurate when simulating the Darcy flow and the non-Darcy flow (Forchheimer flow) within the fractured porous media.


## 2.1. Darcy’s Law and Forchheimer’s Law

Considering that the velocity of the fluid in the porous rock matrix is much lower than in the fracture, the fluid flow in the matrix always obeys Darcy’s law. For the single-phase flow and incompressible flow in porous media, Darcy’s law can be expressed as (in the absence of gravity)(1)v+Kmμ∇P=0,where v[LT-1] is the flow velocity; Km[L2] is the intrinsic permeability of porous media; ∇P[ML-2T-2] is the pressure gradient in the porous media.When the flow velocity is low in the fracture, the relationship between the volumetric flow rate and the pressure gradient is linear, and it can be described as the Cubic law:(2)Q=-we312μ∇P=-KfAhμ∇P,where Q[L3T-1] is the volumetric flow rate or discharge; w[L] is the fracture aperture; e[L] is the aperture of the idealized parallel fracture, and it is the equivalent hydraulic aperture eh for the rough fracture; μ[ML-1T-1] is the dynamic viscosity of the fluid; Kf[L2] is the intrinsic permeability defined as Kf=e2/12; Ah is the cross-section area. When the slope of the flow rate and the pressure gradient no longer coincide, the flow would enter into the non-Darcy regime, and its constitutive relationship can be described as Forchheimer’s law:(3)-∇P=AQ+BQ2,where ∇P[ML-2T-2] is the pressure gradient along the flow direction; A[ML-5T-1] and B[ML-8] are the coefficients, which describe the energy losses of the flow caused by viscosity and inertia, respectively. A and B in (3) are commonly written as follows:(4)A=μkAh=12μweh3(5)B=βρAh2=βρw2eh2,where β[L-1] is called the non-Darcy coefficient, whose value is related to the aperture and roughness of the fracture. Reynolds number Re is defined as the ratio of inertial force to viscous force, while Re can be calculated as follows:(6)Re=ρQμw=ρvehμ,where ρ[ML-3] is the density of fluid.


## 2.2. The Pore Pressure Cohesive Element

The pore pressure cohesive element (PPCE, 6-node element) is developed from a cohesive element (4-node element), which can simulate the fluid flow in fractured rock [33]. Chen et al. [34] and Carrier and Granet [35] study hydraulic fractures with the PPCE, indicating the PPCE can simulate the fluid flow in fractured rock precisely. As shown in Figure 1, in the two-dimensional case, nodes 3, 4, 5, and 6 have not only a freedom degree of 8 (the default pore pressure of ABAQUS is 8), but also the freedom degree of 1 and 2 (the freedom degree of the displacement), whereas nodes 1 and 2 only have a freedom degree of 8. In the simulation, nodes 1, 3, and 5 (2, 4, and 6) have different pore pressure values, but six nodes have the same coordinates when the displacement is excluded, so that the entire unit looks like a line unit. The fluid flow in the PPCE includes the tangential flow within the fracture and the normal flow across the fracture, as shown in Figure 1, while Figure 2 shows the meshing examples for 2D intersecting fractures with the PPCE in ABAQUS.Figure 1
Fluid flow within the PPCE in 2D.Figure 2
The meshing example for 2D intersecting fractures with the PPCE.In ABAQUS, the default constitutive relationship of tangential flow is described as (2), and the constitutive relationship of normal flow is described as follows:(7)Q=cPi-Po,where c[M-1L2T] is the fluid leak-off coefficients at the top and bottom element surfaces; Pi and Po[ML-1T-2] are the midface pore pressure and the pore pressure on the sides, respectively.In the flow exchange of fractures and the matrix, the entire basin must follow the law of mass conservation. As shown in Figure3, m1 and m2 are porous media and f is the fracture. Vector n is a vector normal to f. In m1, m2, for single-phase incompressible fluids, it can be expressed as follows:(8)divv=0inm1,m2.In the flow exchange of fractures and the matrix, the law of mass conservation(9)divQf=v1·n1+v2·n2onf,where Qf=∫-d/2d/2vf,τdn, vf,τ being the tangential flow rate in the fracture; v1·n1, v2·n2 represent the flow rate into the fissure and out of the fracture, respectively.Figure 3
Fractured media.


## 2.3. The Analysis Method Based on the PPCE for Simulation of Non-Darcy in Fractured Rock

For the standard calculation in linear FEM, when taking the theory of flow for example, the relationship can be written as follows:(10)Q=KP,where Q is the volumetric flow rate of the node, K is hydraulic conductivity of the matrix, and P is pore pressure of the node.If the relationship between the flow rate and pressure is linear, then the intrinsic permeability matrix is a constant matrix. However, for the non-Darcy flow, the problem is complicated. The Forchheimer flow rateQ can be obtained by solving (4):(11)Q=2-∇PA+A2+4B-∇PThen,(12)K=2A+A2+4B-∇P.From (12), it is obvious that the hydraulic conductivity K is no longer a constant one as pore pressure changes. Figure 4 abstractly describes the nonlinear relationship between the flow rate and pressure and the “slope” of a certain point on the curve is the intrinsic permeability matrix K. There are many ways of solving this kind of nonlinear problem, such as the iterative method, the incremental method, and the incremental iteration method in FEM.Figure 4
Relationship between the flow rate and pore pressure.As for the non-Darcy flow, the constitutive default of tangential flow can be adapted with a subroutine (nodes 1 and 2 in Figure1), and then the case of the non-Darcy flow can be simulated in fractured media. For a line element with 2 nodes, (10) can be written as follows:(13)Q1Q2=1LK-K-KKP1P2,where Q1 and Q2 are the flow rate of the node,L is the length of the line element,K is hydraulic conductivity, and P1 and P2 are pore pressure of the node.In this study, the primary incremental method is used for the calculation of the non-Darcy flow. As shown in Figure5, specific steps are as follows. Each increment (ΔPl or ΔQl) is small enough in the calculation. At the first increment, the quadratic term in (3) is ignored, and then Pl-1 and Ql-1 can be obtained for the calculation of the intrinsic permeability matrix in the next increase. Repeat the above steps until the non-Darcy flow is steady. At the end of each increment, there is an error justification, Δε<ε. If the error is beyond the threshold, the increment should be smaller for precision.Figure 5
Relationship between the flow rate and pore pressure in the incremental method.


## 2.4. Validation

To validate the correctness and effectiveness of the PPCE and its subroutine, several cases with the Darcy flow and the Forchheimer flow are calculated to compare with previous studies.


### 2.4.1. Fractured Media within the Darcy Flow in the Fracture

In this case, the steady state of fractured media is simulated. The dimensions and boundaries are shown in Figure6. The length of the fracture is L=8m, Kf=0.08 m2, and Km=8×10-6 m2.Figure 6
2D Fractured media with applied pressure at top and bottom.The discretization and pressure fields of this field obtained by ABAQUS-PPCE are compared with DFM (the discrete fracture model) and FM-Meshfree (fracture mapping with mesh-free), which are proposed by Lamb [36], as shown in Figures 7 and 8. To be specific, Figure 7 shows the approaches of four kinds of algorithms for the same model. The mesh nodes should be set on the fracture for the FM-Meshfree model while nodes and element edges need to be placed on the fracture for ABAQUS-PPCE and DFM. It can tell from Figure 8 that the pore pressure contours obtained by four algorithms are consistent, and Figure 9 shows the profile of pore pressure for a vertical section taken through the center of the fractured media for comparison. It can be found that ABAQUS-PPCE’s pore pressure curve and other curves are consistent with the result calculated in Lamb [36], which can validate that ABAQUS-PPCE has its effectiveness in fractured media within the Darcy flow in the fracture simulation.Figure 7
Domain discretization of each method.
(a)
ABAQUS-PPCE (b)
DFM (c)
FM-Meshfree (d)
FM-Meshfree (irregular nodal distribution)Figure 8
Pore pressure fields obtained by ABAQUS-PPCE and other methods.
(a)
ABAQUS-PPCE (b)
DFM (c)
FM-Meshfree (d)
FM-Meshfree (irregular nodal distribution)Figure 9
The pressure profile for a vertical section taken through the center of the fractured media.


### 2.4.2. Free Surface Flow in the Fractured Rock Dam (The Darcy Flow)

In this case, ABAQUS-PPCE is used to simulate one of the most common problems in geotechnical engineering: unconfined flow analysis. The dam model’s dimensions and water elevations of the upstream and the downstream are shown in Figure10(a). To evaluate the effects of the fracture on free surface flow, a fracture with the aperture e=0.001m and the length L=2.5m is embedded in the dam as shown in Figure 10(b). The results obtained by ABAQUS-PPCE and Ren et al. [20] are compared in Figures 11 and 12. When the experiment result calculated by ABAQUS-PPCE is slightly higher than others (about 0.4 m), it can be found that the curves of free surface flow coincide with each other. Thus, it is reasonable to conclude that ABAQUS-PPCE is also correct and efficient in solving unconfined flow problems.Figure 10
Dimensions and boundaries of the dam model.
(a)
Porous media (b)
Fractured mediaFigure 11
Pore pressure fields obtained by ABAQUS-PPCE.
(a)
Porous media (b)
Fractured mediaFigure 12
Comparison of the results obtained by different methods.


### 2.4.3. Fractured Media within Forchheimer Flow in the Fracture

Frih et al. [30, 32] validated the correctness of the interface model by being compared with the standard model, both in the case of the Darcy flow and in the case of the Forchheimer flow. The dimensions and boundaries are shown in Figure 13, where Kf=10-9m2, e=0.01 m (ignoring cubic law), β=10000m-1, and Km=10-12m2. The steady flow results of pore pressure and velocity simulated by ABAQUS-PPCE are shown in Figure 14. The comparison of tangential velocities between ABAQUS-PPCE and the interface model is given in Figure 15. According to the velocity along the fracture in Figure 15 and the recalculation equation (6), it can be seen that Re is about 3000. Combined with previous studies, the flow in this fracture is turbulent. It is clear that both the Darcy flow and Forchheimer flow of ABAQUS-PPCE results are close to the interface model. Besides, Figure 15 shows that the flow velocity in the fracture is affected by the pore pressure in porous media, which indicates it has the role of a barrier. Obviously, this situation cannot be achieved by the discrete media model, so the dual-media model has its advantages.Figure 13
Dimensions and boundaries of the model with one fracture.Figure 14
(a) The 3D pore pressure field and (b) velocity for the Forchheimer flow case.
(a) (b)Figure 15
Darcy and Forchheimer flow velocity in the fracture given by the interface model and ABAQUS-PPCE.Above all, by comparing ABAQUS-PPCE’s results with other existing examples, it can be found that ABAQUS-PPCE is effective and accurate when simulating the Darcy flow and the non-Darcy flow (Forchheimer flow) within the fractured porous media.


## 2.4.1. Fractured Media within the Darcy Flow in the Fracture

In this case, the steady state of fractured media is simulated. The dimensions and boundaries are shown in Figure6. The length of the fracture is L=8m, Kf=0.08 m2, and Km=8×10-6 m2.Figure 6
2D Fractured media with applied pressure at top and bottom.The discretization and pressure fields of this field obtained by ABAQUS-PPCE are compared with DFM (the discrete fracture model) and FM-Meshfree (fracture mapping with mesh-free), which are proposed by Lamb [36], as shown in Figures 7 and 8. To be specific, Figure 7 shows the approaches of four kinds of algorithms for the same model. The mesh nodes should be set on the fracture for the FM-Meshfree model while nodes and element edges need to be placed on the fracture for ABAQUS-PPCE and DFM. It can tell from Figure 8 that the pore pressure contours obtained by four algorithms are consistent, and Figure 9 shows the profile of pore pressure for a vertical section taken through the center of the fractured media for comparison. It can be found that ABAQUS-PPCE’s pore pressure curve and other curves are consistent with the result calculated in Lamb [36], which can validate that ABAQUS-PPCE has its effectiveness in fractured media within the Darcy flow in the fracture simulation.Figure 7
Domain discretization of each method.
(a)
ABAQUS-PPCE (b)
DFM (c)
FM-Meshfree (d)
FM-Meshfree (irregular nodal distribution)Figure 8
Pore pressure fields obtained by ABAQUS-PPCE and other methods.
(a)
ABAQUS-PPCE (b)
DFM (c)
FM-Meshfree (d)
FM-Meshfree (irregular nodal distribution)Figure 9
The pressure profile for a vertical section taken through the center of the fractured media.


## 2.4.2. Free Surface Flow in the Fractured Rock Dam (The Darcy Flow)

In this case, ABAQUS-PPCE is used to simulate one of the most common problems in geotechnical engineering: unconfined flow analysis. The dam model’s dimensions and water elevations of the upstream and the downstream are shown in Figure10(a). To evaluate the effects of the fracture on free surface flow, a fracture with the aperture e=0.001m and the length L=2.5m is embedded in the dam as shown in Figure 10(b). The results obtained by ABAQUS-PPCE and Ren et al. [20] are compared in Figures 11 and 12. When the experiment result calculated by ABAQUS-PPCE is slightly higher than others (about 0.4 m), it can be found that the curves of free surface flow coincide with each other. Thus, it is reasonable to conclude that ABAQUS-PPCE is also correct and efficient in solving unconfined flow problems.Figure 10
Dimensions and boundaries of the dam model.
(a)
Porous media (b)
Fractured mediaFigure 11
Pore pressure fields obtained by ABAQUS-PPCE.
(a)
Porous media (b)
Fractured mediaFigure 12
Comparison of the results obtained by different methods.


## 2.4.3. Fractured Media within Forchheimer Flow in the Fracture

Frih et al. [30, 32] validated the correctness of the interface model by being compared with the standard model, both in the case of the Darcy flow and in the case of the Forchheimer flow. The dimensions and boundaries are shown in Figure 13, where Kf=10-9m2, e=0.01 m (ignoring cubic law), β=10000m-1, and Km=10-12m2. The steady flow results of pore pressure and velocity simulated by ABAQUS-PPCE are shown in Figure 14. The comparison of tangential velocities between ABAQUS-PPCE and the interface model is given in Figure 15. According to the velocity along the fracture in Figure 15 and the recalculation equation (6), it can be seen that Re is about 3000. Combined with previous studies, the flow in this fracture is turbulent. It is clear that both the Darcy flow and Forchheimer flow of ABAQUS-PPCE results are close to the interface model. Besides, Figure 15 shows that the flow velocity in the fracture is affected by the pore pressure in porous media, which indicates it has the role of a barrier. Obviously, this situation cannot be achieved by the discrete media model, so the dual-media model has its advantages.Figure 13
Dimensions and boundaries of the model with one fracture.Figure 14
(a) The 3D pore pressure field and (b) velocity for the Forchheimer flow case.
(a) (b)Figure 15
Darcy and Forchheimer flow velocity in the fracture given by the interface model and ABAQUS-PPCE.Above all, by comparing ABAQUS-PPCE’s results with other existing examples, it can be found that ABAQUS-PPCE is effective and accurate when simulating the Darcy flow and the non-Darcy flow (Forchheimer flow) within the fractured porous media.


## 3. The Stochastic Fracture Network and Darcy and Non-Darcy Flow Simulations of Fractured Rock

In this section, a fundamental question in the field of flow in fractured rock, REV, and equivalent hydraulic properties is researched. Firstly, fracture networks are generated by Monte Carlo Simulation technique. Then REV and equivalent hydraulic properties are studies in the cases of Darcy and non-Darcy flows. Different criteria are adopted to determine whether there is inconsistency in the REV in the case of the Darcy and non-Darcy flow.


### 3.1. Generation of the Stochastic Fracture Network

Table1 shows the necessary parameters of the fracture network for this study, and the intrinsic permeability of the rock matrix is 3×10-15 m2. The program based on the Monte Carlo Simulation technique [37] is used for the generation, and the basic parameters are input. In order to eliminate the boundary effects, ten fracture networks which are large enough are first generated in the program, and then smaller fracture networks are cut by the boundaries of the “analysis domain,” as shown in Figure 16. From each of the ten large “analysis domains,” fracture network models are extracted with various sizes from 1 m × 1 m to 10 m × 10 m for calculation. In total, 100 models are generated and simulated. To justify whether the calculated intrinsic permeability has a tensor quality at a specific REV size, the fracture network models are rotated at an interval of 15° in the clockwise direction for the calculation. Figure 16 only shows the angles of 0°, 30°, 60°, and 90° models and totally 72 models are generated for this purpose.Table 1
Fracture parameters used for the generation of the fracture network.
Set Density (m-2) Dip (°) Length (m) Aperture (μm) Distribution Standard deviation Mean Distribution Standard deviation Mean 1 0.0231 Normal 2.08 108.39 Log-normal 2.22 14.51 85 2 0.0361 Normal 2.97 142.05 Log-normal 1.56 10.67 85 3 0.0247 Normal 2.58 151.91 Log-normal 2.26 13.8 85 4 0.0231 Normal 2.08 71.61 Log-normal 2.22 14.51 85 5 0.0361 Normal 2.97 37.95 Log-normal 1.56 10.67 85 6 0.0247 Normal 2.58 28.09 Log-normal 2.26 13.8 85Figure 16
Generation of the fracture network and the rotated models.The generated fracture network is imported into ABAQUS for modeling. The element type of the fracture is COH2D4P, and the rock matrix is CPE6MP, which is a triangle element to facilitate the process. A precise model in the side length of 2 m and its meshes are shown in Figure17.Figure 17
A certain fracture network model inside the length of 2 m and its meshes.


### 3.2. Calculation of the Equivalent Intrinsic Permeability

#### 3.2.1. The Darcy Flow

When calculating the fractured rock flow, one assumption is that the fracture is the parallel smooth plane, which means the roughness of the breach is ignored. Figure18 shows the general boundary conditions for the calculation of two types of constant pressure and two types of linear pressures in each direction. Then the flow rates in the x-direction and y-direction can be calculated, and the complete components of the 2D intrinsic permeability tensors of the models, kxx, kyy, kxy, and kyx can be obtained using the following:(14)Qi=Akijμ∂P∂xj,where Qi is the flow rate crossing the boundary along with i-axis, A is the cross-section of the model, kij is the intrinsic permeability tensor, and  ∂P/∂xj is pressure gradient along the j direction.Figure 18
Hydraulic boundary conditions for calculation (P1 and P2 represent different hydraulic pressure).The values of1/kα of each direction would be plotted on a polar diagram. If the values of 1/kα can be fitted into an ellipse, then an intrinsic permeability tensor at a given side length scale may exist. The purpose of checking the possible tensor quality, such as an average intrinsic permeability matrix kij-, was calculated at this range using the following (take α=0 for reference axes):(15)kij-=1N∑1Nkpqααipαjq,where N is the number of rotations, αip and αjq are the directional cosines, and kpqα is the calculated intrinsic permeability in each rotated model.


#### 3.2.2. The Non-Darcy Flow

Before calculating the Darcy flow in a fractured network, the hydraulic conductivity of each fracture needs to be calculated by the Cubic law based on the equivalent hydraulic aperture (equal to the aperture in a smooth parallel plate model) of each fracture, and then the inherent permeability tensork can be calculated. Like the Darcy flow calculation, before calculating the non-Darcy flow, it is also necessary to determine the non-Darcy flow quadratic coefficient B (the cubic law determines the first-order coefficient A) and then the intrinsic permeability tensor k, and non-Darcy coefficient tensor β of the entire fracture network can be determined. With the quadratic coefficient B of a single fracture and the non-Darcy coefficient β, it can be identified by using empirical equations or experiments [15, 24, 38, 39]. Since all the splits in this paper are assumed to be smooth parallel plates with an aperture of 85 um (Table 1), all single breaches in the fracture network have the same coefficients A and B, and therefore they have the same Forchheimer’s equation. As this simple fracture model has a fixed aperture, to obtain the quadratic coefficient B and the non-Darcy coefficient β accurately, the method of using numerical experiments to determine the parameters is adopted in this paper.In this paper, CFD calculation software FLUENT is used for nonlinear simulation of the smooth parallel plate model with lengthL = 100 mm and aperture=85 μm. Non-Darcy coefficient β is obtained when the aperture e=85 μm, which is used to determine the REV size in the case of the non-Darcy flow. In the study, a rectangular grid is utilized to discretize the model. The unit length is 0.1 mm and the aperture is 8.5 um. The governing equation of hydrodynamics is the Navier-Stokes equation. For this paper, the nonlinear partial differential equations of the mass conservation equation (continuity equation) and the momentum conservation equation (momentum equation) of isothermal and incompressible Newtonian fluids can be written as(16)∇v=0(17)∂ρv∂t+ρv∇v-∇μ∇v=-∇P,where t[T] is time. Due to the assumption that effects of gravity are neglected, the physical term is ignored in (17). The turbulence model is calculated using a Realizable k-ε model. The left side of the model is set as the pressure inlet boundary, while the right side of the model is set as the pressure outlet boundary, and the exit boundary is set to zero pressure. The upper and lower boundaries of the fracture model are set as boundary conditions with no fluid flow and no slip. Water density and dynamic viscosity take the corresponding value of 20°C. Results are shown in Figure 19. Fit the numerical data with (3), and B=1.1869E13 kg/m−8 can be obtained with R2=0.999. Thus, (5) can be used to calculate β, and β=85.75 m−1 and all the non-Darcy flow parameters are determined.Figure 19
Numerical simulation results of the smooth parallel fracture.For the non-Darcy flow of the fractured rock mass, consistent with the non-Darcy flow of a single fracture, the slope between the flow rate and the pressure gradient decreases as the pressure gradient increases. Therefore, to determine the REV and non-Darcy coefficient tensors in the non-Darcy flow of the fractured rock mass, different hydraulic boundaries need to be applied in the same model. According to the research of Chen et al. [15] and Liu et al. [28], five different pressure gradients are applied (∂P/∂xj=1×105, 2×105, 3×105, 4×105, and 5×105). Although the five pressure gradients can ensure that the flow regime in the fracture network is the non-Darcy flow in general, due to the complexity of the fracture network, different single fractures in the fracture network may have different flow regimes, and especially the connectivity of some branch fractures is poor (such as the fracture in the red circle in Figure 17). Since a section is in the matrix of the rock mass, the flow regime in the fracture should be the Darcy state. However, many scholars [40–44] think that (3) could probably be applied over the entire range of the flow rate/velocity, including that it reduces linear Darcy’s flow at a sufficiently low flow rate/velocity. Therefore, the flow regime of all fractures in the fracture network can be described by Forchheimer’s equation.Ignoring the “preferential flow” effect [45], totally 500 models are generated. Then (18) is adopted to fit the simulation data, and, therefore, the equivalent intrinsic permeability tensor kij and non-Darcy coefficient tensor βij can be obtained.(18)∂P∂xi=μkijAQj+βijρA2QQj,where βij is the non-Darcy coefficient tensor; Q is a modulus of the flow rate throughout the entire region. The process of solving a model is as follows. Take the model and the boundary shown in Figure 18 as an example. According to the boundary shown in Figure 18, similar to the expansion of Darcy’s equation (14), (18) can be expanded into four linear equations, as shown in the following:(19)∂P1x∂x=μkxxAQ1x+μkxyAQ1y+βxxρA2Q1Q1x+βxyρA2Q1Q1y0=μkyxAQ1x+μkyyAQ1y+βyxρA2Q1Q1x+βyyρA2Q1Q1y0=μkxxAQ2x+μkxyAQ2y+βxxρA2Q2Q2x+βxyρA2Q2Q2y∂P1y∂y=μkyxAQ2x+μkyyAQ2y+βyxρA2Q2Q2x+βyyρA2Q2Q2y,where Q1x is the flow rate of x-direction; Q1y is the flow rate of y-direction; then Q1=Q1x2+Q1y2, the pressure gradient of x-direction is ∂P1x/∂x=P1-P2/L; the pressure gradient of y-direction is ∂P1y/∂y=0; the first two equations of (19) can be obtained; when there is pressure difference in y-direction, the flow rate of x-direction is Q2x. When Q2=Q2x2+Q2y2, the pressure gradient of x-direction is ∂P1x/∂x=0; the pressure gradient of y-direction is ∂P1y/∂y=P1-P2/L; then the last 2 equations are obtained. In the simulation of two-dimensional flow, the side length of the square is L, then A=L. Thus, four equations can be obtained under a pressure gradient, but there are 4 inherent permeability tensors k and 4 non-Darcy coefficient tensor β. In this paper, the fluid flow under five pressure gradients was simulated, forming 45 linear equations. Solve the overdetermined system of 8 coefficients, and then solve the overdetermined equations by the least square method, and the tensor of intrinsic permeability k and non-Darcy coefficient β of the model of the determined size and the rotation angle can be calculated.Bachmat [46], Hsu and Cheng [47], and Wang et al. [48] have proved and applied the theory of the tensor of non-Darcy coefficient β. For the determination of the equivalent intrinsic permeability tensor and the intrinsic permeability ellipse, the same method in the Darcy flow is carried out. The only difference is that, as for the non-Darcy flow, the non-Darcy coefficient tensor also needs to be determined. First, the values of 1/βα of each direction would be plotted on a polar diagram. Then (20) is used to calculate the average non-Darcy coefficient βij- (take the α=0 for reference axes):(20)βij-=1N∑1Nβpqααipαjq,where βpqα is the calculated non-Darcy coefficient in each rotated model.


## 3.1. Generation of the Stochastic Fracture Network

Table1 shows the necessary parameters of the fracture network for this study, and the intrinsic permeability of the rock matrix is 3×10-15 m2. The program based on the Monte Carlo Simulation technique [37] is used for the generation, and the basic parameters are input. In order to eliminate the boundary effects, ten fracture networks which are large enough are first generated in the program, and then smaller fracture networks are cut by the boundaries of the “analysis domain,” as shown in Figure 16. From each of the ten large “analysis domains,” fracture network models are extracted with various sizes from 1 m × 1 m to 10 m × 10 m for calculation. In total, 100 models are generated and simulated. To justify whether the calculated intrinsic permeability has a tensor quality at a specific REV size, the fracture network models are rotated at an interval of 15° in the clockwise direction for the calculation. Figure 16 only shows the angles of 0°, 30°, 60°, and 90° models and totally 72 models are generated for this purpose.Table 1
Fracture parameters used for the generation of the fracture network.
Set Density (m-2) Dip (°) Length (m) Aperture (μm) Distribution Standard deviation Mean Distribution Standard deviation Mean 1 0.0231 Normal 2.08 108.39 Log-normal 2.22 14.51 85 2 0.0361 Normal 2.97 142.05 Log-normal 1.56 10.67 85 3 0.0247 Normal 2.58 151.91 Log-normal 2.26 13.8 85 4 0.0231 Normal 2.08 71.61 Log-normal 2.22 14.51 85 5 0.0361 Normal 2.97 37.95 Log-normal 1.56 10.67 85 6 0.0247 Normal 2.58 28.09 Log-normal 2.26 13.8 85Figure 16
Generation of the fracture network and the rotated models.The generated fracture network is imported into ABAQUS for modeling. The element type of the fracture is COH2D4P, and the rock matrix is CPE6MP, which is a triangle element to facilitate the process. A precise model in the side length of 2 m and its meshes are shown in Figure17.Figure 17
A certain fracture network model inside the length of 2 m and its meshes.


## 3.2. Calculation of the Equivalent Intrinsic Permeability

### 3.2.1. The Darcy Flow

When calculating the fractured rock flow, one assumption is that the fracture is the parallel smooth plane, which means the roughness of the breach is ignored. Figure18 shows the general boundary conditions for the calculation of two types of constant pressure and two types of linear pressures in each direction. Then the flow rates in the x-direction and y-direction can be calculated, and the complete components of the 2D intrinsic permeability tensors of the models, kxx, kyy, kxy, and kyx can be obtained using the following:(14)Qi=Akijμ∂P∂xj,where Qi is the flow rate crossing the boundary along with i-axis, A is the cross-section of the model, kij is the intrinsic permeability tensor, and  ∂P/∂xj is pressure gradient along the j direction.Figure 18
Hydraulic boundary conditions for calculation (P1 and P2 represent different hydraulic pressure).The values of1/kα of each direction would be plotted on a polar diagram. If the values of 1/kα can be fitted into an ellipse, then an intrinsic permeability tensor at a given side length scale may exist. The purpose of checking the possible tensor quality, such as an average intrinsic permeability matrix kij-, was calculated at this range using the following (take α=0 for reference axes):(15)kij-=1N∑1Nkpqααipαjq,where N is the number of rotations, αip and αjq are the directional cosines, and kpqα is the calculated intrinsic permeability in each rotated model.


### 3.2.2. The Non-Darcy Flow

Before calculating the Darcy flow in a fractured network, the hydraulic conductivity of each fracture needs to be calculated by the Cubic law based on the equivalent hydraulic aperture (equal to the aperture in a smooth parallel plate model) of each fracture, and then the inherent permeability tensork can be calculated. Like the Darcy flow calculation, before calculating the non-Darcy flow, it is also necessary to determine the non-Darcy flow quadratic coefficient B (the cubic law determines the first-order coefficient A) and then the intrinsic permeability tensor k, and non-Darcy coefficient tensor β of the entire fracture network can be determined. With the quadratic coefficient B of a single fracture and the non-Darcy coefficient β, it can be identified by using empirical equations or experiments [15, 24, 38, 39]. Since all the splits in this paper are assumed to be smooth parallel plates with an aperture of 85 um (Table 1), all single breaches in the fracture network have the same coefficients A and B, and therefore they have the same Forchheimer’s equation. As this simple fracture model has a fixed aperture, to obtain the quadratic coefficient B and the non-Darcy coefficient β accurately, the method of using numerical experiments to determine the parameters is adopted in this paper.In this paper, CFD calculation software FLUENT is used for nonlinear simulation of the smooth parallel plate model with lengthL = 100 mm and aperture=85 μm. Non-Darcy coefficient β is obtained when the aperture e=85 μm, which is used to determine the REV size in the case of the non-Darcy flow. In the study, a rectangular grid is utilized to discretize the model. The unit length is 0.1 mm and the aperture is 8.5 um. The governing equation of hydrodynamics is the Navier-Stokes equation. For this paper, the nonlinear partial differential equations of the mass conservation equation (continuity equation) and the momentum conservation equation (momentum equation) of isothermal and incompressible Newtonian fluids can be written as(16)∇v=0(17)∂ρv∂t+ρv∇v-∇μ∇v=-∇P,where t[T] is time. Due to the assumption that effects of gravity are neglected, the physical term is ignored in (17). The turbulence model is calculated using a Realizable k-ε model. The left side of the model is set as the pressure inlet boundary, while the right side of the model is set as the pressure outlet boundary, and the exit boundary is set to zero pressure. The upper and lower boundaries of the fracture model are set as boundary conditions with no fluid flow and no slip. Water density and dynamic viscosity take the corresponding value of 20°C. Results are shown in Figure 19. Fit the numerical data with (3), and B=1.1869E13 kg/m−8 can be obtained with R2=0.999. Thus, (5) can be used to calculate β, and β=85.75 m−1 and all the non-Darcy flow parameters are determined.Figure 19
Numerical simulation results of the smooth parallel fracture.For the non-Darcy flow of the fractured rock mass, consistent with the non-Darcy flow of a single fracture, the slope between the flow rate and the pressure gradient decreases as the pressure gradient increases. Therefore, to determine the REV and non-Darcy coefficient tensors in the non-Darcy flow of the fractured rock mass, different hydraulic boundaries need to be applied in the same model. According to the research of Chen et al. [15] and Liu et al. [28], five different pressure gradients are applied (∂P/∂xj=1×105, 2×105, 3×105, 4×105, and 5×105). Although the five pressure gradients can ensure that the flow regime in the fracture network is the non-Darcy flow in general, due to the complexity of the fracture network, different single fractures in the fracture network may have different flow regimes, and especially the connectivity of some branch fractures is poor (such as the fracture in the red circle in Figure 17). Since a section is in the matrix of the rock mass, the flow regime in the fracture should be the Darcy state. However, many scholars [40–44] think that (3) could probably be applied over the entire range of the flow rate/velocity, including that it reduces linear Darcy’s flow at a sufficiently low flow rate/velocity. Therefore, the flow regime of all fractures in the fracture network can be described by Forchheimer’s equation.Ignoring the “preferential flow” effect [45], totally 500 models are generated. Then (18) is adopted to fit the simulation data, and, therefore, the equivalent intrinsic permeability tensor kij and non-Darcy coefficient tensor βij can be obtained.(18)∂P∂xi=μkijAQj+βijρA2QQj,where βij is the non-Darcy coefficient tensor; Q is a modulus of the flow rate throughout the entire region. The process of solving a model is as follows. Take the model and the boundary shown in Figure 18 as an example. According to the boundary shown in Figure 18, similar to the expansion of Darcy’s equation (14), (18) can be expanded into four linear equations, as shown in the following:(19)∂P1x∂x=μkxxAQ1x+μkxyAQ1y+βxxρA2Q1Q1x+βxyρA2Q1Q1y0=μkyxAQ1x+μkyyAQ1y+βyxρA2Q1Q1x+βyyρA2Q1Q1y0=μkxxAQ2x+μkxyAQ2y+βxxρA2Q2Q2x+βxyρA2Q2Q2y∂P1y∂y=μkyxAQ2x+μkyyAQ2y+βyxρA2Q2Q2x+βyyρA2Q2Q2y,where Q1x is the flow rate of x-direction; Q1y is the flow rate of y-direction; then Q1=Q1x2+Q1y2, the pressure gradient of x-direction is ∂P1x/∂x=P1-P2/L; the pressure gradient of y-direction is ∂P1y/∂y=0; the first two equations of (19) can be obtained; when there is pressure difference in y-direction, the flow rate of x-direction is Q2x. When Q2=Q2x2+Q2y2, the pressure gradient of x-direction is ∂P1x/∂x=0; the pressure gradient of y-direction is ∂P1y/∂y=P1-P2/L; then the last 2 equations are obtained. In the simulation of two-dimensional flow, the side length of the square is L, then A=L. Thus, four equations can be obtained under a pressure gradient, but there are 4 inherent permeability tensors k and 4 non-Darcy coefficient tensor β. In this paper, the fluid flow under five pressure gradients was simulated, forming 45 linear equations. Solve the overdetermined system of 8 coefficients, and then solve the overdetermined equations by the least square method, and the tensor of intrinsic permeability k and non-Darcy coefficient β of the model of the determined size and the rotation angle can be calculated.Bachmat [46], Hsu and Cheng [47], and Wang et al. [48] have proved and applied the theory of the tensor of non-Darcy coefficient β. For the determination of the equivalent intrinsic permeability tensor and the intrinsic permeability ellipse, the same method in the Darcy flow is carried out. The only difference is that, as for the non-Darcy flow, the non-Darcy coefficient tensor also needs to be determined. First, the values of 1/βα of each direction would be plotted on a polar diagram. Then (20) is used to calculate the average non-Darcy coefficient βij- (take the α=0 for reference axes):(20)βij-=1N∑1Nβpqααipαjq,where βpqα is the calculated non-Darcy coefficient in each rotated model.


## 3.2.1. The Darcy Flow

When calculating the fractured rock flow, one assumption is that the fracture is the parallel smooth plane, which means the roughness of the breach is ignored. Figure18 shows the general boundary conditions for the calculation of two types of constant pressure and two types of linear pressures in each direction. Then the flow rates in the x-direction and y-direction can be calculated, and the complete components of the 2D intrinsic permeability tensors of the models, kxx, kyy, kxy, and kyx can be obtained using the following:(14)Qi=Akijμ∂P∂xj,where Qi is the flow rate crossing the boundary along with i-axis, A is the cross-section of the model, kij is the intrinsic permeability tensor, and  ∂P/∂xj is pressure gradient along the j direction.Figure 18
Hydraulic boundary conditions for calculation (P1 and P2 represent different hydraulic pressure).The values of1/kα of each direction would be plotted on a polar diagram. If the values of 1/kα can be fitted into an ellipse, then an intrinsic permeability tensor at a given side length scale may exist. The purpose of checking the possible tensor quality, such as an average intrinsic permeability matrix kij-, was calculated at this range using the following (take α=0 for reference axes):(15)kij-=1N∑1Nkpqααipαjq,where N is the number of rotations, αip and αjq are the directional cosines, and kpqα is the calculated intrinsic permeability in each rotated model.


## 3.2.2. The Non-Darcy Flow

Before calculating the Darcy flow in a fractured network, the hydraulic conductivity of each fracture needs to be calculated by the Cubic law based on the equivalent hydraulic aperture (equal to the aperture in a smooth parallel plate model) of each fracture, and then the inherent permeability tensork can be calculated. Like the Darcy flow calculation, before calculating the non-Darcy flow, it is also necessary to determine the non-Darcy flow quadratic coefficient B (the cubic law determines the first-order coefficient A) and then the intrinsic permeability tensor k, and non-Darcy coefficient tensor β of the entire fracture network can be determined. With the quadratic coefficient B of a single fracture and the non-Darcy coefficient β, it can be identified by using empirical equations or experiments [15, 24, 38, 39]. Since all the splits in this paper are assumed to be smooth parallel plates with an aperture of 85 um (Table 1), all single breaches in the fracture network have the same coefficients A and B, and therefore they have the same Forchheimer’s equation. As this simple fracture model has a fixed aperture, to obtain the quadratic coefficient B and the non-Darcy coefficient β accurately, the method of using numerical experiments to determine the parameters is adopted in this paper.In this paper, CFD calculation software FLUENT is used for nonlinear simulation of the smooth parallel plate model with lengthL = 100 mm and aperture=85 μm. Non-Darcy coefficient β is obtained when the aperture e=85 μm, which is used to determine the REV size in the case of the non-Darcy flow. In the study, a rectangular grid is utilized to discretize the model. The unit length is 0.1 mm and the aperture is 8.5 um. The governing equation of hydrodynamics is the Navier-Stokes equation. For this paper, the nonlinear partial differential equations of the mass conservation equation (continuity equation) and the momentum conservation equation (momentum equation) of isothermal and incompressible Newtonian fluids can be written as(16)∇v=0(17)∂ρv∂t+ρv∇v-∇μ∇v=-∇P,where t[T] is time. Due to the assumption that effects of gravity are neglected, the physical term is ignored in (17). The turbulence model is calculated using a Realizable k-ε model. The left side of the model is set as the pressure inlet boundary, while the right side of the model is set as the pressure outlet boundary, and the exit boundary is set to zero pressure. The upper and lower boundaries of the fracture model are set as boundary conditions with no fluid flow and no slip. Water density and dynamic viscosity take the corresponding value of 20°C. Results are shown in Figure 19. Fit the numerical data with (3), and B=1.1869E13 kg/m−8 can be obtained with R2=0.999. Thus, (5) can be used to calculate β, and β=85.75 m−1 and all the non-Darcy flow parameters are determined.Figure 19
Numerical simulation results of the smooth parallel fracture.For the non-Darcy flow of the fractured rock mass, consistent with the non-Darcy flow of a single fracture, the slope between the flow rate and the pressure gradient decreases as the pressure gradient increases. Therefore, to determine the REV and non-Darcy coefficient tensors in the non-Darcy flow of the fractured rock mass, different hydraulic boundaries need to be applied in the same model. According to the research of Chen et al. [15] and Liu et al. [28], five different pressure gradients are applied (∂P/∂xj=1×105, 2×105, 3×105, 4×105, and 5×105). Although the five pressure gradients can ensure that the flow regime in the fracture network is the non-Darcy flow in general, due to the complexity of the fracture network, different single fractures in the fracture network may have different flow regimes, and especially the connectivity of some branch fractures is poor (such as the fracture in the red circle in Figure 17). Since a section is in the matrix of the rock mass, the flow regime in the fracture should be the Darcy state. However, many scholars [40–44] think that (3) could probably be applied over the entire range of the flow rate/velocity, including that it reduces linear Darcy’s flow at a sufficiently low flow rate/velocity. Therefore, the flow regime of all fractures in the fracture network can be described by Forchheimer’s equation.Ignoring the “preferential flow” effect [45], totally 500 models are generated. Then (18) is adopted to fit the simulation data, and, therefore, the equivalent intrinsic permeability tensor kij and non-Darcy coefficient tensor βij can be obtained.(18)∂P∂xi=μkijAQj+βijρA2QQj,where βij is the non-Darcy coefficient tensor; Q is a modulus of the flow rate throughout the entire region. The process of solving a model is as follows. Take the model and the boundary shown in Figure 18 as an example. According to the boundary shown in Figure 18, similar to the expansion of Darcy’s equation (14), (18) can be expanded into four linear equations, as shown in the following:(19)∂P1x∂x=μkxxAQ1x+μkxyAQ1y+βxxρA2Q1Q1x+βxyρA2Q1Q1y0=μkyxAQ1x+μkyyAQ1y+βyxρA2Q1Q1x+βyyρA2Q1Q1y0=μkxxAQ2x+μkxyAQ2y+βxxρA2Q2Q2x+βxyρA2Q2Q2y∂P1y∂y=μkyxAQ2x+μkyyAQ2y+βyxρA2Q2Q2x+βyyρA2Q2Q2y,where Q1x is the flow rate of x-direction; Q1y is the flow rate of y-direction; then Q1=Q1x2+Q1y2, the pressure gradient of x-direction is ∂P1x/∂x=P1-P2/L; the pressure gradient of y-direction is ∂P1y/∂y=0; the first two equations of (19) can be obtained; when there is pressure difference in y-direction, the flow rate of x-direction is Q2x. When Q2=Q2x2+Q2y2, the pressure gradient of x-direction is ∂P1x/∂x=0; the pressure gradient of y-direction is ∂P1y/∂y=P1-P2/L; then the last 2 equations are obtained. In the simulation of two-dimensional flow, the side length of the square is L, then A=L. Thus, four equations can be obtained under a pressure gradient, but there are 4 inherent permeability tensors k and 4 non-Darcy coefficient tensor β. In this paper, the fluid flow under five pressure gradients was simulated, forming 45 linear equations. Solve the overdetermined system of 8 coefficients, and then solve the overdetermined equations by the least square method, and the tensor of intrinsic permeability k and non-Darcy coefficient β of the model of the determined size and the rotation angle can be calculated.Bachmat [46], Hsu and Cheng [47], and Wang et al. [48] have proved and applied the theory of the tensor of non-Darcy coefficient β. For the determination of the equivalent intrinsic permeability tensor and the intrinsic permeability ellipse, the same method in the Darcy flow is carried out. The only difference is that, as for the non-Darcy flow, the non-Darcy coefficient tensor also needs to be determined. First, the values of 1/βα of each direction would be plotted on a polar diagram. Then (20) is used to calculate the average non-Darcy coefficient βij- (take the α=0 for reference axes):(20)βij-=1N∑1Nβpqααipαjq,where βpqα is the calculated non-Darcy coefficient in each rotated model.


## 4. Determination of the REV

In this section, firstly, the numerical results of the Darcy flow about size effect on the intrinsic permeability are presented. Then the statistical results of the non-Darcy flow are presented with size effect on inherent permeability and the non-Darcy coefficient. The graphic presentation of results is shown more clearly, with direct results of intrinsic permeability components (kxx, kyy, kxy, and kyx) and non-Darcy coefficient components (βxx, βyy, βxy, and βyx), including their simulation results, mean, maximum, and minimum. The coefficient of variation (CV) is used for the determination of the REV. For comparison of the REV scale, the values of “CV” with the increase of the model size are also plotted on the graph. Finally, for the evaluation of the intrinsic permeability tensor and the non-Darcy coefficient tensor, calculated and fitted ellipses are drawn at a certain model size.


### 4.1. Results of the Darcy Flow

Figure20 shows the pore pressure field at the model side length of 9 m and the rotation angle of 0° in the Darcy flow, from which the pore pressure variation of the flow in fractured rock can be obviously seen. Calculated values of the minimum, maximum, and mean of the intrinsic permeability components of 10 fracture network realizations of different sizes in the Darcy flow are shown in Figure 21. The results show that the range of variation of each intrinsic permeability components, kxx, kyy, kxy, and kyx, becomes smaller with the increase of the model side length, which indicates the existence of the REV. The “CV” values of the diagonal components of the intrinsic permeability tensor kxx and kyy are shown in Figure 22, which shows the values of “CV” both larger than 30% at the side length of 1 m, and it means that the difference of each model is somewhat big. As the side length of 6 m is exceeded, kxx and kyy never change significantly. The “CV” values are both below 20%. In addition to that, the “CV” values are both below 10% with the side length of 8 m. Besides, the values of kxy and kyx are almost close to zero with the increase in the model size, which indicates that a nearly symmetric intrinsic permeability matrix would present. Comparing with the ellipse at the side length of 1 m, an ideal ellipse can be reached at the side length of 8 m (as shown in Figure 27, the dashed line represents average 1/k, and the solid line represents the calculated 1/k).Figure 20
The pore pressure field at the side length of 9 m and the rotation angle of 0°.Figure 21
The results, Max, Min, and mean of calculated intrinsic permeability components in the Darcy flow. (a)kxx, (b) kyy, (c) kxy, and  kyx.
(a) (b) (c)Figure 22
The values of coefficients of variation ofkxx  and  kyy with increasing side lengths in the Darcy flow.


### 4.2. Results of the Non-Darcy Flow

Figure23 shows the curves of calculated results, Min, Max, and mean of intrinsic permeability components, kxx, kyy, kxy, and kyx, which are scattered when the side length is small. With the increase in the model side length, the intrinsic permeability components tend to be steady, indicating the existence of the REV in the non-Darcy flow. The non-Darcy coefficient component (βxx, βyy, βxy, and βyx) variation curves are shown in Figure 24. Like the intrinsic permeability component curves, the non-Darcy coefficient component curves tend to be steady with the increase in the side length.Figure 23
The results, Max, Min, and mean of calculated intrinsic permeability components in the non-Darcy flow. (a)kxx, (b) kyy, (c) kxy,  and  kyx.
(a) (b) (c)Figure 24
The results, Max, Min, and mean of calculated non-Darcy coefficient components in the non-Darcy flow. (a)βxx, (b) βyy, (c) βxy,  and  βyx.
(a) (b) (c)Figures25 and 26 show the values of “CV” of kxx, kyy, βxx, and βyy with the increase in the model side length. As for kxx and kyy, at the side length of 1 m, the values of “CV” are beyond 30%. The values of “CV” tend to decrease with the increase of the model side length and are below 20% at the side length of 6 m. At the side length of 8 m, the values of “CV” are below 10%. As for βxx and βyy, the values of “CV” are both beyond 40% at the side length of 8 m. With the increase of the model side length, the values of “CV” decrease and are below 20% at the side length of 7 m. The values of “CV” are below 10% at the side length of 9 m.Figure 25
The values of coefficients of variation ofkxx  and  kyy with increasing side lengths in the non-Darcy flow.Figure 26
The values of coefficients of variation ofβxx  and  βyy with increasing side lengths in the non-Darcy flow.Figure 27
The equivalent intrinsic permeability ellipse with the increasing model size in the Darcy flow.The intrinsic permeability tensor ellipses of the non-Darcy flow are drawn in Figure28, from which it can be found that the tendency of approaching an ellipse by the simulation data becomes stronger as the side length of the model increases. As for the non-Darcy coefficient tensor ellipse, an ideal ellipse can be achieved at the side length of 7 m as shown in Figure 29, which means that the non-Darcy coefficient has the quality of the tensor.Figure 28
The equivalent intrinsic permeability ellipse with the increasing model size in the non-Darcy flow.Figure 29
The equivalent non-Darcy coefficient ellipse with the increasing model size in the non-Darcy flow.


### 4.3. Scale Effect and the REV Size

According to the calculated results above, the following conclusions are mainly drawn. (I) As for the Darcy flow, if an acceptable value of “CV” of 20% is set, then the side length of 6 m can be approximated for the REV size. If the value of “CV” of 10% is accepted, then the side length of 8 m can be approximated by the REV size. (II) As for the non-Darcy flow, two parameters (the intrinsic permeabilityk and the non-Darcy coefficient β) are necessary to determine the REV size. For the intrinsic permeability k, if the value of “CV” of 20% is accepted, then the side length of 6 m can be approximated by the REV size. Besides, if the value of “CV” of 10% is accepted, then the side length of 8 m can be approximated by the REV size. For the non-Darcy coefficient β, if the value of “CV” of 20% is accepted, then the side length of 7 m can be approximated by the REV size, and if the value of “CV” of 10% is accepted, then the side length of 9 m can be approximated by the REV size.Besides, it is worth noting that, for the flow of Darcy and non-Darcy, the variation of the CV value of the equivalent permeability coefficientk is consistent, or, in other words, both flow regimes share the same REV. When the shapes of the intrinsic permeability ellipse are compared, they are almost identical. Thus, it is proved that the intrinsic permeability k is a physical quantity that reflects the inherent properties of fractured rock and has nothing to do with flow regimes.Above all, what we can conclude is that if the value of “CV” of 20% is accepted,REVDarcy=6 m and REVnon-Darcy=7 m. If the value of “CV” of 10% is accepted, REVDarcy=8 m and REVnon-Darcy=9 m. Therefore, the REV size of the non-Darcy flow is a little larger than the REV size of the Darcy flow.


## 4.1. Results of the Darcy Flow

Figure20 shows the pore pressure field at the model side length of 9 m and the rotation angle of 0° in the Darcy flow, from which the pore pressure variation of the flow in fractured rock can be obviously seen. Calculated values of the minimum, maximum, and mean of the intrinsic permeability components of 10 fracture network realizations of different sizes in the Darcy flow are shown in Figure 21. The results show that the range of variation of each intrinsic permeability components, kxx, kyy, kxy, and kyx, becomes smaller with the increase of the model side length, which indicates the existence of the REV. The “CV” values of the diagonal components of the intrinsic permeability tensor kxx and kyy are shown in Figure 22, which shows the values of “CV” both larger than 30% at the side length of 1 m, and it means that the difference of each model is somewhat big. As the side length of 6 m is exceeded, kxx and kyy never change significantly. The “CV” values are both below 20%. In addition to that, the “CV” values are both below 10% with the side length of 8 m. Besides, the values of kxy and kyx are almost close to zero with the increase in the model size, which indicates that a nearly symmetric intrinsic permeability matrix would present. Comparing with the ellipse at the side length of 1 m, an ideal ellipse can be reached at the side length of 8 m (as shown in Figure 27, the dashed line represents average 1/k, and the solid line represents the calculated 1/k).Figure 20
The pore pressure field at the side length of 9 m and the rotation angle of 0°.Figure 21
The results, Max, Min, and mean of calculated intrinsic permeability components in the Darcy flow. (a)kxx, (b) kyy, (c) kxy, and  kyx.
(a) (b) (c)Figure 22
The values of coefficients of variation ofkxx  and  kyy with increasing side lengths in the Darcy flow.


## 4.2. Results of the Non-Darcy Flow

Figure23 shows the curves of calculated results, Min, Max, and mean of intrinsic permeability components, kxx, kyy, kxy, and kyx, which are scattered when the side length is small. With the increase in the model side length, the intrinsic permeability components tend to be steady, indicating the existence of the REV in the non-Darcy flow. The non-Darcy coefficient component (βxx, βyy, βxy, and βyx) variation curves are shown in Figure 24. Like the intrinsic permeability component curves, the non-Darcy coefficient component curves tend to be steady with the increase in the side length.Figure 23
The results, Max, Min, and mean of calculated intrinsic permeability components in the non-Darcy flow. (a)kxx, (b) kyy, (c) kxy,  and  kyx.
(a) (b) (c)Figure 24
The results, Max, Min, and mean of calculated non-Darcy coefficient components in the non-Darcy flow. (a)βxx, (b) βyy, (c) βxy,  and  βyx.
(a) (b) (c)Figures25 and 26 show the values of “CV” of kxx, kyy, βxx, and βyy with the increase in the model side length. As for kxx and kyy, at the side length of 1 m, the values of “CV” are beyond 30%. The values of “CV” tend to decrease with the increase of the model side length and are below 20% at the side length of 6 m. At the side length of 8 m, the values of “CV” are below 10%. As for βxx and βyy, the values of “CV” are both beyond 40% at the side length of 8 m. With the increase of the model side length, the values of “CV” decrease and are below 20% at the side length of 7 m. The values of “CV” are below 10% at the side length of 9 m.Figure 25
The values of coefficients of variation ofkxx  and  kyy with increasing side lengths in the non-Darcy flow.Figure 26
The values of coefficients of variation ofβxx  and  βyy with increasing side lengths in the non-Darcy flow.Figure 27
The equivalent intrinsic permeability ellipse with the increasing model size in the Darcy flow.The intrinsic permeability tensor ellipses of the non-Darcy flow are drawn in Figure28, from which it can be found that the tendency of approaching an ellipse by the simulation data becomes stronger as the side length of the model increases. As for the non-Darcy coefficient tensor ellipse, an ideal ellipse can be achieved at the side length of 7 m as shown in Figure 29, which means that the non-Darcy coefficient has the quality of the tensor.Figure 28
The equivalent intrinsic permeability ellipse with the increasing model size in the non-Darcy flow.Figure 29
The equivalent non-Darcy coefficient ellipse with the increasing model size in the non-Darcy flow.


## 4.3. Scale Effect and the REV Size

According to the calculated results above, the following conclusions are mainly drawn. (I) As for the Darcy flow, if an acceptable value of “CV” of 20% is set, then the side length of 6 m can be approximated for the REV size. If the value of “CV” of 10% is accepted, then the side length of 8 m can be approximated by the REV size. (II) As for the non-Darcy flow, two parameters (the intrinsic permeabilityk and the non-Darcy coefficient β) are necessary to determine the REV size. For the intrinsic permeability k, if the value of “CV” of 20% is accepted, then the side length of 6 m can be approximated by the REV size. Besides, if the value of “CV” of 10% is accepted, then the side length of 8 m can be approximated by the REV size. For the non-Darcy coefficient β, if the value of “CV” of 20% is accepted, then the side length of 7 m can be approximated by the REV size, and if the value of “CV” of 10% is accepted, then the side length of 9 m can be approximated by the REV size.Besides, it is worth noting that, for the flow of Darcy and non-Darcy, the variation of the CV value of the equivalent permeability coefficientk is consistent, or, in other words, both flow regimes share the same REV. When the shapes of the intrinsic permeability ellipse are compared, they are almost identical. Thus, it is proved that the intrinsic permeability k is a physical quantity that reflects the inherent properties of fractured rock and has nothing to do with flow regimes.Above all, what we can conclude is that if the value of “CV” of 20% is accepted,REVDarcy=6 m and REVnon-Darcy=7 m. If the value of “CV” of 10% is accepted, REVDarcy=8 m and REVnon-Darcy=9 m. Therefore, the REV size of the non-Darcy flow is a little larger than the REV size of the Darcy flow.


## 5. Discussions and Conclusions

Above all, the pore pressure cohesive element (PPCE) with a subroutine is developed, and the simulation of the non-Darcy flow in the fracture is realized. Then the simulation method of the Darcy and non-Darcy flow in the dual-porosity model, like fractured rock, is built. Main conclusions are as follows.(I) The PPCE can be inserted to simulate the Darcy flow in the dual-porosity model, like fractured rock, in ABAQUS. With the subroutine, the PPCE can also simulate the non-Darcy flow in fractured rock. The calculated results are compared with other existing results of the dual-porosity simulation to validate its correctness and effectiveness. A new and more convenient method is provided to simulate the flow in the dual-porosity model.(II) The Monte Carlo Simulation technique is adopted for the generation of the fractured network. Based on the dual-porosity model with the PPCE, generated fracture networks with different side lengths and rotation angles are simulated in both Darcy and non-Darcy flow. On the one hand, the intrinsic permeabilityk calculated by Darcy’s law and Forchheimer’s law has identical convergence, and the determined REVs with the same criterion are identical. On the other hand, the value of “CV” of non-Darcy coefficient β has poor convergence, indicating that the REVs determined in the Darcy flow and the non-Darcy flow are inconsistent. Thus, the REV determined in the Darcy flow, the traditional way, fails to suit the case of the non-Darcy flow.(III) According to the study on the determination of the REV in the Darcy flow and the non-Darcy flow, the REV determined in the non-Darcy flow is a little larger than that in the Darcy flow. Thus, it is appropriate to propose that the REV should be determined by the non-Darcy coefficientβ, which can satisfy not only the case of the Darcy flow, but also the case of the non-Darcy flow.In this paper, based on the dual-medium model, the Monte Carlo fractured rock mass is studied tentatively using the non-Darcy coefficient to determine the REV value under the non-Darcy flow. The REV of the Darcy and non-Darcy flow is compared. Some new understandings have been found. For example, the REV of the non-Darcy flow is inconsistent with the REV of the Darcy flow, and the REV of the non-Darcy flow is larger than the REV of the Darcy flow.The inconsistency between the REVs mainly results from the REV of the non-Darcy flow determined by non-Darcy coefficientβ, and the REV of the Darcy flow is determined by inherent permeability k, when the non-Darcy coefficient of porous media is considered to be only related to intrinsic permeability k [49–51]. Then the REV determined by non-Darcy coefficient β and the REV determined by inherent permeability k should be consistent, but more studies have shown that the non-Darcy coefficient is not only related to the intrinsic permeability k of media. Evans and Civan [52], Geertsma [53], and Macdonald et al. [54] found that non-Darcy coefficient β is also related to the effective porosity of the pore medium; Li and Engler [55] and Veyskarami et al. [56] proposed that, for isotropic media, several correlations relate β to petrophysical parameters (in particular ϕ,k, and the tortuosity τ); Knupp and Lage [57] also pointed out that β is related to the inherent permeability k and inertia coefficient cF, where cF represents the microform resistance exerted by the solid porous substrate; that is to say, it depends on the geometry of the solids in each basic representative volume of the porous medium. It is noteworthy that cF and k are independent parameters. Therefore, it is not only related to the intrinsic permeability k but also related to other physical properties of the medium, which in turn leads to the inconsistency between the REV determined by β and the REV defined by k, which may also be the reason for the disparity of the REV between the Darcy flow and non-Darcy flow in the fractured rock mass.The conclusion that the REV of the non-Darcy flow is larger than that of the Darcy flow is only based on the results shown in the model calculation in this paper. Whether this conclusion is universal and whether this difference can be quantitatively described and which physical parameters of the fractured rock mass influence this outcome need to be further explored and studied in the future work.Besides, it should also be mentioned that the PPCE has a wide range of applications in flow in fractured rock and soil, and the method proposed in this paper can also be applied in the simulation of on-site engineering. What is more, when it comes to the intersecting fractures in the non-Darcy flow, the “preference flow” is ignored in this study. It is suggested that similar studies involve this factor. For more work on the REV and its intrinsic permeability tensor of the fractured rock mass in the non-Darcy flow, factors like roughness, stress, and aperture also should be taken into consideration.



---

*Source: 7535927-2018-04-22.xml*