Numerical Simulation of Multi-Fracture Propagation Based on the Extended Finite Element Method

Abstract

Multi-stage, multi-cluster fracturing in horizontal wells is widely used as one of the most effective methods for unconventional reservoir transformation. This study is based on the extended finite element method and establishes a multi-hydraulic fracturing propagation model that couples rock damage, stress, and fluid flow, and the influence of horizontal stress difference and cluster spacing on fracture propagation is quantitatively analyzed. The simulation results show that changes in horizontal stress differences and inter-cluster spacing have a significant impact on the final propagation morphology of hydraulic fractures, and the change of the fracture initiation sequence forms different stress shadow areas, which in turn affects the propagation morphology of the fractures. When two fractures simultaneously propagate, they will eventually form a “repulsive” deviation, and a smaller stress difference and a decrease in inter-cluster spacing will lead to a more significant deviation of the fracture. Specifically, when the horizontal stress difference is 4 MPa and the cluster spacing is 6 m, the offset of the fracture tip along the direction of minimum horizontal principal stress is about 1.6 m, compared to the initial perforation position. When two fractures propagate sequentially, the fractures do not significantly deviate and propagate along the direction of maximum horizontal principal stress. When fractures propagate sequentially, the stress difference has little effect on the morphology of the fracture, but changes in inter-cluster spacing will significantly affect the length of the fracture. This study quantifies the effect of inter-fracture interference on fracture propagation morphology, providing guidance for optimizing the construction parameters of multi-stage hydraulic fracturing.

1. Introduction

Hydraulic fracturing is an important stimulation measure for the exploration and development of some unconventional reservoirs, such as shale reservoirs and tight sandstone reservoirs [1,2]. A single-plane fracture is not enough to realize the economic exploitation of unconventional oil and gas resources, and the complex fracture network is a necessary condition to improve the oil and gas production in low-permeability reservoirs, which is usually realized by multi-stage hydraulic fracturing of horizontal wells. In order to improve the operation efficiency, the usual construction method is to form several perforation clusters and create multiple hydraulic fractures in the same section of the reservoir through multi-stage hydraulic fracturing of horizontal wells, thereby increasing the volume of reservoir reconstruction [3,4,5]. However, according to some scholars’ statistical analysis of production data, less than half of the perforation clusters contribute to the oil and gas production within a section, which means that most perforation clusters do not achieve the expected effect of increasing production [6,7]. The real increase in oil and gas production comes from the combination of unconventional well structures and large-scale hydraulic fracturing. Single-stage fracturing has evolved into multi-stage production enhancement, and single-well-fracturing has evolved into simultaneous fracturing of branch wells to increase reservoir control and improve well productivity. This has become a commonly used measure to increase production in the industry today [8,9,10].

Therefore, there is an urgent need to improve the efficiency of production and reduce extraction costs; however, there are many factors that affect the propagation of hydraulic fractures: (1) The physical properties of the reservoir near the wellbore are heterogeneous, including natural faults, in situ stress, and rock yield strength. (2) Interactions between different hydraulic fractures within the same horizontal wellbore may inhibit further propagation of some fractures, especially if a fracture is subjected to additional stress from adjacent fractures. This phenomenon is known as the stress shadow effect, which can occur even in homogeneous reservoirs [11,12,13]. (3) The dynamic fluid distribution between hydraulic fractures is related to wellbore hydraulics and perforation characteristics, and is influenced by wellbore friction, perforation friction, and hydraulic fracture propagation control.

Regarding the required number of fractures or spacing between fractures, it is important to note that the optimal production should be achieved through a cost-effective and efficient fracture network. Roussel and Sharma have studied this issue from the perspectives of production and geomechanics [14]. A study was conducted on the distance between hydraulic fractures from a production perspective. Yu and Sepehrnoori [15] observed that fractures too close to each other did not significantly increase production. According to their research results, fracture spacing, length, and well spacing can be optimized under certain porosity, permeability, and fracture conductivity conditions. However, this study did not consider the impact of geological factors of the reservoir on determining the optimal distance between fractures. Many studies have shown that geological and geo-mechanical information about rocks and their variations in the reservoir is also important, as stress primarily controls the initiation and propagation of fractures [16,17]. Stress shadow mainly refers to the phenomenon of local high stress in the direction perpendicular to the fracture surface near the fracture center, which will lead to the reorientation of the direction of the maximum stress in the stress-affected area, and then make the subsequent fracture propagation deviate or even parallel to the wellbore axis. Therefore, it is necessary to optimize the fracture spacing to obtain the maximum number of fractures perpendicular to the wellbore [18].

Currently, there are many numerical methods used for hydraulic fracturing analysis, among which the finite element method (FEM) is the most widely used. The FEM is essentially a numerical method for the mechanics of continuous media. In order to use the finite element method for analyzing hydraulic fracturing and other discontinuous problems, improvements need to be made to the traditional finite element method. Methods for improving the finite element method for non-continuous problems can be divided into two categories: unfixed-mesh methods and fixed-mesh methods [19,20]. Physics-based hydraulic fracturing models usually include different coupled components [21,22]: rock deformation and fracture opening and closing, fluid flow within the fracture, and fracture initiation and propagation. Many scholars have already carried out comprehensive research on the propagation of fractures through the use of numerical simulation methods. Sendon [23] established an analytical model of induced stress around a single planar fracture, which became the cornerstone for subsequent research on fracture interference. Bunger and Peirce conducted research on the propagation of winged fractures and the interaction between them in isotropic media, proposing simple measures to promote synchronized fracture propagation [24]. Zhang et al. studied the deflection and propagation of hydraulic fractures encountering planar bedding planes using boundary element and finite difference methods. The results showed that the deflection of fractures and fluid invasion into weak interfaces mainly depend on the local stress and deformation state at the intersection [25]. Most of the existing studies on fracture interactions have focused on static descriptions, with less attention paid to changes in the magnitude and orientation of the inter-fracture stress during dynamic fracture propagation. Furthermore, there has been little research on the impact of the order of fracture initiation on competitive fracture propagation, and a lack of quantitative exploration of the stress shadow effect on fracture propagation.

In this study, we employ the extended finite element method (XFEM) to simulate hydraulic fracture propagation, which is a type of fixed-grid method. The key idea of XFEM is to enrich the finite element interpolation space with additional discontinuous functions, so that the displacement jumps associated with the fracture can be modeled within the element [26,27]. Different from the existing research, based on the bilinear T-S criterion, this study uses stiffness degradation to describe the damage evolution process of fractures, and realizes the coupling process of fluid flow and solid damage deformation to simulate the expansion of two clusters of different fractures in the same section of horizontal wells. Using the advantages of finite element software simulation, the dynamic change process of the stress shadow size and influence area during fracture propagation is intuitively characterized. The fracture length and width data are extracted by the plug-in program, and the influence of cluster spacing, the horizontal stress difference, and the fracture initiation sequence on fracture propagation morphology is emphatically explored.

2. Materials and Methods

2.1. Fluid–Solid Coupling Simulation

During hydraulic fracturing in saturated porous media, the fluid acting on the fracture surface causes deformation of the solid phase in the porous media. The change in rock pore pressure caused by the injection of fracturing fluid leads to changes in fracture morphology and permeability. Therefore, the process of fracture propagation in geological formations is a dynamic coupling process of viscous fluid flow and rock deformation [4,28,29]. In this study, we investigated the extension of hydraulic fractures and the deformation of the rock matrix during hydraulic fracturing. A homogeneous and isotropic 2D elastic medium was chosen as the region for hydraulic fracture propagation. The rock deformation was described by linear elasticity theory, and the equilibrium equation was satisfied by:

$$\left\{\begin{array}{l}{\sigma }_{ij,j}+f_i=0\hfill \\ {\epsilon }_{ij}=\left(u_{j,i}+u_{i,j}\right)/2\hfill \\ {\sigma }_{ij}=C_{ijkl}{\epsilon }_{kt}\hfill \end{array}\right.$$

(1)

where, ${\sigma }_{ij}$ is the stress tensor, $f_i$ is the volume force on the rock mass, ${\epsilon }_{ij}$ is the strain tensor, $u$ is the medium displacement, and $C$ is the elastic tensor.

To simulate fluid flow in hydraulic fractures, the porous media continuity equation was applied in weak form by imposing pore pressure at each node to simulate the flow of fracturing fluid within the porous medium:

$$\frac{\mathrm{d}}{\mathrm{d}t}\left({{\displaystyle \int }}_V{\rho }_w{\phi }_w\mathrm{d}V\right)+{{\displaystyle \int }}_S{\rho }_{\mathrm{w}}{\phi }_{w}n{q}_w\mathrm{d}S=0$$

(2)

where ${\rho }_w$ is the density of the fluid, kg/m³$, {\phi }_w$ is the matrix porosity, $q_w$ is the average flow velocity of the fluid in the matrix, m³/s, and $n$ is the outer normal direction of the surface S.

The flow equation in the above equation satisfies the Darcy equation:

$$q_w=-\frac{1}{n_{w}g{\rho }_w}{k}^{\prime }\cdot \left(\nabla p_w-{\rho }_{w}g\right)$$

(3)

In this study, the concept of the permeability coefficient was used to characterize the fluid permeability, which is expressed as follows:

$$k^{\prime }=\frac{k\rho g}{\mu }$$

(4)

where, $k^{\prime }$ is the permeability coefficient, m/s, k is the permeability, m², $\mu$ is the fluid viscosity, Pa·s, $n_w$ is the ratio of fluid volume to total volume, ∇P_w is the pressure gradient in the direction of fracturing fluid flow, and $g$ is the free-fall acceleration, m/s².

2.2. The Mathematical Model of Extended Finite Element

2.2.1. Extended Finite Element Method

In the traditional finite element method, the problem domain is discretized into a series of small elements, and the solution is approximated using polynomial functions within each element. This method performs well in dealing with simple geometries and continuous media. However, when there are local discontinuities such as fractures, it requires increasing the mesh refinement to capture these discontinuities, leading to increased computational costs. The extended finite element method is based on the classical finite element method and incorporates additional generalized degrees of freedom to represent local discontinuities. This allows simulations to be performed on relatively coarse grids without the need for mesh refinement around fractures [26,30].

The extended finite element method improves the interpolation shape function within the element. The method incorporates the use of the Heaviside step function, $H_{\left(x\right)}$, and the asymptotic fracture tip function, $F_{\alpha }\left(x\right)$, to account for the discontinuity of the fracture surfaces, where $H_{\left(x\right)}$ is used to characterize the intermittent displacement field:

$$H(x)=\mathrm{sgn}(\phi (x))=\left\{\begin{array}{l}1,\phi (x)>0\hfill \\ -1,\phi (x)<0\hfill \end{array}\right.$$

(5)

$F_{\alpha }\left(x\right)$ is used to characterize the singular displacement field and describe the singularity of the fracture tip stress. The expression in polar coordinates is:

$$F_{\alpha }(x)=\left[\sqrt{r}\sin \frac{\theta }{2},\sqrt{r}\cos \frac{\theta }{2},\sqrt{r}\sin \theta \sin \frac{\theta }{2}\right.\left.\sqrt{r}\sin \frac{\theta }{2}\cos \frac{\theta }{2}\right],\alpha =1,2,3,4$$

(6)

Therefore, its displacement vector function, u, which characterizes the overall division characteristics, is:

$$u=\sum _{I=1}^N{N}_I(x)\left[u_I+H(x)a_I+\sum _{\alpha =1}^4{F}_{\alpha }(x)b_I^{\alpha }\right]$$

(7)

where $N_1\left(x\right)$ is the ordinary nodal displacement form function, $u_I$ is the continuous part of the displacement solution, $a_I$ and $b_I^{\alpha }$ are the nodal extended degree of freedom vectors, $H(x)$ is the intermittent jump function of the fracture surface, and $F_{\alpha }(x)$ is the fracture tip stress asymptotic function [31,32]. The first part of the equation applies to all nodes within the model, the second part applies to nodes whose functional shape is formed by internal fractures, and the third part is used only for nodes whose shape function is separated by the fracture tip [33]. As shown in Figure 1, when an element is undamaged, each virtual node lies entirely on the primary node; when the element is cracked, the fracture element is divided into two parts, each of which consists of several virtual nodes and primary nodes dependent on the fracture direction. Each virtual node is unbound from its origin and can be moved [34].

The expressions for the interpolation of the pressure field and the displacement field have the same form:

$$p=\sum _{I=1}^N{N}_l(x)\left[p_I+\phi (x)p_I+\sum _{\alpha =1}^4{P}_{\alpha }(x)p_I^{\alpha }\right]$$

(8)

where, $p_I$, $p_I^{\alpha }$ are the node pressure expansion degrees of freedom, and $\phi (x)$ and $P_{\alpha }(x)$ are the enhancement functions of the pressure nodes.

2.2.2. Fracture Propagation and Extension Criteria Based on Extended Finite Element

In this study, the extended finite element method was used to simulate the expansion of hydraulic fractures. The tensile-separation criterion based on the damage mechanics of cohesive units was used for the initiation and extension of hydraulic fractures [35]. As shown in Figure 2, the first part involves determining the effective displacement at complete damage, ${\delta }_m^f$, relative to the effective displacement at the onset of damage, ${\delta }_m^0$, or the energy dissipation due to failure, $G_C$. The second part mainly defines the evolution of the damage parameter, $D$, between initial failure and complete failure. This parameter can be specified by a linear or exponential softening law, which is directly provided in a table that shows the relationship between effective displacement and damage. When using the linear displacement-based extension criterion, it can be expressed as:

$$D=\frac{{\delta }_{\mathrm{m}}^{\mathrm{f}}\left({\delta }_{\mathrm{m}}^{\max }-{\delta }_{\mathrm{m}}^0\right)}{{\delta }_{\mathrm{m}}^{\max }\left({\delta }_{\mathrm{m}}^{\mathrm{f}}-{\delta }_{\mathrm{m}}^0\right)}$$

(9)

where ${\delta }_m^{max}$ is the maximum displacement of the unit, ${\delta }_m^f$ is the displacement at which the unit is opened, and ${\delta }_m^0$ is the displacement at which the unit starts to be damaged.

The constitutive model of rock before damage is linear elasticity, and the element that degrades to complete failure after damage. In this study, the maximum principal stress criterion was used as the fracture initiation criterion, that is:

$$f=\frac{{\sigma }_{\max }}{{\sigma }_{\max c}}$$

(10)

where, ${\sigma }_{\max c}$ is the maximum critical stress to which the rock is subjected, MPa, and ${\sigma }_{\max }$ is the maximum principal stress to which the rock is subjected, MPa.

In this study, the type of damage evolution for hydraulic fracture propagation was selected as displacement, and mode-independent was used as the mixed-mode behavior. When the initial damage reached 0.001 mm, the element was considered completely broken. During the process of hydraulic fracture propagation, the formation’s permeability and porosity will change with the change of the effective stress of the porous medium in the formation, so the coupling relationship between the formation stress field and the seepage flow field must be considered. According to the principle of virtual work, the stress balance equation can be obtained as follows:

$${{\displaystyle \int }}_V\sigma \delta \epsilon dV={{\displaystyle \int }}_{S}t\delta vdS+{{\displaystyle \int }}_V\widehat{f}\delta vdV$$

(11)

where $\delta v$ is the virtual velocity, m/s, $\delta \epsilon$ is the virtual variation rate, s⁻¹, $\widehat{f}$ is the volume force per unit volume, N/m³, t is the external surface force per unit area, N/m³, and $\sigma$ is the total stress in the porous medium of the formation, Pa.

Based on the principle of mass conservation, the continuity equation of fluid media can be obtained as follows:

$${{\displaystyle \int }}_V\delta v\frac{1}{J}\frac{\mathrm{d}}{\mathrm{d}l}\left(J{\rho }_{\mathrm{w}}{n}_{\mathrm{w}}\right)\mathrm{d}V+{{\displaystyle \int }}_V\delta v\frac{\partial }{\partial x}\left({\rho }_{\mathrm{w}}{n}_{\mathrm{w}}{v}_{\mathrm{w}}\right)\mathrm{d}V=0$$

(12)

where, $J$ is the rate of change of the formation pore volume, $n_{\mathrm{w}}$ is the ratio of formation liquid volume to total volume, ${\rho }_{\mathrm{w}}$ is the density of liquid in the formation pore, kg/m³, and $v_{\mathrm{w}}$ is the flow velocity of the formation pore liquid, m/s.

2.3. Stress Field Superimposed on Fracture Propagation of Multiple Clusters in Horizontally Oriented Well

During the fracturing process of the horizontal well, the formation of horizontal wellbore will significantly affect the original in situ stress field. Therefore, the injection of fracturing fluid will cause the pore pressure around the wellbore to increase and generate additional stress due to the difference between the pore pressure of the wellbore and that of the formation. In situ stress is the natural stress that exists in a rock mass, and it is regularly distributed in a three-dimensional state. There are many factors that affect the change of in situ stress, among which the obvious ones are wellbore pressure, the original tectonic stress component, and the change caused by fracturing fluid seepage [36]. For the process of propagation of multi-cluster fractures in the horizontal well, besides considering the stress disturbance caused by the propagation of different fractures, the significant temperature difference of the reservoir at different well depths and the thermal effect of the drilling fluid and fracturing fluid injected on the wellbore and fractures also bring additional stress effects. Figure 3 shows the multi-fracture stress superposition model, which regards the rock as an infinitesimal deformation, porous elastic body, according to the principle of stress superposition, the disturbed stress field of horizontal well-fracturing is the sum of the initial in situ stress, wellbore pressure, fracturing fluid seepage, fracture-induced stress, and additional stress caused by thermal effects, that is:

$${\sigma }_{\mathrm{sum}}={\sigma }_{\mathrm{ini}}+{\sigma }_{\mathrm{w}}+{\sigma }_1+{\sigma }_{\mathrm{f}}+{\sigma }_T$$

(13)

where, ${\sigma }_{\mathrm{sum}}$ is the total fracture disturbance stress, ${\sigma }_{\mathrm{ini}}$ is the initial ground stress, ${\sigma }_{\mathrm{w}}$ is the stress around the wellbore after drilling, ${\sigma }_1$ is the fracture fluid percolation stress, ${\sigma }_{\mathrm{f}}$ is the fracture-induced stress, and ${\sigma }_T$ is the thermal stress.

3. Model Construction and Validation

3.1. Model Construction

In the simulation of simultaneous expansion of multiple clusters of fractures in horizontal wells, a total of three fractures with equal spacing were set in the horizontal wells, where the direction of the horizontal wellbore was the direction of the minimum horizontal principal stress to ensure that the formed hydraulic fractures were transverse fractures. The size of the numerical simulation model was 100 m × 100 m, and the related geo-mechanical parameters and fracturing construction parameters are shown in

Table 1. Mechanical parameters of stratigraphic rocks.

Variables	Unit	Value of Reservoir Rock
Young’s modulus	GPa	15
Poisson’s ratio	/	0.25
Permeability coefficient	m/s	1 × 10−7
Initial pore ratio	/	0.1
Filtration loss factor	m/(Pa·s)	1 × 10−14
Tensile strength	MPa	6

.

The simulation adopted the super-hydrostatic pressure system, with a fracturing fluid flow rate and viscosity of 0.12 m³/min and 1 mPa·s, respectively. The perforation length was 1 m, and the pumping time of each hydraulic fracture was maintained for 100 s. To simplify the calculation, the horizontal wellbore was located at the bottom of the model, and the perforation direction was aligned with the maximum horizontal principal stress direction. The established geometric model is shown in Figure 4. A constant-pressure water fracturing system was applied, and the boundary pore pressure was set to 0 MPa and kept constant during the fracturing process. The reservoir rock matrix was represented by a structured, standardized, linear quadrilateral plane strain element type known as CPE4P, and the hydraulic fracture propagation was simulated with two-node linear truss elements (T2D2).

In the XFEM, fracture is modeled independent of mesh configuration and element type. This means that no remeshing is required, and discontinuity, including fractures, need not be aligned with element boundaries [37]. However, different mesh configurations have an impact on the simulation convergence and effectiveness. Therefore, we conducted an investigation on various mesh sizes and selected four different mesh sizes for evaluation: 0.3 m, 0.5 m, 0.8 m, and 1 m, respectively. The relevant mechanical parameters are referenced from Table 1. The cluster spacing was 6 m, the horizontal differential stress was 4 MPa, and the fracturing fluid displacement and viscosity were 0.12 m³/min and 1 mPa·s, respectively. The length of the perforation was 1 m, and the injection time was set to 100 s.

Figure 5 displays the morphology of fractures after simultaneous initiation and propagation of dual-cluster fractures extracted using the plugin under different mesh sizes. As observed from Figure 5, it can be inferred that changing the mesh size had minimal impact on the final morphology of the fractures after propagation. Two cases with the largest difference in half-fracture height were selected for comparison. When the mesh size was 0.3 m, the half-fracture height on the right side was 15.37 m, while with a mesh size of 0.8 m, the half-fracture height was 14.91 m. The difference between the two cases was only 0.46 m. Considering that smaller mesh sizes are suitable for smaller initial perforation lengths, it is important to balance the computational efficiency and simulation accuracy. Finer mesh configurations tend to significantly increase the computation time and can lead to convergence issues [38]. Taking into account both the computational time and the simulation effectiveness, a mesh size of 0.5 m was selected as the optimal choice for the fracture propagation numerical model in this study. This mesh size strikes a balance between capturing essential details and maintaining an acceptable computational efficiency, ensuring a reasonable compromise between accuracy and computational resources.

3.2. Model Verification

The interaction among multiple fractures plays a crucial role in determining the geometric shape of the fractures during their propagation. In this study, we conducted a comparative analysis of the dual-cluster fracture propagation model for horizontal wells proposed by Wu et al. [39]. The wellbore and initial perforation locations in their model were similar to those depicted in Figure 4. The horizontal wellbore was aligned parallel to the minimum horizontal principal stress direction, while the initial perforation locations were positioned along the maximum horizontal principal stress direction. The model input parameters are presented in

Table 2. Input parameters (after Wu et al., 2012 [39]).

Variables	Unit	Value
Young’s modulus	psi	4.35 × 106
Poisson’s ratio	/	0.35
Maximum horizontal stress	psi	6903
Minimum horizontal stress	psi	6773
Injection rate	bbl/min	40
Fluid viscosity	cp	1
Distance between initiation points	ft	33

, and the simulation results are shown in Figure 6. Figure 6a illustrates the fracture opening after expansion in the current model, while Figure 6b displays a comparison of the fracture propagation paths between the two models. In a horizontal well, parallel fractures may deviate from each other due to the stress shadow effect. Once the fractures open, they exert additional stresses on the surrounding rock and neighboring fractures. This can lead to local variations in the direction of the horizontal principal stress and deviations in the fracture paths from the planar geometry.

From Figure 6, it can be observed that the fracture propagation paths of the two models exhibited a high degree of consistency, with only slight deviations at the fracture tips. The analysis suggests that these deviations might be attributed to the different boundary effects caused by the variations in model sizes. As the fractures propagated, they encountered different resistances at the fracture tips, leading to slight differences in the final morphology of the fractures. Based on the above analysis, it can be concluded that the simulation results of the two models mentioned earlier exhibited good agreement. This confirmed the effectiveness and accuracy of the model used in this study.

4. Numerical Simulation Results

To investigate the competition of fracture propagation under simultaneous initiation of two fractures, two perforated intervals with a perforation depth of 1 m were placed horizontally along the direction of the minimum horizontal principal stress in the study area as the initial fractures. The injection time was set to 100 s and the injection rate was chosen as 0.002 m/s. The two perforated intervals concurrently injected fluid. This paper also explores the interference effect of the stress disturbance caused by the propagation of the first fracture on the subsequent initiation fracture under the condition of two fractures initiating at different times. In order to simulate the stress shadow effect induced by fracture propagation under reservoir conditions, two fractures were initiated at different times with a constant injection rate of 0.002 m/s and a duration of 100 s each, followed by a 600 s pressure depletion after the first fracture was completed. The simulation of simultaneous and sequential fracture initiation revealed how stress field and fracture spacing influenced the complexity and diversity of fracture propagation patterns. To investigate the effects of stress difference and fracture spacing on fracture propagation, we performed simulations with different values of these parameters. We varied the stress difference from 2 to 8 MPa and the fracture spacing from 6 to 12 m. Here, we present the results and analysis of how these factors influenced fracture propagation.

4.1. Horizontal Differential Principal Stress

In the simulations under different reservoir stress conditions, we set different values of minimum and minimum horizontal principal stress, such that the horizontal principal stress difference was 2, 4, 6, or 8 MPa. The overburden stress was set to 15 MPa.

Table 3. Reservoir ground stress conditions and cluster spacing.

Case	Vertical Stress/ MPa	Minimum Horizontal Principal Stress/MPa	Maximum Horizontal Principal Stress/MPa	Horizontal Principal Stress Difference/MPa	Cluster Spacing/m
Case 1	15	6	8	2	10
Case 2	15	6	10	4	10
Case 3	15	6	12	6	10
Case 4	15	6	14	8	10

shows the reservoir rock mechanical parameters for this area. The other basic rock mechanical parameters are presented in Table 1. The perforation cluster spacing was kept at 10 m. As shown in Figure 7, the distribution of reservoir pore pressure underwent changes as the stress difference increased. Furthermore, the deviation angles of the fracture morphology gradually decreased as a result, and stress concentration phenomena appeared at the tips of the fractures. By programming and extracting the nodal information of the two-dimensional fractures, the post-propagation coordinate information of the fractures can be obtained. Based on Figure 8, the morphology of the fractures exhibited consistency under different stress difference conditions when the two fractures simultaneously initiated and propagated. This consistency is attributed to the initial perforation cluster being symmetrically located in the model and having consistent boundary conditions. Under a stress difference of 2 MPa, the deviation of the fracture tip was significant, deflecting approximately 1.02 m along the direction of the minimum horizontal principal stress, as compared to the initial initiation point. As the stress difference increased, the variation in the fracture morphology notably decreased. Under the stress difference conditions of 6 MPa and 8 MPa, the fracture morphology was approximately consistent.

In contrast, for simulation calculations involving the sequential initiation of two fractures, the interference from the stress shadow effect generated by the first fracture after expansion may affect the propagation of the second fracture. Thus, the morphology of the fractures will exhibit significant differences compared to the case where both fractures simultaneously initiated. As illustrated in Figure 9, the change in the horizontal principal stress difference had little influence on the final morphology of the fracture. However, due to the stress interference from the previously initiated fracture on the subsequently initiated fracture, the average half-length of the later-initiated fracture decreased by approximately 6 m under the four different stress difference conditions. In contrast, the length of the previously initiated fracture did not decrease compared to the simultaneous initiation condition. These findings suggest that stress interference between fractures can have a detrimental impact on reservoir fracturing and enhanced oil recovery by reducing the hydraulic fracturing effective area.

Figure 10 shows the variation of the minimum horizontal principal stress vector under the condition of a stress difference of 2 MPa. Due to the default setting of the simulator, the vector represents compressive stress as a negative value. From the figure, it can be observed that after the first hydraulic fracture extended, the surrounding stress field was disturbed, resulting in an elliptical-shaped stress shadow area that restricted the propagation of the second fracture. Figure 11 shows the distribution of the minimum horizontal principal stress at the ends of the two sequentially propagating fractures under the four different stress conditions. It can be observed from the figure that the minimum horizontal principal stress at the ends of the initial fracture reached its maximum value, resulting in a high-stress area at the tip of the fracture. When the second fracture propagated, it eventually stopped due to the fact that the pore pressure inside the fracture was lower than the fracture pressure of the element.

The subsequent fracture initiation pressure was also affected by stress interference. Under the condition of staggered fracture initiation, the difference in the pore pressure of the fracture elements was not significant among the four stress difference conditions. Taking the stress difference of 2 MPa as an example, Figure 12 shows the pore pressure changes with time at the perforation locations of the two fractures under the stress difference of 2 MPa. It can be seen from the figure that as time increased, the pore pressure of the two fractures first increased, and then decreased to a certain level and remained stable. The first fracture started to propagate in 0–100 s, followed by the pressure-relief stage of the first fracture from 100 to 700 s, and then the second fracture started to propagate. The initiation pressure of the first fracture was around 36.31 MPa, while the subsequent fracture initiation pressure increased to 43.08 MPa. This suggests that the propagation of the first fracture hindered the propagation of the second fracture, and ultimately, the pore pressure in both fractures tended to be consistent.

4.2. Perforation Cluster Spacing

The variation of cluster spacing also had a significant impact on the morphology of the fracture propagation. In this study, the effect of different perforation cluster spacing on fracture propagation was investigated. Cluster spacings of 6 m, 8 m, 10 m, and 12 m were set, respectively (

Table 4. Mechanical parameters of reservoir rocks and cluster spacing.

Case	Vertical Stress /MPa	Minimum Horizontal Principal Stress /MPa	Maximum Horizontal Principal Stress /MPa	Cluster Spacing /m
Case 5	15	6	10	6
Case 6	15	6	10	8
Case 7	15	6	10	10
Case 8	15	6	10	12

). The mechanical parameters of the reservoir rock are still based on Table 1, with the minimum horizontal principal stress set to 6 MPa, the maximum horizontal principal stress set to 10 MPa, and the overlying stress set to 15 MPa.

Figure 13 shows the distribution of fracture pore pressure after simultaneous initiation and propagation of fractures under different cluster-spacing conditions. It can be seen from the figure that the maximum pore pressure was generated at the fracture tip. Combined with the fracture morphology in Figure 14, it can be seen that as the cluster spacing decreased, the inter-fracture deflection angle gradually increased, and the end of the hydraulic fracture with a cluster spacing of 6 m deviated from the initial perforation cluster position by about 1.6 m. It should be noted that, in contrast to the simulation results under the condition of the changing stress difference (Figure 8), the two simultaneously initiated fractures under the conditions of 8 m and 12 m cluster spacing showed obvious asymmetry. The reason for this is that, under the distance condition of the shot-hole cluster, the initial shot-hole positions of the two fractures were not in the middle of the model, and the influence of the boundary effect led to greater resistance to the propagation of the second fracture.

Similarly, simulations were conducted to study the sequential initiation and propagation of fractures under different cluster-spacing conditions. The obtained fracture patterns are shown in Figure 15, which exhibited significant differences in the propagation of the two fractures. The later-initiated fracture was subjected to interference from the first fracture, resulting in a significant decrease in its half-length, but with a smaller deflection. This indicates that fracture deflection is more susceptible to interference when multiple fractures are competing for propagation at the same time. Moreover, the length of the fractures significantly varied during sequential propagation, and the interference became more pronounced as the cluster spacing decreased.

When the cluster spacing was 6 m, the half-length of the second fracture was reduced by about 11.5 m. By referring to the minimum horizontal principal stress vector map of the fracture propagation process under the cluster spacing of 6 m, as shown in Figure 16, it can be inferred that the propagation of the first fracture resulted in significant stress interference, which changed the magnitude and direction of the minimum horizontal principal stress and created stress shadow effects around the fracture, where the minimum horizontal principal stress around the fracture was higher than the original minimum horizontal principal stress, resulting in the stress shadow effect around the fracture. When the second fracture propagated, the same pumping displacement generated the same pore pressure at the fracture tip. However, due to the presence of high-stress zones, the fracture propagation was hindered, and the half-length of the fracture significantly decreased. In addition, for the case of sequential initiation, in the simulation examples of cases 5–8, the average initiation pressure of the second fracture was lower than that of the first fracture, with a difference of about 8 MPa. The trend of the relationship between the initial perforation pore pressure and time was similar to that shown in Figure 12.

5. Conclusions

In this study, the bilinear T-S criterion and stiffness degradation were used to describe the damage evolution of fractures, and a coupled model of fluid flow and solid damage deformation was constructed. Based on the extended finite element model, the propagation simulation of multiple cluster fractures during the horizontal well-fracturing transformation process was realized, and the following main conclusions were obtained.

(1) During the process of multi-cluster fracture propagation and extension in horizontal wells, the stress difference and the size change of inter-cluster clusters had a greater impact on the final morphology of fractures and the area of reservoir transformation. The fracture initiation sequence also directly affected the final fracture morphology. When the two fractures simultaneously initiated, a more obvious fracture deflection occurred; when the fractures initiated at different times, the fracture length of the later-initiated fracture was significantly reduced due to the stress interference caused by the earlier-propagated fracture.

(2) When simultaneously initiating, changing the stress difference and inter-cluster spacing formed two relatively symmetrical deflected fractures. The smaller the minimum horizontal principal stress difference and the inter-cluster spacing, the more obvious the fracture deflection. When the horizontal stress difference was 4 MPa and the inter-cluster spacing was 6 m, the fracture tip deviated from the initial perforation position by about 1.6 m along the direction of the minimum horizontal principal stress.

(3) When the two fractures separately initiated, the fractures basically extended along the direction of the maximum horizontal principal stress. Due to the propagation of the first fracture, the size and direction of the surrounding stress field were changed, and a high-stress zone was induced around the fracture, forming an elliptical stress shadow area, which hindered the propagation of subsequent fractures. The variation in the stress difference had a minimal impact on the fracture length. Under four different stress difference conditions, the length of the subsequent fracture decreased by approximately 6 m. On the other hand, the change in cluster spacing had a significant influence on the fracture length. When the cluster spacing was 6 m, the half-length of the fracture decreased the most. The subsequently expanding fracture had a reduction of approximately 11.5 m compared to the initially expanding fracture.

(4) The pore pressure at the perforation location also experienced significant changes with variations in the fracturing sequence. When the two fractures were simultaneously initiated, the fracture initiation pressure for both fractures was around 36 MPa. However, in the case of sequential fracturing, the initiation pressure for the subsequent fracture increased by 7 MPa, reaching 43.08 MPa. This is because the expansion of the first fracture altered the original stress field, resulting in an increased resistance for subsequent fracture initiation. Eventually, the pore pressure at the perforation locations of both fractures stabilized at approximately 13 MPa.