Simulation of Key Influencing Factors of Hydraulic Fracturing Fracture Propagation in a Shale Reservoir Based on the Displacement Discontinuity Method (DDM)

Abstract

In the process of the large-scale hydraulic fracturing of a shale gas field in the Weiyuan area of Sichuan province, the quantitative description and evaluation of hydraulic fracture expansion morphology and the three-dimensional distribution law are the key points of evaluation of block fracturing transformation effect. Many scholars have used the finite element method, discrete element method, grid-free method and other numerical simulation methods to quantitatively characterize hydraulic fractures, but there are often the problems that the indoor physical simulation results are much different from the actual results and the accuracy of most quantitative studies is poor. Considering rock mechanics parameters and based on the displacement discontinuity method (DDM), a single-stage multi-cluster fracture propagation model of horizontal well was established. The effects of Young’s modulus, Poisson’s ratio, the in situ stress difference, the approximation angle, the perforation cluster number and the perforation spacing on the formation of complex fracture networks and on the geometrical parameters of hydraulic fractures were simulated. The research results can provide theoretical reference and practical guidance for the optimization of large-scale fracturing parameters and the quantitative post-fracturing evaluation of horizontal wells in unconventional reservoirs such as shale gas reservoirs.

1. Introduction

The exploration and development of shale gas in China is still in the early stage of rapid development. The shale gas reservoir structure is complex. However, as a common unconventional reservoir development technology, hydraulic fracturing has been widely used in Weiyuan, Changning and other areas [1]. The morphology and geometry parameters of hydraulic fracturing are mainly affected by natural fractures, rock mechanic parameters of the reservoir, perforation spacing, the number of perforation clusters and other operation parameters [1,2,3]. When natural fractures develop, the probability of communication between hydraulic fractures and natural fractures increases, and reservoirs have the potential to form complex fracture networks. The larger the Young’s modulus, the less likely the rock is to undergo strain under the same stress, and the simpler the crack morphology formed [4,5]. The heterogeneity of the reservoir and the mutual interference between fractures are the root causes of uneven fracture propagation, and the number of perforation clusters and the perforation spacing have important effects on the distribution of the induced stress field [1,6]. In the process of multi-cluster fracturing in horizontal wells, the negative effects of the uneven expansion of fractures include the inability to achieve optimal reservoir transformation and the expansion of advantageous fractures, which leads to serious inter-well interference caused by communication with adjacent wellbore or hydraulic fractures [7,8,9]. In addition, due to the different properties and filling methods of different natural fracture cements, the cementation strength of natural fractures is different, so the cementation strength of natural fractures also has a great influence on the opening of fractures. As an important factor, the approaching angle also affects the opening of natural cracks.

When extracting unconventional oil and gas resources by hydraulic fracturing, the numerical simulation method has become the best method to study the fracture growth pattern in complex media. With the improvement in research theory and computing power, the numerical simulation of hydraulic fracturing has changed from the early simple two-dimensional model to the complex quasi-three-dimensional model. In recent decades, a series of numerical simulation calculation methods such as finite element method, boundary element method, discrete element method, meshless method and finite difference method have appeared [10]. Cheng et al. used the finite element method to deal with the elastic and flow equations, studied the problem of multi-layer hydraulic crack propagation in a plane homogeneous longitudinal direction, and proposed a finite element solution for plane arbitrarily shaped cracks perpendicular to the interface of multi-layer media [11,12,13]. However, this method has limitations in simulating the crack intersection and non-planar propagation. By establishing a finite element model of the dynamic competitive propagation of multiple fractures, Kresse et al. proposed an optimization method of perforation cluster spacing to ensure the balanced propagation of fractures by strengthening stress interference between fractures [14]. Zou et al. applied the boundary element method to discretized wellbore wall and fracture trace, solving the problem of the initiation and propagation of hydraulic fractures [15]. The model can simulate the non-plane propagation of hydraulic fractures on a horizontal plane, which is of great significance for the study of the hydraulic fracture propagation mechanism. Above, there are many studies on multi-layer hydraulic fracture propagation, non-plane hydraulic fracture propagation on the horizontal plane, induced stress field and geostress field and other factors influencing fracture propagation, fluid–structure coupling, and dynamic fracture network simulation of horizontal wells. However, there are two defects in the sensitivity analysis of these studies on related influencing factors: ① Indoor experiments are limited by indoor conditions and core size, resulting in significant deviations between physical simulation results and actual conditions; ② Most studies only analyze the corresponding qualitative trends or have poor accuracy in quantitative characterization.

In this work, a fracturing model is established using the key core technology research project software Frsmart 2.0 platform of CNPC. The fracture model simulates fracture propagation, fluid flow and proppant migration in a fully coupled fluid mechanics manner. By considering the parameters of rock mechanics and based on the displacement discontinuity method, a single-segment multi-cluster fracture propagation model for horizontal wells is established. The effects of Young’s modulus, Poisson’s ratio, in situ stress difference, approximation angle, perforation cluster number and perforation spacing on the formation of complex fracture networks and on the geometrical parameters of hydraulic fractures are simulated. Contrary to traditional numerical simulation methods such as finite element, the DDM method is used in this work, and the calculation accuracy was higher in the quantitative analysis of the sensitivity of related influencing factors. Its main principles include the following: ① the dimension of the calculation model is reduced by one dimension, so the calculation amount is reduced. The number of grids is reduced, and the input data preparation is simple. Under the same conditions, boundary elements tend to be more accurate than finite element and finite difference methods. ② DDM unit adopts a 1/2 displacement mode to simulate the singular stress field of the crack tip well, and the integral of its coefficient is precisely obtained, so it has high precision. This study provides theoretical reference and practical guidance for the optimization of large-scale fracturing parameters and quantitative post-fracturing evaluation of horizontal wells in unconventional reservoirs such as shale gas.

2. The Establishment of Hydraulic Fracturing Model

In this work, a simulation of hydraulic fracture propagation based on DDM was established to solve the problem of a strange stress field at the fracture tip, and sensitivity quantification characterization was, respectively, conducted on the reservoir and construction factors that affect the expansion of multiple fractures. The high-precision simulation results provide theoretical reference and practical guidance for the optimization of large-scale fracturing parameters and quantitative evaluation of hydraulic fracturing in shale reservoirs. The fracture simulation adopted in this work is a non-planar 3D fracture model, which follows the following assumptions:

(1)

All hydraulic fractures are vertical (vertical stress equal to the maximum principal stress or at least greater than the minimum horizontal principal stress).

(2)

Hydraulic cracks can turn in the horizontal plane.

(3)

There are parallel formations within each fracturing stage.

(4)

The fracture opening follows a linear elastic fracture mechanics theory.

In the process of the initiation and expansion of hydraulic fractures, the fracture stress interferes, and the fracture growth path will be distorted. If hydraulic fractures meet natural fractures in the process of expansion, the hydraulic fractures may directly penetrate the natural fractures or be captured by natural fractures. Fracture stress interference and fracture propagation intersection behavior should meet the corresponding criteria. Under the premise of meeting the above assumptions, displacement discontinuity method (DDM) is used to describe the fracture initiation and propagation intersection behavior of a hydraulic fracture, and the fracture stress interference and fracture propagation intersection behavior criteria are restricted. It is used to simulate the formation process of fracture network and the sensitivity analysis of each related variable.

The basic equations of DDM and the basic criteria of crack propagation are as follows (the flow chart as Figure 1).

2.1. Displacement Discontinuity Method

The relative displacement between two surfaces on a fracture surface is called displacement discontinuity, and its magnitude is the displacement discontinuity quantity. The displacement discontinuity method takes displacement discontinuity as the basic variable and deduces the stress value and displacement value of any point calculated based on displacement discontinuity, which is called the basic solution [5,14].

The deformation of a crack is described by the three-dimensional displacement discontinuity method:

$$\begin{array}{c}-{\sigma }_{sl}^i={\displaystyle {\displaystyle \sum }_{K=1}^N}{K}_{sl,sl}^{ik}{D}_{sl}^k+{\displaystyle {\displaystyle \sum }_{K=1}^N}{K}_{sl,sl}^{ik}{D}_{sh}^k+{\displaystyle {\displaystyle \sum }_{K=1}^N}{K}_{sl,nn}^{ik}{D}_{nn}^k\\ -{\sigma }_{sh}^i={\displaystyle {\displaystyle \sum }_{K=1}^N}{K}_{sh,sl}^{ik}{D}_{sl}^k+{\displaystyle {\displaystyle \sum }_{K=1}^N}{K}_{sh,sh}^{ik}{D}_{sh}^k+{\displaystyle {\displaystyle \sum }_{K=1}^N}{K}_{sh,nn}^{ik}{D}_{nn}^k\\ p-{\sigma }_{nn}^i={\displaystyle {\displaystyle \sum }_{K=1}^N}{K}_{nn,sl}^{ik}{D}_{sl}^k+{\displaystyle {\displaystyle \sum }_{K=1}^N}{K}_{nn,sh}^{ik}{D}_{sh}^k+{\displaystyle {\displaystyle \sum }_{K=1}^N}{K}_{nn,nn}^{ik}{D}_{nn}^k\end{array}$$

(1)

where D_nn—normal displacement discontinuity surface; D_sl, D_sh—a surface with discontinuous shear displacement in the direction of crack length and crack height; P—fluid pressure; σ_nn, σ_sl, σ_sh—effective formation normal and shear stresses, which can be obtained from the distant stress field by the superposition principle. The nine matrices Ks are boundary influence matrices related to the element size, position, orientation, and elastic constant. Superscripts i and k—element index; N—the total number of elements.

Fluid flow and proppant migration in fractures are controlled by the volume conservation equation:

$$\begin{array}{cc}\frac{\partial w}{\partial t}=& 𝛻\cdot \left[w(1-c)v^f+wc{v}^p\right]+\delta (x,y)q_0-q_l\\ & \frac{\partial (cw)}{\partial t}=𝛻\cdot \left(cw{v}^p\right)+\delta (x,y)c_0{q}_0\end{array}$$

(2)

where T—time; w—crack width, D_nn—normal displacement discontinuity surface; c—volume concentration of the proppant; q_o—the injection fluid rate; C₀—the specified volume concentration of proppant in the fluid injection; δ(x, y)—the Dirac function; and q_l corresponds to the frac fluid filtration velocity of the one-dimensional Carter filtration model:

$$q_l=\frac{2C_L}{\sqrt{t-\tau \left(x,y\right)}}$$

(3)

where τ(x, y)—the time the surface position is first exposed to fracturing fluid; C_L—Carter’s filtration coefficient.

The average frac fluid volume velocity (V^f), average proppant volume velocity (V^p), and proppant settling velocity (V^s) are, respectively, expressed as

$$V^f=\frac{w^2}{12{\mu }^f}{Q}^f\left(c,\frac{w}{a}\right)\left(𝛻p-{\rho }^{f}g\right)$$

(4)

$$V^p=Q^p\left(c,\frac{w}{a}\right)\left(v^f+v^s\right)$$

(5)

$$V^s=\frac{a^2}{12{\mu }^f}\left({\rho }^p+{\rho }^f\right)g$$

(6)

where p—fluid pressure, g—gravity vector, a—proppant radius, μ^f—apparent viscosity of the fracturing fluid, ρ^f—fracturing fluid density, v^s—proppant settling velocity, and ρ^p—proppant density.

The following function is used in a planar three-dimensional fracture to reflect the effect of proppant concentration on the viscosity of the fracturing fluid. This function can be used to describe the transition from the Poisson flow state to the Darcy flow state when the proppant concentration reaches its maximum value.

$$Q^f(c,\frac{w}{a})={\left(1-\frac{c}{c_{\max }}\right)}^{\beta }+\frac{a^2}{w^2}\frac{c}{c_{\max }}\stackrel{-}{D}$$

(7)

where β = 1.5, $\overline{D}=8{(1-c\alpha c_{\max }}_{\max }$), α = 4.1, and $c_{\max }$ = 0.585.

When the sand plug or proppant concentration reaches the maximum, the function value is 0; otherwise it is 1. The proppant function expression is different; all this can be obtained from the experimental data.

The equation of perforating pressure drop:

$$P_r=\frac{8\rho }{{\pi }^2{C}_d^2{d}_p}\left(\frac{q_{ri}}{n}\right)$$

(8)

where C_d is the dimensionless flow coefficient; d_p is the perforation diameter; q_ri is the volumetric injection rate at a specific injection location; n is the number of holes at a specific perforation; ρ is the mud density, and $\rho =\left(1-c\right) {\rho }^f+c{\rho }^p$. It was assumed that the pressure at the wellbore perforation is equal everywhere.

2.2. The Fracture Propagation Length under Load Increment

Through the analysis of energy angle, we found that the structure will release energy continuously during the process of crack propagation, and only when the energy is greater than the energy required to maintain the crack propagation can the crack continue to expand. In most cases, after the same load increment is applied, the length of the crack expansion will be different due to the different energy released by the crack expansion. Since the equivalent stress intensity factor K_e is proportional to the energy release rate G_e, the following steps can be used to calculate the crack growth length after each load increment is applied [16,17]:

(1)

Apply an additional incremental load on the basis of the existing load, and finally obtain the equivalent stress intensity factor K_I and K_II by calculating the stress intensity factor K_e at the crack tip under the current total load.

(2)

Conduct judgment on crack propagation. If the calculated equivalent stress intensity factor K_e is less than the critical stress intensity factor K_c, then the crack will not expand under the current load increment and need to go back to step (1). If K_e is greater than or equal to K_c, then the crack begins to expand. Based on the calculated crack expansion angle θ₁, the original crack extends a short length counterclockwise along the θ₁ direction.

(3)

Under the condition that the current total load remains unchanged, the stress field after the fracture extension for a short length is calculated, and the equivalent stress intensity factor K_e at the new crack tip is obtained.

(4)

Conduct judgment on crack propagation. If the equivalent stress intensity factor K_e is greater than or equal to the critical stress intensity factor K_c, the crack will continue to expand. According to the calculated new crack propagation angle θ₂, rotate θ₂ counterclockwise along the new crack tip of the previous section and extend it a further distance, as shown in Figure 2, then return to step (3). If K_e is less than K_c, then under the action of this incremental load, the crack has already expanded, so return to step (1).

2.3. Fracture Intersection and Expansion Criteria

The stress field distribution when the hydraulic fracture extends to the intersection with the natural fracture is shown in Figure 3. $\beta$—the approaching angle between the hydraulic fracture and the natural fracture, ${\tau }_{\beta }$ and ${\sigma }_{\beta ,n}$—the shear stress and normal stress acting on the surface of the natural fracture, and ${\sigma }_H$ and ${\sigma }_h$—the maximum and minimum horizontal principal stress.

When the maximum principal stress is greater than the tensile strength of the rock, the hydraulic fracture will expand forward through the natural fracture under the induction of the maximum horizontal principal stress, and the critical conditions are

$${\sigma }_1=T_0$$

(9)

$$\left|{\tau }_{\beta }\right|<C-{\mu }_{\mathrm{f}}{\sigma }_{\beta ,n}=0$$

(10)

where C is the natural crack interface bonding force, and ${\mu }_f$ is the natural crack friction factor.

The maximum circumferential stress criterion is used to judge the expansion behavior of cracks. The mutual integral method is used to calculate the stress intensity factor K_I and K_II at the crack tip and obtain the equivalent stress intensity factor $K_e$. By comparing the equivalent stress intensity factor $K_e$ in the fracture direction with the fracture toughness of the rock mass $K_{IC}$, the crack expansion can be judged. The critical conditions are

$$\left\{\begin{array}{c}{K}_e=\cos \frac{\theta }{2}\left(K_{\mathrm{I}}{\cos }^2\frac{\theta }{2}-\frac{3K_{\mathrm{II}}}{2}\sin \theta \right)\\ K_e\ge K_{IC}\end{array}\right.$$

(11)

The formula for calculating the crack initiation angle is

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

(12)

2.4. Fracture Propagation-Induced Stress Field Superposition Criterion

When calculating fracture propagation by the displacement discontinuity method, it is necessary to divide the fracture into multiple line elements. The total displacement and stress of the operation object are calculated by superposition of all discontinuous elements. The influence of discontinuous components on the line element i ($D_s^i$ and $D_n^i$) on the displacement and stress at any point in the operand can be calculated by the following equation [15,18].

The expression of the displacement and stress of line element j to the midpoint of line element i is as follows:

$$\begin{array}{l}{u}_s^j=D_s^j\left[(1-2\nu )\sin \theta {\overline{F}}_2+2(1-\nu )\cos {\overline{F}}_3-\overline{y}(\sin \theta {\overline{F}}_4+\cos \theta {\overline{F}}_5)\right]\\ +D_n^j\left[-(1-2\nu )\cos \theta {\overline{F}}_2+2(1-\nu )\sin \theta {\overline{F}}_3-\overline{y}(\cos \theta {\overline{F}}_4+\sin \theta {\overline{F}}_5)\right]\end{array}$$

(13)

$$\begin{array}{l}{u}_n^j=D_s^i\left[(1-2\nu )\cos \theta {\overline{F}}_2+2(1-\nu )\cos \theta {\overline{F}}_3-\overline{y}(\sin \theta {\overline{F}}_4+\cos \theta {\overline{F}}_5)\right]\\ +D_n^i\left[(1-2\nu )\sin \theta {\overline{F}}_2+2(1-\nu )\cos \theta {\overline{F}}_3+\overline{y}(\sin \theta {\overline{F}}_4+\cos \theta {\overline{F}}_5)\right]\end{array}$$

(14)

$$\begin{array}{l}{\sigma }_s^j=2G{D}_s^i\left[\sin 2\theta {\overline{F}}_4-\cos 2\theta {\overline{F}}_5+\overline{y}(\sin 2\theta {\overline{F}}_6+\cos 2\theta {\overline{F}}_7)\right]\\ +2G{D}_{\mathrm{s}}^i\left[-\overline{y}(\cos \theta {\overline{F}}_6+\cos \theta {\overline{F}}_7)\right]\end{array}$$

(15)

$$\begin{array}{l}{\sigma }_s^j=2G{D}_s^i\left[\sin 2\theta {\overline{F}}_4-\sin 2\theta {\overline{F}}_5+\overline{y}(\cos 2\theta {\overline{F}}_6+\sin 2\theta {\overline{F}}_7)\right]\\ +2G{D}_{\mathrm{s}}^i\left[-{\overline{F}}_5-\overline{y}(\sin 2\theta {\overline{F}}_6+\cos 2\theta {\overline{F}}_7)\right]\end{array}$$

(16)

2.5. Model Accuracy Verification

In order to verify the accuracy of the hydraulic fracturing model, as shown in Figure 4, fracturing data from several horizontal wells in Weiyuan X well area were selected, and data from 10 samples were input into the model for calculation. Based on the comprehensive filtration coefficient of the formation and the total volume of fracturing fluid injected into the formation, the volume of fracturing fluid used for fracturing, namely SRV, can be obtained. By simulating single-well SRVs with a different pumping flow rate at the same operation time, when the displacement is low, the net pressure inside the joint is insufficient to activate natural fractures. As the pumping flow rate increases, the net pressure inside the joint continues to rise. When the net pressure is increased to the required net pressure to activate the natural fracture, the natural fracture around the main fracture will be activated and opened, forming a complex fracture system. At this time, the SRV increases sharply. It can be seen from the sample data that when the displacement increases to 8 m³/min, hydraulic fractures begin to induce a large number of secondary fractures to open around them. Moreover, the numerical simulation results of the selected 10 sample wells have a high fit with the SRV of the sample wells in the field, which verifies the accuracy of the hydraulic fracturing model.

2.6. Model Parameter Settings

It is assumed that the reservoir has three layers: upper, middle and lower; the upper and lower layers are interlayers, the middle reservoir is 40 m thick, the fixed fracture height is the reservoir thickness, and the model is 2000 m long and 500 m wide. The basic parameter settings of fracking model are shown in

Table 1. Hydraulic fracturing model parameter settings.

Parameter	Top Interlayer	Middle Interlayer	Bottom Interlayer
Top depth (m)	2850	2900	2940
Bottom depth (m)	2900	2940	2990
Rock type	Shale
Porosity	0.01
Permeability (mD)	0.01
Young’s modulus (GPa)	20
Poisson’s ratio	0.25
Maximum horizontal principal stress (MPa)	39	43	44
Minimum horizontal principal stress (MPa)	40	44	45
Vertical stress (MPa)	41	45	46
Fracturing fluid	Slicwater
Propping agent	40/70 mesh quartz sand
Hole diameter (mm)	10
Single cluster perforation density (/m)	16
Flow rate (m3/min)	12
Filtration coefficient (m/s0.5)	2 × 10−5
Fracture toughness (MPa·m0.5)	1

.

Here, a coefficient of difference in inflow liquid volume was introduced to quantitatively evaluate the difference in inflow liquid volume into each cluster crack. This coefficient is defined by the standard deviation of the distribution of inflow liquid volume in each cluster:

$$S_d=\sqrt{\frac{1}{N_f-1}{\displaystyle \sum _{i=1}^{N_f}{\left(\frac{v_i}{v_t}-\frac{1}{N_f}\right)}^2}}$$

(17)

where V_t—total liquid volume, m³; S_d—the differential coefficient of the inflow of each cluster of liquid in the section, without a dimension.

The size of this coefficient can characterize the degree of difference in the amount of liquid in each cluster fracture during multi-cluster fracturing, and the amount of liquid in each cluster fracture determines the expansion form of each cluster fracture. The greater the difference in the amount of liquid in the cluster, the more uneven the crack expansion. Generally, 4.4 is taken as the limit, and if S_d is lower than 4.4, it is considered that the amount of liquid flowing into each cluster fracture is uniform, and the fracture expansion morphology is more uniform.

3. Simulation of Influencing Factors of Hydraulic Fracturing Fracture Parameters

There are two main types of factors affecting the formation of a complex fracture network by hydraulic fracturing, internal reservoir factors and external operation factors. The internal reservoir factors are mainly Young’s modulus, Poisson’s ratio, the in situ stress difference, the natural fracture approaching angle, etc. The external operation factors are mainly divided into the perforation cluster number, perforation spacing, pumping flow rate, fracturing fluid viscosity and so on. Only when the above two kinds of factors are considered together can a complex fracture network system be formed to carry out sufficient volumetric fracturing transformation of the reservoir. The following will carry out the quantitative research on the sensitivity of the relevant parameters from these two aspects.

3.1. Reservoir Influencing Factors

3.1.1. Young’s Modulus

Young’s modulus (E) is one of the important parameters affecting the propagation of hydraulic fractures, which reflects the difficulty of rock deformation under the action of external forces. The higher the brittleness of the rock, the easier it is to be pressed open. We simulated crack morphology using a single segment with three clusters of perforations and a cluster spacing of 30 m, with Young’s modulus set at 10 GPa, 20 GPa, 30 GPa, and 40 GPa. The simulation results are shown in Figure 5. With the same basic parameter settings, as the Young’s modulus of the reservoir increased, the length and width of each cluster of fractures in a single-stage three-cluster fracturing increased. The fracture lengths were 169.8 m, 189.7 m, 209.7 m, and 229.7 m, respectively. The deflection angle of the two sides of the fractures increased, and the middle fracture was squeezed by the induced stress of the two sides of the fractures, with a width smaller than that of the two sides and almost no deflection. According to the simulation results, under the condition of a certain fracture height, the integral area of each curve and coordinate was the fracture bottom area under the control of the same fracture height. That is, in Figure 6, although the elastic modulus was different, the scale of fracture 1 and 2 was the same. From the geometric meaning of the simulation results, it can be shown that fractures formed by reservoirs with a low elastic modulus are usually short and wide, while the fractures formed by reservoirs with a high elastic modulus are usually longer and narrower.

3.1.2. Poisson’s Ratio

Poisson’s ratio, which is the ratio of the transverse strain to axial strain of rock, has no effect on the final shape of hydraulic fracture, but Poisson’s ratio has a certain effect [19,20]. Generally speaking, the smaller the Poisson’s ratio is, the more likely the hydraulic fracture formed in the half length direction is to deform in unit time when the rock is subjected to shear force, that is, the faster the expansion speed. The fracture morphology was simulated by using a single section with three clusters of perforations, a cluster spacing of 30 m, and Poisson’s ratio of 0.1, 0.2, 0.3, 0.4, as shown in Figure 7. Under the same setting of basic parameters, with the increase in reservoir Poisson’s ratio, the change in the final fracture morphology was not significant, but the simulation end time was 30 s, 35 s, 38 s, 40 s, respectively. By extracting the fracture pressure values of four fracturing models in the simulation process and drawing them into curves (Figure 8), it can be seen from the curve change trend that with the increase in Poisson’s ratio, the fracture pressure of the shale reservoir will also increase, and the difficulty of rock cracking will also increase. In other words, under the same operation conditions, reservoirs with low Poisson’s ratio will show a certain lag in fracture expansion; that is, the smaller the Poisson’s ratio, the faster the fracture expansion speed.

3.1.3. In Situ Stress Difference and Approximation Angle

The crack approach angle refers to the angle between a newly generated crack and a natural crack when they intersect. Generally speaking, the smaller the geostress difference and approximation angle, the easier hydraulic fractures are to be captured by natural fractures. However, the larger the geostress difference and approximation angle, the easier hydraulic fractures are to extend forward through natural fractures [19,20]. According to the maximum likelihood diagram of well block X, the natural fractures in well block X in the Weiyuan area of Sichuan Basin were mainly low-angle (30°), and the distribution of in situ stress field had strong heterogeneity (Figure 9). Therefore, this work considered the comprehensive effect of in situ stress difference and fracture approach angle, and conducted orthogonal numerical simulation.

A single segment and four cluster horizontal well fracturing model with a cluster spacing of 30 m was established. The natural fracture strength was randomly distributed at 3 MPa and the standard deviation was 2 MPa to simulate the heterogeneity of specific reservoirs in the Weiyuan area. The natural fracture density was 0.002 pieces/m². By changing the combination of the approximation angle of natural fractures and the horizontal principal stress difference, orthogonal numerical simulations were conducted when the natural fracture azimuth angle was 30°, 45°, 60°, and the ground stress difference was 0, 5, and 10 MPa.

The simulation results show that when the approach angle is 30 ° and the geostress difference is 0 MPa, almost all hydraulic fractures are captured by natural fractures (Figure 10). When the hydraulic fractures extend to the natural fractures, they will extend along the weak bonding surface of the natural fractures [19,20], thereby activating the natural fractures. As the horizontal stress difference and approach angle increase, the hydraulic fractures tend to extend forward through the natural fractures, and the number of transverse fractures formed gradually increases. When the main stress difference is 10 MPa and the approach angle is 60°, almost all natural cracks are penetrated, and the formed crack network is similar to a grid, forming a complex fracture network.

The local stress difference gradually increased from 0 to 20 MPa, and the coefficient of difference in liquid inflow of the three different positions of the perforation cluster within the section increased from 4 to 8.2 (Figure 11). This indicates that as the geostress difference increases, the difference in liquid inflow into each cluster crack gradually increases. Specifically, the amount of liquid inflow into the outer crack is more than that into the inner crack, and the development of the inner crack is suppressed, with the size of the crack being smaller than that of the outer crack. After the local stress difference exceeds 12 MPa, the coefficient of difference in the liquid inflow of the perforation cluster at three locations within the section reaches 8.2, which is not conducive to the uniform expansion of cracks and makes it difficult to form a complex crack network.

3.2. Operation Factors

3.2.1. Number of Perforating Clusters

Changing the number of perforation clusters will change the morphology, density and length of the fracture, and affect. With the same amount of injected fluid, as the number of fracturing clusters increases, the size of a single fracture will continue to shrink (Figure 12). When conducting single section and two clusters of perforations, during the extension process of the two cracks, due to mutual interference of stress between the cracks, they will gradually deflect towards the direction away from each other. The shear forces on the two cracks are basically equal, so the deflection angle is also the same [19,20]. When single-stage three-cluster perforating was carried out, the middle crack was inhibited by the induced stress generated by the expansion of cracks on both sides, resulting in its length and width being smaller than the cracks on both sides. At the same time, there was almost no deflection due to the equal and opposite shear forces on both sides. When four clusters of fractures were fractured within a segment, the middle two clusters of fractures are severely disturbed by stress, and the half length of the fractures decreased from 190 m in a single segment with three clusters of perforations to 150 m. Compared to three clusters of perforations, the width of the middle fractures was also smaller; the pattern of five cluster perforations within the segment was consistent with the previous plan.

Figure 13 depicts the inflow of liquid into each crack cluster under different cluster numbers in a single segment. When the number of clusters in a single segment was less than 4, the interference of stress shadow between clusters was small, and the cracks expanded uniformly. When the number of clusters exceeded 4, the difference in liquid content between each cluster reached 4%, and the cracks began to expand unevenly. The outer cracks were deflected and became larger due to induced stress interference. Meanwhile, the simulation results in Figure 14 show that both the average control area of the mesh and the average deflection angle of the crack increased with the increase in the number of clusters, and the increased amplitude continued to increase. The number of clusters increased from 4 to 5, and the maximum increase in the control area of the mesh was 62.3%; the maximum increase in average deflection angle was 59%.

3.2.2. Perforation Spacing

The determination of cluster spacing is influenced by two factors: inter-cluster stress interference and seepage interference. When the distance between clusters was large, the interference of stress shadows between clusters was relatively small, and the crack propagation in each cluster of a single segment was more uniform. Due to the competition between clusters, there were differences in the amount of liquid entering each cluster crack. Specifically, the amount of liquid entering the outer crack was greater than that of the inner crack. However, the inner crack, due to the shear stress on both sides, did not deflect, but the crack size was smaller than that of the outer crack. Due to the more fracturing fluid entering the outer crack, its size was larger, resulting in an imbalance in the shear force on both sides of the outer crack, which can cause the crack to deflect [15]. Other factors being equal, we set the cluster spacing as 20 m, 30 m, 40 m and 50 m, and conducted single-section three-cluster perforation (Figure 15). When the cluster spacing was 20 m, the deflection angle of hydraulic fractures was relatively large, and the expansion of hydraulic fractures was mainly affected by induced stress. As the cluster spacing increased, the deflection angle of outer fractures decreased, and the fractures tended to be straight. The expansion of hydraulic fractures was mainly controlled by the original geostress field, and the fractures expanded along the direction of the maximum horizontal principal stress. At the same time, with the increase in the cluster spacing, the length of the hydraulic fracture (209.7 m) had little overall change, while the width gradually decreased and the height became larger (Figure 16). This is because in the basic setting of this simulation, the vertical fracture toughness was small: 0.5 MPa·m^0.5. In addition, the formation filtration coefficient was 2 × 10⁻⁵ m/s^0.5, which made the fracturing fluid filtration loss small and was basically used for fracture-making. In addition, the small vertical fracture toughness made the fracture easier to expand vertically.

3.2.3. Fracturing Fluid Viscosity

In large-scale hydraulic fracturing operations, the viscosity of fracturing fluid has an important influence on the complexity of fracture propagation. The lower the viscosity of the fracturing fluid, the higher the complexity of the fracture network. When high-viscosity fracturing fluid was used, there was obvious main fracture expansion, and the hydraulic fracture hardly intersected with the natural fracture, which was more likely to form a single fracture, while the low viscosity fracturing fluid was more likely to form a complex fracture network. Other parameters remained unchanged. The viscosity of the fracturing fluid was set as 5 MPa·s, 10 MPa·s, 20 MPa·s and 40 MPa·s, respectively (Figure 17). The density of natural fractures was 0.004 pieces/m², the length was 40 m and the azimuth angle was 30°. According to the simulation results, when the viscosity of fracturing fluid was low, secondary fractures were more developed. Hydraulic fracturing and natural fractures not only formed longitudinal fractures along the direction of maximum principal stress, but also transversed fractures perpendicular to the wellbore. When the viscosity of the fracturing fluid was high, the development trend of the main fracture gradually appeared, the hydraulic fracture was relatively straight, and the interaction with the natural fracture was less. The fracturing fluid viscosity had little influence on the hydraulic fracture width, and the change curve of influence on length is shown in Figure 18.

3.2.4. Pumping Flow Rate

With the increase in the pumping flow rate, the net pressure in the hydraulic fracture also increased; the smaller the pumping flow rate, the more fracturing fluid loss per unit time. When the displacement was low, the net pressure in the fracture was not enough to activate the natural fracture. However, with the increase in displacement, the net pressure in the seam gradually increases. When the net pressure increases enough to activate the natural fractures, the natural fractures are activated and interact with the hydraulic fractures to form a complex fracture network. Other parameters remain unchanged. Set the total injection volume of fracturing to 500 m³, the pumping flow rate to 3 m³/min, 8 m³/min, 12 m³/min, 16 m³/min, fracturing fluid viscosity to 10 MPa·s, natural fracture density to 0.004 pieces/m², length to 40 m, azimuth angle to 45° (Figure 19). According to the simulation results, it can be seen that when the pumping flow rate is 3 m³/min, the fracturing fluid displacement is less than the formation filtration rate during the pumping flow rate, and the fracturing fluid filtration loss is serious, and a large-scale fracture network cannot be formed. When the pressure fluid displacement increases to 12 m³/min, the optimal sewing mesh size has been formed. When the pressure fluid displacement continues to increase to 16 m³/min, the increase in seam mesh control area is minimal, but the operation cost increases more. The changes in crack length and width are shown in Figure 20. The larger the displacement, the longer and wider the crack. Therefore, the optimal pumping flow rate in Weiyuan area is 12 m³/min.

3.2.5. Fracturing Fluid Consumption

The amount of fracturing fluid determines the scale of fractures, and optimizing and researching the amount of fracturing fluid can provide a basis for the construction design of large-section and multi-cluster horizontal Wells in Weiyuan area of Sichuan Province. A horizontal well model with four clusters in a single section and cluster spacing of 30 m was set, other parameters remain unchanged and the fracturing fluid dosage was set as 2000 m³, 3000 m³ and 4000 m³, to study the influence of fracturing fluid dosage on fracture scale. According to Figure 21 of the simulation results, the length and width of fractures increased linearly with the increase in fracturing from 2000 m³ to 4000 m³. Considering that the well spacing of horizontal Wells in Weiyuan area of Sichuan is generally 400 m, excessive hydraulic fractures usually cause inter-well interference, it is recommended that the amount of fracturing fluid in a single stage should be controlled within 3000 m³ (Figure 22).

4. Conclusions

(1)

Both rock modulus and Poisson’s ratio are important indicators of rock brittleness, and are also the main parameters for formulating reservoir reconstruction schemes. The larger the Young’s modulus and Poisson’s ratio are, the more brittle the rock is, and the easier it is to crack during hydraulic fracturing. In the reservoir with higher Poisson’s ratio, the propagation rate of hydraulic fracture is faster.

(2)

When the in situ stress difference and approximation angle are small, the hydraulic fractures will not only induce the natural fractures to extend in the direction of the maximum principal stress, but also form transverse fractures, and finally form a grid-like fracture network. The low-angle (30°) natural fractures in Weiyuan X well area, Sichuan Province have the potential to form complex fracture networks. At the same time, the strong heterogeneity of geostress difference distribution brings many uncertain factors to the formation of fracture network. Therefore, when selecting the well location, it is necessary to combine the maximum likelihood chart and preferentially select the area with small in situ stress difference.

(3)

When the number of perforating clusters in a single stage is more than four clusters, the interfracture interference is more serious. When the cluster spacing is less than 30 m, the fracture shape is distorted, which is the premise of forming a complex fracture network.

(4)

The viscosity of fracturing fluid has no obvious influence on the width of fracture expansion, but has great influence on the length of fracture. The higher the viscosity of the fracturing fluid, the more obvious the development trend of the main fracture. When the fracturing fluid with low viscosity (10 MPa·s) is used, the more obvious the development trend of the main fracture is or does not show the development trend of the main fracture, and more secondary fractures are activated to form a complex fracture network with the main fracture.

(5)

The construction displacement is related to the net pressure inside the fracture, and only enough net pressure inside the fracture can open the natural fracture. When the construction displacement is 3 m³/min, the filtration loss is more serious, and when it is above 12 m³/min, the fracture network can reach the best scale. It is suggested that the average construction displacement of large-scale hydraulic fracturing operations in Weiyuan shale gas Wells is 12 m/min.

(6)

The amount of fracturing fluid determines the scale of the fracture. In order to avoid inter-well interference, it is recommended that the maximum amount of fracturing fluid used in single section of horizontal well construction in Weiyuan area should be 3000 m³.