A Novel Numerical Model for Fluid Flow in 3D Fractured Porous Media Based on an Equivalent Matrix-Fracture Network

Yao, Chi; He, Chen; Yang, Jianhua; Jiang, Qinghui; Huang, Jinsong; Zhou, Chuangbing

doi:https://doi.org/10.1155/2019/9736729

Geofluids

On this page

Abstract Introduction Conclusion Data Availability Conflicts of Interest Acknowledgments References Copyright Related Articles

Research Article | Open Access

Volume 2019 | Article ID 9736729 | https://doi.org/10.1155/2019/9736729

A Novel Numerical Model for Fluid Flow in 3D Fractured Porous Media Based on an Equivalent Matrix-Fracture Network

Chi Yao,¹Chen He,¹Jianhua Yang,¹Qinghui Jiang,¹Jinsong Huang,¹and Chuangbing Zhou¹

Academic Editor: Baojun Bai

Received01 Jul 2018

Accepted10 Oct 2018

Published03 Jan 2019

Abstract

An original 3D numerical approach for fluid flow in fractured porous media is proposed. The whole research domain is discretized by the Delaunay tetrahedron based on the concept of node saturation. Tetrahedral blocks are impermeable, and fluid only flows through the interconnected interfaces between blocks. Fractures and the porous matrix are replaced by the triangular interface network, which is the so-called equivalent matrix-fracture network (EMFN). In this way, the three-dimensional seepage problem becomes a two-dimensional problem. The finite element method is used to solve the steady-state flow problem. The big finding is that the ratio of the macroconductivity of the whole interface network to the local conductivity of an interface is linearly related to the cubic root of the number of nodes used for mesh generation. A formula is presented to describe this relationship. With this formula, we can make sure that the EMFN produces the same macroscopic hydraulic conductivity as the intact rock. The approach is applied in a series of numerical tests to demonstrate its efficiency. Effects of the hydraulic aperture of fracture and connectivity of the fracture network on the effective hydraulic conductivity of fractured rock masses are systematically investigated.

1. Introduction

The research for seepage flow in the subsurface system is of great significance in many geological and environmental engineering applications, for example, petroleum and gas exploitation, greenhouse gas storage, and nuclear waste treatment. Most rocks in nature are cut by fracture or fracture network, and the fractures have significant influence on the deformation and seepage characteristics of rocks [1–4]. It is well known that fluid tends to flow through the easiest path to the downstream. For fractured rock, the conductivity of fractures is often several orders of magnitude higher than that of the porous matrix, which means that fractures and fracture networks can strongly influence or even dominate the conductivity and flow paths of rock masses. Moreover, fluid exchange between the porous matrix and fractures often results in complex hydraulic behaviors, such as seepage flow in aquifers and reservoirs. It is generally hard to solve seepage problems in such complicated materials with analytical methods. On the other hand, many numerical simulation methods, such as finite element method [5, 6], finite volume method [7, 8], boundary element method [9], and numerical manifold method [10, 11], provide efficient tools for the analysis of fluid flow in fractured porous media and can help us better predict the hydraulic properties of underground engineering.

There are three major kinds of model for fluid flow simulation in fractured porous media which based on equivalent continuum/porous media (EPM), discrete fracture network (DFN), or fractured porous media (FPM). The EPM-based model assumes that fluid flow obeys Darcy’s law, in which fractures and matrix are approximated as a continuum and the classical continuous seepage theory is used for analysis [12–14]. This kind of model provides a phenomenological description of flow in fractured porous media in large dimension. However, it is limited by whether the REV exists and cannot reflect exactly the local conductivity anisotropy induced by fractures and the preferential flow effect within fracture networks. DFN-based models assume that rock blocks are impermeable and flow movement only takes place in fracture networks [15–17]. The geometry and topology of fracture networks are generally described explicitly by this kind of approach. DFN-based models are applicable only when the fracture network is well-connected, and the conductivity of the rock matrix is negligible. Methods based on the FPM concept consider that fluid flows simultaneously within fractures and the rock matrix [7, 11, 18–21]. This kind of approach inherits the advantages of EPM and DFN. It not only describes the geometry of the fracture network but also takes into consideration the contribution from the matrix and the fluid exchange between matrix and fractures. However, this method has difficulties in numerical mesh generation due to geometry complexity of fracture networks, and the computational cost is generally very high. The three kinds of model mentioned previously have their advantages in fluid seepage simulation, respectively, among which the FPM-based models are usually more comprehensive, more consistent with physical reality, and thus could provide more details.

To overcome the drawbacks of three previously mentioned approaches, Yao et al. [22] proposed a novel equivalent matrix-fracture network (EMFN) approach, in which flows in the rock matrix and fractures are represented by the flow in an equivalent fracture network. However, Yao et al. [22] only considered 2D problems where the Voronoi diagram was used for discrete 2D domains and the one-dimensional equivalent fracture network was used to study the conductivity evolution of fractured rock masses. It is well known that the 2D fracture network tends to underestimate the conductivity compared with those based on 3D description since the connectivity of the 2D model is inferior to the 3D model [5, 23, 24]. As a result, 3D models are necessary for the prediction of hydraulic properties of fractured porous media.

In this study, the EMFN approach is further developed to model 3D fluid flows in fractured porous media. The porous matrix and fracture network are replaced by an equivalent matrix-fracture network (EMFN). Delaunay tetrahedrons are used to construct the 3D fracture network. It is shown that Delaunay tetrahedrons are efficient and able to handle irregular boundary geometry and complex fracture distributions. The proposed approach is applied to study the conductivity of complex fractured rock masses. Results show that there is a unique relationship between the macroscopic conductivity of the whole interface network and the microscopic conductivity of interfaces, the ratio of which is linearly related to the cubic root of the number of nodes used for mesh generation. With this relationship, the so-called equivalent matrix-fracture network (EMFN) can produce the same conductivity of intact rock, which lays the foundation for the proposed model. The effectiveness of this approach is then demonstrated by several examples.

2. The Equivalent Matrix-Fracture Network Model

2.1. The Numerical Model

In this approach, rock masses and fractures are discretized into a triangular interface network in the 3D domain. Fluid flow occurs on the network and follows Darcy’s law. Firstly, a 2D model is used to illustrate the concept of the equivalent matrix-fracture network model. A 2D schematic of one streamline (solid red line) in the model is depicted in Figure 1(a), where the triangular interfaces are reduced to line segments. A 3D conceptual diagram of the interface network is also shown in Figure 1(b). The triangular interfaces between tetrahedrons are the paths through which fluid flows. Rock matrix and fractures are represented by interface networks in the same mesh. The fracture thickness is not considered. The only difference between the matrix interface and the fracture interface is that they are assigned different transmissivities while building the equivalent equation. In this way, the pore structure of the porous media and the fracture network is replaced by an equivalent interface network, which is called the 3D equivalent matrix-fracture network (EMFN). The core concept of the EMFN model is to use a sufficiently complex and chaotic discrete and permeable interface network to fill a three-dimensional volume, making the effective permeability of the interface network equivalent to the continuous media. Then, the fractures and matrix are simultaneously represented as lower-dimensional domains (2D surfaces) in the volume, and the Galerkin-based finite element method is used to solve the flow problem.

(a)

(b)

Assuming that the transmissivity of an individual interface is , the continuity equation of the triangular interface under steady-state flow conditions in the local coordinate (Figure 2(a)) is expressed as where h is the hydraulic head, 1 m of which is equivalent to 9.81 kPa. Suppose an intersection line is the common edge of two triangular fracture surfaces (Figure 2(b)). Neglecting the storage capacity of the porous media, for any node on this line, the equation of mass conservation is where is the number of interfaces passing through this node, and is the flow rate from this node to the interface ,

(a)

(b)

Illustrated in Figure 2(b), is the direction vector that is perpendicular to the intersection line and parallel to the plane . The flow problem in the whole domain also needs to satisfy the hydraulic head boundary conditions and the flux boundary conditions.

The finite element method is used to solve the steady-state flow problem in the EMFN model. The hydraulic head of any node on a triangular element can be represented by the hydraulic head at three vertices (named ) surrounding this element: where , , and are the shape functions. After the finite element approximation, the total equivalence equation of the steady-state flow is built up: where where , , and are local coordinates and and are, respectively, the unknown hydraulic head vector and the vector calculated from the known boundary condition. The solution of the whole seepage flow field can be obtained by solving (5).

2.2. Three-Dimensional Mesh Generator for Complex Fracture Network Based on Delaunay Tessellation

In the meshing process of the research domain with irregular boundary geometry and incorporating complex fracture distribution, auto-unstructured mesh generation techniques, such as Delaunay triangulation (in 2D) and tetrahedralization (in 3D) [7, 25], are relatively appropriate and efficient. In this paper, the 3D research domain is discretized by Delaunay tetrahedrons (Figure 1(b)). Delaunay tetrahedralization is a fast and relatively high-quality procedure for automatic spatial subdivision. Once the spatial node set and the geometric boundary condition have been given, the unique tetrahedral mesh can be generated. The node set exists as vertexes of tetrahedrons in the mesh.

A random node saturation procedure is introduced in the node set generation illustrated in Figure 3. A parameter is used to define the minimum distance between any two nodes in the mesh, then the node keeps being put one by one randomly into the domain until no other node meets the condition. Then, the open-source Delaunay tetrahedralization code, TetGen [26], is applied to discretize the domain based on the node set. The fracture geometry can be described by a special Delaunay constrained facet and edge generation strategy. The process of the mesh generation program is as follows: Step 1:Define a research domain (in this paper, the domain is always a cube with side length as shown in Figure 1(b)), and then nodes with a random 3D coordinate were put into this domain one by one. Some of the nodes are put in firstly because they are used to form the fracture planes, intersection lines, and domain boundaries; the others are used to form the rock matrix and filled the rest of the domain.Step 2:Generate a strict Delaunay tetrahedral mesh using the existing node set.Step 3:Check the integrity of the fracture planes. If it is found that the area of some planes is missing, call the local node infill algorithm (described later in this section) and jump back to the Step 2.Step 4:Output the mesh data.

It needs to be pointed out that, in Step 1, almost all node insertion processes follow the node saturation principle. However, due to the presence of intersecting cracks with a very small dihedral angle, the nodes near their intersection line may break this principle and have a smaller . Fortunately, the number of abnormal nodes is far smaller than the total number of nodes when the fracture network is not quite dense, which means that the effect of the abnormal nodes on the calculation of rock conductivity is negligibly small. This problem needs to be avoided through a graded-mesh-based approach in the future.

Here, we introduce the local node infill algorithm in Step 3. For a node set (including the fracture geometry), the mesh generated by Delaunay’s method is unique (the situation of 5 or more nodes on one spherical surface can be avoided by disturbing the coordinates of nodes). Unfortunately, some of the tetrahedrons adopted may pass through these fracture planes, resulting in the loss of the area of some fracture planes. A schematic fracture plane is illustrated in Figure 4(a), which lost a triangular area. The reason is shown in Figure 4(b): three tetrahedrons (pqab, pqbc, and pqca) passed through the plane, causing the plane to lose a triangular area abc, where pq is the common edge of three tetrahedrons.

(a)

(b)

(c)

In rigid Delaunay tetrahedralization, the circumscribed sphere of any tetrahedron includes no other nodes besides four vertexes of the tetrahedron. So, one can appropriately put some nodes into the domain appropriately to break these undesirable tetrahedrons and re-form new tetrahedrons that do not pass through fracture planes. In fact, these nodes are the intersection nodes of fracture planes with the common edge of abnormal tetrahedrons, for example, the intersection node of pq with the plane abc in Figure 4(c). The final mesh is illustrated in Figure 5, where the green interface represents the fracture planes and the blue interface represents the equivalent matrix (for clarity, some interfaces close to the boundary are neglected).

2.3. Relationship between Mesh Parameter and Macroscopic Conductivity of EDFN Model

By adopting the random node saturation method, the density of the node set is determined through . In fact, there is a linear relationship between and the total node number of the cubic domain with side length :

Formula (7) is derived from 10 groups of numerical tests, in which , and each group contains 8 specimens. Numerical results are shown in Figure 6; the solid line is the result of a linear fit. Moreover, the coefficient of variation of of one group is, respectively, 1.10%, 0.12%, and 0.01%, when equals 10, 50, and 200.

In the EMFN model, the pore system of the rock matrix is replaced by an interconnected equivalent matrix-fracture network, which is constructed by an assembly of triangular interfaces. The hydraulic conductivity of each equivalent matrix-fracture is used to ensure this system is hydraulically equivalent to the rock matrix. In this paper, all are assumed to be identical. To verify the effectiveness of this method and to study the mesh dependency of the EMFN model at the same time, a group of numerical experiments is designed. There are 6 groups (each group contains 8 specimens) of meshes without any fracture, and the only difference between groups is . The research domain is cubic and . A pair of opposite boundaries are hydraulic head boundaries, and the others are impervious boundaries. According to Darcy’s law,

Here, . Defining dimensionless quantity to evaluate the transmission capacity of the research domain relating to , where the subscript represents the direction (-, -, or -direction) in the global coordinate system. The results are listed in Table 1, where is the average number of nodes, is the average value of one group in the direction, and is the coefficient of variation of , , and in a group. One can see that the Delaunay mesh generated from a randomly and uniformly distributed set of nodes with identical on each equivalent fracture has almost the same in three dimensions, which means the equivalent permeability of the matrix-fracture network generated by the process in Section 2.2 is almost isotropic. Further, the variation coefficient of conductivity becomes smaller with the decrease in .

From Figure 7, one can see that and have a linear relationship:

Substituting (10) into (9) results in

According to Darcy’s theorem, macroscopic flow rate and macroscopic conductivity coefficient have the following relationship:

It can be derived from (11) and (12) that

When using the mesh density parameter as an argument, (13) can be rewritten as

At present, the relationship between the macroscopic conductivity coefficient and the microscopic parameter of the EMFN model has been established.

2.4. A Modified Method for Complex Fracture Network

The function describes the relationship between macroscopic hydraulic parameters and the microscopic parameters in EMFN. It is obtained from data fitting for the mesh that does not include any fracture, which means that the nodes used to generate the Delaunay-based mesh satisfy the “random node saturation procedure” in any part of the research domain. In other words, the density of chaos-oriented interfaces is homogeneous, and the influence of a slight disturbance in node density is negligible. However, when it comes to complex fracture networks (such as those shown in Figure 8, where the green disks are fractures and the blue cube is the research domain), the effect of additional refinement nodes cannot be ignored because they are too many. Where the refinement nodes are added, the density of the interface will be higher than that of other areas where they do not contain refinement nodes. Refinement nodes exist only in the fracture plane. Under this feature, according to the formation rule of Delaunay’s tetrahedron and the EMFN model, the equivalent permeability of the matrix will be affected only within the distance of one tetrahedron around the refinement node, and the transmissivity within the fracture plane will not be affected. The denser the interface network is, the higher the matrix network equivalent permeability is.

(a)

(b)

(c)

In light of the problem, a practical strategy to reduce the influence of refinement nodes on the equivalent permeability of the matrix network is developed. We mark all the refinement nodes added by Step 3 introduced in Section 2.2 and then find out in all the matrix interfaces that at least one of their vertex nodes is the refinement node. Finally, the conductivity of these interfaces is reduced because they are denser than we expected. Through a series of tests, a pragmatic reduction coefficient, 0.55, was determined. That is, if is the transmissivity of matrix interfaces used in (14), then the transmissivity of these dense matrix interfaces is 0.45 . We cannot prove mathematically that the reduction of transmissivity is correct, but it works.

A series of experiments were carried out to verify the effectiveness of the reduction strategy. We designed a research domain that , and in which the centers of fractures are randomly distributed, and their orientation is random and isotropic with periodic boundary conditions (i.e., the example in Figure 8). The radius of fractures is identical and equal to 0.25 L. Under these conditions, samples with different fracture densities were tested. Fracture densities are assessed using the number of fractures inside the research domain. In this test, all the transmissivities of interfaces are identical and equal to . Thus, although fractures are indeed described in the mesh, the fracture interfaces are no different from the matrix interfaces. So the influence of nonhomogeneous interface density and the effect of transmissivity reduction strategy can be easily known through steady-state flow simulation. We define the hydraulic conductivity of the matrix calculated by (14) as and the hydraulic conductivity of fractured media obtained through EMFN as . As shown in Figure 9, the black squares are under the reduction strategy, which does not change significantly with the increase in fracture number (more specifically, the rise of refinement node). However, the results (red spots) obtained from the samples that not have a reduction procedure are obviously larger. The percentages of refinement nodes in total nodes are also noted in Figure 9.

2.5. Verification of the EMFN Model

In the EMFN model, either fracture interfaces or the equivalent matrix interfaces are all treated in the same way of the equivalent-fracture network, in which the “fracture interfaces” and the “matrix interfaces” are given independently; for fracture it is , and for matrix it is . Unfortunately, for flow problems in the discrete fracture porous media, it is hard to find an analytic solution of conductivity even for very simple conditions [27]. However, we can use the superposition principle in conductivity estimation to estimate the conductivity in some well-designed conditions. Two examples are given in this section. The first one is to test the conductivity of the rock a single rectangular fracture passing through the domain; the other one is to examine the conductivity tensor of the rock which has a disk fracture embedded with varying orientations in the centroid of the domain.

When a rectangular fracture is parallel to the hydraulic gradient (sketching is depicted in Figure 10), its contribution to flow rate in the direction of is easy to get by where is the flow rate through fracture and is the trace length of the fracture. In this example, , , , and . Substituting (15) into (9) results in the -direction as . The transmission capacity of the matrix can be calculated by (9), (12), and (14) as . The result of simulation using the EMFN model is .

(a)

(b)

The second example is illustrated in Figure 11(a): a disk fracture (30 m in diameter) is embedded in the cubic domain, , , and . The boundary condition is depicted in Figure 11(c), and a pair of boundaries perpendicular to the -axis are given, respectively, constant hydraulic head and ; others are corresponding linear variable hydraulic head boundaries. The center of the disk fracture coincides with the centroid of the cubic domain; the normal vector of the fracture is always perpendicular to the -axis. is defined as the angle formed by the -axis and . The conductivity tensors when equal to 0°, 15°, 30°, 45°, 60°, 75°, and 90° are calculated.

(a)

(b)

(c)

The conductivity tensor is the function of and can be expressed as

When , where are, respectively, the coefficient of conductivity in the directions. Components in (16) can be written as

Substituting adopted from the numerical simulation when (here , ) into (16), (17), and (18), the relationships between and components of are received. is defined as the coefficient of conductivity simulated in the same conditions but without any fracture (in this case, and equal to , which means the fracture has no improvement on the conductivity in the direction parallel to ). The comparison between the numerical result using the EMFN model and the analytical result calculated by (16), (17), and (18) is plotted in Figure 12. One can see that the numerical results obtained are very close to the theoretical prediction values, which proves the validity of the EMFN method proposed in this paper.

3. Effects of Fracture Geometry on the Effective Conductivity

Different from general models based on the DFN model, the EMFN model is very suitable for studying the effects of the sparse fracture network on the effective conductivity of rock masses in which the contribution of the matrix to the macroscopic conductivity cannot be neglected. In this section, the effect of the hydraulic aperture, the relative position between fractures, and the connectivity of the fracture network on the conductivity of rock mass are studied by the EMFN model.

3.1. Hydraulic Aperture of Fracture

In this example, a single disk fracture is embedded in the cubic domain (). The fracture is parallel to the xoy-plane and has the same centroid as the cube (illustrated in Figure 13). The diameter of the disk is set as, respectively, 40 m, 50 m, and 60 m. Mesh . is defined as the coefficient of conductivity in the direction, and is the same coefficient but in the condition that no fracture was embedded. In this paper, the fracture is simplified to the disk plane. According to the cubic law, the aperture of fracture is proportional to the transmissivity in the plane: where and are, respectively, the coefficient of kinematic viscosity of water and the gravity acceleration. A series of simulations are conducted by varying the fracture transmissivity and the diameter of fracture. Results are depicted in Figure 14, where the ordinate represents the relative hydraulic conductivity and the abscissa is the ratio of fracture transmissivity and the transmissivity of the equivalent matrix-fracture network . It can be observed that hydraulic conductivity increases until grows to 10⁵, after which effective conductivity in the direction will remain constant and the increase in has no more impact on .

3.2. Interaction between Two Fractures

From Section 2.5, we know that a single fracture has almost no contribution to the macroscopic conductivity of rock in the same direction of its normal vector. However, the interaction of multiple fractures may have a complex impact on seepage behaviors of rock. Illustrated in Figure 15, two disk fractures (, ) are embedded in the research domain. For the mesh, and . Fracture 1 is put at the center of the cubic domain, and its normal vector is parallel to the -axis; for fracture 2, the normal vector is parallel to the -axis, and the center of the fracture moves along the centerline of the cube in the -direction. Thus, we can use the value of the -coordinate () of the centroid node of the fracture to describe its spatial position. We propose to test the coefficient of conductivity in three directions (parallel to three axes), which are . The change in relative conductivity in three directions is shown in Figure 16. One can observe that as fracture 2 moves, trends of , , and are completely different.

From the test in Section 3.1, we obtained that when there is only one single fracture (, , and ). In this section, also tends to 1.16 in the direction when fracture 2 is moving away from fracture 1 and, notably, when fracture 2 is at the center of the cubic domain. Figure 17 illustrates the contours of the hydraulic head at the intersecting surface when the direction of is parallel to the -axis, which corresponds to the black squares in Figure 16. Obviously, the hydraulic head field in the cubic domain is symmetrical about this section. One can observe that the hydraulic head field is significantly influenced by the appearance of fracture 2, which makes the difference of . In fact, it can be seen from the distribution of the hydraulic head field on the fractures that although the individual fracture 2 does not contribute to the conductivity of rock mass in the direction, it helps fracture 1 to obtain more flow recharge, making much larger than the one when only fracture 1 is embedded. In particular, the hydraulic head field in Figures 17(a) and 17(d) are almost the same. In Figure 17(d), there is almost no hydraulic head change in fracture 2, so it barely helps the fluid flow in fracture 1, which explains the same value of in Figure 16 at and . Moreover, in Figure 17, we can see that the change of the hydraulic head in fractures is much slower than in the rock matrix; in other words, it consumes much less energy of fluid flowing through the fracture than the rock matrix in the same distance, which is reasonable.

(a)

(b)

(c)

(d)

3.3. Rock Mass with Multiple Fracture Clusters

Illustrated in Figure 18 is a cubic domain with 3 fracture clusters embedded, and the distribution parameters are listed in Table 2. Following Baecher et al.’s model [28], fractures are supposed to be disk-shaped. The orientation of the normal vector of fractures obeys the Fisher distribution in each cluster, where is the dispersion parameter about the mean direction of fracture normal vectors [29]. The centroid position and size of fractures follow, respectively, uniform distribution and negative exponential distribution. Here, , , and ; the hydraulic head upstream and downstream boundaries are, respectively, 100 m and 0 m.

(a)

(b)

(c)

(d)

To show the hydraulic head field clearly, the equivalent matrix-fractures are hidden. Figure 18(a) is the hydraulic head in fractures. Obviously, the fracture network connects the upstream and downstream boundaries. Figure 18(c) is the hydraulic head contour in the rock matrix; the black thick solid line is the intersecting line formed by the contour plane and macroscopic fractures. To compare with the previous example, the geometric structure of the fracture network in Figures 18(b) and 18(d) is identical to Figures 18(a) and 18(c), but the fractures intersecting with plane are all eliminated, which means that the fracture network does not penetrate the domain in the -direction. One can see that the connectivity of the fracture network significantly influences the hydraulic head field in fractured porous media. When fractures cut between each other and form a good percolation pathway, the distribution of the hydraulic gradient is relatively large along the path of fractures (the left side in Figure 18). However, for the nonpercolating situation (see Figure 18(b), the fracture did not cross the research domain in the -direction); the change of the hydraulic head occurs mainly in the matrix (see Figure 18(d)). The relative conductivity under different is depicted in Figure 19 (the data of dark solid lines is taken from the example shown in Figure 18, and other lighter ones are obtained from ten parallel tests under the same fracture distributions). It can be seen from Figure 19 that the relative conductivity of fractured rock mass increases with the increase in . The difference is when ; still increases and has a linear relationship with in the percolating one, but becomes almost a constant in the nonpercolating one. Corresponding to the situation of single fracture in Section 3.1, the fracture network that does not penetrate the entire porous media which can enhance the conductivity of the whole domain, but the conductivity of the rock matrix limits its upper limit.

4. Conclusion and Perspective

In this study, an original method for simulating fluid flow in three-dimensional fractured porous media is proposed. The porous matrix and fractures are hydraulically replaced by an equivalent matrix-fracture network (EMFN). The finite element method is used to solve steady-state seepage flow in discretized triangular mesh. This method can naturally describe geometrical properties of each fracture and the contribution on equivalent conductivity of the fracture network and porous matrix.

In the proposed method, two major obstacles have been overcome: one is the 3D mesh generation techniques based on the constrained Delaunay algorithm that can describe the geometry and topology of the fracture network and the rock matrix surrounding it; the other one is the relationship between macroscopic hydraulic parameters and the microscopic parameters in EMFN which ensures that the equivalent conductivity of fractured porous media is correctly simulated. The validity of the original method is verified through two simple examples from the perspective of equivalent conductivity and conductivity tensor. Moreover, the code is subsequently applied to study the effects of a single fracture aperture and the relative position of two fractures on the equivalent conductivity of fractured porous media. It is found that isolated fracture and fracture clusters surrounded by the porous media can also significantly influence the hydraulic head field distribution in its vicinity and increase the equivalent conductivity as those percolating ones. This effect depends on the properties of fractures individually and their relative position, which means the effective conductivity of fractured porous media cannot be obtained by simply accumulating the contributions from the porous matrix and each fracture. Finally, the code is used to simulate the steady-state flow in rock with 3 fracture clusters, and one controlled test is constructed. The results show that for fractured porous media with similar geometrical statistical parameters, the contribution of the matrix to equivalent conductivity becomes smaller as the conductivity of fractures increases. But for those where the fracture network does not penetrate the entire media, the equivalent conductivity will eventually reach an upper limit and will not increase as the conductivity of the fracture network increases.

Up to now, the original method has been successfully used in steady-state flow simulation in fractured porous media. However, there are still several areas that need to be improved in the future. The density of the mesh is required to be uniformly distributed in space in the current method, which is computationally demanding when the fracture network is sparse and could cause the higher evaluation of matrix conductivity when the fracture network is too dense. Therefore, it is important to develop a high-quality variable density mesh generation technology. Moreover, the method can be further developed for hydraulic-mechanical coupling problems based on 3D DFN representation (as is done in the previous work for the 2D problem [22]), and for the stability analyses for the fractured rock [30–32].

Data Availability

The data used to support the findings of this study are included within the article.

Conflicts of Interest

The authors declare that there is no conflict of interest regarding the publication of this paper.

Acknowledgments

This work was supported by the National Natural Science Foundation of China (CN) (grant numbers U1765207, 41762020, and 51509104).

References

N. Chen, J. Kemeny, Q. Jiang, and Z. Pan, “Automatic extraction of blocks from 3D point clouds of fractured rock,” Computers & Geosciences, vol. 109, pp. 149–161, 2017.
View at: Publisher Site | Google Scholar
Q. Jiang and C. Zhou, “A rigorous solution for the stability of polyhedral rock blocks,” Computers and Geotechnics, vol. 90, pp. 190–201, 2017.
View at: Publisher Site | Google Scholar
Q. Jiang, Y. Sun, B. Yi, T. Li, and F. Xiong, “Inverse analysis for geomaterial parameter identification using Pareto multiobjective optimization,” International Journal for Numerical and Analytical Methods in Geomechanics, vol. 42, no. 14, pp. 1698–1718, 2018.
View at: Publisher Site | Google Scholar
Q. Jiang and C. Zhou, “A rigorous method for three-dimensional asymmetrical slope stability analysis,” Canadian Geotechnical Journal, vol. 55, no. 4, pp. 495–513, 2018.
View at: Publisher Site | Google Scholar
N. Huang, Y. Jiang, B. Li, and R. Liu, “A numerical method for simulating fluid flow through 3-D fracture networks,” Journal of Natural Gas Science and Engineering, vol. 33, pp. 1271–1281, 2016.
View at: Publisher Site | Google Scholar
J. Noorishad and M. Mehran, “An upstream finite element method for solution of transient transport equation in fractured porous media,” Water Resources Research, vol. 18, no. 3, pp. 588–596, 1982.
View at: Publisher Site | Google Scholar
I. I. Bogdanov, V. V. Mourzenko, J. F. Thovert, and P. M. Adler, “Effective permeability of fractured porous media in steady state flow,” Water Resources Research, vol. 39, no. 1, 2003.
View at: Publisher Site | Google Scholar
Y. Caillabet, P. Fabrie, P. Landereau, B. Noetinger, and M. Quintard, “Implementation of a finite-volume method for the determination of effective parameters in fissured porous media,” Numerical Methods for Partial Differential Equations, vol. 16, no. 2, pp. 237–263, 2000.
View at: Publisher Site | Google Scholar
A. Teimoori, Z. Chen, S. S. Rahman, and T. Tran, “Effective permeability calculation using boundary element method in naturally fractured reservoirs,” Petroleum Science and Technology, vol. 23, no. 5-6, pp. 693–709, 2005.
View at: Publisher Site | Google Scholar
M. Hu, Y. Wang, and J. Rutqvist, “Fully coupled hydro-mechanical numerical manifold modeling of porous rock with dominant fractures,” Acta Geotechnica, vol. 12, no. 2, pp. 231–252, 2017.
View at: Publisher Site | Google Scholar
Q. H. Zhang, S. Z. Lin, Z. Q. Xie, and H. D. Su, “Fractured porous medium flow analysis using numerical manifold method with independent covers,” Journal of Hydrology, vol. 542, pp. 790–808, 2016.
View at: Publisher Site | Google Scholar
A. M. Hitchmough, M. S. Riley, A. W. Herbert, and J. H. Tellam, “Estimating the hydraulic properties of the fracture network in a sandstone aquifer,” Journal of Contaminant Hydrology, vol. 93, no. 1–4, pp. 38–57, 2007.
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,” Water Resources Research, vol. 18, no. 3, pp. 645–658, 1982.
View at: Publisher Site | Google Scholar
J. F. Pankow, R. L. Johnson, J. P. Hewetson, and J. A. Cherry, “An evaluation of contaminant migration patterns at two waste disposal sites on fractured porous media in terms of the equivalent porous medium (EPM) model,” Journal of Contaminant Hydrology, vol. 1, no. 1-2, pp. 65–76, 1986.
View at: Publisher Site | Google Scholar
N. Koudina, R. Gonzalez Garcia, J. F. Thovert, and P. M. Adler, “Permeability of three-dimensional fracture networks,” Physical Review E, vol. 57, no. 4, pp. 4466–4479, 1998.
View at: Publisher Site | Google Scholar
C. Yao, Q. H. Jiang, W. Wei, and C. B. Zhou, “The variational inequality formulation for unconfined seepage through three-dimensional dense fracture networks,” Science China Technological Sciences, vol. 56, no. 5, pp. 1241–1247, 2013.
View at: Publisher Site | Google Scholar
Q. H. Zhang and J. M. Yin, “Solution of two key issues in arbitrary three-dimensional discrete fracture network flow models,” Journal of Hydrology, vol. 514, pp. 281–296, 2014.
View at: Publisher Site | Google Scholar
M. Belayneh, S. Geiger, and S. K. Matthäi, “Numerical simulation of water injection into layered fractured carbonate reservoir analogs,” AAPG Bulletin, vol. 90, no. 10, pp. 1473–1493, 2006.
View at: Publisher Site | Google Scholar
I. I. Bogdanov, V. V. Mourzenko, J. F. Thovert, and P. M. Adler, “Effective permeability of fractured porous media with power-law distribution of fracture sizes,” Physical Review E, vol. 76, no. 3, article 036309, 2007.
View at: Publisher Site | Google Scholar
V. V. Mourzenko, I. I. Bogdanov, J. F. Thovert, and P. M. Adler, “Three-dimensional numerical simulation of single-phase transient compressible flows and well-tests in fractured formations,” Mathematics and Computers in Simulation, vol. 81, no. 10, pp. 2270–2281, 2011.
View at: Publisher Site | Google Scholar
V. Reichenberger, H. Jakobs, P. Bastian, and R. Helmig, “A mixed-dimensional finite volume method for two-phase flow in fractured porous media,” Advances in Water Resources, vol. 29, no. 7, pp. 1020–1036, 2006.
View at: Publisher Site | Google Scholar
C. Yao, Q. H. Jiang, and J. F. Shao, “A numerical analysis of permeability evolution in rocks with multiple fractures,” Transport in Porous Media, vol. 108, no. 2, pp. 289–311, 2015.
View at: Publisher Site | Google Scholar
N. Huang, Y. Jiang, R. Liu, and B. Li, “Estimation of permeability of 3-D discrete fracture networks: an alternative possibility based on trace map analysis,” Engineering Geology, vol. 226, pp. 12–19, 2017.
View at: Publisher Site | Google Scholar
P. S. Lang, A. Paluszny, and R. W. Zimmerman, “Permeability tensor of three-dimensional fractured porous rock and a comparison to trace map predictions,” Journal of Geophysical Research: Solid Earth, vol. 119, no. 8, pp. 6288–6307, 2014.
View at: Publisher Site | Google Scholar
T. Taniguchi and E. Fillion, “Numerical experiments for 3-dimensional flow analysis in a fractured rock with porous matrix,” Advances in Water Resources, vol. 19, no. 2, pp. 97–107, 1996.
View at: Publisher Site | Google Scholar
H. Si, “TetGen, a Delaunay-based quality tetrahedral mesh generator,” ACM Transactions on Mathematical Software, vol. 41, no. 2, pp. 1–36, 2015.
View at: Publisher Site | Google Scholar
M. F. Lough, S. H. Lee, and J. Kamath, “An efficient boundary integral formulation for flow through fractured porous media,” Journal of Computational Physics, vol. 143, no. 2, pp. 462–483, 1998.
View at: Publisher Site | Google Scholar
G. B. Baecher, N. A. Lanney, and H. H. Einstein, “Statistical description of rock properties and sampling,” in The 18th U.S. Symposium on Rock Mechanics (USRMS), Golden, CO, USA, June 1977.
View at: Google Scholar
M. Khamforoush, K. Shams, J. F. Thovert, and P. M. Adler, “Permeability and percolation of anisotropic three-dimensional fracture networks,” Physical Review E, vol. 77, no. 5, article 056307, 2008.
View at: Publisher Site | Google Scholar
B. Wang, C. Yao, J. Yang, and S. Jiang, “Numerical simulation of macro-meso mechanical behaviours of sandstone containing a single open fissure under uniaxial compression,” European Journal of Environmental and Civil Engineering, vol. 22, Supplement 1, pp. s99–s113, 2018.
View at: Publisher Site | Google Scholar
J. H. Yang, C. Yao, Q. H. Jiang, W. B. Lu, and S. H. Jiang, “2D numerical analysis of rock damage induced by dynamic in-situ stress redistribution and blast loading in underground blasting excavation,” Tunnelling and Underground Space Technology, vol. 70, pp. 221–232, 2017.
View at: Publisher Site | Google Scholar
C. Yao, J. F. Shao, Q. H. Jiang, and C. B. Zhou, “Numerical study of excavation induced fractures using an extended rigid block spring method,” Computers and Geotechnics, vol. 85, pp. 368–383, 2017.
View at: Publisher Site | Google Scholar

Copyright

Copyright © 2019 Chi Yao 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

1865

Downloads

1291

Citations