Abstract
The formation of the fracture network in shale hydraulic fracturing is the key to the successful development of shale gas. In order to analyze the mechanism of hydraulic fracturing fracture propagation in cemented fractured formations, a numerical simulation about fracture behavior in cemented joints was conducted based firstly on the block discrete element. And the critical pressure of three fracture propagation modes under the intersection of hydraulic fracturing fracture and closed natural fracture is derived, and the parameter analysis is carried out by univariate analysis and the response surface method (RSM). The results show that at a low intersecting angle, hydraulic fractures will turn and move forward at the same time, forming intersecting fractures. At medium angles, the cracks only turn. At high angles, the crack will expand directly forward without turning. In conclusion, low-angle intersecting fractures are more likely to form complex fracture networks, followed by medium-angle intersecting fractures, and high-angle intersecting fractures have more difficulty in forming fracture networks. The research results have important theoretical guiding significance for the hydraulic fracturing design.
1. Introduction
The porosity and permeability of shale are extremely low. To achieve the effective development of shale gas, fractures are required to provide flow channels for gas. Therefore, hydraulic fracturing is applied to increase fractures in shale. Natural weakness planes contribute to the formation of a complex hydraulic fracture network. Fortunately, many shale outcrops, cores, and image logs show that there are weakness planes including joints, faults, bed-parallel fractures, early compacted fractures, and fractures associated with concretions in shale reservoir [1–3]. There is also some evidence of natural fractures as shown in Figure 1. Those outcrops from north China provide proofs of cemented fractures and natural fracture networks in rocks (Figures 1(a) and 1(b)), while the cores from the Changning-Weiyuan Shale provide proofs of open or cemented natural fractures in deep formation (Figures 1(c) and 1(d)). Those weakness planes including open fractures and cemented cracks are heterogeneous. They can reduce rock strength and may divert hydraulic fracture propagation. Therefore, the interaction between the weakness planes and hydraulic fractures is significant for the formation of a complex fracture network.
(a) Cemented fractures of outcrops
(b) Natural fractures network of outcrops
(c) Natural fractures of cores (depth 3635 m)
(d) Cemented fractures of cores (depth 1531 m)
The experimental study primarily focuses on investigating the behaviors of a hydraulic fracture intersecting with a natural fracture. Hydraulic fracturing tests at a laboratory scale have been carried out at different approximation angles and differential stresses, which were run for a condition where the hydraulic fracture opened existing fractures to determine if more constructive crossing interaction could be obtained. For hydraulic fractures, when the angle of approach angle was 30 degrees, existing fractures tend to open and prevent the induced fracture from crossing. At an angle of 60 and 90, hydraulic fractures tend to cross existing fractures when the differential stress was high enough. Based on the results, an elastic solution has been used as a basis for a hydraulic/natural fracture interaction criterion for fracture opening [4]. But the shear slip of the discontinuous surface was ignored in Blanton’s criterion which only considers the condition that the fractures opened when the fluid pressure overcame the normal stress on the discontinuous surfaces. Therefore, Warpinski and Teufel gave the critical pressure of the shear slip on the fracture surfaces [5]. Then, a series of experimental investigations designed to assess the conditions required for the crossing were presented. And a stress solution of the linear elastic fracture mechanics near a fracture tip was used to calculate the criterion of the hydraulic/natural fracture interaction [6]. Besides, other related experiments about the natural/hydraulic fracture interaction were performed with a certain purpose. A criterion has been developed for fractures crossing frictional interfaces with cohesion at nonorthogonal intersection angles. Laboratory experiments were conducted to validate the criterion [7]. Bahorich et al. performed 9 hydraulic fracturing experiments in gypsum cement blocks. And the cemented natural fractures were represented by plaster discontinuities [8]. And a series of semicircular bend tests on Marcellus Shale core samples containing calcite-filled natural fractures (veins) were performed to investigate the influence of weak planes on hydraulic fracture propagation [9].
Additional numerical analyses have also been performed to develop propagation criteria for determining propagation direction when a hydraulic fracture encounters natural fractures. For example, the discrete element method [1, 10, 11], extended finite element method [12–14], cohesive zoned method [1, 15], phase field method [16], and displacement discontinuity method [17] were implemented to perform the numerical simulation which analyzes the crossing or arresting behavior of hydraulic and natural fractures. The discrete element method is one of the effective methods for hydraulic fracturing. Five typical coal models were established to simulate hydraulic fracturing in the coal seam based on the two-dimensional particle flow code (PFC2D) by Wang et al. [10]. Hydraulic fracture propagation and proppant transport in a two-layer formation with stress drop were studied by Zhang and Dontsov [18].
Previous studies mainly focus on the influence of damaged natural fractures on hydraulic fractures. The jointed surface of the weakness planes is not cracked due to cementation, but its strength is weaker. When the fracturing fluid reaches the weakness plane, the initiation and propagation of induced fractures are likely along the weakness plane. Therefore, the interaction between the weakness plane and the hydraulic fracture has an important influence on fracture propagation.
The purpose of this paper is to comprehensively consider the expansion path and morphology of hydraulic fractures intersecting with natural closed fractures under the shear and tensile action of fracture surfaces. The expansion patterns of hydraulic fractures after they meet with intersecting closed fractures are of great significance for the formation of complex fracture networks after fracturing. In this paper, a numerical simulation about fracture behavior in cemented joints was conducted based firstly on block discrete element. Then, the failure modes of hydraulic fractures in fractured formation are discussed, and the failure conditions under different failure modes are given. Based on the different failure conditions of the intersection, the fracture morphology with two boundaries after fracture failure at different intersecting angles is discussed. Besides, the influence of parameters on fracture morphology was analyzed by univariate analysis and response surface method (RSM).
2. Numerical Model
2.1. Block Discrete Element Method [18, 19]
The block discrete element method is a three-dimensional numerical modeling for advanced geotechnical analysis. The block discrete element method simulates the response of discontinuous media (such as jointed rock or masonry bricks) that is subject to either static or dynamic loading. The discontinuous material is represented as an assemblage of discrete blocks. Models may contain a mix of rigid or deformable blocks. Deformable blocks are defined by a continuum mesh of finite-difference zones, with each zone behaving according to a prescribed linear or nonlinear stress-strain law. The relative motion of the discontinuities is also governed by linear or nonlinear force-displacement relations for movement in both the normal and shear directions. Joint models and properties can be assigned separately to individual, or sets of, discontinuities.
The block discrete element method is based on the distinct element method (DEM) for discontinuum modeling. It differs from particle-based methods in its ability to represent a zero-initial porosity condition, as well as interlocked irregular block shapes that provide resistance to block rotation (moments) after contact breakage [20]. The block discrete element method consists of the continuum and discrete model simulations. Continuum blocks need not be represented by particles. Therefore, the calculation amount of this method is less than that of particle-based methods. This method can deal with larger models with lower equipment requirements than particle-based methods.
2.1.1. The Constitutive Model of Joints
The thickness of the joint is far less than the scale of its plane, so its deformation characteristics are described by the stress-displacement relationship. The constitutive relation of the structural plane mainly studies the relationship between the stress of the structural plane and its normal deformation and tangential deformation. There are two stresses on the structure surface: normal stress and shear stress , and two corresponding displacements: normal displacement and tangent displacement . The stress-displacement relation matrix is expressed as where , which is the normal stiffness coefficient and represents the effect of normal displacement on normal stress; , which is the shear stiffness coefficient and represents the effect of shear displacement on shear stress; , which is the dilation stiffness coefficient and represents the effect of shear displacement on normal stress; and , which represents the effect of normal displacement on shear stress.
The Goodman element is adopted in this calculation; i.e., and are set as 0 and and are respectively used to describe the normal deformation and tangential deformation of the structural surface.
The Mohr-Coulomb equation for shear strength of the structure plane is where and and are the cohesive force and friction angle of the structural plane, respectively.
2.1.2. Contact Friction Joint Model
The contact friction joint model and the coulomb slip model are used in this calculation. It is assumed that the normal stress increment (positive pressure) and the shear vector increment between blocks are proportional to the normal displacement and the tangential displacement increment , respectively, in the elastic stage. The contact stiffness of the normal and tangential springs is and , respectively. Assuming that there is no tension in the joint and Coulomb’s law is satisfied, then the relation can be expressed as where and are the normal and tangential components of contact force , respectively; and are the normal and tangential stiffness of joints, and are the normal and tangential relative displacements of joints, respectively; and are the friction coefficient and cohesive force of joint material, respectively; and is the length of the contact surface.
2.1.3. Fluid Flow in Crack [10, 21]
It is assumed that the fracture surface is impermeable and fluid flows only within the fracture surface. The flow of fluid in the crack satisfies the cube law. Using the modified cube law, where is the fluid flow rate per crack width, is the equivalent hydraulic aperture, is the density of the fluid, is the fluid dynamic viscosity, is the hydraulic gradient, and is the hydraulic conductivity.
In the hydromechanics coupling of fractures, the influence of mechanical deformation on fracture permeability is mainly manifested as the change of fracture aperture. In the elastic stage, the crack aperture is expressed as an equation related to the effective stress.
In the plastic stage, considering the expansion caused by a joint slip, the crack opening can be expressed as where is initial fracture aperture, is effective stress, and is a coefficient about the influence of stiffness and is set as 1 in this study.
2.1.4. Motion Equation of Blocks
The deformation of blocks is nonnegligible. The block is divided into tetrahedral units. The vertices of a tetrahedral element are called grid difference points. The motion equation is established at each node as follows: where is the outside surface, is the mass on the grid point, is the acceleration of gravity, and is the resultant force of an external force applied to a node.
The node force is 0 in the balance state. Otherwise, there is an acceleration of the node according Newton’s second law.
For every time step, strain and rotation are related with displacement. Their forms are as follows:
The incremental method is applied in calculation. It is not limited to small strain problems. Therefore, the constitution model of deformation blocks takes incremental form. The equation is as follows: where and are Lamé coefficients, and is the Kronecker function.
2.2. Numerical Model Description
Based on the block discrete element theory, the propagation behavior of hydraulic fractures in rock mass with cemented joints is simulated by the following model. The model including two cutters and one vertical joint that can simulate the behavior of fracture crossing and fracture containment is shown in Figure 2. As shown in Figure 2(b), the intersection angles between vertical and two cutters are and . The size of the model is 20 m in 3 dimensions. A 3D model of rock with joints and elements is shown in Figure 2(c). The distance between the two cutter centers is 4 m. And the parameters are shown in Table 1.
(a)
(b)
(c)
In this simulation, the and are set with the same angle. 30°, 45°, and 90° are chosen to investigate the effect of the intersection angle on crack propagation. And the vertical joint is divided into 3 parts by two cutters. The fluid-driven crack is always contained in the middle part of the vertical joint. The simulation results are shown in the next section.
3. Results and Discussion
3.1. Fracture Propagation in Rocks with Cemented Joints
The time-dependent injection pressure was recorded. The result is shown in Figure 3. As shown in the figure, the pressure increases rapidly at the initial period. The crack was initialized when it reached maximum pressure at 56.4 MPa. Then, there is a vibration period after crack open. At last, the pressure tended to become smooth and steady. This is consistent with the trend of actual monitoring data on site, which means that the simulation results are suitable for reference research.
One of the main purposes of numerical simulation is to investigate the behavior of a fluid-driven crack in the rock with cemented joints. In other words, the shape of open fractures after fracking is one of the most important goals. The shape is described by the aperture of fractures in this study. The results are plotted in Figure 4. Figures 4(a)–4(c) is the spatial distribution of fractures with different angles between the vertical joint and two cutters. As shown in the pictures, with different intersection angles, the final shape is varied. According to a predecessor’s studies, it is easy to understand that the crack extends straight forward with a high intersection angle, while the crack is always captured by the meeting natural fracture but would not stretch forward when the intersection angle is large enough. However, at the situation that two angles are 30°, it penetrates again. Besides, it is also captured by natural fractures. Why can the fluid-driven crack penetrate again? This is determined by the stress field around the fractures. In Section 3.2, the fracture propagation criterion based on elastic mechanics is discussed.
(a)
(b)
(c)
(d)
(e)
(f)
Another noteworthy feature different from the 2D model simulation is that fracture height is restrained by two cutters. In the 2D model, the fluid-driven crack must extend forward in their own form, penetrating or captured, when approaching two joints. However, in the 3D model, it is not obligatory for the fluid-driven crack extending forward as in the 2D. Figure 4(d)–4(f) is the aperture of vertical joints with different angles at 30°, 60°, and 90°, respectively. As shown in the figures, no matter what the angles, the vertical propagation is restricted by the two cutters. And it extends left and right. The fracture containment problem is another problem still to be solved.
3.2. Critical Pressure for the Failure of the Weakness Plane
There are many weak surface structures, such as bedding and fractures, in shale due to the action of geological tectonic stress. It is more likely to be damaged in fracturing for the reason that the strength of the weakness plane is low. The fracture will be preferentially expanded along the weakness plane. Therefore, the intersection of the hydraulic fracture and the weakness plane is the key problem of fracture propagation in the formation with the weakness plane.
The intersection model of the hydraulic fracture and the natural fracture is shown in Figure 5. Suppose there are two intersecting closed natural fractures, NF1 and NF2, at an approach angle of . The assumptions in this paper are as follows. The hydraulic fracture extends along the fracture plane NF2 parallel to the maximum principal stress. The hydraulic fracture tip reaches the intersection of the two fractures but has not yet passed the intersection point. NF1 and NF2 have certain tensile and shear strength before cracking, and the tensile and the shear strength disappear after cracking. To simplify the calculation, this paper selects the same method as the W&T criterion and assumes that the fracture is blunted at the interface, so that the singularity of the hydraulic fracture tip can be ignored. More detailed information can be found in [5].
The failure modes of NF1 and NF2 under fluid action are different with different angles of fracture and principal stress. In general, the failure form of NF1 may be a tensile failure or a shear failure. After tensile failure, cracks will open directly. However, cracks do not necessarily open after shear failure. It also requires that the fluid pressure can overcome the normal stress on the fracture surface and prop up the fracture. The failure mode of NF2 without shear stress on the fracture surface can only be tensile failure. The failure forms are summarized in Table 2. Different failure forms should meet the following conditions.
is the fluid pressure at the hydraulic fracture tip. is the critical pressure to shear failure of NF1. is the opening pressure of NF1 without bond strength. is the critical tensile failure pressure of NF1 considering the bond strength. is the critical tensile failure pressure of NF2 considering the bond strength.
Once the critical conditions are determined, the failure forms of fracture propagation can be determined after obtaining the failure conditions of different forms of fractures. Therefore, the critical pressure of different fracture failure forms will be calculated in the following part of this section.
According to Warpinski and Teufel, the critical fluid pressure in the crack during shear failure of NF1 can be calculated. where is the normal stress on the fracture surface, is the shear stress on the fracture surface, is the maximum principal stress, is the minimum principal stress, is the approach angle, is the inherent shear strength of the interface, and is the coefficient of friction.
The strength of the weakness plane disappears after cracking. When the fluid pressure inside the crack is greater than the normal stress on the fracture surface, the fracture surface opens. Therefore, the fluid pressure in the crack after shear failure is
Tensile failure occurs when the fluid pressure in NF1 is greater than the sum of normal stress and tensile strength of the fracture surface. At this point, the critical pressure of fluid in the joint is where is the tensile strength of the cemented fracture.
Tensile failure occurs when the fluid pressure of fracture NF2 is greater than the sum of normal stress on the surface and tensile strength of the fracture. At this time, the critical fluid pressure of tensile failure for fracture NF2 is
The pressure conditions for crack opening can be determined by the above formulas. When the natural fracture NF1 opens, the hydraulic fracture is captured by fracture NF1. When fracture NF2 opens, the hydraulic fracture passes through fracture NF1. The critical pressure of fracture failure at different angles of NF1 and NF2 is calculated to determine the fracture expansion form at the intersection point.
3.3. The Behavior of Hydraulic Fracture When Approaching Intersecting Cemented Joints
The critical pressure of different failure forms of cracks can be calculated by the formulas (11)–(14) and Table 2. The parameters selected in this article are shown in Table 3. And the results are shown in Figure 6.
Obviously, for NF1, the critical pressure of shear failure decreases first and then increases with the increase of the intersecting angle. The critical pressure of tensile failure and opening after shear failure increase with the increase of the intersecting angle. However, the critical pressure of the NF2 tensile failure is independent of the intersection angle.
At a small angle (0° and 10.67°), as the pressure in the crack rises, fluid pressure first reached . However, NF1 is not broken by the shear at this time, so it cannot be opened. The pressure then reached and . At this point, the two values are close; both NF1 and NF2 undergo tensile failure. Hydraulic fractures extend to both NF1 and NF2, forming cross fractures. At the angle range of 10.67° to 15°, the fracture fluid pressure in the first achieves Popen, and no fractures are failures now. Then, the pressure continued to rise to and shear failure occurred in NF1. At this point, NF1 shear failure and cracking occurred and hydraulic fracture only extended in NF1. When the fracture intersection angle is 15° to 37.7°, the fluid pressure first reaches , and then NF1 shear failure occurs. At this point, the pressure has not reached the opening pressure . As the pressure rises to , NF1, which has been shear damaged, cracks. When the intersection angle is 37.7° to 68°, the shear failure of NF1 occurs first, but it does not open, because the fluid pressure at this time is less than . With the increase of pressure, the fluid pressure reaches and the tensile failure of NF2 occurs. Since NF1 is not open, hydraulic cracks will only propagate along NF2. When the intersection angle is between 68 and 90 degrees, the fluid pressure first reaches and the hydraulic fracture will only extend in NF2. It is noteworthy that the values of and are close to each other in the range of 37.7 degrees, at which time NF1 and NF2 will open simultaneously.
To sum up, with the change of the intersecting angle, the expansion morphology of intersecting cracks can be divided into three categories. At low angles (0-10.67), NF1 and NF2 open at the same time, forming a cross-state. At the middle angle (10.67-37.7), shear cracking occurred in NF1 and no cracking occurred in NF2. The fracture morphology shows turning. At high angles (37.7-90), NF2 cracks and NF1 does not open. The fracture shape is a single fracture extending forward. In summary, the intersection angle can be divided into three parts according to the fracture propagation pattern, namely, low angle, medium angle, and high angle. The boundary between the low angle and the middle angle is defined as the L-M boundary (10.67 in Figure 6). And the boundary between the medium angle and the high angle is defined as the M-H boundary (37.7 in Figure 6). It should be noted that is close to at the M-H boundary, so the NF1 shear failure and the NF2 tensile failure may occur at the same time nearby. In addition, the failure mechanism of the crack is different even if the crack propagates in the same form. The concrete failure form should be judged according to its failure mechanism.
From Figure 6, it can be seen that in the low-angle region or near the M-H boundary, the crack number changes from 1 to 3. In the middle-angle area, the crack number is changed from the original 1 to 2. In both cases, the number of hydraulic cracks crossing the intersection point is higher than that before crossing. Whether the number of cracks increases or not is of great significance to the formation of a fracture network after fracturing. Therefore, it is more advantageous to form a complex fracture network at low and medium angles. When the angle is high, the cracks extend forward directly and the number of cracks does not increase, so it is not conducive to the formation of the fracture network. In summary, it is most advantageous to form a fracture network in the low intersection angle and angle boundary range, followed by a medium angle, and difficult to form an effective fracture network at a high angle. If the range of favorable angles can be widened, the effect of fracturing can be improved. Therefore, this paper carries out parameter analysis to clarify the influence of different parameters on the L-M boundary and M-H boundary.
4. Parameter Analysis
4.1. Effect of Minimum Principal Stress
From (11)–(14), we can see that has a separate term in each formula. When the difference of in situ stress is constant, if only the minimum principal stress is changed, the critical pressure will change the same value and their difference will not change. Therefore, the change of the minimum horizontal principal stress will not cause the change of the angle range. Therefore, in the case of constant in situ stress difference, the influence of the minimum principal stress on the angle range division of fracture propagation can be neglected. However, the increase of the minimum principal stress will lead to an increase of the critical pressure, so the increase of the minimum principal stress will increase the difficulty of crack opening.
4.2. The Effect of Inherent Shear Strength of the Interface
The inherent shear strength influences the shear behavior of cracks. The higher the inherent shear strength is, the greater the shear strength of cracks is. According to the propagation criterion of intersecting cemented cracks, the inherent shear strength only affects and the two are positively correlated.
The two angle boundaries vary with as shown in Figure 7. It can be seen that with the increase of , the L-M boundary value increases gradually, while the M-H boundary keeps constant. Figure 7(a) shows that when is small, the L-M boundary is close to 0, which is due to the shear failure of NF1 caused by the decrease of , leading to the left movement of the intersection of and . When is too large, the shear strength of NF1 is too high and the shear failure of NF1 will not occur. As can be seen from Figure 7(b), when , and do not intersect and the middle-angle range disappears. In short, as increases, the L-M boundary moves to the right and the low-angle range expands. The fracture network after fracturing is the most complex in the low-angle range. Therefore, the increase of can improve the complexity of the fracture network after fracturing.
(a)
(b)
4.3. The Effect of the Friction Coefficient
The influence of the friction coefficient on each boundary is shown in Figure 8(a). As can be seen from the figure, the friction coefficient only affects the L-M boundary. With the increase of , the L-M boundary moves to the left and the range of the middle angle becomes wider. The changes of critical pressures under different conditions are shown in Figure 8(b). As shown in the figure, as increased, the curve increased in the middle part while it decreased in the two ends. This causes the left intersection point of curve and to move to the left, resulting in the decrease of the L-M value.
(a)
(b)
4.4. The Effect of the Differential Stress
As mentioned in Section 4.1, the change of minimum principal stress has no effect on the boundary ranges. It is the differential stresses that affect the boundary range. The critical pressures (Figures 9(a)–9(c)) under different differential stress conditions (0 MPa, 2 MPa, and 5 MPa) were calculated, and the change curves of the L-M boundary and the M-H boundary with differential stress under 4-10 MPa were given (Figure 9(d)).
(a)
(b)
(c)
(d)
Figure 9(a) is the curve of each critical pressure when the differential stress is 0. It can be seen from the figure that under the condition of uniform stress, the intersection angle has no influence on the critical pressure. Moreover, and were equal and less than . Under this condition, both NF1 and NF2 will undergo tensile failure at the same time, forming cross-cracks. This indicates that under the condition of uniform stress, it will be advantageous to form complex fracture networks. Figure 9(b) shows the critical pressure curve when the stress difference is 2 MPa. It can be seen from the figure that under the condition of nonuniform stress field, the critical conditions of NF1 will change with the angle of intersection. The was lower than that under uniform stress. This indicates that the increase of differential stress reduces the difficulty of the shear failure of NF1. Under the condition of lower differential stress, as was greater than and , NF1 had not yet undergone shear failure. Therefore, in the case of low angle, NF1 and NF2 suffer from tensile failure at the same time, while in the case of high angle, only NF2 suffers from tensile failure. Figure 9(c) shows the critical pressure curves when the differential stress is 5 MPa. According to the figure, with the increase of the different stress, the critical financial pressure curve will become a classical three-interval model. Figure 9(d) shows the boundary change with the change of differential stress. As can be seen from the figure, both the L-M boundary and the M-H boundary decrease with the increase of differential stress. Since the angle range larger than the M-H boundary cannot form an effective fracture network, it indicates that it is more difficult to form a fracture network under the condition of high differential stress. In addition, when the boundary is less than the L-M boundary, the number of final cracks formed is the largest. The reduction of the L-M boundary will lead to the reduction of the complexity of the joint network. In conclusion, under the condition of uniform stress, a complex joint network will be formed. Low differential stress conditions are conducive to the formation of complex fracture networks, while high stress conditions are not conducive to the formation of fracture networks. The complexity of the fracture network formed under high differential stress is reduced.
4.5. The Effect of the Tensile Strength of Fractures
The tensile strength of cracks affects their tensile cracking. The greater the tensile strength, the higher the critical pressure for tensile cracking. The influence of the fracture tensile strength on each critical pressure is shown in Figure 10. Figure 10(a) shows the change of the boundary with the tensile strength. It can be seen from the figure that with the increase of the tensile strength, the L-M boundary decreases slightly, the M-H boundary increases, and the range of the middle angle expands. This indicates that the increase of the fracture tensile strength is conducive to the formation of the fracture network. According to Figure 10(b), as the fracture strength increased, the values of and increased, while and remained constant. As the difficulty of fracture tensile failure increases, shear failure is more likely to occur under the same pressure. Shear failure is the main cause of fracture steering. Therefore, the appropriate increase of fracture tensile strength is conducive to fracture diversion. But increased tensile strength also makes fracturing more difficult.
(a)
(b)
4.6. Multivariate Regression Analysis Based On the Response Surface Method
As mentioned in Section 3.3, the premise of forming a complex fracture network is the increase of the number of fractures. When the angle is less than the L-M boundary, the number of fractures changes from 1 to 3, and in the middle-angle range, the number of fractures changes from 1 to 2, both of which are conducive to the formation of the complex fracture network. So it makes sense to study the boundary. The above research can analyze the influence of a single factor on the two kinds of the boundary but cannot comprehensively evaluate the comprehensive influence of multiple factors. The response surface method can be used to study the dependence of the response variable on independent variables. If the number of independent variables is small (no more than 3), then the response surface method is one of the best methods. The response surface method can not only study the influence of individual variables on response variables but also study the influence of interaction between independent variables on response variables. The general form of the linear regression equation fitted by the response surface method is where is the number of variables, is the variable, and is the coefficient to be solved.
The selection and combination of independent variables are called the experimental design. The central composite design (CCD) is a common method for experimental design. The design tables of two experiments are shown in Tables 4 and 5. Specific variable selection will be based on different response variables.
4.6.1. L-M Boundary
According to the previous analysis, for the L-M boundary, the influence of the tensile strength of joints is negligible. Therefore, in the analysis of the L-M boundary, three independent variables, , , and , were selected—A: principal stress difference; B: inherent shear strength, ; and C: coefficient of friction, . The selection of variables and the experimental design are shown in Tables 4 and 5. Then, the RSM for the L-M boundary was conducted. The regression equation is shown in equation (16). And the analysis of variance is listed in Table 6.
The model -value of 10080.89 implies the model is significant. There is only a 0.01% chance that an -value this large could occur due to noise. Values of “” less than 0.0500 indicate that model terms are significant. And the significance order is .
The factor effect graph shows the linear effect of changing the level of a single factor. It is constructed by predicting the responses for the low (-1) and high (+1) levels of a factor. As shown in Figure 11, the L-M boundary is negatively correlated with principal stress difference and coefficient of friction and positively correlated with inherent shear strength. The result is consistent with the univariate analysis above.
(a)
(b)
(c)
4.6.2. M-H Boundary
From the univariate analysis, it is known that the M-H boundary is not affected by and . Therefore, in the analysis of the M-H boundary, only two independent variables and are used for analysis. They are : principal stress difference, and : tensile strength of cemented joint, . The selection of variables and the experimental design are shown in Tables 4 and 5. Then, the RSM for the L-M boundary was conducted. The regression equation is shown in equation (17). And the analysis of variance is listed in Table 7.
The model -value of 35479.45 implies the model is significant. There is only a 0.01% chance that an -value this large could occur due to noise. Values of “” less than 0.0500 indicate model terms are significant. In this case , , , , and are significant model terms. And the significance order is .
The factor effect graph of variances for the M-H boundary is shown in Figure 12. The M-H boundary is negatively correlated with principal stress difference and positively correlated with the tensile strength of the cemented joint. The result is consistent with the univariate analysis above.
5. Conclusion
In order to analyze the mechanism of hydraulic fracturing fracture propagation in formations with cemented joints, the critical pressure of three fracture propagation modes under the intersection of hydraulic fracturing fracture and closed natural fracture is derived, and the parameter analysis is carried out. The main conclusions are as follows: (1)The intersection angle can be divided into three parts according to the fracture propagation pattern, namely, low angle, medium angle, and high angle. Three types of behavior for a fracture are a cross-state at a low angle, turning at a middle angle, and extending forward at a high angle. The first two conditions are conducive to the formation of complex fracture networks(2)Univariate analysis was performed on each variable by the control variable method. The result shows that the increase of and can improve the complexity of the fracture network after fracturing. And low differential stress conditions are conducive to the formation of complex fracture networks. Besides, the appropriate increase of the fracture tensile strength is conducive to fracture diversion(3)The relations between the two boundaries and different independent variables are obtained, and their regression formulas are obtained by fitting by RSM
Data Availability
The data and code for numerical simulation and boundary calculation used to support the findings of this study are available from the corresponding author upon request.
Conflicts of Interest
The authors declared that they have no conflicts of interest to this work.
Acknowledgments
We would also like to thank Liuke Huang for his help in the calculation of the block discrete element method. The authors gratefully acknowledge the financial support given by the China National Science and Technology Major Project (Grant No. 2017ZX05037001).