Regular Article
Numerical simulation of horizontal well network fracturing in glutenite reservoir
^{1}
School of Petroleum Engineering, China University of Petroleum (East China), Qingdao, Shandong
266580, China
^{2}
Downhole Service Company, Chuanqing Drilling Engineering Company, Chengdu, Sichuan
610051, China
^{3}
Haimo America, Inc, Houston, Texas
77042, USA
^{*} Corresponding author: limingzhong_upc@hotmail.com
Received:
25
October
2017
Accepted:
7
August
2018
Hydraulic fracturing is an economically effective technology developing the glutenite reservoirs, which have far stronger heterogeneity than the conventional sandstone reservoir. According to the field production experience of Shengli Oilfield, horizontalwell fracturing is more likely to develop a complex fractured network, which improves the stimulated volume of reservoir effectively. But the clear mechanism of horizontalwell hydraulic fracture propagation in the glutenite reservoirs is still not obtained, thus it is difficult to effectively carry out the design of fracturing plan. Based on the characteristics of the glutenite reservoirs, a coupled FlowStressDamage (FSD) model of hydraulic fracture propagation is established. The numerical simulation of fracturing expansion in the horizontal well of the glutenite reservoir is conducted. It is shown that a square meshlike fracture network is developed near the horizontal well in the reservoir with lower stress difference, in which fracture is more prone to propagate along the direction of the minimum principal stress as well. High fracturing fluids injection displacement and high fracturing fluid viscosity lead to the rise of static pressure of the fracture, which results in the rise of fracture complexity, and greater probability to deflect when encountering gravels. As the perforation density increases, the microfractures generated at each perforation gather together faster, and the range of the stimulated reservoir is also relatively large. For reservoirs with high gravel content, the complexity of fracture network and the effect of fracture communication are obviously increased, and the range of fracture deflection is relatively large. In the case of the same gravel distribution, the higher the tensile strength of the gravel, the greater fracture tortuosity and diversion was observed. In this paper, a simulation method of horizontal well fracture network propagation in the reservoirs is introduced, and the result provides the theoretical support for fracture network morphology prediction and plan design of hydraulic fracturing in the glutenite reservoir.
© M. Qi et al., published by IFP Energies nouvelles, 2018
This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
1 Introduction
The glutenite reservoirs are generally characterized by lowporosity, low permeability and deep burial. The core analysis of glutenite reservoir shows that the pore space is intergranular and the radius of pore throat is much smaller than sandstone reservoir (Meng et al., 2010). With all these reservoir characteristics, hydraulic fracturing is inevitable in glutenite field development. Horizontal well fracturing process can effectively utilize the rock characteristics of glutenite reservoirs and the interference effect between microfractures generated from each perforation, forming a more complex fracture network, thus increasing the stimulated volume of the reservoir and improving the field oil production. According to the field application experience of the glutenite block in Shengli Oilfield, China, after the horizontal fracturing, the daily well production is 13.7t, much more than 1.2t with conventional straight well fracturing treatment, the utilization rate of reserves is increased to 96.7% as well (Zhao, 2014), which shows that the horizontal well fracturing technology has a good field application prospect in reconstruction of glutenite reservoirs.
In the past several years, a lot of effort has been put in research on the numerical simulation of hydraulic fracture propagation, and several models for predicting the fracture morphology, including PK, PKN, KGD, P3D, and PL3D, etc. have been put forward (Detournay, 1990). In the research on hydraulic fracture propagation, many coupling terms are considered, including the fracture plane deformation caused by fluids pressure, fluidsolid coupling in fluids flow inside the fracture, influence of reservoir heterogeneity, etc. (Guo et al., 2015b, 2016; Olson and Taleghani, 2009; Rui et al., 2018). Currently, numerical simulation of hydraulic fracturing for formation with heterogeneous media and random fracture propagation directions is the hotspot of research (Delorme et al., 2013a, 2013b; Li et al., 2011; Weng et al., 2011). The numerical simulation of complex fracture propagation mainly involves the following five methods: (1) conventional finite element method: fracture propagation is simulated by the principle of damage mechanics. One is cohesive element method (Paul et al., 2017), which does not apply to random fracture propagation. In another method, all elements are assumed to have the same type. According to the judgment criteria, the dynamic fracture propagation can be simulated by changing the material properties of the elements (Faivre et al., 2016). Although this method is affected by the size of elements with lower precision, it is easy to simulate the random fracture propagation; (2) extended finite element method (Belytschko and Black, 2015): this method has a high precision in simulating fracture propagation. When there are a large amount of other heterogeneous media (natural fractures), the results may be difficult to converge; (3) boundary element method (Olson, 2008): it has certain advantages in simulating random fracture propagation, but may have problems in treating the load and flow on fracture surface; (4) meshless method (Moës et al., 2015): this is a new approach to the simulation of fracture propagation, and easy to use for random fracture propagation. However, theoretical maturity and calculation efficiency should be further improved; (5) analytic model (Weng et al., 2011): the simulation precision is lowered due to the use of analytic model and formula for determining the direction of fracture propagation.
A lot of effort has been put in research on the hydraulic fracture propagation in the glutenite reservoirs (Zhao et al., 2007, 2015). Li et al. adopted the finite element method (FEM) to simulate the dynamic process of rock failure in the glutenite reservoirs, conducted research on influence of horizontal geostress difference, gravel size and content on the fracture morphology, and identified five fracturing modes (termination, branching, deflection, penetration, and attraction) of fracture encountering the gravel based on numerical simulation, obtaining some valuable conclusions (Li et al., 2013). Ma et al. conducted the physical experiment in six natural glutenite rock specimens, and quantified the influence of gravel size, geostress difference, and fracturing fluids type on the fracture morphology, width, and breakdown pressure (Ma et al., 2017). Meng et al. applied the largesize true triaxial hydraulic fracturing system to simulate the hydraulic fracture propagation in the artificial specimens, simulated the fracture propagation in the conditions of presence or nonpresence of gravel, different geostress difference, and different gravel sizes, and quantified the influence of gravel on fracture propagation by combining with fracturing treatment curves (Meng et al., 2010). Liu et al. applied the technology of high resolution microcomputed tomography in the glutenite specimen from the hydraulic fracturing true triaxial physical experiment, obtaining the visual description of fracture propagation in the glutenite specimen (Liu et al., 2004). Based on CDEM method, and referring to the result of the physical experiment by Liu et al, Ju et al. established the 3D model of the glutenite specimen, and conducted the influence of heterogeneity and principal stress difference on the position of fracture initiation and propagation (Ju et al., 2016). In the above study, the rule of fracture propagation in the glutenite reservoirs was simulated in terms of experiment and numerical simulation, but the scale and duration of simulation were limited, and there was little research on influence of gravel and matrix rock mechanical parameters, fracturing fluids displacement and viscosity, gravel size, etc. on the propagation of hydraulic fracture.
In this paper, a damage mechanics method is introduced to establish a 2D coupled FlowStressDamage (FSD) model of hydraulic fracture propagation in the glutenite reservoirs, and the model is resolved through programming. The influence of reservoir principal stress difference, hydraulic fracturing treatment parameters, gravel content, tensile strength, etc. on the hydraulic fracture propagation, is studied, aiming at capturing the profound understanding of horizontal well hydraulic fracture network propagation in the glutenite reservoirs.
2 Mathematical model of fracture propagation
Compared with other hydraulic fracturing numerical simulation models, the FSD coupled model can simulate the progressive failure process of materials. Its calculation method is based on the finite element and statistical damage theory, taking the heterogeneity of the material properties into consideration, and reflecting in the simulation result (Li et al., 2012; Yang et al., 2011). In the process of simulating fracturing, failure treatment was carried out for units that satisfy a given strength criterion, thus achieving the simulation of heterogeneous material fracturing. The FSD model can be applied not only to simulate the dynamic fracturing process of various reservoirs (including but not limited to shale, Guo et al., 2017; glutenite, hard rock, Zhu et al., 2015 and brittle rock, Li et al., 2011; etc.), but also to the study of fracturing in deviated (Zhu et al., 2015) and horizontal wells (Guo et al., 2015a).
Based on the classical Biot fluidsolid coupling equation, the damage mechanics method was used to establish the 2D fracture propagation model with seepagestressdamage coupling (Biot, 1941, 2004). By selfcompiled FORTRAN program, the finite element simulation of hydraulic fracturing was carried out. Through the redefinition of the mechanical parameters of the meshes in the finite element model, the effect of gravel in the sand conglomerate reservoir on the fracture expansion is simulated. The ABAQUS software was used to open the grid file and select grids that need to be set to gravel, and the independent parameter assignment software was used to assign the gravel grid separately. Finally, the other grids were assigned in the computation program. The software “TECPLOT” was used for data analysis and visual processing.
Rock deformation is modeled with the theory of linear elasticity. Coupling of fluid flow and geomechanical deformations are based on the poroelastic theory developed by Biot. The fluid flow inside the fracture is usually simplified to flow along a channel with a very narrow opening, which is represented by a nonlinear partial differential equation that relates the fluid flow velocity with the fracture width and pressure gradient along the fracture. The fracture propagation process is usually considered in the framework of Linear Elastic Fracture Mechanics (LEFM) theory (Tian et al., 2001). To simplify the processes, we used different type of element with their own properties in our model. Fluid element represents the induced hydraulic fracture. It has very high permeability and very low value of Young’s modulus and Poisson’s ratio. Gravel element represents the gravel inside the reservoir which has higher tensile strength and Young’s modulus than normal formation element and has lower permeability and porosity. We have used formation element with formation properties.
2.1 Seepage equation (Chen et al., 1995)
(1)where P is the pore pressure (MPa); μ is the fluid viscosity; C _{ t } is the total compression factor (C _{ f } + C _{ l }), C _{ f } is the rock compressibility, C _{ l } is the fluid compressibility; t is the time; K is the permeability (×10^{−3} μm^{2}). This equation is the mathematical equation for unstable seepage of elastic, porous, singlephase microcompressive fluid.
2.2 Deformation field equation
(1) Equilibrium equation (Yang et al., 2002)
Suppose the deformation is linear elastic. According to the effective stress principle, the equation characterizing the coupling between the stress field of solid and pressure field of fluid can be obtained:(2)where is the tensor of effective stress; P is the pore pressure (MPa); f _{ i } is the external load (MPa); δ _{ ij } is the Kronecker constant; α is the effective stress coefficient ranging between 0 and 1.
(2) Geometric equation(3)where ε _{ ij } is strain tensor; u is displacement.
(3) Constitutive equation (Rahman et al., 2009)
The stress and strain of elastic material are linearly related:(4)
The constitutive relation of damage mechanics is expressed as follows:(5)where D is the damage variable:(6)where E _{ c } is the elastic modulus of the damage material; E is the elastic modulus of the damage material.
(4) Heterogeneity of rock medium
In order to describe the heterogeneity of the reservoir rock, suppose its mechanical properties obey the Weibull distribution, which is defined by the following distribution density function (Weibull, 1951):(7)where σ _{ c } is the mechanical parameter for the elements of rock medium (elastic modulus, Poisson’s ratio). σ _{0} is the parameter related to the average of all element parameters; m is uniformity coefficient, which reflects the uniformity of rock medium.
(5) Permeabilitydamage coupling equation
When the maximum tensile principal stress of the element reaches the tensile strength of the rock (σ _{3} ≤ −σ _{ t }), tensile damage occurs (see criteria for tensile failure). The constitutive equation for rocks with elastic brittleness and residual stress is shown as follows:(8)
Permeability k is given by:(9)where σ _{t} is uniaxial tensile strength; σ _{tr} is residual strength at initial tensile damage; ε _{t0} is the threshold of strain under tensile damage according to criteria for uniaxial tension. When the strain under uniaxial tension reaches ε _{t0}, the element begins to undergo damage, but does not totally lose stress bearing ability. When the maximum tensile principle strain reaches the specific ultimate strain ε _{tu}, the element is considered to completely lose the stress bearing ability. The element is completely damaged and undergoes tensile failure, i.e. D = 1. At this time, in the permeability equation, the effective stress coefficient α is equal to 1. ξ (ξ > 1) is the coefficient of sudden jump of permeability. It characterizes the times by which the permeability is increased before and after damage of element under the same stress state. β is coupling parameter, which can be determined by stressstrain permeability test.
Under the action of compressive or shear stress, the MohrCoulomb yield criterion can be used as the second criterion of damage threshold (criterion of shear failure):(10)where σ _{1} and σ _{3} are maximum and minimum principal stress; φ is inner frictional angle; f _{ c } is uniaxial compressive strength.
The constitutive equation under uniaxial compressive stress is written as follows:(11)
Permeability k is expressed as(12)where f _{cr} is residual stress; ε _{c0} is the maximum compressive principal strain when the maximum compressive principal stress reaches the level of uniaxial compressive strength.
Combining equations (1)–(12) with initial conditions and boundary conditions, the overall control equation is obtained. The discrete control equations in finite element method (Eqs. (1) and (2)) are used to iteratively solve the two processes of fluid seepage and rock deformation.
3 Model validation
To validate the rationality of numerical model, a comparison of numerical simulation with the result of physical simulation of hydraulic fracturing in the natural glutenite specimens from Ma et al. (2017) is conducted, with same experimental parameters, and similar gravel sizes and distribution morphology. The accuracy of the model is validated by determining whether the fracture propagation morphology and propagation rule of the fracture tip encountering the gravel in the numerical simulation is the same as those in the physical simulation. According to the parameters listed in Table 1, the numerical simulation is operated with Specimen I with larger gravel grain and Specimen II with smaller gravel grain.
Experimental data in physical simulation and numerical model.
As shown in Figures 1a and 1b, in Specimen I with higher content of pebble gravel with sizes <20 mm, the hydraulic fracture morphology is dominated by the singlewing fracture along the horizontal maximum principal stress, with fracture propagation tendency along the interface between gravel and sandstone, but the gravel has little influence on the global fracture morphology, which is due to small gravel size. The result of numerical simulation in Figure 2 (Specimen I) shows that though the fracture deflects locally after encountering the gravel, the fracture generally propagates along the maximum horizontal geostress, forming the relatively straight fracture plane.
Fig. 1. Result of physical simulation carried by Ma et al. 
Fig. 2. Result of validating the software simulation. 
In Specimen II with higher content of gravel with sizes between 40 and 100 mm (Figs. 1c and 1d), the fluids are injected where the gravel is located, and the result of fracturing is shown in Figure 1d. Compared with Specimen I, Specimen II has more complex fracture propagation morphology, forming a series of fracturing morphologies such as fracture reorientation and branched fractures, which is affected by largesize gravel. The fracture propagation trajectory in the numerical simulation is shown in Figure 2. Specimen II displays more complex morphology than Specimen I, and after encountering the largesize gravel, the branched fractures are generated inside the gravel, forming complex fracture morphologies, such as deflection penetrating the gravel, and propagation around the gravel.
Comparison shows that the fracture propagation from the result of numerical simulation is basically consistent with that of physical simulation not only in terms of fracture global propagation morphology but also in concrete details, e.g., propagation of fracture encountering gravel (penetrating gravel, propagation around gravel, and attraction of gravel to the fracture) is displayed similarly in both simulations. Thus, the numerical model has higher reliability, and it is applicable in the simulation of hydraulic fracture propagation morphology in the glutenite reservoirs.
4 Results and analysis
The schematic of the simulating model is shown in Figure 3. Both the X and Y direction of the model is of the length of 2 m, where the X direction is the horizontal minimum stress direction, and the Y direction is the horizontal maximum stress direction. The basic model uses a perforation density of 16 perforations per meter, and to minimize the effect of the boundary effect, the perforations were placed symmetrically in the middle of the model. The gravel size is 40–100 mm, evenly distributed in the model. The basic model parameters are mainly derived from logging and core analysis of glutenite reservoirs in Yong’an block of Shengli Oilfield, as partly shown in Table 2. Other mechanical parameters of matrix and gravel used in the simulation are the same as those in the core experiments, as shown in Table 1.
Fig. 3. Schematic of the numerical model for horizontal well fracturing. 
Mechanical properties and injection parameter of the basic model.
4.1 Influence of horizontal geostress difference
Hydraulic fracture always propagates along the maximum horizontal geostress (σ _{ν} > σ _{H} > σ _{h}), thus the straight double wings more likely occur in the formation with a higher difference between maximum horizontal geostress and minimum horizontal geostress (Δσ = σ _{H} − σ _{h}) (Guo et al., 2015b). Under the influence of the gravels, the hydraulic fracture pattern shows greater complexity than other reservoirs especially for those with a lower geostress difference. For the glutenite reservoir with horizontal well fracturing, the stress state of the reservoir leads to different fracture network morphology, as shown in Figure 4.
Fig. 4. Influence of geostress difference on propagation trajectory of fracture network in the glutenite reservoir (σ _{ν} > σ _{H} > σ _{h}, Δσ = σ _{H} − σ _{h}). 
In Model 4(a) with a stress difference of 5 MPa, main fractures are generated at perforations that located at the edge of the section. In the middle of the perforation section, the microfractures generated from the tip of the perforations generally extend along the X direction and converge to each other. The morphology of the fracture network nearby the well is like square mesh, in this paper we call this kind of fracture network square meshlike network. The main fractures that extended from the perforations at edge also produced substantial steering under the influence of gravels and the interference effect of multifracture. In Model 4(b) with a stress difference of 10 MPa, although there are still some microfractures generated from the middle of perforation section propagate towards X direction and intersect with others, the outwardly extending fractures are not confined to those generated by the perforations at the edge. Main fractures that extended outwardly will still occur deflection under the influence of gravel and intersect with each other. The stress difference in Model 4(c) is 15MPa, which can be seen different from the other two models. Although square meshlike fracture network still exists in the middle of the perforation section, the fractures generated by the upper perforations on the left side of the model are converged into the main fracture with a bigger width and extend along the oblique 45degree angle, and the lower part form into the main fracture that extends straightly. The area covered by the fracture network in Model 4(c) is relatively smaller than the other two models.
According to the previous analysis, when the stress difference is small, fractures generated from the middle of perforation section will expand along the direction of the minimum principal stress under the influence of the multifracture disturbance and the influence of the gravel. When the geostress difference is large, the deflection effect of the main fracture is not that obvious. With the comparison of these models, the number of main fractures and the stimulated reservoir area in Model 3(b) with the stress difference of 10 MPa is the largest, so the application effect of horizontal well fracturing is better in reservoirs that have an appropriate geostress difference.
4.2 Influence of fracturing fluids injection displacement
The fracturing fluids injection displacement is a treatment parameter requiring deliberate concern during the hydraulic fracturing plan design in the glutenite reservoir. In the process of hydraulic fracturing, the fracturing fluids injection displacement influences the hydrostatic pressure inside the fracture and the penetration force on gravel by fracture tips reaching the gravel, and thus it influences the gravel failure, and propagation morphology and length of the fracture. Also, high displacement is conducive to the formation of more complex hydraulic fracturing network. The study of Olson has shown that the complexity of the fracture path is closely related to the high relative net pressure (low horizontal stress difference or high net pressure) (Olson 2008; Olson and Taleghani, 2009). The high net pressure is more related to the induced stress and the interaction between multiple fractures, conducting to the bending of the fracture expansion path, thus forming a complex fracture network.
The expansion of hydraulic fracturing network under three sets of injection displacement is shown in Figure 5. The injection displacement of Model 5(a) is 0.0005 m^{3} s^{−1}. The fracture width in this model is relatively short, and microfractures that extending from different perforations are easy to converge and communicate with each other. But the formation of the main fracture is difficult, and the stimulated area is significantly smaller than the other two models. The injection displacement of Model 4(b) is 0.0009 m^{3} s^{−1}, which leads to a much more complicated fracture network due to the rise of the net pressure. The confluence effect of the microfractures generated from the perforations is weakened concerning smaller displacement, and 2–3 microfractures are converged to form an independent main fracture easily, and six of these main fractures can be found in this model. The injection displacement in Model 5(c) is 0.0012 m^{3} s^{−1}. In Model 5(c), the fracture tip frequently deflects after encountering gravels, resulting in a higher tortuosity of the fracture than the other two sets of models. The fracture network morphology in the middle of the perforation no longer appears in a square meshlike shape, and the deviation of the main fractures is much more obvious, especially which generated from the leftmost perforation. But the number of main fractures in this model is five, one less than Model 5(b), thus affecting the model’s stimulated area. Therefore, with the increase in injection displacement, the fracture deflection is more obvious, but the reservoir stimulation effect is not necessarily better.
Fig. 5. Influence of fracturing fluids injection displacement on propagation trajectory of fracture network in the glutenite reservoir. 
4.3 Influence of fracturing fluids viscosity
The viscosity of the fracturing fluid and the injection displacement affect the relative static pressure of the tip of the hydraulic fracture, which affects the mutual interference between multifractures and the expansion after encountering gravels. Also, the high viscosity fracturing fluids benefit in forming higher and wider fracture and reduce the fracturing fluids leak off effectively. The viscosity of fracturing fluids changes obviously as the formation temperature and pressure vary after the injection into the reservoir. Therefore, it is of significance to study the influence of fracturing fluids viscosity on the fracture propagation (Fig. 6).
Fig. 6. Influence of fracturing fluids viscosity on propagation trajectory of fracture network in the glutenite reservoir. 
In Model 6(a), the viscosity of the fracturing fluid is 20 MPa s. Compared with the other two models, the fractures are more prone to propagate around the gravel instead of penetrating them. Due to low fracture net pressure, no main fracture was formed in the middle of the perforation section, and the fracture deviation is not that significant as the other two models. The fracture network morphology of Model 6(b) and Model 6(c) are roughly similar. Fractures of Model 6(c) with a higher fluid viscosity have better connectivity and deflection degree. But the main fracture generated from the left third perforation converged with other fractures faster than that in Model 6(b), thus affecting the stimulated area. Therefore, fracturing fluid viscosity should not be too low when conducting horizontal well fracturing. But high fracturing fluid viscosity will accelerate the gathering of fractures, which is not conducive to develop a high stimulated area.
4.4 Influence of perforation density
When the distance between the perforations is close, the stress concentration and interference effect caused by other perforations and microfractures make the stress distribution more complicated, thus influencing the extension of the fracture network. The interaction effect of perforations is conducive to communication and gathering of the microfractures, which makes it much easier to generate main fractures. To study the influence of perforation density on the fracture network in the glutenite reservoir, three models with perforation density of 8, 16 and 24 perforations per meter are established. The simulation results are shown in Figure 7.
Fig. 7. Influence of perforation density on propagation trajectory of fracture network in the glutenite reservoir. 
In Model 7(a) with eight perforations, the induced stress and interference effects generated by perforations are relatively small, which result in the independence of each microfracture. Unlike the other two models, microfractures expand radically along the direction of horizontal maximum principal stress (Y direction) at first. While even if communication and gathering occur under the influence of gravels and the interference effect from other fractures during the process of fracture expansion, still the complexity and deviation effect of the fracture network is not compared with the other two models. As the perforation density increases, the interaction between perforations and microfractures gradually increases, which leads to a faster communication and aggregation rate. Microfractures form a square meshlike network as shown in the Model 7(b). Model 7(c) reflects the expansion of the fracture network under high perforation density. The speed of communication and gathering of microfractures is much faster, and a square meshlike network was formed immediately after the injection. The net pressure of the main fracture in Model 7(c) is greater than the others due to the aggregation of more microfractures, therefore result in a more obvious deflection.
4.5 Influence of gravel content
Compared with the conventional sandstone reservoir, the glutenite reservoirs have more complex lithology. For example, rock specimens sampling from the glutenite reservoirs in Yong’an Block, Shengli Oilfield show great heterogeneity, and the gravel sizes in the specimen vary between 1 cm and 20 cm, and with gravel content maximum of 64% and the minimum of 2%. A fracturing model of a horizontal well in sandstone reservoirs was established as a control group, and two models of gravel content 10% and 40% were set up to study the effect of gravel content on the formation of fracture network as shown in Figure 8.
Fig. 8. Influence of gravel content on propagation trajectory of fracture network in the glutenite reservoir. 
Model 8(a) reflects the expansion of the hydraulic fracture network under the influence of interaction between multiple fractures during the horizontal fracturing. Not that all the hydraulic fractures developed at each perforation can expand into the main fracture, but the action of deflection frequently occurs under the effect of interference, finally gathering with other fractures. With the absence of gravel, the main fractures in the model are mostly developed under the direction of maximum principal stress, but deflection caused by the border effect also exists in some fractures. Fracture network shows much more deflection and tortuosity in glutenite reservoir compared with the sandstone group. From the simulation results of Models 8(b) and 8(c), it can be seen that the effect of gravel content on the expansion of the fracture network is mainly reflected in the following tips: the communication of microfractures in Model 8(b) is relatively poor; for two symmetrical perforations, fracture expansion and communication from one side hinder that from the other, therefore no square meshlike network was formed; Model 8(c) shows a better fracture network complexity and deflection, which lead to a larger stimulated area.
4.6 Influence of gravel tensile strength
During the process of hydraulic fracturing, the most important rock failure form is the tensile failure. The fracture generation and propagation during hydraulic fracturing are dominated by the tensile failure in the rock matrix, thus the tensile strength of gravel is an important element influencing the propagation of hydraulic fracture encountering the gravel. A simulation of propagation trajectory of hydraulic fracture network in the glutenite reservoir is operated with the tensile strength of gravel (σ _{tg}) 15 MPa, 20 MPa, and 25 MPa.
It is shown in Figure 9 that the influence of the gravel tensile strength on the fracture network is mainly reflected in fracture deflection. The fracture is more prone to deflect in the model with higher strength gravel. The tensile strength of Model 9(b) is 20 MPa, the overall fracture morphology is similar to that in Model 9(a), but the fracture is more prone to deflect when encountering gravels, which results in a fracture network with higher tortuosity. Model 9(c) with gravel tensile of 25 MPa generates a fracture network that propagates toward the direction of X integral. But it is worth noting that the general fracture length, number of main fractures and stimulation area in Model 9(c) are the least of these three models. Therefore, the effect of horizontal well fracturing is better in the reservoir with lower gravel tensile strength.
Fig. 9. Influence of gravel tensile strength on propagation trajectory of fracture network in the glutenite reservoir. 
5 Conclusion
In this paper, the influencing factors of fracturing development of horizontal well in the glutenite reservoir are studied by numerical simulation. The hydraulic fracturing performance was evaluated based on the complexity of the fracture network, fracture deviation, and the stimulated reservoir area, and the following conclusions were obtained.

Based on the damage mechanics method, a coupled FSD model of hydraulic fracture propagation in the glutenite reservoirs is established, and the accuracy of the model is validated through comparison between results of the physical experiment and numerical simulation.

Square meshlike fracture network was formed in the middle of the perforation section at a geostress difference of 5 MPa. Fracturing deviation effect is not that obvious with a high geostress difference, and the stimulated area is also relatively lower.

The injection displacement raises the relative static pressure of the tip of the hydraulic fracture directly, as well as the viscosity of the fracturing fluid, which affects the complexity of the fracture network. Microfractures generated from the perforations are more easily to converge into main fractures with a higher displacement, and fracturing deflection after penetrating gravel occurs more often. Although fracturing fluid viscosity should not be too low when conducting horizontal well fracturing, high fracturing fluid viscosity will accelerate the gathering of fractures, which is also not conducive to develop the high stimulated area.

The increase of perforation density mainly enforces the interaction between the perforations and the microfractures, which can effectively promote the communication and connection. In the case of 24 perforations per meter, the aggregation rate of microfractures is obviously faster, and the net pressure of the main fracture is greater than the others, therefore result in a more obvious deflection.

The main fractures in the model are mostly developed along the direction of maximum principal stress with the absence of gravel, but deflection caused by interference and border effect also exists. With the gravel content increases, the communication speed of microfractures and complexity of the network increase, as well as the fracture network diversion. The fracture shows great deflection and tortuosity in the reservoir with high gravel tensile strength and is more prone to deflect when encountering gravels, but the stimulated area is smaller instead.
Acknowledgments
The authors would like to acknowledge the financial support of the National Natural Science Foundation of China (Grant No. 51404288), the Fundamental Research Funds for the Central Universities (Grant No. 17CX02077) and the Applied Basic Research Project for Qingdao (Grant No. 171120jch).
References
 Belytschko T., Black T. (2015) Elastic crack growth in finite elements with minimal remeshing, Int. J. Numer. Methods Eng. 45, 5, 601–620. [CrossRef] [Google Scholar]
 Biot M.A. (1941) General theory of threedimensional consolidation, J. Appl. Phys. 12, 2, 155–164. [CrossRef] [Google Scholar]
 Biot M.A. (2004) Theory of elasticity and consolidation for a porous anisotropic solid, J. Appl. Phys. 26, 2, 182–185. [CrossRef] [MathSciNet] [Google Scholar]
 Chen H.Y., Teufel L.W., Lee R.L. (1995) Coupled Fluid Flow and Geomechanics in Reservoir Study – I. Theory and Governing Equations, IEEE Systems Conference, Society of Petroleum Engineers. [Google Scholar]
 Delorme M., Daniel J.M., KadaKloucha C., Khvoenkova N., Schueller S., Souque C. (2013a) Reservoir Stimulation and Induced Microseismic Events on 3D Discrete Fracture Network for Unconventional Reservoirs, Unconventional Resources Technology Conference. [Google Scholar]
 Delorme M., Mota R.O., Khvoenkova N., Fourno A., Nœtinger B. (2013b) A methodology to characterize fractured reservoirs constrained by statistical geological analysis and production: a real field case study, Geol. Soc. London Spec. Publ. 374, SP37414. [Google Scholar]
 Detournay E. (1990) A Poroelastic PKN Hydraulic Fracture Model Based on an Explicit Moving Mesh Algorithm, J. Energy Resour. Technol 112, 4, 224–230. [CrossRef] [Google Scholar]
 Faivre M., Paul B., Golfier F., Giot R., Massin P., Colombo D. (2016) 2D coupled HMXFEM modeling with cohesive zone model and applications to fluiddriven fracture network, Eng. Fract. Mech. 159, 115–143. [CrossRef] [Google Scholar]
 Guo T., Zhang S., Zou Y., Xiao B. (2015a) Numerical simulation of hydraulic fracture propagation in shale gas reservoir, J. Nat. Gas Sci. Eng. 26, 847–856. [CrossRef] [Google Scholar]
 Guo T., Zhang S., Ge H., Wang X., Lei X., Xiao B. (2015b) A new method for evaluation of fracture network formation capacity of rock, Fuel 140, 778–787. [CrossRef] [Google Scholar]
 Guo T., Qu Z., Gong D., Lei X., Liu M. (2016) Numerical simulation of directional propagation of hydraulic fracture guided by vertical multiradial boreholes, Journal of Natural Gas Science & Engineering 35, 175–188. [CrossRef] [Google Scholar]
 Guo T., Li Y., Ding Y., Qu Z., Gai N., Rui Z. (2017) Evaluation of acid fracturing treatments in shale formation, Energy Fuels 31, 10, 10479–10489. [CrossRef] [Google Scholar]
 Ju Y., Liu P., Chen J., Yang Y., Ranjith P.G. (2016) CDEMbased analysis of the 3D initiation and propagation of hydrofracturing cracks in heterogeneous glutenites, J. Nat. Gas Sci. Eng. 35, 614–623. [CrossRef] [Google Scholar]
 Li G., Tang C.A., Li L.C. (2011) Threedimensional micro flowstressdamage (FSD) model and application in hydraulic fracturing in brittle and heterogeneous rocks, Key Eng. Mater. 452–453, 8, 581–584. [Google Scholar]
 Li L., Meng Q., Wang S., Li G., Tang C. (2013) A numerical investigation of the hydraulic fracturing behaviour of conglomerate in glutenite formation, Acta Geotech. 8, 6, 597–618. [CrossRef] [Google Scholar]
 Li L.C., Tang C.A., Li G. (2012) Numerical simulation of 3D hydraulic fracturing based on an improved flowstressdamage model and a parallel FEM technique, Rock Mech. Rock Eng. 45, 5, 801–818. [Google Scholar]
 Liu H.Y., Kou S.Q., Lindqvist P.A. (2004) Numerical simulation of the fracture process in cutting heterogeneous brittle material, Int. J. Rock Mech. Min. Sci. 41, 3, 14–19. [CrossRef] [Google Scholar]
 Ma X., Zou Y., Li N., Chen M., Zhang Y., Liu Z. (2017) Experimental study on the mechanism of hydraulic fracture growth in a glutenite reservoir, J. Struct. Geol. 97, 37–47. [CrossRef] [Google Scholar]
 Meng Q.M., Zhang S.C., Guo X.M., Chen X.H., Zhang Y. (2010) A primary investigation on propagation mechanism for hydraulic fractures in glutenite formation, J. Oil Gas Technol. 32, 4, 119–123. [Google Scholar]
 Moës Nicolas, Dolbow J., Belytschko T. (2015) A finite element method for crack growth without remeshing, Int. J. Numer. Methods Eng. 46, 1, 131–150. [Google Scholar]
 Olson J.E. (2008) Multifracture propagation modeling: Applications to hydraulic fracturing in shales and tight gas sands, American Rock Mechanics Association. [Google Scholar]
 Olson J.E., Taleghani A.D. (2009) Modeling simultaneous growth of multiple hydraulic fractures and their interaction with natural fractures. Society of Petroleum Engineers, DOI: 10.2118/119739MS. [Google Scholar]
 Paul B., Ndeffo M., Massin P., Moës N. (2017) An integration technique for 3D curved cracks and branched discontinuities within the extended finite element method, Finite Elem. Anal. Des. 123, 19–50. [CrossRef] [Google Scholar]
 Rahman M.M., Aghighi M.A., Rahman S.S., Ravoof S.A. (2009) Interaction between induced hydraulic fracture and preexisting natural fracture in a poroelastic environment: Effect of pore pressure change and the orientation of natural fractures, Society of Petroleum Engineers. [Google Scholar]
 Rui Z., Wang X., Zhang Z., Lu J., Chen G., Zhou X., Patil S. (2018) A realistic and integrated model for evaluating oil sands development with steam assisted gravity drainage technology in Canada, Appl. Energy 213, 76–91. [CrossRef] [Google Scholar]
 Tian Y., Tan C., Wan Z., Qi F. (2001) Coupling analysis of seepage and stresses in rock failure process, Chin. J. Geotech. Eng. 23, 4, 489–493. [Google Scholar]
 Weibull W. (1951) A statistical distribution function of wide applicability, J. Appl. Mech. 13, 2, 293–297. [Google Scholar]
 Weng X., Kresse O., Cohen C.E., Wu R., Gu H. (2011) Modeling of Hydraulic Fracture Network Propagation in a Naturally Fractured Formation, Society of Petroleum Engineers. DOI: 10.2118/140253MS. [Google Scholar]
 Yang T.H., Tham L.G., Tang C.A., Leng X.F., Li L.C. (2002) Influence of heterogeneity on hydraulic fracturing in rocks, Chin. J. Geotech. Eng. 24, 6, 724–728. [Google Scholar]
 Yang T.H., Xu T., Liu H.Y., Tang C.A., Shi B.M., Yu Q.X. (2011) Stressdamageflow coupling model and its application to pressure relief coal bed methane in deep coal seam, Int. J. Coal Geol. 86, 4, 357–366. [CrossRef] [Google Scholar]
 Zhao C. (2014) 3D fracturing network optimization techniques for horizontal wells in sandstoneconglomerate formations, Petrol. Drill. Tech. 5, 95–99. [Google Scholar]
 Zhao Y.Z., Qu L.Z., Wang X.Z., Cheng Y.F., Shen H.C. (2007) Simulation experiment on prolongation law of hydraulic fracture for different lithologic formations, J. China Univ. Petrol. 31, 3, 63–66. [Google Scholar]
 Zhao Z., Guo J., Ma S. (2015) The experimental investigation of hydraulic fracture propagation characteristics in glutenite formation, Adv. Mater. Sci. Eng. 2, 1–5. [CrossRef] [Google Scholar]
 Zhu H., Zhao X., Guo J., Jin X., An F., Wang Y., Lai X. (2015) Coupled flowstressdamage simulation of deviatedwellbore fracturing in hardrock, J. Nat. Gas Sci. Eng. 26, 711–724. [CrossRef] [Google Scholar]
All Tables
All Figures
Fig. 1. Result of physical simulation carried by Ma et al. 

In the text 
Fig. 2. Result of validating the software simulation. 

In the text 
Fig. 3. Schematic of the numerical model for horizontal well fracturing. 

In the text 
Fig. 4. Influence of geostress difference on propagation trajectory of fracture network in the glutenite reservoir (σ _{ν} > σ _{H} > σ _{h}, Δσ = σ _{H} − σ _{h}). 

In the text 
Fig. 5. Influence of fracturing fluids injection displacement on propagation trajectory of fracture network in the glutenite reservoir. 

In the text 
Fig. 6. Influence of fracturing fluids viscosity on propagation trajectory of fracture network in the glutenite reservoir. 

In the text 
Fig. 7. Influence of perforation density on propagation trajectory of fracture network in the glutenite reservoir. 

In the text 
Fig. 8. Influence of gravel content on propagation trajectory of fracture network in the glutenite reservoir. 

In the text 
Fig. 9. Influence of gravel tensile strength on propagation trajectory of fracture network in the glutenite reservoir. 

In the text 