Abstract
There has been a growing consensus that preexisting natural fractures play an important role during stimulation. A novel fully coupled hydromechanical model using extended finite element method is proposed. This directly coupled scheme avoids the cumbersome process during calculating the fluid pressure in complicated fracture networks and translating into an equivalent nodal force. Numerical examples are presented to simulate the hydraulic fracture propagation paths for simultaneous multifracture treatments with properly using the stress shadow effects for horizontal wells and to reveal the deformation response and interaction mechanism between hydraulic induced fracture and nonintersected natural fractures at orthotropic and nonorthotropic angles. With the stress shadow effects, the induced hydraulic flexural fracture deflecting to wellbore rather than transverse fracture would be formed during the progress of simultaneous fracturing for a horizontal well. The coupled hydromechanical simulation reveals that the adjacent section to the intersection is opened and the others are closed for orthogonal natural fracture, while the nonorthogonal natural fracture is activated near the intersection firstly and along the whole section with increasing perturbed stresses. The results imply that the induced hydraulic fracture tends to cross orthotropic natural fracture, while it is prior to being arrested by the nonorthotropic natural fracture.
1. Introduction
Unconventional reservoirs, such as shale gas, have an extremely tight matrix, but well-developed natural fractures [1]. Volume stimulation drives these natural fractures to expand continuously and introduce shear slip in brittle rocks, thus forming complicated networks of intersecting natural fractures and artificial fractures that further increase the stimulated volume [2]. Therefore, the volume stimulation method has become one of the key technologies in unconventional oil and gas development [3]. However, simulation hydraulic fracturing in the presence of natural fracture network is a challenging task, owing to the complex interactions between fracturing fluid, rock matrix, and rock interfaces, as well as the mechanical behaviors before, during, and after the intersections between hydraulic fractures and preexisting natural fractures [4–6].
The morphology of complicated fracture networks created by volume fracturing can be estimated by determining the mechanical mechanism and propagation behavior of hydraulically induced fractures before, during, and after intersecting with natural fractures. In laboratory experiments, a true triaxial fracturing simulation system can be used in combination with high-energy CT scanning, acoustic emission, and microscopic X-ray chromatography to characterize the intersection behavior between hydraulically induced fractures and natural fractures during the fracturing of artificial rock or reservoir outcrop rock samples. These methods also reveal the influences of in situ stress state, horizontal principle stress difference, approach angle, natural fracture interface properties, injection fluid amount, and fracturing fluid viscosity [7–11]. Scholars have developed a series of criteria to describe the intersection behavior between hydraulic fractures and natural fractures. For example, Renshaw and Pollard proposed a criterion for the intersection between hydraulic fractures and orthogonal fractures (R-P criterion) [12]. Gu and Weng expanded the R-P criterion to encompass nonorthogonal conditions (G-W criterion) by considering that the approach angles of hydraulic fractures and natural fractures are between 0° and 90° [13]. Meanwhile, Chen et al. considered the dip angle of natural fractures and expanded the R-P criterion to three-dimensional conditions [14]. Dahi-Taleghani and Olson investigated the tensile property and shearing behavior of natural fractures while hydraulic fractures are approaching but before they come into contact [15]. All of these criteria only consider the mechanical behavior caused by the interaction between hydraulic fractures and natural fractures, without accounting for the influence of the hydraulic effect induced by the injection of fracturing fluid.
Volume fracturing presents a complicated hydromechanical coupling problem, with two special conditions. Firstly, there would be fluid flow within the fracture network and matrix rock and fluid exchange or fluid leak between the fracture and the surrounding porous media, because of the pressure difference between the fracture network and the matrix rock during fracturing fluid pumping. Secondly, although the net pressure drives the initiation and propagation of fractures, the increased pore fluid pressure in the matrix surrounding the fractures would impact fracture propagation. There were considerable differences between piston-like displacement, single phase and multiphase leakage behaviors for nonintersected fractures, and matrix rock using extended finite element method (XFEM) [16]. Based on the vital important role of leakage, this paper promoted the multiphase fluid flow model to fully coupled hydromechanical model for fractured porous media to provide a numerical tool for simulating hydraulic fracturing of fractured porous media. This decision is enlightened on the strong discontinuity of displacement at either side of fractures and the weak discontinuity of pore fluid pressure, as XFEM represents the most commonly used method in simulating discontinuous media. This paper also discusses the deformation features and mechanical mechanisms before the intersection between hydraulically induced fractures and orthogonal and nonorthogonal natural fractures.
2. Fully Coupled Hydromechanical Model
The fractured porous media were composed of matrix rock and macro fractures as a whole. Although there were plenty of micro fractures, which show performances similar to those of macro fractures, developed in matrix rock, the fractured porous media were considered to be constructed by an intact matrix and macro fracture networks for simplicity, as displayed in Figure 1.
The strong form, weak form, and the associated discretization form of governing equations for hydromechanical coupling inside both matrix domain and fracture domain are demonstrated in this section.
2.1. Strong Form
Figure 2 shows a bounded control body () with an outer boundary of .
The effective stress equilibrium equation for porous media is written as
According to the law of mass conservation, the continuity equations for the wetting and nonwetting phase fluid flow in the matrix domain are expressed as
Similarly, the continuity equations for the wetting and nonwetting phase fluid flow inside the fracture domain are given by
The displacement and pressure of the wetting and nonwetting phase initially are provided according to the initial conditions
The displacement and pressure of the wetting and nonwetting phase at the boundary are provided by the forced boundary conditions
The stress and mass flow rate of the wetting and nonwetting phase at the boundary are provided by the natural boundary conditions
The fluid pressure is applied on the fracture surface
The pore fluid pressure at the fracture boundary is continuous, with fluid exchange occurring between fractures and the surrounding matrix:
2.2. Weak Form
The weak form of the governing equations can be obtained by combining the weighted residual method with the associated boundary conditions. With regard to fluid flow in fractures, since the fracture width is much smaller than its length, the fluid pressure variation along its width can be neglected (see Figure 3). The areal integral inside the fracture domain can be transformed into the linear integral along the fracture surface . The fluid exchange items or leak between the fractures and the surrounding matrix rock can then be vanished by adding the equivalent integrals of the matrix rock and fractures, which have the same absolute value but opposite signs. Thus, it becomes unnecessary to explicitly express the fluid exchange items or leak [16–18]:
2.3. Discretization Form
XFEM introduces Heaviside step functions and branch functions to enrich the shape function, which make it easy to deal with typical strong discontinuous problems like fractures without grid matching and remeshing. The XFEM displacement field for multiple fractures is [19]
In fractured porous media, the natural fractures are the main channels for seepage. The pore fluid pressure field is continuous at the fracture surface, while the derivative (i.e., flow rate) is discontinuous. This presents a typical weak discontinuous problem. Similar to the interpolation of displacement field, the pressure field can be interpolated as [16, 18]
The discretization form of the fully coupled hydromechanical model for fractured porous media can be obtained by regarding the trial functions and in (9)–(11) as the variation of displacement and pressure interpolation functions, as shown in
The stiffness matrix , coupling matrix , flow matrix , compressibility matrix , saturation matrix , equivalent nodal force vector , and equivalent node flow rate vector are assembled matrixes for each element matrix labeled by the degrees of freedom, where each element matrix is expressed as
Equation (14) is implicitly solved in the time domain using the generalized trapezoidal method. The implicit pressure explicit saturation method (IMPES) is used to treat the multiphase fluid flow. This strong nonlinear multiphase fluid hydromechanical coupling problem can be solved by incorporating the Newmark scheme [20].
3. Numerical Simulation
3.1. Simultaneous Fracturing for Horizontal Wells
The cluster perforation and simultaneous fracturing are key technologies for volume stimulation of horizontal wells. The propped hydraulic fracture induces additional stresses in the surrounding rock. This stress interference between fractures can complicate the propagation of hydraulic fractures to form complex fracture networks. Crosby et al. [21] and Rahman et al. [22] investigated the propagation behavior of two isolated hydraulic fractures induced by the simultaneous fracturing of a horizontal well.
This study simulated Crosby’s DC-8 specimen, as Figure 4 shows. Modeling was performed using material properties and injection parameters similar to those employed during the laboratory test DC-8. The specimen was meshed into uniform 80 × 80 structured bilinear quadrilateral elements, with the propagation step set at a constant value of 0.01 m. The coupling process was neglected during this simulation because the synthetic fracture test material was cast by silica flour and acrylic emulsion to provide the specimen with a negligible leak coefficient, 3.8 × 10−8 , and an ultralow permeability, 8 × 10−5 mD. Therefore, only the first item in (14) was used. Equation (9) states that the fluid pressure inside the fracture, , should be calculated, which presents a very complicated problem, especially for fractures that are not straight. Moreover, the fluid pressure inside the fracture must be transferred into an equivalent nodal force using the linear integral form, which requires integral points on the fracture surface [23]. We set the net pressure inside the fractures at 3 MPa.
The maximum circumferential stress criterion is selected to determine the propagation condition and crack direction. The stress intensity factors of hydraulic fracture were calculated with the interaction integral
The first item of the right side can be calculated easily. The second item of the right side was derived explicitly:
Figure 5 shows comparisons between the XFEM simulation results, experimental results, and HYFRANC finite element simulation results.
The results of XFEM simulation are more consistent with the experimental result in the upper fracture than in the lower fracture. This can be largely attributed to the fact that the lower fracture was initiated preferentially during the experiment. Therefore, the propagation path was relatively straight. However, the two fractures were considered to initiate at the same time in the simulation. Then, the lower fracture was forced to deflect towards the wellbore for stress shadow effect. In comparison, the HYFRANC simulation results [21] and the XFEM simulation results show high consistency with the experimental results. Moreover, the simulation efficiency improved significantly.
The stress distribution at the sixth propagation step (Figure 6) shows that the -axis stress, , between the two fractures increased from 5 MPa to about 7.62 MPa (compressive), the -axis stress, , increased from 7 MPa to about 8.06 MPa (compressive), and the shear stress, , was lower than 0.0002 MPa. In other words, the stress state between the two fractures changed significantly. The maximum horizontal nominal stress direction changed from the -axis to the -axis. Perforation fracturing between two fractures tends to cause the hydraulic fracture to propagate along the -axis, intersect with other fractures, and form a fracture network.
(a) -axis stress (MPa)
(b) -axis stress (MPa)
(c) planar shearing stress (0.1 kPa)
3.2. Coupled Hydromechanical Behavior for Hydraulic Fracturing
Figure 7 shows a 2D fractured porous media control domain. The physical model has a length and width of 150 m, which contains three fractures. Table 1 lists detailed parameters for the model. The hydraulic fracture, which is parallel to the global -axis, has a length of 50 m and an initial width of 1 mm. The orthogonal natural fracture, which is parallel to the global -axis, has a length of 40 m and an initial width of 0.1 mm. The nonorthogonal natural fracture, which lies at an included angle of 45° to the global -axis, has a length of 30 m and an initial width of 0.1 mm. The initial stress is zero, while the pore fluid pressure inside both the fractures and the matrix rock is 50 MPa. The displacements of the upper, lower, right, and left sides are all restricted (), and they are therefore impermeable boundaries. In reality, the length, width, height, and net pressure would be increased with increasing injection time, and the stress intensity factors and induced stresses would be strengthened to force the deformation of natural fracture, subsequently. For simplicity, the fractures were presumed as nonpropagating cracks to investigate the deformation behavior between hydraulic and natural fractures before intersection. However, the net pressure and width of fractures were solved with the fully coupled hydromechanical model. Although this assumption was not coinciding with the reality, the numerical simulation can represent the deformation behavior between hydraulic and natural fractures to some extent.
A uniform 80 × 79 bilinear quadrilateral mesh was used for calculation, with the intention of locating the injection well point at the element boundary, while the fracture location was independent of the generated mesh.
The -axis displacement, , distribution is shown in Figure 8. The -axis displacement, , distribution is shown in Figure 9. The displacement fields demonstrate that the injected fluid initially flows along the first fracture, forcing the fluid pressure inside the fracture to increase. Meanwhile, the fluid filtrated into the surrounding matrix rock, leading to an increase of pore fluid pressure in the matrix rock. The increase in fluid pressure drives the deformation of fractured porous media, generating displacement along the -axis and the -axis. Since the contact condition is not considered here, parts of the other fracture surfaces may come into contact and intersect with each other, while other parts become staggered.
(a) = 60 s
(b) = 600 s
(a) = 60 s
(b) = 600 s
Fracture width represents the key parameter in hydraulic fractures. Equation (12) allows the displacement of the two surfaces of a fracture at any position to be obtained. The fracture width can be calculated by subtracting the positive-side displacement from the negative-side displacement. The width decreases with increasing distance from the injecting well point, and it increases with increasing injection duration, as Figure 10 shows. The fracture width is distributed symmetrically. The orthogonal natural fracture remains closed at first under the effect of stress induced by the left crack tip of a hydraulic fracture. As the fluid pressure inside the hydraulic fracture increases, the crack-tip-induced stress increases and forces the center section to open, while the faraway section remains closed, as shown in Figure 11. This is one of the reasons why a hydraulic fracture intersects a natural fracture at a large intersection angle. In contrast, the stress induced by a hydraulic fracture drives the nonorthogonal natural fracture to open. As Figure 12 shows, only a small part of the initial fracture remains closed, while the whole fracture opens finally.
4. Conclusion
Hydraulic fracturing is an interactive process incorporating both mechanical and hydraulic effects. The hydromechanical coupling model and associated solution scheme proposed under the framework of XFEM can effectively avoid tedious calculations of both the fluid pressure distribution in fractures and the transformation of equivalent nodal forces. This model also considers the influence of pore pressure variation in the rock matrix on fracture propagation.
The mechanical module can simulate the fracture propagation pathway of the simultaneous fracturing for a horizontal well. Under the influence of stress perturbation between hydraulic fractures, not the straight transverse fractures but the flexural fractures deflecting toward the wellbore were created for the staged fracturing of a horizontal well. The stress shadow effect introduced by the propped hydraulic fracture reduces the horizontal stress differential between fractures or even leads to stress inversion. Therefore, it is beneficial to enhance the complexity of fracture networks.
The induced stresses added to the crack tip of the hydraulic fracture keep most of the orthogonally intersected natural fractures closed, with only a small part opening due to the induced stress and fracturing fluid leakage. The stress concentration at the crack tip of the hydraulic fracture fully opens most of the nonorthogonal intersected natural fractures as time passes. The hydraulic fracture tends to traverse the orthogonal natural fracture and continue to propagate, while it is more easily arrested by the nonorthogonal natural fracture.
Nomenclature
| : | Overall stress |
| : | Effective stress |
| : | Fluid density |
| : | Fluid density of wetting phase |
| : | Fluid density of nonwetting phase |
| : | Gravity acceleration vector |
| : | Unit vector |
| : | Biot coefficient |
| ϕ: | Porosity |
| : | Bulk modulus of rock particle |
| : | Bulk modulus of wetting phase |
| : | Bulk modulus of nonwetting phase |
| : | Saturation of wetting phase in matrix rock |
| : | Saturation of nonwetting phase in matrix rock |
| : | Saturation of wetting phase in fracture |
| : | Saturation of nonwetting phase in fracture |
| : | Pore fluid pressure |
| : | Capillary pressure of matrix rock |
| : | Capillary pressure of fracture |
| : | Wetting phase fluid pressure in matrix rock |
| : | Nonwetting phase fluid pressure in matrix rock |
| : | Initial wetting phase fluid pressure in matrix rock |
| : | Initial nonwetting phase fluid pressure in matrix rock |
| : | Wetting phase fluid pressure at the boundary |
| : | Nonwetting phase fluid pressure at the boundary |
| : | Wetting phase fluid pressure in fracture |
| : | Nonwetting phase fluid pressure in fracture |
| : | Permeability tensor of matrix rock |
| : | Fracture permeability |
| : | Relative permeability of wetting phase in matrix rock |
| : | Relative permeability of nonwetting phase in matrix rock |
| : | Relative permeability of wetting phase in fracture |
| : | Relative permeability of nonwetting phase in fracture |
| : | Viscosity of wetting phase |
| : | Viscosity of nonwetting phase |
| : | Ground source sink term of wetting phase |
| : | Bulk coefficient of wetting phase |
| : | Displacement vector |
| : | Boundary displacement vector |
| : | Initial displacement |
| : | Strain tensor |
| : | Boundary force tensor |
| : | Flow rate of wetting phase at the boundary |
| : | Flow rate of nonwetting phase at the boundary |
| : | Unit normal vector at the boundary |
| : | Value difference between two sides of fracture |
| : | Transformation matrix between local coordinate and global coordinate |
| : | Unit normal vector in local coordinate system of fracture surface |
| : | Unit tangential vector in local coordinate system of fracture surface |
| : | Interpolating shape function |
| : | Heaviside function |
| : | Ramp function |
| : | Branch function for displacement |
| : | Junction function |
| : | Signed distance function |
| : | Branch function of pore pressure |
| : | Interpolating function matrix of displacement |
| : | Interpolating function matrix of pore pressure |
| K: | Stiffness matrix |
| L: | Coupling matrix |
| C: | Compressibility matrix |
| S: | Saturation matrix |
| B: | Strain matrix |
| : | Equivalent nodal force vector |
| : | Elasticity matrix of rock |
| : | Equivalent nodal flow rate vector |
| : | Fracture width. |
Conflicts of Interest
The authors declare that there are no conflicts of interest regarding the publication of this paper.
Acknowledgments
The authors are grateful for the research support of the National Natural Science Foundation of China (NSFC, no. 51404207) and the Postdoctoral Program of Postdoctoral Work Centre, Southwest Oil & Gas Field Company (no. 20150304-08).