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

Wang, Yuan; Niu, Yu-long; Feng, Qiang

doi:https://doi.org/10.1155/2018/7535927

Geofluids

On this page

Abstract Introduction Conclusions Conflicts of Interest Acknowledgments References Copyright Related Articles

Research Article | Open Access

Volume 2018 | Article ID 7535927 | https://doi.org/10.1155/2018/7535927

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

Yuan Wang,¹Yu-long Niu,^1,2and Qiang Feng^1,2

Academic Editor: Cornelius Langenbruch

Received21 Sept 2017

Revised10 Feb 2018

Accepted19 Feb 2018

Published22 Apr 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 permeability 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 tensors obtained in the Darcy flow and the non-Darcy flow are basically the same.

1. Introduction

Hydraulic property of the fractured rock mass plays a vital role in on-site engineering, like the underground basement, oil storage, CO₂ 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)where is the flow velocity; is the intrinsic permeability of porous media; 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:where is the volumetric flow rate or discharge; is the fracture aperture; is the aperture of the idealized parallel fracture, and it is the equivalent hydraulic aperture for the rough fracture; is the dynamic viscosity of the fluid; is the intrinsic permeability defined as ; 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:where is the pressure gradient along the flow direction; and are the coefficients, which describe the energy losses of the flow caused by viscosity and inertia, respectively. and in (3) are commonly written as follows:where 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:where 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.

In ABAQUS, the default constitutive relationship of tangential flow is described as (2), and the constitutive relationship of normal flow is described as follows:where is the fluid leak-off coefficients at the top and bottom element surfaces; and 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 Figure 3, and are porous media and is the fracture. Vector is a vector normal to . In , , for single-phase incompressible fluids, it can be expressed as follows:In the flow exchange of fractures and the matrix, the law of mass conservationwhere , being the tangential flow rate in the fracture; , represent the flow rate into the fissure and out of the fracture, respectively.

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:where is the volumetric flow rate of the node, is hydraulic conductivity of the matrix, and 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 rate can be obtained by solving (4):

Then,

From (12), it is obvious that the hydraulic conductivity 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 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.

As for the non-Darcy flow, the constitutive default of tangential flow can be adapted with a subroutine (nodes 1 and 2 in Figure 1), 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:where and are the flow rate of the node, L is the length of the line element, K is hydraulic conductivity, and and 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 Figure 5, specific steps are as follows. Each increment ( or ) is small enough in the calculation. At the first increment, the quadratic term in (3) is ignored, and then and 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.

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 Figure 6. The length of the fracture is m, m², and m².

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.

(a) ABAQUS-PPCE

(b) DFM

(c) FM-Meshfree

(d) FM-Meshfree (irregular nodal distribution)

(a) ABAQUS-PPCE

(b) DFM

(c) FM-Meshfree

(d) FM-Meshfree (irregular nodal distribution)

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 Figure 10(a). To evaluate the effects of the fracture on free surface flow, a fracture with the aperture and the length 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.

(a) Porous media

(b) Fractured media

(a) Porous media

(b) Fractured media

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 , m (ignoring cubic law), , and . 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.

(a)

(b)

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

Table 1 shows the necessary parameters of the fracture network for this study, and the intrinsic permeability of the rock matrix is m². 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.

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 Figure 17.

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. Figure 18 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 -direction and -direction can be calculated, and the complete components of the 2D intrinsic permeability tensors of the models, , , , and can be obtained using the following:where is the flow rate crossing the boundary along with -axis, is the cross-section of the model, is the intrinsic permeability tensor, and is pressure gradient along the direction.

The values of of each direction would be plotted on a polar diagram. If the values of 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 , was calculated at this range using the following (take for reference axes):where is the number of rotations, and are the directional cosines, and 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 tensor 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 (the cubic law determines the first-order coefficient ) and then the intrinsic permeability tensor , and non-Darcy coefficient tensor of the entire fracture network can be determined. With the quadratic coefficient 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 and , and therefore they have the same Forchheimer’s equation. As this simple fracture model has a fixed aperture, to obtain the quadratic coefficient 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 length L = 100 mm and aperture μm. Non-Darcy coefficient is obtained when the aperture μ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 aswhere 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 - 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 kg/m⁻⁸ can be obtained with . Thus, (5) can be used to calculate , and m⁻¹ and all the non-Darcy flow parameters are determined.

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 (, , , , and ). 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 and non-Darcy coefficient tensor can be obtained.where is the non-Darcy coefficient tensor; 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:where is the flow rate of -direction; is the flow rate of -direction; then , the pressure gradient of -direction is ; the pressure gradient of -direction is ; the first two equations of (19) can be obtained; when there is pressure difference in -direction, the flow rate of -direction is . When , the pressure gradient of -direction is ; the pressure gradient of -direction is ; then the last 2 equations are obtained. In the simulation of two-dimensional flow, the side length of the square is , then . Thus, four equations can be obtained under a pressure gradient, but there are 4 inherent permeability tensors 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 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 of each direction would be plotted on a polar diagram. Then (20) is used to calculate the average non-Darcy coefficient (take the for reference axes):where 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 (, , , and ) and non-Darcy coefficient components (, , , and ), 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

Figure 20 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, , , , and , 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 and 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, and 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 and 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 , and the solid line represents the calculated ).

(a)

(b)

(c)

4.2. Results of the Non-Darcy Flow

Figure 23 shows the curves of calculated results, Min, Max, and mean of intrinsic permeability components, , , , and , 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 (, , , and ) 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.

(a)

(b)

(c)

(a)

(b)

(c)

Figures 25 and 26 show the values of “CV” of , , , and with the increase in the model side length. As for and , 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 and , 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.

The intrinsic permeability tensor ellipses of the non-Darcy flow are drawn in Figure 28, 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.

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 permeability and the non-Darcy coefficient β) are necessary to determine the REV size. For the intrinsic permeability , 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 coefficient 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 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, m and m. If the value of “CV” of 10% is accepted, m and 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 permeability 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 , when the non-Darcy coefficient of porous media is considered to be only related to intrinsic permeability [49–51]. Then the REV determined by non-Darcy coefficient and the REV determined by inherent permeability should be consistent, but more studies have shown that the non-Darcy coefficient is not only related to the intrinsic permeability 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 and inertia coefficient , where 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 and are independent parameters. Therefore, it is not only related to the intrinsic permeability 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 , 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.

Symbols

	Flow velocity
	Intrinsic permeability of porous media
	Pressure gradient
	Volumetric flow rate or discharge
	Fracture aperture [L]
	Aperture of the idealized parallel smooth fracture [L]
	Equivalent hydraulic aperture [L]
	Dynamic viscosity of fluid
	Intrinsic permeability of fracture
	Cross-section area
,	Forchheimer coefficients ,
	Non-Darcy coefficient
	Density of fluid []
	Reynold’s number [-]
	Fluid leak-off coefficients []
,	Pore pressure []
	Flow rate along with -axis
	Cross-section of the model []
	Intrinsic permeability tensor []
,	Directional cosines [-]
	Calculated intrinsic permeability in each rotated model []
,	Constant associated with fractured rock mass [-]
	Non-Darcy coefficient tensor
:	Modulus of the flow rate
	Calculated non-Darcy coefficient in each rotated model
	Coefficient of variation [-].

Conflicts of Interest

The authors declare that there are no conflicts of interest regarding the publication of this paper.

Acknowledgments

The authors gratefully acknowledge the support of the National Key R&D Program of China no. 2016YFC0401804.

References

J. Bear, “Dynamics of Fluids in Porous Media,” Soil Science, vol. 120, no. 2, pp. 162-163, 1975.
View at: Publisher Site | Google Scholar
J. C. S. Long, J. S. Remer, C. R. Wilson, and P. A. Witherspoon, “Porous media equivalents for networks of discontinuous fractures,” Resources Research, vol. 18, no. 3, pp. 645–658, 1982.
View at: Google Scholar
M. Oda, “An equivalent model for coupled stress and fluid flow analysis in jointed rock masse,” Water Resources Research, vol. 22, no. 13, pp. 1845–1856, 1986.
View at: Publisher Site | Google Scholar
M. Oda, “Permeability tensor for discontinuous rock masses,” Géotechnique, vol. 35, no. 4, pp. 483–495, 1985.
View at: Publisher Site | Google Scholar
K.-B. Min, L. Jing, and O. Stephansson, “Determining the equivalent permeability tensor for fractured rock masses using a stochastic REV approach: Method and application to the field data from Sellafield, UK,” Hydrogeology Journal, vol. 12, no. 5, pp. 497–510, 2004.
View at: Publisher Site | Google Scholar
A. Baghbanan and L. Jing, “Hydraulic properties of fractured rock masses with correlated fracture length and aperture,” International Journal of Rock Mechanics and Mining Sciences, vol. 44, no. 5, pp. 704–719, 2007.
View at: Publisher Site | Google Scholar
K. B. Min, J. Rutqvist, C.-F. Tsang, and L. Jing, “Stress-dependent permeability of fractured rock masses: a numerical study,” International Journal of Rock Mechanics and Mining Sciences, vol. 41, no. 7, pp. 1191–1210, 2004.
View at: Publisher Site | Google Scholar
A. Baghbanan and L. Jing, “Stress effects on permeability in a fractured rock mass with correlated fracture length and aperture,” International Journal of Rock Mechanics and Mining Sciences, vol. 45, no. 8, pp. 1320–1334, 2008.
View at: Publisher Site | Google Scholar
Y.-H. An and Q. Wang, “Analysis of representative element volume size based on 3D fracture network,” Yantu Lixue/Rock and Soil Mechanics, vol. 33, no. 12, pp. 3775–3780, 2012.
View at: Google Scholar
Y.-F. Chen, S.-H. Hu, R. Hu, and C.-B. Zhou, “Estimating hydraulic conductivity of fractured rocks from high-pressure packer tests with an Izbash's law-based empirical model,” Water Resources Research, vol. 51, no. 4, pp. 2096–2118, 2015.
View at: Publisher Site | Google Scholar
Y. Wang, Flow Flow analysis and flow-stress coupled analysis of fissured rock masses (Doctoral dissertation [Doctoral, thesis], Hohai University, Nanjing, China, 1995.
R.-C. Liu, B. Li, Y.-J. Jiang, and L.-Y. Yu, “Effects of equivalent hydraulic aperture and hydraulic gradient on nonlinear seepage properties of rock mass fracture networks,” Yantu Lixue/Rock and Soil Mechanics, vol. 37, no. 11, pp. 3165–3174, 2016.
View at: Publisher Site | Google Scholar
H. Zhou, M. Xiao, and J.-T. Chen, “Study of anchoring mechanism and analysis of anchoring effect of fully grouted rock anchor in large-scale underground caverns,” Rock and Soil Mechanics, vol. 37, no. 5, pp. 1503–1511, 2016.
View at: Publisher Site | Google Scholar
Z. Zhang and J. Nemcik, “Fluid flow regimes and nonlinear flow characteristics in deformable rock fractures,” Journal of Hydrology, vol. 477, pp. 139–151, 2013.
View at: Publisher Site | Google Scholar
Y.-F. Chen, J.-Q. Zhou, S.-H. Hu, R. Hu, and C.-B. Zhou, “Evaluation of Forchheimer equation coefficients for non-Darcy flow in deformable rough-walled fractures,” Journal of Hydrology, vol. 529, pp. 993–1006, 2015.
View at: Publisher Site | Google Scholar
G. I. Barenblatt, I. P. Zheltov, and I. N. Kochina, “Basic concepts in the theory of seepage of homogeneous liquids in fissured rocks [strata],” Journal of Applied Mathematics and Mechanics, vol. 24, no. 5, pp. 1286–1303, 1960.
View at: Publisher Site | Google Scholar
R. W. Zimmerman, G. Chen, T. Hadgu, and G. S. Bodvarsson, “A numerical dual‐porosity model with semianalytical treatment of fracture/matrix flow,” Water Resources Research, vol. 29, no. 7, pp. 2127–2137, 1993.
View at: Publisher Site | Google Scholar
J. E. Warren and P. J. Root, The behavior of naturally fractured reservoirs, 1963, The behavior of naturally fractured reservoirs.
V. Martin, J. Jaffre, and J. E. Roberts, “Modeling fractures and barriers as interfaces for flow in porous media,” SIAM Journal on Scientific Computing, vol. 26, no. 5, pp. 1667–1691, 2005.
View at: Publisher Site | Google Scholar | MathSciNet
F. Ren, G. Ma, Y. Wang, and L. Fan, “Pipe network model for unconfined seepage analysis in fractured rock masses,” International Journal of Rock Mechanics and Mining Sciences, vol. 88, pp. 183–196, 2016.
View at: Publisher Site | Google Scholar
F. Liu, L. Q. Zhao, P. L. Liu, Z. F. Luo, N. Y. Li, and P. S. Wang, An extended finite element model for fluid flow in fractured porous media, Mathematical Problems in Engineering, 2015.
S.-H. Chen, X.-M. Feng, and S. Isam, “Numerical estimation of REV and permeability tensor for fractured rock masses by composite element method,” International Journal for Numerical and Analytical Methods in Geomechanics, vol. 32, no. 12, pp. 1459–1477, 2008.
View at: Publisher Site | Google Scholar
S. H. Chen and X. M. Feng, “Composite element model for rock mass flow flow,” Journal of Hydrodynamics, vol. 18, no. 2, pp. 219–224, 2006.
View at: Google Scholar
V. Tzelepis, K. N. Moutsopoulos, J. N. E. Papaspyros, and V. A. Tsihrintzis, “Experimental investigation of flow behavior in smooth and rough artificial fractures,” Journal of Hydrology, vol. 521, pp. 108–118, 2015.
View at: Publisher Site | Google Scholar
T. W. Schrauf and D. D. Evans, “Laboratory Studies of Gas Flow Through a Single Natural Fracture,” Water Resources Research, vol. 22, no. 7, pp. 1038–1050, 1986.
View at: Publisher Site | Google Scholar
B. Li, R. Liu, and Y. Jiang, “Influences of hydraulic gradient, surface roughness, intersecting angle, and scale effect on nonlinear flow behavior at single fracture intersections,” Journal of Hydrology, vol. 538, pp. 440–453, 2016.
View at: Publisher Site | Google Scholar
R. W. Zimmerman, A. Al-Yaarubi, C. C. Pain, and C. A. Grattoni, “Non-linear regimes of fluid flow in rock fractures,” International Journal of Rock Mechanics and Mining Sciences, vol. 41, no. 1, pp. 1–7, 2004.
View at: Publisher Site | Google Scholar
R. C. Liu, Y. J. Jiang, B. Li, X. S. Wang, B. S. Xu, and L. Y. Yu, “A study of hydraulic properties of rock fracture networks based on rank-one inverse Broyden quasi-Newton method,” Rock and Soil Mechanics, vol. 37, no. 1, pp. 219–228, 2016.
View at: Google Scholar
R. Liu, B. Li, and Y. Jiang, “Critical hydraulic gradient for nonlinear flow through rock fracture networks: The roles of aperture, surface roughness, and number of intersections,” Advances in Water Resources, vol. 88, pp. 53–65, 2016.
View at: Publisher Site | Google Scholar
N. Frih, J. E. Roberts, and A. Saada, “Modeling fractures as interfaces: A model for Forchheimer fractures,” Computational Geosciences, vol. 12, no. 1, pp. 91–104, 2008.
View at: Publisher Site | Google Scholar
P. Knabner and J. E. Roberts, “Mathematical analysis of a discrete fracture model coupling Darcy flow in the matrix with Darcy-Forchheimer flow in the fracture,” ESAIM: Mathematical Modelling and Numerical Analysis, vol. 48, no. 5, pp. 1451–1472, 2014.
View at: Publisher Site | Google Scholar
N. Frih, J. E. Roberts, and A. Saada, “Un modèle Darcy-Forchheimer pour un écoulement dans un milieu poreux fracturé,” Revue Africaine de la Recherche en Informatique et Mathématiques Appliquées, vol. 5, pp. 129–143, 2006.
View at: Google Scholar
Abaqus, V, “6.14 Documentation, Systemes Simulia Corporation,” 2014.
Z. Chen, A. P. Bunger, X. Zhang, and R. G. Jeffrey, “Cohesive zone finite element-based modeling of hydraulic fractures,” Acta Mechanica Solida Sinica, vol. 22, no. 5, pp. 443–452, 2009.
View at: Publisher Site | Google Scholar
B. Carrier and S. Granet, “Numerical modeling of hydraulic fracture problem in permeable medium using cohesive zone model,” Engineering Fracture Mechanics, vol. 79, pp. 312–328, 2012.
View at: Publisher Site | Google Scholar
A. R. Lamb, “Coupled deformation, fluid flow and fracture propagation in porous media,” 2011.
View at: Google Scholar
D. Feng, Distribution of Cracks' Geometry and Saturated Unsatured Intrinsic permeability of Cracks Soil (Doctoral dissertation [Doctoral, thesis], Hohai University, Nanjing, China, 2013.
J. Q. Zhou, S. H. Hu, S. Fang, Y. F. Chen, and C. B. Zhou, “Nonlinear flow behavior at low Reynolds numbers through rough-walled fractures subjected to normal compressive loading,” International Journal of Rock Mechanics and Mining Sciences, vol. 80, pp. 202–218, 2015.
View at: Google Scholar
Y. Wang, F. Qin, Z. Xia, and X. Ni, “Non-darcy flow model and numerical simulation for predicting water inflow in deep tunnel,” Yanshilixue Yu Gongcheng Xuebao/Chinese Journal of Rock Mechanics and Engineering, vol. 31, no. 9, pp. 1862–1868, 2012.
View at: Google Scholar
P. G. Ranjith and W. Darlington, “Nonlinear single-phase flow in real rock joints,” Water Resources Research, vol. 43, no. 9, Article ID W09502, 2007.
View at: Publisher Site | Google Scholar
C. Cherubini, C. I. Giasi, and N. Pastore, “Bench scale laboratory tests to analyze non-linear flow in fractured media,” Hydrology and Earth System Sciences, vol. 16, no. 8, pp. 2511–2522, 2012.
View at: Publisher Site | Google Scholar
M. Javadi, M. Sharifzadeh, K. Shahriar, and Y. Mitani, “Critical Reynolds number for nonlinear flow through rough-walled fractures: the role of shear processes,” Water Resources Research, vol. 50, no. 2, pp. 1789–1804, 2014.
View at: Publisher Site | Google Scholar
J.-Q. Zhou, S.-H. Hu, Y.-F. Chen, M. Wang, and C.-B. Zhou, “The friction factor in the forchheimer equation for rock fractures,” Rock Mechanics and Rock Engineering, vol. 49, no. 8, pp. 3055–3068, 2016.
View at: Publisher Site | Google Scholar
L. Yu, R. Liu, and Y. Jiang, “A review of critical conditions for the onset of nonlinear fluid flow in rock fractures,” Geofluids, vol. 2017, Article ID 2176932, 17 pages, 2017.
View at: Publisher Site | Google Scholar
G. Kosakowski and B. Berkowitz, “Flow pattern variability in natural fracture intersections,” Geophysical Research Letters, vol. 26, no. 12, pp. 1765–1768, 1999.
View at: Publisher Site | Google Scholar
Y. Bachmat, “Basic transport coefficients as aquifer characteristics. In Sym-Symposium on Hydrology of Fractured Rocks,” in In IAHS Sym-Symposium on Hydrology of Fractured Rocks, Dubrovinik, Croatia, 1965.
View at: Google Scholar
C. T. Hsu and P. Cheng, “Thermal dispersion in a porous medium,” International Journal of Heat and Mass Transfer, vol. 33, no. 8, pp. 1587–1597, 1990.
View at: Publisher Site | Google Scholar
X. Wang, F. Thauvin, and K. K. Mohanty, “Non-Darcy flow through anisotropic porous media,” Chemical Engineering Science, vol. 54, no. 12, pp. 1859–1869, 1999.
View at: Publisher Site | Google Scholar
M. G. Sidiropoulou, K. N. Moutsopoulos, and V. A. Tsihrintzis, “Determination of Forchheimer equation coefficients a and b,” Hydrological Processes, vol. 21, no. 4, pp. 534–554, 2007.
View at: Publisher Site | Google Scholar
B. J. Eck, M. E. Barrett, and R. J. Charbeneau, “Forchheimer flow in gently sloping layers: Application to drainage of porous asphalt,” Water Resources Research, vol. 48, no. 1, Article ID W01530, 2012.
View at: Publisher Site | Google Scholar
E. Ghane, N. R. Fausey, and L. C. Brown, “Non-Darcy flow of water through woodchip media,” Journal of Hydrology, vol. 519, pp. 3400–3409, 2014.
View at: Publisher Site | Google Scholar
R. D. Evans and F. Civan, “Characterization of non-Darcy multiphase flow in petroleum bearing formation. Final report,” Tech. Rep., Oklahoma Univ., Norman, OK (United States). School of Petroleum and Geological Engineering, 1994.
View at: Google Scholar
J. Geertsma, “stimating the Coefficient of Inertial Resistance in Fluid Flow Through Porous Media,” SPE Journal, vol. 14, no. 14, pp. 445–450, 1974.
View at: Publisher Site | Google Scholar
I. F. Macdonald, M. S. El-Sayed, K. Mow, and F. A. L. Dullien, “Flow through porous media-the Ergun Eq. revisited,” Engineering Chemistry Fundamentals, vol. 18, no. 3, pp. 199–208, 1979.
View at: Google Scholar
D. C. Li and T. W. Engler, “Literature review on correlations of the non-Darcy coefficient. In Permian Basin Oil and Gas Recovery Conference,” 2001.
View at: Google Scholar
M. Veyskarami, A. H. Hassani, and M. H. Ghazanfari, “Modeling of non-Darcy flow through anisotropic porous media: Role of pore space profiles,” Chemical Engineering Science, vol. 151, pp. 93–104, 2016.
View at: Publisher Site | Google Scholar
P. M. Knupp and J. L. Lage, “Generalization of the Forchheimer-extended Darcy flow model to the tensor permeability case via a variational principle,” Journal of Fluid Mechanics, vol. 299, pp. 97–104, 1995.
View at: Publisher Site | Google Scholar

Copyright

Copyright © 2018 Yuan Wang et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

PDF Download Citation

Download other formats

Order printed copies

Views

1303

Downloads

1330

Citations