The Impact of Oriented Perforations on Fracture Propagation and Complexity in Hydraulic Fracturing

To better understand the interaction between hydraulic fracture and oriented perforation, a fully coupled finite element method (FEM)-based hydraulic-geomechanical fracture model accommodating gas sorption and damage has been developed. Damage conforms to a maximum stress criterion in tension and to Mohr–Coulomb limits in shear with heterogeneity represented by a Weibull distribution. Fracturing fluid flow, rock deformation and damage, and fracture propagation are collectively represented to study the complexity of hydraulic fracture initiation with perforations present in the near-wellbore region. The model is rigorously validated against experimental observations replicating failure stresses and styles during uniaxial compression and then hydraulic fracturing. The influences of perforation angle, in situ stress state, initial pore pressure, and properties of the fracturing fluid are fully explored. The numerical results show good agreement with experimental observations and the main features of the hydraulic fracturing process in heterogeneous rock are successfully captured. A larger perforation azimuth (angle) from the direction of the maximum principal stress induces a relatively larger curvature of the fracture during hydraulic fracture reorientation. Hydraulic fractures do not always initiate at the oriented perforations and the fractures induced in hydraulic fracturing are not always even and regular. Hydraulic fractures would initiate both around the wellbore and the oriented perforations when the perforation angle is >75◦. For the liquid-based hydraulic fracturing, the critical perforation angle increases from 70◦ to 80◦, with an increase in liquid viscosity from 10−3 Pa·s to 1 Pa·s. While for the gas fracturing, the critical perforation angle remains 62◦ to 63◦. This study is of great significance in further understanding the near-wellbore impacts on hydraulic fracture propagation and complexity.


Introduction
Hydraulic fracturing is commonly used to improve productivity from hydrocarbon reservoirs and has been a key technique in accessing unconventional shale reservoirs.As such, it is recognized as one of the most effective stimulation techniques for enhanced recovery from reservoirs [1,2].Hydraulic fracturing is used not only to create macro-scale fractures but is also complicit in connecting and reactivating pre-existing fractures at all scales to create transport pathways within reservoirs.The transport linkage between preexisting flaws (i.e., perforations, joints, and natural fractures) are often critical to hydrocarbon production from reservoirs.These natural fractures or other flaws can induce complex fracture geometries and flow pathways resulting from both tensile and shear failure [3][4][5].Implicit in this is understanding fracture initiation and propagation in the near-wellbore region as of great significance for hydraulic fracture treatment and injectivity test interpretation [6][7][8].
Since most wells are cased and require to be perforated to access the reservoir, the majority of hydraulic fracturing treatments are conducted through perforations.Perforations play a vital role in the origin of complex fracture geometries [6,[9][10][11][12].In fractured reservoirs, perforations are often the sole hydraulic connection between the fractured reservoir and the wellbore.Hydraulic fractures initiated from preexisting perforations will reorient themselves normal to the direction of the minimum far-field stress as they propagate out from the wellbore [2,4,[13][14][15][16].Near the wellbore, multiple sub-parallel fractures may be generated due to fracture initiation processes such as perforations and natural flaws located at the open-hole wellbore wall.Hydraulic fractures commonly initiate either from perforations located at the cased wellbore or from preexisting natural fractures intersecting an uncased wellbore [17].However, when the perforation angle and the horizontal stress difference are large, dual fractures may be induced along perforations and the maximum stress direction.This means that the fractures formed by hydraulic fracturing from orientated perforations may be very complex.With regard to moving boundary problems in hydraulic fracturing, many model reduction techniques have been proposed to achieve the desired accuracy and computational efficiency [18][19][20].Recently, several efforts have been made to determine an optimal perforation condition utilizing the integrated reduced-order model [21] and 3D fracture-propagation model [22].Near-wellbore effects contain the fracture curvature, pressure enhancement in perforation hole and fracture branching, which are crucial issues regarding stimulated reservoir volume thus production.Therefore, in this paper, phenomena on the hydraulic fracture curving and branching near wellbore are focused.Furthermore, a finite element method (FEM)-based damage model is proposed to simulate the complete evolution of hydraulic fractures initiation and propagation.The impacts of oriented perforations on the near-wellbore fracture complexity are examined in terms of fracture trajectory, horizontal differential stress, fracturing fluid viscosity, and wellbore pressure.
Mechanical and transport properties of the reservoir may be significantly altered during stimulation [23][24][25][26][27][28].In consideration of mentioned above, the impact of perforations and pre-existing features in reservoirs on the evolution of fluid-driven fractures, their interaction with pre-existing features and the resulting complexity of the evolving fracture network are investigated.This response is explored through a novel fully coupled model of rock damage and fluid flow for both slightly compressible and compressible fluids.Damage conforms to a maximum stress criterion in tension and to Mohr-Coulomb limits in shear with heterogeneity represented by a Weibull distribution.Fracturing fluid flow, rock deformation and damage, and fracture propagation are collectively represented to study the complexity of hydraulic fracture initiation with perforations present in the near-wellbore region.Specifically, the main objective of this study is thus to investigate the effects of preexisting perforation angle, horizontal differential stress, fracturing fluid, and initial pore pressure on fracture propagation and complexity.Research findings are applied to optimize perforation design in order to minimize high treatment pressures and to control complexity in the near-wellbore fracture region.

Conceptual Model
The near-wellbore fractures initiated from two symmetrically oriented perforations considered in this study are shown in Figure 1.The 2D hydraulic fracture model considers a cross-section through a vertical well with its axis aligned in the direction of maximum or intermediate principal stress.These two perforations are represented by two preexisting fractures (2D perforations rather than tunnels), which serve as initiation sites for the propagation of longitudinal hydraulic fractures.On the condition that the perforated interval is long, or the well is not cemented and cased, these hydraulic fractures and their propagation within the transverse plane are assumed to propagate under plane-strain conditions [17,29].After initiation, and beyond the influence of the wellbore, the hydraulic fractures will reorient themselves from the perforation to the orientation of the minimum principal stress.However, when the perforation angles are large, more than two hydraulic fractures may initiate in a homogeneous rock instead of just a single zigzag shape fracture, as described in Figure 1.
Processes 2018, 6, x FOR PEER REVIEW 3 of 20 principal stress.However, when the perforation angles are large, more than two hydraulic fractures may initiate in a homogeneous rock instead of just a single zigzag shape fracture, as described in Figure 1.The in situ stress field is defined in terms of the in-plane maximum horizontal stress σH and minimum horizontal stress σh.The preexisting perforation length L is assumed to be very small relative to the wellbore radius r-therefore the impact of preexisting perforation length on hydraulic fracture path can be neglected.The perforation is defined by its azimuth θ with respect to the orientation of maximum horizontal stress σH (x-axis).If the perforation tunnel is oriented in the direction of the maximum principal stress ( 0 θ = ), the hydraulic fractures would propagate directly in this (x-axis) direction.When the tunnels are oblique to this direction ( 0 θ ≠ ) the hydraulic fractures must ultimately reorient themselves parallel to the direction of the maximum horizontal stress when they exit the stress shadows of the near-wellbore region and of the perforation.Beyond such stress-shadows, the hydraulic fractures tend to propagate in the plane of least resistance perpendicular to the minimum horizontal stress.

Formulation of the Conceptual Model
During hydraulic fracturing, fluids will penetrate into the porous medium surrounding both the borehole and the perforation tunnel.Thus, the evolution of the fracture(s) then results from this coupled phenomenon of fluid flow, solid deformation, and damage.The governing equations for this fully coupled hydrological-mechanical damage behavior were developed based on conservation of mass, momentum, and energy.The effects of rock damage on the evolving mechanical and transport parameters are also considered.

Governing Equations for Mechanical Response
The poroelastic constitutive relation for rock deformation is given by Reference [30]: where G is the equivalent shear modulus (Pa), K is the bulk modulus (Pa), p is the fluid pressure (Pa), ij δ is the Kronecker delta, α is the Biot coefficient, and s ε is the adsorption strain (m/m).
The gas adsorption strain may be represented by the Langmuir isotherm, as follows [31]: The in situ stress field is defined in terms of the in-plane maximum horizontal stress σ H and minimum horizontal stress σ h .The preexisting perforation length L is assumed to be very small relative to the wellbore radius r-therefore the impact of preexisting perforation length on hydraulic fracture path can be neglected.The perforation is defined by its azimuth θ with respect to the orientation of maximum horizontal stress σ H (x-axis).If the perforation tunnel is oriented in the direction of the maximum principal stress (θ = 0), the hydraulic fractures would propagate directly in this (x-axis) direction.When the tunnels are oblique to this direction (θ = 0) the hydraulic fractures must ultimately reorient themselves parallel to the direction of the maximum horizontal stress when they exit the stress shadows of the near-wellbore region and of the perforation.Beyond such stress-shadows, the hydraulic fractures tend to propagate in the plane of least resistance perpendicular to the minimum horizontal stress.

Formulation of the Conceptual Model
During hydraulic fracturing, fluids will penetrate into the porous medium surrounding both the borehole and the perforation tunnel.Thus, the evolution of the fracture(s) then results from this coupled phenomenon of fluid flow, solid deformation, and damage.The governing equations for this fully coupled hydrological-mechanical damage behavior were developed based on conservation of mass, momentum, and energy.The effects of rock damage on the evolving mechanical and transport parameters are also considered.

Governing Equations for Mechanical Response
The poroelastic constitutive relation for rock deformation is given by Reference [30]: where G is the equivalent shear modulus (Pa), K is the bulk modulus (Pa), p is the fluid pressure (Pa), δ ij is the Kronecker delta, α is the Biot coefficient, and ε s is the adsorption strain (m/m).The gas adsorption strain may be represented by the Langmuir isotherm, as follows [31]: where ε L is the Langmuir strain constant (m/m) and P L is the Langmuir pressure constant (Pa).From this, a modified Navier-type constitutive deformation [32], in terms of displacement u i (m) and fluid pressure p (Pa), can be expressed as: where the effects of fluid pressure and injection rate on the rock deformation are defined as Equation (3).

Governing Equations for Fluid Flow
For the purpose to study the effect of different fracturing fluids on fracture complexity, fluid flow for both slightly compressible (liquids) and compressible (gases) fluids are presented, respectively.

Slightly Compressible Fluids
In general, the equation of state is defined with regard to the fluid compressibility, c f (Pa −1 ), as: Assuming that the fluid compressibility c f is constant over a prescribed range of pressure, after integration, Equation (4) can be written as: where ρ 0 is the density (kg•m −3 ) at the reference pressure p 0 (Pa).According to Taylor series expansion, this may be expressed as: with an approximate result given by: Similarly, the rock compressibility may be defined as: and again, after integration, this may be written as: where φ 0 is the porosity at p 0 .Similarly, this may be approximated as: Neglecting the effect of gravity, which will be small in relation to fluid pressure gradients, fluid flow may be defined in terms of Darcy's law as: where k is the equivalent permeability of the rock (m 2 ) and µ f is the fluid viscosity (Pa•s).
Processes 2018, 6, 213 5 of 19 According to the mass conservation equation, the governing equation for fluid transport in a porous medium is defined as: where m is fluid mass, ρ is fluid density (kg•m −3 ), q is the Darcy velocity vector and Q s is the source or sink of the compressible fluid defined in terms of mass rate (kg•m −3 •s −1 ).Supposing that the rock is fluid-saturated, the fluid content per unit volume can be expressed as m = ρφ, where φ is the porosity.For slightly compressible fluids the partial derivative of m with respect to time can be written as: The total compressibility may be defined as: Combining Equations ( 4)-( 14) yields the governing equation for slightly compressible fluid flow in the fracture network as: enabling solution for evolving flux or pressure boundary conditions and the evolution of permeability with damage.

Compressible Fluids
For compressible fluid (gas) the compressibility c g is both much larger than for a slightly compressible fluid (liquid) and variable with pressure.In such a case, the pressure-volume-temperature (PVT) relation for a non-ideal or real gas can be written as: where M is the molecular weight (kg/mol), Z is the gas compressibility factor (Pa −1 ), R is the universal gas constant (J•mol −1 •K −1 ), and T is the temperature (K).Assuming that the gas compressibility and viscosity are constant, then the governing relation is defined as: accommodating the pressure dependent compressibility of the fluid.

Governing Equations Accommodating Rock Heterogeneity and Damage Evolution
To characterize heterogeneity at the mesoscopic-level, mechanical parameters such as strength, Young's modulus, and Poisson's ratio may be assigned according to the Weibull distribution.This distribution defines the parameter according to a probability density function as [33]: where u is the parameter of interest (such as strength, Young's modulus, or Poisson's ratio); the scale parameter u 0 is related to the average of the element parameter, and m is the shape parameter of the Weibull distribution function.m is defined as the degree of rock homogeneity and called the homogeneity index.The variations of f (u) with respect to m are shown in Figure 2. Obviously, a higher m value represents a more homogeneous rock.
Processes 2018, 6, x FOR PEER REVIEW 6 of 20 where u is the parameter of interest (such as strength, Young's modulus, or Poisson's ratio); the scale parameter u0 is related to the average of the element parameter, and m is the shape parameter of the Weibull distribution function.m is defined as the degree of rock homogeneity and called the homogeneity index.The variations of ( ) f u with respect to m are shown in Figure 2. Obviously, a higher m value represents a more homogeneous rock.As illustrated in Figure 3, rock damage in tension or shear is initiated when its state of stress (positive for compression) satisfies the maximum tensile stress criterion or the Mohr-Coulomb criterion, respectively as [27,33]: where ft0 and fc0 are uniaxial tensile and compressive strength (Pa), respectively, 1 σ and 3 σ are first and third principal stresses (Pa), respectively, θ is the internal frictional angle (°), and F1 and F2 are two damage threshold functions (Pa).As illustrated in Figure 3, rock damage in tension or shear is initiated when its state of stress (positive for compression) satisfies the maximum tensile stress criterion or the Mohr-Coulomb criterion, respectively as [27,33]: where f t0 and f c0 are uniaxial tensile and compressive strength (Pa), respectively, σ 1 and σ 3 are first and third principal stresses (Pa), respectively, θ is the internal frictional angle ( • ), and F 1 and F 2 are two damage threshold functions (Pa).
Based on the theory of elastic damage, the elastic modulus of the damaged rock is expressed as [33]: where D represents the damage variable, and E and E 0 are the Young's modulus (Pa) of the damaged and the undamaged element, respectively.In the present study, the element, as well as its damage, is assumed isotropic, so the E, E 0 and D parameters are all scalars.Under any stress and initial conditions, the tensile stress criterion is applied preferentially.In other words, the maximum tensile stress criterion is first applied to judge whether the elements are first damaged in tension or not.Only elements that survive this test will be checked for damage in shear using the Mohr-Coulomb criterion.Based on the theory of elastic damage, the elastic modulus of the damaged rock is expressed as [33]: where D represents the damage variable, and E and E0 are the Young's modulus (Pa) of the damaged and the undamaged element, respectively.In the present study, the element, as well as its damage, is assumed isotropic, so the E, E0 and D parameters are all scalars.Under any stress and initial conditions, the tensile stress criterion is applied preferentially.In other words, the maximum tensile stress criterion is first applied to judge whether the elements are first damaged in tension or not.Only elements that survive this test will be checked for damage in shear using the Mohr-Coulomb criterion.
In terms of the damage constitutive law shown in Figure 3, the damage variable for the rock can be calculated as [33]: the maximum tensile and maximum compressive principal strains (m/m) when tensile and shear damage occurs, respectively, and n is a constitutive coefficient specified as 2.0.In Equation ( 19), when F1 < 0 and F2 < 0 the applied stress is insufficient to satisfy the maximum tensile stress criterion and the Mohr-Coulomb failure criterion, respectively.F1 = 0 and dF1 > 0 implies rock damage in the tensile mode when the stress state satisfies the maximum tensile stress criterion and the rock is still under load.F2 = 0 and dF2 > 0 implies rock damage in the shear mode when the stress state satisfies the Mohr-Coulomb failure criterion and the rock remains loaded.In terms of the damage constitutive law shown in Figure 3, the damage variable for the rock can be calculated as [33]: where ε 1 and ε 3 are the major and minor principal strains (m/m), respectively.ε t0 and ε c0 are the maximum tensile and maximum compressive principal strains (m/m) when tensile and shear damage occurs, respectively, and n is a constitutive coefficient specified as 2.0.In Equation (19), when F 1 < 0 and F 2 < 0 the applied stress is insufficient to satisfy the maximum tensile stress criterion and the Mohr-Coulomb failure criterion, respectively.F 1 = 0 and dF 1 > 0 implies rock damage in the tensile mode when the stress state satisfies the maximum tensile stress criterion and the rock is still under load.F 2 = 0 and dF 2 > 0 implies rock damage in the shear mode when the stress state satisfies the Mohr-Coulomb failure criterion and the rock remains loaded.The damage variables of Equation ( 21) remain in the range from 0 to 1.0 regardless of the form or magnitude of damage.However, negative and positive damage magnitudes are respectively adapted for damage in tension and in shear, merely to allow visualization of the two damage modes in post-processed figures.In this regard, tensile damage is represented as negative numbers (−1 ≤ D < 0), while shear damage is represented as positive numbers (0 < D ≤ 1).In the respect of the irreversibility of damage, the damage variable may only increment monotonically from zero during loading (dF 1 > 0 or dF 2 > 0) and remain unchanged during unloading (dF 1 < 0 or dF 2 < 0).In this respect, the damage, defined by Equation (20), reduces the elastic modulus E and the shear modulus G of the rock, via to Equations ( 20) and (21).

Model Validation against Experimental Observations
In this section, the hydraulic fracturing geomechanical model proposed in this study is validated against observed mechanisms of uniaxial compression and hydraulic fracturing [11].

Comparisons of Breakdown Pressure and Fracture Geometry
The relationship between axial stress and axial strain during modeled rupture of a uniaxial specimen is shown in Figure 4 relative to the mechanical parameters recovered from experiment [11].Apparent from Table 1 and Figure 4 is that the experimental data adequately replicate the anticipated stress-strain behavior.Additionally, it should be noted that the simulation of uniaxial compression is conducted to provide basic mechanical properties to inform the numerical simulation of hydraulic fracturing conducted in the same material.A visual comparison between experimental results [11] and the numerical simulations is shown in Figure 5.As can be seen from Figure 5, the fracture pattern recovered from the numerical model is consistent with experimental observations.Hydraulic fractures initiate outward along the perforation and as they propagate away from the near-wellbore towards the region of unaltered in-situ stress they reorient themselves perpendicular to the minimum horizontal stress.
modulus G of the rock, via to Equations ( 20) and (21).

Model Validation against Experimental Observations
In this section, the hydraulic fracturing geomechanical model proposed in this study is validated against observed mechanisms of uniaxial compression and hydraulic fracturing [11].

Comparisons of Breakdown Pressure and Fracture Geometry
The relationship between axial stress and axial strain during modeled rupture of a uniaxial specimen is shown in Figure 4 relative to the mechanical parameters recovered from experiment [11].Apparent from Table 1 and Figure 4 is that the experimental data adequately replicate the anticipated stress-strain behavior.Additionally, it should be noted that the simulation of uniaxial compression is conducted to provide basic mechanical properties to inform the numerical simulation of hydraulic fracturing conducted in the same material.A visual comparison between experimental results [11] and the numerical simulations is shown in Figure 5.As can be seen from Figure 5, the fracture pattern recovered from the numerical model is consistent with experimental observations.Hydraulic fractures initiate outward along the perforation and as they propagate away from the near-wellbore towards the region of unaltered in-situ stress they reorient themselves perpendicular to the minimum horizontal stress.

Geometric Model and Boundary Conditions
The geometry and boundary conditions of the specimen are simplified as shown in Figure 6.The numerical specimen, 300 mm on edge with a 10 mm diameter borehole at the center is subjected to a biaxial horizontal stress.The boundary conditions correspond to a confining pressure σH applied on the top boundary and σh applied on the right boundary with rollers applied along both the left side and the base.There are no-flux conditions on all boundaries except for the borehole wall into which a constant fluid injection rate is applied.The fluid injection rate in the experiments is 2.1 × 10 −9 m 3 •s −1 .The numerical mechanical properties for the simulation are listed in Table 2.

Geometric Model and Boundary Conditions
The geometry and boundary conditions of the specimen are simplified as shown in Figure 6.The numerical specimen, 300 mm on edge with a 10 mm diameter borehole at the center is subjected to a biaxial horizontal stress.The boundary conditions correspond to a confining pressure σ H applied on the top boundary and σ h applied on the right boundary with rollers applied along both the left side and the base.There are no-flux conditions on all boundaries except for the borehole wall into which a constant fluid injection rate is applied.The fluid injection rate in the experiments is 2.1 × 10 −9 m 3 •s −1 .The numerical mechanical properties for the simulation are listed in Table 2.

Geometric Model and Boundary Conditions
The geometry and boundary conditions of the specimen are simplified as shown in Figure 6.The numerical specimen, 300 mm on edge with a 10 mm diameter borehole at the center is subjected to a biaxial horizontal stress.The boundary conditions correspond to a confining pressure σH applied on the top boundary and σh applied on the right boundary with rollers applied along both the left side and the base.There are no-flux conditions on all boundaries except for the borehole wall into which a constant fluid injection rate is applied.The fluid injection rate in the experiments is 2.1 × 10 −9 m 3 •s −1 .The numerical mechanical properties for the simulation are listed in Table 2.

Effects of Preexisting Perforation Orientation
Perforation orientation has a significant effect on the geometry of the hydraulic fracture initiated from the wellbore.The azimuth and length of the preexisting perforation are two important factors that affect hydraulic fracture propagation and pattern.An optimal preexisting perforation should initiate only a single wing fracture with minimum tortuosity at an achievable breakdown pressure [9,34].Numerical simulations for cases with preexisting perforation angles of 0 • , 15 • , 30 • , 45 • , 60 • , and 75 • are conducted to examine the response.As shown in Figures 7 and 8, a larger perforation angle induces a relatively larger curvature during the reorientation of the hydraulic fractures and takes a long distance for the hydraulic fracture to rotate to the direction completely aligned with the direction of the maximum horizontal stress.A larger perforation angle requires a higher breakdown pressure to initiate the hydraulic fracture.In addition, another observation is that hydraulic fractures are initiated both around the wellbore and preexisting fractures when the perforation angle is 75 • .Additionally, hydraulic fractures do not always initiate from preexisting perforations.The fractures formed by hydraulic fracturing with oriented perforations have a complex morphology.Therefore, for perforations in field operations, perforations with a high angle are generally not recommended and should be avoided as much as possible [14].Specifically, there is a significant difference in breakdown pressure between numerical and experimental results, when the perforation angle is 75 • .This may be due to experimental inconsistencies or due to heterogeneity in the sample, since it is expected that the breakdown pressure would increase with an increase of the perforation angle [9,14,17].

Effects of Preexisting Perforation Orientation
Perforation orientation has a significant effect on the geometry of the hydraulic fracture initiated from the wellbore.The azimuth and length of the preexisting perforation are two important factors that affect hydraulic fracture propagation and pattern.An optimal preexisting perforation should initiate only a single wing fracture with minimum tortuosity at an achievable breakdown pressure [9,34].Numerical simulations for cases with preexisting perforation angles of 0°, 15°, 30°, 45°, 60°, and 75° are conducted to examine the response.As shown in Figures 7 and 8, a larger perforation angle induces a relatively larger curvature during the reorientation of the hydraulic fractures and takes a long distance for the hydraulic fracture to rotate to the direction completely aligned with the direction of the maximum horizontal stress.A larger perforation angle requires a higher breakdown pressure to initiate the hydraulic fracture.In addition, another observation is that hydraulic fractures are initiated both around the wellbore and preexisting fractures when the perforation angle is 75°.Additionally, hydraulic fractures do not always initiate from preexisting perforations.The fractures formed by hydraulic fracturing with oriented perforations have a complex morphology.Therefore, for perforations in field operations, perforations with a high angle are generally not recommended and should be avoided as much as possible [14].Specifically, there is a significant difference in breakdown pressure between numerical and experimental results, when the perforation angle is 75°.This may be due to experimental inconsistencies or due to heterogeneity in the sample, since it is expected that the breakdown pressure would increase with an increase of the perforation angle [9,14,17].

Analysis of Near-Wellbore Hydraulic Fracture Complexity
In situ stress, initial pore pressure, and fluid viscosity are important factors that affect the tortuosity and extension pressure of near-wellbore hydraulic fractures.In order to fully understand hydraulic fracture behavior in the field, sensitivity studies are conducted to investigate the effects of in situ horizontal differential stress, initial pore pressure, and fracturing fluid viscosity.

Effects of Horizontal Differential Stress
In order to investigate the effects of different in situ stresses on hydraulic fracture behavior, stress ratios λ of 1, 1.2, 1.4, 1.6, 1.8, and 2 are selected (stress ratios define the ratio between maximum horizontal stress and minimum horizontal stress).The magnitude of σh is held constant at 20 MPa.Stress ratio has a significant effect on both hydraulic facture propagation and reorientation, as shown in Figure 9, where different stress ratios correspond to different fracture morphologies.When the stress ratio is equal to unity (hydrostatic) the hydraulic fractures initiate and propagate along the direction of the oriented perforation.As the stress ratio λ is proportional to the difference between σH and σh, the larger this difference, the larger the propensity for the hydraulic fractures to reorient themselves to the maximum horizontal stress direction [35].A larger stress ratio results in both a smaller curvature during reorientation and occurs at a shorter distance for the fractures to reorient their direction to align with the direction of the maximum horizontal stress.Moreover, based on the results presented in Figure 10 and Table 3, it appears that the stress ratio also has a significant influence on the initiation pressure and breakdown pressure during hydraulic fracturing.It is apparent that there are clear reductions both in the initiation pressure and the breakdown pressure with increasing stress ratio.

Analysis of Near-Wellbore Hydraulic Fracture Complexity
In situ stress, initial pore pressure, and fluid viscosity are important factors that affect the tortuosity and extension pressure of near-wellbore hydraulic fractures.In order to fully understand hydraulic fracture behavior in the field, sensitivity studies are conducted to investigate the effects of in situ horizontal differential stress, initial pore pressure, and fracturing fluid viscosity.

Effects of Horizontal Differential Stress
In order to investigate the effects of different in situ stresses on hydraulic fracture behavior, stress ratios λ of 1, 1.2, 1.4, 1.6, 1.8, and 2 are selected (stress ratios define the ratio between maximum horizontal stress and minimum horizontal stress).The magnitude of σ h is held constant at 20 MPa.Stress ratio has a significant effect on both hydraulic facture propagation and reorientation, as shown in Figure 9, where different stress ratios correspond to different fracture morphologies.When the stress ratio is equal to unity (hydrostatic) the hydraulic fractures initiate and propagate along the direction of the oriented perforation.As the stress ratio λ is proportional to the difference between σ H and σ h , the larger this difference, the larger the propensity for the hydraulic fractures to reorient themselves to the maximum horizontal stress direction [35].A larger stress ratio results in both a smaller curvature during reorientation and occurs at a shorter distance for the fractures to reorient their direction to align with the direction of the maximum horizontal stress.Moreover, based on the results presented in Figure 10 and Table 3, it appears that the stress ratio also has a significant influence on the initiation pressure and breakdown pressure during hydraulic fracturing.It is apparent that there are clear reductions both in the initiation pressure and the breakdown pressure with increasing stress ratio.As mentioned above, a larger perforation angle results in a complicated fracture morphology, which hydraulic fractures are initiated both from perforations and wellbore, as shown in Figure 7 (perforation angle equals 75 • ).In order to avoid to forming such complicated fracture morphology, the critical perforation angle is defined, which is used as the minimum perforation angle to form complex fracture morphology.A series of numerical simulations have been conducted to identify the changes of critical perforation angle with stress ratios and fracturing fluids.The results of critical perforation angle against stress ratio are shown in Figure 11.Under the condition of hydrostatic pressure, the hydraulic fractures initiate and propagate along the direction of the oriented perforation, and thus will not form complex fracture morphology.The larger the stress ratio is, the less the critical perforation angle is, and more easily to form complicated fracture morphology.Such complex fracture morphology can result in difficult transportation of proppants and a significant increase in the required injection pressure for propagating hydraulic fractures, which are an undesirable phenomenon that should be avoided as much as possible in field operations [10,14].
the critical perforation angle is defined, which is used as the minimum perforation angle to form complex fracture morphology.A series of numerical simulations have been conducted to identify the changes of critical perforation angle with stress ratios and fracturing fluids.The results of critical perforation angle against stress ratio are shown in Figure 11.Under the condition of hydrostatic pressure, the hydraulic fractures initiate and propagate along the direction of the oriented perforation, and thus will not form complex fracture morphology.The larger the stress ratio is, the less the critical perforation angle is, and more easily to form complicated fracture morphology.Such complex fracture morphology can result in difficult transportation of proppants and a significant increase in the required injection pressure for propagating hydraulic fractures, which are an undesirable phenomenon that should be avoided as much as possible in field operations [10,14].

Effects of Initial Pore Pressure
With constant total stresses of σh (20 MPa) and σH (32 MPa), the initial pore pressure is varied from 2 MPa to 16 MPa.The purpose of these numerical simulations is to examine the effects of initial pore pressure on initiation pressure and breakdown pressure.An increase of fluid pore pressure can decrease static friction and thereby facilitate fracture propagation on favorably oriented planes in a deviatoric stress field.As shown in Figures 12 and 13, the resulting differential initiation pressure (Pc-P0) and breakdown pressure (Pb-P0) decrease with an increase in the initial pore pressure.This is consistent with the D-C (Detournay-Cheng) criterion [36].To incorporate an effective stress law into the D-C criterion, the geomechanical model proposed in this study correctly describes the relationship between breakdown pressure and the far-field stress in hydraulic fracturing [37].The effects of pore pressure on the initiation pressure are illustrated by the initiation equation based on poroelastic theory [38].The initiation pressure decreases with an increase in the initial pore pressure [39].There are two classical approaches to define initiation pressure in terms of far-field stresses [38,40].These represent behavior for both nonpenetrating injection fluids [40], and for penetrating fluids [38], with initiation pressures given by: Critical perforation angle (º)

Stress ratio
Critical perforation angle vs stress ratio

Effects of Initial Pore Pressure
With constant total stresses of σ h (20 MPa) and σ H (32 MPa), the initial pore pressure is varied from 2 MPa to 16 MPa.The purpose of these numerical simulations is to examine the effects of initial pore pressure on initiation pressure and breakdown pressure.An increase of fluid pore pressure can decrease static friction and thereby facilitate fracture propagation on favorably oriented planes in a deviatoric stress field.As shown in Figures 12 and 13, the resulting differential initiation pressure (P c -P 0 ) and breakdown pressure (P b -P 0 ) decrease with an increase in the initial pore pressure.This is consistent with the D-C (Detournay-Cheng) criterion [36].To incorporate an effective stress law into the D-C criterion, the geomechanical model proposed in this study correctly describes the relationship between breakdown pressure and the far-field stress in hydraulic fracturing [37].The effects of pore pressure on the initiation pressure are illustrated by the initiation equation based on poroelastic theory [38].The initiation pressure decreases with an increase in the initial pore pressure [39].There are two classical approaches to define initiation pressure in terms of far-field stresses [38,40].These represent behavior for both nonpenetrating injection fluids [40], and for penetrating fluids [38], with initiation pressures given by: (22) where P H-W and P H-F are the initiation pressures related to Hubbert-Willis [40] and Haimson-Fairhurst [38] equations, respectively, T is the tensile strength of the rock, and σ h and σ H are the far-field principal stresses.As 0 ≤ φ ≤ α ≤ 1 (φ is rock porosity) and 0 ≤ ν ≤ 0.5 (ν is Poisson' ratio of rock) for rock, then obtain 0 ≤ α(1 − 2ν)/(1 − ν) ≤ 1.Therefore, the initiation pressure for a penetrating fluid is always smaller than (or equal to) that for a nonpenetrating fluid.It can be seen from Figure 12 that the resulting differential initiation pressure (P c -P 0 ) during gas fracturing is significantly different from that for liquid fracturing.The obvious effect of gas fracturing is in the reduction of the initiation pressure.This is in good agreement with Equation (22).In addition, as shown in Figures 12 and 13, for water-based fracturing, the resulting P c -P 0 and P b -P 0 decrease linearly with an increase in the initial pore pressure.However, for gas fracturing, the resulting P c -P 0 and P b -P 0 indicate a nonlinear decrease with an increase in initial pore pressure.The influences of gas penetration complicate the mechanism of the fracturing process.The gas penetration not only alters the pore pressure in the reservoir, but also the gas adsorption-induced damage modifies the mechanical properties of the reservoir rock [27,41,42].According to the Langmuir adsorption isotherm, as shown in Equation ( 2), with an increase in pore pressure, the sorption capacity and volumetric strain simultaneously increase with the adsorptive strain potentially resulting in additional rock damage.This adsorption-induced damage is fully coupled to the gas fracturing model proposed in this paper.As a consequence, both the initiation pressure and breakdown pressure of gas fracturing show nonlinear decreases with increasing pore pressure.
where PH-W and PH-F are the initiation pressures related to Hubbert-Willis [40] and  equations, respectively, T is the tensile strength of the rock, and σh and σH are the far-field principal stresses.As 0 1 φ α ≤ ≤ ≤ (φ is rock porosity) and 0 0.5 ν ≤ ≤ (ν is Poisson' ratio of rock) for rock, then obtain ( ) ( ) Therefore, the initiation pressure for a penetrating fluid is always smaller than (or equal to) that for a nonpenetrating fluid.It can be seen from Figure 12 that the resulting differential initiation pressure (Pc-P0) during gas fracturing is significantly different from that for liquid fracturing.The obvious effect of gas fracturing is in the reduction of the initiation pressure.This is in good agreement with Equation ( 22).In addition, as shown in Figures 12 and 13, for water-based fracturing, the resulting Pc-P0 and Pb-P0 decrease linearly with an increase in the initial pore pressure.However, for gas fracturing, the resulting Pc-P0 and Pb-P0 indicate a nonlinear decrease with an increase in initial pore pressure.The influences of gas penetration complicate the mechanism of the fracturing process.The gas penetration not only alters the pore pressure in the reservoir, but also the gas adsorption-induced damage modifies the mechanical properties of the reservoir rock [27,41,42].According to the Langmuir adsorption isotherm, as shown in Equation ( 2), with an increase in pore pressure, the sorption capacity and volumetric strain simultaneously increase with the adsorptive strain potentially resulting in additional rock damage.This adsorption-induced damage is fully coupled to the gas fracturing model proposed in this paper.As a consequence, both the initiation pressure and breakdown pressure of gas fracturing show nonlinear decreases with increasing pore pressure.

Effects of Fracturing Fluids
The impacts of both slightly compressible (liquid) and compressible (gas) fluids in fracturing performance were examined to better understand fracture behavior during hydraulic fracturing with liquids or gases.Fluid viscosity is an important variable that influences fracture propagation -the higher the fluid viscosity the slower the local flow rate in the fracture and the higher the fluid pressure under the same injection rate [43,44].As shown in Figure 14, when the viscosity of the liquid is high the breakdown pressure is large.Conversely, as the gas has a lower viscosity, the breakdown pressure is also smaller.This fact is illustrated in Figure 14, which shows that the relationship between breakdown pressure and fluid viscosity follows a power-law function.In addition, the higher the fracturing fluid viscosity, the more complex the fracture geometry local to the wellbore.This is because a fluid with relatively higher viscosity has a relatively slower flow rate in the near-wellbore region [45].Figure 15 shows the variations of critical perforation angle with different viscosities of fracturing fluids.For the

Effects of Fracturing Fluids
The impacts of both slightly compressible (liquid) and compressible (gas) fluids in fracturing performance were examined to better understand fracture behavior during hydraulic fracturing with liquids or gases.Fluid viscosity is an important variable that influences fracture propagationthe higher the fluid viscosity the slower the local flow rate in the fracture and the higher the fluid pressure under the same injection rate [43,44].As shown in Figure 14, when the viscosity of the liquid is high the breakdown pressure is large.Conversely, as the gas has a lower viscosity, the breakdown pressure is also smaller.This fact is illustrated in Figure 14, which shows that the relationship between breakdown pressure and fluid viscosity follows a power-law function.In addition, the higher the fracturing fluid viscosity, the more complex the fracture geometry local to the wellbore.This is because a fluid with relatively higher viscosity has a relatively slower flow rate in the near-wellbore region [45].Figure 15 shows the variations of critical perforation angle with different viscosities of fracturing fluids.For the liquid-based hydraulic fracturing, the critical perforation increases with an increase in liquid viscosity.While for the gas fracturing, the critical perforation angle does not change a lot with the changes of gas viscosity.Although it is essential to create a dominant hydraulic fracture for hydrocarbons recovery, one should carefully design the fracturing fluid to avoid damage to the near-wellbore fracture system.

Conclusions
In this study, a fully-coupled geomechanical hydraulic fracture model is proposed under the rubric of FEM-based damage mechanics.This approach accommodates maximum tensile stress theory with shear failure according to Mohr-Coulomb theory and is applied to predict fracture propagation and complexity in heterogeneous rocks.The physical process involves fully coupling

Conclusions
In this study, a fully-coupled geomechanical hydraulic fracture model is proposed under the rubric of FEM-based damage mechanics.This approach accommodates maximum tensile stress theory with shear failure according to Mohr-Coulomb theory and is applied to predict fracture propagation and complexity in heterogeneous rocks.The physical process involves fully coupling fluid flow in the fracture and pressure diffusion into the surrounding porous rock while simultaneously accommodating rock damage, deformation, and fracture propagation.The form of the near-wellbore fractures initiating from two symmetrically-disposed radial perforations are analyzed.The effects of perforation angle, horizontal differential stress, fracturing fluid, and initial pore pressure on hydraulic fracturing are investigated.The following conclusions are derived from the numerical simulations.
A larger perforation angle induces a relatively larger curvature of the fracture during hydraulic fracture reorientation and requires a longer distance for the hydraulic fracture to rotate to the direction that is completely aligned with the direction of the maximum horizontal stress.Hydraulic fractures do not always initiate at the oriented perforations and the fractures induced during hydraulic fracturing are not always even and regular.Hydraulic fractures are initiated both around the wellbore and the oriented perforations when the perforation angle is >75 • .A larger perforation angle may cause a more complex fracture morphology.Therefore, perforations at a high angle with respect to the maximum principal stress are generally not recommended.
Stress ratio has a significant influence on fracture initiation pressure, breakdown pressure, and fracture propagation.Under the condition of hydrostatic pressure (λ = 1), the hydraulic fractures initiate and propagate along the direction of the oriented perforation, and thus will not form complex fracture morphology.The resulting fracture propagation does not always follow the directions provided by the perforations but rather reorients to the direction of the maximum horizontal stress.A larger stress ratio results in reduced fracture initiation and breakdown pressures.A larger stress ratio has a smaller critical perforation angle, thus, more easily to form complex fracture morphology.
The resulting differential initiation pressure and breakdown pressure decrease with an increase in the initial pore pressure.Particularly, the initiation pressures of the numerical simulations are in good agreement with the solutions of the Hubbert-Willis [40] and Haimson-Fairhurst [38] equations.In addition, since the adsorption-induced strain is fully coupled to the gas fracturing model, the behavior during gas fracturing is significantly different to that during water-based fracturing.The gas penetration not only alters the pore pressure in the reservoir but the gas adsorption-induced damage also modifies the mechanical properties of the reservoir rock.
The fracturing fluid influences the initiation pressure, breakdown pressure, and critical perforation angle by affecting the fluid flow behavior and the evolution of the poroelastic stress around the wellbore.The fracture initiation pressure and breakdown pressure increase with an increase in the viscosity.For the liquid-based hydraulic fracturing, breakdown pressure increases 55.6%, and the critical perforation angle increases from 70 • to 80 • , when the viscosity is increased to 1 Pa•s.While for the gas fracturing, breakdown pressure increases 26.8% when the viscosity is increased to 10 −4 Pa•s.Notably, the critical perforation angle, remaining 62 • to 63 • , and does not change a lot with the varies of gas viscosity.The critical perforation angle is sensitive to the viscosity of liquid and is insensitive to the viscosity of gas.
The proposed model can be used to determine an optimal perforation condition in various hydraulic fracturing.Research findings are applied to optimize perforation design in order to minimize high treatment pressures and to control complexity in the near-wellbore fracture region.This study is of great significance in further understanding the near-wellbore impacts on hydraulic fracture propagation and complexity.

Figure 2 .
Figure 2. Distributions of rock properties for different homogeneous indices (Mechanical parameter 0 u is 1).

Figure 2 .
Figure 2. Distributions of rock properties for different homogeneous indices (Mechanical parameter u 0 is 1).

Figure 3 .
Figure 3.The damage constitutive law for rock under uniaxial stress condition.

where 1 ε and 3 ε
are the major and minor principal strains (m/m), respectively.

Figure 3 .
Figure 3.The damage constitutive law for rock under uniaxial stress condition.

Figure 4 .
Figure 4. Simulation of uniaxial compression of a sample compared with experimental observations by Chen et al. [11].

Figure 4 .
Figure 4. Simulation of uniaxial compression of a sample compared with experimental observations by Chen et al. [11].

Figure 5 .
Figure 5.Comparison of model simulation results with experimental results.(a) Experimental fracture pattern from Chen et al. [11].(b) Fracture pattern obtained from numerical simulation.

Figure 6 .
Figure 6.Geometry of the numerical model used to investigate near-wellbore hydraulic fracturing processes.

Figure 5 .
Figure 5.Comparison of model simulation results with experimental results.(a) Experimental fracture pattern from Chen et al. [11].(b) Fracture pattern obtained from numerical simulation.

Figure 5 .
Figure 5.Comparison of model simulation results with experimental results.(a) Experimental fracture pattern from Chen et al. [11].(b) Fracture pattern obtained from numerical simulation.

Figure 6 .Figure 6 .
Figure 6.Geometry of the numerical model used to investigate near-wellbore hydraulic fracturing processes.

Figure 7 .
Figure 7. Numerical results of near-wellbore fracture patterns for different preexisting perforation angles.

Figure 8 .
Figure 8.Comparison of experimental results and numerical results for breakdown pressures.

Figure 9 .
Figure 9. Simulation results of fracture tortuosity/complexity of hydraulic fracturing produced at various stress ratios.

Figure 10 .
Figure 10.Initiation and breakdown pressures with different stress ratios.

Figure 9 . 20 Figure 9 .
Figure 9. Simulation results of fracture tortuosity/complexity of hydraulic fracturing produced at various stress ratios.

Figure 10 .
Figure 10.Initiation and breakdown pressures with different stress ratios.

Figure 10 .
Figure 10.Initiation and breakdown pressures with different stress ratios.

Figure 11 .
Figure 11.The relationship between critical perforation angle and stress ratio.

Figure 11 .
Figure 11.The relationship between critical perforation angle and stress ratio.

Figure 12 .
Figure 12.Relationship between initiation pressure and initial pore pressure.

Figure 12 .
Figure 12.Relationship between initiation pressure and initial pore pressure.

Figure 13 .
Figure 13.Breakdown pressure as a function of initial pore pressure.

Figure 14 .
Figure 14.Relationship between breakdown pressure and fracturing fluid viscosity.

Figure 15 .
Figure 15.The relationship between critical perforation angle and stress ratio.

Figure 15 .
Figure 15.The relationship between critical perforation angle and stress ratio.

Table 1 .
Comparison of experimental data and numerical data.

Table 1 .
Comparison of experimental data and numerical data.

Table 1 .
Comparison of experimental data and numerical data.

Table 2 .
Input parameters for the near-wellbore hydraulic fracturing model.

Table 2 .
Input parameters for the near-wellbore hydraulic fracturing model.

Table 3 .
Initiation pressure and breakdown pressure for different stress ratios.

Table 3 .
Initiation pressure and breakdown pressure for different stress ratios.

Table 3 .
Initiation pressure and breakdown pressure for different stress ratios.