Numerical Simulation of Influencing Factors of Hydraulic Fracture Network Development in Reservoirs with Pre-Existing Fractures

Abstract

The intersection behavior of hydraulic fractures and single natural fractures has been studied in detail; however, in fractured reservoirs, natural fractures are numerous and interlaced and the intersection of hydraulic fractures and multiple natural fractures occurs during the fracturing process. This intersection behavior is more complex and there is a lack of research on this topic at present. In this study, a numerical simulation model of the interaction between hydraulic fractures and a series of natural fractures was established, the main factors that affect the formation scale of a fracture network during the hydraulic fracturing of a fractured reservoir were studied using the numerical simulation method, and the parameters were also studied. The results showed that the natural fracture trend, in situ stress difference, and injection flow rate have an impact on the scale of a fracture network. The larger the in situ stress difference, the smaller the scale of the fracture network, which gradually changes from multiple clusters of fractures to single fractures. The larger the injection flow rate, the larger the scale of the fracture network. In the uniform stress field, the direction of a natural fracture is closer to the direction of principal stress, so the lower the fracture extension pressure, the smaller the scale of the network. On the contrary, the farther away from the principal stress direction, the lower the fracture extension pressure and the higher the extension pressure, the larger the scale of the fracture network.

1. Introduction

For tight reservoirs, hydraulic fracturing is the most commonly used technique to improve the seepage capacity and increase the production of reservoirs [1,2,3,4,5], as the size of the fractures that are formed by hydraulic fracturing directly affects final production [6,7,8]. The methods that are commonly used for the numerical simulation of hydraulic fracturing include the extended finite element method (XFEM), the displacement discontinuity method (DDM), the phase field method, and the cohesive zone method (CZM). Dahi found from field hydraulic fracture diagnostic data that the fractures that are produced by fracturing a formation in which natural fractures are present do not produce ideal symmetrical bifurcations, but rather they deflect and interact with the natural fractures. He divided fracture intersection behavior into three categories: hydraulic fractures that are arrested by natural fractures, hydraulic fractures that cross natural fractures, and hydraulic fractures that expand together with natural fractures [9]. Zhang pointed out that the angle of deflection occurs when a hydraulic fracture intersects and crosses a natural fracture and that the angle of deflection is related to the ratio of the two ground stresses in the plane [10,11,12,13]. Yao developed an energy-based fracture intersection criterion from the Griffith’s criterion and produced a fracture intersection criterion that is related to the stress difference and approach angle [14]. Luo found that as the angle between the hydraulic fracture and the direction of the maximum horizontal stress increases, the likelihood of the hydraulic fracture crossing a natural fracture with the same angle of approach first increases and then decreases, according to XFEM simulations [15]. However, these fracture intersection criteria were studied using intersections of single hydraulic fractures and single natural fractures and did not consider the effects of stress shadowing when the tip of the hydraulic fracture is close to the natural fracture face. To explore the expansion law of hydraulic fractures in multiple natural fracture strata, Li simulated the expansion trajectory of hydraulic fractures when there were multiple natural fractures of the same orientation in the strata using XFEM simulations [16]. Hou simulated the complex morphology of hydraulic fracture extension in shale formations under the disturbance of randomly distributed natural fractures using the displacement discontinuity method [17,18]. Zhang conducted a study using the discrete element method and indicated that the behavior of hydraulic fractures is strongly related to the hydrostatic pressure within the fractures [19]. The driving force for the deflection of the direction of hydraulic fracture extension is mainly the induced stress field that is generated by the natural fractures and the more natural fractures there are, the more complex the hydraulic fracture extension. Hou simulated the hydraulic fracturing process in shale formations that contained randomly distributed laminae and high-angle fractures using 3DEC and showed that it is a competitive extension process [20,21].

The natural fracture directions in these studies were either the same or completely random. In contrast, an observation of outcrops revealed the existence of a special natural fracture distribution pattern, in which natural fractures extend along two directions with each occupying approximately the same general area, in the manner of a checkerboard. This natural fracture distribution pattern is caused by a shift in the direction of ground stress during the geological history of the strata. In this paper, we used the CZM method to simulate hydraulic fracturing in a formation of natural fractures with a checkerboard-liked distribution, studied the effects of differences in ground stress and the fracturing fluid injection rate on the final fracture network size under two fracture orientations, and developed a method to optimize the perforate direction and fracturing fluid injection rate in order to provide guidance for improving the fracturing effects in such reservoirs.

2. Governing Equations

Cohesive elements are a special type of element, which was used in this study to represent the possible path of a hydraulic fracture. Firstly, the model needed to be meshed; then, the originally continuously distributed mesh was discretized and the cohesive element was inserted in between the elements to reconstruct the continuous model. During the process of analysis and simulation, the coupling analysis of multiple physical fields was carried out. The mechanical response of the cohesive unit was described by the critical failure criterion of fracture mechanics and the damage evolution criterion after the failure. Different equations were used to describe the fluid migration, according to the degree of failure and the opening of the cohesive unit (after complete failure). In this way, the simulation was completed.

When simulating the hydraulic fracturing process, the initiation and propagation process of a fracture described the deformation and failure process of the cohesive element. Unlike the description of the stress state of a conventional element, which requires the use of six independent components, the stress state of the cohesive element was reduced to three principal directions in the space among them, $t_n$, $t_s$, and $t_t$, as shown in Figure 1.

The relationship between traction and separation is shown in Figure 2. Before the cohesive element was damaged, the relationship between the traction and separation was linear and could be expressed by Equation (1):

$$\left\{\begin{array}{c}{t}_n\\ t_s\\ t_t\end{array}\right\}=\left[\begin{array}{ccc}{k}_{nn}& & \\ & k_{ss}& \\ & & k_{tt}\end{array}\right]\left\{\begin{array}{c}{\delta }_n\\ {\delta }_s\\ {\delta }_t\end{array}\right\}$$

(1)

where $k_{nn}$, $K_{ss}$, $K_{tt}$ is the stiffness in three principal directions of the cohesive element, Pa/m; $t_n$, $t_s$, $t_t$ is the traction in three principal directions of the cohesive element, Pa; and ${\delta }_n$, ${\delta }_s$, ${\delta }_t$ is the separation in three principal directions of the cohesive element, m.

As the traction increased, when the three directions of traction satisfied the relationship in Equation (2), the cohesive element started to become damaged:

$${\left\{\frac{\langle t_n\rangle }{t_n^0}\right\}}^2+{\left\{\frac{t_s}{t_s^0}\right\}}^2+{\left\{\frac{t_t}{t_t^0}\right\}}^2=1$$

(2)

where the symbol $\langle \rangle$ represents Macaulay brackets, which are used to signify that a pure compressive deformation or stress state does not initiate damage; and $t_n^0$, $t_s^0$, $t_t^0$ is the strength of the cohesive element in three principal directions, Pa.

The damage was a gradual process and was not instantaneous. The relationship between the traction and separation during this period could be represented by Equation (3):

$$\{\begin{array}{l}D=\frac{{\delta }_m^f\left({\delta }^{\max }-{\delta }_m^0\right)}{{\delta }^{\max }\left({\delta }_m^f-{\delta }_m^0\right)}\\ t_n=\{\begin{array}{ll}(1-D)\overline{t_n}& \left(\overline{t_n}>0\right)\\ \overline{t_n}& \left(\overline{t_n}<0\right)\end{array}\\ t_s=(1-D)\overline{t_s}\\ t_t=(1-D)\overline{t_t}\end{array}$$

(3)

where $D$ is the damage variable; and $\overline{t_n}$, $\overline{t_s}$, $\overline{t_t}$ is the predicted traction that was calculated using Equation (1). When $D$ is equal to 1, the cohesive element is completely destroyed and fractures are formed.

The forming of fractures had a great influence on the fluid transport within the formation. The flow within the porous media was Darcy flow, as shown by Equation (4):

$$q_r=-\frac{k}{{\mu }_w}\nabla p$$

(4)

where $q_r$ is the flow velocity in porous media, m/s; $k$ is the permeability of the formation, m²; ${\mu }_w$ is the viscosity of the fluid, $\mathrm{Pa}\cdot \mathrm{s}$; and $\nabla p$ is the gradient of the pore pressure.

The flow inside the fracture was simplified into a tangential flow along the fracture extension direction and a normal flow that was perpendicular to the fracture direction. Assuming that the fluid was a Newtonian fluid and neglecting the pressure gradient normal to the fracture, the controlling equation of the fluid inside the fracture could be expressed as Equation (5):

$$\{\begin{array}{l}{q}_h=-\frac{d^3}{12{\mu }_w}\nabla p_h\\ q_t=C_t\left(p-p_t\right)\\ q_b=C_b\left(p-p_b\right)\end{array}$$

(5)

where $q_h$ is tangential flow rate in the fracture, m/s; $\nabla p_h$ is the tangential fluid pressure gradient in the fracture, Pa/m; $d$ is the fracture opening, m; $q_t$, $q_b$ is the flow rate of the fluid in the fracture leakoff to both sides of the fracture wall, m/s; $p_t$, $p_b$ is the pore pressure on the outside of the fracture wall, Pa; $C_t$, $C_b$ is the leakoff coefficient, $\mathrm{m}/\left(\mathrm{s}\cdot \mathrm{Pa}\right)$.

3. Numerical Model for Hydraulic Fracturing Simulation

When there is only one fracture in a formation, the fracture expands along the horizontal direction of the maximum ground stress. However, there was a series of natural fractures in the actual formation that was studied here, and the presence of the fractures created an additional stress field. The more natural fractures there are, the more complex the overall stress field of the formation becomes, which makes the paths of the hydraulic fractures during the fracturing process unpredictable.

In this paper, a numerical simulation model was established, which had a length and width of 10 m and a thickness of 0.1 m. Then, two groups of natural fracture networks were randomly generated. Each group consisted of 54 fractures of different lengths: the first group had 27 fractures with the angle of −20° from the horizontal and 27 fractures with the angle of 80° from the horizontal, while the second group had 27 fractures with the angle of 40° from the horizontal and 27 fractures with the angle of 140° from the horizontal. In order to prevent the local enrichment of the randomly generated natural fractures in the plane from affecting the simulation results, a large sample of random fractures (10,000 fractures) was first generated using the random fracture generation and then equally spaced sampling was performed on them to ensure the uniform distribution of the natural fractures without losing randomness. The schematic diagram of the two sets of natural fracture networks is shown in Figure 3 below.

The fracturing fluid that was used in the simulation was water; its viscosity was assumed to be constant and was taken as 1. The water was injected at a constant flow rate into the center of the model during the fracturing process. The maximum horizontal in situ stress was in the lateral direction of the model and the minimum horizontal in situ stress was in the longitudinal direction of the model. The other parameters that were used in the model are shown in

Table 1. The simulation parameters.

Main Parameter	Value
Young’s modulus (GPa)	40
Poisson’s ratio	0.16
Density (kg/cm3)	2300
Horizontal maximum in situ stress (MPa)	66
Horizontal minimum in situ stress (MPa)	61
Vertical in situ stress (MPa)	78
Pore pressure (MPa)	30
Permeability (mD)	1
Porosity	0.01
Hydraulic fracture tensile strength (MPa)	6
Nature fracture tensile strength (MPa)	0
Hydraulic fracture shear strength (MPa)	6
Nature fracture shear strength (MPa)	0
Leakoff coefficient (×10−11 $\mathrm{m}/\left(\mathrm{s}\cdot \mathrm{Pa}\right)$)	5

.

4. Results and Discussion

4.1. Fracture Pressure Characteristics

Fracturing simulations were carried out using 0 MPa, 3 MPa, and 5 MPa as the plane stress difference. The pore pressure at the injection points of the natural fracture networks of group A and group B varied with time during the fracturing process, as shown in Figure 4, Figure 5, Figure 6 and Figure 7.

The shape of the curves under different ground stress differences was basically the same, as shown in Figure 4. All of them showed a pressure increase at the beginning that was followed by a sudden decrease due to the formation of the fracture, after which the fracture extended and curve became smoother. The pressure of the formation of a fracture was about 20 MPa under different ground stress differences, between which there were no significant differences. A fracture pressure of 3 MPa during the last 90 s of the fracture extension stage showed a small decrease, followed by a gradual increase. On the contrary, a fracture pressure of 5 MPa showed a small increase that was followed by a gradual decrease and a pressure of 0 MPa increased uniformly throughout the extension stage.

The displacement contour plot of the natural fracture network of group A at different moments under a ground stress difference of 0 MPa is shown in Figure 5. The figure shows that the fractures mainly extended along the longitudinal direction and a lateral displacement was produced near the injection point due to the opening of the longitudinal fractures, which caused a lateral fracture that extended to the right under the combined effects of the high pore pressure. Natural fractures bifurcate as they extend due to natural fracture influence. Since the ground stress difference was 0 MPa, the resistance to fracture extension was the same in all directions; so, the final fracture morphology was strongly related to the fracture direction that was formed by the initial rupture and was secondly influenced by the natural fracture distribution. The main fracture in the longitudinal direction and the orientation of the fracture in the transverse direction, which are shown in the figure, both well confirmed the influence of the orientation of the natural fractures on the overall orientation of the fracture network that was formed by the hydraulic fracturing.

The pore pressure of the injection point in the natural fracture network of group B changed with time, as shown in Figure 6. Compared to group A, the pore pressure at the water injection point of group B changed more during the fracturing process. Pressure drops also appeared after the initial fracture of the formation, but the pressure increased again immediately and then directly entered the fracture extension stage. The pressure curves under different stress differences in the fracture extension stage were quite different. When the stress difference was 0 MPa, the curve obviously had greater fluctuations. When the fracturing simulation with a stress difference of 3 MPa was carried out to 85 s, the pressure suddenly decreased by about 3 MPa and then gradually increased. This phenomenon also occurred in the fracturing simulation with a stress difference of 0 MPa, but it occurred earlier because the larger the stress difference, the easier the natural fracture is to shear and slide. Then, the fracture opened in the process of pore pressure increase, which resulted in a pressure drop that was similar to that of the formation of a fracture. When the stress difference was 5 MPa, the pore pressure at the injection point rose steadily during the whole fracture propagation process without fluctuation because the natural fractures in the critical state had opened in advance under the large in situ stress difference and the unopened natural fractures could not shear or slide under that formation pressure. Compared to the simulation results of group A, the pore pressure of group B under different in situ stress differences gradually increased, which showed that the trend of the natural fractures in the formation had an impact on the fracture extension pressure.

The displacement contour plot of group B is shown in Figure 7. Similar to the results from group A, the displacement at the water injection point was the largest, up to 4 mm, and there was an obvious shear phenomenon near to the fracture. Shear fractures tend to be self-supporting and enhance connectivity between fractures, which has an important effect on reservoir reconstruction.

4.2. Effects of In-Situ Stress Difference

To evaluate the effects of ground stress difference on the fracturing results of the two sets of natural fracture networks, the ground stress difference was taken to be 0 MPa, 3 MPa, and 5 MPa with a constant water injection rate of 0.005 m³/s.

The fracturing simulation results of group A and group B under different in situ stress differences are shown in Figure 8 below. The simulation results of both groups showed that with an increase in stress difference, the scale of the fracture network gradually changed from a complex multi-cluster of fractures into a simple single fracture, which was more obvious in group A. With a stress difference increase from 0 MPa to 3 MPa, the fractures in group B decreased in length but the network still contained a series of clusters. However, the fractures in group A changed from three clusters to a single fracture, which significantly decreased the area of the fractured zone. By comparing group A and group B under the same in situ stress difference, it was found that more natural fractures in group B sheared open under the influence of the in situ stress difference and the fracture network that was formed by the fracturing process was also more complex. This showed that the influence of in situ stress difference on the hydraulic fracturing of a fractured reservoir is also related to the trend of the natural fractures. When the angle between the natural fracture strike and the maximum in situ stress direction was large or small, the scale of the fracture network was small, while when the fracture strike was close to 45°, the scale of the fracture network was large.

Figure 9 shows the distribution of the fracture criteria from Equation (2) during the fracturing of group A and group B under different stress differences. The results showed that with the increase in in situ stress difference, more and more natural fractures reached the failure conditions; however, because there was not enough fluid in the fractures, the fractures closed. Under the same stress difference, the number of natural fractures in group B that met the fracture criterion was significantly higher, which further reflected that when the natural fracture strike was different, the influence of the in situ stress difference on whether it was damaged was also different.

4.3. Effects of Fracturing Fluid Inject Rate

In order to compare the fracturing effects of the two groups of natural fracture networks under different fracturing fluid injection rates with a stress difference of 0 MPa, the hydraulic fracturing simulations of 500 s were carried out with inject rates of 0.001 m³/s, 0.005 m³/s, and 0.01 m³/s. The simulation results are shown in Figure 10 below.

The simulation results showed that with the increase in flow rate, the scale of the fracture network increased, but this was not obvious. At the same time, some natural fractures that were far away from the main fractures also opened.

For group A, new clusters appeared on the left after increasing the flow rate, which was caused by the intersection of natural fractures. Increasing the flow rate did not produce additional hydraulic fractures in group A. In group B, however, increasing the flow rate not only opened more natural fractures but also produced additional hydraulic fractures, which then formed some small clusters and better improved the fracturing effect. This showed that for reservoirs with different trends of natural fractures, a greater flow rate does not necessarily produce a better fracturing effect. When the natural fracture strikes were close to the principal stress direction, an increase in the flow rate was of little significance.

The pore pressure distribution of group A and group B under different injection flow rates is shown in Figure 11 below. The results showed that under the same injection rate, the pore pressure of the formations around the fracture networks of group A and group B was about the same, but the pore pressure in the fractures was different. For example, when the injection flow rate was 0.01 m³/s, the pore pressure of the formation around the fracture network in group A was about 13 MPa, while that in group B was 16 MPa and the pressure of the fracture tip was as high as 24 MPa. This difference showed that the direction of the natural fractures also had a great influence on the extension pressure of the hydraulic fractures. In group A, when the hydraulic fractures met natural fractures, it was mainly the hydraulic fractures that intersected with the natural fractures. At this time, the direction of the hydraulic fractures was consistent with that of the natural fracture and the natural fractures did not produce any resistance to the expansion of the hydraulic fractures. However, in group B, the hydraulic fractures either intersected vertically with the natural fractures or turned because of the natural fractures. When the hydraulic fractures intersected with the open natural fractures, the pressure of the fluid in the natural fractures was equal to the resistance of the hydraulic fracture propagation. When the hydraulic fractures turned due to the natural fractures, it showed that the induced stress field of the natural fractures had an influence on the hydraulic fractures. All of this shows that the natural fractures in group B seriously affected the expansion of the hydraulic fractures; so, the pressure in group B was higher than that in group A, but at the same time, the scale of fracture network that was formed by group B was also larger.

5. Conclusions

(1).

With the increase in the in situ stress difference from 0 MPa to 3 MPa, the fracture that was formed by the fracturing process gradually changed from a complex multi-cluster fracture network into a simple single fracture and the fracturing effect decreased gradually. This change was especially obvious when the natural fracture direction was close to the principal stress direction;

(2).

The greater the flow rate of the injection point, the higher the fluid pressure in the fracture; however, there was no obvious change in pore pressure in the formation around the fracture. Meanwhile, the greater the flow rate of the injection point, the larger the scale of the fracture network that was formed;

(3).

The natural fracture direction had an influence on the fracture network scale and the fracture extension pressure during the fracturing process. In the uniform stress field, while the natural fracture trend was closer to the direction of principal stress, the scale of the fracture network was smaller and the extension pressure was lower.