Geofluids
Volume 2019 (2019), Article ID 9402392, 16 pages
https://doi.org/10.1155/2019/9402392

Research Article

Numerical Simulation of Artificial Fracture Propagation in Shale Gas Reservoirs Based on FPS-Cohesive Finite Element Method

Xiaoqiang Liu,¹ Zhanqing Qu,¹ Tiankui Guo,¹ Ying Sun,¹ Zhifeng Shi,² Luyang Chen,³ and Yunlong Li¹

¹College of Petroleum Engineering, China University of Petroleum (East China), Qingdao 266580, China
²The Fifth Oil Production Plant, PetroChina Huabei Oilfield Company, Xinji 052360, China
³Chunliang Oil Production Plant, Shengli Oilfield, Sinopec, Binzhou 256500, China

Correspondence should be addressed to Zhanqing Qu; quzhq@upc.edu.cn and Tiankui Guo; guotiankui@126.com

Received 18 September 2018; Accepted 9 December 2018; Published 17 February 2019

Academic Editor: Shengnan Nancy Chen

Copyright © 2019 Xiaoqiang Liu et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Abstract.

The simulation of hydraulic fracturing by the conventional ABAQUS cohesive finite element method requires a preset fracture propagation path, which restricts its application to the hydraulic fracturing simulation of a naturally fractured reservoir under full coupling. Based on the further development of a cohesive finite element, a new dual-attribute element of pore fluid/stress element and cohesive element (PFS-Cohesive) method for a rock matrix is put forward to realize the simulation of an artificial fracture propagating along the arbitrary path. The effect of a single spontaneous fracture, two intersected natural fractures, and multiple intersected spontaneous fractures on the expansion of an artificial fracture is analyzed by this method. Numerical simulation results show that the in situ stress, approaching angle between the artificial fracture and natural fracture, and natural fracture cementation strength have a significant influence on the propagation morphology of the fracture. When two intersected natural fractures exist, the second one will inhibit the propagation of artificial fractures along the small angle of the first natural fractures. Under different in situ stress differences, the length as well as aperture of the hydraulic fracture in a rock matrix increases with the development of cementation superiority of natural fractures. And with the increasing of in situ horizontal stress differences, the length of the artificial fracture in a rock matrix decreases, while the aperture increases. The numerical simulation result of the influence of a single natural fracture on the propagation of an artificial fracture is in agreement with that of the experiment, which proves the accuracy of the PFS-Cohesive FEM for simulating hydraulic fracturing in shale gas reservoirs.

1. Introduction

As an important unconventional fossil hydrogen energy source, shale gas reservoirs have been widely developed in recent years. Due to its ultralow permeability, large-scale hydraulic fracturing has become a necessary technology for economic exploitation [1]. Shale reservoirs generally develop weak bedding surfaces and natural fractures. Artificial fractures will intersect with natural fractures to get a complex fracture net during the process of hydraulic fracturing [2–4]. A correct understanding of the intersection mechanism between multinatural fractures and artificial fractures is of great importance in guiding the hydraulic fracturing operation and increasing the stimulated reservoir volume.

Many scholars have studied the propagation pattern of an artificial fracture intersecting with a single natural fracture by physical modeling and numerical modeling [5, 6]. It is generally believed that there are three situations of the intersection between an artificial fracture and a single natural fracture [7–11]: the artificial fracture intersects with the natural fracture, the natural fracture intercepts, and the artificial fracture is deflected into the natural fracture. Rock samples from shale gas reservoirs were scanned by microseismic monitoring technology, and the results showed that there are multi-intersected natural fractures in the rock samples [12–14]. Zhao et al. analyzed the tendency of a propagating angle and the impacts of propagating pressure of a complex fracture network influenced by natural fractures [15]. Hou et al. proposed propagation area evaluation of artificial fracture networks in shale gas reservoirs and analyzed the possibility that the artificial fracture communicates with the natural fractures [16]. However, few scholars have researched the effects of multi-intersected natural fractures on artificial fracture propagation. And there is also no clear expression of an artificial fracture propagation pattern in shale reservoirs with multi-intersected natural fractures.

The previous numerical study of the intersection of an artificial fracture with a natural fracture mainly adopts linear elastic mechanics, but this method ignores the fluid-solid coupling during hydraulic fracturing. The finite element method (FEM) and the extended finite element method (XFEM) and the discontinuous displacement methods based on the boundary element (DDM) and the discrete fracture network model (DFN) are gradually proposed [17–19]. Due to the complexity of dealing with natural fractures and mesh division, the above methods are seldom used in the large-scale hydraulic fracturing simulation with multinatural fractures under full coupling. In recent years, the ABAQUS cohesive finite element method is introduced to simulate the intersection of artificial fractures with natural fractures and some scholars have carried out related research. Guo et al. [20] used the cohesive finite element methods to build up a two-dimensional model for the first time. The fracture propagation path was presented to analyze the influence of the approaching angle, horizontal stress difference, viscosity, and injection rate of fracturing fluid on the propagation pattern of artificial fractures intersecting with a single natural fracture. Gonzalez et al. [21] used the cohesive finite element method to study the role of natural fracture cementation strength on the propagation of artificial fractures. Four triangular cohesive elements were taken to ensure the continuity of pressure in the cohesive element. Haddad et al. [22] established a three-dimensional model of the intersection of a single artificial fracture and a single natural fracture. A pore fluid/stress element was used for a rock matrix, and a cohesive element was used for an artificial fracture and a natural fracture. Xu et al. [23] studied the impact of an approaching angle and initiation pressure on the propagation of artificial fractures and the opening of a natural fracture by cohesive element method. Zhang et al. [24] used the cohesive finite element method to study the influence of stress, fracturing parameters, and fracture geometry on the stress interference between an artificial fracture and a natural fracture. Although the multifracture interference was involved in the model, it still used the preset fracture propagation path.

In general, the effect of an intersected natural fracture on artificial fracture propagation is not yet clear. And the conventional ABAQUS cohesive finite element simulating hydraulic fracturing requires a preset fracture propagation path, which restricts its application to the hydraulic fracturing simulation. And the cohesive finite element is only used for simulation of an artificial fracture intersecting with a single natural fracture. While intersected natural fractures are widely distributed in shale reservoirs, the effective connection of a multi-intersected natural fracture, as well as the free propagation of artificial fractures in a rock matrix, has become necessary for numerical simulation of large-scale hydraulic fracturing in shale reservoirs. In this paper, a new dual-attribute element of pore fluid/stress element and cohesive element (PFS-Cohesive) method for a rock matrix is presented based on the further development of a cohesive finite element, which can realize artificial fracture propagation along the arbitrary path. The effect of a multi-intersected natural fracture on artificial fracture propagation is analyzed by this new method. The results obtained can be used to explain the expansion mechanism of artificial fractures in shale reservoirs.

2. Computational Model

Rock stress equilibrium equation, continuity equation for fluid flow in porous media, and fracture flow model are proposed in this paper, and the effects of matrix poroelastic deformation [25–27] are considered in numerical simulation. The coupling of rock stress equilibrium equation and fluid flow continuity equation in porous media is realized through the effective stress concept. The fluid flow continuity equation in porous media is linked to fracture flow model by leak-off coefficient. Natural fractures are simulated by a cohesive element, and the dual-attribute element of a pore fluid/stress element and cohesive element (PFS-Cohesive) method is developed by Fortran for a rock matrix. Fracturing fluid and liquid in rock pores has influence on rock matrix stress, porosity, and permeability. And matrix poroelastic deformation also affects the fracturing liquid flow in turn.

2.1. Fluid-Solid Coupling Equation

The definition of effective stress of a rock matrix in saturated single-phase liquid is proposed as [27]

where is the effective stress (Pa), is the total stress (Pa), is the Biot coefficient, is the matrix fluid pressure (Pa), and is the second-order identify tensor. The definition of the Biot coefficient is

where is the bulk modulus of a rock (Pa) and is the bulk modulus of mineral grains (Pa).

The relationship between stress and strain could be expressed as incremental forms as

where is the elastic-plastic matrix and is the rock compression strain caused by pore fluid pressure, which can be written as

where is the rock pore pressure (Pa) and .

The rock stress equilibrium equation could be expressed by the principle of virtual work. The virtual work of the rock is equivalent to the virtual work produced by the force (including body force and surface fore) acting on the rock.

where is the virtual displacement (m).

Combining formulas (1)–(4), formula (5) can be written as

After some manipulation, formula (6) can be changed as [26, 27]

where is an elastic-plastic matrix with a time step, is saturation, is for the rock pore pressures (Pa), and

The flow of liquid in porous media conforms to Darcy flow. Darcy’s law of porous medium flow is written as

where is the liquid flow velocity (m/s), is permeability (×10⁻¹⁵ μm²), is pressure difference (Pa), and is liquid density (kg/m³).

The mass conservation for fluid flow in a rock is as follows [28]:

where is saturation of liquid in a porous rock and is fluid volume modulus (GPa).

2.2. Fracture Flow Model

The flow of liquid in a fracture includes tangential flow and normal flow, shown in Figure 1. The tangential flow mode is used to describe the fluid flow along the direction of the fracture, and a normal flow way is worked to simulate fracturing fluid leak off into the rock matrix. The liquid is supposed to be incompressible Newtonian fluid, and tangential flow is governed by lubricating equation [29]

where is tangential fluid flow (m/s), is tangential permeability, is tangential flow fluid pressure gradient (Pa/m), and is fracture aperture (m). Based on Reynold’s equation, the tangential permeability is defined as

where is viscosity of liquid in the fracture (Pa·s).

Figure 1: Fluid flow in a cohesive element.

Normal flow indicates leakage of fracturing fluid to porous medium and is simulated by defining the leakage coefficient of a rock matrix, as shown in Figure 2. The normal flow equation can be written as [29]

where and are flow rates of the top surface and bottom surface, respectively, (m/s), and are the top leakage coefficient and bottom leakage coefficient, respectively, and are the pore pressure at the top surface and bottom surface, respectively, (Pa), and is the fluid pressure in the fracture (Pa).

Figure 2: The leakage coefficient is shown as a permeable layer.

2.3. Fracture Initiation Criteria and Damage Evolution Law

ABAQUS provides a variety of fracture initiation criteria based on stress and strain. The secondary stress criteria are adopted in this model. When the sum of the ratio of the stress in a normal direction and tangential direction to the corresponding critical stress reaches 1, the fracture starts to crack. The fracture initiation criteria can be written as [30]

where is normal stress of the element in a normal direction (Pa), and are shear stresses of the element in two tangent directions (Pa), is critical stress in a normal direction when the element fails, also defined as the tensile strength of the rock (Pa), and are critical stresses in two tangent directions when the element fails, also defined as the shear strength of the rock (Pa), and indicates that no damage occurs when the element is subjected to pressure stress, defined as

Damage evolution refers to the energy needs for further damage of a rock after initial crack. This paper introduces the Benzegagh-Kenane criteria as a damage law of artificial fracture growth during the hydraulic fracturing process, which can be represented as [31]

where and are energy correlation rates in directions and , constants depend on material properties, is critical fracture energy release rates for composite fractures, and and are critical energy release rates in tangential and normal directions.

2.4. The Dual-Attribute Element of the PFS-Cohesive Method

The cohesive element is embedded in the pore fluid/stress element, and the dual-attribute element of pore fluid/stress element and cohesive element (PFS-Cohesive) method for a rock matrix is presented, shown in Figure 3.

Figure 3: Cohesive finite element equipped with pore fluid/stress finite element (PFS-Cohesive) in a rock matrix.

The dual element of the rock matrix consists of 6 nodes, of which there are two displacement degrees of freedom (u1-u8) in no. 1-4 nodes. The no. 5 and no. 6 nodes have pore pressure freedom (P1, P2) only, and the injected fracturing fluid pressure and flow rate are dispersed in the middle layer of no. 5 and no. 6. The pore fluid/stress element is used to simulate the property of the reservoir rock, such as the permeability of formation, Young’s modulus, and Poisson’s ratio. The cohesive element is used to simulate the occurrence and propagation of fracture, as well as liquid flow in the fracture. The two element types share the same nodes of nos. 1-4.

2.5. Model Description

Three models are established by using the PFS-Cohesive method (shown in Figure 4). The natural fracture is simulated by a cohesive element, and the rock matrix is simulated by the dual-attribute element of the PFS-Cohesive method. The size of the first model is , and the natural fracture length is 9 m. The effect of a single natural fracture on the artificial fracture propagation morphology is analyzed by using this model. The main purpose of this model is to verify the accuracy of the PFS-Cohesive method by comparing with experimental results. The size of the second model is . The length of the long natural fracture is 9 m, and the length of the short natural fracture is 5 m. This model is used to analyze artificial fracture propagation with two intersected natural fractures. The size of the third model is . The natural fracture length is between 2 m and 9 m, and natural fractures are intersected or discrete in 60° with maximum in situ horizontal stress. The third model is used to analyze the influence mechanism of multiple natural fractures on artificial fracture propagation morphology at reservoir scales. The size of the element has a significant influence on the fracture propagation pattern. If the size of the element is too large, the fracture will change its propagation path abruptly along the edge of the element, which will affect the accuracy of the model. While if the size of the element is too small, it has little effect on fracture propagation but will reduce the calculation efficiency of the model. In order to ensure the calculation accuracy of the model and improve the calculation efficiency, local mesh refinement is applied to the fracture propagation area. After many practices, in the first model and the second model, the spacing of grid nodes of the fracture propagation area is 0.1 m and the spacing of grid nodes of the model boundary is 0.5 m, which can ensure the accuracy of model results. In the third model, due to the large scale of model and the wide distribution of natural fractures, after many practices, the spacing of grid nodes of the fracture propagation area is 0.2 m and the spacing of grid nodes of the model boundary is 0.5 m, which can also ensure the accuracy of the model. The basic parameters of the three models are shown in Table 1.

(a) Single natural fracture

(b) Two intersected natural fractures

(c) Multi-intersected or discrete natural fractures

Figure 4: Schematic model diagram.

Table 1: Basic model parameters.

Parameter
Value
Parameter
Value

Maximum in situ horizontal stress
36 MPa
Minimum in situ horizontal stress
26/31/36 MPa

Overburden stress
40 MPa
Initial pore pressure
24 MPa

Young’s modulus
10 GPa
Poisson’s ratio
0.24

Porosity
0.02
Permeability
0.06 mD

Tensile strength of the rock
2.5 MPa
Tensile strength of the natural fracture
2/1.7/1.5 MPa

Shear strength of the rock
2 MPa
Shear strength of the natural fracture
1.5/1.3/1 MPa

Normal energy release rate of the rock
28 N/mm
Normal energy release rate of the natural fracture
26/24/22 N/mm

Tangential energy release rate of the rock
28 N/mm
Tangential energy release rate of the natural fracture
26/24/22 N/mm

Filtration coefficient

Approaching angle
30°/45°/60°/75°

Injection rate
5 m³/min·m
Viscosity
15 mPa•s

3. Result Analysis of the Three Models

3.1. The First Model: A Single Natural Fracture

The approaching angel is defined as the acute angle between the artificial fracture and natural fracture. According to the first model, the effect of a single natural fracture on artificial fracture propagation morphology with approaching angles of 30°, 45°, 60°, and 75° is studied, when the in situ stress level differences are 0 MPa, 5 MPa, and 10 MPa (shown in Figure 5).

Figure 5: Influence of the approaching angle and in situ horizontal stress difference on artificial fracture propagation morphology.

Under the same approaching angle, a single natural fracture has different effects on the expansion of an artificial fracture at different in situ horizontal stress difference levels. When the approaching angle reaches 30° and the in situ stress difference is less than 10 MPa, the artificial fracture extends along the right side of the natural fracture (shown in Figures 5(a)–5(c)). When the approaching angle is 45° and the in situ horizontal stress difference is less than 5 MPa, the artificial fracture propagates along the right wing of the natural fracture (shown in Figures 5(d) and 5(e)), while when the horizontal stress difference is 10 MPa, the artificial fractures cross the natural fractures (shown in Figure 5(f)). When the approaching angles are 60° and 75° and the value of the horizontal stress difference is 0 MPa, artificial fractures extend along the two wings of the natural fracture (shown in Figures 5(g) and 5(j)). As the in situ stress difference reaches 5 MPa, the artificial fracture is branched into two parts. One branch propagates along the natural fracture. The other branch crosses the natural fracture and extends along the direction of maximum in situ stress (shown in Figures 5(h) and 5(k)). As the in situ stress difference continues to increase to 10 MPa, the artificial fracture directly crosses through the natural fracture (shown in Figures 5(i) and 5(l)).

Under the same in situ horizontal stress difference, the effects of the single natural fracture on artificial fracture propagation morphology are different with different approaching angles. As the stress difference is 0 MPa and the approaching angles are 30° and 45°, the artificial fracture propagates along the right side of the natural fracture, while the left wing of the natural fracture is only partly opened (shown in Figures 5(a) and 5(d)). When the approaching angle increases to 60° and 75°, the left wing and right wing of the natural fracture are opened by fracturing fluid and the artificial fracture extends along the natural fractures on both wings (showed in Figures 5(g) and 5(j)). When the in situ horizontal stress difference is between 5 MPa and 10 MPa and the approaching angles are 30°, 45°, and 60°, the artificial fracture extends along the right wing of the natural fracture if the artificial fracture is arrested by the natural fracture (shown in Figures 5(b), 5(c), 5(e), and 5(h)), while artificial fracture extends along the two wings of the natural fracture with the approaching angle increasing to 75° (shown in Figure 5(k)).

Figure 6 shows that the in situ horizontal stress in the direction perpendicular to the natural fracture causes the natural fracture to be in the compressive stress state. The condition for opening the natural fracture is that the injected liquid pressure is larger than the sum of fracture cohesion and the normal stress. When the artificial fracture and natural fracture intersect at different approaching angles, the branch of the natural fracture in the OB direction is easier to be opened, while the branch of the natural fracture in the OA direction is relatively difficult to be opened. This is mainly affected by the additional resistance effect of the small angle. The fracturing fluid has the smallest flow resistance along the large angle (∠COB) of the OB direction. However, the flow direction of fracturing fluid along the small angle (∠COA) is opposite to the initial flow direction, which will produce additional resistance. Additional pressure is required for fracture initiation and propagation along the branch in the OA direction.

Figure 6: Additional resistance effect of a small angle.

3.2. The Second Model: Two Intersected Natural Fractures

The influence of two intersected natural fractures on artificial fracture propagation with approaching angles of 30°, 45°, 60°, and 75° is studied, when the in situ horizontal stress differences are 0 MPa, 5 MPa, and 10 MPa. The results are shown in Figures 7–10.

Figure 7: Approaching angle of 30°.

When the approaching angle is 30°, the effect of the two intersected natural fractures on the propagation of an artificial fracture with different horizontal stresses is shown in Figure 7. Phase 1 is the fracture propagation pattern when the artificial fracture intersects with the first natural fracture. The artificial fracture extends along the natural fracture when the in situ horizontal stress differences are 0 MPa, 5 MPa, and 10 MPa separately, and affected by the additional resistance effect of the small angle, the artificial fracture only extends along the right wing of the first natural fracture (shown in Figures 7(a), 7(d), and 7(g)). Phase 2 is the fracture propagation pattern when the artificial fracture intersects with the second natural fracture. The artificial fracture propagates along the large angle between the fracture propagation direction and the branch of the second natural fracture (shown in Figures 7(b), 7(e), and 7(h)). Phase 3 is the final fracture propagation pattern with in situ stress differences of 0 MPa, 5 MPa, and 10 MPa. The artificial fracture extends along the large angle direction when it intersects with the first natural fracture and the second natural fracture (shown in Figures 7(c), 7(f), and 7(i)).

Comparing with Figure 5(a)–5(c), it can be seen that the existence of the second natural fracture changes the propagation path of the artificial fracture along the first natural fracture. When the in situ horizontal stress difference is 0 MPa, the second natural fracture inhibits the propagation of the artificial fracture along the left wing of the first natural fracture.

Figure 8 shows the influence of the two intersected natural fractures on the artificial fracture extension pattern with different horizontal stresses when the approaching angle is 45°. When the in situ horizontal stress is less than 5 MPa, the artificial fracture deflects along the first natural fracture and the second fracture and extends along the direction of the large angle branch of the natural fracture (shown in Figures 8(c) and 8(f)). As the horizontal stress difference reaches 10 MPa, the artificial fracture passes through the first natural fracture directly and it is not affected by the second natural fracture (shown in Figure 8(i)).

Figure 8: Approaching angle of 45°.

The effect of two intersected natural cracks on artificial fracture propagation in different in situ horizontal stress differences with the approaching angle of 60° is shown in Figure 9. When the in situ stress level difference is 0 MPa, the artificial fracture extends along the two wings of the first natural fracture and it extends faster in the direction of the right wing (shown in Figures 9(a) and 9(b)). Then, the artificial fracture extends in the direction of the large angle between the second natural fracture and the artificial fracture (showed in Figure 9(c)). The artificial fracture only extends along the right wing of the first natural fracture in phase 1 when the in situ horizontal stress difference is 5 MPa (shown in Figure 9(d)), and the final pattern of fracture propagation is shown in Figure 9(f). Compared with the single natural fracture of the first model (Figure 5(h)), the existence of the second natural fracture restrains the artificial fracture from crossing the first natural fracture directly and makes it completely divert and extends along the natural fracture. As the in situ stress difference reaches 10 MPa, the artificial fracture passes the first natural fracture directly (shown in Figure 9(i)).

Figure 9: Approaching angle of 60°.

The effect of the two intersected natural fractures on artificial fracture propagation under different in situ stresses with the approaching angle of 75° is shown in Figure 10. When the in situ horizontal stress difference is 0 MPa, the artificial fracture extends along the two wings of the first natural fracture. After intersecting with the second natural fracture, one branch of artificial fracture crosses the second natural fractures and another branch of the artificial fracture extends along the second natural fracture (shown in Figure 10(c)). As the in situ horizontal stress difference is 5 MPa, the artificial fracture crosses the first natural fracture. And the two wings of the first natural fracture are partly cracked as well (shown in Figure 10(f)). When the formation stress difference is 10 MPa, the artificial fracture crosses the first natural fracture directly (shown in Figure 10(i)).

Figure 10: Approaching angle of 75°.

According to the results of the second model, the influences of the two intersected natural fractures on artificial fracture propagation with different approaching angles and in situ stress differences are summarized as shown in Table 2. Under the condition of a small in situ stress difference, the artificial fracture is easy to extend along the natural fracture, and when the in situ stress difference increases to a certain value, the artificial fracture will cross the natural fracture directly. With a small in situ stress difference, the smaller the approaching angle is, the easier that the artificial fracture extends along one branch of the natural fracture, resulting in a single-fracture pattern. With the increase of the approaching angle, the artificial fracture will extend along the two wings of the natural fracture, which increases the complexity of the fracture pattern. However, in the case of high in situ stress difference, artificial fractures tend to cross the first natural fracture directly with the increase of the approaching angle. And the second natural fracture intersected with the first one cannot be cracked, which is unbeneficial in forming a complex fracture network. In all, the smaller the in situ stress difference is and the larger the approaching angle is, the easier it is to form a relatively complex fracture network.

Table 2: The influence of two intersected natural fractures on artificial fracture propagation.

Approaching angle (°)
In situ stress difference () (MPa)
The 1st natural fracture
The 2nd natural fracture

30
0
Partly deflected
Partly deflected

5
Partly deflected
Partly deflected

10
Partly deflected
Partly deflected

45
0
Deflected
Partly deflected

5
Partly deflected
Partly deflected

10
Crossed

60
0
Deflected
Partly deflected

5
Partly deflected
Partly deflected

10
Crossed

75
0
Partly deflected

5

10
Crossed

The state of being partly deflected indicates that the artificial fracture can only propagate along one wing of the natural fracture. The state of being deflected indicates that the artificial fracture will propagate along the two wings of the natural fracture. The state of being crossed indicates that the artificial fracture will cross the natural fracture directly.

3.3. The Third Model: Multi-Intersected Natural Fractures

3.3.1. In Situ Horizontal Stress Difference

The third model is used to study the propagation of an artificial fracture under the effect of multi-intersected natural fractures at reservoir scales. The in situ horizontal stress differences are 0 MPa, 5 MPa, and 10 MPa, and the other parameters are based on Table 1. The results are shown in Figure 11.

(a) 0 MPa

(b) 5 MPa

(c) 10 MPa

Figure 11: The effects of in situ horizontal stress difference on fracture propagation.

When the in situ pressure level difference is 0 MPa, the artificial fracture almost extends along the direction of the natural fracture, forming a main fracture that is similar with the strike of the natural fracture. At the same time, the artificial fracture is arrested with a part of the natural fractures, forming some branch fractures. As the in situ horizontal stress difference is 5 MPa, the propagation morphology of the artificial fracture is influenced by the in situ horizontal stress difference and natural fracture synchronously. The artificial fracture extends along the natural fracture at a part of intersections, but the main propagation path of the artificial fracture is deflected to the direction of the maximum in situ horizontal stress. As the in situ horizontal stress difference reaches 10 MPa, the artificial fracture extends along the maximum in situ horizontal stress unaffected by the natural fracture.

3.3.2. Cementation Strength of the Natural Fracture

According to Table 1, by changing the tensile strength, shear strength, and energy release rate of the natural fracture, three cementation strengths of the natural fracture (strong cementation, medium cementation, and weak cementation) are defined. The influence of the different natural fracture cementation strengths (strong cementation, medium cementation, and weak cementation) on artificial fracture propagation morphology is analyzed.

When the in situ stress differences are 0 MPa, 5 MPa, and 10 MPa, the relationship between the aperture and the length of the artificial fracture with three different natural fracture cementation strengths is showed in Figure 12. The length and aperture of the artificial fracture in the rock matrix increase with the increase of natural fracture cementation strength under the condition of different in situ horizontal stress differences. The length of the artificial fracture in the rock matrix decreases with the raise of in situ stress difference, while the aperture of artificial fracture in the rock matrix increases. This is mainly because the weaker the cementation strength of natural fracture is, the more easily the artificial fracture extends along the natural fracture. Fracturing fluid is more likely to enter a weak cemented natural fracture, resulting in an increase in the length of fractured natural fractures and the decrease in artificial fracture length in the rock matrix.

(a) 0 MPa

(b) 5 MPa

(c) 10 MPa

Figure 12: The influence of cementation strength of a natural fracture on artificial fracture propagation.

Figure 13 also confirms that the length of the fractured natural fracture increases with the decrease of the cementation strength of the natural fracture. And the artificial fracture is more likely to cross the natural fractures directly with the raise of in situ stress difference, resulting in a significant reduction in the length of the fractured natural fracture.

Figure 13: The influence of cementation strength of a natural fracture on natural fracture propagation.

4. Model Verification

In order to affirm the accuracy of the PFS-Cohesive method for numerical simulation of fracture intersection, triaxial experiments are carried out with different approaching angles and in situ horizontal stress differences. Figure 14 shows the influence of the approaching angles of 45°, 60°, and 90° on the artificial fracture propagation morphology when the in situ horizontal stress difference is 3 MPa. The pattern of the artificial fracture intersecting with the natural fracture is deflected, partly deflected, crossed, and partly deflected with the approaching angles of 45°, 60°, and 90°.

Figure 14: The influence of the approaching angle on artificial fracture propagation.

Many other scholars also studied the propagation morphology of an artificial fracture intersecting with a single natural fracture [32–34]. The experimental results of Zhou et al. [32] are most representative and are widely used in numerical simulation verification. Zhou et al. also used the triaxial hydraulic fracturing experimental test analysis of an artificial fracture propagation pattern intersecting with a single natural fracture with horizontal stress differences of 3 MPa, 5 MPa, 7 MPa, and 10 MPa and with approaching angles of 30°, 60°, and 90°. The experimental results of Zhou et al. [32] and us are summarized, and a chart (shown in Figure 15) is proposed, which is used to compare with the numerical simulation results of the first model to affirm the accuracy of this novel PFS-Cohesive method for the simulation of fracture intersection.

Figure 15: Result comparison between numerical simulation and physical simulation.

Figure 15 shows that the results of numerical simulation are highly similar with the experimental results, except for the abnormal point corresponding to the approaching angle of 30° and in situ horizontal stress difference of 10 MPa. It proves the accuracy of the PFS-Cohesive finite element method for simulating fracture intersection.

Wang et al. [35] used the extended finite element method (XFEM) to simulate the intersection of an artificial fracture with a single natural fracture. The same criterion on the basis of the energy release rate method is also adopted in his model to determine the cross/arrest behavior between the artificial fracture and natural fracture. The effect of a single natural fracture on artificial fracture propagation morphology with approaching angles of 90°, 75°, 60°, and 45° is studied, when the in situ stress level difference is between 0.69 MPa and 10.35 MPa. The simulation results of the extended finite element method (Wang et al. [35]) and FPS-Cohesive finite element method are compared as shown in Figure 16. The comparison results show that the PFS-Cohesive results are highly similar with the XFEM results, which further proves the accuracy of the PFS-Cohesive finite element method for simulating fracture intersection.

Figure 16: Result comparison between PFS-Cohesive method and XFEM method.

5. Conclusions

(1)The dual-attribute element of a pore fluid/stress element and cohesive element (PFS-Cohesive) method for a rock matrix is presented based on the further development of a cohesive finite element. It solves the limitation of the ABAQUS cohesive finite element method in simulating hydraulic fracturing requiring a preset fracture propagation path. This new method can realize the simulation of artificial fracture propagation along any arbitrary path. And the accuracy of this method is verified by experimental results and XFEM method results(2)The in situ horizontal stress difference, approaching angle, and cementation of the natural fracture have significant effects on the propagation morphology of an artificial fracture. In general, when the artificial fracture intersects with the natural fracture, a smaller in situ stress difference, a greater approaching angle, and a weaker cementation strength of the natural fracture are favorable to form a relatively complex fracture network(3)When two intersected natural fractures exist, the second natural fractures will inhibit an artificial fracture extending along the direction of a small angle branch of the first natural fracture. After artificial fracture propagating along the direction of the large angle branch of the first natural fracture, it is easier to continue to propagate along the direction of the large angle branch of the second natural fracture, resulting in only one wing of the first natural fracture being cracked(4)When the in situ horizontal stress difference is 0 MPa, the artificial fracture almost extends along the natural fracture, forming a main fracture consistent with the strike of the natural fracture. The artificial fracture tends to deflect to the maximum in situ horizontal stress direction with the raise of in situ horizontal stress difference. When the in situ horizontal stress difference is 10 MPa, the artificial fracture tends to extend along the direction of the maximum in situ horizontal stress and it is almost unaffected by the natural fracture(5)The length and aperture of an artificial fracture in a rock matrix increases with the increase of natural fracture cementation strength. And with the rising of in situ horizontal stress difference, the artificial fracture length in a rock matrix decreases, while the aperture of the artificial fracture increases

Data Availability

The data used to support the findings of this study are available from the corresponding author upon request.

Conflicts of Interest

We have no competing interests.

Acknowledgments

The authors would like to acknowledge the financial support of the National Natural Science Foundation of China (Grant no. 51874338) and express their gratitude to the Fundamental Research Funds for the Central Universities (Grant no. 17CX02077), the Applied Basic Research Project of Qingdao Province (Grant no. 17-1-1-20-jch), and the Innovation Funding Program of the China University of Petroleum (East China) (YCX2018010).

References

S. Jiang, X. L. Tang, S. X. Long et al., “Reservoir quality, gas accumulation and completion quality assessment of Silurian Longmaxi marine shale gas play in the Sichuan Basin, China,” Journal of Natural Gas Science and Engineering, vol. 39, pp. 203–215, 2017.

S. C. Zhang, T. K. Guo, T. Zhou, Y. S. Zou, and S. R. Mu, “Fracture propagation mechanism experiment of hydraulic fracturing in natural shale,” Acta Petrolei Sinica, vol. 35, no. 3, pp. 496–503, 2014.

S. Heng, C. H. Yang, Y. T. Guo, C. Y. Wang, and L. Wang, “Influence of bedding planes on hydraulic fracture propagation in shale formations,” Chinese Journal of Rock Mechanics & Engineering, vol. 34, no. 2, pp. 228–237, 2015.

X. S. Shang, Y. H. Ding, L. F. Yang, Y. J. Lu, X. M. Yan, and Y. H. Wang, “Network fracturing technique in shale reservoirs based on weak discontinuity and fracture intersection,” Natural Gas Geoscience, vol. 27, no. 10, pp. 1883–1891, 2016.

C. H. Ahn, R. Dilmore, and J. Y. Wang, “Modeling of hydraulic fracture propagation in shale gas reservoirs: a three-dimensional, two-phase model,” Journal of Energy Resources Technology, vol. 139, no. 1, article 012903, p. 13, 2017.

D. S. Zhou, P. Zheng, J. Peng, and P. He, “Induced stress and interaction of fractures during hydraulic fracturing in shale formation,” Journal of Energy Resources Technology, vol. 137, no. 6, article 062902, p. 6, 2015.

M. M. Rahman, “A fully coupled numerical poroelastic model to investigate interaction between induced hydraulic fracture and pre-existing natural fracture in a naturally fractured reservoir: potential application in tight gas and geothermal reservoirs,” in SPE Annual Technical Conference and Exhibition, New Orleans, LA, USA, October 2009.

J. E. Olson, B. Bahorich, and J. Holder, “Examining hydraulic fracture: natural fracture interaction in hydrostone block experiments,” in SPE Hydraulic Fracturing Technology Conference, The Woodlands, TX, USA, February 2012.

R. Keshavarzi and S. Mohammadi, “A new approach for numerical modeling of hydraulic fracture propagation in naturally fractured reservoirs,” in SPE/EAGE European Unconventional Resources Conference & Exhibition-From Potential to Production, Vienna, Austria, March 2012.

G. Ma, F. Zhang, X. Liu, and D. Feng, “Research on the influence of natural fractgures on hydraulic fractures propagation in coal and rock mass,” Journal of Henan University (Natural Science), vol. 35, no. 2, pp. 178–182, 2016.

J. Xie, H. Huang, H. Ma et al., “Numerical investigation of effect of natural fractures on hydraulic-fracture propagation in unconventional reservoirs,” Journal of Natural Gas Science and Engineering, vol. 54, pp. 143–153, 2018.

A. Dahi-Taleghani and J. E. Olson, “Numerical modeling of multistranded hydraulic fracture propagation: accounting for the interaction between induced and natural fractures,” SPE Journal, vol. 16, no. 3, pp. 575–581, 2011.

R. Keshavarzi, S. Mohammadi, and H. Bayesteh, “Hydraulic fracture propagation in unconventional reservoirs: the role of natural fracture,” in 46th US Rock Mechanics/Geomechanics Symposium, Chicago, IL, USA, June 2012.

S. C. Maxwell, C. Waltman, N. R. Warpinski, M. J. Mayerhofer, and N. Boroumand, “Imaging seismic deformation induced by hydraulic fracture complexity,” SPE Reservoir Evaluation & Engineering, vol. 12, no. 1, pp. 48–52, 2009.

J. Zhao, Y. M. Li, S. Wang, Y. S. Jiang, and L. H. Zhang, “Simulation of a complex fracture influenced by natural fractures,” Natural Gas Industry, vol. 34, no. 1, pp. 68–73, 2014.

B. Hou, M. Chen, Z. Li, Y. Wang, and C. Diao, “Propagation area evaluation of hydraulic fracture networks in shale gas reservoirs,” Petroleum Exploration and Development, vol. 41, no. 6, pp. 833–838, 2014.

P. Wen and T. Fan, “The discontinuity displacement method applied to three-dimensional coplanar crack problem for any boundary-value condition,” Engineering Fracture Mechanics, vol. 48, no. 5, pp. 691–702, 1994.

N. Moës, J. Dolbow, and T. Belytschko, “A finite element method for crack growth without remeshing,” International Journal for Numerical Methods in Engineering, vol. 46, no. 1, pp. 131–150, 1999.

T. Wang, Z. L. Liu, and Z. Zhuo, “Numerical simulation of interaction between hydraulic fractures and natural fractures,” in The 24th annual conference of Beijing mechanics Society, pp. 754–756, Beijing, 2017.

J. Guo, X. Zhao, H. Zhu, X. Zhang, and R. Pan, “Numerical simulation of interaction of hydraulic fracture and natural fracture based on the cohesive zone finite element method,” Journal of Natural Gas Science and Engineering, vol. 25, no. 4, pp. 180–188, 2015.

M. Gonzalez-Chavez, A. Dahi Taleghani, and J. E. Olson, “A cohesive model for modeling hydraulic fractures in naturally fractured formations,” in SPE Hydraulic Fracturing Technology Conference, The Woodlands, TX, USA, February 2015.

M. Haddad, J. Du, and S. Vidal-Gilbert, “Integration of dynamic microseismic data with a true 3D modeling of hydraulic fracture propagation in vaca muerta shale,” in SPE Hydraulic Fracturing Technology Conference, The Woodlands, TX, USA, February 2016.

W. L. Xu, S. F. Xue, X. Zhang, C. Li, and F. Sun, “Numerical simulation study on propagation regularity of hydraulic fracture in natural fractured reservoirs,” Contemporary Chemical Industry, vol. 46, no. 7, pp. 1371–1374, 2017.

R. Zhang, G. S. Li, and H. Y. Zhu, “Dynamic propagation of multiple horizontal fractures and mutual interference between induced stresses,” Journal of Southwest Petroleum University (Science & Technology Edition), vol. 39, no. 1, pp. 91–99, 2017.

B. Carrier and S. Granet, “Numerical modeling of hydraulic fracture problem in permeable medium using cohesive zone model,” Engineering Fracture Mechanics, vol. 79, pp. 312–328, 2012.

S. Salimzadeh, A. Paluszny, and R. W. Zimmerman, “Three-dimensional poroelastic effects during hydraulic fracturing in permeable rocks,” International Journal of Solids and Structures, vol. 108, pp. 153–163, 2017.

M. A. Biot, “General theory of three-dimensional consolidation,” Journal of Applied Physics, vol. 12, no. 2, pp. 155–164, 1941.

O. C. Zienkiewicz and R. L. Taylor, The Finite Element Method: An Introduction with Differential Equations, Elsevier, Burlington, 2005.

T. J. Boone and A. R. Ingraffea, “A numerical procedure for simulation of hydraulically driven fracture propagation in poroelastic media,” International Journal for Numerical and Analytical Methods in Geomechanics, vol. 14, no. 1, pp. 27–47, 1990.

H. Y. Wang, “Poro-elasto-plastic modeling of complex hydraulic fracture propagation: simultaneous multi-fracturing and producing well interference,” Acta Mechanica, vol. 227, no. 2, pp. 507–525, 2016.

M. L. Benzeggagh and M. Kenane, “Measurement of mixed-mode delamination fracture toughness of unidirectional glass/epoxy composites with mixed-mode bending apparatus,” Composites Science and Technology, vol. 56, no. 4, pp. 439–449, 1996.

J. Zhou, M. Chen, Y. Jin, and G. Q. Zhang, “Analysis of fracture propagation behavior and fracture geometry using a tri-axial fracturing system in naturally fractured reservoirs,” International Journal of Rock Mechanics and Mining Sciences, vol. 45, no. 7, pp. 1143–1152, 2008.

F. Yao, M. Chen, X. D. Wu, and G. Q. Zhang, “Physical simulation of hydraulic fracture propagation in natural fractured formations,” Oil Drilling & Production Technology, vol. 30, no. 3, pp. 83–86, 2008.

W. Cheng, Y. Jin, M. Chen, T. Xu, Y. Zhang, and C. Diao, “A criterion for identifying hydraulic fractures crossing natural fractures in 3D space,” Petroleum Exploration and Development, vol. 41, no. 3, pp. 371–376, 2014.

X. L. Wang, F. Shi, C. Liu, D. T. Lu, H. Liu, and H. A. Wu, “Extended finite element simulation of fracture network propagation in formation containing frictional and cemented natural fractures,” Journal of Natural Gas Science and Engineering, vol. 50, no. 1, pp. 309–324, 2018.