Regular Article
A comparative study on simulating flowinduced fracture deformation in subsurface media by means of extended FEM and FVM
^{1}
School of Chemical Engineering and Technology, Xi’an Jiaotong University, 710049 Xi’an, PR China
^{2}
School of Mechanical Engineering, Beijing Key Laboratory of Pipeline Critical Technology and Equipment for Deepwater Oil & Gas Development, Beijing Institute of Petrochemical Technology, 102617 Beijing, PR China
^{*} Corresponding authors: yang.fs@xjtu.edu.cn; yubobox@vip.163.com
Received:
25
January
2020
Accepted:
15
May
2020
Accurate and efficient simulation on the fluid flow and deformation in porous media is of increasing importance in a diverse range of engineering fields. At present, there are only several methods can be used to simulate the deformation of fractured porous media. It is very important to know their application scopes, advantages, and disadvantages for solving the practical problems correctly. Therefore, in this paper, we compared two numerical simulation methods for flowinduced fracture deformation in porous media. One is the Extended Finite Element Method (XFEM), which is based on the classical finite element method and can simulate strong or weak discontinuous problems. The other is developed within the finitevolume framework, termed Extended Finite Volume Method (XFVM). We designed three test cases, including single fracture, cross fractures and eight discrete fractures, to investigate the accuracy and efficiency of XFEM and XFVM. The reference solutions were provided by the commercial software, COMSOL, where the standard finite element method is implemented. The research findings showed that the accuracy of the XFEM was slightly higher than that of the XFVM, but the latter was more efficient. These results are likely to be useful in decision making regarding choice of solving methods for the multifield coupling problem in fractured porous media.
© T. Li et al., published by IFP Energies nouvelles, 2020
This is an Open Access article distributed under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
1 Introduction
Accurate and effective simulation of coupled flow and geomechanics (also known as coupled HydroMechanical [HM]) plays a vital role in studying various geophysical and engineering problems, such as hydrocarbon reservoirs exploitation, CO_{2} 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 linearelastic theory and Darcy’s law, the ThreeDimensional (3D) consolidation theory was established by Biot [6]. One of the first solutions for the 3D consolidation theory was presented and the nonmonotonic porewater 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 cellcentered 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 poroelastoplasticity 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 nonconforming grids to discrete the computational domain [18–21]. Yan et al. [20] developed a stabilized XFEM scheme to address the displacement oscillation along macrofracture 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 XFEMbased method for modeling hydraulic fracture propagation. To save the CPU time for solving largescale 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 volumebased simulation method for fracture deformation [22–26]. Based on the cellcentered finitevolume approach, Ucar et al. [22] proposed a multipoint 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 flowinduced 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 singlephase 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 lowerdimensional representation of the fracture system [30].
Based on the EDFM, the flow equation expressed by equation (4) can be separated into two parts:
Matrix:
Fracture:
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 cellcentered Finite Volume Method (FVM) is used. The TwoPoint 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 TwoDimensional (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 matrixmatrix, matrixfracture, and fracturefracture. and represent the matrix and fracture pressure calculated at the last time step.
Fig. 1 The schematic of the Embedded Discrete Fracture Model (EDFM). The red points represent the matrix grids, which storage fluid pressure; the blue points denote the fracture grids, which storage fracture pressure and displacement for the XFVM; the bright green points storage the matrix displacement for the XFVM and XFEM; e, w, n, s represent the interface of control volume; E, W, N, S represent the internal interface in a control volume. 
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:
The Dirichlet boundary condition is defined as on , where is the prescribed displacement on the boundary .
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 .
Fig. 2 A twodimensional elastic body with an internal fracture. 
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 ydirections:(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.
Fig. 3 Computational domain and boundary conditions. The red lines represent the discrete fractures. The red point at the left corner denotes the injection well, and the blue points at the right corner represent the projection well. (a) Single fracture; (b) cross fractures; (c) eight discrete fractures. 
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.
Pressure boundary conditions: the pressure on the left and right boundaries is 10^{7} Pa and 10^{5} Pa, respectively; the noflow boundary condition is applied on the top and bottom sides; the initial pressure is zero.
Displacement boundary conditions: the fixed boundary conditions in the x and ydirections are applied on the four sides.
Computational parameters for these three test cases.
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 ydirections for XFEM and XFVM are very close to that of COMSOL. It can be observed from the ydisplacement 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.
Fig. 4 Computational grids used in the COMSOL. (a) 11 400 elements; (b) 16 130 elements; (c) 38 914 elements. 
Fig. 5 Computational grids used in XFEM and XFVM. (a) Matrix: 10 201, fracture: 43; (b) matrix: 10 201, fracture; 84; (c) matrix: 10 201, fracture: 232. 
Fig. 6 Results comparison of fluid pressure between COMSOL and FVM. 
Fig. 7 Comparisons of displacement fields among COMSOL, XFEM, and XFVM for single fracture. The first row represents the displacement fields in the xdirection, the second row represents the displacement fields in the ydirection. Left: results of COMSOL; middle: results of XFEM; right: results of 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%.
Fig. 8 Comparison of the ydirection displacement field on the fracture upper surface among COMSOL, XFEM, and XFVM. 
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^{−3} 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^{7} Pa. The production well is located in the upper right corner, and the production pressure is 10^{5} Pa. The noflow boundary conditions are applied on the four sides. The left and right boundaries are fixed in the xdirection. The upper and bottom boundaries are fixed in the ydirection. 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 ydirections 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 ydirections 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 xdirection 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%.
Fig. 9 Results comparison of fluid pressure between COMSOL and FVM. 
Fig. 10 Comparisons of displacement fields among COMSOL, XFEM, and FVM for cross fractures. The first row represents the displacement fields in the xdirection, the second row represents the displacement fields in the ydirection. Left: results of COMSOL; middle: results of XFEM; right: results of XFVM. 
Fig. 11 Comparisons of xdirection displacement field on the upper surface of horizontal fracture among COMSOL, XFEM, and XFVM. 
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.
Fig. 12 Results comparison of fluid pressure between COMSOL and FVM. 
Fig. 13 Comparison of displacement fields among COMSOL, XFEM, and FVM for complex fractures. The first row represents the displacement fields in the xdirection; the second row represents the displacement fields in the ydirection. Left: results of COMSOL; middle: results of XFEM; right: results of XFVM. 
Fig. 14 Comparisons of ydirection displacement field on the lower surface of horizontal fracture at the bottom of the domain among COMSOL, XFEM, and 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 2norm 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).
Fig. 15 Sparsity matrix pattern for the three cases. Left: XFEM; right: XFVM. 
Condition number of the discrete matrix for XFEM and XFVM.
The CPU time for the XFEM and XFVM.
5 Conclusion
In this study, we discuss two numerical methods for simulating flowinduced 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:
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.
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.
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 flowinduced deformation problems. For example, how to more efficiently simulate the deformation of complex fractures system for the XFEM and keep the discrete matrix not illconditioned, 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
 Kim J., Tchelepi H.A., Juanes R. (2011) Stability and convergence of sequential methods for coupled flow and geomechanics: Fixedstress and fixedstrain splits, Comput. Meth. Appl. Mech. Eng. 200, 13, 1591–1606. [CrossRef] [Google Scholar]
 Huang Z.Q., Winterfeld P.H., Xiong Y., Wu Y.S., Yao J. (2015) Parallel simulation of fullycoupled thermalhydromechanical processes in CO_{2} leakage through fluiddriven fracture zones, Int. J. Greenhouse Gas Control 34, 39–51. [CrossRef] [Google Scholar]
 Li S., Li X., Zhang D. (2016) A fully coupled thermohydromechanical, threedimensional model for hydraulic stimulation treatments, J. Nat. Gas Sci. Eng. 34, 64–84. [Google Scholar]
 Jin L., Zoback M.D. (2017) Fully coupled nonlinear fluid flow and poroelasticity in arbitrarily fractured porous media: A hybriddimensional computational model, J. Geophys. Res.: Solid Earth 122, 10, 7626–7658. [CrossRef] [Google Scholar]
 Terzaghi K. (1943) Theoretical soil mechanics, John Wiley and Sons, New York, NY. [Google Scholar]
 Biot M.A. (1941) General theory of threedimensional consolidation, J. Appl. Phys. 12, 2, 155. [Google Scholar]
 Cheng H.D., Abousleiman Y., Detournay E., Cui L., Roegiers J.C. (1996) Mandel’s problem revisited, Géotechnique 46, 2, 187–195. [CrossRef] [Google Scholar]
 Murad M.A., Loula A.F.D. (1994) On stability and convergence of finite element approximations of biot’s consolidation problem, Int. J. Numer. Meth. Eng. 37, 4, 645–667. [CrossRef] [MathSciNet] [Google Scholar]
 Jha B., Juanes R. (2007) A locally conservative finite element framework for the simulation of coupled flow and reservoir geomechanics, Acta Geotech. 2, 3, 139–153. [CrossRef] [Google Scholar]
 Deng Y., Sun D., Liang Y., Yu B., Wang P. (2017) Implementation of the ideal algorithm for complex steadystate incompressible fluid flow problems in openfoam, Int. Commun. Heat Mass Trans. 88, 63–73. [CrossRef] [Google Scholar]
 Li T.Y., Gao Y.Q., Han D.X., Yang F.S., Yu B. (2020) A novel POD reducedorder model based on EDFM for steadystate and transient heat transfer in fractured geothermal reservoir, Int. J. Heat Mass Trans. 146, 2020, 118783. [CrossRef] [Google Scholar]
 Vasilyeva M., Chung E.T., Efendiev Y., Kim J. (2018) Constrained energy minimization based upscaling for coupled flow and mechanics, J. Comput. Phys. 376, 660–674. doi: 10.1016/j.jcp.2018.09.054. [Google Scholar]
 Asadi R., AtaieAshtiani B., Simmons C.T. (2014) Finite volume coupling strategies for the solution of a Biot consolidation model, Comput. Geotech. 55, 494–505. [Google Scholar]
 Nordbotten M.J. (2014) Finite volume hydromechanical simulation in porous media, Water Res. Res. 50, 5, 4379–4394. [CrossRef] [Google Scholar]
 Tang T., Hededal O., Cardiff P. (2015) On finite volume method implementation of poroelastoplasticity soil model, Int. J. Numer. Anal. Meth. Geomech. 39, 13, 1410–1430. [CrossRef] [Google Scholar]
 Nilsen H.M., Nordbotten J., Raynaud X. (2018) Comparison between cellcentered and nodalbased discretization schemes for linear elasticity, Comput. Geosci. 22, 1, 233–260. [Google Scholar]
 Suliman R., Oxtoby O.F., Malan A.G., Kok S. (2014) An enhanced finite volume method to model 2D linear elastic structures, Appl. Math. Model. 38, 7–8, 2265–2279. [Google Scholar]
 Belytschko T., Black T. (1999) Elastic crack growth in finite elements with minimal remeshing, Int. J. Numer. Meth. Eng. 45, 5, 601–620. [CrossRef] [Google Scholar]
 Bordas S., Nguyen P.V., Dunant C., Guidoum A., NguyenDang H. (2007) An extended finite element library, Int. J. Numer. Meth. Eng. 71, 6, 703–732. [CrossRef] [Google Scholar]
 Yan X., Huang Z., Yao J., Li Y., Fan D., Zhang K. (2018) An efficient hydromechanical model for coupled multiporosity and discrete fracture porous media, Comput. Mech. 62, 943–962. [Google Scholar]
 Shi F., Wang X., Liu C., Liu H., Wu H. (2017) An xfembased method with reduction technique for modeling hydraulic fracture propagation in formations containing frictional natural fractures, Eng. Fract. Mech. 173, 64–90. [Google Scholar]
 Ucar E., Keilegavelen E., Berre I., Nordbotten J.M. (2018) A finitevolume discretization for deformation of fractured media, Comput. Geosci. 22, 4, 993–1007. [Google Scholar]
 Deb R., Jenny P. (2014) Modeling of failure along predefined planes in fractured reservoirs, in: Proceedings, ThirtyNinth Workshop Geothermal Reservoir Engineering, Stanford University, Stanford, CA. [Google Scholar]
 Deb R., Jenny P. (2017) Finite volumebased modeling of flowinduced shear failure along fracture manifolds, Int. J. Numer. Anal. Meth. Geomech. 41, 3, 1922–1942. [CrossRef] [Google Scholar]
 Deb R., Jenny P. (2017) Modeling of shear failure in fractured reservoirs with a porous matrix, Comput. Geosci. 21, 5–6, 1119–1134. [Google Scholar]
 Deb R. (2018) Numerical modeling of fluid injection induced shear failure, tensile opening and flowmechanics coupling, PhD Thesis, ETH, Zurich. [Google Scholar]
 Karvounis D.C., Jenny P. (2016) Adaptive hierarchical fracture model for enhanced geothermal systems, Multiscale Model. Simul. 14, 1, 207–231. [Google Scholar]
 Lee S.H., Lough M.F., Jensen C.L. (2001) Hierarchical modeling of flow in naturally fractured formations with multiple length scales, Water Resour. Res. 37, 3, 443–455. [Google Scholar]
 Moinfar A., Sepehrnoori K., Johns R.T., Varavei A. (2013) Coupled geomechanics and flow simulation for an embedded discrete fracture model, Society of Petroleum Engineers. [Google Scholar]
 Gunnar J., Miller S.A. (2017) On the role of thermal stresses during hydraulic stimulation of geothermal reservoirs, Geofluids 2017, 1–15. [Google Scholar]
 Gunnar J., Valley B., Miller S.A. (2018) THERMAIDA matlab package for thermohydraulic modeling and fracture stability analysis in fractured reservoirs, arXiv:1806.10942 [physics.compph] [Google Scholar]
 Shah S., Møyner O., Tene M., Lie K.A., Hajibeygi H. (2016) The multiscale restriction smoothed basis method for fractured porous media (FMsRSB), J. Comput. Phys. 318, 36–57. [Google Scholar]
 Peaceman D.W. (1983) Interpretation of wellblock pressures in numerical reservoir simulation with nonsquare gridblocks and anisotropic permeability, Soc. Pet. Eng. J. 8, 3, 183–194. [Google Scholar]
All Tables
All Figures
Fig. 1 The schematic of the Embedded Discrete Fracture Model (EDFM). The red points represent the matrix grids, which storage fluid pressure; the blue points denote the fracture grids, which storage fracture pressure and displacement for the XFVM; the bright green points storage the matrix displacement for the XFVM and XFEM; e, w, n, s represent the interface of control volume; E, W, N, S represent the internal interface in a control volume. 

In the text 
Fig. 2 A twodimensional elastic body with an internal fracture. 

In the text 
Fig. 3 Computational domain and boundary conditions. The red lines represent the discrete fractures. The red point at the left corner denotes the injection well, and the blue points at the right corner represent the projection well. (a) Single fracture; (b) cross fractures; (c) eight discrete fractures. 

In the text 
Fig. 4 Computational grids used in the COMSOL. (a) 11 400 elements; (b) 16 130 elements; (c) 38 914 elements. 

In the text 
Fig. 5 Computational grids used in XFEM and XFVM. (a) Matrix: 10 201, fracture: 43; (b) matrix: 10 201, fracture; 84; (c) matrix: 10 201, fracture: 232. 

In the text 
Fig. 6 Results comparison of fluid pressure between COMSOL and FVM. 

In the text 
Fig. 7 Comparisons of displacement fields among COMSOL, XFEM, and XFVM for single fracture. The first row represents the displacement fields in the xdirection, the second row represents the displacement fields in the ydirection. Left: results of COMSOL; middle: results of XFEM; right: results of XFVM. 

In the text 
Fig. 8 Comparison of the ydirection displacement field on the fracture upper surface among COMSOL, XFEM, and XFVM. 

In the text 
Fig. 9 Results comparison of fluid pressure between COMSOL and FVM. 

In the text 
Fig. 10 Comparisons of displacement fields among COMSOL, XFEM, and FVM for cross fractures. The first row represents the displacement fields in the xdirection, the second row represents the displacement fields in the ydirection. Left: results of COMSOL; middle: results of XFEM; right: results of XFVM. 

In the text 
Fig. 11 Comparisons of xdirection displacement field on the upper surface of horizontal fracture among COMSOL, XFEM, and XFVM. 

In the text 
Fig. 12 Results comparison of fluid pressure between COMSOL and FVM. 

In the text 
Fig. 13 Comparison of displacement fields among COMSOL, XFEM, and FVM for complex fractures. The first row represents the displacement fields in the xdirection; the second row represents the displacement fields in the ydirection. Left: results of COMSOL; middle: results of XFEM; right: results of XFVM. 

In the text 
Fig. 14 Comparisons of ydirection displacement field on the lower surface of horizontal fracture at the bottom of the domain among COMSOL, XFEM, and XFVM. 

In the text 
Fig. 15 Sparsity matrix pattern for the three cases. Left: XFEM; right: XFVM. 

In the text 