Abstract
Reservoir simulation is critical to the design of reservoir development plan and has been extensively used. However, it is challengeable for to simulate the production process for fractured reservoirs, because of the fracture geometry and the fracture deformation. Specifically, three problems need to be solved. First, there is a lack of mathematical models that can predict the fracture deformation with acceptable precision. Second, the fracture deformation is stress-dependent; therefore, a geocoupled equation should be used to quantify the stress change, but the solution is extremely expensive. Third, the fracture geometries pose great challenges to traditional gridding techniques. This paper proposes a new geocoupling simulation method that is capable of modeling the complex fracture geometry as well as the fracture mechanical behavior. The geomechanical effects and the reservoir production performance are modeled through an implicit geocoupled model, which is developed based on the poro-mechanics theory. The fixed stress strategy is used to solve the geocoupled equations. Moreover, a comprehensive fracture modeling method is proposed, in terms of the fracture deformation model and the fracture gridding technique to model the fracture effects. Ultimately, this method is used to analyze two field-scale cases. The results demonstrate that this method exhibits good practicability and has practical significance for fractured reservoir development.
1. Introduction
In recent years, more and more attention has been attracted to fractured reservoirs. Actually, most unconventional reservoirs can be classified as fractured reservoirs. In fractured reservoir development, fracture is the major flow channel. Understanding the fracture effects on production is important, but is also extremely difficult. The reasons are:
1.1. Fracture Geometry
Fracture geometry is critical to reservoir simulation but is difficult to be integrated into the simulation algorithm [1]. On one hand, the geometry is irregular; on the other hand, the combined form of fracture networks is varied. As the properties of fractures are very different from those of the matrix, fractured reservoirs exhibit a high degree of heterogeneity. Describing the geometric features of fractures is challenging because the fractures exhibit various length scales and complex shapes. In a simulation, as the evolution of mechanical unknowns are not as dramatic as that of flow unknowns, the fracture geometry could be simplified in geomechanics simulation, but in flow simulation, the fracture geometry should be modeled explicitly. Traditionally, reservoir simulations use the dual-porosity dual permeability (DPDK) model [2] to model the fracture effects. DPDK treats the fractures and the matrix as two overlapped interacting continua. In a DPDK grid element, fractures are distributed in a predetermined pattern, and the transmissibility between fracture and matrix is derived accordingly. As the fractures are idealized in DPDK, the transmissibility is misrepresented in the case of large-scale fractures, making the application of DPDK limited in fractured reservoir simulation. Different from DPDK, both the discrete fracture model [3] (DFM) and the embedded discrete fracture model [4, 5] (EDFM) are capable of representing real fracture geometries. DFM is the most direct method because the fractures are modeled using grid elements. Flow behaviors in fractures and between fractures and matrix are calculated using the grid transmissibility. The nature of DFM brings both advantages and disadvantages for DFM. The advantages are as follows: (1) the flux calculation method is simple and (2) the accuracy is high. The disadvantage is that the fidelity of fracture representation depends on the gridding technique; currently, the “X” shaped fractures could not be modelled. The limitation of DFM leads the development of the embedded discrete fracture model (EDFM). EDFM models fractures with virtual 2D polygon plates, which are inserted into the background grid. As the computational domains of the fractures and the matrix are different, the complicated gridding is not necessary in EDFM. The flow mass transfer between fractures and matrix is calculated using the non-neighboring connections (NNC). Three types of NNC are required in EDFM; they are as follows: (1) NNC between a fracture cell and its neighboring matrix grid block; (2) NNC between two cells of an individual fracture; and (3) NNC between two intersecting fractures. The detailed procedure using EDFM modeling fracture flow could be found in the work of Li [6].
1.2. Fracture Mechanic Behavior
In reservoir geomechanics simulation, the fracture mechanical deformation is modeled with a tailored constitutive model which is a superposition of fracture and matrix properties. The matrix can be treated as an elastic-plastic block, and the deformation could be calculated according to the elastic-plastic mechanics. However, fracture deformation is a complicated question because of the different stress-strain relationships between normal and shearing deformations [7]. For normal deformation, the constitutive relationship is monotonic and nonlinear. Researchers have proposed different empirical constitutive models [8–10] to describe the relationship. The most commonly used models are Shehata’s model [11] and Bandis’ model [12]. But for shear deformation, the constitutive relationships between shear stress and shear strain are not monotonic due to peak-shear strength [13]. Specifically, the fracture pre-peak deformation is usually assumed to be elastic, whereas the post-peak deformation is plastic [14]. Among the extant models, Barton’s empirical model [15] is the most frequently utilized due to its simplicity and capability of predicting the shearing deformation of fractures with acceptable precision. Moreover, some theoretical models are based on the solid mechanics are proposed, and related works can be seen in Desai [16] and Amadei [17]. The most widely used theoretical model is Plesa’s model [18], which is based on the theory of plasticity and contact mechanics. Currently, existing fracture constitutive models are still immature. Empirical models provide reasonable results in a predetermined range of stresses, but since the principles of energy conservation are not involved in such models, some nonphysical results may be produced when complex mechanical deformation processes are involved; moreover, because the model establishment is not based on a strict mathematical derivation, the robustness of empirical models cannot be guaranteed. Even through the theoretical model is more rigorous, the existing mechanics theory is not capable of modeling fracture behaviors in all aspects.
1.3. Geomechanics Simulation
Because of the importance of reservoir deformation, 3D mechanical calculation is necessary for fractured reservoir simulation. Two types of simulation methods, i.e., the continuum method and the discrete method, can be used. The finite difference method (FDM) discretizes the mechanical equations by replacing the partial derivatives with the numerical differences defined on neighboring grid points directly. Under preassigned initial/boundary conditions, the mechanical deformation can be obtained by solving the discrete equation. However, because of the complex geometry of the reservoir numerical domain, FDM may not be a suitable method. The finite element method (FEM) [19] which was first proposed in the early 1960s offers another choice. The advantages of FEM are its flexibility in modeling reservoir geometries and the convenience of dealing with complicated boundary conditions. To solve a problem, FEM also discretizes the target reservoir into many regular-shaped elements. The difference is that FEM does not discrete the equation directly, but set up the formulation in a weak form, and replaces the initial function with a set of trial functions for each subdomain. The discrete process for a FEM equation follows the law of energy conservation, so its results are more physical. In recent years, to improve the capacity of FEM in fractured reservoir applications, many corresponding developments have been proposed, such as the Generalized/eXtend Finite Element Method (GFEM), which is developed based on the partition of unity principle [20], the Virtual Element Method [21] (VEM), and the Boundary Element Method [22] (BEM). As reservoir simulation is for macroscale problems, the block-type discrete methods are also suitable. The most well-known and representative discrete method is the distinct element method (DEM) [23]. In the original framework of DEM, the solution strategy is explicit. Shi and Goodman further proposed an implicit solution strategy for DEM, which is also named the [24] discontinuous deformation analysis (DDA).
To integrate the geomechanical behaviors, the geomechanics equation and the flow equation should be coupled, and the coupling process would lead to a complex equation system. Three candidate solution schemes exist [25]: the loosely coupled scheme, the fully coupled scheme, and the sequentially coupled scheme. In loosely coupled scheme, the coupling is explicit. The geomechanical equation is solved every certain time step. Although this method saves computational costs compared to other coupling techniques, it is not rigorous and fails to capture some physical phenomena. In contrast, both fully coupled scheme and sequentially coupled scheme are implicit. Fully coupled scheme solves the governing equations of fluid flow and mechanics simultaneously at every time step. It is also unconditionally stable when the coupled problem is well-posed, but the computational cost can be prohibitively expensive. Sequential implicit scheme partitions the coupled problem into two subproblems, and the subproblems could be solved by different methods. The solution is identical to that calculated from the fully coupled scheme, but the computational cost is much less. The fixed stress splitting strategy [26] is an ideal sequential implicit scheme. It solves the flow sub-problem first with fixed volumetric stress. Kim et al. have shown that fixed-stress splitting exhibits excellent convergence properties for both single-phase and multiphase flow problems [27]. From a pragmatic perspective view, fixed stress splitting is preferred to other schemes in geocoupled reservoir simulations.
This paper proposes a new geocoupled model for fractured reservoir simulation. The core idea of this model is that it emphasizes the effects of fractures and geomechanics simultaneously. The remainder of this paper is organized as follows: Section 2 presents the mathematical model for compositional flow and geomechanics. Section 3 describes the fracture model, including the fracture constitutive model, and the method to estimate fracture permeability. Section 4 proposes the fracture gridding technique, which includes the grid system and fracture modeling method. Section 5 first validates the accuracy of this model with two benchmark cases and then presents two large-scale cases to illustrate the applicability of this model. This work is concluded in Section 6 with a discussion of the advantages of this model.
2. Governing Equations
We treat that the solid computational domain and the fluid computational domain are overlapped. The governing equation of fluid flow is set up based on the theory of mass balance, and the governing equation of geomechanics is based on the theory of momentum conservation. Under the assumption of (1) quasi-static, (2) linear elastic, and (3) infinitesimal transformation, the governing equation for geomechanics can be written as where is the divergence operator; is the rank-4 elastic drained bulk modulus; is the strain tensor; is the Biot coefficient; is the equivalent pore pressure; , , and denote the saturation of water, gas, and oil, respectively; , , and denote the pressure of water, gas, and oil phase, respectively; is the rank-2 identity tensor; is the true porosity; is the fluid density; is the rock density; denotes gravity; and is the displacement.
Related studies show that the phase transition behavior of underground fluid [28, 29] is very complex; therefore, we use the fully compositional flow model as the flow governing equations. The mass-balance equation for hydrocarbon component is expressed in terms of the Eulerian derivative as where and are the mole fractions of component in the oil and the gas phase, respectively; is the volumetric strain, which is the trace of the strain tensor ; and are the viscosity of oil and gas, respectively; and are the relative permeability of oil and gas, respectively; is the absolute permeability; and is the mole flow rate of component from the grid to wells and/or boundaries. Owning to the Euler’s form and the infinitesimal transformation, both the flow equation and the geomechanics equation can be discretized on static grids. From the consideration of efficiency and stability, fixed stress splitting and iteration is an ideal solution method for our model. This is also proven in the benchmark tests hereinafter. Li et, al. [3] proposed flow equation in terms of coupling terms for the fixed-stress splitting, which can be expressed as where is the coupling coefficient and subscripts and denote a component and a variable, respectively. For component , they are where .
For the water phase, the coupling coefficients are as follows: where and .
When two hydrocarbon phases coexist, the equation system is not closed. Equal fugacity equations and linear constraints for compositions and saturations should be added to the equation system to constrain the variable set:
Equations (3) to (19) describe the geocoupled compositional flow equation system, where the variable set is . The Peng-Robinson (PR) equation of state and the L.B.C. model are used to calculate the fluid property.
3. Fracture Deformation Model
3.1. Fracture Constitutive Model
We choose the Barton-Bandis model to estimate fracture’s deformation because (1) Barton’s model can model the deformation on normal and shear deformation simultaneously and (2) the experimental stress range of the Barton-Bandis model could cover that of our problem. The definition of a fracture’s compliance follows the work of Bandis et al. [12] and is expressed in matrix form as where the normal and shear stiffness are defined as follows: where is the normal stiffness; is the initial normal stiffness; is the effective normal stress; is the maximum closure; is the shear stiffness; is the maximum shear stiffness; and is the dimensionless compressive strength of rock. and are related to the joint roughness coefficient (JRC), joint compressive strength (JCS), and unstressed joint aperture :
The definition of fracture compliance matrix describes the mechanical property of a fracture in its local coordinate system . To construct a pseudo-continuum, the fracture compliance matrix should be transformed into the global coordinate system [x, y, z], which accords with the coordinates of the solid and fluid zone. Equation (24) constitutes the final formulation of the transformation matrix :
For a fractured rock mass, the compliance is where is the resulting compliance of the fractured rock mass; is the compliance of the intact rock matrix, which is assumed to be linear-elastic and is the inverse of the element stiffness matrix; is the number of fractures; and is the fracture spacing.
3.2. Estimation of Fracture Permeability
In the ideal situation, in which a fracture only consists of two smooth and parallel plates, the hydraulic aperture is equal to the mechanical aperture, and thus the fracture permeability could be evaluated through the cubic law with respect to the hydraulic aperture as follows: where the subscript denotes the reference state; is the fracture permeability; and is the hydraulic aperture. Due to the roughness of the two facing walls, the distribution of local apertures in a fracture is not uniform. For validity of the cubic law, this model employs Barton’s empirical equation (Equation (27)) [30] to convert the mechanical aperture to the hydraulic aperture:
Consequently, the key to updating fracture permeability is estimating the mechanical aperture , which can be expressed in terms of normal and shearing deformations: where is the initial mechanical fracture aperture under zero effective stress; is the change of fracture aperture due to the perturbation of normal effective stress; and is the shear induced dilation. Following Barton [31], can be calculated as follows:
The shear induced dilation is calculated using Asadollahi’s model [32], in which the shear stress-dilation relationship comprises four stages: (1) linear elastic stage; (2) hardening stage; (3) softening stage; and (4) residual stage. The shear induced dilation is formulated as where is the shear displacement of fracture; is the peak shear displacement; and is the dilation displacement. and are both functions of joint roughness coefficient (JRC) and joint compressive strength (JCS):
It should be noted that the fracture permeability model is not unique in the geocoupled simulator. Equations (26) to (32) mostly follow the work of Barton [31], but other choices could exist. One advantage of the pseudo-continuum method is the flexibility in dealing with fracture geometries. Indeed, it is possible to use any empirical or tabular constitutive models for fractures.
4. Fracture Gridding Technique
For reservoirs with preexisting fractures, it is reasonable to treat the fractured block as a continuum, because the mechanical behavior of a fracture constitutes prior knowledge. By keeping the strain energy conservative, a fractured rock block can be represented as a pseudo-continuum [7], which exhibits the same mechanical behavior as the discontinuity. In particular, the mechanical property of the pseudo-continuum (fractured block) is obtained by summing the compliances of fractures and the intact rock.
In the flow problem, EDFM is used to model large-scale fractures, and DPDK is used to model small-scale fractures. The fracture deformation has two influences on the simulation process: the cumulative term calculation and the permeability updating. In the cumulative term calculation, the fracture deformation affects the calculation of material derivative and the volume strain distribution in DPDK model. In the calculation of flow term, the fracture deformation changes the permeability of fractures and then affects the calculation of fracture conductivity. In a geomechanical simulation, the constitutive properties of fractures directly affect the constitutive properties of fractured blocks. Based on the principle of strain energy balance, this section proposes the geomechanical discrete method of fractures, including the construction of equivalent continuous elements and the distribution method of fracture element deformation. The former is used to model the influence of fracture constitutive properties on element properties. The latter is used to distribute the volume strain and calculate the fracture conductivity.
The construction of an equivalent element is based on the conservation of strain energy. For an element, the energy conservation equation is where is the volume and the footnotes f, m, and t denote the fracture, matrix, and the equivalent element, respectively. is the strain of the equivalent element; is the stress tensor of the equivalent element; and , and are the compliance of equivalent element, fracture, matrix, respectively. Hence, the can be calculated as where and represent large scale and small scale, respectively. It can be seen that the flexibility of the quasi-continuous element is equal to the weighted sum of the volume of the crack and the matrix. This equation assumes that the permeability of fracture and matrix are described in the same coordinate system, but in reality, the mechanical properties of fracture and matrix are often established in different coordinate systems, so it is necessary to combine the coordinate systems of the two. According to the space conversion principle of vector, the coordinate conversion matrix is where is the angle between local axis and global axis . Then, we can get the construction equation for the equivalent element: where is the number of large fractures in the element and is the number of small fractures. In DPDK, small fractures are assumed to be uniformly distributed in the element at the same angle; they are equivalently treated as a change in porosity. In this case, the second item at the right end of the equation can be replaced by the following formula:
In flow simulation, fracture and matrix are calculated separately. Therefore, it is necessary to decompose the deformation for fracture and matrix according to the deformation of the equivalent element, respectively. The deformation decomposition has two functions in the simulation. The first is to provide volume strain for the generation of the cumulative term in DPDK model. The second is to provide the stress state of fractures for the calculation of fracture permeability, so as to calculate the opening of fractures. The deformation distribution of equivalent elements is also based on the strain energy balance theory. The total deformation of the equivalent continuous element is the sum of the deformation of fracture and matrix:
According to the energy conservation theory, the strain energy of an equivalent element meets the following condition in 1D situation:
Then, the strain energy balance equation can be extended to
It can be concluded that when the quasi-continuous volume strain is known, the volume strain of crack and matrix are calculated as follows:
It is obvious that
Therefore, the strain of crack and matrix can be calculated as follows:
The stress calculation method under three-dimensional state of a fracture is where is the stiffness matrix of a fractured element in current state.
The governing equation system is linearized through the Newton-Raphson method. The coupling problem is solved through the fixed-stress splitting strategy, which solve the flow problem first with a fixed volumetric stress in one iteration step. The fracture information is coupled with the governing equation system explicitly and is updated at the beginning of each time step. In each coupling step, the flow problem is solved first with the volumetric stress fixed; then, the geomechanical problem is solved with a frozen pressure. In the intermediate step, the displacement does not need to be calculated, because the mechanics problem is quasi-static.
5. Numerical Results
We show the validity and applicability of this model through some cases in this section. Case 1 is the Mandel’s problem, which is widely used to validate the accuracy of geocoupled simulators. Case 2 is a reservoir with two fractures, which is set up to validate the accuracy of this method in dealing with fracture problems. Case 3 and Case 4 are two field-scale cases; they are set up to demonstrate the applicability of this method.
Case 1. Mandel’s problem. The physical model of this case is an infinitely long rectangular specimen sandwiched between two frictionless, rigid plates. Two lateral sides of the specimen are drained and traction-free. Initially, a force acts on the top surface, causing the Mandel-Cryer effect. The geometry set for Mandel’s problem refers to Figure 1. The benchmark solution can be found in Goulet’ work [33], and the numerical model is built using the parameters listed in Table 1. Because of the symmetry of the specimen, we only generate the numerical domain for a quarter of the specimen. The numerical solutions for pore pressure along the -direction at various times are shown in Figure 2. It can be seen that the Mandel-Cryer effect is captured by this model. The solution from this model and the analytical solution are in good agreement.
Case 2. Comparison with CMG. Case 2 is a water-flooding problem in a rectangular reservoir. The reservoir size is . The traction is set to 32.9 MPa to face I+, J+, and K+, and the fixed displacement boundary condition is acted to face I-, J-, and K-. All sides are no-flow boundaries. The physical properties are listed in Table 1. Initially, the fluid pressure is 300 Bar, and the reservoir is oil-saturated. The grids are generated as shown in Figure 3. We set the initial oil production rate to 1589 m3/d and the water injection rate to 476.96 m3/d. We choose the CMG [34] as the benchmark simulator, which employs the explicitly coupled method to model the geocoupled problems. Figure 4 shows the comparison of the oil production rate and the water production rate. The result is outputted once when the solution reaches convergence. When the solution is not converged, the time step would be reduced, and the simulation would be rolled back. It can be seen that the results calculated by this model are very close to those of CMG, but because of the implicitly coupled strategy, the convergence of this method is better.
Case 3. Tight oil production.
Case 3 is a tight oil production problem, which contains a horizontal well that is 15 m long, with eight hydraulic fractures striking at different orientations. The stages 2, 4, and 6 are closed. The fractures are modeled using the EDFM method. The reservoir is discretized into 15084 grid cells, as shown in Figure 5. The geomechanical domain and the reservoir domain are identical. The top surface of the over-burden is traction-free, and the bottom of the under-burden is fixed. The vertical stress gradient is 20.5 kPa/m. The ratio between the maximum horizontal stress to the vertical stress is 0.9, and the ratio between the minimum horizontal stress to the vertical stress is 0.8. Initially, the fluid pressure is 155 Bar. Other physical properties of the case are listed in Table 2. The well is working under a constant bottom hole pressure (BHP) of 50 Bar.
Figure 6 shows the simulation results; the simulation results of the initial state, the 1440th day, and the 2550th day are shown from left to right, from which we can see that the geomechanics affects the reservoir to a considerable degree. As the major flow channel, the fractures are critical to production. The volumetric strain and fluid pressure evolve mostly around the stimulated reservoir volume. Because of the press evolve, the difference between the maximum and minimum principal stress changes may affect the subsequent re-fracturing effect.
Case 4. Tight oil production: CO2 Huff-N-Puff.
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
(i)
Case 4 is designed to show the capability of this method in simulating the coupled geomechanics and compositional flow process in the presence of phase transition. As shown in Figure 7, the flow domain is a fractured reservoir with one horizontal well and two hydraulic fractures. Initially, the reservoir is oil-saturated, and the fluid composition is shown in Table 3. The reservoir fluid consists of fourteen pseudo hydrocarbon components, as shown in Table 2. The initial reservoir pressure is 250 Bar. The boundary conditions and physical properties used in this example are the same as those in Case 3. The schedule is from 0 to 1000 days, and the well is set to production with the bore pressure of 200 Bar; then, we start the Huff-N-Puff operation. The operation is implemented for 5 cycles; each cycle is 90 days, of which the first 30 days are the gas injection stage; and the daily gas injection is 800 m3. After that, soak the well for 20 days and open for 40 days. Figure 8 compares the cumulative oil production under different conditions. It can be seen from the figure that the effect of geomechanics predicted from the geocoupled model is lower than that of the component model. After 1000 days of production, the production capacity of production wells was seriously insufficient. After the implementation of the Huff-N-Puff operation, the production capacity was effectively improved and the final cumulative production was increased by about 7.7%.
The main mechanism of the Huff-N-Puff operation is that it improves the fluidity of fluid by injecting light components while increasing the reservoir pressure. In the first 1000 days, the composition of the fluid is evenly distributed. The heterogeneous distribution is caused by the injection of CO2. Therefore, the reconstruction degree of reservoir fluid can be tracked according to the molar concentration of CO2. Figure 9 shows the carbon dioxide distribution simulated from the component model and the geocoupled model, respectively. From the fluid geocoupled model, it can be seen that after the CO2 injection, because the transverse permeability is much larger than the longitudinal permeability, it first improves the fluid of the perforated layer. After well soaking, due to its low molar mass, the component distribution is affected by buoyancy, resulting in its funnel-shaped distribution in the latter situation. In the geomechanical model, due to fracture closure behavior, the permeability of the fractured reservoir is low, and the fluid pressure is high. Therefore, the CO2 distribution predicted in the geocoupled model is smaller than that in the component model, and the concentration of CO2 in the distribution area is higher than that simulated in the component model.
6. Conclusion
In this work, a new geocoupled fractured reservoir model is proposed, which has the following advantages: (1)Modeling the dynamic geomechanics field. Because of the implicit fully coupled model, the geomechanics field can be calculated and outputted in every time step. The geomechanics numerical field and the reservoir numerical field are identical(2)Modeling the fracture deformation behaviors. The fracture treatment is a systematic work. In this work, the fracture mechanical behaviors are modeled through Barton’s model, and a novel fracture gridding technique, which is based on strain energy conservation, is proposed to integrate the fracture model into the simulation
We presented some numerical cases to demonstrate the accuracy and applicability of this model. Moreover, numerical cases show that this method is capable of dealing with some complicated problems, such as simulating the geocoupling process in the case of EDFM and DFM and coupling compositional flow and geomechanics. The practicality of this method makes it applicable for more complicated field problems.
Data Availability
The data used to support the findings of this study are available from the corresponding author upon request.
Conflicts of Interest
The authors declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.
Acknowledgments
This work is partially funded by the National Natural Science Foundation of China (Grant No. 51974356).