# Continuum Damage Modeling of Hydraulic Fracture from Perforations in Horizontal Wells

Authors: Feng Sun; Peng Jia; Shifeng Xue
Journal: Mathematical Problems in Engineering
 (2019)

---

## 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.


---

## Body

## 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(1)σij,j+fi=σij,j′-αp,i+fi=0Combined with the constitutive relationship of Hooke’s law, the elastic deformation equation of rock is obtained as(2)E21+νui,jj+E21+ν1-2νuj,ji-αp,i+fi=0where σij is the stress tensor, Pa; σij′ is the effective stress tensor, Pa; tension is taken as positive; fi is the body force in solid skeleton; α is the Biot coefficient; p is pore pressure, Pa; u is displacement, m; E is Young’s modulus, GPa; ν is Poisson’s ratio.The mass balance equation of fluid flow in deformable porous media is written as(3)S∂p∂t-∇·q+α∂εv∂t=QfAccording to Darcy’s law for fluid flow, the flow rateq is described as(4)q=-kμ∇p+ρfg∇Zwhere S is the storage coefficient, Pa−1; εv is the volumetric strain of rock mass; Qf is the source-sink term of the seepage process, s−1; k is rock permeability, μm2; μ is fluid viscosity, Pa∙s; ρf is fluid density, kg/m3; Z is elevation, m.Domain of the problem is illustrated as in Figure1. Consider a porous continuum in the domain Ω bounded by surface Γ, the initial conditions and boundary conditions are stated as follows.Figure 1
Schematic of the problem domain.(i) Initial conditions(5)ui0=u0,p0=p0inΩ(ii) Boundary conditions(6)ui=u-tonΓuσij·n=TitonΓσp=p-tonΓp-kμ∇p=q-tonΓqwhere u0, p0 are the displacement and pressure of initial domain, u-(t), p-(t) are the displacement and pressure of boundary, Ti(t) is the boundary force, n is the unit normal vector perpendicular to the boundary surface, q-(t) is the fluid flux, Γu, Γσ, Γp, and Γq 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].(7)Fσ1=σ1-ft0where F(σ1) is the yield function of tensile failure, Pa; σ1 is the maximum principal stress, Pa; ft0 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:(8)E=1-DE0where E and E0 represent the Young’s modulus of the damaged elements and undamaged elements and D is the damage variable.Damage evolution of the element in tensional state can be expressed as [26] (9)D=0εt<εt01-λεt0εtεt≥εt0where εt=∑i=13εi+εi/22 is the equivalent principal strain for a tensile state, εi (i=1,2,3) are the three principal strains, εt0 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(10)E1-D21+νui,jj+E1-D21+ν1-2νuj,ji-αp,i+fi=0Because of the damage, permeability is correlated to damage effect within the failure elements. The permeability can be described as [23]: (11)k=ξk0exp⁡-ησ-′0<D≤1where k0 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.


## 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(1)σij,j+fi=σij,j′-αp,i+fi=0Combined with the constitutive relationship of Hooke’s law, the elastic deformation equation of rock is obtained as(2)E21+νui,jj+E21+ν1-2νuj,ji-αp,i+fi=0where σij is the stress tensor, Pa; σij′ is the effective stress tensor, Pa; tension is taken as positive; fi is the body force in solid skeleton; α is the Biot coefficient; p is pore pressure, Pa; u is displacement, m; E is Young’s modulus, GPa; ν is Poisson’s ratio.The mass balance equation of fluid flow in deformable porous media is written as(3)S∂p∂t-∇·q+α∂εv∂t=QfAccording to Darcy’s law for fluid flow, the flow rateq is described as(4)q=-kμ∇p+ρfg∇Zwhere S is the storage coefficient, Pa−1; εv is the volumetric strain of rock mass; Qf is the source-sink term of the seepage process, s−1; k is rock permeability, μm2; μ is fluid viscosity, Pa∙s; ρf is fluid density, kg/m3; Z is elevation, m.Domain of the problem is illustrated as in Figure1. Consider a porous continuum in the domain Ω bounded by surface Γ, the initial conditions and boundary conditions are stated as follows.Figure 1
Schematic of the problem domain.(i) Initial conditions(5)ui0=u0,p0=p0inΩ(ii) Boundary conditions(6)ui=u-tonΓuσij·n=TitonΓσp=p-tonΓp-kμ∇p=q-tonΓqwhere u0, p0 are the displacement and pressure of initial domain, u-(t), p-(t) are the displacement and pressure of boundary, Ti(t) is the boundary force, n is the unit normal vector perpendicular to the boundary surface, q-(t) is the fluid flux, Γu, Γσ, Γp, and Γq 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].(7)Fσ1=σ1-ft0where F(σ1) is the yield function of tensile failure, Pa; σ1 is the maximum principal stress, Pa; ft0 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:(8)E=1-DE0where E and E0 represent the Young’s modulus of the damaged elements and undamaged elements and D is the damage variable.Damage evolution of the element in tensional state can be expressed as [26] (9)D=0εt<εt01-λεt0εtεt≥εt0where εt=∑i=13εi+εi/22 is the equivalent principal strain for a tensile state, εi (i=1,2,3) are the three principal strains, εt0 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(10)E1-D21+νui,jj+E1-D21+ν1-2νuj,ji-αp,i+fi=0Because of the damage, permeability is correlated to damage effect within the failure elements. The permeability can be described as [23]: (11)k=ξk0exp⁡-ησ-′0<D≤1where k0 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:(12)Ku+Qp=fu(13)Sp˙+Hp+QTu˙=fqwhere u, p are the vectors of displacement and pressure and u˙, p˙ are corresponding time derivatives. The matrix expressions are listed as follows.(14)K=∫ΩBTDBdΩQ=∫ΩαBTmNpdΩfu=∫ΩNuTfdΩ+∫ΓσNuTTdΓS=∫ΩNpTSNpdΩH=∫Ω∇NpTku∇NpdΩfp=∫Γqe-NpTq-dΓ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 (15)Kut+△t=Kut+QΔpt+△fu(16)S+θΔtHpt+△t=S-1-θΔtHpt-QTΔut+Δtfqwhere ut+△t, pt+△t represent the solution at the current time step, ut, pt is the solution at last time step, θ is an integration parameter in the interval 0≤θ≤1. In the following numerical example, we choose θ=1.The displacement convergence criterion is used to judge the convergence of incremental solution(17)△ε=∑i=1node△ut2ut2≤10-5where △ε is the iterative tolerance of the displacement increment, △ut=ut+Δt-ut is the displacement increment at current time step, and ut 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 pressurep is computed first from (16). Then the displacement, strain, and stress are computed from (15) from updated p 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 K and H 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.Figure 2
Flow chart of the iterative calculations in coupled FE model.


### 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 p0 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 t=0. For solution algorithm verification, we choose Youth’s modulus E=1.0 GPa, Possion’s ration ν=0.25, permeability k=1.0 um2, fluid viscosity μ=1.0 × 10-3Pa∙s, Biot coefficient α=1.0, and initial pressure p0=1.0 MPa.Figure 3
Geometry for Terzaghi’s consolidation problem.As is shown in Figure4, 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.Figure 4
Comparison of numerical results and analytical solution of Terzaghi’s problem. (a) Pore pressure and (b) vertical displacement at the top.
(a) (b)


## 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:(12)Ku+Qp=fu(13)Sp˙+Hp+QTu˙=fqwhere u, p are the vectors of displacement and pressure and u˙, p˙ are corresponding time derivatives. The matrix expressions are listed as follows.(14)K=∫ΩBTDBdΩQ=∫ΩαBTmNpdΩfu=∫ΩNuTfdΩ+∫ΓσNuTTdΓS=∫ΩNpTSNpdΩH=∫Ω∇NpTku∇NpdΩfp=∫Γqe-NpTq-dΓ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 (15)Kut+△t=Kut+QΔpt+△fu(16)S+θΔtHpt+△t=S-1-θΔtHpt-QTΔut+Δtfqwhere ut+△t, pt+△t represent the solution at the current time step, ut, pt is the solution at last time step, θ is an integration parameter in the interval 0≤θ≤1. In the following numerical example, we choose θ=1.The displacement convergence criterion is used to judge the convergence of incremental solution(17)△ε=∑i=1node△ut2ut2≤10-5where △ε is the iterative tolerance of the displacement increment, △ut=ut+Δt-ut is the displacement increment at current time step, and ut 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 pressurep is computed first from (16). Then the displacement, strain, and stress are computed from (15) from updated p 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 K and H 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.Figure 2
Flow chart of the iterative calculations in coupled FE model.


## 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 p0 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 t=0. For solution algorithm verification, we choose Youth’s modulus E=1.0 GPa, Possion’s ration ν=0.25, permeability k=1.0 um2, fluid viscosity μ=1.0 × 10-3Pa∙s, Biot coefficient α=1.0, and initial pressure p0=1.0 MPa.Figure 3
Geometry for Terzaghi’s consolidation problem.As is shown in Figure4, 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.Figure 4
Comparison of numerical results and analytical solution of Terzaghi’s problem. (a) Pore pressure and (b) vertical displacement at the top.
(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σh 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 σH. 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σv=43.0 MPa, maximum horizontal stress σH=38.0 MPa, and minimum horizontal stress σh=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 p=18.0MPa is applied as the initial condition. No-flow boundary condition is applied on all sides.Figure 5
Schematic of the model. (a) Model geometry, (b) spiral perforation model, (c) wellbore angle diagram, and (d) perforation angle diagram.For an open-hole horizontal well, it was assumed that the same bottom-hole pressure (BHP)pp loads were applied on the wellbore and perforation surfaces. For a cased horizontal well, the BHP pp 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 pw applied to the borehole boundary in a cased-and-cemented wellbore has been proposed [8, 13]:(18)pw=TFpp-cσh+cσhTF=1+vs/Es21-vs/Ro2-Ri2Ri21+v/E+1+vs/EsRi2+1-2vsRo2/Ro2-Ri2·Ro2rw2where Es, vs are Young’s modulus and Poisson’s ratio of the steel, TF is the transmission factor, c is the in-situ stress factor, Ri, Ro are the inner and outer radii of the casing, and rw is the inner radius of wellbore.The reservoir properties and well parameters are reported in Table1.Table 1
Reservoir properties and well parameters.
Parameters Values Young’s modulus of rockE0 (GPa) 21.5 Poisson’s ratio of rockν 0.23 Biot coefficientα 0.65 Permeabilityk (m2) 3.0×10−15 Tensile strengthft0(MPa) 3.0 Residual strength coefficientλ 0.15 Mutation coefficientξ 10.0 Coupling coefficientη 0.025 Wellbore diameterrw(mm) 178 External diameter of the casingRo(mm) 178 Internal diameter of the casingRi(mm) 164 Young’s modulus of casingEs(GPa) 210 Poisson’s ratio of casingvs 0.21 In situ stress factorc 0.3Figure6 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 (pf) 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 p-(t) 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.Figure 6
Typical fracturing curve of the fracturing treatment.


### 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 Figure7. 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%.Figure 7
Comparison of tangential stress in open-hole horizontal well. (a) Tangential stress along the wellbore and (b) tangential stress along the perforation tunnel.
(a) (b)Figure7(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%.Figure8 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.Figure 8
Comparison of fracture geometry in near wellbore of spiral perforation. (a) Fracture geometry in the FEM model and (b) fracture geometry in experiment [17].
(a) (b)


## 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σh 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 σH. 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σv=43.0 MPa, maximum horizontal stress σH=38.0 MPa, and minimum horizontal stress σh=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 p=18.0MPa is applied as the initial condition. No-flow boundary condition is applied on all sides.Figure 5
Schematic of the model. (a) Model geometry, (b) spiral perforation model, (c) wellbore angle diagram, and (d) perforation angle diagram.For an open-hole horizontal well, it was assumed that the same bottom-hole pressure (BHP)pp loads were applied on the wellbore and perforation surfaces. For a cased horizontal well, the BHP pp 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 pw applied to the borehole boundary in a cased-and-cemented wellbore has been proposed [8, 13]:(18)pw=TFpp-cσh+cσhTF=1+vs/Es21-vs/Ro2-Ri2Ri21+v/E+1+vs/EsRi2+1-2vsRo2/Ro2-Ri2·Ro2rw2where Es, vs are Young’s modulus and Poisson’s ratio of the steel, TF is the transmission factor, c is the in-situ stress factor, Ri, Ro are the inner and outer radii of the casing, and rw is the inner radius of wellbore.The reservoir properties and well parameters are reported in Table1.Table 1
Reservoir properties and well parameters.
Parameters Values Young’s modulus of rockE0 (GPa) 21.5 Poisson’s ratio of rockν 0.23 Biot coefficientα 0.65 Permeabilityk (m2) 3.0×10−15 Tensile strengthft0(MPa) 3.0 Residual strength coefficientλ 0.15 Mutation coefficientξ 10.0 Coupling coefficientη 0.025 Wellbore diameterrw(mm) 178 External diameter of the casingRo(mm) 178 Internal diameter of the casingRi(mm) 164 Young’s modulus of casingEs(GPa) 210 Poisson’s ratio of casingvs 0.21 In situ stress factorc 0.3Figure6 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 (pf) 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 p-(t) 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.Figure 6
Typical fracturing curve of the fracturing treatment.


## 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 Figure7. 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%.Figure 7
Comparison of tangential stress in open-hole horizontal well. (a) Tangential stress along the wellbore and (b) tangential stress along the perforation tunnel.
(a) (b)Figure7(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%.Figure8 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.Figure 8
Comparison of fracture geometry in near wellbore of spiral perforation. (a) Fracture geometry in the FEM model and (b) fracture geometry in experiment [17].
(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 Figure9, 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)).Figure 9
Initial fracture geometry for open-hole perforation with different angleβ. (a) β=0°, (b) β=30°, (c) β=60°, and (d) β=90°.
(a) (b) (c) (d)For cased-cement horizontal well, the initiated fractures are all transverse to the wellbore (Figure10). 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)).Figure 10
Initial fracture geometry for case-cement perforation with different angleβ. (a) β=0°, (b) β=30°, (c) β=60°, and (d) β=90°.
(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 σv. The fracture initial pressure of the tunnel is mainly affected by the maximum horizontal stress σH and minimum horizontal stress σh. 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 σv and minimum horizontal stress σh.As shown in Figure11, 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 pf grows up to its maximum at β=90°, the pf difference between various perforation angles for case-cement horizontal well is 3.5MPa.Figure 11
Effects of deflection angle on perforation initiation pressure.


### 5.2. The Initial Fracture Characters of Horizontal Wells

Figure12 presents the maximum stress principal change of spiral perforations during the pressurization period. With the BHP pp 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)).Figure 12
The maximum principal stress change of spiral perforation tunnels.Figure 13
Damage evolution of spiral perforation during the pressurization period. (a)pp=41.5 MPa, (b) pp =43.5 MPa, and (c) pp =45.0 MPa.
(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. Aspp 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 (Figures8(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 Figure14. 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 pp=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.Figure 14
The maximum principal stress change of in-plane perforation tunnels.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.Figure 15
Damage evolution of in-plane perforation during the pressurization period. (a)pp =39.5MPa, (b) pp =41.5MPa, and (c) pp =45.0MPa.
(a) (b) (c)


### 5.3. The Near-Wellbore Fracture Geometry of Horizontal Wells

Figure16 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.Figure 16
Comparison of near-wellbore fracture geometry with different perforation strategies. (a) Open-hole horizontal well, (b) case-cement horizontal well, and (c) in-plane perforation horizontal well.
(a) (b) (c)Figure16(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.


## 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 Figure9, 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)).Figure 9
Initial fracture geometry for open-hole perforation with different angleβ. (a) β=0°, (b) β=30°, (c) β=60°, and (d) β=90°.
(a) (b) (c) (d)For cased-cement horizontal well, the initiated fractures are all transverse to the wellbore (Figure10). 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)).Figure 10
Initial fracture geometry for case-cement perforation with different angleβ. (a) β=0°, (b) β=30°, (c) β=60°, and (d) β=90°.
(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 σv. The fracture initial pressure of the tunnel is mainly affected by the maximum horizontal stress σH and minimum horizontal stress σh. 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 σv and minimum horizontal stress σh.As shown in Figure11, 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 pf grows up to its maximum at β=90°, the pf difference between various perforation angles for case-cement horizontal well is 3.5MPa.Figure 11
Effects of deflection angle on perforation initiation pressure.


## 5.2. The Initial Fracture Characters of Horizontal Wells

Figure12 presents the maximum stress principal change of spiral perforations during the pressurization period. With the BHP pp 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)).Figure 12
The maximum principal stress change of spiral perforation tunnels.Figure 13
Damage evolution of spiral perforation during the pressurization period. (a)pp=41.5 MPa, (b) pp =43.5 MPa, and (c) pp =45.0 MPa.
(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. Aspp 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 (Figures8(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 Figure14. 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 pp=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.Figure 14
The maximum principal stress change of in-plane perforation tunnels.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.Figure 15
Damage evolution of in-plane perforation during the pressurization period. (a)pp =39.5MPa, (b) pp =41.5MPa, and (c) pp =45.0MPa.
(a) (b) (c)


## 5.3. The Near-Wellbore Fracture Geometry of Horizontal Wells

Figure16 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.Figure 16
Comparison of near-wellbore fracture geometry with different perforation strategies. (a) Open-hole horizontal well, (b) case-cement horizontal well, and (c) in-plane perforation horizontal well.
(a) (b) (c)Figure16(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.



---

*Source: 9304961-2019-05-27.xml*