Abstract
Aiming at the problem of quantifying the composition of complex fracture network and fracture activation after fracturing, a complex fracture network model of weak plane with prefabricated structure is established based on finite element method, global embedded cohesive zone model (CZM), and real shale outcrop. Considering the influence of fully coupled stress/fluid, the effects of weak plane azimuth, horizontal stress difference, fracturing fluid viscosity, and injection rate on fracture network composition, geometry, and stimulated reservoir volume (SRV) are studied. The concept of fracture relative activation rate, which can quantitatively analyze fracture network composition and fracture activation, is proposed. The results show that the fracture geometry of the two kinds of conjugated shale after fracturing is controlled by the most weak mechanical plane, and the fracture network is, respectively, axisymmetric and centrosymmetric. The fracture network is composed of the weak plane fractures with the dominant free gas transport and the matrix microfractures with the dominant adsorbed gas transport. The effects of weak plane azimuth, horizontal stress difference, and fracturing fluid viscosity on SRV length are not monotonous, while the increase of the azimuth, horizontal stress difference, fracturing fluid viscosity, and the reduction of appropriate injection rate will lead to the increase of SRV width. The effect of horizontal stress difference on the relative activation rate of fractures does not change monotonically, while the increase of weak plane azimuth, fracturing fluid viscosity, and injection rate will lead to the increase of the relative activation rate of matrix microfractures, the increase of the total length of activated fractures, and the decrease of the relative activation rate of weak plane fractures.
1. Introduction
Multistage and cluster fracturing of horizontal wells is a key technology for the United States to realize the shale gas revolution and change the global energy trade pattern. Although this technology has been adopted on a large scale in China, there are still some unsolved problems in the composition of complex fracture network after fracturing, fine analysis of the activation rate of different types of fractures, and the description of fracture geometry [1]. Shale has structural anisotropy due to a large number of local faults, bedding planes, and natural fractures (NFs), as shown in Figures 1(a) and 1(b). The developed structural weak plane will have a great impact on the initiation, propagation, and geometric shape of hydraulic fractures (HFs) in shale. Therefore, it is of great significance to study the composition and complex geometry of hydraulic fractures in shale reservoir under different structural weak plane distribution for the efficient development of shale reservoir [2, 3].
Affected by bedding and natural fracture development, the geometric shape of hydraulic fracture in shale is quite different from that in conventional reservoir. The hydraulic fractures in the conventional reservoir are symmetrically distributed with double wings, while the hydraulic fractures in the shale reservoir can expand and turn along the matrix, bedding, and natural fractures under the coupling action of hydraulics and mechanics, showing a complex and unpredictable fracture network [4]. At present, the numerical simulation methods for the propagation mechanism of complex fracture networks mainly include displacement discontinuous method, finite element method, extended finite element method, discrete element method, cohesive element method, and phase field method. The cohesive element method is widely used in civil engineering, mechanical engineering, and material science [5]. It was Barenblatt who first proposed the concept of cohesive zone, and he first used cohesive element method to simulate crack propagation in brittle materials [6]. Mokryakov [7] successively simulated the fracturing process in soft rock by using the cohesive element method and the hydraulic fracturing expansion behavior considering the influence of fracture tip and plastic inside rock. Guo et al. [8] established the intersection model of hydraulic fracture and natural fracture based on cohesive element method and discussed the influence of different factors on fracture intersection. Guo et al. [9] used this method to simulate the extension behavior of hydraulic fractures in layered shale. With the development of global embedded cohesive element technology, a new cohesive zone model is gradually formed based on viscosity element method to simulate cracks with arbitrary paths. Wang [10] established a hydraulic fracture extension model considering fluid structure coupling by using extended finite element and cohesive zone model and studied the important effects of inelastic deformation on propagation pressure and fracture geometry. Dahi Taleghani et al. [11] first obtained the main parameters to ensure the accuracy of traction separation criterion through indoor experiments and then simulated the intersection process of hydraulic fracture and natural fracture using cohesive zone model. Wang [12] also established a global embedded cohesive unit model to simulate complex fracture propagation based on the cohesive zone model and considered that under the influence of natural fracture density, size, distribution, stress, and construction parameters, only complex fractures rather than connected fracture networks may be formed in the fractured reservoir. Yu et al. [13] established a finite element model of coupled fluid flow and deformation based on the adaptive embedded cohesive zone model and simulated the propagation behavior of fractures in rock matrix and cross natural fracture network. Li [14] analyzed the influence of fracturing temporary plugging technology on fracture complexity based on the modified cohesive zone model. Compared with other numerical simulation methods, the cohesive zone model can eliminate the singularity at the crack tip and overcome the disadvantage of nonconvergence of linear elastic fracture mechanics. For the reservoir with weak plane of shale structure, using the cohesive zone model to globally embed the cohesive unit in the weak plane and matrix, and set the corresponding mechanical parameters, respectively, the propagation dynamics of hydraulic fractures in the weak plane structure and matrix can be accurately simulated.
Based on the real shale outcrop with two different weak plane characteristics and the global embedded cohesive zone model, a shale model with ideal prefabricated weak plane structure is established. Using this model, the geometric shape, fracture network composition, and relative activation of hydraulic fractures in structurally anisotropic shale are studied. Figure 1(c) shows the shale model with conjugate cross distribution of weak plane, which is composed of weak plane cohesive unit (Figure 1(f)), matrix cohesive unit (Figure 1(e)), and matrix finite element (Figure 1(d)). Although the global embedding of cohesive element will increase the overall calculation cost, it can more accurately capture the intersection, passage, blocking, and capture of hydraulic fractures in matrix and weak plane structures. Through the simulation of the activation behavior of weak plane fractures and matrix fractures and the study of their internal mechanism, the research results of this paper have certain guidance and reference value for optimizing hydraulic fracturing design, fine evaluation of postfracturing reconstruction effect, and productivity prediction.
2. Mathematical Model
2.1. Fracture Propagation Model
Figure 2(a) shows the normal and tangential fluid flows in the fracture. Tangential flow promotes fracture propagation, and normal flow represents fracturing fluid filtration into the formation. According to Newtonian fluid rheology theory, fluid is considered incompressible. The tangential flow in the crack is controlled by the lubrication equation, which is derived from Poiseuille’s law [15]: where is the volume flow through the fracture section, is the fracture width, is the viscosity of fracturing fluid, and represents the fluid pressure gradient along the fracture direction.
The fluid normal filtration loss in the fracture can be expressed as [16] where and are the filtration coefficient at crack top and bottom; and , respectively, represent pore pressures in matrix units adjacent to cohesive units (as shown in Figure 2(b)); and is the filtration loss of fracturing fluid.
According to the principle of mass conservation, the mass conservation equation in cracks can be expressed as where represents the fracturing fluid injection rate and is the Dirac delta function.
The fracture opening is determined by the properties of cohesive elements, fluid properties, pore pressure, fluid pressure in the fracture, stress distribution, and damage criteria and can be determined by the displacement of the top and bottom of the fracture wall [17]: where and , respectively, represent the displacements on the top wall and bottom wall of cracks and is the unit normal vector on crack wall.
2.2. Law of Traction Separation
In this paper, the secondary stress failure criterion is used to judge the initiation and propagation of cracks. The secondary stress failure criterion can determine the initiation and propagation of cracks by three stress components and the peak values of three stress components. The criterion can be expressed as where , , and represent the normal stress and first and second shear stress components, respectively, and , , and , respectively, represent the tensile strength and shear strength of the complete cohesive unit.
The stress component of the traction separation law is also affected by the damage variable when the bond element fails [9]: where is the nominal stress vector and is the stress component of undamaged displacement predicted by elastic traction separation law.
As shown in Figure 3, bilinear cohesion law is adopted in this paper to describe the relationship between traction force and displacement [18]. In the above formula, dimensionless represents the damage degree of cohesive unit during fracturing, and its value can be determined according to the traction separation law [12]: where and represent the displacement component when the traction of cohesive unit reaches its limit and and , respectively, are defined as where , , and , respectively, represent the normal displacement component and the first and second shear displacement component; indicates the tensile strength or shear strength of the material; and is the fracture energy of mixed mode. Based on Benzeggagh-Kenane fracture criterion [11], damage evolution of cohesive units in the process of fracture expansion can be determined, which can be defined as where , , and represent the work done by the traction force in three directions and , , and represent the fracture energy in normal direction and two shear directions, respectively.
2.3. Fluid Flow and Geomechanical Coupling Model
Based on the volumetric virtual work principle, the stress balance equation can be written as [19–21] where is the effective stress, indicates pore pressure, is the virtual strain rate matrix, represents virtual velocity, and represent the surface traction force per unit area and the physical force per unit volume, respectively, and represents the identity matrix.
The continuity equation of fluid seepage in reservoir can be expressed as [22, 23] where is the volume change rate of porous medium, indicates fluid density, and indicates the porosity of the reservoir matrix.
3. Advantages of the Model
As shown in Figure 4(a), it is the grid after the global embedding of the cohesive zone model. One-dimensional cohesive elements are embedded between all four node grid elements (the red-dotted line between grid elements indicates cohesive elements). The embedded global cohesive unit forms the potential extension path of the crack. This method is similar to the finite element discrete element method used by Zhang et al. [24] and Blanton and Q. Wang et al. [25, 26] to simulate hydraulic fracturing behavior. Figure 4(b) shows the conventional cohesive element model. Under this condition, the crack can only extend along the set path. Compared with the conventional cohesive element model, the model in this paper can capture the extension direction of cracks more accurately. Figure 4(c) shows the composition of the fluid pressure node at the intersection in the global embedded cohesive element. When cohesive elements intersect, four cohesive elements (CE1~4) use a fluid pressure node at the intersection. After the secondary development with ABAQUS, the hydraulic fracturing model is established by using the finite element method, and then, the cohesive element is embedded globally to form the potential multidegrees of freedom of fracture propagation.
4. Results Discussion and Analysis
4.1. Shale Model with Conjugate Distribution of Structural Weak Plane
As shown in Figures 1(a) and 1(b), the structural weak planes in the two shale outcrops are relatively developed, and the horizontal weak plane and longitudinal weak plane are staggered with each other at a certain angle in a conjugate distribution form. According to the experiment, the mechanical properties of the two weak planes in the shale outcrop in Figure 1(a) are similar, while the mechanical properties of the horizontal weak plane and the longitudinal weak plane in the shale outcrop in Figure 1(b) are different. Based on the structural distribution characteristics and mechanical characteristics of the weak plane of the two shale outcrops, two prefabricated ideal distribution weak plane models are established, as shown in Figure 5. Figure 5(a) is a weak plane model based on the shale outcrop characteristics in Figure 1(a). All weak planes in the model have the same mechanical properties. Figure 5(b) is a conjugate weak plane model composed of weak plane I and weak plane II based on the characteristics of shale outcrop in Figure 1(b), in which weak plane I and weak plane II have different mechanical properties. The approaching angle is the angle between the weak plane direction and the horizontal maximum principal stress direction. The boundary condition of the model is that the displacement along the direction perpendicular to the boundary is fixed at 0, and the pore pressure on the boundary is set as a constant, as shown in Figure 5(a). The main input parameters in the simulation are shown in Table 1.
(a) Model- I
(b) Model- II
4.2. Fracture Geometry Based on Model I and Model II
Based on model I, under the conditions of an approximation angle of 45°, a stress difference of 5 MPa, and an injection rate of 0.1 m3/s, the geometric forms of hydraulic fractures in conjugated weak plane shale at different times were simulated, as shown in Figure 6. The parameters of weak plane cohesive unit in simulation are the same as those of weak plane I cohesive unit in Table 2. In order to accurately observe the geometric shape of the crack, the opening distribution of the damaged cohesive unit after concealing the matrix unit is shown in Figure 7, which is the real hydraulic fracture network formed after compression. The color band values in Figures 6 and 7 and all subsequent figures represent the crack opening (PFOEN), and the unit is m. It can be seen from Figure 6 that under the above initial stress distribution and construction parameters, the main hydraulic fracture formed after fracturing mainly expands along the weak plane structure with weak mechanical properties. Figure 7 shows that during fracturing, while the hydraulic fracture expands along the weak plane structure, it also compresses some cohesive units in the matrix block near the weak plane to form matrix microfractures. Therefore, the real fracture network formed by fracturing in shale with weak plane structure is composed of weak plane fractures and matrix microfractures, and the fracture geometry is in the shape of axisymmetric network. From the perspective of production, the activated matrix type microfractures can become a new flow channel of shale gas existing on the surface of kerogen and pores in the form of adsorption, while the activated weak plane type fractures become the main flow channel for the accumulation and transmission of all desorbed and free shale gas. After stimulation and reconstruction, the proportion of two types of activated fractures will have an important impact on later production.
(a) 2 s
(b) 5 s
(c) 10 s
(a) 2 s
(b) 5 s
(c) 10 s
Based on model II, under the conditions of approximation angle of 45°, stress difference of 5 MPa, and injection rate of 0.1 m3/s, the distribution of hydraulic main fracture opening at different times was simulated, as shown in Figure 6. Figure 8 shows the opening distribution of weak plane fracture and matrix microfracture. In the simulation, the main mechanical parameters of matrix cohesive unit and prefabricated weak plane I and weak plane II cohesive units are shown in Table 2. It can be seen that the shear and tensile resistance of cohesive unit in weak plane II is relatively weak; that is, it is easier to be destroyed and form hydraulic fractures. The analysis of Figures 9 and 8 shows that, under the influence of the mechanical properties of the weak plane, the fracture network tends to expand along the weak plane II with weak mechanical properties and finally forms a center-symmetric fracture network composed of weak plane type fractures and matrix type microfractures. Among the factors affecting the mesh geometry, the structural weak plane with the weakest mechanical properties will play a dominant role.
(a) 2 s
(b) 5 s
(c) 10 s
(a) 2 s
(b) 5 s
(c) 10 s
4.3. Multivariate Analysis Based on Model I
4.3.1. Influence of Approaching Angle on Fracture Geometry
In order to quantitatively analyze the impact of the approach angle on the shape of the fracture network and the fracture activation rate, this paper uses the changes of the major axis and minor axis of the postfracturing reservoir reconstruction volume (SRV) (it is considered that the SRV after the actual reconstruction is an ellipse) to analyze the shape of the formed fracture network (as shown in Figure 10(b)). Then, the relative activation rate of fracture is introduced to characterize the relative activation rate of hydraulic fracture in weak surface structure and matrix. The relative activation rate of cracks can be expressed by the ratio of the activated crack length in the weak plane structure or matrix to the total activated crack length: where and indicate the relative activation degree of weak surface type and matrix type microcracks, respectively; and represents the total number of activated cohesive units in the weak plane structure and matrix, respectively; and denote the number of the weak plane structure and the activated cohesive unit in the matrix, respectively; and and represent the length of the and cohesive units in the weak plane and matrix, respectively, m.
(a)
(b)
Figures 11(a) and 11(b), respectively, show the influence of approaching angle on the relative activation of SRV long axis, short axis, weak plane type fractures, and matrix type microcracks. As can be seen from Figure 11, with the increase of approaching angle, the short axis of SRV gradually increases, while the long axis of SRV first decreases and then increases. The reason why the long axis of SRV increases when the approaching angle is 75° is that under the influence of stress difference, after exceeding a certain degree, the failure conditions of weak plane cohesive unit and matrix cohesive unit are similar, and hydraulic fractures can expand along weak plane structure or matrix to form main fractures. With the increase of approaching angle, the total length of fractures increases first and then decreases, while the relative activation rates of weak plane fractures decreases gradually, while the relative activation rates of matrix microfractures increases gradually. The results show that there is an optimal approaching angle to maximize the total length of reservoir fractures under different approaching angles, and the increase of approaching angle can increase the activation rates and complexity of matrix microfractures and reduce the activation rates of weak plane fractures.
(a)
(b)
4.3.2. Influence of Horizontal Stress Difference on Crack Propagation
Figures 12(a) and 12(b) show the effects of horizontal stress difference on the relative activation rates of SRV long axis, short axis, weak plane type, and matrix type microcracks, respectively. It can be seen from the figure that the fracture network formed under different stress differences is still composed of weak plane cracks and matrix microcracks, of which the main hydraulic cracks are composed of weak plane cracks. With the increase of horizontal stress difference, the short axis of SRV decreases gradually, and the length of SRV decreases first and then increases. The reason why the long axis of SRV decreases at 5 MPa is that when the horizontal stress difference increases, the weak plane cohesive unit is relatively more difficult to destroy, and more fluid is used to activate the matrix microcracks. As can be seen from Figure 12(b), with the increase of stress difference, the total length of activated fractures and the relative activation rates of matrix microfractures increase first and then decrease, while the relative activation rates of weak plane fractures decrease first and then increase. In general, the horizontal stress difference has an important effect on the geometry of fracture network and the relative activation rate of fracture. Small or large horizontal stress difference is not conducive to the formation of matrix microcracks, and small horizontal stress difference is more conducive to the formation of weak plane fractures. Under different horizontal stress differences, there exists an intermediate horizontal stress difference value to maximize the total length of activated fracture.
(a)
(b)
4.3.3. Influence of Fracturing Fluid Viscosity on Fracture Propagation
The influence of viscosity on length and width of SRV and fracture relative activation rate is shown in Figure 13. The results show that the fracture fluid viscosity has a significant effect on the composition and geometry of fracture network. With the increase of fracturing fluid viscosity, the width of SRV and the relative activation rate of weak plane fracture gradually decrease, while the length of SRV firstly decreases and then increases. The opening of weak plane fracture, the total length of active fracture, and the relative activation rate of matrix microfracture gradually increase. These results indicate that the change of fracturing fluid viscosity can significantly change the ratio of weak plane type fractures to matrix type microfractures and affect the expected productivity. The fracturing fluid with smaller viscosity has stronger deep penetration ability, which can make the fracturing sweep range larger. The fracturing fluid with higher viscosity has stronger reconstruction ability near the wellbore. The shale matrix near the injection point is fully broken by forming matrix microfractures to form a channel for the outflow of adsorbed gas. The increase in fracturing fluid viscosity essentially reduces the efficiency of flow, increases the resistance of fluid flow, and greatly increases the fluid pressure in the fracture near wellbore area. The increase of fluid pressure in the fracture will produce stronger stress shadowing effect, resulting in more shear-dominated matrix microfractures at the end of the main fracture, which is also the reason why the length of SRV increases at 30 mPa.s.
(a)
(b)
4.3.4. Effect of Injection Rate on Crack Propagation
Injection rate is also an important parameter to control the geometry of fracture network. Based on model I, under the condition of stress difference of 5 MPa and fracturing fluid viscosity of 1 mPa.s, the effects of different injection rate on SRV and fracture relative activation rate were simulated (as shown in Figure 10). The results show that with the increase of injection rate, the extension range, SRV length, and width of weak plane crack and the relative activation rate of weak plane crack gradually decrease, while the opening of weak plane crack, the total length of activated crack, and the relative activation rate of matrix microcrack gradually increase. Therefore, properly increasing the fracturing fluid injection rate and viscosity will help to activate more matrix microfractures and break the matrix blocks near the wellbore. Properly reducing the fracturing fluid injection rate and viscosity will help to activate more weak plane fractures and realize deep penetration by using the developed weak plane structure.
5. Conclusion
Based on the global embedded cohesive zone model and the characteristics of real shale outcrop, a hydraulic fracturing fracture network simulation model coupled with fluid/stress is established for the shale with conjugate cross weak plane structure distribution. The effects of weak plane approaching angle, horizontal stress difference, fracturing fluid viscosity, and injection rate on fracture network geometry and SRV are comprehensively studied. The concept of fracture relative activation rate was proposed for quantitative analysis of fracture network composition and fracture activation behavior. The results show the following: (1) the fracture geometry after fracturing of shale is dominated by the weak plane structure with the weakest mechanical properties, which shows an axially symmetric network and a centrosymmetric network, respectively. The fracture network consists of weak plane fracture and matrix microfracture. (2) The weak plane approaching angle, horizontal stress difference, and fracturing fluid viscosity do not show monotonic effects on SRV length. The increase of approaching angle, horizontal stress difference, fracturing fluid viscosity, and appropriate decrease of injection rate will lead to the increase of SRV width. (3) The horizontal stress difference has a nonmonotonic effect on the relative activation rate of fractures. The increase of weak plane approaching angle, fracturing fluid viscosity, and injection rate will lead to the gradual increase of the relative activation rate and total activated fracture length of matrix microfractures and the decrease of the relative activation rate of weak plane fractures. (4) The fundamental reason for the change of fracture geometry, total length of active fracture, and relative activation rate of fracture is the comprehensive influence of stress shadow effect, fluid and geological stress coupling effect, and cohesive unit damage criterion under different geological parameters and different construction parameters. The research in this paper has important reference significance for hydraulic fracturing design, postfracturing productivity simulation, and prediction of shale with weak plane structure.
Data Availability
No data were used to support this study.
Conflicts of Interest
The authors declare that they have no conflicts of interest.