Simulation of Fracture Morphology during Sequential Fracturing

Abstract

During hydraulic fracturing, the aperture of hydraulic fractures will shrink by the in-situ stress, but will not fully close because of the existence of proppant inside the fracture. In previous studies, few people noticed the existence of proppant, which has resulted in the inaccuracy of simulation results. In this study, based on the boundary element method, a numerical simulation model for sequential fracturing was established, which respectively considered the influence of proppant in staged fracturing and zipper fracturing. In addition, the influence mechanism of proppant on fracture morphology is then revealed. Simulation results show that the residual aperture of the previous hydraulic fracture, which was produced by proppant, may increase with the increase of proppant stiffness and fracture spacing and may also be shrunk by the dynamic propagation of subsequent hydraulic fracture. However, the residual aperture will rebound after hydraulic fracturing construction is finished. The shrinkage and rebound values of residual aperture of hydraulic fracture are usually less than 1 mm. In addition, at the same time, the residual aperture of previous hydraulic fracture may also influence the propagation of subsequent hydraulic fracture. These influences are represented by the bend of fractures in multistage fracturing and the intersection in zipper fracturing. With the increase of well spacing, the influence degree of residual aperture on subsequent fracture propagation is reduced. The previous hydraulic fracture cannot have a significant effect on the deflection of subsequent hydraulic fracture when fracture spacing is between 10 and 30 m. The above research has important guiding significance for controlling fracture morphology in hydraulic fracturing.

1. Introduction

Unconventional oil and gas resources are important components of traditional fossil fuels [1]. Studies have shown that human consumption of oil and natural gas remains dominant. However, due to the low permeability and porosity of the unconventional reservoir, it cannot achieve economic productivity without technological measures of increasing recovery efficiency. Waterflooding is one of these methods. However, waterflooding has the problem of salt precipitation and permeability reduction. Salt precipitation may result in permeability reduction during waterflooding under high pressure and high-temperature conditions [2,3] in addition to permeability reduction due to asphaltene precipitation in the near-wellbore region of oil wells [4]. In order to reduce the production costs and improve the recovery efficiency of oil and gas, various fracturing techniques have been proposed, like the horizontal well staged fracturing technology [5], zipper fracturing technology [6], synchronous fracturing technology [7], refracturing technology [8], etc. The common point of the above technologies is that multiple fractures propagate sequentially and interfere with each other. Especially in recent years, with the proposal of close cutting fracturing technology and infill wells conception, the interference situation between fracturing wells and hydraulic fractures (HFs) has received obsessive attention [9,10]. Revealing the disturbance between HFs and fracturing wells can help to improve the effectiveness of reservoir fracturing, ultimately enhancing oil and gas recovery [11]. Therefore, it is necessary to study the effects of multi-fracture interactions on fracture trajectories and deformation.

Due to the limitation of sample size, it is different to display the mutual interference during fracture propagation in physical experiments. Thus, numerical simulation technology is an important way to study this scientific problem. Zheng [12] simulated the well interference on hydraulic fracture propagation during zipper fracturing, proposing that asymmetric fracture propagations occur in multi-well fracturing, and the lateral growths of interior fractures are suppressed due to the intense inter-well stress interference. Guo [13] investigated the effect of six factors on the fracture propagation. These factors include plasticity and brittleness, horizontal stress difference, cluster spacing, fracturing fluid viscosity, pumping rate, and cracking sequence. The results showed that, in simultaneous fracturing, the fracture network is formed spontaneously. With the same brittleness of the reservoir, fracture in zipper fracturing always propagates slightly longer than that in simultaneous fracturing. The results provide theoretical support for both zipper and synchronous fracturing modes, which helps design well completion and fracturing operation parameters. Zheng [14,15] proposed simulation of the propagation of multiple closely spaced fractures, and simulation results show that closely spaced perforation clusters in a stage tend to be unevenly and asymmetrically initiated, and a coalesce situation likely occurs when fracture space is too small. Zhang [16] verified this mechanism by a physical experiment. In multiple vertical wells, 3D fractured areas produced in simultaneous fracturing are much smaller than those produced in sequential fracturing [17]. Saberhosseini [18] used an extended finite-element model by ABAQUS to study the multiple fracture growth in modified zipper fracturing, proposing that we should pay attention to the flow rates and the injection to prevent multiple hydraulic fracture deviation and collision. Zhou [6] simulated the periodic variation law of induced stress in modified zipper fracturing, and proposed that arranging a perforating position through fracturing units can effectively increase the reservoir reconstruction area. In addition, at the same time, Sukumar [19] quantified the degree of pressure communication between adjacent reservoirs based on field data for the first time. Some three-dimensional and fully coupled models were proposed to study the interaction between multi-fractures [20,21,22].

However, in numerical simulation, the interference model between closed fracture and dynamic propagation fracture still has some deficiency. Whether in simultaneous or sequential fracturing simulation, the reverse support of proppant on compression fractures is often ignored [15,23,24,25]. Based on the boundary element method, Li [26] proposed nonlinear normal stiffness and studied the propagation of rock joint under compression. However, the rock joint model was proposed based on a natural fracture; it is not suitable for HFs. In general, hydraulic fractures are supported by two forces during fracturing: pressure from the injected fluid and reactive force from proppant. The above rock joint model did not consider the initial aperture of HF caused by hydraulic fracturing. Cheng [27] studied the simultaneous and sequential fracturing while considering the nonlinear closure of the fracture. However, in his simulation model, the shrinkage of HF aperture is calculated on the basis of the contraction in the previous step, which cannot reflect the superposition of proppant counterforce. In this paper, a mathematical model is established, taking the comprehensive change of the fracture aperture as unknown in order to reflect the superposition effect of the counter-support force. In addition to this, there are still some fracture interaction mechanisms that have not been revealed, such as fracture morphology when normal stiffness is considered in zipper fracturing.

In conclusion, the aim of this paper is to establish a multi-fracture propagation interaction model considering fracture plane support and study the interaction mechanism between hydraulic fractures. The rest of this paper is organized into the following four sections. Section 2 is the construction of the computational model and relevant validation. Numerical simulation analysis and discussions are presented in Section 3 to study the aperture and trajectory of HFs, the induced stress distribution situation. Section 4 presents some conclusions.

2. Numerical Model

2.1. Model Assumption

In reality, the propagation of hydraulic fracture is the combined effect of multiple physical fields. At present, it is difficult to establish a simulation model coupling all physical fields. Thus, in order to simplify the calculation model, and focus our research on fracture morphology and stress interference in the horizontal well, the following assumptions are introduced to establish the numerical model: (1) Isothermal conditions, the influence of temperature on fracture deformation during fracture propagation is ignored; (2) The model is homogeneous and isotropic, that is, do not consider the heterogeneity of rock mechanical properties and abnormal pressure caused by geological tectonic movement; (3) The target formation belongs to linearly elastic medium, that is, the load and deformation obey Hooke’s law [28]; (4) Fracture belongs to the KGD model, that is, HF is a vertical fracture with constant height, and the horizontal plane satisfies the plane strain condition for a short fracture relative to height [29]; (5) The study domain is infinite and ignores the abnormal pressure caused by geological tectonic movement; (6) The fracture propagation process follows the theory of linear elastic fracture mechanics [30], that is, the linear theory of elasticity is used to analyze the crack, and some characteristic parameters, like stress intensity factor, are used to judge the crack growth; (7) The fluid within the fracture is incompressible and belongs to the Poiseuille flow between two parallel plates [31].

2.2. Discontinuous Displacement Method

The displacement discontinuity method (DDM) is one of the boundary element methods (BEM) and was first proposed by Crouch [32]. In this paper, we use DDM to simulate the deformation of fracture and induced stress field. Separate HF into N segments of discontinuous microelements. When a discrete fracture microelement is subjected to field stress in the discrete state, relative displacement between the upper and lower surfaces will occur. In general, the shear stress and normal stress in each fracture element are as follows [33]:

$$\{\begin{array}{c}{\sigma }_s^i={\displaystyle \sum }_{j=1}^N{A}_{ss}^{i,j}{D}_s^j+{\displaystyle \sum }_{j=1}^N{A}_{sn}^{i,j}{D}_n^j \\ {\sigma }_n^i={\displaystyle \sum }_{j=1}^N{A}_{ns}^{i,j}{D}_s^j+{\displaystyle \sum }_{j=1}^N{A}_{nn}^{i,j}{D}_n^j \end{array}\left(i,j=1,2,\dots N\right)$$

(1)

where $D_s,D_n$ is the normal and shear relative displacements of the discontinuous microelement in the s–n coordinate system, m; ${\mathrm{A}}_{ss}^{i,j},{\mathrm{A}}_{sn}^{i,j},{\mathrm{A}}_{nn}^{i,j},{\mathrm{A}}_{ns}^{i,j}$ are the influence coefficients of unit i under the discontinuous displacement of unit j, Pa/m; ${\sigma }_s^i, {\sigma }_n^i$ are the shear stress and normal stress in the fracture element, respectively, Pa. The impact coefficients are as follows:

$$\begin{array}{c}{A}_{ss}^{i,j}=2G[-f_{{}^{\prime }\overline{xy}}\sin 2\mathrm{\gamma }-f_{{}^{\prime }\overline{yy}}\cos 2\mathrm{\gamma }-\overline{y}(f_{{}^{\prime }\overline{xyy}}\sin 2\mathrm{\gamma }-f_{{}^{\prime }\overline{yyy}}\cos 2\mathrm{\gamma }\left)\right]\\ A_{sn}^{i,j}=2G[-y(f_{{}^{\prime }\overline{xyy}}\cos 2\mathrm{\gamma }+f_{{}^{\prime }\overline{yyy}}\sin 2\mathrm{\gamma }\left)\right]\\ A_{nn}^{i,j}=2G[2f_{{}^{\prime }\overline{xy}}\sin 2\mathrm{\gamma }-f_{{}^{\prime }\overline{yy}}\sin 2\mathrm{\gamma }-\overline{y}(f_{{}^{\prime }\overline{xyy}}\cos 2\mathrm{\gamma }-f_{{}^{\prime }\overline{yyy}}\sin 2\mathrm{\gamma }\left)\right]\\ A_{ns}^{i,j}=2G[-f_{{}^{\prime }\overline{yy}}+\overline{y}(f_{{}^{\prime }x\overline{yy}}\sin 2\mathrm{\gamma }-f_{{}^{\prime }\overline{yyy}}\cos 2\mathrm{\gamma }\left)\right] \end{array}$$

(2)

where f_′x, f_′y, f_′xx, f_′xy, f_′xyy, f_′yyy, are the derivatives of function f(x,y), respectively; ν is the Poisson ratio; G is the shear modulus, MPa. For a detailed description of rock deformation, please refer to our previous studies [15].

Put Equation (1) into matrix form:

$$\left[\begin{array}{cc}\overline{A_{ss}}& \overline{A_{sn}}\\ \overline{A_{ns}}& \overline{A_{nn}}\end{array}\right]\left[\begin{array}{c}\overline{D_s}\\ \overline{D_n}\end{array}\right]=\left[\begin{array}{c}\overline{{\sigma }_s}\\ \overline{{\sigma }_n}\end{array}\right]$$

(3)

By eliminating $\overline{D_s}$ in Equation (3), the width of HF ($\overline{D_n}$) can be directly solved:

$$\overline{A_{ns}}\cdot {\overline{A_{ss}}}^{\left(-1\right)}\cdot \overline{{\sigma }_s}-\overline{A_{ns}}\cdot {\overline{A_{ss}}}^{\left(-1\right)}\cdot \overline{A_{sn}}\cdot \overline{D_n}+\overline{A_{nn}}\cdot \overline{D_n}-\overline{{\sigma }_n}=0$$

(4)

where:

$$\begin{array}{c}\overline{A_{ss}}=\left[\begin{array}{ccc}\begin{array}{c}\begin{array}{cc}{A}_{ss}^{11}& A_{ss}^{12}\end{array}\\ \begin{array}{cc}{A}_{ss}^{21}& A_{ss}^{22}\end{array}\end{array}& \cdots & \begin{array}{c}{A}_{ss}^{1N}\\ A_{ss}^{2N}\end{array}\\ ⋮& \ddots & ⋮\\ \begin{array}{cc}{A}_{ss}^{N1}& A_{ss}^{N2}\end{array}& \cdots & A_{ss}^{NN}\end{array}\right] \overline{A_{sn}}=\left[\begin{array}{ccc}\begin{array}{c}\begin{array}{cc}{A}_{sn}^{11}& A_{sn}^{12}\end{array}\\ \begin{array}{cc}{A}_{sn}^{21}& A_{sn}^{22}\end{array}\end{array}& \cdots & \begin{array}{c}{A}_{sn}^{1N}\\ A_{sn}^{2N}\end{array}\\ ⋮& \ddots & ⋮\\ \begin{array}{cc}{A}_{sn}^{N1}& A_{sn}^{N2}\end{array}& \cdots & A_{sn}^{NN}\end{array}\right]\overline{D_s}=\left[\begin{array}{c}{D}_s^1\\ \begin{array}{c}{D}_s^2\\ ⋮\end{array}\\ D_s^N\end{array}\right]\\ \overline{A_{ns}}=\left[\begin{array}{ccc}\begin{array}{c}\begin{array}{cc}{A}_{ns}^{11}& A_{ns}^{12}\end{array}\\ \begin{array}{cc}{A}_{ns}^{21}& A_{ns}^{22}\end{array}\end{array}& \cdots & \begin{array}{c}{A}_{ns}^{1N}\\ A_{ns}^{2N}\end{array}\\ ⋮& \ddots & ⋮\\ \begin{array}{cc}{A}_{ns}^{N1}& A_{ns}^{N2}\end{array}& \cdots & A_{ns}^{NN}\end{array}\right]\overline{A_{nn}}=\left[\begin{array}{ccc}\begin{array}{c}\begin{array}{cc}{A}_{nn}^{11}& A_{nn}^{12}\end{array}\\ \begin{array}{cc}{A}_{nn}^{21}& A_{nn}^{22}\end{array}\end{array}& \cdots & \begin{array}{c}{A}_{nn}^{1N}\\ A_{nn}^{2N}\end{array}\\ ⋮& \ddots & ⋮\\ \begin{array}{cc}{A}_{nn}^{N1}& A_{nn}^{N2}\end{array}& \cdots & A_{nn}^{NN}\end{array}\right]\overline{D_n}=\left[\begin{array}{c}{D}_n^1\\ \begin{array}{c}{D}_n^2\\ ⋮\end{array}\\ D_n^N\end{array}\right]\\ \overline{{\sigma }_s}=\left[\begin{array}{c}-\frac{1}{2}\left({\sigma }_h-{\sigma }_H\right)\sin 2{\mathsf{\theta }}_1\\ \begin{array}{c}-\frac{1}{2}\left({\sigma }_h-{\sigma }_H\right)\sin 2{\mathsf{\theta }}_2\\ ⋮\end{array}\\ -\frac{1}{2}\left({\sigma }_h-{\sigma }_H\right)\sin 2{\mathsf{\theta }}_{\mathrm{N}}\end{array}\right] \overline{{\sigma }_n}=\left[\begin{array}{c}{p}_{f1}-{\sigma }_h{\sin }^2{\mathsf{\theta }}_1-{\sigma }_H{\cos }^2{\mathsf{\theta }}_1\\ \begin{array}{c}{p}_{f2}-{\sigma }_h{\sin }^2{\mathsf{\theta }}_2-{\sigma }_H{\cos }^2{\mathsf{\theta }}_2\\ ⋮\end{array}\\ p_{fN}-{\sigma }_h{\sin }^2{\mathsf{\theta }}_{\mathrm{N}}-{\sigma }_H{\cos }^2{\mathsf{\theta }}_{\mathrm{N}}\end{array}\right]\end{array}$$

where P_f1 is the injection pressure, Pa. θ is the angle between discrete element and maximum horizontal principal stress, °.

Then, the specific stresses and displacements of each point in the reservoir caused by j-th microelements can be obtained as follows [34]:

$$\{\begin{array}{c}{\mathrm{u}}^j{}_x=D_{\bar{x}}[2(1-v)f_{{}^{\prime }\bar{y}}-y{f}_{{}^{\prime }\overline{xx}}]+D_{\bar{y}}[-(1-2v)f_{{}^{\prime }\bar{x}}-y{f}_{{}^{\prime }\overline{xy}}]\\ {\mathrm{u}}^j{}_y=D_{\bar{x}}[2(1-v)f_{{}^{\prime }\bar{x}}-y{f}_{{}^{\prime }\overline{xy}}]+D_{\bar{y}}[2(1-v)f_{{}^{\prime }\bar{y}}-y{f}_{{}^{\prime }\overline{yy}}]\\ {{\sigma }^j}_{xx}=2G{D}_{\bar{x}}(2f_{{}^{\prime }\overline{xy}}+y{f}_{{}^{\prime }\overline{xyy}})+2G{D}_{\bar{y}}(f_{{}^{\prime }\overline{yy}}+y{f}_{{}^{\prime }\overline{yyy}})\\ {{\sigma }^j}_{yy}=2G{D}_{\bar{x}}(-y{f}_{{}^{\prime }\overline{xyy}})+2G{D}_{\bar{y}}(f_{{}^{\prime }\overline{yy}}+y{f}_{{}^{\prime }\overline{yyy}})\\ {{\sigma }^j}_{xy}=2G{D}_{\bar{x}}(f_{{}^{\prime }\overline{yy}}+y{f}_{{}^{\prime }\overline{xyy}})+2G{D}_y(-y{f}_{{}^{\prime }\overline{yyy}})\end{array}$$

(5)

where ${\sigma }^j{}_{xx},{\sigma }^j{}_{yy},{\sigma }^j{}_{xy},{\mathrm{u}}^j{}_x,{\mathrm{u}}^j{}_y$ are the induced stress and displacements.

2.3. Deformation of Hydraulic Fracture

Whether there is staged fracturing in a single well or zipper fracturing in a double well, it can be considered that hydraulic fractures are opened sequentially. It is known from microseismic monitoring [35] and experiments [16] that, in hydraulic fracturing, hydraulic fracture clusters will intersect and eventually form a main hydraulic fracture zone. Thus, in order to simplify the computational model, we assume that each fracturing stage includes one perforation cluster and creates only one fracture at a time, starting from the toe to the heel. The above hypothesis can be found in the literature [36].

Under the in-situ stress, the previous opening HF after fracturing construction has the tendency to close. The existence of proppant can make HF maintain a certain aperture to improve the oil and gas seepage capacity in the fracture. The opening degree is related to the mechanical properties of the rock. Assume that HFs after hydraulic fracturing are fully filled with proppant. Goodman [37] first proposed the normal stiffness, K_n, and tangential stiffness, K_s, which were defined as the ratio of the increment of normal stress and shear stress, respectively, to the related increment of displacement. Zhou [38] tested the stiffness with stress wave propagation. In this paper, we introduce the support stiffness coefficients $K_s$ and $K_n$ to simulate the reverse supporting force from proppant on the fracture surface [33,36,39]. Its expression is as follows:

$$\{\begin{array}{c}{\left({\mathsf{\sigma }}_n\right)}_{stiff}=-K_n\mathsf{\Delta }{D}_n\\ {\left({\mathsf{\sigma }}_s\right)}_{stiff}=-K_s\mathsf{\Delta }{D}_s\end{array}$$

(6)

Each HF is divided in N_i (i represent different fracturing stage) fracture elements. At the initial stage of fracturing, there is only one fracture, and the fracture surface is affected by both in-situ stress and fluid pressure inside the fracture. Considering that the model size is infinite, the stress equilibrium equation on each fracture element is given by:

$$\{\begin{array}{c}{\displaystyle \sum }_{j=1}^{N_1}\left(A_{ss}^{i,j}{\left[D_s^j\right]}_{max}^1+A_{sn}^{i,j}{\left[D_n^j\right]}_{max}^1\right)={\mathsf{\sigma }}_{\mathrm{s}}^{\mathrm{i}}\\ {\displaystyle \sum }_{j=1}^{N_1}\left(A_{ns}^{i,j}{\left[D_s^j\right]}_{max}^1+A_{sn}^{i,j}{\left[D_n^j\right]}_{max}^1\right)={\mathsf{\sigma }}_{\mathrm{n}}^{\mathrm{i}}\end{array}\left(i=1,2\dots N_1\right)$$

(7)

The boundary influence coefficients $A_{ss}^{i,j}$, etc., in the above are obtained from Equation (2). The total stress at element i is obtained by adding the in-situ stress and the injection pressure in fracture:

$$\{\begin{array}{c}{\mathsf{\sigma }}_s^{\mathrm{i}}={\left({\mathsf{\sigma }}_s^{\mathrm{i}}\right)}_p+\left({\mathsf{\sigma }}_s^{\mathrm{i}}\right){ }_0^{\infty }\\ {\mathsf{\sigma }}_{\mathrm{n}}^{\mathrm{i}}{\mathsf{\sigma }}_{\mathrm{n}}^{\mathrm{i}}={\left({\mathsf{\sigma }}_{\mathrm{n}}^{\mathrm{i}}\right)}_p+\left({\mathsf{\sigma }}_{\mathrm{n}}^{\mathrm{i}}\right){ }_0^{\infty }\end{array} \left(i=1,2\dots N_1\right)$$

(8)

where $\left({\mathsf{\sigma }}_s^{\mathrm{i}}{\mathsf{\sigma }}_{\mathrm{s}}^{\mathrm{i}}\right){ }_0^{\infty }, \left({\mathsf{\sigma }}_n^{\mathrm{i}}\right){ }_0^{\infty }$ are far-field stress, Pa; ${\left({\mathsf{\sigma }}_s^{\mathrm{i}}\right)}_p, {\left({\mathsf{\sigma }}_n^{\mathrm{i}}\right)}_p$ are stress inside of the fracture, Pa; the s, n are the local coordinate system about the ith element.

Generally, the stress components are as follows:

$$\{\begin{array}{c}{\left({\mathsf{\sigma }}_{\mathrm{s}}^{\mathrm{i}}\right)}_p=0; \\ {\left({\mathsf{\sigma }}_{\mathrm{n}}^{\mathrm{i}}\right)}_p=P_f;\end{array} \{\begin{array}{c}\left({\mathsf{\sigma }}_{\mathrm{s}}^{\mathrm{i}}\right){ }_0^{\infty }=-\frac{1}{2}\left({\mathsf{\sigma }}_{\mathrm{h}}-{\mathsf{\sigma }}_{\mathrm{H}}\right)\sin 2{\mathsf{\beta }}^{\mathrm{i}}; \\ \left({\mathsf{\sigma }}_{\mathrm{n}}^{\mathrm{i}}\right){ }_0^{\infty }=-{\mathsf{\sigma }}_{\mathrm{h}}{\sin }^2{\mathsf{\theta }}_{\mathrm{i}}-{\mathsf{\sigma }}_{\mathrm{H}}{\cos }^2{\mathsf{\beta }}^{\mathrm{i}};\end{array}$$

(9)

During this hydraulic fracturing stage, the discontinuous displacements (${\left[D_s^j\right]}_0, {\left[D_n^j\right]}_0$) of the first HF can be obtained by Equations (6)–(8). Stop injection and shut off the pump when the first HF reaches the preset length. In addition, now under the action of in-situ stress, the aperture of HF may be shrunk compared with the pump not being shut off. However, due to the existence of proppant, the fracture surface will generate counter support force from the inside to the outside. The total stress at element i is obtained by adding the in-situ stress, the injection pressure in fracture, and the reaction force from proppant. Thus, the fracture aperture after the end of first fracturing stage is as follows:

$$\{\begin{array}{c}{\sum }_{j=1}^{N_1}(A_{ss}^{i,j}{\left[D_s^j\right]}_{\mathit{max}}^1+A_{sn}^{i,j}{\left[D_n^j\right]}_{\mathit{max}}^1)+{\sum }_{j=1}^{N_1}(A_{ss}^{i,j}\mathsf{\Delta }{D}_s^j+A_{sn}^{i,j}\mathsf{\Delta }{D}_n^j)\\ =-\frac{1}{2}({\mathsf{\sigma }}_{\mathrm{h}}-{\mathsf{\sigma }}_{\mathrm{H}})\sin 2{\mathsf{\theta }}_{\mathrm{i}}-K_s\mathsf{\Delta }{D}_s^i\\ {\sum }_{j=1}^{N_1}(A_{ns}^{i,j}{\left[D_s^j\right]}_{\mathit{max}}^1+A_{nn}^{i,j}{\left[D_n^j\right]}_{\mathit{max}}^1)+{\sum }_{j=1}^{N_1}(A_{ns}^{i,j}\mathsf{\Delta }{D}_s^j+A_{nn}^{i,j}\mathsf{\Delta }{D}_n^j)\\ =-{\mathsf{\sigma }}_{\mathrm{h}}{\sin }^2{\mathsf{\theta }}_{\mathrm{i}}-{\mathsf{\sigma }}_{\mathrm{H}}{\cos }^2{\mathsf{\theta }}_{\mathrm{i}}-K_n\mathsf{\Delta }{D}_n^i\end{array}(i=1,2\dots N_1)$$

(10)

where [D]_max represents the maximum discontinuous displacement of the node when fracture surface opens, m. Introduce the total stress and deformations of element:

$$\{\begin{array}{c}{\left[D_s^j\right]}^{{}^{\prime }}={\left[D_s^j\right]}_{max}^1+\mathsf{\Delta }{D}_s^j;\\ {\left[D_n^j\right]}^{{}^{\prime }}={\left[D_n^j\right]}_{max}^1+\mathsf{\Delta }{D}_n^j;\end{array} \{\begin{array}{c}{\left({\mathsf{\sigma }}_{\mathrm{s}}^{\mathrm{i}}\right)}^{{}^{\prime }}=\left({\mathsf{\sigma }}_{\mathrm{s}}^{\mathrm{i}}\right){ }_0^{\infty }+{\left({\mathsf{\sigma }}_{\mathrm{s}}^{\mathrm{i}}\right)}_{stiff}\\ {\left({\mathsf{\sigma }}_{\mathrm{n}}^{\mathrm{i}}\right)}^{{}^{\prime }}=\left({\mathsf{\sigma }}_{\mathrm{n}}^{\mathrm{i}}\right){ }_0^{\infty }+{\left({\mathsf{\sigma }}_{\mathrm{n}}^{\mathrm{i}}\right)}_{stiff}\end{array}$$

(11)

where ${\left({\mathsf{\sigma }}_{\mathrm{s}}^{\mathrm{i}}\right)}_{stiff}, {\left({\mathsf{\sigma }}_{\mathrm{n}}^{\mathrm{i}}\right)}_{stiff}$ are what we will call the induced stress, Pa; $\mathsf{\Delta }{D}_s^j$, $\mathsf{\Delta }{D}_n^j$ are the induced displacement, m.

Then, the deformation equation is as below:

$$\{\begin{array}{c}{\displaystyle \sum }_{j=1}^{N_1}\left(A_{ss}^{i,j}{\left[D_s^j\right]}^{{}^{\prime }}+A_{sn}^{i,j}{\left[D_n^j\right]}^{{}^{\prime }}\right)={\left({\mathsf{\sigma }}_{\mathrm{s}}^{\mathrm{i}}\right)}^{{}^{\prime }}\\ {\displaystyle \sum }_{j=1}^{N_1}\left(A_{ns}^{i,j}{\left[D_s^j\right]}^{{}^{\prime }}+A_{sn}^{i,j}{\left[D_n^j\right]}^{{}^{\prime }}\right)={\left({\mathsf{\sigma }}_{\mathrm{n}}^{\mathrm{i}}\right)}^{{}^{\prime }}\end{array}$$

(12)

Now, the next stage fracturing is started again and the injection is pumped on until the second HF propagated to the preset length. During this stage, due to the existence of stress interaction, the previous HF will affect the subsequent HF and vice versa. In addition, at the same time, the previous aperture will also change. Then, the HF satisfies the following equation:

$$\{\begin{array}{c}{\sum }_{j=1}^{N_1}(A_{ss}^{i,j}{\left[D_s^j\right]}_{\mathit{max}}^1+A_{sn}^{i,j}{\left[D_n^j\right]}_{\mathit{max}}^1)+{\sum }_{j=1}^{N_1}(A_{sS}^{i,j}\mathsf{\Delta }{D}_s^j+A_{sn}^{i,j}\mathsf{\Delta }{D}_n^j)\\ +{\sum }_{j=N_1+1}^{N_1+N_1}[A_{ss}^{i,j}{\left[D_s^j\right]}_{\mathit{max}}^2+A_{sn}^{i,j}{\left[D_n^j\right]}_{\mathit{max}}^2]={\left({\mathsf{\sigma }}_{\mathrm{s}}^{\mathrm{i}}\right)}^{{}^{″}}\\ {\sum }_{j=1}^{N_1}(A_{ns}^{i,j}{\left[D_s^j\right]}_{\mathit{max}}^1+A_{sn}^{i,j}{\left[D_n^j\right]}_{\mathit{max}}^1)+{\sum }_{j=1}^{N_1}(A_{ns}^{i,j}\mathsf{\Delta }{D}_s^j+A_{nn}^{i,j}\mathsf{\Delta }{D}_n^j)\\ +{\sum }_{j=N_1+1}^{N_1+N_1}[A_{ns}^{i,j}{\left[D_s^j\right]}_{\mathit{max}}^2+A_{nn}^{i,j}{\left[D_n^j\right]}_{\mathit{max}}^2]={\left({\mathsf{\sigma }}_{\mathrm{n}}^{\mathrm{i}}\right)}^{{}^{″}}\end{array}(i=1,2\dots N_1+N_2)$$

(13)

For propped HFs, the expression is similar to Formula (11):

$$\{\begin{array}{c}{\left({\mathsf{\sigma }}_{\mathrm{s}}^{\mathrm{i}}\right)}^{{}^{″}}=\left({\mathsf{\sigma }}_{\mathrm{s}}^{\mathrm{i}}\right){ }_0^{\infty }+{\left({\mathsf{\sigma }}_{\mathrm{s}}^{\mathrm{i}}\right)}_{stiff}\\ {\left({\mathsf{\sigma }}_{\mathrm{n}}^{\mathrm{i}}\right)}^{{}^{″}}=\left({\mathsf{\sigma }}_{\mathrm{n}}^{\mathrm{i}}\right){ }_0^{\infty }+{\left({\mathsf{\sigma }}_{\mathrm{n}}^{\mathrm{i}}\right)}_{stiff}\end{array} \left(i=1, 2, \dots , N_1\right)$$

(14)

For unpropped HFs, the expression on the right side of the formula is similar to Formula (8) and is as below:

$$\{\begin{array}{c}{\left({\mathsf{\sigma }}_s^{\mathrm{i}}\right)}^{{}^{″}}={\left({\mathsf{\sigma }}_s^{\mathrm{i}}\right)}_p+\left({\mathsf{\sigma }}_s^{\mathrm{i}}\right){ }_0^{\infty }\\ {\left({\mathsf{\sigma }}_n^{\mathrm{i}}\right)}^{{}^{″}}={\left({\mathsf{\sigma }}_n^{\mathrm{i}}\right)}_p+\left({\mathsf{\sigma }}_n^{\mathrm{i}}\right){ }_0^{\infty }\end{array} \left(i=N_1+1, N_1+2, \dots , N_1+N_2\right)$$

(15)

To simplify the theoretical model, a total discontinuous displacement value is introduced:

$$\{\begin{array}{c}{\left[D_{ℝ}^j\right]}^{{}^{″}}={\left[D_{ℝ}^j\right]}_{max}^2+\mathrm{\Delta }{D}_{\mathbb{R}}^j;\left(j=1, 2, \dots , N_1\right) \\ {\left[D_{ℝ}^j\right]}^{″}=D_{ℝ}^j;\left(i=N_1+1, N_1+2, \dots , N_1+N_2\right) \end{array}\left(ℝ=s,n\right)$$

(16)

Then, Formula (13) can be abbreviated as the following form:

$$\{\begin{array}{c}{\displaystyle \sum }_{j=1}^{N_1+N_2}\left(A_{ss}^{i,j}{\left[D_s^j\right]}^{″}+A_{sn}^{i,j}{\left[D_n^j\right]}^{″}\right)={\left({\mathsf{\sigma }}_{\mathrm{s}}^{\mathrm{i}}\right)}^{″}\\ {\displaystyle \sum }_{j=1}^{N_1+N_2}\left(A_{ns}^{i,j}{\left[D_s^j\right]}^{″}+A_{sn}^{i,j}{\left[D_n^j\right]}^{″}\right)={\left({\mathsf{\sigma }}_{\mathrm{n}}^{\mathrm{i}}\right)}^{″}\end{array} \left(i=1, 2, \dots , N_1+N_2\right)$$

(17)

At this step, the unknown quantities are the variation of the propped fracture and the maximum width of the unpropped fracture. Now, shot off the second fracturing well, the formula about the aperture of HFs after reaching equilibrium is as below:

$$\{\begin{array}{c}{\sum }_{j=1}^{N_1+N_2}(A_{SS}^{i,j}{\left[D_s^j\right]}_{\mathit{max}}^{1,2}+A_{sn}^{i,j}{\left[D_n^j\right]}_{\mathit{max}}^{1,2})+{\sum }_{j=1}^{N_1+N_2}(A_{ss}^{i,j}\mathsf{\Delta }{D}_s^j+A_{sn}^{i,j}\mathsf{\Delta }{D}_n^j)\\ =-\frac{1}{2}({\mathsf{\sigma }}_{\mathrm{h}}-{\mathsf{\sigma }}_{\mathrm{H}})\sin 2{\mathsf{\theta }}_{\mathrm{i}}-K_s\mathsf{\Delta }{D}_s^i\\ {\sum }_{j=1}^{N_1+N_2}(A_{ns}^{i,j}{\left[D_s^j\right]}_{\mathit{max}}^{1,2}+A_{sn}^{i,j}{\left[D_n^j\right]}_{\mathit{max}}^{1,2})+{\sum }_{j=1}^{N_1+N_2}(A_{ns}^{i,j}\mathsf{\Delta }{D}_s^j+A_{nn}^{i,j}\mathsf{\Delta }{D}_n^j)\\ =-{\mathsf{\sigma }}_{\mathrm{h}}{\sin }^2{\mathsf{\theta }}_{\mathrm{i}}-{\mathsf{\sigma }}_{\mathrm{H}}{\cos }^2{\mathsf{\theta }}_{\mathrm{i}}-K_n\mathsf{\Delta }{D}_n^i\end{array}(i=1,2\dots N_1)$$

(18)

Thus, it can be seen that there are two forms about fracture deformation in sequential fracturing: one is to close the injection well, and there are only supporting fractures in the domain. Supposing that there are M fractures, the fracture deformation equation is as follows:

$$\{\begin{array}{c}{\displaystyle \sum _{j=1}^{\mathbb{N}}}(A_{ss}^{i,j}{\left[D_s^j\right]}_{\mathit{max}}^{1,2\dots M}+A_{sn}^{i,j}{\left[D_n^j\right]}_{\mathit{max}}^{1,2\dots M})+{\displaystyle \sum _{j=1}^{\mathbb{N}}}(A_{ss}^{i,j}\mathsf{\Delta }{D}_s^j+A_{sn}^{i,j}\mathsf{\Delta }{D}_n^j)={\left({\mathsf{\sigma }}_s^{\mathrm{i}}\right)}^{″}\\ {\displaystyle \sum _{j=1}^{\mathbb{N}}}(A_{ns}^{i,j}{\left[D_s^j\right]}_{\mathit{max}}^{1,2\dots M}+A_{sn}^{i,j}{\left[D_n^j\right]}_{\mathit{max}}^{1,2\dots M})+{\displaystyle \sum _{j=1}^{\mathbb{N}}}(A_{ns}^{i,j}\mathsf{\Delta }{D}_s^j+A_{nn}^{i,j}\mathsf{\Delta }{D}_n^j)={\left({\mathsf{\sigma }}_{\mathrm{n}}^{\mathrm{i}}\right)}^{″}\end{array}\left(\begin{array}{c}i=1,2\dots N; N=N_1+\dots +N_M\end{array}\right)$$

(19)

$$\{\begin{array}{c}{\left({\mathsf{\sigma }}_{\mathrm{s}}^{\mathrm{i}}\right)}^{″}=-\frac{1}{2}\left({\mathsf{\sigma }}_{\mathrm{h}}-{\mathsf{\sigma }}_{\mathrm{H}}\right)\sin 2{\mathsf{\theta }}_{\mathrm{i}}-K_s\mathsf{\Delta }{D}_s^i \\ {\left({\mathsf{\sigma }}_{\mathrm{n}}^{\mathrm{i}}\right)}^{″}=-{\mathsf{\sigma }}_{\mathrm{h}}{\sin }^2{\mathsf{\theta }}_{\mathrm{i}}-{\mathsf{\sigma }}_{\mathrm{H}}{\cos }^2{\mathsf{\theta }}_{\mathrm{i}}-K_n\mathsf{\Delta }{D}_n^i\end{array} \left(\begin{array}{c}i=1,2\dots \mathbb{N}; \mathbb{N}=N_1+\dots +N_M\end{array}\right)$$

(20)

The other is that: the previous M-1 HFs are under the compression of in-situ stress and induce stress, and the Mth HF is the current dynamic propagation fracture.

$$\{\begin{array}{c}{\displaystyle \sum _{j=1}^{\mathbb{N}-N_M}}(A_{ss}^{i,j}{\left[D_s^j\right]}_{\mathit{max}}^{1,2\dots M-1}+A_{sn}^{i,j}{\left[D_n^j\right]}_{\mathit{max}}^{1,2\dots M-1})+{\displaystyle \sum _{j=1}^{\mathbb{N}-N_M}}(A_{ss}^{i,j}\mathsf{\Delta }{D}_s^j+A_{sn}^{i,j}\mathsf{\Delta }{D}_n^j)\\ +{\displaystyle \sum _{j=\mathbb{N}-N_M+1}^{\mathbb{N}}}[A_{SS}^{i,j}{\left[D_s^j\right]}_{\mathit{max}}^{1,2\dots M-1}+A_{sn}^{i,j}{\left[D_n^j\right]}_{\mathit{max}}^{1,2\dots M-1}]={\left({\mathsf{\sigma }}_{\mathrm{s}}^{\mathrm{i}}\right)}^{″}\\ {\displaystyle \sum _{j=1}^{\mathbb{N}-N_M}}(A_{ns}^{i,j}{\left[D_s^j\right]}_{\mathit{max}}^{1,2\dots M-1}+A_{sn}^{i,j}{\left[D_n^j\right]}_{\mathit{max}}^{1,2\dots M-1})+{\displaystyle \sum _{j=1}^{\mathbb{N}-N_M}}(A_{ns}^{i,j}\mathsf{\Delta }{D}_s^j+A_{nn}^{i,j}\mathsf{\Delta }{D}_n^j)\\ +{\displaystyle \sum _{j=\mathbb{N}-N_M+1}^{\mathbb{N}}}[A_{ns}^{i,j}{\left[D_s^j\right]}_{\mathit{max}}^M+A_{nn}^{i,j}{\left[D_n^j\right]}_{\mathit{max}}^M]={\left({\mathsf{\sigma }}_{\mathrm{n}}^{\mathrm{i}}\right)}^{″}\end{array}\left(\begin{array}{c}i=1,2\dots \mathbb{N}; \mathbb{N}=N_1+\dots +N_M\end{array}\right)$$

(21)

$$\{\begin{array}{c}{\left({\mathsf{\sigma }}_{\mathrm{s}}^{\mathrm{i}}\right)}^{″}=-\frac{1}{2}\left({\mathsf{\sigma }}_{\mathrm{h}}-{\mathsf{\sigma }}_{\mathrm{H}}\right)\sin 2{\mathsf{\theta }}_{\mathrm{i}}-K_s\mathsf{\Delta }{D}_s^i \\ {\left({\mathsf{\sigma }}_{\mathrm{n}}^{\mathrm{i}}\right)}^{″}=-{\mathsf{\sigma }}_{\mathrm{h}}{\sin }^2{\mathsf{\theta }}_{\mathrm{i}}-{\mathsf{\sigma }}_{\mathrm{H}}{\cos }^2{\mathsf{\theta }}_{\mathrm{i}}-K_n\mathsf{\Delta }{D}_n^i\end{array} \left(i=1, 2, \dots , \mathbb{N}-N_M\right)$$

(22)

$$\{\begin{array}{c}{\left({\mathsf{\sigma }}_{\mathrm{s}}^{\mathrm{i}}\right)}^{″}=-\frac{1}{2}\left({\mathsf{\sigma }}_{\mathrm{h}}-{\mathsf{\sigma }}_{\mathrm{H}}\right)\sin 2{\mathsf{\theta }}_{\mathrm{i}} \\ {\left({\mathsf{\sigma }}_{\mathrm{n}}^{\mathrm{i}}\right)}^{″}=P_f-{\mathsf{\sigma }}_{\mathrm{h}}{\sin }^2{\mathsf{\theta }}_{\mathrm{i}}-{\mathsf{\sigma }}_{\mathrm{H}}{\cos }^2{\mathsf{\theta }}_{\mathrm{i}}\end{array} \left(i=\mathbb{N}-N_M+1, N_1+2, \dots , \mathbb{N}\right)$$

(23)

2.4. The Initiation and Propagation of Fracture

The stress intensity factor (SIF) is an important parameter to describe the field stress of the fracture tip, which is calculated by discontinuous displacement value on the fracture surface. The SIF in point a of the fracture tip can be shown [40]:

$$\{\begin{array}{c}{K}_{\mathsf{I}}=\frac{0.806E\sqrt{\pi }}{4\left(1-{\upsilon }^2\right)\sqrt{2a}}{D}_n\\ K_{\mathsf{I}\mathsf{I}}=\frac{0.806E\sqrt{\pi }}{4\left(1-{\upsilon }^2\right)\sqrt{2a}}{D}_s\end{array}$$

(24)

where $K_{\mathsf{I}}$ and $K_{\mathsf{I}\mathsf{I}}$ is the stress intensity factor of different types of cracks, MPa∙m^0.5; E is Young’s modulus, MPa; ν is Poisson’s ratio, and a is the half-length of crack, m.

At the same time, the maximum circumferential stress (MCS) criterion is introduced to judge the initiation and propagation direction of fracture [41]. The fracture deflection angle (θ₀) can be obtained by solving the first-order partial derivative for the circumferential stress:

$${\theta }_0=\arcsin \left[\frac{K_{\mathsf{I}\mathsf{I}}\left(K_{\mathsf{I}}+3\sqrt{K_{\mathsf{I}}{}^2+8K_{\mathsf{I}\mathsf{I}}{}^2}\right)}{K_{\mathsf{I}}{}^2+9K_{\mathsf{I}\mathsf{I}}{}^2}\right]$$

(25)

The mixed mode of the stress intensity factor is used to identify the initiation of HF and NFs:

$$K_e=\cos \frac{{\theta }_0}{2}\left[K_{\mathsf{I}}{\cos }^2\frac{\theta }{2}-\frac{3}{2}{K}_{\mathsf{I}\mathsf{I}}\sin {\theta }_0\right]$$

(26)

When $K_{IC}>{\mathrm{K}}_e$, the fracture starts to propagate, where K_IC is the fracture toughness, MPa∙m^1/2.

2.5. Flowchart

The flow chart is shown in Figure 1. The solution procedure is as follows: ① Input initial parameters, such as maximum horizontal stress, minimum horizontal stress, Young’s modulus, Poisson’s ratio, fracture toughness, fracture number, fracture spacing, injection pressure, perforation depth, and perforating angle; ② Hydraulic fracture discretization and judge well condition; ③ Calculate the D_n and D_s of HF until the accuracy requirements are met by Equations (19)–(23) based on different well condition; ④ Check the iterative step, if the iterative step is greater than the preset value, calculate induced stress (Equation (5)) and exit the calculation, otherwise judge whether HF initiation occurs (Equation (26)); ⑤ If initiation occurs, calculate the deflection angle of the element (Equation (25)), increase fracture units and iterative step, and then go to step ②, else go to step ② directly.

2.6. Model Validation

Firstly, we verify the induced stress caused by fracture deformation. As the calculation model meets the plane strain condition [30], the analytical solution of the induced field stress around the fracture can be expressed as [42]:

$$\{\begin{array}{c}{\mathsf{\sigma }}_{xx}=-P_f\frac{r}{a}{\left(\frac{a^2}{r_1{r}_2}\right)}^{\frac{3}{2}}\sin \mathsf{\theta }\sin \left[\frac{3}{2}({\mathsf{\theta }}_1+{\mathsf{\theta }}_2)\right]-P_f[\frac{r}{r_1{r}_2}\cos -\left(\mathsf{\theta }\frac{{\mathsf{\theta }}_1+{\mathsf{\theta }}_2}{2}\right)-1]\\ {\mathsf{\sigma }}_{yy}=P_f\frac{r}{a}{\left(\frac{a^2}{r_1{r}_2}\right)}^{\frac{3}{2}}\sin \mathsf{\theta }\sin \left[\frac{3}{2}({\mathsf{\theta }}_1+{\mathsf{\theta }}_2)\right]-P_f[\frac{r}{r_1{r}_2}\cos (\mathsf{\theta }-\frac{{\mathsf{\theta }}_1+{\mathsf{\theta }}_2}{2})-1]\\ {\mathsf{\sigma }}_{xy}=-P_f\frac{r}{a}{\left(\frac{a^2}{r_1{r}_2}\right)}^{\frac{3}{2}}\sin \mathsf{\theta }\cos \left[\frac{3}{2}({\mathsf{\theta }}_1+{\mathsf{\theta }}_2)\right]\end{array}$$

(27)

where P_f is the injection pressure, Pa; r, r₁, r₂ are respectively the distance between fracture and point Q, m; θ, θ₁, θ₂ are the angle between the x-axis and the line distance, °. The specific stress model is shown in Figure 2a:

Initial simulation parameters: injection pressure (P_f) is −3 MPa, the half-length of crack (a) is 1 m, the coordinates of point Q: x = 0.5 m, and y increases from 0 m to 10 m. Calculate the induced stress at point Q under a different vertical distance (coordinate y) from the fracture plane. The normal stress and shear stress obtained by simulation and calculation decrease with the increase of the distance from the fracture plane, and the numerical solution is basically consistent with the analytical solution. In Figure 2b, the solid line represents analytical results, and the dotted line represents the simulation results.

Secondly, we verify the stress intensity factor at the fracture tip. There is not only tensile failure, but also shear failure in the fracture, and the calculation model of the stress intensity factor at the fracture tip under the simple model is constructed (Figure 3). Assuming that there is a crack with a length of 2a in an infinitely uniform thin plate, the angle between fracture and the uniaxial tensile force P is θ. The calculation formula of the analytical solution of the stress intensity factor at the crack tip obtained by the Westergaard stress function method and the stress intensity factors of ${\mathrm{K}}_{\mathsf{I}}$ and ${\mathrm{K}}_{\mathsf{I}\mathsf{I}}$ can be expressed as:

$$\{\begin{array}{c}{\mathrm{K}}_{\mathsf{I}}=P\sqrt{\pi a}·{\sin }^2\mathsf{\theta } \\ {\mathrm{K}}_{\mathsf{I}\mathsf{I}}=P\sqrt{\pi a}·\sin \mathsf{\theta }·\cos \mathsf{\theta }\end{array}$$

(28)

In this paper, we calculated the stress intensity factors of ${\mathrm{K}}_{\mathsf{I}}$ and ${\mathrm{K}}_{\mathsf{I}\mathsf{I}}$ at the tip of the fracture, where the analytical results about ${\mathrm{K}}_{\mathsf{I}}$ and ${\mathrm{K}}_{\mathsf{I}\mathsf{I}}$ are respectively 1.3293 MPa∙m^1/2 and 0.7765 MPa∙m^1/2. In addition, at the same time, we provided the corresponding simulation results with different numbers of discrete elements. By comparing the simulation results with the analytical results, the relative error of our calculation model is obtained:

δ = (K_anal − K_simu)/K_anal

(29)

where δ is the relative error, K_anal is analytical results, MPa∙m^1/2, K_simu is simulation results, MPa∙m^1/2.

It can be seen from

Table 1. Relative errors of I-type and II-type stress intensity factors at the crack tip.

Element Number	$K_I/\left(MPa·m^{1/2}\right)$ (Analytical Results, 1.3293)		$K_{II}/\left(MPa·m^{1/2}\right)$ (Analytical Results, 0.7765)
	The Relative Error	$Simulation Results$	The Relative Error	$Simulation Results$
r	5.0%	1.2622	6.1%	0.7288
4	2.1%	1.3016	3.2%	0.7515
6	1.1%	1.3151	2.1%	0.7593
10	0.2%	1.3261	1.4%	0.7657

, as the number of discrete elements increases, the relative error decreases. The calculation accuracy can be guaranteed by increasing the number of discrete elements in numerical simulation. In the validation, a = 1 m, P_f = 1 MPa, θ = 60°.

3. Results and Discussion

The sequential fracturing model mentioned above is suitable for solving the mutual interference when the fractures are opened sequentially. In the following, we will talk about the fracture deformation and stress distribution in two kinds of sequential fracturing technology: one is about the staged fracturing technology in a single horizontal well, and the other is about sequential fracturing technology in a double horizontal well (zipper fracturing). The general range of values for parameters in numerical simulation can refer to our previous research [12]. The initial simulation parameters in this paper are as shown in

Table 2. Initial simulation parameters.

Young’s modulus (MPa)	19,830	Minimum horizontal stress (MPa)	−20
Poisson’s ratio	0.261	Fracture toughness (MPa·m1/2)	2.5
Maximum horizontal stress (MPa)	−24	Injection pressure (MPa)	−20
Injection angle (°)	90	Fracturing cluster number	3
Fracture spacing(m)	40	Perforating depth (m)	1

.

In order to clearly display the propping effect of proppant on the fracture surface, we simulated the variation law of HF aperture during sequential fracturing of three HFs in a well, and revealed the responding mechanism of induced stress under proppant action. At the same time, we simulated the trajectory morphology under fracture interference during zipper fracturing. The specific physical model is as shown in Figure 4. Inject the fracturing liquid into HF 1 first, stop injection when HF 1 reaches the target length and then begin the fracturing construction on HF 2. Complete the fracturing of three HFs in this order.

3.1. Deformation of Hydraulic Fracture and Induced Stress in a Single Well

The induced stress after hydraulic fracturing is mainly generated by the residual aperture of the fracture. Previous studies revealed that the induced compression stress around HF is always greater than anywhere else. Thus, in order to clearly reveal the interference mechanism from the previous hydraulic fracture completed by fracturing operation, we simulated the propagation of HF during staged fracturing first. Figure 5 shows the simulation results when hydraulic fracturing conduction is finished and three HFs reach equilibrium under in-situ stress.

From Figure 5, we can see that the deflection of hydraulic fractures in dynamic propagation may be affected by the HFs that have been fractured, and is away from the propped HF. Although the three fractures have the same injection pressure and are reached to balance under the same in-situ stress, induced stress produced by HF 1 is bigger than HF2, while HF 2 is bigger than HF 3. This is because, during the fracturing of HF 2 and 3, induced stress generated by propped fractures already exists, which leads to a gradual decrease in the maximum fracture aperture during subsequent hydraulic fracturing. The reduction of maximum HF aperture directly results in a reduction in proppant injection volume and residual aperture of HF.

In contrast, the fracturing HF may also affect the aperture of fractured HF. In order to quantitatively evaluate the aperture of HF 1 influenced by different fracturing stages, we conducted the following simulation, as shown in Figure 6.

Due to fracturing fluid filtration in HF, previously injected fracturing fluid cannot support HF to maintain initial maximum aperture. Comparing aperture of HF 1 under the situation of pump on and off of HF 1 in Figure 6a, we can see that, under the formation of compressive stress, the aperture of HF may have some retraction after hydraulic fracturing but will not be closed completely under proppant support. At the same time, if there is another HF that is under fracturing construction (HF 2 or 3), the induced stress generated by dynamic extension HF may compress HF 1 further (Figure 6a). However, from Figure 6b, we can see that, compared with HF 2 being pumped on, the aperture of HF 1 will have a weeny rebound after HF 2 is pumped off. The fracture residual width will decrease with an increasing number of stages, but the decrease is small and below 0.2 mm.

3.2. Sensitivity Analysis

We have analyzed the interaction mechanism of HFs in sequential fracturing above. However, the interaction degree is affected by many factors; therefore, a series of sensitivity analysis is carried out—the first is about the stiffness of proppant. We also quantitatively studied the residual aperture of HF under different stiffness of bracing member (Figure 7a). The induced stress caused by residual aperture of HF was shown in Figure 7b.

Simulation results showed that the support effect acting on the fracture surface is not linear. When the stiffness is less than 1 × 10¹⁰, the deformation of HF aperture is obvious. Conversely, when the stiffness is greater than 5 × 10¹⁰, the variation of residual HF aperture is not obvious. The greater the stiffness, the greater the influence of residual deformation on the in-situ stress after the fracture surface is closed.

In order to reveal the influence degree of this induced stress by residual aperture on subsequent HF deformation, based on the physical model as Figure 4a and Table 2, a numerical simulation study was carried out. Figure 8 is the simulation results about HF 2 with different stiffness and fracture spacing. To make the simulation results easy to compare, we did not display HF 1 temporarily.

From Figure 8a, we can see that the greater stiffness and the smaller fracturing spacing, the easier it is to cause subsequent fractures to be deflected. Bigger stiffness can keep larger residual aperture of HFs, and further produce bigger residual induced stress. This induced stress will deflect other fracture trajectories. However, the interference will decrease with the increase of fracture spacing (Figure 8b).

3.3. Deformation of HF and Distribution of Induced Stress under Zipper Fracturing

We have evaluated the influence of residual HF aperture in the above in staged fracturing when considering the reverse supporting force by proppant. Sequential fracturing not only exists in a single well, but also in multi-well fracturing. In this part, based on the sequential fracturing calculation model considering the supporting force of proppant in double wells, we simulated the propagation morphology of HF and characteristic of induced stress distribution in zipper fracturing. The physical model is as shown in Figure 4b, and the simulation results are shown in Figure 9. Well spacing is 70 m.

Simulation results in Figure 9 showed that, under the interference of adjacent fractures, hydraulic fractures tend to intersect among the two horizontal wells. In addition, at the same time, we should not ignore the existence of residual aperture of HF after fracturing. These residual apertures produced by proppant will exist for a long time and influence the propagation of subsequent HF. A series of sensitivity analyses were conducted to further investigate the influence mechanism of residual apertures on subsequent hydraulic fracture propagation.

From Figure 10, we can see that, with the decrease of well spacing and fracture spacing, HF 2 is more likely to deflect to HF 1. Setting a reasonable well spacing and fracturing spacing is beneficial for preventing the intersection of HF in adjacent wells. Compared to fracture spacing, the deflection of HF 2 is influenced by well spacing.

4. Conclusions

In this paper, a hydraulic fracture propagation model was established, which considered the influence of residual fracture aperture on subsequent hydraulic fracture propagation. The interaction of hydraulic fractures during sequential fracturing is studied by a series of numerical simulations:

(1)

As large quantities of proppant are injected into hydraulic fracture during hydraulic fracturing, fractures will not completely close after the fracturing operation is completed. This residual aperture caused by proppant may produce induced stress and change the distribution of in-situ stress. Induced stress by residual aperture gradually decreases with the increase of vertical distance from the fracture plane and the decrease of residual aperture of fracture.

(2)

The residual aperture will also influence the propagation and maximum aperture of subsequent fracture. When the fracture spacing in sequence fracturing is closer, the residual aperture will inhibit the opening degree of the subsequent fracture, which in turn affects the injection volume of the proppant. During sequence fracturing, fractures tend to exclusion and turn away in staged fracturing, on the contrary, which tends to approach and intersect in zipper fracturing.

(3)

Subsequent fracturing in turn compresses the previously cracked fracture, resulting in a further reduction in residual aperture, and after the fracture construction is completed, the previously pressurized fracture aperture is rebound. As the number of hydraulic fracture increases, the residual aperture of the previously pressed fracture gradually decreases. However, the fluctuation of fracture aperture mentioned above is small and less than 0.2 mm.

(4)

Sensitivity analysis shows that, in staged fracturing, the smaller the fracturing spacing, the more likely subsequent fractures are to be deflected, while in zipper fracturing, the effect of fracture spacing is not obvious. Well spacing can obviously influence the deflection of subsequent fracture in zipper fracture. With the increase of stiffness, the residual aperture of the hydraulic fracture increases, and the subsequent fractures are more likely to be deflected.

There are still some shortcomings in this paper. For example, the proppant stiffness is still fixed, which cannot reflect the nonlinear propping effect of proppant on fracture surfaces. At the same time, the fracturing fluid filtration process in the fracture is not considered. In the future, we will focus on optimizing these deficiencies and coupling them with existing models.