A Fully Coupled Hydro-Mechanical Approach for Multi-Fracture Propagation Simulations

Abstract

Hydraulic fracturing is a complex nonlinear hydro-mechanical coupled process. Accurate numerical simulation is of great significance for reducing fracturing costs and improving reservoir development benefits. The aim of this paper is to propose an efficient numerical simulation method for the fracturing-to-production problem under a unified framework that has good convergence and accuracy. A hydro-mechanical coupled fracturing model (HMFM) is established for poroelastic media saturated with a compressible fluid, and the local characteristics of the physical field are fully considered. Each fracture is explicitly characterized using the discrete fracture model (DFM), which can better reflect the physical characteristics near fractures. Based on the extended finite element method (XFEM) and the Newton–Raphson method, a fully coupled approach named Unified Extended Finite Element (UXFEM) is developed, which can solve the nonlinear system of equations that describe the solution under a unified framework. UXFEM can accurately capture the local physical characteristics of different physical fields on the orthogonal structured grids. It realizes the grid-fracture decoupling, and fractures can propagate in any direction, which shows greater flexibility in simulating fracture propagation. The fully coupled approach can better reflect the essential relationship between pressure, stress, and fracture, which is beneficial to studying hydro-mechanical coupled problems. To validate the UXFEM, UXFEM is compared with the classical KGD model, analytic solution, and COMSOL solution. Finally, based on UXFEM, the interference phenomenon and fracturing-to-production study are carried out to prove the broad practical application prospect of this new fully coupled approach.

1. Introduction

Fracturing, as an effective technical means, can significantly improve the development efficiency of low-permeability and ultra-low-permeability reservoirs, such as shale and tight sandstone. The reservoir often develops various discontinuous structures such as natural fractures, faults, and caves, resulting in significant multi-scale characteristics of the reservoir [1,2,3]. Moreover, all physical processes in the subsurface are highly nonlinear and complex multiphysics problems [4,5], which poses considerable challenges to studying fracture propagation and other related research.

As the most traditional research method, experiments have played a considerable role in promoting the study of fracture propagation [2,3,6,7,8,9,10]. However, the limitations of physical simulation are also prominent, especially for complex multiphysics coupling problems. The experimental cost is relatively high, and the results are often obtained at the laboratory scale, which cannot accurately reflect the underground situation. Furthermore, some complex problems are often unable to be experimentally studied. In this context, numerical simulation technology is applied to the design before fracturing, monitoring during fracturing, and post-fracturing dynamic production studies. Many scholars have conducted extensive and meaningful research on the fracture propagation problem in petroleum engineering and proposed some numerical simulation methods. Those methods mainly include the finite element method (FEM) and its derivative methods, the boundary element method (BEM), and the discrete element method (DEM).

FEM is a flexible, effective, and widely used numerical method [11,12,13,14,15,16,17]. Fractures are highly coupled with grids, and fractures must be set along the boundary of mesh grids. Therefore, fractures must propagate along the grid boundary or continuously reconstruct mesh grids, which requires tremendous calculation. At the same time, the shape function of FEM is continuous. It is often difficult to accurately describe various discontinuities when describing the fluid pressure and the solid displacement fields. Thus, it is often necessary to use local mesh refinement to achieve high-accuracy calculations, which also increases the calculation burden. The characteristics of fracture geometry and physical field are challenging to characterize efficiently and accurately in traditional FEM, and the limitations of FEM are significant.

Some FEM-based extension methods have been proposed successively, such as node-splitting FEM [18,19], generalized FEM [20,21], and XFEM [22,23,24,25,26,27]. Node-splitting FEM allows FEM to describe the fracture width simply. Still, the fracture propagation path must follow the grid boundary, which cannot accurately describe the fracture propagation and interaction. XFEM is a method based on the partition of unity method (PUM), which uses enrichment functions to capture the physical characteristics of fracture walls and fracture tips. XFEM is an efficient method for solving discontinuous problems, which can well solve strong and weak discontinuity problems. The core idea is to capture various discontinuities with the help of enrichment functions constructed based on analytical solutions or asymptotic analytical solutions, and fractures and grids are independent of each other. XFEM, based on the discrete fracture model, can explicitly characterize each fracture. Many scholars have constructed enrichment functions to describe fracture surfaces, fracture tips, and intersecting fractures based on physical field characteristics and analytical solution characteristics, which can efficiently deal with complex fracture mechanical behavior. However, most of the current XFEM models cannot accurately consider the fluid flow and exchange between fractures and matrix, usually ignoring the fluid leak-off term during fracturing stimulation or replacing it with an empirical leak-off equation, which cannot meet the simulation requirements. Moreover, XFEM also has considerable troubles in the numerical simulation of complex fracture flow, and it is difficult to achieve a fully coupled approach under a unified numerical framework.

BEM [28,29,30,31,32,33,34] is a dimensionality reduction research method based on Green’s formula. It can calculate fracture aperture and stress more accurately and conveniently than FEM and is suitable for simulating fracture propagation with complex topologies. However, this method faces two main problems. One is that it cannot accurately describe the influence of fractures on fluid/solid physical fields. The other one is that it is difficult to consider the leak-off during fracture propagation [35,36].

DEM [37] is a numerical simulation method specially used to solve the problem of discontinuous media. It can deal with fracture intersecting, branching, merging, and kinking problems. However, calibrating/updating particle properties remains a complex technical challenge. The particle number needs to be large enough to achieve sufficient accuracy, which prevents this technique from being widely used in large-scale models [36].

Different methods have different advantages and different defects. In order to synthesize their advantages and compensate for each other’s weaknesses, some hybrid numerical simulation methods are also widely used. Recently, Li et al. [36,38,39,40] developed a meaningful thermal/hydro-mechanical (THM) model. They solved it using the hybrid numerical simulation method, realizing the two-dimensional and three-dimensional simulation of the construction of a complex fracture network. Settgast et al. [41] developed a fully coupled finite element/finite volume approach for simulating field-scale hydraulically driven fractures in three dimensions. Guo et al. [42] adopted the mixed finite element/displacement discontinuity method to solve for the spatial–temporal evolutions of pore pressure and in situ stress because of parent-well production and injection and models the fracture propagation during infill-well completion based on updated heterogeneous in situ stresses. Recently, a coupled simulation strategy combining the embedded discrete fracture method (EDFM) and the XFEM has been developed to simulate the fluid-driven fracture propagation process in porous media. EDFM and XFEM are used to simulate fracture-related fluid mechanics and solid mechanics, respectively, with information exchanged under the iterative numerical coupling scheme, and it realizes two-dimensional and three-dimensional hydraulic fracturing fracture propagation simulation [43,44,45]. Liu et al. [46] developed a hydro-mechanical model for non-planar hydraulic fracture propagation in ductile formations, which is solved by the hybrid extended finite element/finite volume method. Zhang et al. [47] combined the extended finite element/phase-field method to solve the discontinuous and continuous hydraulic fracturing formulations.

However, there are still many insurmountable problems in developing realistic simulation tools for the hydraulic fracturing process. The problems of numerical simulation of fracture propagation are mainly reflected in the contradiction between calculation efficiency and calculation accuracy, reflected explicitly in several aspects: (1) There are strict requirements for grid division; (2) The fracture propagation path is not arbitrary, and the fracture morphology of simulation results is distorted; (3) Matrix flow, fracture flow, and solid deformation coupling are challenging, and the leak-off term during the fracturing process is often ignored or simplified; (4) The nonlinearity is strong, and the fully coupled approach is scarce; (5) There are still risks and challenges in the convergence and stability of numerical calculation. The reasons for this contradiction are complex and diverse, such as rock heterogeneity, unclear mechanism of fracture propagation, and aggravation of heterogeneity caused by the intersection of fractures and artificial stimulation measures. In recent years, some new methods and models have been proposed and achieved good results, such as Peridynamics [48,49,50,51].

This paper uses the discrete fracture model (DFM) to characterize each fracture [52,53,54,55,56]. A fully coupled numerical approach for solving nonlinear hydro-mechanical coupled models is established based on XFEM. In this method, the solid rock deformation obeys quasi-static linear elasticity, characterized by the stress balance equation. The matrix flow follows Darcy’s law, and the fracture fluid flow obeys the cubic law. In the derivation, this paper uses the normal flow velocity discontinuity term in the weak form of matrix and fracture flow equations to automatically characterize the fluid leak-off term. This paper uses XFEM to compile a fully coupled UXFEM solver for the abovementioned HMFM. This approach solves the fluid pressure field and the solid displacement field under the same framework, realizing the unification and full coupling of the model and solving.

According to the different physical characteristics of solid deformation and fluid flow, UXFEM adopts different enrichment functions and successfully realizes the numerical simulation of hydro-mechanical fracture propagation. The advantages of XFEM are entirely inherited in UXFEM. Fractures are decoupled from mesh grids, and fractures are allowed to propagate in any direction and can be deflected at any angle without mesh reconstruction or other special processing. In this paper, the maximum circumferential tensile stress criterion is used as the fracture propagation criterion, the displacement-based method is used to calculate the stress intensity factors (SIFs), and fractures always propagate along the direction of the maximum circumferential stress. When UXFEM solves nonlinear mathematical models, it adopts Newton–Raphson iterative linearization to deal with nonlinear systems. With the advantages of XFEM, UXFEM can achieve the balance between calculation efficiency and accuracy to a certain extent, which has practical significance for the actual engineering-scale fracturing simulation research.

2. Unified Extended Finite Element Method

2.1. Fracture Description

This paper uses pair of orthogonal level set functions $f(x),\mathrm{g}(x)$ to describe each fracture explicitly, and some enrichment functions are constructed based on level set functions. The fracture surface is defined as $f(x)=0$, and the fracture tip is defined as $f(x)=0,g(x)=0$.

The core of XFEM is to capture the discontinuity of the physical field by using enrichment functions with discontinuous properties and enriched degree of freedoms (DOF) on the basis of FEM, which makes the description of the physical field characteristics independent of grids, bringing a lot of conveniences. The fracture is a strong discontinuity for the displacement field, while the fracture is a weak discontinuity for the fluid pressure field. That means the pressure on the two sides of the fracture is continuous, but the pressure derivative is discontinuous. In this paper, the set of all nodes is denoted as $N_{all}$, and the enriched nodes, including surface nodes and tip nodes, are denoted as $N_s$ and $N_t$, respectively (Figure 1).

2.2. Solid Displacement Field Approximation and Enrichment Functions

In the XFEM framework, the displacement approximation is written as:

$$u(\mathit{x})={\displaystyle \sum _{i\in N_{all}}{N}_i(\mathit{x})u_i}+{\displaystyle \sum _{j\in N_s}{N}_j(\mathit{x})(H(f(\mathit{x}))-H((f({\mathit{x}}_j)))}{a}_j+{\displaystyle \sum _{q\in N_t}{N}_q(\mathit{x}){\displaystyle \sum _{l=1}^4(F_l(\mathit{x})-F_l({\mathit{x}}_q))}}{b}_q^l$$

(1)

where $N_i(\mathit{x}),N_j(\mathit{x}),N_q(\mathit{x})$ are standard FEM shape function; $u_i$ is standard displacement degrees of freedom (DDOF) for nodes $N_{all}$; $a_j^n,b_{ql}^m$ are added enriched DDOFs for nodes $N_s$ and $N_t$ respectively; $H(\cdot )$ is the Heaviside step function; $F_l(\cdot )$ is fracture tip displacement enrichment function.

For expression brevity, we combine the last two terms in Equation (1) and rewrite Equation (1) as:

$$u(\mathit{x})={\displaystyle \sum _{i\in N_{all}}{N}_i(\mathit{x})u_i}+{\displaystyle \sum _{j\in N_{enr}}{\Psi }_j(\mathit{x}){\tilde{u}}_j}=\mathit{N}\cdot \mathit{U}$$

(2)

where $N_{enr}=\left\{N_s,N_t\times 4\right\}$; ${\Psi }_j$ and ${\tilde{u}}_j$ denote enrichment functions and its added enriched DDOFs; $\mathit{N}$ is the combination of standard FEM shape function $N_i$ and enrichment function ${\Psi }_j$; and $\mathit{U}$ is the combination of standard DDOFs and added enriched DDOFs.

For the element completely penetrated by the fracture, the displacement field is discontinuous on both sides, which is a strong discontinuity. The Heaviside step function is used to enrich the surface nodes.

$$H(f(\mathit{x}))=\left\{\begin{array}{l}1,f(\mathit{x})\ge 0\hfill \\ 0,f(\mathit{x})<0\hfill \end{array}\right.$$

(3)

For fracture tips, the stress has the singularity of $\mathrm{O}(r^{-\frac{1}{2}})$. For isotropic linear elastic materials, tip nodes can be enriched by the fracture tip displacement enrichment function [22,57,58]:

$$[F_l(r,\theta ),l=1,2,3,4]=[\sqrt{r}\sin \frac{\theta }{2}\sqrt{r}\cos \frac{\theta }{2}\sqrt{r}\sin \frac{\theta }{2}\sin \theta \sqrt{r}\cos \frac{\theta }{2}\sin \theta ]$$

(4)

where $r,\theta$ are polar coordinates with origin at the fracture tip and $x_1$ axis oriented into the body and parallel to the fracture surfaces (Figure 2).

2.3. Fluid Pressure Field Approximation and Enrichment Functions

In the XFEM framework, the pressure approximation is written as:

$$p(\mathit{x})={\displaystyle \sum _{i\in N_{all}}{N}_i(\mathit{x})p_i}+{\displaystyle \sum _{a\in N_s}{N}_a(\mathit{x})({\varphi }_s(x)-{\varphi }_s(x_a))}{\tilde{p}}_a+{\displaystyle \sum _{b\in N_t}{N}_b(\mathit{x})({\varphi }_t(x)-{\varphi }_t(x_b))}{\tilde{p}}_b^m$$

(5)

where $N_i(\mathit{x}),N_a(\mathit{x}),N_b(\mathit{x})$ are standard FEM shape function; $p_i$ is standard pressure degrees of freedom (PDOFs) for node $N_{all}$; $a_j^n,b_{ql}^m$ are added enriched PDOFs for nodes $N_s$ and $N_t$ respectively; ${\varphi }_s(\cdot )$ is the modified level set absolute value function; ${\varphi }_t(\cdot )$ is fracture tip pressure enrichment function.

For expression brevity, we combine the last two terms in Equation (5) and rewrite Equation (5) as:

$$p(\mathit{x})={\displaystyle \sum _{i\in N_{all}}{N}_i(\mathit{x})p_i}+{\displaystyle \sum _{j\in N_{enr}}{\Phi }_j(\mathit{x}){\tilde{p}}_j}=\mathit{H}\cdot \mathit{P}$$

(6)

where $N_{enr}=\left\{N_s,N_t\right\}$; ${\Psi }_j$ and ${\tilde{u}}_j$ denote enrichment functions and their added enriched PDOFs; $\mathit{H}$ is the combination of standard FEM shape function $N_i$ and enrichment function ${\mathrm{\Phi }}_j$; and $\mathit{P}$ is the combination of standard PDOFs and added enriched PDOFs.

For the fracture surface, the pressure is continuous, but the pressure gradient is discontinuous. Moës [59] uses the modified level set absolute value function to capture the fluid pressure characteristics, avoiding the appearance of blending elements.

$${\varphi }_s(\mathit{x})={\displaystyle \sum _j\left|f({\mathit{x}}_j)\right|N_j(\mathit{x})-\left|{\displaystyle \sum _{j}f({\mathit{x}}_j)N_j(\mathit{x})}\right|}$$

(7)

According to the study by Chen et al. [60], the pressure gradient has the singularity of $\mathrm{O}(r^{-\frac{1}{2}})$. A new enrichment function was constructed based on the asymptotic analytic solution [61,62,63].

$${\varphi }_t(\mathit{x})=\sqrt{r}\cos \frac{\theta }{2}$$

(8)

2.4. XFEM Discretization

This paper only considers the quasi-static process with infinitesimal strain in porous media, which is homogeneous, isotropic, and linearly elastic. The porous medium domain $\mathrm{\Omega }$, with boundary $\mathrm{\Gamma },{\mathrm{\Gamma }}_t,{\mathrm{\Gamma }}_u$ (Figure 3). The fracture ${\mathrm{\Gamma }}_f$ consisting of two surfaces is embedded in the domain $\mathrm{\Omega }$. In this paper, fluid flow and matrix deformation are coupled based on poroelasticity [64], matrix and fracture flow are coupled based on discontinuous flow on the fracture surface, and fracture deformation and fracture flow are coupled using the cubic law. The coupled governing equations are Equation (A4) and (A21) in Appendix A. For the detailed formula derivation, please refer to the Appendix A.

According to the basic theory of XFEM, the fracture opening can be expressed as:

$$w_f=〚\mathit{u}〛\cdot {\mathit{n}}_f=\left({\displaystyle \sum _{i\in N_s}{N}_i(x)a_i}+2{\displaystyle \sum _{q\in N_t}\sqrt{r}}{N}_q(x)b_q\right)\cdot {\mathit{n}}_f=〚N〛\cdot \mathit{U}\cdot {\mathit{n}}_f$$

(9)

where $〚\mathit{u}〛$ denotes the displacement jump across the fracture; $a_i\mathrm{and}{b}_q$ are discontinuous enriched DDOFs; $〚N〛$ denotes the discontinuous shape function matrix.

By substituting the displacement and pressure approximation expressions into the weak form of the governing equations, the discrete calculation format can be obtained as follows:

$$\left[\begin{array}{cccc}0& 0& 0& 0\\ 0& 0& 0& 0\\ M_{uu}& M_{u\tilde{u}}& M_{pp}& M_{p\tilde{p}}\\ M_{\tilde{u}u}& M_{\tilde{u}\tilde{u}}& M_{\tilde{p}p}& M_{\tilde{p}\tilde{p}}\end{array}\right]\left[\begin{array}{c}\dot{U}\\ \dot{\tilde{U}}\\ \dot{P}\\ \dot{\tilde{P}}\end{array}\right]+\left[\begin{array}{cccc}{K}_{uu}& K_{u\tilde{u}}& K_{pp}& K_{p\tilde{p}}\\ K_{\tilde{u}u}& K_{\tilde{u}\tilde{u}}& K_{\tilde{p}p}& K_{\tilde{p}\tilde{p}}\\ 0& 0& K{K}_{pp}& K{K}_{pp}\\ 0& 0& K{K}_{pp}& K{K}_{pp}\end{array}\right]\left[\begin{array}{c}U\\ \tilde{U}\\ P\\ \tilde{P}\end{array}\right]=\left[\begin{array}{c}{F}_1\\ {\tilde{F}}_1\\ F_2\\ {\tilde{F}}_2\end{array}\right]$$

(10)

Equation (10) can be simplified to the following form:

$$\left[\begin{array}{cc}0& 0\\ {\mathit{M}}_u& {\mathit{M}}_p\end{array}\right]\left[\begin{array}{c}\dot{\mathit{U}}\\ \dot{\mathit{P}}\end{array}\right]+\left[\begin{array}{cc}{\mathit{K}}_u& {\mathit{K}}_p\\ 0& \mathit{K}{\mathit{K}}_p\end{array}\right]\left[\begin{array}{c}\mathit{U}\\ \mathit{P}\end{array}\right]=\left[\begin{array}{c}{\mathit{F}}_1\\ {\mathit{F}}_2\end{array}\right]$$

(11)

where

$${\mathit{K}}_u={\displaystyle {\int }_{\mathrm{\Omega }}{\mathit{B}}^T:\mathit{C}:\mathit{B}d\mathrm{\Omega }}$$

$${\mathit{K}}_p=-{\displaystyle {\int }_{\mathrm{\Omega }}{\mathit{B}}^T\alpha {\varsigma }^T{\mathit{H}}^{T}d\mathrm{\Omega }-{\displaystyle {\int }_{{\mathrm{\Gamma }}_f}〚N〛{\mathit{n}}_f{\mathit{H}}^{T}d\mathrm{\Gamma }}}$$

$${\mathit{F}}_1={\displaystyle {\int }_{{\mathrm{\Gamma }}_t}\mathit{N}\cdot \overline{\mathit{\sigma }}d\mathrm{\Gamma }}-{\displaystyle {\int }_{\mathrm{\Omega }}\mathit{N}\cdot ({\sigma }_0+\alpha p_0\mathit{I})d\mathrm{\Omega }}$$

$${\mathit{M}}_u={\displaystyle {\int }_{\mathrm{\Omega }}\alpha \rho \mathit{H}\varsigma \mathit{B}d\mathrm{\Omega }+{\displaystyle {\int }_{{\mathrm{\Gamma }}_f}\rho \mathit{H}{\mathit{n}}_f{〚N〛}^{T}d\mathrm{\Gamma }}}$$

$${\mathit{M}}_p={\displaystyle {\int }_{\mathrm{\Omega }}\rho Q^{-1}\mathit{H}{\mathit{H}}^T}d\mathrm{\Omega }+{\displaystyle {\int }_{{\mathrm{\Gamma }}_f}\frac{\rho w}{K_l}\mathit{H}{\mathit{H}}^T}d\mathrm{\Gamma }$$

$$\mathit{K}{\mathit{K}}_p={\displaystyle {\int }_{\mathrm{\Omega }}\frac{\rho k_m}{\mu }}{\mathit{D}}^T\mathit{D}d\mathrm{\Omega }+{\displaystyle \underset{{\mathrm{\Gamma }}_f}{\int }\frac{\rho w^3}{12\mu }{\mathit{D}}_s{\mathit{D}}_s{}^{T}d\mathrm{\Gamma }}$$

$${\mathit{F}}_2={\displaystyle {\int }_{{\mathrm{\Gamma }}_f}H{Q}_{m}d\mathrm{\Gamma }}-{\displaystyle {\int }_{{\mathrm{\Gamma }}_q}\mathit{H}\overline{q}d\mathrm{\Gamma }}$$

$$\varsigma =[1,1,0]$$

$$\mathit{N}={[{\mathit{N}}^u{\mathit{N}}^{\tilde{u}}]}^T,\mathit{B}=[\nabla {\mathit{N}}^u\nabla {\mathit{N}}^{\tilde{u}}]$$

$$\mathit{H}={[{\mathit{N}}^p{\mathit{N}}^{\tilde{p}}]}^T,\mathit{D}=[\nabla {\mathit{N}}^p\nabla {\mathit{N}}^{\tilde{p}}]$$

where ${\mathit{H}}_s$ is the total shape function at the point source, ${\mathit{D}}_s$ is the directional derivative of ${\mathit{H}}_s$ along the tangential direction of the fracture, i.e., ${\mathit{D}}_s=\partial {\mathit{H}}_s/\partial s$. The superscripts “$u$” and “$\tilde{u}$” correspond to standard DDOFs and enriched DDOFs, respectively, and the superscripts “$p$” and “$\tilde{p}$” correspond to standard PDOFs and enriched PDOFs, respectively. Equation (11) is a set of coupled nonlinear equations, which are solved by the Newton–Rapson method.

3. Fracture Propagation and Solution Strategy

3.1. Fracture Propagation Criterion

In practice, the hydraulic fracture usually propagates in a mixed mode. The maximum circumferential stress criterion is adopted as the fracture propagation criterion. It is assumed that the fracture propagates when the effective stress intensity factor $K_e$ along that direction reaches the fracture toughness $K_{IC}$, and it will deflect by an angle of [65]:

$$\theta =2\arctan \left(\frac{-2K_{II}/K_I}{1+\sqrt{1+8{(K_{II}/K_I)}^2}}\right)$$

(12)

The effective stress intensity factor $K_e$ is expressed as:

$$K_e=\cos \frac{\theta }{2}\left(K_I{\cos }^2\frac{\theta }{2}-1.5K_{II}\sin \theta \right)\ge K_{IC}$$

(13)

The displacement-based approach is used to calculate the SIFs. $K_I$ and $K_{II}$ are determined by the discontinuous displacement of the fracture tip element. The relational equation is [66].

$$\left\{\begin{array}{c}{K}_I=\frac{0.806\sqrt{\pi }E{D}_n}{4(1-{\nu }^2)\sqrt{r}}\\ K_{II}=\frac{0.806\sqrt{\pi }E{D}_s}{4(1-{\nu }^2)\sqrt{r}}\end{array}\right.$$

(14)

where $E$ is Young’s modulus, $\nu$ is Poisson’s ratio, $D_n$ and $D_s$ are normal and tangential discontinuous displacements of fractures. It should be noted that the calculation models in this paper are all under the plane strain condition.

3.2. HMFM Fracturing Simulation Process

Fracture propagation is a highly nonlinear hydro-mechanical coupled process. In this paper, the fully coupled approach is implemented. The fully coupled approach forms a single large system of equations that solve for all of the displacement/pressure unknowns at once. Figure 4 summarizes the fully coupled implementation of the solution of the HMFM.

4. Results and Discussions

4.1. Validation 1: Comparison between UXFEM and KGD Model

The KGD model is one of the most commonly used fracturing models. The viscosity-dominated hydraulic pressure is calculated using the approximate solution of the KGD model derived by Detournay [67].

Table 1. Parameters for validation 1.

Parameters	Value	Unit
Young modulus	60	GPa
Poisson ratio	0.25	/
Fracture toughness	1.2	MPa·m0.5
Fluid viscosity	50	mPa·s
Injection rate	0.01	m2/s
Convergence tolerance	0.01	/

lists the parameters used by the model. The fracture length and width at the wellbore position are compared with those calculated by UXFEM and plane strain KGD model, as shown in Figure 5. It can be seen that the evolution of fracture length and width calculated by UXFEM in this paper can well fit the KGD model. The calculation error increases with the fracture length because the infinite reservoir assumption cannot be satisfied with the fracture propagation.

4.2. Validation 2: Comparison of KI and KII in Numerical Simulation and Analytical Solution

Rice [68] deduces the analytical expressions of the SIFs at the fracture tip under hydrostatic pressure and horizontal stress. For the single fracture propagation in an infinite domain, the SIFs are expressed as:

$$\left\{\begin{array}{l}{K}_I=\sqrt{\pi L_f}\left[p-\left({\sigma }_H{\sin }^2\beta +{\sigma }_h{\cos }^2\beta \right)\right]\hfill \\ K_{\mathrm{II}}=\frac{\sqrt{\pi L_f}}{2}\left[{\sigma }_H-{\sigma }_h\right]\sin 2\beta \hfill \end{array}\right.$$

(15)

where $L_f$ is the half length of the fracture; $p$ is hydrostatic pressure in the fracture; S is the inclination angle of fracture and maximum horizontal principal stress. ${\sigma }_H$ and ${\sigma }_h$ are the maximum horizontal stress and the minimum horizontal stress, respectively.

In order to approximate the infinite formation and eliminate the influence of the boundary, the model size is set to 200 × 200 m, the initial half-length of the fracture is 5 m, and the domain is discretized into structured grids (159 × 159).

Table 2. Parameters for validation 2.

Parameters	Value	Unit
Young’s modulus	60	GPa
Poisson’s ratio	0.25	/
Maximum horizontal in-situ stress	60	MPa
Minimum horizontal in-situ stress	55	MPa
Hydrostatic pressure applied to the fracture surface	65	MPa
Half-length of the fracture	5	m

lists the parameters used by the model. Calculate the SIFs (K_I and K_II) of different inclination angles, and the results are compared, as shown in Figure 6. By comparison, the SIFs calculated by UXFEM in this paper can well fit the analytical solution. The maximum error of KI is not more than 2%, which proves the method’s accuracy.

4.3. Validation 3: Comparison between UXFEM and COMSOL Multiphysics^® 5.6

The purpose of this section is to show the accuracy and feasibility of our fully coupled approach. We consider a quasi-static mud loss case of a shale oil reservoir. The simulation parameters are shown in

Table 3. Parameters for Validation 3.

Parameters	Value	Unit
Initial reservoir pressure	10	MPa
Pressure at injection point	25	MPa
Maximum horizontal in-situ stress	20	MPa
Minimum horizontal in-situ stress	15	MPa
Young’s modulus	40	GPa
Poisson’s ratio	0.2	/
Matrix permeability	1 × 10−18	m2
Initial matrix porosity	0.15	/
Bulk modulus of the solid	50	GPa
Bulk modulus of the fluid	2.5	GPa
Fluid density	1000	kg/m3
Fluid viscosity	10	mPa·s
Biot coefficient	0.85	/

. The geometry (Figure 7) is built to carry out the simulation process. The domain is discretized into 25,860 unstructured grids in COMSOL and 3481 structured grids in UXFEM. Figure 8 shows the displacement field, effective stress component ${\mathit{\sigma }}_y^{\prime }$, and pressure field calculated by COMSOL and UXFEM. In Figure 8, the results calculated by UXFEM and COMSOL are in good agreement. From Figure 8a,d, the displacement at the injection point is the largest, implying that the fracture width is the largest at the injection point. Figure 8b,e shows that there is an obvious stress concentration at the fracture tip. Figure 8c,f shows that the pressure is higher along the fracture, and a small low-pressure area appears around the tips. Figure 9 shows that fracture widths calculated by UXFEM and COMSOL can be matched. Calculating the L₂ error norm:

$${‖F-F^h‖}_{L_2}={\left[{\displaystyle \underset{\mathrm{\Omega }}{\int }(F-F^h)d\mathrm{\Omega }}\right]}^{0.5}$$

(16)

where $F^h$ is the COMSOL solution and $F$ is the UXFEM solution.

The displacement L₂ error, effective stress component ${\mathit{\sigma }}_y^{\prime }$ L₂ error, and pressure L₂ error are 0.0083%, 0.041%, and 0.078%, respectively. A very good agreement between the results proves the effectiveness and accuracy of this fully coupled approach.

4.4. Case Study1: Propagation of Three Parallel Fractures

A 200 m × 200 m two-dimensional model is established to study the fracture deflection propagation caused by the stress shadow effect. The initial half-length of the fracture was 5 m, the distance between the three fractures was 14 m or 21 m, and the midpoints of the three fractures are the injection points. The injection rate at injection points is 0.01 m²/s. Other parameters are listed in

Table 4. Parameters for case study.

Parameters	Value	Unit
Initial reservoir pressure	50	MPa
Maximum horizontal in situ stress	60	MPa
Minimum horizontal in situ stress	55	MPa
Matrix permeability	1 × 10−17	m2
Initial matrix porosity	0.1	/
Initial fracture length	10	m
Injection rate at injection points	0.01	m2/s
Bulk modulus of the fluid	2.5	GPa
Fluid density	1000	kg/m3
Fluid viscosity	50	mPa·s
Fracture toughness $K_{IC}$	5	MPa·m0.5
Biot coefficient	0.85	/

.

The final morphology of the simulated fracture, displacement field, and effective stress component ${\mathit{\sigma }}_y^{\prime }$ are shown in Figure 10. For the model of three fractures, the middle fracture is squeezed and suppressed by the same on both sides. The middle fracture almost propagates along a straight line, and its length is significantly shorter than the other two fractures. Due to the stress interference between fractures, the side fractures deflect in the opposite direction, away from the middle fracture [69]. Fractures will actively choose to propagate along the path with the least resistance. As the fracture spacing is smaller, the fracture deflection angle is larger, indicating that the fracture spacing is an important factor affecting the stress interference intensity.

4.5. Case Study2: Propagation in Multiple En Échelon Fractures

En échelon fractures are common in geology and are caused by the mechanical interaction between their near-tip stress fields [36,70]. To study the evolution of En échelon fractures, a 200 × 200 m geometry model is built. The model contains three prefabricated fractures, with an initial fracture length of 10 m and a vertical spacing of 16 m. The initial parameters of the three fractures are completely the same. The injection point is located at the midpoint of the fracture (The red points in Figure 11), and the injection rate is set to 0.01 m²/s. Other parameters of the model are the same as those in 5.4 (Table 4). See Figure 11 for the calculation results. It can be seen that after the fractures overlap, the propagation of the middle fracture is restrained, and the tendency of mutual attraction between adjacent fractures is becoming more and more obvious. Enhancement or impediment of the growth ahead of an individual fracture front leads to asymmetric fracture geometry [36].

The stress interference of parallel fractures has been studied previously, and the fractures tend to mutually repel each other. However, for such parallel fractures with incomplete left-right overlap, the stress shadow effect will cause the stress direction to change, and the overlapping tips will appear as a phenomenon of mutual attraction. In the process of hydraulic fracturing, we hope that hydraulic fractures can form complex shapes. However, this phenomenon may lead to fracture collision and reduce the fracturing effect, which needs to be avoided. Therefore, the optimization of the overlapping length and spacing of fractures is an important job.

4.6. Case Study3: Multi-Well and Multi-Cluster Simultaneous Fracturing and Production

The problem considered here is the fracturing-production problem of 7 hydraulic fractures generated by two horizontal wells (#1 and #3) and one vertical well (#2), numbered 1 to 7, respectively. The schematic diagram of the initial fractures and the model is shown in Figure 12. The size of the model is 200 m × 100 m. Two horizontal wells are along the straight-line x = 60 m and x = 140 m, respectively. The vertical well is located in the middle of the model. Table 4 lists the detailed parameters.

This study aims to analyze the interaction between evolved multi-fractures and local stress-field-pressure fields and help explore the interference phenomenon’s influence on fracture morphology and production. Taking well #1 as an example, it can be seen that the fractures are subject to obvious stress interference from Figure 13 and Figure 14, and the stress interference is manifested in two points:

(1)

The fractures will undergo obvious diversion. For well #1, Frac. 1–3 have different degrees of deflection. The deflection of the two Frac. 1 and 3 are more prominent. In the early stage, the middle fracture Frac. 2 does not deflect obviously due to the interference of superimposed stress. In the later stage of fracturing, Frac. 2 and 4 attract each other, showing En échelon fracture morphology [69].

(2)

There are noticeable differences in fracture length and width. The length of Frac. 2 and 6 of the two horizontal wells is shorter, followed by Frac. 3 and 5, and the longest is Frac. 1 and 7. Frac. 4 cannot propagate far enough and even the front edge of the fracture tends to close. Although Frac. 1 and 7 have propagated long enough, one side of the fracture tends to close under the effect of stress shadow. These phenomena are very likely to lead to the reduction of fracturing efficiency and cost waste.

The complex morphology of fractures is the result of stress interference. Studies have found that 25–30% of the fractures do not bring production [71,72], indicating that 25–30% of the fractures are not opened or closed after being opened, which is also the effect of stress interference. Therefore, stress interference has both positive and negative effects. Numerical simulation can be used to find the best design for well placement, perforation, and fracturing so that the fractures can be fully opened and stabilized to maximize the benefits of fracturing.

Change the bottom hole condition to a constant pressure of 25 MP, and proppant instead of the hydraulic pressure supports the fracture. Assuming that the fracture has high conductivity and obeys Darcy’s law, it is only necessary to change the fracture flow coefficient $\frac{w^2}{12\mu }$ in Equation (A4) into the hydraulic fracture permeability. The stress-sensitive parameters are brought into the governing equation for calculation during production. For the estimation method of stress-sensitive parameters, see the study of Lu et al. [73]. The pressure change in the production process is shown in Figure 15, which shows that there is almost no production interference between fractures and wells within 30 days, and obvious interference can be observed after 300 days. The calculated cumulative production is shown in Figure 16, which clearly shows the rising pattern of the cumulative production. The early production increases rapidly because the reservoir near the well is fully stimulated. However, the reservoir’s extremely low permeability makes it difficult to maintain rapid production.

Four calculation examples with different well-spacing (40 m, 60 m, 80 m, and 100 m) are simulated for comparison to verify the effect of well-spacing on fracture morphology and production. After 60 s injection, the fracture morphology and displacement field are as Figure 17. The final morphology of fractures under the stress interference is different with different well-spacing. Then, realize conversion from hydraulic fracturing to production, as shown in Figure 18, and the cumulative production curve after 2000 days is plotted in Figure 19. Figure 18 shows that small well-spacing will lead to slow pressure drop near boundaries, resulting in a low production rate. Thus, the cumulative production of small well-spacing is smaller than that of large well-spacing, which can be explained by the small well-spacing leading to a poor fracturing effect. Therefore, reasonable well-spacing is a crucial factor affecting the formation of a complex fracture network and efficient production.

5. Conclusions

Aiming at solving the fracturing-to-production problem under the unified framework, a fully coupled approach UXFEM for solving the HMFM under the same framework is established based on the extended finite element method and Newton–Rapson method. This technology inherits all the excellent characteristics of XFEM and can fully reflect the mechanical mechanism of fracture propagation. Meanwhile, it also avoids some potential problems caused by the iterative coupling method, such as convergence reduction. The effectiveness and accuracy of UXFEM in HMFM research are validated by comparing it with the results of classical KGD models, analytical solutions, and COMSOL Multiphysics^® 5.6. To study the interference problem during the fracturing process, three examples are used to illustrate the effect of stress shadow on fracturing. At the same time, after simultaneous fracturing, this paper also successfully realized the fracturing-to-production analysis. The research shows that hydraulic fractures are interfered with by the stress between wells and fractures, which plays an essential role in forming the final fracture network morphology and cumulative production. Proper well location design, perforation plan, and production plan are the keys to determining whether the stress interference can be used reasonably to increase the fracturing effects and reservoir development efficiency. This research has certain practical significance for fracturing design and production prediction.