- Email: [email protected]
A simulation method based on energy criterion for network fracturing in shale gas reservoirs
Journal of Natural Gas Science and Engineering 52 (2018) 295–303
Contents lists available at ScienceDirect
Journal of Natural Gas Science and Engineering journal homepage: www.elsevier.com/locate/jngse
A simulation method based on energy criterion for network fracturing in shale gas reservoirs
T
Haifeng Zhaoa,∗, Xiaohua Wanga, Wei Wangb, Erfei Mua a b
School of Petroleum Engineering, China University of Petroleum, Beijing, 102249, China Hancheng Branch, PetroChina Coalbed Methane Co., Ltd., Hancheng, Shanxi, 715409, China
A R T I C L E I N F O
A B S T R A C T
Keywords: Shale gas Network fracturing Energy criterion Stimulated reservoir volume Fluid-solid coupling
The fracture network in shale gas reservoirs has a plurality of extending frontiers and the propagation paths of fractures are difficult to predict in advance. When the network fracturing is simulated by the extended finite element method (XFEM) based on stress intensity factors, the computation is arduous and time-consuming, which limits the application of XFEM in field fracturing design. Therefore, a simulating method of crack propagation based on energy criterion is proposed. Based on the principle of energy balance, the energy criterion for controlling crack frontiers extending is established by the fluid injection energy, the rock strain energy, the surface energy, the seismic energy and the friction loss energy. The criterion only needs to calculate the energy distribution in the crack frontiers to judge whether the fracture could expand, avoiding the fussy calculation of the stress intensity factor. Based on the proposed energy criterion, the fluid solid coupling effect and the fractures interaction, an energy method is established for network fracturing simulation in shale gas horizontal well. The results show that the fracture network morphology obtained by simulation is basically consistent with micro-seismic monitoring result, and the computational efficiency of this simulation method is 10000 times higher than that of XFEM. The energy method provides an idea and basis for the simulation and design of shale reservoir fracturing.
1. Introduction The horizontal well with segmented volume fracturing is the main technology to exploit the shale gas effectively. The complex network is formed during the hydraulic fracturing process for shale reservoirs via the micro-seismic monitoring. Thus, the investigation of the formation mechanism and simulation of fracture network is essential (Wu et al., 2012; Zou et al., 2015; Fisher et al., 2004, 2005) and the criterion of rock failure and fracture extension is the key step. The calculation complexity and time consumption of the fracture network models are significantly different when different criterions of fracture propagation are used in numerical simulations. Therefore, the criterion study of rock failure and fracture propagation for shale reservoir fracturing plays a vital role in characterizing the formation of the fracture network. Most of the current fracture propagation models are on the basis of stress intensity factor (SIF) (Irwin, 1957). Some complex fracture network models such as the wire-mesh model (Xu et al., 2009), the unconventional fracture model (Weng et al., 2011) and the discrete network model (Meyer, 1989, Meyer and Bazan, 2011) are established. These models are solved by different numerical algorithms, such as
∗
boundary element (Olson and Taleghani, 2009; Hou et al., 2015), extended finite element (Taleghani, 2010; Zeng et al., 2015) and discrete element (Meyer, 1989; Meyer and Bazan, 2011; Wang et al., 2015). The fracture criterions used in the mentioned models include the maximum energy release rate criterion, the minimum strain energy density criterion and the maximum tensile stress criterion (Zhao, 2004), and the theoretical foundation of these criterions is the SIF calculation based on the fracture toughness criterion. However, the fracture network has many propagation frontiers and the fracture geometry is very complex, making it difficult to calculate the stress distribution accurately. And the establishment of the integral path to calculate SIF is sometimes impossible (Martin, 2000). In addition, the calculation of network simulation is very time-wasting based on the fracture toughness criterion, which deviates from the engineering demand. Therefore, we attempt to study hydraulic crack propagation from the viewpoint of energy balance. Salamon (1984) established the energy balance law for rock failure processes during the mining course and developed methodologies for calculating energy terms in rock mechanics. Yao et al. (2018) assumed that the fractures extend forward when hydraulic energy is greater than the sum of strain energy and
Corresponding author. E-mail address: [email protected] (H. Zhao).
https://doi.org/10.1016/j.jngse.2018.02.003 Received 31 May 2017; Received in revised form 1 February 2018; Accepted 4 February 2018 Available online 07 February 2018 1875-5100/ © 2018 Elsevier B.V. All rights reserved.
Journal of Natural Gas Science and Engineering 52 (2018) 295–303
H. Zhao et al.
surface energy, but the friction dissipation energy was not taken into account when using the Griffith stability criterion. Khademian (2016); Khademian et al. (2017) focused on the radiated seismic energy during fracturing. They calculated strain energy, friction work, dissipative plastic work, and radiated seismic energy for analyzing rock failures in compression or shear. Forasmuch as, the energy criterion controlling the extension of fracture frontiers is established by balancing hydraulic fracturing injection energy, strain energy, surface energy, seismic energy and friction loss energy based on the above theories. Then an energy method for fracturing simulation of shale gas horizontal well is proposed combined with the energy criterion and the fluid-solid coupling method used to calculate fracturing fluid flow rate and pressure distribution in the fracture network. The energy method avoids the complex calculation of SIF in the crack frontiers of fracture network and provides an idea and basis for the fracturing simulation and design in shale reservoirs.
n
⎛
N
⎞
n
⎠
i=1
∑ Ek = (q − qlti) ⎜σmin + St + ∑ ∑ Δpij ⎟ ∑ ti i=1 j=0
⎝
(1)
The whole surface energy of fracture network is n
U=
N
∑ ∑ γs hij lij
(2)
i=1 j=0
The strain energy of rock is n
N
∑ Es = ∑ ∑ σmin hij lij wij (3)
i=1 j=0
The frictional loss energy is n
∑ Ef
= (q − qlti)
i=1
2. The energy mechanism and propagation criterion of fracture network formation
n
N
∑ ti ∑ ∑ Δpij i=1 j=0
(4)
where i represents the i-th time step in the iterative calculation; j ren presents the j-th level fracture; q is the pump rate, m3/s; T = ∑i = 1 ti is the injection time of fracturing fluid, s; σmin is the minimum principal stress, MPa; St is the tensile strength, MPa; Δpij is the pressure drop in the fracture, MPa; γs is the rock surface energy density, J/m2; hij is the fracture height, m; lij is the fracture length, m; wij is the average fracture width, wij = (wi − 1, j + wi, j )/2 , m,. According to previous studies, the reservoir rock is in the state of drastic change of local ground stress field. The I-II mixed fracture usually occurs during the process of hydraulic fracturing. For the twodimensional plane, the rock surface energy density (Nuismer, 1974; Aimene and Nairn, 2014) of the I-II type mixed fractures can be expressed as Eq. (5).
2.1. Energy mechanism of fracture network formation The hydraulic fracturing is a dynamic coupling process of rock deformation and fracturing fluid flow. On one hand, the micro-cracks occur inside the rock under the pressure of the fracturing fluid, and gradually propagate through the rock. Ultimately, the macro-cracks are formed. On the other hand, the fracturing fluid flows in the cracks and supports a certain width of cracks. With the high displacement injection of the fracturing fluid, the fluid pressure increases instantaneously, forcing the cracks to expand continuously. The hydraulic fracture would encounter natural fractures during the expansion process of fracture. When the hydraulic energy supplied by the fracturing fluid is greater than the energy required by natural crack propagation at both crack tips, the crack branches would be formed. Then the branch cracks continue to propagate. The branch cracks will form branches once again if the fractures run into the other natural fracture. In this way, the fracture network consisting of hydraulic cracks and natural fractures is eventually formed. Hydraulic fractures not only form branches with natural fractures, but also produce bifurcation, which finally produces the crack system like a root. In the shale reservoir, the complex fracture network shown in Fig. 1 is composed of natural fractures, hydraulic fractures and induced joints.
γs =
2 2 + KIIC (KIC ) 2E′
(5)
In the case of plane stress: E′ = E ; plane strain: E′ = E /(1 − υ2) ; KIC and KIIC are I and II rock fracture toughness respectively, which could be measured by laboratory experiments.
3. Dynamic mechanism and numerical solution of fracture network propagation 3.1. Dynamic mechanism of fracture network propagation
2.2. Energy criterion for crack propagation Dynamic mechanism of fracture network expansion was studied from the perspective of fracture dynamic mechanics (Zhao et al., 2012). When hydraulic fracture encounters the natural fractures, it is considered that the formation of complex fracture network needs to simultaneously meet the extending conditions at both ends of the natural fractures. The pressure required for the simultaneous propagation of both ends of the natural fracture is higher than the pressure required by the quasi-static expansion of the previous fracture. The fracture dynamics is used to determine the critical pump rate which could achieve the simultaneous expansion of branch cracks. The mechanism of the network expansion is studied from the perspective of energy in this paper. When hydraulic fracturing is carried out with large and variable pump rate, there is obvious dynamic shock wave in the stratum. A large number of micro seismic events is captured by the micro seismic monitoring technology. In the meantime, the stress state of rock and the tip position of crack vary with time rapidly. The inertial terms of the system balance equation can't be neglected when the fluid and rock are taken as a whole. Therefore, the crack propagation should be treated as a dynamical problem in this case.
During the hydraulic fracturing process, the energy in the cracks mainly contains the following parts (Fig. 2): the hydraulic energy of fracturing fluid (Ek); surface energy (U) to form new cracks; the rock strain energy (Es) which is equivalent to the work done by the fluid pressure; hydraulic friction loss energy (Ef) including the internal friction loss of fracturing fluid and friction loss between fracturing fluid and cracks wall; seismic energy (Ew) that monitored by micro-seismic during the fracturing process. Since micro-seismic energy is very small compared to the total hydraulic energy, it is ignored in our study. On the basis of the formation of the complex fracture network after fracturing in shale reservoir, without considering the bifurcation of the fracture itself, an energy criterion for the expansion of the whole fracture network under the condition of fracture classification is given. When ∑ Ek > ∑ (U + Es + Ef ) , the hydraulic fractures begin to propagate; when ∑ Ek < ∑ (U + Es + Ef ) , the hydraulic fractures cease to expand. The physical quantities in the energy criterion are calculated by Eqs. (1)–(4). Total injected hydraulic energy is 296
Journal of Natural Gas Science and Engineering 52 (2018) 295–303
H. Zhao et al.
Fig. 1. The propagation process of the fracture network.
3.2. The mathematical model of fracture network propagation
ρ
3.2.1. The governing equation
D→ υij Dt
pij + η (S ) ∇2→ υij + 2D⋅∇η = −∇→
(6)
where D is the strain rate tensor, D = (∇υ + ∇υT )/2 ; ρ is fluid density, kg/m3; η is the fluid viscosity, Pa·s; → υij is the fluid velocity vector; S = 2tr (D 2) . In the view of the actual situation of crack propagation, Eq. (6) is simplified. The pressure drop in the direction of fracture length and fracture height is represented by Eq. (7).
1) The fluid migration equation in the fracture The flow of fracturing fluid in a single fracture is simplified as a flow between two smooth parallel plates. The flow between the parallel plates follows the Navier-Stokes equation. It is assumed that the fracture fluid is incompressible liquid and the effect of mass force on the flow field is neglected. The general equation of the whole fracture network system is represented by Eq. (6) (Yang et al., 2002; Liu et al., 2003).
∂pij ∂x
=
∂ (υx )ij ⎤ ∂pij ∂ (υz )ij ⎤ ∂ ⎡ ∂ ⎡ , η (S ) η (S ) = ⎥ ⎢ y z y ∂y ⎢ ∂ ∂ ∂ ∂y ⎥ ⎣ ⎦ ⎣ ⎦
(7)
For power-law fluids, its viscosity follows the relation given by Eq.
Fig. 2. Diagram of energy distribution during the fracturing.
297
Journal of Natural Gas Science and Engineering 52 (2018) 295–303
H. Zhao et al.
For two-dimensional flow, the viscosity function η (S ) can be simplified:
Table 1 The specific parameters of the example. Parameters
Value
Parameters
Value
Overburden stress Maximum principal stress Minimum principal stress
63.5 MPa 58.5 MPa
Fracture toughness I Fracture toughness II
0.8 MPa m0.5 1.14 MPa m0.5
52.5 MPa
52.98J/m2
Young's modulus of rock Poisson's ratio Tensile strength
34.5 GPa 0.24 5 MPa
Rock strain energy density Fracturing fluid viscosity Average pump rate
⎜
⎟
⎜
⎟
(9)
The vector equation (6) is expanded and simplified. The volume flow rate of the fracturing fluid along the X and Z directions are expressed by fluid pressure pij and crack width wij (Kozicki et al., 1966; Nolte, 1979):
10 mPa s 13 m3/min
⎧ qx = − S
∂pij ∂x
⎨ q = −S ∂pij ∂z ⎩ z
(8).
η = K ′⋅S
2 (n ′− 1)/2
2
∂ (υz )ij ⎞ ⎤ ⎡ ∂ (υx )ij ⎞ η = K′⎢⎛ +⎛ ⎥ y ∂ ⎠ ⎝ ∂y ⎠ ⎦ ⎣⎝
(n ′− 1)/2
(8)
where S =
1 2n ′ K ′− n ′ 2n ′ + 1
Fig. 3. Fracture network shape at different time.
298
(10) wij
( ) 2
2n ′+ 1 n′ ⎡
⎢ ⎣
2
n ′− 1 2 − 2n ′
( ) + ( ) ⎤⎥⎦ ∂pij ∂x
∂pij ∂z
.
Journal of Natural Gas Science and Engineering 52 (2018) 295–303
H. Zhao et al.
Fig. 4. The energy variation over time during the fracturing.
∂Aij (x , t ) ∂ ⎛ ∂pij ⎞ ∂ ⎛ ∂pij ⎞ + ⎜S ⎟ + ⎜S ⎟ = ∂t ∂z ⎝ ∂z ⎠ ∂x ⎝ ∂x ⎠
Table 2 The geometric results of a single segment after fracturing. Geometric parameters
simulation results Micro-seismic interpretation results
Fracturenetwork length (m)
Fracturenetwork width (m)
365.34 334.53
92 101.41
Fracture height (m)
75 75.54
SRV/ (106m3)
1.32 1.53
Maximum width at main seam (mm)
∂ ∂x
5.4 -
=
where qij ∂ (q )ij
is
the
(12)
1 w (2n′+ 1)/ n′ 2n ′ K ′− n′ (n′+ 1)/ n′ ⎡ 2 ⎢ ⎨ 2n ′ + 1
⎣
⎧
2n ′
⎨ 2n ′ + 1 ⎩ ∂Aij (x , t ) ∂t
1 w
K ′− n′ +
(2n′+ 1)/ n′
2(n′+ 1)/ n′
⎡ ⎢ ⎣
2
n′− 1 2 − 2n′ ∂p ⎫ ij ⎤ + ∂x ⎬ ⎥
( ) ( )⎦ ∂pij
+
∂x
2
∂pij ∂z
( ) + ( ) ⎤⎦⎥ ∂pij ∂x
∂pij
⎭
n′− 1 2 − 2n′ ∂p ⎫ ij
∂z
∂z
⎬ ⎭
2cl hij
(13)
t − τ (x , z )
Therefore, the relationship between fluid pressure pij and width wij at different positions in the fracture network is established.
The continuity equation of the fracturing fluid inside the fracture is (Amini et al., 2015; Carter, 1957):
∂t
⎧ ⎩
∂ ∂z
∂Aij (x , t )
t − τ (x , z )
The above equations are expanded:
2) The continuity equation
∇⋅qij + (ql )ij +
2cl hij
3) The governing equation of rock deformation
=0
injection
(11) volume
rate
of
fracturing
Assuming that the rock deformation is elastic, the displacement governing equation of reservoir rock is expressed as (Chun and Ghassemi, 2012; Yew and Weng, 2015):
fluid,
∂ (q )ij
∇⋅qij = ∂xx + ∂zz ; qx and qz are volume flow rate in the X and Z directions respectively; (ql )ij is the filtration loss of fracturing fluid, (ql )ij = 2cl hij / t − τ (x , z ) ; cl is the fluid leakoff coefficient; τ (x , z ) is the time which the leak begins at the position (x, z) of crack; Aij (x , t ) is the cross-sectional area at time t, Aij (x , t ) =
hij 2 − hij 2
∫
G∇2 ui ′ +
G uk ′ , k ′ i ′ + Fi ′ = 0 (i′′k′ = 1,2) 1 − 2υ
(14)
Eq. (11) is expanded to obtain the width equation, it is represented by Eq. (15) (Bui, 1977; Adachi et al., 2007):
wij (x , t ) dz ; wij is the
fracture width. Substituting Eq. (10) into Eq. (11), the two-dimensional migration equation of the fracturing fluid in the fracture can be obtained:
(σn )ij − pij =
G 4π (1 − ν )
Fig. 5. Fracture morphology obtained by XFEM software.
299
∫Ω
ij
∂ ⎛ 1 ⎞ ∂wij ⎤ ⎡ ∂ ⎛ 1 ⎞ ∂wij ⎢ ∂x ⎜ l ⎟ ∂x ′ + ∂z ⎜ l ⎟ ∂z′ ⎥ dx ′dz′ ⎝ ij ⎠ ⎣ ⎝ ij ⎠ ⎦
(15)
Journal of Natural Gas Science and Engineering 52 (2018) 295–303
H. Zhao et al.
Fig. 6. The maximum width variation over time of the two methods.
Fig. 7. The results comparison of the two methods with the actual results.
For any branch of fracture network, the unit normal vector of the crack surface is:
→ n = [cos βj sin θj , cos βj cos θj , sin βj]
⎡ σn = → n⎢ ⎣
σh σH
⎤ →T 2 2 2 2 2 ⎥ n = cos βj cos θj σh + cos βj sin θj σH + sin βj σv σv ⎦ (17)
(16)
where Fi ′ is the external loads, here representing the resultant forces of fluid pressure and ground stress; ui ′ is the displacement of rock nodes in the direction i′, which is the crack width of the Y direction; θj is the angle between natural fracture and hydraulic fracture of the grade j, (°);
The normal stress acting on the crack surface is expressed as:
300
Journal of Natural Gas Science and Engineering 52 (2018) 295–303
H. Zhao et al.
βj is the dip angle of natural fracture of the grade j, (°); σh is the minimum horizontal principal stress, MPa; σH is the maximum horizontal principal stress, MPa; σv is the overburden pressure, MPa.
Finally, the dynamic process of network expansion is solved by the energy method. 4. Example analysis and discussion
3.2.2. The boundary conditions and initial conditions The width of the crack tip of each crack in the fracture network is zero, that is, w = 0 on Ωij ; The fracturing fluid flow at the crack seam is equal to the con∂p struction pump rate, that is, − S ∂n = −q0 on (Ωij ) p ; The fracturing fluid flow at the crack tip of each crack in the frac∂p ture network is zero, that is, − S ∂n = 0 on Ωij .
Combined with the geological and constructional parameters derived from the X2-HF well of the FL shale blocks in China, the formation process of complex fracture network is simulated by using the established method of the fracture network. The specific parameters of the example are shown in Table 1. The horizontal section of the well has been staged and clustered with sand adding fracturing technology, 15 stages, 3 clusters per stage, average cluster interval of 30m. The slickwater is selected as the main fracturing fluid. The natural fractures are randomly distributed, which can be entered manually. The distribution of natural fractures in subsurface reservoir is complicated and the distribution of different geological blocks is different. Before the numerical simulation, the parameters of natural fractures should be selected according to the geological data or the history data of the fractured wells. The pumping time should be determined according to the site construction requirements.
3.2.3. The interaction between hydraulic fracture and natural fracture Zhao et al. (2012) have discussed the interaction between hydraulic fracture and natural fracture. In order to analyze the complex fracture network formed by fracturing in the shale reservoir from the point of mathematics and physics, the theoretical deduction in this paper only considers the simultaneous expansion of hydraulic fractures along both ends of natural fractures. In this case, the fluid pressure at the intersection of the hydraulic fracture and the natural fracture needs to overcome the frictional resistance of the fracturing fluid from the intersection point to the end of the natural fracture and meet the condition of the end rupture of the natural fracture. We assumed fractures propagate along the maximum horizontal principal stress direction. When the hydraulic fracture encounters the natural fracture, it will propagate along natural fracture and turn to the maximum stress direction at the end of natural fracture. The turning radius of the fracture is neglected.
4.1. Simulation of fracture network expansion with the energy method 4.1.1. Simulation of fracture network dynamic expansion The computer program is used to simulate the morphology of fracture network propagation at different time. The dynamic forming process of the fracture network can be obtained. The injection time of fracturing fluid is determined to be 5min, 10min, 15min and 20min respectively. The other geological parameters and constructional parameters are selected according to the above data. The simulation results are shown in Fig. 3. From Fig. 3, the hydraulic cracks forward to communicate with the natural cracks. With the continuous injection of fracturing fluid, natural cracks will open and continue to extend from the both ends of the cracks, and the new hydraulic cracks will be formed. These new cracks spread forward to communicate the new natural cracks. Finally, a dendritic fracture network will be formed.
3.2.4. The flux distribution after hydraulic fractures interacting with natural fractures In order to directly study the propagation process of hydraulic fracture, the hydraulic fracture extends along the natural fracture when the interaction occurs. When the hydraulic fracture reaches the end of the natural crack, it will continue to expand along a certain angle. The concept of fracture classification is introduced, and the relationship between the flux (qN ) of the grade N branch crack and the pump rate (Q) is (Zhao et al., 2012): qN = Q/2N + 1.
4.1.2. Energy allocation during the expansion of fracture network In order to analyze the energy distribution with the injection time during the hydraulic fracturing process, the injection time of fracturing fluid is selected to be 10min. The calculation results of each part of energy are obtained from the software (Fig. 4). The rock strain energy accounts for 34.76% of the total energy input. The friction loss energy consumption in the fracture network accounts for 10.4% of the total input energy. The surface energy required to form the fracture network accounts for 0.015% of the total input energy. Therefore, the input energy is mainly used to overcome the work of ground stress in the large-scale fracturing process of the brittle shale reservoir. We concludes that the deeper the depth of the reservoir, the greater the energy consumed for the fracturing process. Additionally, the energy consumed by the friction resistance is relatively large when the fracturing fluid flows in the crack. Moreover, it can be seen that most of the energy remains at the moment. The cracks expand by the driving force of the residual energy. When the residual energy is not enough to overcome the energy required for the crack propagation, the crack will stop extending.
3.3. The analysis of numerical solution During the fracturing process in the field, the pump rate of fracturing fluid remains unchanged, and the fluid-solid coupling is achieved by the fluid pressure and fracture width of the fracture network. The specific method is to calculate the pressure distribution in fracture for each time step, and then calculate the width, flow distribution and energy distribution of fracture network. Then the energy criterion is used to determine whether the crack propagates in the time step. The width of the fracture network is recalculated at the end of the calculation. Each time step is calculated, and it is necessary to judge whether hydraulic fracture interacts with natural fracture. During the process of calculation, the width, flow distribution, pressure distribution, crack shape and other physical quantities are stored in real time. Therefore, fluid pressure inside the fracture affects the distribution of fracture width; in turn, the crack width changes the fracture shape and the fluid pressure distribution inside the crack. Then the fluid solid coupling between the rock and the fracturing fluid is achieved by the fluid pressure in the fracture network. The fracture network extension model established above is a large set of nonlinear equations. The solution of the model can only rely on numerical methods. In order to verify the model established in this paper and the rationality of the proposed energy method, we only consider the two-dimensional situation of the expansion of the fracture network when compiling the calculation program for the moment. And we assume that the crack height is constant and the filtration loss of fracturing fluid is also neglected.
4.1.3. Stimulated reservoir volume The stimulated reservoir volume (SRV) refers to the product of the area and thickness of the reservoir which is activated or communicated during fracturing (Mayerhofer et al., 2010). From the above analysis, it's known that the shape of fracture network formed by volume fracturing in shale reservoir is tree-type ellipsoid. The formula for SRV calculation is given: 301
Journal of Natural Gas Science and Engineering 52 (2018) 295–303
H. Zhao et al. M
SRV =
4
∑ ⎛3π × k=1
⎝
Wk
Lk × Hk ⎞ 2 × 2 2 ⎠
is considered. And the fluid filtration loss into rock mass is also ignored. Therefore, the two simplifications should be verified deeply in the next work.
(18)
where Wk is the maximum width of the fracture network, m; Lk is the maximum length of the fracture network, m; Hk is the maximum height of the fracture network, m; M is the cluster number of single stage fracture. Eq. (18) can be used to calculate the SRV of a single segment. The simulated results and the data observed by field micro-seismic monitoring are shown in Table 2. By contrast, the error of the SRV is about 13.7%. It demonstrates the validity of the energy criterion and the model established in this paper.
Acknowledgments This work was supported by National Natural Science Foundation of China (NO.11672333) and National Major Science and Technology Project of China (NO.2017ZX05009-003). The authors gratefully appreciate their supports. References
4.2. Comparison with the XFEM method
Adachi, J., Siebrits, E., Peirce, A., Desroches, J., 2007. Computer simulation of hydraulic fractures. Int. J. Rock Mech. Min. Sci. 44 (5), 739–754. Aimene, Y.E., Nairn, J.A., 2014. Modeling multiple hydraulic fractures interacting with natural fractures using the material point method. In: SPE/EAGE European Unconventional Resources Conference and Exhibition. Society of Petroleum Engineers. Amini, K., Soliman, M.Y., House, W.V., 2015. A three-dimensional thermal model for hydraulic fracturing. In: SPE Technical Conference and Exhibition. Bui, H.D., 1977. An integral equation method for solving the problem of a plane crack of arbitrary shape. J. Mech. Phys. Solid. 25, 29–39. Carter, R.D., 1957. Derivation of the General Equation for Estimating the Extent of the Fractured Area. Appendix I of “Optimum Fluid Characteristics for Fracture Extension,” Drilling and Production Practice, GC Howard and CR Fast. American Petroleum Institute, New York, New York, USA, pp. 261–269. Chun, K.H., Ghassemi, A., 2012. Fracture propagation under poroelastic loading. In: 46th US Rock Mechanics/Geomechanics Symposium. American Rock Mechanics Association. Fisher, M.K., Heinze, J.R., Harris, C.D., et al., 2004. Optimizing horizontal completion techniques in the Barnett shale using microseismic fracture mapping. In: SPE Annual Technical Conference and Exhibition. Society of Petroleum Engineers. Fisher, M.K., Wright, C.A., Davidson, B.M., et al., 2005. Integrating fracture mapping technologies to improve stimulations in the Barnett shale. SPE Prod. Facl 20 (2), 85–93. Hou, B., Chen, M., Zhang, B.W., 2015. Propagation of multiple hydraulic fractures in fractured shale reservoir. Chin. J. Geo. Eng 37 (6), 1041–1046. Irwin, G.R., 1957. Analysis of stresses and strains near the end of a crack traversing a plate. J. Appl. Mech. 24, 361–364. Khademian, Z., 2016. Studies of seismicity generated by unstable failures around circular excavations. In: 50th US Rock Mechanics/Geomechanics Symposium. American Rock Mechanics Association. Khademian, Z., Nakagawa, M., Garvey, R., et al., 2017. Role of fluid injection pressure in inducing seismicity. In: Stanford Geothermal Workshop, Stanford, California. Kozicki, W., Chou, C.H., Tiu, C., 1966. Non-Newtonian flow in ducts of arbitrary crosssectional shape. Chem. Eng. Sci. 21 (8), 665–679. Liu, J.J., Feng, X.T., Pei, G.H., 2003. Study on mathematical model of three dimensional hydraulic fracturing. Chin. J. Rock Mech. Eng. 22 (12), 2042–2046. Martin, A.N., 2000. Crack tip plasticity: a different approach to modeling fracture propagation in soft formations. In: SPE Annual Technical Conference and Exhibition. Society of Petroleum Engineers. Mayerhofer, M.J., Lolon, E., Warpinski, N.R., 2010. What is stimulated reservoir volume? SPE Prod. Oper. 25 (1), 89–98. Meyer, B.R., 1989. Three-dimensional hydraulic fracturing simulation on personal computers: theory and comparison studies. In: SPE Eastern Regional Meeting. Society of Petroleum Engineers. Meyer, B.R., Bazan, L.W., 2011. A discrete fracture network model for hydraulically induced fractures - theory, parametric and case studies. In: SPE Hydraulic Fracturing Technology Conference. Society of Petroleum Engineers. Nolte, K.G., 1979. Determination of fracture parameters from fracturing pressure decline. In: SPE Annual Technical Conference and Exhibition. Society of Petroleum Engineers. Nuismer, R.J., 1974. An energy release rate criterion for mixed mode fracture. Int. J. Fract. 11 (2), 245–250. Olson, J.E., Taleghani, A.D., 2009. Modeling simultaneous growth of multiple hydraulic fractures and their interaction with natural fractures. In: SPE Hydraulic Fracturing Technology Conference. Society of Petroleum Engineers. Salamon, M.D.G., 1984. Energy considerations in rock mechanics: fundamental results. J. S. Afr. Inst. Min. Metall 84 (8), 233–246. Taleghani, A.D., 2010. Fracture Reinitiation as a possible branching mechanism during hydraulic fracturing. In: 44th US Rock Mechanics Symposium and 5th US-Canada Rock Mechanics Symposium. American Rock Mechanics Association. Wang, L.X., Tang, D.H., Li, S.H., et al., 2015. Numerical simulation of hydraulic fracturing by a mixed method in two dimensions. Chin. J. Theor. Appl. Mech. 47 (6), 973–983. Weng, X., Kresse, O., Cohen, C.E., et al., 2011. Modeling of hydraulic fracture network propagation in a naturally fractured formation. SPE Prod. Oper. 26 (4), 368–380. Wu, Q., Xu, Y., Wang, X.Q., et al., 2012. Volume fracturing technology of unconventional reservoirs: connotation, optimization design and implementation. Petrol. Explor. Dev. 39 (3), 352–358. Xu, W., Thiercelin, M.J., Walton, I.C., 2009. Characterization of hydraulically-induced shale fracture network using an analytical/semi-analytical model. In: SPE Annual Technical Conference and Exhibition. Society of Petroleum Engineers.
In order to compare the calculation efficiency by using the energy method and the XFEM method to simulate crack propagation, the injection time of the fracturing fluid is selected to be 10min and other parameters are determined according to the above field data. The two different methods are applied to numerically simulate the crack propagation under the same conditions. One method is the simulation program compiled by the above established model, and another method is the finite element software of ABAQUS. The simulation results of fracture morphology are shown in Figs. 3(b) and 5. In Fig. 5, SDEG represents the damage degree of cohesive element and element is totally damaged when SDEG is equal to 1. The length of fracture network obtained by the developed model is 104.92m and the operation time is 5s. However, the crack length is 98.7m and the operation time is 9.5h calculated from ABAQUS software. The time difference of the two methods is 4 orders of magnitude, but the length error is 4.62%. The maximum opening width and the pressure needed for crack extension (Fig. 6 ∼7) are obtained by two different methods, which is extracted from the software backstage. The fracture width obtained from ABAQUS is a litter larger (4%–10%) than the width calculated from the method what we developed. The possible reason is that we assumed fluid can't filtrate into the shale, meaning that fluid only flows in the fracture network. From Fig. 7, the pressure results calculated by the two methods of the developed model and the ABAQUS software are close to the actual bottom hole pressure. The error is less than 1%. Analysis is that the possible reason for the error lies in the stress shadowing effects, and the interaction between fractures as they grow, the dissipated plastic work due to the slip along newly formed or preexisting fractures are neglected. In the energy method, the fracturing fluid flux and pressure are unified, then the expansion of fracture network is directly studied. The stress has little effect on the calculation results of the fracture propagation. However, in the finite element calculation, stress or displacement constraints need to be applied to the boundary when considering the influence of the formation factor. During the process of crack propagation, the stress changes at any moment. The stress value of each part of the reservoir needs to be recalculated when the cracks extend each step forward. The crack propagation is calculated on the basis of the stress, thus its calculation result is significantly influenced by the stress. Therefore, the calculation process of XFEM is complex and the computational efficiency is low. 5. Conclusion Based on the energy balance and the fluid-solid coupling theory during the hydraulic fracturing, an energy method for the simulation of network fracturing in shale gas horizontal wells is established. The dynamic expansion process of the fracture network is obtained by applying the energy method. The calculated results agree well with the actual data, indicating the rationality of the energy method. Simultaneously, the calculation efficiency of the energy method is more than 10000 times higher than that of XFEM. In this paper, only twodimensional fracture network expansion (the crack height is constant) 302
Journal of Natural Gas Science and Engineering 52 (2018) 295–303
H. Zhao et al. Yang, X.F., Chen, M., Liu, X.S., 2002. Theoretical study of the fluid field within the hydraulic fracture. J. Southwest Pet. Inst. 24 (01), 87–90. Yao, Y., Wang, W., Keer, L.M., 2018. An energy based analytical method to predict the influence of natural fractures on hydraulic fracture propagation. Eng. Fract. Mech. 189 (1), 232–245. Yew, C.H., Weng, X.W., 2015. Mechanics of Hydraulic Fracturing, second ed. Gulf Publishing Company, Houston, pp. 8–10. Zeng, Q.D., Yao, J., Sun, Z.X., 2015. Numerical modeling of fracture network propagation
in shale reservoirs. Chin. J. Theor. Appl. Mech. 47 (5), 994–999. Zhao, J.S., 2004. Fracture Mechanics and Fracture Physics. Huazhong University of Science and Technology Press, Wuhan, pp. 45–62. Zhao, H.F., Chen, M., Jin, Y., 2012. Rock fracture kinetics of the fracture mesh system in shale gas reservoirs. Petrol. Explor. Dev. 39 (4), 465–470. Zou, C.N., Yang, Z., Zhu, R.K., et al., 2015. Progress in China's unconventional oil and gas exploration and development and theoretical technologies. Acta Geol. Sin. 89 (6), 979–1007.
303
Contents lists available at ScienceDirect
Journal of Natural Gas Science and Engineering journal homepage: www.elsevier.com/locate/jngse
A simulation method based on energy criterion for network fracturing in shale gas reservoirs
T
Haifeng Zhaoa,∗, Xiaohua Wanga, Wei Wangb, Erfei Mua a b
School of Petroleum Engineering, China University of Petroleum, Beijing, 102249, China Hancheng Branch, PetroChina Coalbed Methane Co., Ltd., Hancheng, Shanxi, 715409, China
A R T I C L E I N F O
A B S T R A C T
Keywords: Shale gas Network fracturing Energy criterion Stimulated reservoir volume Fluid-solid coupling
The fracture network in shale gas reservoirs has a plurality of extending frontiers and the propagation paths of fractures are difficult to predict in advance. When the network fracturing is simulated by the extended finite element method (XFEM) based on stress intensity factors, the computation is arduous and time-consuming, which limits the application of XFEM in field fracturing design. Therefore, a simulating method of crack propagation based on energy criterion is proposed. Based on the principle of energy balance, the energy criterion for controlling crack frontiers extending is established by the fluid injection energy, the rock strain energy, the surface energy, the seismic energy and the friction loss energy. The criterion only needs to calculate the energy distribution in the crack frontiers to judge whether the fracture could expand, avoiding the fussy calculation of the stress intensity factor. Based on the proposed energy criterion, the fluid solid coupling effect and the fractures interaction, an energy method is established for network fracturing simulation in shale gas horizontal well. The results show that the fracture network morphology obtained by simulation is basically consistent with micro-seismic monitoring result, and the computational efficiency of this simulation method is 10000 times higher than that of XFEM. The energy method provides an idea and basis for the simulation and design of shale reservoir fracturing.
1. Introduction The horizontal well with segmented volume fracturing is the main technology to exploit the shale gas effectively. The complex network is formed during the hydraulic fracturing process for shale reservoirs via the micro-seismic monitoring. Thus, the investigation of the formation mechanism and simulation of fracture network is essential (Wu et al., 2012; Zou et al., 2015; Fisher et al., 2004, 2005) and the criterion of rock failure and fracture extension is the key step. The calculation complexity and time consumption of the fracture network models are significantly different when different criterions of fracture propagation are used in numerical simulations. Therefore, the criterion study of rock failure and fracture propagation for shale reservoir fracturing plays a vital role in characterizing the formation of the fracture network. Most of the current fracture propagation models are on the basis of stress intensity factor (SIF) (Irwin, 1957). Some complex fracture network models such as the wire-mesh model (Xu et al., 2009), the unconventional fracture model (Weng et al., 2011) and the discrete network model (Meyer, 1989, Meyer and Bazan, 2011) are established. These models are solved by different numerical algorithms, such as
∗
boundary element (Olson and Taleghani, 2009; Hou et al., 2015), extended finite element (Taleghani, 2010; Zeng et al., 2015) and discrete element (Meyer, 1989; Meyer and Bazan, 2011; Wang et al., 2015). The fracture criterions used in the mentioned models include the maximum energy release rate criterion, the minimum strain energy density criterion and the maximum tensile stress criterion (Zhao, 2004), and the theoretical foundation of these criterions is the SIF calculation based on the fracture toughness criterion. However, the fracture network has many propagation frontiers and the fracture geometry is very complex, making it difficult to calculate the stress distribution accurately. And the establishment of the integral path to calculate SIF is sometimes impossible (Martin, 2000). In addition, the calculation of network simulation is very time-wasting based on the fracture toughness criterion, which deviates from the engineering demand. Therefore, we attempt to study hydraulic crack propagation from the viewpoint of energy balance. Salamon (1984) established the energy balance law for rock failure processes during the mining course and developed methodologies for calculating energy terms in rock mechanics. Yao et al. (2018) assumed that the fractures extend forward when hydraulic energy is greater than the sum of strain energy and
Corresponding author. E-mail address: [email protected] (H. Zhao).
https://doi.org/10.1016/j.jngse.2018.02.003 Received 31 May 2017; Received in revised form 1 February 2018; Accepted 4 February 2018 Available online 07 February 2018 1875-5100/ © 2018 Elsevier B.V. All rights reserved.
Journal of Natural Gas Science and Engineering 52 (2018) 295–303
H. Zhao et al.
surface energy, but the friction dissipation energy was not taken into account when using the Griffith stability criterion. Khademian (2016); Khademian et al. (2017) focused on the radiated seismic energy during fracturing. They calculated strain energy, friction work, dissipative plastic work, and radiated seismic energy for analyzing rock failures in compression or shear. Forasmuch as, the energy criterion controlling the extension of fracture frontiers is established by balancing hydraulic fracturing injection energy, strain energy, surface energy, seismic energy and friction loss energy based on the above theories. Then an energy method for fracturing simulation of shale gas horizontal well is proposed combined with the energy criterion and the fluid-solid coupling method used to calculate fracturing fluid flow rate and pressure distribution in the fracture network. The energy method avoids the complex calculation of SIF in the crack frontiers of fracture network and provides an idea and basis for the fracturing simulation and design in shale reservoirs.
n
⎛
N
⎞
n
⎠
i=1
∑ Ek = (q − qlti) ⎜σmin + St + ∑ ∑ Δpij ⎟ ∑ ti i=1 j=0
⎝
(1)
The whole surface energy of fracture network is n
U=
N
∑ ∑ γs hij lij
(2)
i=1 j=0
The strain energy of rock is n
N
∑ Es = ∑ ∑ σmin hij lij wij (3)
i=1 j=0
The frictional loss energy is n
∑ Ef
= (q − qlti)
i=1
2. The energy mechanism and propagation criterion of fracture network formation
n
N
∑ ti ∑ ∑ Δpij i=1 j=0
(4)
where i represents the i-th time step in the iterative calculation; j ren presents the j-th level fracture; q is the pump rate, m3/s; T = ∑i = 1 ti is the injection time of fracturing fluid, s; σmin is the minimum principal stress, MPa; St is the tensile strength, MPa; Δpij is the pressure drop in the fracture, MPa; γs is the rock surface energy density, J/m2; hij is the fracture height, m; lij is the fracture length, m; wij is the average fracture width, wij = (wi − 1, j + wi, j )/2 , m,. According to previous studies, the reservoir rock is in the state of drastic change of local ground stress field. The I-II mixed fracture usually occurs during the process of hydraulic fracturing. For the twodimensional plane, the rock surface energy density (Nuismer, 1974; Aimene and Nairn, 2014) of the I-II type mixed fractures can be expressed as Eq. (5).
2.1. Energy mechanism of fracture network formation The hydraulic fracturing is a dynamic coupling process of rock deformation and fracturing fluid flow. On one hand, the micro-cracks occur inside the rock under the pressure of the fracturing fluid, and gradually propagate through the rock. Ultimately, the macro-cracks are formed. On the other hand, the fracturing fluid flows in the cracks and supports a certain width of cracks. With the high displacement injection of the fracturing fluid, the fluid pressure increases instantaneously, forcing the cracks to expand continuously. The hydraulic fracture would encounter natural fractures during the expansion process of fracture. When the hydraulic energy supplied by the fracturing fluid is greater than the energy required by natural crack propagation at both crack tips, the crack branches would be formed. Then the branch cracks continue to propagate. The branch cracks will form branches once again if the fractures run into the other natural fracture. In this way, the fracture network consisting of hydraulic cracks and natural fractures is eventually formed. Hydraulic fractures not only form branches with natural fractures, but also produce bifurcation, which finally produces the crack system like a root. In the shale reservoir, the complex fracture network shown in Fig. 1 is composed of natural fractures, hydraulic fractures and induced joints.
γs =
2 2 + KIIC (KIC ) 2E′
(5)
In the case of plane stress: E′ = E ; plane strain: E′ = E /(1 − υ2) ; KIC and KIIC are I and II rock fracture toughness respectively, which could be measured by laboratory experiments.
3. Dynamic mechanism and numerical solution of fracture network propagation 3.1. Dynamic mechanism of fracture network propagation
2.2. Energy criterion for crack propagation Dynamic mechanism of fracture network expansion was studied from the perspective of fracture dynamic mechanics (Zhao et al., 2012). When hydraulic fracture encounters the natural fractures, it is considered that the formation of complex fracture network needs to simultaneously meet the extending conditions at both ends of the natural fractures. The pressure required for the simultaneous propagation of both ends of the natural fracture is higher than the pressure required by the quasi-static expansion of the previous fracture. The fracture dynamics is used to determine the critical pump rate which could achieve the simultaneous expansion of branch cracks. The mechanism of the network expansion is studied from the perspective of energy in this paper. When hydraulic fracturing is carried out with large and variable pump rate, there is obvious dynamic shock wave in the stratum. A large number of micro seismic events is captured by the micro seismic monitoring technology. In the meantime, the stress state of rock and the tip position of crack vary with time rapidly. The inertial terms of the system balance equation can't be neglected when the fluid and rock are taken as a whole. Therefore, the crack propagation should be treated as a dynamical problem in this case.
During the hydraulic fracturing process, the energy in the cracks mainly contains the following parts (Fig. 2): the hydraulic energy of fracturing fluid (Ek); surface energy (U) to form new cracks; the rock strain energy (Es) which is equivalent to the work done by the fluid pressure; hydraulic friction loss energy (Ef) including the internal friction loss of fracturing fluid and friction loss between fracturing fluid and cracks wall; seismic energy (Ew) that monitored by micro-seismic during the fracturing process. Since micro-seismic energy is very small compared to the total hydraulic energy, it is ignored in our study. On the basis of the formation of the complex fracture network after fracturing in shale reservoir, without considering the bifurcation of the fracture itself, an energy criterion for the expansion of the whole fracture network under the condition of fracture classification is given. When ∑ Ek > ∑ (U + Es + Ef ) , the hydraulic fractures begin to propagate; when ∑ Ek < ∑ (U + Es + Ef ) , the hydraulic fractures cease to expand. The physical quantities in the energy criterion are calculated by Eqs. (1)–(4). Total injected hydraulic energy is 296
Journal of Natural Gas Science and Engineering 52 (2018) 295–303
H. Zhao et al.
Fig. 1. The propagation process of the fracture network.
3.2. The mathematical model of fracture network propagation
ρ
3.2.1. The governing equation
D→ υij Dt
pij + η (S ) ∇2→ υij + 2D⋅∇η = −∇→
(6)
where D is the strain rate tensor, D = (∇υ + ∇υT )/2 ; ρ is fluid density, kg/m3; η is the fluid viscosity, Pa·s; → υij is the fluid velocity vector; S = 2tr (D 2) . In the view of the actual situation of crack propagation, Eq. (6) is simplified. The pressure drop in the direction of fracture length and fracture height is represented by Eq. (7).
1) The fluid migration equation in the fracture The flow of fracturing fluid in a single fracture is simplified as a flow between two smooth parallel plates. The flow between the parallel plates follows the Navier-Stokes equation. It is assumed that the fracture fluid is incompressible liquid and the effect of mass force on the flow field is neglected. The general equation of the whole fracture network system is represented by Eq. (6) (Yang et al., 2002; Liu et al., 2003).
∂pij ∂x
=
∂ (υx )ij ⎤ ∂pij ∂ (υz )ij ⎤ ∂ ⎡ ∂ ⎡ , η (S ) η (S ) = ⎥ ⎢ y z y ∂y ⎢ ∂ ∂ ∂ ∂y ⎥ ⎣ ⎦ ⎣ ⎦
(7)
For power-law fluids, its viscosity follows the relation given by Eq.
Fig. 2. Diagram of energy distribution during the fracturing.
297
Journal of Natural Gas Science and Engineering 52 (2018) 295–303
H. Zhao et al.
For two-dimensional flow, the viscosity function η (S ) can be simplified:
Table 1 The specific parameters of the example. Parameters
Value
Parameters
Value
Overburden stress Maximum principal stress Minimum principal stress
63.5 MPa 58.5 MPa
Fracture toughness I Fracture toughness II
0.8 MPa m0.5 1.14 MPa m0.5
52.5 MPa
52.98J/m2
Young's modulus of rock Poisson's ratio Tensile strength
34.5 GPa 0.24 5 MPa
Rock strain energy density Fracturing fluid viscosity Average pump rate
⎜
⎟
⎜
⎟
(9)
The vector equation (6) is expanded and simplified. The volume flow rate of the fracturing fluid along the X and Z directions are expressed by fluid pressure pij and crack width wij (Kozicki et al., 1966; Nolte, 1979):
10 mPa s 13 m3/min
⎧ qx = − S
∂pij ∂x
⎨ q = −S ∂pij ∂z ⎩ z
(8).
η = K ′⋅S
2 (n ′− 1)/2
2
∂ (υz )ij ⎞ ⎤ ⎡ ∂ (υx )ij ⎞ η = K′⎢⎛ +⎛ ⎥ y ∂ ⎠ ⎝ ∂y ⎠ ⎦ ⎣⎝
(n ′− 1)/2
(8)
where S =
1 2n ′ K ′− n ′ 2n ′ + 1
Fig. 3. Fracture network shape at different time.
298
(10) wij
( ) 2
2n ′+ 1 n′ ⎡
⎢ ⎣
2
n ′− 1 2 − 2n ′
( ) + ( ) ⎤⎥⎦ ∂pij ∂x
∂pij ∂z
.
Journal of Natural Gas Science and Engineering 52 (2018) 295–303
H. Zhao et al.
Fig. 4. The energy variation over time during the fracturing.
∂Aij (x , t ) ∂ ⎛ ∂pij ⎞ ∂ ⎛ ∂pij ⎞ + ⎜S ⎟ + ⎜S ⎟ = ∂t ∂z ⎝ ∂z ⎠ ∂x ⎝ ∂x ⎠
Table 2 The geometric results of a single segment after fracturing. Geometric parameters
simulation results Micro-seismic interpretation results
Fracturenetwork length (m)
Fracturenetwork width (m)
365.34 334.53
92 101.41
Fracture height (m)
75 75.54
SRV/ (106m3)
1.32 1.53
Maximum width at main seam (mm)
∂ ∂x
5.4 -
=
where qij ∂ (q )ij
is
the
(12)
1 w (2n′+ 1)/ n′ 2n ′ K ′− n′ (n′+ 1)/ n′ ⎡ 2 ⎢ ⎨ 2n ′ + 1
⎣
⎧
2n ′
⎨ 2n ′ + 1 ⎩ ∂Aij (x , t ) ∂t
1 w
K ′− n′ +
(2n′+ 1)/ n′
2(n′+ 1)/ n′
⎡ ⎢ ⎣
2
n′− 1 2 − 2n′ ∂p ⎫ ij ⎤ + ∂x ⎬ ⎥
( ) ( )⎦ ∂pij
+
∂x
2
∂pij ∂z
( ) + ( ) ⎤⎦⎥ ∂pij ∂x
∂pij
⎭
n′− 1 2 − 2n′ ∂p ⎫ ij
∂z
∂z
⎬ ⎭
2cl hij
(13)
t − τ (x , z )
Therefore, the relationship between fluid pressure pij and width wij at different positions in the fracture network is established.
The continuity equation of the fracturing fluid inside the fracture is (Amini et al., 2015; Carter, 1957):
∂t
⎧ ⎩
∂ ∂z
∂Aij (x , t )
t − τ (x , z )
The above equations are expanded:
2) The continuity equation
∇⋅qij + (ql )ij +
2cl hij
3) The governing equation of rock deformation
=0
injection
(11) volume
rate
of
fracturing
Assuming that the rock deformation is elastic, the displacement governing equation of reservoir rock is expressed as (Chun and Ghassemi, 2012; Yew and Weng, 2015):
fluid,
∂ (q )ij
∇⋅qij = ∂xx + ∂zz ; qx and qz are volume flow rate in the X and Z directions respectively; (ql )ij is the filtration loss of fracturing fluid, (ql )ij = 2cl hij / t − τ (x , z ) ; cl is the fluid leakoff coefficient; τ (x , z ) is the time which the leak begins at the position (x, z) of crack; Aij (x , t ) is the cross-sectional area at time t, Aij (x , t ) =
hij 2 − hij 2
∫
G∇2 ui ′ +
G uk ′ , k ′ i ′ + Fi ′ = 0 (i′′k′ = 1,2) 1 − 2υ
(14)
Eq. (11) is expanded to obtain the width equation, it is represented by Eq. (15) (Bui, 1977; Adachi et al., 2007):
wij (x , t ) dz ; wij is the
fracture width. Substituting Eq. (10) into Eq. (11), the two-dimensional migration equation of the fracturing fluid in the fracture can be obtained:
(σn )ij − pij =
G 4π (1 − ν )
Fig. 5. Fracture morphology obtained by XFEM software.
299
∫Ω
ij
∂ ⎛ 1 ⎞ ∂wij ⎤ ⎡ ∂ ⎛ 1 ⎞ ∂wij ⎢ ∂x ⎜ l ⎟ ∂x ′ + ∂z ⎜ l ⎟ ∂z′ ⎥ dx ′dz′ ⎝ ij ⎠ ⎣ ⎝ ij ⎠ ⎦
(15)
Journal of Natural Gas Science and Engineering 52 (2018) 295–303
H. Zhao et al.
Fig. 6. The maximum width variation over time of the two methods.
Fig. 7. The results comparison of the two methods with the actual results.
For any branch of fracture network, the unit normal vector of the crack surface is:
→ n = [cos βj sin θj , cos βj cos θj , sin βj]
⎡ σn = → n⎢ ⎣
σh σH
⎤ →T 2 2 2 2 2 ⎥ n = cos βj cos θj σh + cos βj sin θj σH + sin βj σv σv ⎦ (17)
(16)
where Fi ′ is the external loads, here representing the resultant forces of fluid pressure and ground stress; ui ′ is the displacement of rock nodes in the direction i′, which is the crack width of the Y direction; θj is the angle between natural fracture and hydraulic fracture of the grade j, (°);
The normal stress acting on the crack surface is expressed as:
300
Journal of Natural Gas Science and Engineering 52 (2018) 295–303
H. Zhao et al.
βj is the dip angle of natural fracture of the grade j, (°); σh is the minimum horizontal principal stress, MPa; σH is the maximum horizontal principal stress, MPa; σv is the overburden pressure, MPa.
Finally, the dynamic process of network expansion is solved by the energy method. 4. Example analysis and discussion
3.2.2. The boundary conditions and initial conditions The width of the crack tip of each crack in the fracture network is zero, that is, w = 0 on Ωij ; The fracturing fluid flow at the crack seam is equal to the con∂p struction pump rate, that is, − S ∂n = −q0 on (Ωij ) p ; The fracturing fluid flow at the crack tip of each crack in the frac∂p ture network is zero, that is, − S ∂n = 0 on Ωij .
Combined with the geological and constructional parameters derived from the X2-HF well of the FL shale blocks in China, the formation process of complex fracture network is simulated by using the established method of the fracture network. The specific parameters of the example are shown in Table 1. The horizontal section of the well has been staged and clustered with sand adding fracturing technology, 15 stages, 3 clusters per stage, average cluster interval of 30m. The slickwater is selected as the main fracturing fluid. The natural fractures are randomly distributed, which can be entered manually. The distribution of natural fractures in subsurface reservoir is complicated and the distribution of different geological blocks is different. Before the numerical simulation, the parameters of natural fractures should be selected according to the geological data or the history data of the fractured wells. The pumping time should be determined according to the site construction requirements.
3.2.3. The interaction between hydraulic fracture and natural fracture Zhao et al. (2012) have discussed the interaction between hydraulic fracture and natural fracture. In order to analyze the complex fracture network formed by fracturing in the shale reservoir from the point of mathematics and physics, the theoretical deduction in this paper only considers the simultaneous expansion of hydraulic fractures along both ends of natural fractures. In this case, the fluid pressure at the intersection of the hydraulic fracture and the natural fracture needs to overcome the frictional resistance of the fracturing fluid from the intersection point to the end of the natural fracture and meet the condition of the end rupture of the natural fracture. We assumed fractures propagate along the maximum horizontal principal stress direction. When the hydraulic fracture encounters the natural fracture, it will propagate along natural fracture and turn to the maximum stress direction at the end of natural fracture. The turning radius of the fracture is neglected.
4.1. Simulation of fracture network expansion with the energy method 4.1.1. Simulation of fracture network dynamic expansion The computer program is used to simulate the morphology of fracture network propagation at different time. The dynamic forming process of the fracture network can be obtained. The injection time of fracturing fluid is determined to be 5min, 10min, 15min and 20min respectively. The other geological parameters and constructional parameters are selected according to the above data. The simulation results are shown in Fig. 3. From Fig. 3, the hydraulic cracks forward to communicate with the natural cracks. With the continuous injection of fracturing fluid, natural cracks will open and continue to extend from the both ends of the cracks, and the new hydraulic cracks will be formed. These new cracks spread forward to communicate the new natural cracks. Finally, a dendritic fracture network will be formed.
3.2.4. The flux distribution after hydraulic fractures interacting with natural fractures In order to directly study the propagation process of hydraulic fracture, the hydraulic fracture extends along the natural fracture when the interaction occurs. When the hydraulic fracture reaches the end of the natural crack, it will continue to expand along a certain angle. The concept of fracture classification is introduced, and the relationship between the flux (qN ) of the grade N branch crack and the pump rate (Q) is (Zhao et al., 2012): qN = Q/2N + 1.
4.1.2. Energy allocation during the expansion of fracture network In order to analyze the energy distribution with the injection time during the hydraulic fracturing process, the injection time of fracturing fluid is selected to be 10min. The calculation results of each part of energy are obtained from the software (Fig. 4). The rock strain energy accounts for 34.76% of the total energy input. The friction loss energy consumption in the fracture network accounts for 10.4% of the total input energy. The surface energy required to form the fracture network accounts for 0.015% of the total input energy. Therefore, the input energy is mainly used to overcome the work of ground stress in the large-scale fracturing process of the brittle shale reservoir. We concludes that the deeper the depth of the reservoir, the greater the energy consumed for the fracturing process. Additionally, the energy consumed by the friction resistance is relatively large when the fracturing fluid flows in the crack. Moreover, it can be seen that most of the energy remains at the moment. The cracks expand by the driving force of the residual energy. When the residual energy is not enough to overcome the energy required for the crack propagation, the crack will stop extending.
3.3. The analysis of numerical solution During the fracturing process in the field, the pump rate of fracturing fluid remains unchanged, and the fluid-solid coupling is achieved by the fluid pressure and fracture width of the fracture network. The specific method is to calculate the pressure distribution in fracture for each time step, and then calculate the width, flow distribution and energy distribution of fracture network. Then the energy criterion is used to determine whether the crack propagates in the time step. The width of the fracture network is recalculated at the end of the calculation. Each time step is calculated, and it is necessary to judge whether hydraulic fracture interacts with natural fracture. During the process of calculation, the width, flow distribution, pressure distribution, crack shape and other physical quantities are stored in real time. Therefore, fluid pressure inside the fracture affects the distribution of fracture width; in turn, the crack width changes the fracture shape and the fluid pressure distribution inside the crack. Then the fluid solid coupling between the rock and the fracturing fluid is achieved by the fluid pressure in the fracture network. The fracture network extension model established above is a large set of nonlinear equations. The solution of the model can only rely on numerical methods. In order to verify the model established in this paper and the rationality of the proposed energy method, we only consider the two-dimensional situation of the expansion of the fracture network when compiling the calculation program for the moment. And we assume that the crack height is constant and the filtration loss of fracturing fluid is also neglected.
4.1.3. Stimulated reservoir volume The stimulated reservoir volume (SRV) refers to the product of the area and thickness of the reservoir which is activated or communicated during fracturing (Mayerhofer et al., 2010). From the above analysis, it's known that the shape of fracture network formed by volume fracturing in shale reservoir is tree-type ellipsoid. The formula for SRV calculation is given: 301
Journal of Natural Gas Science and Engineering 52 (2018) 295–303
H. Zhao et al. M
SRV =
4
∑ ⎛3π × k=1
⎝
Wk
Lk × Hk ⎞ 2 × 2 2 ⎠
is considered. And the fluid filtration loss into rock mass is also ignored. Therefore, the two simplifications should be verified deeply in the next work.
(18)
where Wk is the maximum width of the fracture network, m; Lk is the maximum length of the fracture network, m; Hk is the maximum height of the fracture network, m; M is the cluster number of single stage fracture. Eq. (18) can be used to calculate the SRV of a single segment. The simulated results and the data observed by field micro-seismic monitoring are shown in Table 2. By contrast, the error of the SRV is about 13.7%. It demonstrates the validity of the energy criterion and the model established in this paper.
Acknowledgments This work was supported by National Natural Science Foundation of China (NO.11672333) and National Major Science and Technology Project of China (NO.2017ZX05009-003). The authors gratefully appreciate their supports. References
4.2. Comparison with the XFEM method
Adachi, J., Siebrits, E., Peirce, A., Desroches, J., 2007. Computer simulation of hydraulic fractures. Int. J. Rock Mech. Min. Sci. 44 (5), 739–754. Aimene, Y.E., Nairn, J.A., 2014. Modeling multiple hydraulic fractures interacting with natural fractures using the material point method. In: SPE/EAGE European Unconventional Resources Conference and Exhibition. Society of Petroleum Engineers. Amini, K., Soliman, M.Y., House, W.V., 2015. A three-dimensional thermal model for hydraulic fracturing. In: SPE Technical Conference and Exhibition. Bui, H.D., 1977. An integral equation method for solving the problem of a plane crack of arbitrary shape. J. Mech. Phys. Solid. 25, 29–39. Carter, R.D., 1957. Derivation of the General Equation for Estimating the Extent of the Fractured Area. Appendix I of “Optimum Fluid Characteristics for Fracture Extension,” Drilling and Production Practice, GC Howard and CR Fast. American Petroleum Institute, New York, New York, USA, pp. 261–269. Chun, K.H., Ghassemi, A., 2012. Fracture propagation under poroelastic loading. In: 46th US Rock Mechanics/Geomechanics Symposium. American Rock Mechanics Association. Fisher, M.K., Heinze, J.R., Harris, C.D., et al., 2004. Optimizing horizontal completion techniques in the Barnett shale using microseismic fracture mapping. In: SPE Annual Technical Conference and Exhibition. Society of Petroleum Engineers. Fisher, M.K., Wright, C.A., Davidson, B.M., et al., 2005. Integrating fracture mapping technologies to improve stimulations in the Barnett shale. SPE Prod. Facl 20 (2), 85–93. Hou, B., Chen, M., Zhang, B.W., 2015. Propagation of multiple hydraulic fractures in fractured shale reservoir. Chin. J. Geo. Eng 37 (6), 1041–1046. Irwin, G.R., 1957. Analysis of stresses and strains near the end of a crack traversing a plate. J. Appl. Mech. 24, 361–364. Khademian, Z., 2016. Studies of seismicity generated by unstable failures around circular excavations. In: 50th US Rock Mechanics/Geomechanics Symposium. American Rock Mechanics Association. Khademian, Z., Nakagawa, M., Garvey, R., et al., 2017. Role of fluid injection pressure in inducing seismicity. In: Stanford Geothermal Workshop, Stanford, California. Kozicki, W., Chou, C.H., Tiu, C., 1966. Non-Newtonian flow in ducts of arbitrary crosssectional shape. Chem. Eng. Sci. 21 (8), 665–679. Liu, J.J., Feng, X.T., Pei, G.H., 2003. Study on mathematical model of three dimensional hydraulic fracturing. Chin. J. Rock Mech. Eng. 22 (12), 2042–2046. Martin, A.N., 2000. Crack tip plasticity: a different approach to modeling fracture propagation in soft formations. In: SPE Annual Technical Conference and Exhibition. Society of Petroleum Engineers. Mayerhofer, M.J., Lolon, E., Warpinski, N.R., 2010. What is stimulated reservoir volume? SPE Prod. Oper. 25 (1), 89–98. Meyer, B.R., 1989. Three-dimensional hydraulic fracturing simulation on personal computers: theory and comparison studies. In: SPE Eastern Regional Meeting. Society of Petroleum Engineers. Meyer, B.R., Bazan, L.W., 2011. A discrete fracture network model for hydraulically induced fractures - theory, parametric and case studies. In: SPE Hydraulic Fracturing Technology Conference. Society of Petroleum Engineers. Nolte, K.G., 1979. Determination of fracture parameters from fracturing pressure decline. In: SPE Annual Technical Conference and Exhibition. Society of Petroleum Engineers. Nuismer, R.J., 1974. An energy release rate criterion for mixed mode fracture. Int. J. Fract. 11 (2), 245–250. Olson, J.E., Taleghani, A.D., 2009. Modeling simultaneous growth of multiple hydraulic fractures and their interaction with natural fractures. In: SPE Hydraulic Fracturing Technology Conference. Society of Petroleum Engineers. Salamon, M.D.G., 1984. Energy considerations in rock mechanics: fundamental results. J. S. Afr. Inst. Min. Metall 84 (8), 233–246. Taleghani, A.D., 2010. Fracture Reinitiation as a possible branching mechanism during hydraulic fracturing. In: 44th US Rock Mechanics Symposium and 5th US-Canada Rock Mechanics Symposium. American Rock Mechanics Association. Wang, L.X., Tang, D.H., Li, S.H., et al., 2015. Numerical simulation of hydraulic fracturing by a mixed method in two dimensions. Chin. J. Theor. Appl. Mech. 47 (6), 973–983. Weng, X., Kresse, O., Cohen, C.E., et al., 2011. Modeling of hydraulic fracture network propagation in a naturally fractured formation. SPE Prod. Oper. 26 (4), 368–380. Wu, Q., Xu, Y., Wang, X.Q., et al., 2012. Volume fracturing technology of unconventional reservoirs: connotation, optimization design and implementation. Petrol. Explor. Dev. 39 (3), 352–358. Xu, W., Thiercelin, M.J., Walton, I.C., 2009. Characterization of hydraulically-induced shale fracture network using an analytical/semi-analytical model. In: SPE Annual Technical Conference and Exhibition. Society of Petroleum Engineers.
In order to compare the calculation efficiency by using the energy method and the XFEM method to simulate crack propagation, the injection time of the fracturing fluid is selected to be 10min and other parameters are determined according to the above field data. The two different methods are applied to numerically simulate the crack propagation under the same conditions. One method is the simulation program compiled by the above established model, and another method is the finite element software of ABAQUS. The simulation results of fracture morphology are shown in Figs. 3(b) and 5. In Fig. 5, SDEG represents the damage degree of cohesive element and element is totally damaged when SDEG is equal to 1. The length of fracture network obtained by the developed model is 104.92m and the operation time is 5s. However, the crack length is 98.7m and the operation time is 9.5h calculated from ABAQUS software. The time difference of the two methods is 4 orders of magnitude, but the length error is 4.62%. The maximum opening width and the pressure needed for crack extension (Fig. 6 ∼7) are obtained by two different methods, which is extracted from the software backstage. The fracture width obtained from ABAQUS is a litter larger (4%–10%) than the width calculated from the method what we developed. The possible reason is that we assumed fluid can't filtrate into the shale, meaning that fluid only flows in the fracture network. From Fig. 7, the pressure results calculated by the two methods of the developed model and the ABAQUS software are close to the actual bottom hole pressure. The error is less than 1%. Analysis is that the possible reason for the error lies in the stress shadowing effects, and the interaction between fractures as they grow, the dissipated plastic work due to the slip along newly formed or preexisting fractures are neglected. In the energy method, the fracturing fluid flux and pressure are unified, then the expansion of fracture network is directly studied. The stress has little effect on the calculation results of the fracture propagation. However, in the finite element calculation, stress or displacement constraints need to be applied to the boundary when considering the influence of the formation factor. During the process of crack propagation, the stress changes at any moment. The stress value of each part of the reservoir needs to be recalculated when the cracks extend each step forward. The crack propagation is calculated on the basis of the stress, thus its calculation result is significantly influenced by the stress. Therefore, the calculation process of XFEM is complex and the computational efficiency is low. 5. Conclusion Based on the energy balance and the fluid-solid coupling theory during the hydraulic fracturing, an energy method for the simulation of network fracturing in shale gas horizontal wells is established. The dynamic expansion process of the fracture network is obtained by applying the energy method. The calculated results agree well with the actual data, indicating the rationality of the energy method. Simultaneously, the calculation efficiency of the energy method is more than 10000 times higher than that of XFEM. In this paper, only twodimensional fracture network expansion (the crack height is constant) 302
Journal of Natural Gas Science and Engineering 52 (2018) 295–303
H. Zhao et al. Yang, X.F., Chen, M., Liu, X.S., 2002. Theoretical study of the fluid field within the hydraulic fracture. J. Southwest Pet. Inst. 24 (01), 87–90. Yao, Y., Wang, W., Keer, L.M., 2018. An energy based analytical method to predict the influence of natural fractures on hydraulic fracture propagation. Eng. Fract. Mech. 189 (1), 232–245. Yew, C.H., Weng, X.W., 2015. Mechanics of Hydraulic Fracturing, second ed. Gulf Publishing Company, Houston, pp. 8–10. Zeng, Q.D., Yao, J., Sun, Z.X., 2015. Numerical modeling of fracture network propagation
in shale reservoirs. Chin. J. Theor. Appl. Mech. 47 (5), 994–999. Zhao, J.S., 2004. Fracture Mechanics and Fracture Physics. Huazhong University of Science and Technology Press, Wuhan, pp. 45–62. Zhao, H.F., Chen, M., Jin, Y., 2012. Rock fracture kinetics of the fracture mesh system in shale gas reservoirs. Petrol. Explor. Dev. 39 (4), 465–470. Zou, C.N., Yang, Z., Zhu, R.K., et al., 2015. Progress in China's unconventional oil and gas exploration and development and theoretical technologies. Acta Geol. Sin. 89 (6), 979–1007.
303