Abstract
This paper proposes a numerical model for the fluid flow in fractured porous media with the extended finite element method. The governing equations account for the fluid flow in the porous medium and the discrete natural fractures, as well as the fluid exchange between the fracture and the porous medium surrounding the fracture. The pore fluid pressure is continuous, while its derivatives are discontinuous on both sides of these high conductivity fractures. The pressure field is enriched by the absolute signed distance and appropriate asymptotic functions to capture the discontinuities in derivatives. The most important advantage of this method is that the domain can be partitioned as nonmatching grid without considering the presence of fractures. Arbitrarily multiple, kinking, branching, and intersecting fractures can be treated with the new approach. In particular, for propagating fractures, such as hydraulic fracturing or network volume fracturing in fissured reservoirs, this method can process the complex fluid leak-off behavior without remeshing. Numerical examples are presented to demonstrate the capability of the proposed method in saturated fractured porous media.
1. Introduction
Flow models estimating the flow in fractured porous media mainly include the equivalent continuum model [1], dual continuum model [2], discrete fracture model [3], and discrete fracture network model [4, 5]. These models fall into two categories: continuum models and discrete fracture models. In the equivalent continuum model, a representative elementary volume (REV) is required, the fractures are assumed to distribute regularly, and the equivalent permeability is hard to determine [6], which make it unavailable while several large fractures existed such as hydraulic fractures, complex fractures network, and big faults. Dual continuum model cannot describe the flow in fracture because the fractures are not treated explicitly in it [7, 8]. Discrete fracture network model is more real and gains worldwide concern recently because the fractures are treated explicitly and the effect of one crack on the whole flow could be considered [9]. Yao Jun et al. [10] presented a two-dimensional two-phase finite element flow model based on explicit treating of discrete fracture network. In the model, matrix domain is partitioned using triangle finite element method (FEM) grids while large fracture is partitioned using one-dimensional line elements. The grids are required to match the explicit fractures and nodes are needed on both fracture intersections and crack tips to better solve such static fracture problems. Moreover, tedious remeshing is necessary for propagating fractures [11]. The extended finite element method (XFEM) is an effective approach to avoid remeshing and is now widely used to solve discontinuity problem [12–14]. Khoei et al. [15] simulated the coupled thermohydromechanical (THM) model for impermeable discontinuities in saturated porous medium with XFEM. Mohammadnejad and Khoei [16] presented a coupled hydromechanical (HM) model with combined XFEM and cohesive crack. Lamb [17] and Lamb et al. [18] combined XFEM with the dual-porosity and dual-permeability model to describe the fluid flow, deformation, and fracture propagation in fractured porous medium. By virtue of XFEM and lower-dimensional interface elements, Watanabe et al. [19] researched coupled HM problems in discrete fractures porous media systems. Based on the methodology of XFEM, in this paper, the pressure field is enriched by the absolute signed distance and appropriate asymptotic functions to develop an XFEM model for fluid flow in fractured porous medium. The results demonstrate that the XFEM is an efficient method for simulating fluid flow in fractured porous medium with nonmatching grids, especially when the fissure is propagating, such as hydraulically driven fractures.
2. Governing Equation
The strong form as well as the associated weak form of governing equations for fluid flow inside both matrix domain and fracture domain is demonstrated in the section.
2.1. Strong Form
According to the law of mass conservation, the continuity equation for flow in the matrix domain is written as
Similarly, the continuity equation for flow inside the fracture domain is given bywhere is the Biot constant, is the porosity of matrix rock, is the bulk modulus of solid phase, Pa, is the bulk modulus of pore fluid, Pa, is the pore fluid pressure in matrix, Pa, is time, s, is the permeability tensor of matrix, m2, is the pore fluid viscosity, Pa·s, is the source/sink term on the ground, m3/(m3·s), is the bulk coefficient, is the fluid pressure in fracture, Pa, and is the permeability of fracture, m2, calculated according to the following formula [20]:where is a morphological parameter accounting for the difference between real fracture and ideal parallel fracture, ranging from 1.04 to 1.65. is the mechanical width of fracture, m.
2.2. Weak Form
In order to deduce the weak form of governing equations, a two-dimensional domain bounded by the boundary is considered as shown in Figure 1. Natural fracture with high permeability inside is regarded as a one-dimensional discontinuous line, because the fracture width is much less than the length. and are used to represent two faces of the fracture, respectively.
The initial and boundary conditions are as follows:where is the unit normal vector of the boundary, is fluid density, kg/m3, and is the flow rate on the boundary, kg/(m2·s).
The weak forms of equations for flow in matrix and fracture are derived by weighted residual methodin which is the fracture domain. The integrals over are performed in the local Cartesian coordinate system , whose axes and axes are consonant with tangential direction and normal direction, respectively.
With regard to the flow inside the fracture, both fluid pressure and its interpolating function are assumed to be uniform along fracture width [16], as shown in Figure 2.
Thus, the first integral of (6) can be rewritten as [16]
Because of the supposed uniform distribution of fluid pressure on the axes, the derivative does not vary with ; thus the second integral is expressed as [16]
Substituting the constituents of (6) and rearranging it, the weak form of continuity equation for flow inside fracture is gained. Adding it to the weak form for flow in matrix region, the fluid exchange term can be eliminated as (9) shows. Thus it is not necessary to express fluid exchange explicitly:
3. Pressure Enrichment with XFEM
For permeable natural fracture, the pressure field is continuous at the fracture surface, while the derivative of pressure (i.e., flow rate) is discontinuous. The enrichment of the discontinuity in derivative is constructed as follows [21].
For a node with a bisected support, the absolute signed distance function is applied to realize the enrichment so as to meet the condition of continuous pressure and discontinuous derivative:
For a support which is slit by the discontinuity, the branch function is applied to realize the enrichment:
Its derivate with respect to isin which is the intersection angle between local coordinate system for fracture and global coordinate system.
The first branch function and its derivative for pore pressure are shown in Figure 3. On the surface of fracture (i.e., and ), the interpolation function is continuous while its derivative is discontinuous. The other two interpolation functions share the same characteristics.
(a) The first branch function
(b) Its derivative for pore pressure
For the nodes cut by two intersecting discontinuities, the enrichment is realized through the product of two absolute signed distance functions related to the two discontinuities:
According to the fundamental XFEM, the local enrichment function and extra degrees of freedom are introduced for interpolation of fluid pressure field. Then the pressure after enrichment can be expressed aswhere is the node interpolation function, is the ramp function, is the set of all nodes, is the set of nodes enriched by absolute signed distance function for fracture is the set of nodes enriched by appropriate asymptotic function for fracture tip is the set of nodes enriched by intersecting function for cross point , is the freedom degree of regular nodes (-freedom degree), is the freedom degree of extra nodes in bisected support (-freedom degree), is the freedom degree of extra nodes in slit support (-freedom degree), and is the freedom degree of extra nodes in intersecting support (-freedom degree).
4. Discretization Equation of XFEM
Regarding trial function as the variation of pore pressure interpolation function, the XFEM discretization equation of (9) in space is then gained considering the irrelevance of the values of , , , and :where flow matrix , compressibility matrix , flow rate vector at the equivalent node, the unit tangent vector of local coordinate system, and coordinate transform matrix are, respectively, expressed as follows:
Equation (15) can be rewritten asin which is the sum of matrix of flow in matrix region and matrix of flow in fracture, is the sum of matrix of rock compressibility and matrix of fracture compressibility, is the vector of pore pressure field, and is the vector of flow rate at the equivalent node.
In the time domain, (17) can be discretized by the use of Newmark method as follows [22]:in which parameter satisfies . Only when , the solution is unconditionally stable, when , it is the central difference in the time domain, and when , it is fully implicit.
5. Numerical Examples
5.1. Steady State Flow
For a two-dimensional fractured domain with an independent fracture shown in Figure 4, Lamb [17] compared the calculated results of discrete fracture model (DFM), fracture mapping with mesh free (FM-Mfree), and fracture mapping with finite element method (FM-FEM) and verified the correctness of fracture mapping method. The domain discretization, pressure field, and pressure profile for a vertical section taken from the center of the domain for each method are, respectively, demonstrated in Figures 5, 6, and 7.
(a) DFM
(b) FM-FEM
(c) FM-Mfree
(d) XFEM
(a) DFM
(b) FM-FEM
(c) FM-Mfree
(d) XFEM
It is observed in Figure 5 that the mesh edges need to be set on the fracture and nodes need to be set on the fracture tip for DFM model, and the mesh node needs to be set on the fracture for FM-Mfree model. However, the FM-FEM and XFEM models can simulate the fluid flow in fractured porous medium with nonmatching grid. The permeability of elements with fracture is processed equivalently for FM-FEM model, while the permeability of fracture is calculated only with cubic law or others for XFEM model. The red square represents the enriched node with a bisected support, hollow red circle represents the branch function enriched node with a silt support, and solid red circle represents the improved branch-function-enriched node for XFEM, as shown in Figure 5(d) [23].
As can be seen from Figures 6 and 7, the calculated pressure distribution and pressure profile are basically the same for 4 methods with different cell number and node number. From this, the XFEM approach proposed in this paper is verified.
Because of the introduction of enrichment functions with discontinuous derivate and extra degrees of freedom for node, it is not necessary to consider the location of fracture in the process of domain discretization, which makes it more convenient when treating complex network fractures. In particular, when fracture propagation is concerned such as that during hydraulic fracturing process, there is no need to remesh and transform nodal pressure filed, with obvious merit.
Adding one more fracture in the domain shown in Figure 4, it becomes intersecting fracture with the rest of parameters and boundary conditions unchanged. Lamb [17] computed the pressure distributions in the existence of intersecting fractures by use of DFM and FM-FEM, respectively. Comparing them in Figure 8 with the result using XFEM, it is recognized that the result of XFEM has a good consistency with that of DFM and FM-FEM, which confirms that the fractured-porous-media flow model proposed in this paper can describe not only the flow process in both fracture and matrix but also process independent and intersecting fractures. Thus, it builds a solid foundation for coupled HM simulation and analysis of fracture propagation during volume stimulation under the concept of XFEM.
(a) DFM
(b) FM-FEM
(c) XFEM
5.2. Transient State Flow
Consider a 40 meters long and 50 meters wide horizontal reservoir with several 1.5-millimeters-wide fractures as shown in Figure 9. The permeability and porosity of matrix rock is 100 × 10−15 m2 and 0.18, respectively. The initial pore fluid pressure is 20 MPa. Both the lower boundary and the upper boundary are impermeable. The left boundary has a constant pressure of 25 MPa, while the right boundary has a constant pressure of 20 MPa. The bulk modulus of water phase and solid rock grains is 2.0 GPa and 134.6 GPa, respectively. The viscosity of water phase is 5 mPa·s.
Applying 39 × 49 structured gird, the distributions of pore fluid pressure at different moments ( s, s, s, and s) are demonstrated in Figure 10 after calculating with XFEM. Before pressure spreading to fracture ( s), the isobaric lines appear to be parallel, which is in the stage of linear displacement (Figure 10(a)). During pressure encountering fracture ( s), the isobaric line bents toward tip direction near the fracture tip, and the fracture tip is presented as a point/sink source at the moment (Figures 10(b)~10(d)). After pressure sweeping over the whole fractured region, the pressure drop of whole fractured region is very small, and pressure loss occurs mainly in the matrix (Figures 10(c) and 10(d)).
(a) s
(b) s
(c) s
(d) s
6. Conclusions and Discussions
(1)Both absolute signed distance function and appropriate asymptotic function are continuous but discontinuous in derivatives, which enables them to be enrichment functions during processing pressure field across permeable fractures.(2)Based on the explicit expression of real distribution of fractures, locally enrichment with XFEM was used to develop fluid flow in fractured porous medium. The model is able to elaborately simulate the flow in reservoir with arbitrary independent, kinking, branching, intersecting, and complex network fractures.(3)Domain discretization is independent of fracture location with XFEM flow model. Moreover, remeshing is unnecessary when dealing with dynamic propagating fractures, which enables it to be used for simulating hydraulic fracturing with coupled HM model.
Conflict of Interests
The authors declare that there is no conflict of interests regarding the publication of this paper.
Acknowledgment
The authors are grateful for the research support of the National Natural Science Foundation of China (NSFC, no. 51404207).