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-15 μm2), is pressure difference (Pa), and is liquid density (kg/m3).
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).
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).
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.
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
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).
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.
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.
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)).
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)).
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)).
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.
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
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 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.
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°.
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 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.
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).