Abstract
The perforations play an important role in determining the near-wellbore fracture geometry during the stimulation phase. To address the impact of perforations on fracture geometry in horizontal wells, a 3D coupled hydromechanical finite element model is developed and employed. Based on the theory of continuum damage mechanics, scalar damage variable governs the degradation of the stiffness of the solid. Damage affects the crack element modeling that is used to consider crack behavior and construct a crack-tracking algorithm to simulate propagation. The model was validated against the analytical solutions and perforation fracture experiments. The results indicate that perforation can be used to control the fracturing pressure and propagation behavior of the initial fracture, which has a further effect on the fracture geometry of near-wellbore region in horizontal wells. Optimizing perforation parameters can direct the propagation of the initial fracture toward the preferred fracture plane. The results demonstrate an improved capability to depict the 3D near-wellbore fracture geometry and fracture propagation with a continuum damage model. The model enables the optimization of orientations and perforation parameters, so that most efficient perforating completions can be designed for hydraulic fracture stimulation.
1. Introduction
Horizontal drilling and multiple-fracture completions of horizontal wells are effective techniques used today to extract hydrocarbons economically from unconventional reservoirs (shale oil, shale gas, tight gas) [1, 2]. Perforations provide the fluid conduit between the horizontal wellbore and the formation. The choice of perforation parameters, such as shot density, shot phasing, interval length, and gun orientation, has a significant effect on the geometry of the hydraulic fracture [3–6]. Optimized perforation numbers and location lead to significant improvements in well performance [7, 8].
To improve hydraulic fracturing efficiency, the perforation tunnel is suggested to align with the preferred fracture plane to reduce the tortuosity of near-wellbore zone [9]. The oriented perforation technique is used to minimize flow restrictions and friction pressures during fracturing for vertical wellbores. The fracture tortuosity in near wellbore of the horizontal well is undesirable, as they increase the average treating pressure and potentially impact the wellbore-to-fracture connectivity during production [10].
Gaining a better understanding of fracture initiation and propagation from perforations is essential for efficient hydraulic fracture treatment. Many researchers have analyzed the stress state near the well and perforation with different analytical and numerical modeling approaches. For the analytic approaches [11, 12], the wellbore and perforation cavities are treated separately, each being considered as an infinitely long cylinder. These cylinders are considered elastic and they are loaded with identical internal pressures. The stress state of a perforation is derived by summing the results of two subproblems of pressurized hollow cylinders. Numerical approaches have been applied to analyze the complicated stressed-state of the perforated wellbore and the fracture-initiation process in 3D. The boundary element method (BEM) was employed to investigate the problem of fracture initiation and perforation failure, and 3D models were built to evaluate the fracture-initiation pressure and predict the location and direction of the initial crack [13–15]. The dual boundary element method was employed to investigate the 3D nonplanar fractures evolution from a cavity [16].The finite element method (FEM) was used to obtain the 3D stress state around the horizontal wellbore by considering the completion process, and postfailure behavior model was adopted to simulate fracture propagation [17]. The 3D lattice based method (DEM) was employed to simulate complex hydraulic fracture initiation from perforation and the fracture development process [18].
The above numerical models have mainly focused on predicting the location and direction of initial crack. They rarely provide information on fracture propagation and coalescence among multiple perforations. In this study, a continuum damage mechanics (CDM) based constitutive model has been developed to accurately compute the coupled hydromechanical stress state near the perforation and wellbore and predict complex fracture geometry due to damage (initiation, propagation and coalescence of fractures) in the 3D problem statement. The FEM modeling approach developed in this paper was validated against the results of analytical solutions and perforation fracture experiments. The fracture-initiation pressure of this paper agrees well with the field-measured data. The visualization of complex fracture growth near a wellbore allows one to comprehend the physical process and optimize the perforation strategy.
2. The Governing Equations
To address the impact of perforations on hydraulic fracture initiation and complexity in horizontal wells, numerical simulation of the complicated near-well fracturing process is implemented in two parts. The first part mainly accounts for the fully coupled process of hydraulic fracturing in near-well region, i.e., the interaction between fluid flow and deformation (fracture), based on poroelasticity [19–22]. The second part accounts for postfailure behavior modeling, such as fracture propagation and coalescence, based on the continuum damage mechanics [23, 24].
2.1. Coupled Hydromechanical Models
Based on the theory of poroelasticity [25], the equilibrium equation of porous media with fluid saturated can be expressed as
Combined with the constitutive relationship of Hooke’s law, the elastic deformation equation of rock is obtained aswhere is the stress tensor, Pa; is the effective stress tensor, Pa; tension is taken as positive; is the body force in solid skeleton; is the Biot coefficient; is pore pressure, Pa; is displacement, m; is Young’s modulus, GPa; is Poisson’s ratio.
The mass balance equation of fluid flow in deformable porous media is written as
According to Darcy’s law for fluid flow, the flow rate is described aswhere is the storage coefficient, Pa−1; is the volumetric strain of rock mass; is the source-sink term of the seepage process, s−1; is rock permeability, μm2; is fluid viscosity, Pa∙s; is fluid density, kg/m3; is elevation, m.
Domain of the problem is illustrated as in Figure 1. Consider a porous continuum in the domain bounded by surface , the initial conditions and boundary conditions are stated as follows.
(i) Initial conditions
(ii) Boundary conditionswhere , are the displacement and pressure of initial domain, , are the displacement and pressure of boundary, is the boundary force, is the unit normal vector perpendicular to the boundary surface, is the fluid flux, , , , and are the bound surface of the displacement, stress, pore pressure, and flux, respectively.
2.2. Continuum Damage Evolution
During hydraulic fracturing process, the damage mechanisms in porous rocks are categorized based upon stress state, hydraulic pressure, and environmental conditions. The maximum tensile stress criterion is used for the determination of the location and direction of the initial crack [11, 13].where is the yield function of tensile failure, Pa; is the maximum principal stress, Pa; is the tensile strength of rock, Pa.
The stress distribution, hydraulic fracture initiation, and propagation of perforation tips and near wellbore are analyzed, based on the failure element concept. When the stress components are obtained, the fracture-initiation criterion (7) can recognize the failure point in the scale of element level. In continuum damage mechanics (CDM), damage can be considered as a macroscopic state variable that affects stiffness degradation of the material and can be used to modify the stiffness matrix for a failure element.
Information on failure element, time, and position is recorded in a dataset at each time step. Based on the damage mechanics theory, the Young’s modulus of the element degrades gradually during damage progress [23], and the damage evolution law of the damaged element is defined as follows:where and represent the Young’s modulus of the damaged elements and undamaged elements and is the damage variable.
Damage evolution of the element in tensional state can be expressed as [26] where is the equivalent principal strain for a tensile state, (=1,2,3) are the three principal strains, is the elastic limit of tensile strain, is the residual strength coefficient.
Considering the continuum damage evolution, the constitutive equation (2) can be further written as
Because of the damage, permeability is correlated to damage effect within the failure elements. The permeability can be described as [23]: where is the initial permeability, is a coupling coefficient, is the effective average stress, is the mutation coefficient, which reflects the permeability mutation in the damage element.
3. Numerical Solution
3.1. Finite Element Formulations
Using the Galerkin weighted residual method and the Green-Gauss theorem [27], the force equilibrium equation (1) and the continuity equation for fluid flow (3) can be converted into the matrix form after finite element discretization:where , are the vectors of displacement and pressure and , are corresponding time derivatives. The matrix expressions are listed as follows.
In order to integrate equation (13) with respect to time, the linear interpolation in time can be discretized with standard finite difference methods, and the incremental forms of the above FE equations can be written as where , represent the solution at the current time step, , is the solution at last time step, is an integration parameter in the interval . In the following numerical example, we choose .
The displacement convergence criterion is used to judge the convergence of incremental solutionwhere is the iterative tolerance of the displacement increment, is the displacement increment at current time step, and is the displacement at last time step.
3.2. Computational Procedures
Using FEPG software as a platform, an incremental finite element procedure is developed using the continuum damage evolution criteria, and the iterative decoupling method is used to solve the above equations of geomechanics and reservoir flow variables.
Given a set of initial values, the fluid pressure is computed first from (16). Then the displacement, strain, and stress are computed from (15) from updated and boundary conditions. The coupling iteration is controlled by a convergence criterion, which is based on (17).
With loading history changing, the stress concentration region gradually reaches the fracture-initiation criterion (7). Then, this region’s elements are defined as damage elements. The damage elements’ matrices and are determined using (9) and (11). In this way, the initiation and growth of fracture can be described using continuum damage mechanics. The flow chart of the iteration process is shown in Figure 2.
3.3. Verification of the Solution Algorithm
The coupled hydromechanical solution algorithm is verified against the Terzaghi’s consolidation problem [28, 29]. This problem consists of a fluid-saturated column of soil layer with loading on top. Drainage is only allowed through the top surface of the soil column. Figure 3 shows the geometry and boundary conditions of this problem, and the load is applied instantaneously at time =0. For solution algorithm verification, we choose Youth’s modulus =1.0 GPa, Possion’s ration =0.25, permeability =1.0 um2, fluid viscosity =1.0 × 10-3 Pa∙s, Biot coefficient =1.0, and initial pressure =1.0 MPa.
As is shown in Figure 4, a good matching between the analytical solution [28, 29] and the numerical result is obtained for both pore pressure along the vertical direction at different time and vertical displacement at the top versus time.
(a)
(b)
4. Computational Model and Validations
4.1. Model Configuration
Field examples of horizontal wells in Changqing oilfield (Xi’an, China) are introduced. Numerical simulations have been performed to analyze the developed fracture geometry for traditional spiral perforation and in-plane perforation in a horizontal well. The numerical results were compared with analytical solutions and experimental geometry.
The wellbore is assumed to be located at the center of a block of 3.0 m × 5.0 m × 5.0 m dimensions, and the near-wellbore block of 1.0 m × 1.5 m × 1.5 m is selected to local encryption of element mesh. The horizontal borehole axis aligned with the minimum horizontal stress direction. The model has one perforation cluster, which includes six perforation tunnels, and the phase angle of perforation tunnels is 60°. We assume the direction of two middle perforations is consistent with the maximum horizontal stress . The perforation tunnel length is 400 mm, and the diameter of perforation tunnel is 20 mm.
The applied far-field in-situ stresses: vertical stress =43.0 MPa, maximum horizontal stress =38.0 MPa, and minimum horizontal stress =34.0 MPa are applied on three sides, as shown in Figure 5(a). The other three sides of the model are applied as roller boundary conditions. For fluid phase, the pore pressure of the reservoir =18.0MPa is applied as the initial condition. No-flow boundary condition is applied on all sides.
For an open-hole horizontal well, it was assumed that the same bottom-hole pressure (BHP) loads were applied on the wellbore and perforation surfaces. For a cased horizontal well, the BHP acting on the inside of the casing wall is partially transmitted to the rock face, since the stiffness of casing is much higher than the rock. An approximate approach to the evaluation of pressure applied to the borehole boundary in a cased-and-cemented wellbore has been proposed [8, 13]:where , are Young’s modulus and Poisson’s ratio of the steel, is the transmission factor, is the in-situ stress factor, , are the inner and outer radii of the casing, and is the inner radius of wellbore.
The reservoir properties and well parameters are reported in Table 1.
Figure 6 presents a typical field fracturing records of surface injection pressure and injection rate. The increasing injection rate is up to 8.0-11.0 m3/min. From the pump pressure curve, the average breakdown pressure is 26.9MPa. This data can be converted to the fracture initial pressure () of 43.1MPa, with the addition of a hydraulic head of 18.0MPa and the friction loss of almost 1.8 MPa. In this case, a varying pressure boundary condition is applied on the inner face of perforation tunnels and wellbore. According to the field data of Figure 6, the pressurization rate is settled as 0.05 MPa/s in this case.
4.2. Model Validations
The FEM model was validated against the results of analytical model and perforation fracture experiments. The tangential stress distribution in the open-hole horizontal wellbore is shown as Figure 7. We compare the FEM results with the analytical solution (AS) of tangential stress along the wellbore and perforation tunnel [11, 12]. Figure 7(a) presents that the FEM results are close to the AS results with a difference of less than 3.5%.
(a)
(b)
Figure 7(b) compares the tangential stress along the perforation tunnel with change of perforation deflection angle β. For β=0°, the maximum tangential stress is at tangential angle of perforation tunnel =0°. While for β=90°, the maximum tangential stress is at =90°. The difference between FEM results and AS solutions of the tangential stress along the perforation tunnel is less than 6.5%.
Figure 8 compares the damage fracture geometry in near wellbore of horizontal well between numerical simulation results and that of the perforation fracture experiments [17]. One can see from the numerical results that fracture extends through each perforation as they spiral around the casing and form the fracture geometry with spiral distortion. The laboratory experiment simulated a horizontal well with the wellbore axis along the minimum stress direction, and the fracture geometry is also spiral around the casing. The above results show that the developed FEM model can be used to analyze the stress and fracture geometry in the near-well region.
(a)
(b)
5. Results and Discussions
5.1. Sensitivity of Fracture Initiation to Perforation Orientation
The effect of perforation orientation on initial fracture geometry is analyzed through various angle β relative to the vertical direction.
With the increase of perforation deflection angle, as shown in Figure 9, two possible scenarios of fracture initiation are observed for open-hole perforation horizontal well. With the comparatively small angle (β=0°, β=30°), the zones for rock failure concentrate on the junction site of perforation and the borehole wall. The fracture is along the longitudinal direction of the wellbore (Figures 9(a) and 9(b)). With the increase of the perforation angle (β=60°, β=90°), the initial rupture position is located at the middle and toe end of the perforation tunnel, forming a transverse rupture perpendicular to the horizontal wellbore axis (Figures 9(c) and 9(d)).
(a)
(b)
(c)
(d)
For cased-cement horizontal well, the initiated fractures are all transverse to the wellbore (Figure 10). For β≤60°, fractures initiate from the perforation base (Figures 10(a), 10(b), and 10(c)). While β continues increasing, the initial fracture position is extended from base to the middle and tip of the tunnel (Figure 10(d)).
(a)
(b)
(c)
(d)
According to the numerical results, the initiation of the fractures is sensitive to perforation orientation. The change of stress conditions controlling the fractures initiation can explain this phenomenon.
When the perforation angle β=0°, as shown in the computational model (Figure 5(d)), the axis of the tunnel is along the direction of vertical stress . The fracture initial pressure of the tunnel is mainly affected by the maximum horizontal stress and minimum horizontal stress . As the perforation angle β=90°, the axis of the tunnel is along the direction of the maximum horizontal stress. The fracture initial pressure is dominated by the vertical stress and minimum horizontal stress .
As shown in Figure 11, the difference of fracture pressure between open-hole horizontal wells reaches 6.5 MPa. For the case-cement horizontal well, the cementation effect of casing and cement ring reduced the pressure load applied on the near borehole rock [13]. Figure 11 shows that the casing and cement layer result in decrease of the fracture initial pressure 4.0~7.5 MPa than open-hole horizontal well. With the increase of perforation angle β, the fracture initial pressure grows up to its maximum at β=90°, the difference between various perforation angles for case-cement horizontal well is 3.5MPa.
5.2. The Initial Fracture Characters of Horizontal Wells
Figure 12 presents the maximum stress principal change of spiral perforations during the pressurization period. With the BHP applied on the perforation tunnel reaching 41.5 MPa, the maximum principal stress of element on the middle tunnel reaches 3.02 MPa and surpasses the tensile strength of the rock. The initial fracture appears at the interface of wellbore and middle perforation tunnel, leading to fracture perpendicular to the wellbore (Figure 13(a)).
(a)
(b)
(c)
Based on continuum damage mechanics, the failure elements’ matrices of stiffness and seepage are updated. The postfailure and growth of fracture can be described by the damage-driven crack element modeling. As increases to 43.5 MPa, the maximum principal stress on the side perforation tunnels is also beyond the tensile strength and gives initiation of cracks extending along the perforation axis (Figure 13(b)). With the pressure increase, the damaged zone extends and links up among multiperforations (Figure 13(c)).
The stress interference between the spiral perforation channels controls the direction of fracture propagation in local region of perforation, inducing the fracture tortuosity in near wellbore. According to the FEM results, the fracturing pressure of different perforations ranges between 41.5 MPa and 43.5 MPa, which is in a good agreement with the field fracturing recording FIP 43.1MPa.
The fracture tortuosity (Figures 8(a) and 13(c)) in near-wellbore region of the horizontal well is undesirable, as they increase the average treating pressure and the flow restriction.
The in-plane perforation technique was proposed to change conventional spiral distribution of firing holes and adopt special charge distribution mode for horizontal wells [9]. The in-plane perforation technique designs perforation axis aiming at the preferred fracture plane, and the stress interference of multiperforations increases.
The maximum principal stress change of in-plane perforation is calculated and shown in Figure 14. In this case, the side perforations incline angle γ set as 25°. The maximum principal stress of element on the middle tunnel surpasses the tensile strength of rock at BHP =39.5 MPa. This value is 2.0 MPa lower than the conventional spiral perforation case, which helps reduce the FIP of hydraulic fracture. According to the field fracturing recording of Changqing oilfields, the average fracturing pressure for in-plane perforations horizontal wells is 39.2MPa. This data is consistent well with the calculation results of this paper.
By adjusting the tunnel inclined angle γ, in-plane perforation can also control the initial fracture directions of the perforation. As shown in Figure 15, the fracture of side perforations initiates and extends along the perforation axis, leading to the coalescence of fractures.
(a)
(b)
(c)
5.3. The Near-Wellbore Fracture Geometry of Horizontal Wells
Figure 16 depicts the damage evolution process around the horizontal wellbore for different perforation strategies. For open-hole horizontal well, the crack of multiperforations extends along the axial and vertical directions of the wellbore and form a complex near-wellbore rupture pattern, as illustrated in Figure 16(a). The fracture geometry of case-cement horizontal well is spiral around the casing with spiral distortion (see Figure 16(b)), which increases the average treating pressure and affects the wellbore-to-fracture connectivity.
(a)
(b)
(c)
Figure 16(c) presents the near-wellbore fracture geometry of in-plane perforation horizontal well. The crack extends along the inclined perforation and deflects to the damage zones of the middle tunnels. Two preferred fan shaped fracture plane perpendicular to the axial of the casing. The in-plane perforations guide the initial fracture trajectory toward the preferred transverse fracture plane. This kind of fracturing morphology would avoid the tortuous near-wellbore fracture geometry beneficial to wellbore-to-fracture connectivity.
The above sensitivity study performed revealed that perforation is an effective strategy to control the fracturing pressure and the initial fracture position with the variation of angle between the direction of the perforation channel and the preferred fracture plane. A complex and tortuous near-well fracture geometry tends to be induced in conventional spiral perforations. With optimized direction of jet flow of the in-plane perforation, the initial fracture propagation trajectory directed toward the preferred fracture plane and reduced the tortuosity of near-wellbore. This will improve the efficiency of hydraulic fracturing stimulation.
6. Conclusions
(1)A three-dimensional coupled hydromechanical fracturing model is established to investigate the propagation of near-wellbore fracture geometry, together with a continuum-based damage mechanics constitutive relationship.(2)Finite element program is developed and employed to analyze the impact of perforations on fracture geometry in horizontal wells. The model was validated against the analytical solutions and experimental results, the capability of modeling 3D near-wellbore fracture geometry and fracture propagations demonstrated.(3)Perforation angle variations result in variable initiation pressure and complex fracturing pattern for horizontal wells. By increasing the stress interference effect among perforations and changing the direction, it benefits to lower the fracture initial pressure.(4)Perforation has a further effect on the fracture geometry of near-wellbore region in horizontal wells. Optimizing perforation parameters can direct the propagation of the initial fracture toward the preferred transverse fracture plane benefiting wellbore-to-fracture connectivity. Consequently, it improves the efficiency of hydraulic fracturing stimulation.
Data Availability
The data used to support the findings of this study are available from the corresponding author upon request.
Conflicts of Interest
The authors declare that there are no potential conflicts of interest regarding the publication of this article.
Acknowledgments
This work was supported by the National Natural Science Foundation of China (NSFC) (No. 51304230), the Fundamental Research Funds for the Central Universities (16CX05001A), and the China Scholarship Council Fund (201706455025).