© T. Li et al., published by IFP Energies nouvelles, 2020

1 Introduction

Accurate and effective simulation of coupled flow and geomechanics (also known as coupled Hydro-Mechanical [HM]) plays a vital role in studying various geophysical and engineering problems, such as hydrocarbon reservoirs exploitation, CO₂ sequestration, and geothermal extraction, etc. [1–4]. Extensive theoretical and numerical study has been carried out for coupled flow and mechanical problems in the past several decades. Terzaghi first proposed the coupled process of flow and mechanics in geological media as a consolidation phenomenon [5]. Subsequently, based on the linear-elastic theory and Darcy’s law, the Three-Dimensional (3D) consolidation theory was established by Biot [6]. One of the first solutions for the 3D consolidation theory was presented and the non-monotonic pore-water pressure response was demonstrated by Mandel [7]. For some special cases, the analytic method is simple and effective [5–7]. However, for complex problems, numerical simulation is often commonly used [1–4, 8–11].

In terms of numerical simulation, there are two different families of approaches for discretization and solution of governing equations, which are Finite Element Method (FEM) and Finite Volume Method (FVM). FEM is mainly applied to solve the problems of solid mechanics deformation [8, 9]. FVM is mainly used to solve the problems of flow and heat transfer [10, 11]. Therefore, the mixed solution method (FVM + FEM) is generally used to solve the fluid and mechanical coupling problems in subsurface porous media [12]. Now, the FVM has been extended to the field of solid mechanics and achieved great success [13–16]. Suliman et al. [17] proposed an enhanced FVM for the purpose of modeling linear elastic structures undergoing bending. Compared with the traditional FEM, the advantage of the model was verified and a rigorous error analysis was conducted. Asadi et al. [13] solved a Biot consolidation model with discontinuous coefficients by FVM. The results demonstrated that the finite volume scheme led to a local mass conservation method, which can remove pressure oscillations and yielded higher accuracy for the flow and mechanics parameters. Nordbotten [14] proposed a cell-centered FVM to discretize coupled flow and deformation of porous media. The presented approach was validated by comparison to analytical solutions. Based on the Biot’s consolidation theory and Mohr–Coulomb constitutive relation, Tang et al. [15] presented a new FVM code for poro-elasto-plasticity soil model. The results showed that the implemented solver was capable and efficient at predicting reasonable soil responses with pore pressure coupling under different boundary conditions.

In addition to the deformation of the rock mass, one of the difficulties that cannot be ignored in the simulation of fluid mechanical coupling in porous media is the fracture system. In general, the permeability and porosity of underground porous media are very low, and it is difficult for the working fluid to pass through without hydraulic fracturing or chemical stimulation. Therefore, characterization and simulation of the fracture system is also a key issue for better understanding of flow and mechanical behavior in deformable porous media. In fact, the traditional FEM can deal with the flow and deformation of fractured porous media. However, this will introduce a large number of unstructured grids. To avoid the generation of the complex grids, the Extended Finite Element Method (XFEM) is proposed, which can allow for using non-conforming grids to discrete the computational domain [18–21]. Yan et al. [20] developed a stabilized XFEM scheme to address the displacement oscillation along macro-fracture boundaries. The mixed space discretization and modified fixed stress sequential implicit methods were used to solve the coupling model. The model and results were verified by comparison to COMSOL. Shi et al. [21] proposed an XFEM-based method for modeling hydraulic fracture propagation. To save the CPU time for solving large-scale linear equation systems, a reduction technique was adopted in a tightly coupled model.

In addition to FEM and XFEM, many efforts have been paid for establishing a finite volume-based simulation method for fracture deformation [22–26]. Based on the cell-centered finite-volume approach, Ucar et al. [22] proposed a multi-point stress approximation method to discretize coupled flow and mechanical deformation in fractured porous media. The convergence of the method was examined for several benchmark problems. Recently, Deb et al. [23–26] presented an Extended Finite Volume Method (XFVM) for modeling of shear failure in fractured reservoirs. Similar to the XFEM, the enrichment functions were introduced to simulate the displacement discontinuity of fracture surfaces.

From the above review, according to grid type and algorithm, we can find that there are three main methods to solve the problem of coupled fluid and geomechanics. The first one is, based on the unstructured grids, the flow and mechanical equations are solved together by standard FEM or FVM. The second one is the mixed solving technology, which is the flow equations are solved by the FVM, and the mechanical equations are solved by the XFEM. By doing this, the computational domain can be discretized by structure grids. The last one is the flow field is still obtained by the FVM, and the deformation field is calculated by the XFVM. The advantage of this method is that we can still use structured grid generation technology.

Although a large number of models and algorithms have been proposed for the simulation of fluid and mechanical coupling process, there are no quantitative comparisons among them. It is very important for scientific research to know the advantages and disadvantages of each algorithm and model. Therefore, in this paper, we will investigate the performance of XFEM and XFVM on the simulation of flow-induced fracture deformation in porous media. To avoid the influence of the fluid field on these two methods, the flow equations are discretized and solved by the standard FVM. Then, the fluid pressure is provided to calculate the mechanical deformation.

The outline of this paper is as follows. The mathematical models including flow equation and mechanical equilibrium equation are presented in Section 2. The numerical implementation and discretization of those equations by traditional FVM and extended FEM and FVM are given in Section 3. The numerical comparisons for three test cases and error estimation are presented in Section 4. Finally, the conclusions and remarks are given in Section 5.

2 Mathematical models

2.1 Flow equations

For the sake of simplicity, we only consider the single-phase flow problem in the fractured porous media. Combined with Darcy’s law, the mass balance equation can be obtained as follow:(1)where t is the time; ρ is the fluid density; ϕ is the porosity; is the permeability tensor; μ is the fluid viscosity; p is the fluid pressure; g is the gravitational acceleration; z is the vertical distance; q is the source term.

Considering the slight compressibility of rock and fluid, the first term of LHS in equation (1) needs to be further deduced:(2)

Introducing the compressible coefficients of rock and fluid, , , equation (2) reads:(3)

Assuming small spatial variations for the density, i.e. , therefore, the final expression of equation (1) can be obtained:(4)

In order to consider the influence of fractures on the fluid flow, the Embedded Discrete Fracture Model (EDFM) is adopted [27, 28]. According to the implementation process of EDFM, the fractured media can be distinctly separated into a fracture system and a damaged matrix [29]. Therefore, the fracture and matrix are computationally independent except for the transfer function. Besides, this model ignores the gradient change of fracture pressure in the normal direction, which allows for a lower-dimensional representation of the fracture system [30].

Based on the EDFM, the flow equation expressed by equation (4) can be separated into two parts:

• Matrix:

(5)

• Fracture:

(6)where superscripts m and fr represent matrix and fracture media, respectively; and are the transfer functions [31].

2.2 Mechanical equations

The balance of linear momentum in the porous matrix is given by:(7)where is the total stress; b is the body force vector.

Note that we use the convention that tensile stress is positive. Relation between the effective stress and the pore pressure is given as:(8)where is the total stress; is the effective stress; is the identical tensor; α is the Biot coefficient.

Relation between the stress and the stain is defined as:(9)where is the volumetric strain; λ and G are the Lame’s coefficients.

Relation between the stain and the displacement is given as:(10)(11) where is the displacement vector, .

According to Cauchy stress formula, the tensile and shear stress on the fracture surface are given as:(12)(13)where n represents the normal vector to the fracture surface.

3 Numerical implementation and basic discretization

For the flow equations, the classical cell-centered Finite Volume Method (FVM) is used. The Two-Point Flux Approximation (TPFA) scheme for pressure gradient term is adopted. The implicit scheme for approximation of time is used. For the mechanical equations, we used two solving algorithms, which are XFEM and XFVM.

3.1 Discretization of flow equations

Consider the Two-Dimensional (2D) domain shown in Figure 1. We have the following approximation for the flow equations:(14)(15)where and represent the control volume of matrix and fracture segment, respectively; , , , and represent the transmissivity of matrix-matrix, matrix-fracture, and fracture-fracture. and represent the matrix and fracture pressure calculated at the last time step.

Obviously, equations (14) and (15) need to be solved in coupling. We can further reshape these two equations into a matrix form [32]:(16)where superscript w represents the external well terms, which can be implemented by the Peaceman formula [33]. If no well is drilled into the matrix and fracture domain, the ‘w’ terms will be zero.

3.2 Discretization of mechanical equations

3.2.1 Extended Finite Element Method (XFEM)

To derive the XFEM formulation for an equilibrium body with a strong discontinuity, consider a 2D domain with an internal fracture denoted by shown in Figure 2. The essential boundary conditions are defined as follows:

1. The Dirichlet boundary condition is defined as on , where is the prescribed displacement on the boundary .

2. The Neumann boundary conditions are defined as on , where is the prescribed stress on the boundary , and . For the free boundary, we have . Moreover, the internal boundary condition is defined as .

Therefore, the weak form of equilibrium equation (7) can be obtained multiplying by the test functions and integrating over the domain as:(17)where fracture aperture variation .

For XFEM, to discrete the integral equation (17), we need to add some special enrichment functions to the conventional finite element shape functions. Based on the type of discontinuity and the location of fractures, the most commonly used enrichment functions for the strong discontinuity are Heaviside step function H(x), tip branch function and Junction function [19–21]. Therefore, the enriched displacement field can be defined as(18)where , , are the node enrichment variables.

By substituting the XFEM approximation (Eq. (18)) into the weak form of the equilibrium equation (17), the discrete system can be finally obtained:(19)where K represents the global stiffness matrix, which is assembled by the element stiffness :(20)

The strain matrix B is defined as follows:(20a)(20b)(20c)(20d)

In equation (19), R is the global external force vector, which is assembled by the element force vector :(21)The node force vector can be defined as follows:(21a)(21b)(21c)(21d)

3.2.2 Extended Finite Volume Method (XFVM)

Similar to the XFEM, the element displacement field of the XFVM is defined as [26]:(22)where and represent the normal and shear displacement vector on the fracture surface.

As shown in Figure 1, the integration of equation (7) over the domain can be given by(23)

Using the divergence theorem, equation (23) can be converted to the integration over the boundary :(24)

For convenience, the equation (24) can be further converted to the form of stress components in x and y-directions:(25a)(25b)

Combining the equations (9)–(11), we can obtain the force balance equations in the form of displacement:(26a)(26b)

The key idea of XFVM is how to calculate the integral displacement derivative in equations (26a) and (26b). Integration of any function F along the internal interface is defined as:(27)

As shown in equation (22), the contributions of integral displacement derivative of equation (27) consist of two parts, the continuous and discontinuous degrees of freedom:(28)where,(28a)(28b)(28c)

Due to the equation (28) introducing the additional degrees of freedom, we need to couple the force constraint equations of fractures, as shown in equations (12) and (13). Therefore, integrating these two equations over a fracture segment L, we can obtain the following force constraint equations in terms of displacement:(29)(30)

Similar to the equation (28), the contributions of integral displacement derivative still include two parts:(31)where,(31a)(31b)(31c)

For the fracture intersection, we need to add an interaction term between interacting fracture segments. For details, one can refer to [26]. Finally, we can obtain the discrete form of force equilibrium equation:(32)where and represent the coefficient matrix; is the force vector.

4 Numerical examples and results

Three test cases are designed to study the accuracy and efficiency of the XFEM and XFVM. The computational domain and boundary conditions are shown in Figure 3. To study the accuracy of these three cases, the numerical results of XFEM and XFVM are compared to the commercial software (COMSOL), where a standard FEM solver is implemented. Then, the pattern of discrete sparse matrix and computational efficiency for the XFEM and XFVM are discussed.

4.1 Single fracture

In this test case, a single horizontal fracture is predefined in the computational domain, as shown in Figure 3a. The length of the square matrix and fracture is 100 m and 4 m, respectively. The other computational parameters are given in Table 1. The corresponding boundary conditions are as follows.

1. Pressure boundary conditions: the pressure on the left and right boundaries is 10⁷ Pa and 10⁵ Pa, respectively; the no-flow boundary condition is applied on the top and bottom sides; the initial pressure is zero.

2. Displacement boundary conditions: the fixed boundary conditions in the x and y-directions are applied on the four sides.

Figure 4a shows the free triangle grids generated by the COMSOL. The total numbers of the domain and boundary elements are 11 400 and 352, respectively. Figure 5a shows the structure grids adopted in XFEM and XFVM. There are 101 × 101 mesh points in the matrix area and 43 grids for fracture. The pressure distributions obtained by the COMSOL and FVM are depicted in Figure 6. First, because of the stronger conductivity of the fracture, the internal fluid velocity is faster and the pressure gradient is smaller than that of the matrix. Second, we can see that the two groups of results are in good agreement in both matrix and fracture regions. Therefore, the solution of the flow equation based on the FVM is first verified. Figure 7 shows the comparison of displacement fields between XFEM and XFVM. The results of displacement in x and y-directions for XFEM and XFVM are very close to that of COMSOL. It can be observed from the y-displacement that the fracture surface shows a trend of closing under the action of pore pressure. In order to compare the difference between them quantitatively, we define a Maximum Relative Error (MRE)(33)where represents the results obtained by COMSOL; denotes the results calculated by XFEM and XFVM.

Figure 8 shows the displacement comparisons of fracture upper surface among these three groups of results. According to the equation (33), we can obtain the MRE along the fracture surface is 4.54% between XFEM and COMSOL, and the MRE between XFVM and COMSOL is 3.41%.

4.2 Two intersecting fractures

As shown in Figure 3b, a cross fracture is considered in the domain. These two fractures are 40 m in length and the aperture is 10⁻³ m. The intersection coordinate is (50 m, 50 m). The injection well is located at the bottom left corner, and the injection pressure is 10⁷ Pa. The production well is located in the upper right corner, and the production pressure is 10⁵ Pa. The no-flow boundary conditions are applied on the four sides. The left and right boundaries are fixed in the x-direction. The upper and bottom boundaries are fixed in the y-direction. Initial displacement and fluid pressure are zero. The other computational parameters are given in Table 1. As shown in Figure 4b, the total number of domain elements employed in COMSOL is 16 130, and the total number of boundary elements is 504. The total number of matrix and fracture grids is 10 201 and 84, as shown in Figure 5b.

The pressure field distributions and comparisons are depicted in Figure 9. We can see that, the pressure gradient is larger near injection well and production well, and smaller near fractures. The results of FVM are in good agreement with that of COMSOL. Figure 10 shows the comparison of the results of displacement in x and y-directions between COMSOL, XFEM, and XFVM for this cross fracture. We can see that, under the action of pore pressure, the displacement fields in x and y-directions show a trend of symmetrical distribution. Moreover, the solution of XFEM and XFVM is almost the same as that of COMSOL. Figure 11 shows the comparisons of x-direction displacement on the upper surface of horizontal fracture among COMSOL, XFEM, and XFVM. As can be seen from Figure 11, the results of displacement fields obtained by XFEM and XFVM are similar, but the results of XFEM are closer to that COMSOL than XFVM. The MRE between XFEM and COMSOL is 7.97%, and for the XFVM, it is 9.67%.

4.3 Eight discrete fractures

In this test case, we consider eight discrete fractures included in a 2D porous media. The pressure and displacement boundary conditions are the same with case 2. The computational grids adopted in COMSOL are shown in Figure 3c. The total number of elements is 38 914, and the total degrees of freedom is 98 292. On the other hand, the grids used in XFEM and XFVM are shown in Figure 4c. The number of matrix grids is the same with case 1 and case 2. The total numbers of fracture grids are 232. The other related parameters are presented in Table 1. Figure 12 shows the comparison of the results of fluid pressure between COMSOL and FVM. Because of the existence of discrete fractures, the distribution of the pressure field presents certain anisotropy. It can be seen that these two groups of results are very close. Figure 13 shows the comparison of displacement fields among COMSOL, XFEM, and FVM. Because the distribution of fracture is not very regular, the distribution of the displacement field also presents anisotropy. At the same time, we can clearly view that the results of XFEM and XFVM are almost the same as the reference solutions. We choose the horizontal fracture at the bottom of the computational area for accuracy comparison, as shown in Figure 14. Compared with COMSOL, the MRE for the XFEM is 3.03% and 7.78% for XFVM.

In order to study the sparsity of the discrete matrix obtained by these two methods, we define the following condition numbers.(34)where represents the 2-norm of matrix K.

We can see that from Figure 15, as the number of fractures increases, the increase of variables of non diagonal block in matrix. Although the number of variables of the sparse matrix obtained by XFEM is more than the XFVM, the condition number of the sparse matrix obtained by XFVM is larger than the XFEM, as shown in Table 2. The condition number of the sparse matrix in cross fractures case is larger than the uncrossed cases. The effect could be explained by the fact that more degrees of freedom need to be introduced into the cross element for XFEM, and the influence between different fracture segments needs to be considered for XFVM. Table 3 shows the CPU times for these three test cases. We can see that XFVM is more efficient than XFEM. The reason is that XFEM mainly calculates the integrand function by numerical integration shown in equations (20) and (21). However, XFVM can accurately calculate the integral value of the integrand function at the control volume interface given in equations (28) and (31).

5 Conclusion

In this study, we discuss two numerical methods for simulating flow-induced fracture deformation in porous media. First, the flow equations are established based on the embedded discrete fracture method and solved by using the traditional finite volume method. Then, the mechanical equilibrium equations are discretized by two methods, XFEM and XFVM. Through three test examples, we study the accuracy and efficiency of these two methods in comparison with the standard finite element method implemented in COMSOL. The main conclusions and remarks can be drawn:

1. In the first test case, the accuracy of the numerical calculation of XFEM and XFVM is not much different. In the following two test examples, the accuracy of XFEM is slightly higher than that of XFVM. Compared to the COMSOL, the maximum relative error among these three cases for XFEM is 7.97% and 9.67% for XFVM. Similar to the XFEM, XFVM also introduces the enrichment function to simulate the displacement discontinuity of the fracture surface. However, this method ignores the influence of the fracture tip, so it may lose some accuracy.

2. Although the condition number of the discrete matrix obtained by XFVM is much larger than that obtained by XFEM, its computational efficiency is 4–5 times higher than that of XFEM. There are two main reasons. One is XFVM obtains the integral value of the integrated function by direct offline integration, while XFEM obtains that by online Gauss integration, which will consume a lot of CPU time. The other is XFEM needs to add a lot of extra degrees of freedom (at least four additional degrees of freedom per element) to elements with fractures passing through. However, for the XFVM, only two additional degrees of freedom (tangential and normal displacement components of fractures) need to be introduced to the damaged elements.

3. This study is only a preliminary discussion of these two methods. There are still many key challenges for XFEM and XFVM, which should be addressed in the numerical simulation of flow-induced deformation problems. For example, how to more efficiently simulate the deformation of complex fractures system for the XFEM and keep the discrete matrix not ill-conditioned, and how to more accurately simulate the deformation of cross fractures for the XFVM under the action of external force.

Acknowledgments

The study is supported by the National Natural Science Foundation of China (Nos. 51936001 and 51706021), the Project of Construction of Innovative Teams and Teacher Career Development for Universities and Colleges under Beijing Municipality (No. IDHT20170507), the Program of Great Wall Scholar (No. CIT&TCD20180313), and the Beijing Youth Talent Support Program (No. CIT&TCD201804037).

References

All Tables

All Figures