A coupled hydraulic-mechanical-damage geotechnical model for simulation of fracture propagation in geological media during hydraulic fracturing

Accepted Manuscript A coupled hydraulic-mechanical-damage geotechnical model for simulation of fracture propagation in geological media during hydraulic fracturing Tianjiao Li, Lianchong Li, Chun'an Tang, Zilin Zhang, Ming Li, Liaoyuan Zhang, Aishan Li PII:

S0920-4105(18)30972-0

DOI:

https://doi.org/10.1016/j.petrol.2018.10.104

Reference:

PETROL 5467

To appear in:

Journal of Petroleum Science and Engineering

Received Date: 27 April 2018 Revised Date:

21 October 2018

Accepted Date: 30 October 2018

Please cite this article as: Li, T., Li, L., Tang, Chun'., Zhang, Z., Li, M., Zhang, L., Li, A., A coupled hydraulic-mechanical-damage geotechnical model for simulation of fracture propagation in geological media during hydraulic fracturing, Journal of Petroleum Science and Engineering (2018), doi: https:// doi.org/10.1016/j.petrol.2018.10.104. This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

ACCEPTED MANUSCRIPT A Coupled Hydraulic-Mechanical-Damage Geotechnical Model for Simulation of Fracture Propagation in Geological Media During Hydraulic Fracturing Tianjiao Lia, Lianchong Lib,*, Chun’an Tanga, Zilin Zhangc, Ming Lic, Liaoyuan Zhangc, Aishan Lic a. School of Civil Engineering, Dalian University of Technology, Dalian 116024, China b. School of Resources and Civil Engineering, Northeastern University, Shenyang 110819, China c. Oil Production Technology Research Institute, Shengli Oilfield Branch Company, Dongying 257000, China

RI PT

Corresponding author. E-mail address: [email protected] (L. Li).

Abstract

This paper proposed a coupled hydraulic-mechanical-damage (HMD) geotechnical model for simulating the entire failure process in geological media. In this model, elastic damage mechanics

SC

describes the damage process of elements. Darcy’s Law and unsteady seepage flow equation describe fluid flow in elements, where compression of fluid and rock are considered. The coupled process is fulfilled by Biot’s consolidation theory describing how pore-pressure affects the stress field and the

M AN U

stress/aperture/mesh-dependent permeability that represents the influence of stress on seepage field. The HMD model is implemented in the calculation module of Rock-Failure-Process-Analysis on Petroleum-problems (RFPA-Petrol) codes. The RFPA suite of software codes is based on the idea that the heterogeneity leads to non-linearity and causes progressive failure behavior observed in brittle rock. The capability of RFPA-Petrol simulator is demonstrated by a typical coupled problem related to double fractures propagation during hydraulic fracturing. This simulation successfully reproduces the

TE D

asynchronous initiation, asymmetric growth and non-equal deflection of fractures in the numerical model. Then this simulator was used to investigate how fracture spacing and horizontal stress ratio influence propagation and reorientation of three hydraulic fractures from a horizontal well. At last, a horizontal well drilled by SHENGLI oilfield is simulated to certificate RFPA-Petrol’s ability on

EP

field-scale models. It is found that larger fracture spacing and higher initial fracture height will form longer and wider fractured zones. HMD-based RFPA-Petrol is suitable for simulation of laboratory

AC C

scale and field scale fracture propagation problems.

Keywords

HMD, numerical simulation, unsteady flow, hydraulic fracturing, stress shadow effect

1.

Introduction The productivity of oil and gas wells has been enhanced significantly in low-permeability

reservoirs, along with more technical and numerical developments in hydraulic fracturing (Li et al., 2011; Rahman et al., 2002), since the initial theoretical model for hydraulically induced fractures was proposed by Hubbert and Wills (Hubbert and Willis, 1972). However, the numerical simulation of coupled hydraulic-mechanical processes is a complicated problem with non-linear and non-local 1 / 24

ACCEPTED MANUSCRIPT properties enhanced by the discontinuities of the crack surface (Adachi et al., 2007; Lecampion, 2009; Sepehri et al., 2015). Throughout the last several years, except the Finite Element Method (FEM), several new numerical methods have been introduced to simulate the entire fracturing process without remeshing (MA et al., 2011; Sharafisafa and Nazem, 2014; Yang et al., 2016; Zheng et al., 2014), among which are the Extended Finite Element Method (XFEM) (Belytschko and Black, 1999),

RI PT

Discrete Element Method (DEM) ) (Cundall, 1971), Numerical Manifold Method (NMM) (Shi, 1992) and Rock Failure Process Analysis (RFPA) Method (Tang, 1997). These methods have been proven successful at simulating fracture propagation.

The extended finite element method (XFEM), is a numerical method for modeling cracks and crack growth within a standard finite element framework (Zheng and Luo, 2015). It enables a large

SC

class of crack growth problems except severely curved cracks to be solved without remeshing, which results in faster simulation of the fracture propagation (Belytschko and Black, 1999; Sepehri et al.,

M AN U

2015). It has been improved by embedding other methods. Sukumar et al. (Sukumar et al., 2003) combined the XFEM with the fast-marching method (FMM) for modeling three-dimensional planar crack growth. Khoei et al. (Khoei et al., 2016) added an equivalent continuum model to XFEM for simulating fluid flow through the deformable and porous media with multiple length scales of fractures. Lecamipon (Lecampion, 2009) took an internal fluid pressure inside fractures into consideration when applying the XFEM to hydraulic fracture problems. Furthermore, Li et al. (Li et al., 2015) integrated

TE D

the XFEM with the dual-permeability method (DPM) to present a multi-scale extended finite element model and investigate the multi-scale problem in shale gas reservoirs. The Discrete (or distinct) Element Method (DEM), is based on the generation of discontinuous particles interacting with each other (Marina et al., 2015). The discontinuities between particles can

EP

exhibit non-linear mechanical behaviors, including slip and opening (Damjanac and Cundall, 2016). Some of the applications of DEM to brittle materials problems were conducted (Hazeghian and Soroush, 2016; Marina et al., 2015; Prochazka, 2004; Tan et al., 2009). Sharafisafa and Nazem

AC C

(Sharafisafa and Nazem, 2014) used DEM to simulate fracture propagation and coalescence in fractured rock masses. Zeeb and Konietzky (Zeeb and Konietzky, 2015) applied DEM code 3DECTMto investigate stress shadow and fracture interaction in an Enhanced Geothermal System (EGS)under hydraulic fracturing. Deng et al. (Deng et al., 2014) deployed 3D DEM to simulate the shale–proppant interactions and evolution of fracture aperture during hydraulic fracturing. Damjanac and Cundall (Damjanac and Cundall, 2016) used fully coupled hydraulic-mechanical models to investigate hydraulic fracturing in naturally fractured rocks. Shimizu et al. (Shimizu et al., 2011) developed the flow-coupled DEM code to discuss the influence of the fluid viscosity and the particle size distribution during hydraulic fracturing in competent rock. The numerical manifold method (NMM), is originated from the topological manifold and 2 / 24

ACCEPTED MANUSCRIPT differential manifold (Ma et al., 2009). Moreover, its advantage is modeling discontinuous problems more naturally (Ma et al., 2009; Yang et al., 2016). This method has been improved by many researchers. Hu et al. (Hu et al., 2016) used the NMM model with a non-conforming mesh and a two-cover-mesh system to analyze flow in porous media with complexly intersecting fractures. He et al. (He et al., 2016) incorporated the hybrid crack element (HCE) method into the NMM to obtain the

RI PT

stress intensity factors and simulate crack propagation. Wu et al. (Wu et al., 2017)proposed a developed micro-mechanical based NMM which could closely mimic the deformation and failure characteristics of rock to investigate the rock macroscopic response and fracture processes. Zhang et al. (Zhang et al., 2015) utilized a loose coupling of NMM and the simulation of steady flow in a crack network to study the hydraulic fracturing process. Yang et al. (Yang et al., 2016) developed the NMM to analyze the

SC

direction of 3D fracture propagation. Zhang et al. (Zhang et al., 2016) developed a 3D fractured porous medium flow model by NMM to simulate fluid flow in arbitrarily and complexly fractured rock mass.

M AN U

In this paper, FEM based 3D simulator RFPA-Petrol is developed to simulate hydraulic fractures. It is based on the idea that the heterogeneity leads to macro non-linearity and causes progressive failure behavior observed in brittle rock (Tang, 1997; Tang et al., 1998). Therefore, global non-linear behaviors seen in brittle rock can be simulated with brittle-elastic elements differing in strength and physical constants depending on the heterogeneity index of the rock materials (Tang, 1997). When the pre-described strength is reached, the assumed parameters of the failed elements vary to describing its

TE D

weak and soft mechanical properties (Tang, 1997). During the progressive failure process, failed elements connect to form a fracture, and no sliding interface or additional "gap" element is needed (Tang, 1997). RFPA has excellent characteristics. It does not require remeshing after each time step or adding any artificial boundary element when the fracture propagates. It results in lower computation

EP

time and is well-suited for modeling fractures. RFPA’s capabilities of investigating hydraulic fracture problems have been demonstrated by many authors (He et al., 2017; Li et al., 2016; Li et al., 2016; Li

2.

AC C

et al., 2013; Li et al., 2017; Wang et al., 2009). Simulation methodology

2.1 Basic assumptions

RFPA with FEM capability was implemented as a basic simulator for geomechanical analysis. In

RFPA-Petrol, the eight-node isoparametric hexahedral element is used as the basic element (Li et al., 2013). In this study, it is developed with the following assumptions: A fractured porous rock is assumed to be an inhomogeneous, isothermal, deformable fluid-saturated porous medium. Its heterogeneity is considered by assuming that the properties, such as Young’s modulus and the strength properties of each element, conform to the Weibull distribution (specified by the Weibull distribution parameters) (Li et al., 2013; Tang et al., 1998). At the elemental 3 / 24

ACCEPTED MANUSCRIPT scale, the rock mass is treated as an elastobrittle material with a residual strength (Li et al., 2012; Li et al., 2013). Its mechanical behavior is modeled by an elastic damage constitutive law, and the residual strain/deformation upon unloading is not considered (Li et al., 2012; Li et al., 2013). Additionally, the coupled process between stress/strain and fluid flow in the deforming rock mass is governed by Biot’s consolidation theory (Biot, 1941).

RI PT

The compressible fracturing fluid filling in the porous rock is governed by Darcy’s law. The isotropic conditions are considered for the hydraulic behavior at elemental scale. The hydraulic behavior responds to changes in stress field by permeability variation, i.e., the permeability of an element varies as a function of the stress state during elastic deformation(Li et al., 2016), increases according to a deformation-dependent law when elements are damaged and obey assumption of pipe

SC

laminar flow for failure elements. The differences of permeability induced by the compressibility of

2.2 Static mechanical equations

M AN U

water are neglected.

In the RFPA-Petrol model, the elastic mechanics can be used to describe rock deformation at the elemental scale before damage (Zhou and Hou, 2013). So, the constitutive relations of an isotropic linear poroelastic medium based on continuum mechanics can be expressed as follows:

2ν G tr ( ε ) I − α pI 1 − 2ν

TE D

σ = 2Gε +

(2-1)

The equations of equilibrium and the strain–displacement relations can be written in the general form:

∇σ + f = 0

AC C

EP

and

ε=

1 (( ∇ u ) T + ∇ u ) 2

(2-2)

(2-3)

2.3 Fluid flow equations

In the case of elements with isotropically permeability, the fluid conductivity equation is governed

by (Freeze and Cherry, 1980):

ρ gκ 2 ∂H ∇ H −q = S µ ∂t

(2-4)

S is storativity, also referred to as storage coefficient. It is the volume of water that an aquifer releases from or takes into storage per unit surface area per unit change in the water head (Fader, 1974). It is determined by the compressibility of matrix and fluid (Freeze and Cherry, 1980; Sneddon and Elliot, 1946):

4 / 24

ACCEPTED MANUSCRIPT S = ρ g ( cM + n β )

(2-5)

The compressibility of the matrix can be considered by experimentally obtained pressure-void ratio curve or approximately determined by the vertical uniaxial compressibility as the following expression (Ferronato et al., 2010): cM = ( λ +2G ) =

(1 + ν )(1 − 2ν ) 1 1 −ν

(2-6)

E

RI PT

-1

2.4 Constitutive law of elements

Mechanical properties and permeability characteristics of rock component are different. They lead

SC

to great spatial heterogeneity in mesoscopic properties of rock (Li et al., 2013). In RFPA-Petrol modeling, the material is composed of elements whose material parameters (such as Young’s modulus

2013; Liang et al., 2004), as:

f ( x) =

M AN U

and strength) is assumed to follow a Weibull function with threshold values (Li et al., 2016; Li et al.,

m x    x0  x0 

m −1

 x  exp  −   x0 

m

(2-7)

According to this function, the main physical parameters inputted in RFPA-Petrol are determined.

TE D

2.5 Damage evolution of elements

An Element assumed to be elastobrittle with a residual strength has only one failure mode. The compressive stress is supposed to be positive. Its stress-strain relationship is shown in Fig.1. It is considered to failure in a tensile mode when the minimum principal strain exceeds the experimentally

EP

determined uniaxial tensile strain at failure or the minimum principal stress exceeds the uniaxial tensile strength, and fails in a shear mode when the shear stress satisfies the Mohr-Coulomb failure criterion

AC C

(Rankine, 2013). The criterions (including the maximum strain criterion, the maximum stress criterion, and the Mohr-Coulomb failure criterion) are expressed as flowing in equations:

or

ε3 ≤

σ t0 E

σ3 ≤ σt0

(2-8)

(2-9)

or

σ1 − σ 3

1 + sin φ ≥ σ c0 1 − sin φ

(2-10)

When an element’s state satisfies any one of these criterions, it begins to accumulate damage. In elastic damage mechanics, the elastic modulus of the element gradually declines with damage. Based on the strain equivalence hypothesis, the elastic modulus of the damaged element is given by: 5 / 24

ACCEPTED MANUSCRIPT Ed = (1 − D ) E0

(2-11)

2.5.1 Damage evolution of an element in the tensile state The damage variable of an element in a tensile mode under multiaxial stress states can be

ε > εt 0

 0   σ D = 1 − rt  ε E0  1 

ε tu < ε ≤ ε t 0 ε ≤ ε tu

SC

ε tu =ηε t 0

RI PT

expressed as (Li et al., 2013; Li et al., 2017; Li et al., 2018; Liang et al., 2004):

(2-12)

(2-13)

where ε is the equivalent principal strain. It is defined as follows (Li et al., 2017; Zhou et al., 1998):

−ε12 + −ε 22 + −ε 32

M AN U

ε =−

where ε 1 , ε 2 and ε 3 are the three principal strains, and

x x = 0

(2-14)

is a function defined as follows:

x≥0 x<0

(2-15)

TE D

2.5.2 Damage evolution of an element in the compression state The damage variable of an element in a shear mode under multiaxial stress states can be expressed as (Li et al., 2013; Li et al., 2017; Li et al., 2018; Liang et al., 2004):

EP

 0  D= σ rc 1 − ε E 1 0 

ε1 < ε c 0 ε c 0 ≤ ε1

(2-16)

AC C

In RFPA-Petrol, the state of elements can be classified by damage variable into three kinds that are elastic elements (D=0), damaged elements (0
ACCEPTED MANUSCRIPT smooth parallel sides (Li et al., 2013), which means that the permeability of a damaged element is isotropic as shown in Fig.4 (Yuan and Harrison, 2005). In order to achieve a simple approximation for hydraulic fracture aperture, the fracture whose length is assumed to be much greater than width may be regarded as a confined-height fracture. Also, consider that the fracture is being opened by an internal pressure that equal to the closure stress which

RI PT

holds the fracture closed plus a net pressure (Economides and Michael, 2007). Under these assumptions and conditions, fracture aperture is represented as a linear function of net pressure (Economides and Michael, 2007; Sneddon and Elliott, 1946): w=

2 pnet h f (1 − υ 2 ) E

(2-17)

SC

The net pressure could be represented by the changes in the normal stress on the fracture wall (Zhou and Hou, 2013), so it is determined by the pressure inside the fracture and the minimum in-situ stress as follows:

2.7 Permeability variation of elements

M AN U

pnet = p f − σ c

(2-18)

The permeability of an element varies as functions of the stress state during elastic deformation, increases according to an aperture-dependent law when elements are damaged and obey assumption of

TE D

pipe laminar flow for full failure elements.

2.7.1 Stress-dependent permeability of elastic elements

EP

During elastic deformation, the reduction of strength and stiffness lead to an increase in permeability of element by causing dilatant deformation of rock (Li et al., 2016; Li et al., 2013; Zhu et al., 2006). The empirical law of permeability variation for an elastic element depending on the stress

AC C

history can be described by (Louis, 1975; Tang et al., 2002)

κ e = κ 0e

σ  − b  ii −α P   3 

(2-19)

In RFPA-Petrol, Eq. (2-19) applies to the stress-dependent permeability of elastic elements (Li et al., 2013).

2.7.2 Aperture-dependent permeability of damaged elements For elements fail in shear mode, the rearrangement of grains and pores due to shear slip affects permeability slightly (Crawford and Yale, 2002). Therefore, the shear strain’s influence on permeability is neglected, and the permeability of damaged elements undergo shear obeys Eq. (2-19). When an element fails in tensile mode and begins to accumulate damage, fluid flow through 7 / 24

ACCEPTED MANUSCRIPT damage channels is simplified as laminar flow between parallel plates where the volumetric flow per unit width normal to the direction of flow varies as the cube of the fracture aperture (Klimczak et al., 2010; Oron and Berkowitz, 1998; Witherspoon et al., 1980). It is found to be valid whether the fracture surfaces are held open or are being closed under stress, and it is fairly accurate for large apertures where fracture roughness can be neglected (Sisavath et al., 2003; Witherspoon et al., 1980). So, we

RI PT

assume that a damaged rock element (when D>0) may be represented hydraulically as a unit of rock containing a fracture with smooth parallel sides (Li et al., 2013). The cubic law gives the unit discharge in a fracture by the following:

q0 =

d 03 ρ g ρg J f = κd ⋅ ⋅ A⋅ J f 12µ µ

(2-20)

SC

So, permeability for a damaged element in tensile state uniquely defined by fracture aperture can be expressed as follows:

d0 2 12

M AN U

κd =

(2-21)

d0 which is the distance between two parallel sides is equal to the aperture of the element. Then, d0 could be calculated by:

do ≈ w =

2 pnet h f (1 − ν 2 ) E

(2-22)

TE D

2.7.3 mesh-dependent permeability of full failure elements

For complete failure elements in the tensile state, fluid flowing through it is conceptualized by using the assumption of pipe laminar flow. Therefore, the permeability of a complete failure element

EP

can be determined as follows (Chen et al., 2004):

κf =

df2 32

(2-23)

The equivalent diameter df can be calculated by side length of the element by df2=4a2/π, where a is the

AC C

mesh size.

In hydraulic fracturing, most of the failure elements are in the tensile state. Eq. (2-21) and Eq.

(2-23) describe the permeability variation for damaged elements and complete failure elements separately.

2.8 Numerical implementation of HMD model in RFPA-Petrol The above equations for unsteady seepage, static mechanics, and elastic damage mechanics are utilized in the geomechanical simulator RFPA-Petrol. Fig.5 illustrates the numerical analysis process and computation scheme of this simulator. This program consists of three main modules which are modeling module, calculation module, and damage analysis module. Among them, the calculation 8 / 24

ACCEPTED MANUSCRIPT module is written in FORTRAN language. The pre- and post- processing modules, i.e., modeling module and damage analysis module, are implemented by Visual C++ language in RFPA (Li et al., 2016).

3.

The application of RFPA-Petrol on parallel multiple fracture analysis

RI PT

3.1 A basic verification of the capability of RFPA-Petrol with analytical solutions In this section, a thick-wall cylinder was employed as the basic verification model to validate the capabilities of the proposed program RFPA-Petrol. Based on the plane strain assumption, the displacement in height was set as 0, and the model can be simplified as shown in Fig. 6. This model

SC

with an inside diameter of 5 mm and an outside diameter of 25 mm was under a constant inner fluid pressure of 1.609MPa at the borehole wall (Zhan et al., 2009). Other parameters related to this validation example is in Table1. Based on neglecting the convective heat flux, a steady-state

M AN U

calculation was conducted to investigate the variation of stresses in this cylinder.

The stress state at any point in the cross-section under pore-pressure is determined by the difference of internal pore-pressure and external pressure rather than the external pressure itself. Therefore, the radial and tangential stresses induced by hydraulic pressure at any point are formulated as follows (Zhan et al., 2009):

(1 + υ )( ra rb ) 2 ( rb2 − ra2 )

r r  rb  ln r ln r ln r Sk  2 a + 2 b − 2 a  r ra rb  

TE D

σr =

2

2 ( rb2 − ra2 )

AC C

EP

σθ =

(1 + υ )( ra rb )

2

rb r r   ln r ln r ln r Sk  − 2 a + 2 b − 2 a  r ra rb   Sk =

     

   − 1 − υ Sk 2   

Pb − Pa r ln b ra

(3-1)

(3-2)

(3-3)

The stress maps obtained from numerical simulation are illustrated in Fig.7. It describes that

pore-pressure and minimum principal stress are mainly concentrated within the zones around the hole. They emerge gradient distribution along the radius. While the maximum principal stress firstly increased then decreased with the radius. Based on the plane strain assumption, the maximum principal stress equal to the radial stress and the minimum principal stress equal to the tangential stress. Fig.8 shows the analytically and numerically obtained radial and tangential stress variation with radius at θ= 0. The analytical stress values are calculated from Eq. (3-1) to Eq. (3-3). As can be seen from the line graph, the tangential stress is largest in magnitude on the wall of the hole. Then it declines sharply and reaches 0.5MPa at the edge of the model. The simulated radial and tangential stress are 9 / 24

ACCEPTED MANUSCRIPT compared with the analytical results and exhibited good agreement in most regions. As RFPA-Petrol uses eight-node isoparametric hexahedral elements to mesh, the tangential stress has a relative error of 2% when r equals 2.5mm. 3.2 A verification of the capability of RFPA-Petrol in multiple fracture propagation with laboratory

RI PT

results A model demonstrating the capability of RFPA-Petrol in simulating processes of multiple fracture propagation was built based on some hydraulic fracturing experiments (Rahman et al., 2002). The laboratory tests were designed to investigate transverse fractures propagation from two parallel circular

SC

perforations under hydraulic fracturing. Dimensions and parameters in the simulation as shown in Table 2 were selected to be similar to the laboratory test DC-3 (Rahman et al., 2002) for comparing the numerical results with the laboratory results.

M AN U

In order to evaluate elastic modulus and strength that would be used in the simulation of multiple hydraulic fracturing, a uniaxial compression numerical simulation was built. The numerical input elastic modulus, strength and degree of material homogeneity were calculated by Eq. (2-7). As shown in Fig.9 the model measuring 50×50×100mm was meshed by 50×50×100 elements with a side length of 1 mm. The top uniaxial displacement loading rate was 0.003mm/step.

Fig.10 shows the strain-stress curve for m=4. Table 3 presents macroscopic physical parameters,

TE D

microscopic physical parameters and the physical parameters obtained from the numerical test. Macroscopic physical parameters are compared with physical output parameters. It is found that the absolute error of elastic modulus and strength are 1.3% and 2.9% respectively. In sum, the numerical input physical parameters are suitable for simulating multiple hydraulic fracturing.

EP

The laboratory hydraulic fracturing sample measuring 400×400×400mm was confined in a triaxial loading setup where the principal stresses of 7MPa vertical, 6MPa maximum horizontal, and 5MPa

AC C

minimum horizontal stresses were applied. Two perforations were simulated by two parallel plastic discs are perpendicular to the direction of minimum horizontal stress and are separated by dx = 0, dy=0, dz= 30mm (Rahman et al., 2002). Fig.11 illustrates the cutaway view of simplified laboratory sample and cross-section of the

numerical model. Fig.11(b) shows that there is a grouted hole between two perforations. It was filled by grout, so it was not considered in the simulated model as shown in Fig.11(c). As mesh size is 4mm in simulation, dz must be a multiple of 4mm. dz= 32mm was chosen in the simulation. In this model, the y-axis is the direction of the vertical stress. In the simulation, the injection rates in two perforations were same at the beginning. Then injection into primary fracture was ceased and resumed with the change of primary fracture propagation state and secondary fracture propagation state. Injection into perforations was both shut in 10 / 24

ACCEPTED MANUSCRIPT after the 2000s. The injection rate for the numerical test is illustrated in Fig.12(b). Fig.12 shows that the numerical pressure responses resemble the laboratory results especially in variation tendency although the discrepancy between them is still considerable. The RFPA-Petrol simulation successfully obtains the second peak of primary fracture injection pressure without special artificial treatment. Fig.12 shows the pressure response to fractures in the entire hydraulic fracturing

RI PT

process. The injection pressures of two perforations both increases slightly from 0MPa to 0.8MPa between 0s to 500s. Then they have a dramatic increase throughout 500s to 920s, rising from 0.8MPa to 21.19MPa and 11.9MPa respectively. During this period, the difference in pressure of two perforations grows gradually. Then the injection pressure in one of the perforations reaches the breakdown pressure and primary fracture initiates. After peaking at 21.19MPa when t=920s, the

SC

pressure of primary fracture had fallen back to 15.8MPa by 1360s due to the decrease in primary injection rate and fracture propagation. Just before it, the pressure of another perforation reaches the

M AN U

breakdown pressure which equaled to 21.9MPa and the secondary fracture initiates. Then there is a rapidly decrease for injection pressure in secondary fracture from 21.9MPa to 13.9MPa between the 1160s and 2000s. However, injection pressure in primary fracture increases to 17.09MPa and reaches its second peak at the 1560s and has a sharp decrease from 17.09MPa to 15.2MPa throughout the 1560s to 2000s. During this period, the difference between primary fracture curve and secondary fracture curve becomes small, but they do not overlap because the grouted hole connecting wellbores are not

TE D

simulated in the numerical model and region between two pre-existing fractures are not damaged and connected. After t=2000s, the injection pressure of primary fracture and secondary pressure decrease gradually and reach 10MPa and 6MPa respectively which are the extension pressure at the end of the calculation. The pressure in the primary fracture has a sharp drop because it reaches the model

EP

boundary.

To sum up, the variation tendency of numerical injection pressure curves is similar to the experiment curves. The breakdown pressure of the secondary fracture is larger than that of primary

AC C

fracture as the secondary fracture initiates in the presence of the primary fracture. The second peak of primary fracture curve is lower than the breakdown pressures of both fractures. The injection pressures of them declined quickly after fracture initiation due to fracture growth and subsequent leak-off. Fracturing two adjacent perforations that are parallel to each other produces a complex pressure

field between and near them. Table 4 shows that the minimum stress concentrating on the two perforations contributes to an increase in the magnitude of tensile stress acting on them before primary fractures initiating. Then the primary fracture initiates. If the injection of it did not shut down at that time, it was hard for the secondary fracture to propagate. After the injection into the primary fracture is ceased, the tensile stress in the other perforation begins to accumulate, and the secondary fracture cracks finally. 11 / 24

ACCEPTED MANUSCRIPT Table 4 illustrates the evolution of the minimum stress and pore-pressure in the cubic numerical model. It is cut by x=0.2m, y=0.2m, and z=0.1m in order to show the results clearly. The moments of outputting images are marked as the key points (KP) in Fig.12. From Point A to Point B, pore-pressure is small, and the stress field is determined by the in-situ stresses. Then stress reversal caused by the increasing pore-pressure appears in the area near two perforations between Point B and Point C. At

RI PT

Point C, pore-pressure in the two perforations are not equal anymore. The perforation with larger pore-pressure starts to crack as the primary hydraulic fracture. Its altered local distribution of stresses significantly. The minimum stress in the area between two perforations decreases while the minimum stress on the lower side of primary fracture increases. Then the injection of lower perforation is shut down. However, the primary fracture does not stop propagating caused by stress concentration at the

SC

crack tip. Though Point C to Point D, the minimum stress in the upper perforation decreases quickly because of the increasing pore-pressure. At Point D it reaches the tensile strength and the upper

M AN U

perforation cracks as the secondary fracture. It also causes stress redistribution resembling the primary fracture. Four high compressive stress zones which lead to the reorientation of fractures are shown clearly. From Point D to Point E, although pore-pressures in both perforations’ declines, the tensile stresses concentrate on the area near crack tips of two fractures. Therefore, they keep propagation. However, the propagation velocity of secondary fracture is faster due to larger pore-pressure in it. At Point E, the pore-pressure in lower perforation reaches the value of the other perforation quickly

TE D

because of the reopening of injection in it. Moreover, the minimum stress at its crack tip decreases resulting in the improvement of its propagation rate and reaches the magnitude of tensile strength in upper perforation at Point F. During Point F to Point H, the minimum stress of the elements near the fracture except elements near crack tips increases gradually and nearly changes from negative to

EP

positive. The primary fracture reaches the model boundary at Point H, and the calculation is stopped. During the numerical cubic model cracking, two fractures grow perpendicular to the z-axis and deviate from the opposing fracture planes instead of propagating perpendicular to the direction of

AC C

z-axis which is the direction of far-field minimum horizontal stress. This phenomenon is caused by the stress shadow effect (the increase in the minimum horizontal stress by adjacent hydraulic fracture (Rios et al., 2013)). The time-dependent fracture geometry and aperture are illustrated in Fig.13. It is shown

that the two fractures are characterized by highly asymmetric propagation. They form V-shape spalled zones that are away from each other. As shown in Fig.13, the numerical fracture geometry is not on the symmetry of plane x-y during the entire damage process. The differences between primary fracture and secondary fracture are showed. The average initial deflection angle of primary fracture is approximately 20 degrees, whereas secondary fracture’s average deflection angle is about 40 degrees. Compared with primary fracture, the curvatures of the secondary fracture are smaller. This is caused by redistribution of stresses after the primary fracture initiating which strongly interferes the other 12 / 24

ACCEPTED MANUSCRIPT perforation. Moreover, comparing the right part with the left part in Fig.14, both the deflection angle of primary fracture and deflection angle of secondary in the left side are larger than that in the right side. It is shown that the larger deflection angle of primary fractures is 29 degrees, and the larger deflection angle of secondary fracture is 45 degrees. Fracture propagation paths on the 3D numerical model cross section are compared to available

RI PT

photographs gained from laboratory test DC-3 sample (Rahman et al., 2002). In order to compare experimentally obtained deflection angles with simulation deflection angles clearly, the picture of the laboratory model cross section is rotated clockwise by 180 degrees. Fig.14 clearly shows the similarities and differences between laboratory results and numerical results. The simulation reproduces the propagation path obtained from the experiment. The average deflection angle of

SC

secondary fracture calculated by RFPA-Petrol is close to the experiment result with a relative error of 2%. However, there is a significant error of deflection angle of primary fracture in -x side between

M AN U

numerical result and experiment result, which can be explained by the following reason. The pressure of primary fracture declined slower after damaged in numerical simulation than that in the laboratory. In sum, RFPA-Petrol can simulate the dynamic propagation of fracture. During the calculation, the leak-off was considered naturally without additional equations. Comparing the fracture geometry obtained from RFPA-Petrol simulation and experiment, it is shown that RFPA-Petrol can predict the deflection of double fractures and reflect their asymmetry during multiple fracturing stimulation.

TE D

This research has attempted to simulate the propagation of parallel fractures from horizontal wells. Turning to a discussion of the results under the condition that the distance between two perforation plates is approximately 1.5 times the half-length of perforations. First, the breakdown pressure of the primary fracture is slightly larger than the secondary fracture. Second, the secondary fracture has one

EP

peak while the primary fracture has two peaks because of shutting down and restarting its injection. Third, the second peak value of injection pressure of primary fracture is smaller than its first peak and the value of secondary fracture injection pressure at that time. Forth, the two fractures are characterized

AC C

by highly asymmetric propagation. Compared with primary fracture, the curvatures of the secondary fracture are larger. Overall, this simulation shows agreement with laboratory results under a similar range of conditions.

Simulated fracture propagation paths (Bryant et al., 2015; Rahman et al., 2002; Wang, 2016; Wu

et al., 2012) that incorporate the HMD model is compared with other numerical results of multiple fracture propagation in Fig.15. The maximum horizontal in-situ stress is in the vertical direction. Although the dimensions of model c, d and e are different from Rahman et al.’s laboratory tests (Rahman et al., 2002) and the numerical model in this research, they simulated the same process that is

stress interference between two parallel fractures under hydraulic fracturing. Comparing the HMD fracture geometry with other numerical results, it is found that HMD results show strong asymmetry as 13 / 24

ACCEPTED MANUSCRIPT the laboratory results due to the consideration of heterogeneity and accumulation of damage. However, the discrepancy between numerically obtained results and experimentally obtained results cannot be ignored. It is mainly due to mesh sensitivity, the imprecise empirical failure element factors, and model simplification. The RFPA-Petrol simulation slightly overestimates the failure zone as the element size which is not small enough to draw aperture and calculate permeability after damage.

RI PT

Injection pressure vs. injection period curve for the numerical test has a slower decline after initiation of fractures because the confined-height assumption of aperture leads to the underestimated of damage permeability. The neglect of injection tubes may also be a factor. Even so, the simulation provides supplementary information on the principal stress distribution, fracture-induced stress redistribution, pore-pressure distribution, and fracture aperture map that readily influenced by equipment error or

SC

cannot be observed directly during laboratory test. Based on the discussion, it is found that this study has presented high reliability in the simulation and contributed significant advancement to the

M AN U

numerical modeling of fracture propagation. The RFPA-Petrol simulator may have an enhanced capability of simulating more complex hydraulic-fracturing problems.

3.3 The analysis of multiple hydraulic fractures propagation from a single horizontal well The HMD model simulator was applied to investigate the propagation of multiple hydraulic fractures with different fracture spacing and different horizontal stress ratio. The basic material

TE D

properties and boundary conditions are similar to the model in part 3.2 (this paper). The initial hydraulic fractures are three circles whose radius are 20mm and parallel to the xy plane. The partial section of the numerical model built in the simulator is shown in Fig.16. In this model, the y-direction

EP

is same as the vertical stress direction, and the z-direction is the minimum horizontal stress direction. 3.3.1 Multiple hydraulic fractures with different fracture spacing

AC C

Seven numerical models were built to explore the impact of fracture spacing on the propagation of the middle fracture. The fracture spacing difference of neighboring models is 12 mm. The specific fracture spacing is listed in Table 5. Results show that the upper and lower fractures tend to propagate away from each other in every

model. Fig.17 shows maximum lengths of the continuous fractured zone in the x-direction and the y-direction of the middle fracture in numerical model 1~7. For simplicity, we call them as fracture length in the x-direction and fracture length in the y-direction. As other researchers’ results (Pan et al., 2014; Peirce and Bunger, 2015) showed, the length of the middle fracture after hydraulic fracturing increases with the fracture spacing. However, it is also shown that the direction of local minimum stress near the middle fracture changed as the two curves in the Fig.17 cross at (2.36,3.44). Before the intersection, the hydraulic fracture length in the y-direction is longer than that in the x-direction. 14 / 24

ACCEPTED MANUSCRIPT However, the fracture length in the y-direction declines sharply after fracture spacing is 40 mm because of the reorientation of local principal stresses near the middle fracture. At the same time, the length of fracture in the x-direction increases dramatically and exceeds the fracture length in the y-direction which means that the middle fracture propagates more easily in the x-direction than the y-direction after the intersection of two lines. These results cannot be gained from normal 2D simulations.

RI PT

Fig.18 shows the breakdown pressure (or peak pressure) of the middle fracture in model 1 to model 7. Initially, the peak pressure remains at a low level and declines between dz=16 mm to dz=40 mm. Then the breakdown pressure increases sharply as the direction of maximum principal stress turns. However, it has a decline from dz=64 mm to dz=88mm which means the fracture propagates more easily with the increase of fracture spacing. It is found that an outlier is present when dz/r equals to 2.6

SC

in both Fig.17 and Fig.18. It may be induced by the reorientation of local principal stresses near the middle fracture.

M AN U

Fig.19b illustrates that the middle fracture nearly propagates parallel to z-direction and merges into the other fractures when the fracture spacing is 52mm, i.e. dz/r=2.6. Two more minimum principal stress field figures of model 1 and model 7 are shown as Fig.19a and Fig.19c. They have the smallest and largest fracture spacing in this research. When the spacing is pretty small, the middle fracture will not crack, and the zone between the upper fracture and the lower fracture will be damaged slightly. When the spacing is large enough, three fractures propagate almost simultaneously. When fracture

TE D

spacing is larger than 16 mm and smaller than 88 mm, the middle fracture tends to coalesce into the other fractures and only two dominating hydraulic fracture propagate. 3.3.2 Multiple hydraulic fractures under different horizontal stress ratios

EP

Models with three parallel fractures whose spacings are 88 mm are simulated under different horizontal maximum stresses in x-direction which are 3 MPa, 4 MPa, 5 MPa, 6 MPa, 8 MPa, 10 MPa,

AC C

12 MPa, and 15 MPa. The corresponding horizontal differential stress and horizontal stress ratio are listed in Table 6.

It is found that the horizontal stress ratio has a significant influence on the propagation paths of

fractures. Fig.20 shows the length of the continuous fractured zone in the x-direction and the y-direction of the middle fracture. Both of the two curves are divided by three zones: the stress in the x-direction is the minimum in-situ stress, the stress in the x-direction is the middle in-situ stress and stress in the x-direction is the maximum in-situ stress. In the first zone, the fracture lengths in the y-direction are larger than that in the x-direction because the stress x is the minimum in-situ stress and the fracture propagate perpendicular to the minimum principal stress. In the second zone, the fracture length in the y-direction is still larger. Both of the two curves have an increasing tendency during this period as the stress in the z-direction is the minimum in-situ stress, and its direction is perpendicular to 15 / 24

ACCEPTED MANUSCRIPT the initial fracture planes. In the third zone, the stress in the x-direction exceeds the vertical stress and becomes the maximum in-situ stress, so the fracture length in the x-direction is longer than that in the y-direction. Meanwhile, it is found that the two curves decline sharply after σH/σh equal to 1.6. That is, when the maximum horizontal stress is larger than the vertical stress and minimum horizontal stress, the fracture length of the middle fracture decrease with the increase of horizontal stress ratio.

RI PT

The breakdown pressures of the middle fracture in different models are shown in Fig.21. In general, the breakdown pressure has an increasing tendency. However, it has slight fluctuations between σH/σh=1 and σH/σh=2. After σH/σh equal to 2, it grows linearly which means it is harder for the middle fracture to crack when the horizontal stress ratio is large.

Fig.22 shows the partial section of the propagation path inside the numerical models. Three

SC

typical models whose horizontal stress ratios are in the first zone, the second zone and the third zone separately are compared in this figure. In model 9, the deflection angle of upper fracture and lower

M AN U

fracture are larger than the other two models’. The middle fracture propagates along the initial fracture plane firstly due to the concentration of stress near the initial middle fracture. Then its propagation path turns to parallel to z-direction and merges into the other fractures at last. In model 11, three initial fractures propagate perpendicular to the vertical direction that is the direction of minimum principal stress and have small deflection angles because of the stress interference. In model 13, the deflection angles are the smallest compared to model 9 and model 11.

TE D

Numerical simulations of three parallel fractures propagation under hydraulic fracturing show that fracture spacing and horizontal stress ratio have a significant effect on fracture geometry and local stress field of the middle fracture. If all of the fractures are supposed to form dominating hydraulic fractures, the fracture spacing and horizontal stress ratio should be large enough.

Application of RFPA-Petrol on a field-scale problem

EP

4.

Simulations of horizontal well Fan154-P1 were established to verify RFPA-Petrol’s ability on

AC C

field scale problems. The horizontal well was drilled in block Fan154 in Shandong province of China by SHENGLI oilfield (Huang, 2012; Liu et al., 2012; Ren et al., 2011). This reservoir is mainly composed of fine-grained lithic arkose with 41.9% quartz and 35.1% feldspar in the lithic fragment, with 13-14% clay, 23-39% mixture of illite and smectite, 23-25% illite, 17-29% kaolinite, 8-31% chlorite in sandstone (Huang, 2012). The dimension of the reservoir model is 300m×100m×68m. Its mesh size, initial conditions, and boundary conditions are listed in Table 7. The reservoir profile and conceptual model is illustrated as Fig.23(Ren et al., 2011). It has four layers including one oil layer, one water layer, and two barriers. Their material parameters are shown in Table 8(Huang, 2012; Liu et al., 2012; Niu, 2017; Ren et al., 2011). In this simulation, the water layer and the oil layer share the same material parameters. In field fracturing operation, the horizontal well was parallel to the minimum horizontal stress and 16 / 24

ACCEPTED MANUSCRIPT perforations were vertical to it. It was drilled 4 meters below the top of oil layer to avoid fracturing the water layer. Open hole ball drop fracturing was chosen to stimulate this reservoir. This fracturing method induced a reduction of perforations’ initial length near the toe of the horizontal well. In this simulation, model-s1 and model-s2 were built with fracture spacing 100m and 30m. Initial fracture length and spacing used in the simulation is shown in Fig.24. The original length of longest

fractures F2, F2-1, F2-2, F2-3, and F2-4 are 2m. The others are 1m. Simulated

injection

pressure

and

flux

and

field

RI PT

perforations, i.e., half-height of initial fractures F1, F1-1 and F1-2 are 3m. Half-height of initial

measured

data

are

shown

in

Fig.25(Unconventional Project Department of Shengli Oilfield, 2013). In order to compare them clearly, we leave out the period of initial pressure test and move the curve of field measured pressure and flux

SC

40min in the negative direction of t. Simulated pressure data obtained from F1 in model-s1 and F1-1 in model-s2. Model-s1 and field fracturing operation have a similar fracture spacing, so we compared

M AN U

their curves firstly. The breakdown pressure of model-s1 is slightly higher than field measured pressure because the radius of perforation in the simulation is larger than that in field fracturing operation which results in the reduction of stress concentration. The pressure of model-s1 declined sharply after fracture initiating due to the mesh size is large and leads to the increase of leak-off fluid volume. The closure pressure of model-s1 is larger because the injection is not shut down after t=40min. Then, comparing the simulated pressure of model-s1 and model-s2, as heterogeneity is considered and physical element

model-s1 and model-s2.

TE D

parameters near F1 and F1-1 are different, there are errors between the breakdown/closure pressure of

The height and length of each fracture at t=80min are listed in Table 9. Data of model-s1 are discussed firstly. It shows that fracture F1 has the largest fractured zone in model-s1. That is, fracture

EP

with larger initial fracture height will form a longer and higher fractured zone. This conclusion is also suitable for model-s2. Average fracture length of F1-1 and F1-2 is larger than average fracture length of F2-1 to F2-4. Both of them are larger than that of F3-1 to F3-4.

AC C

Then fractures which have a same initial height in model-s1 and model-s2 are compared. It is

found that the height of F1-1 and F1-2 are almost identical with F1, while the length of F1-1 is smaller than F1. This may result from the high compressive zone between the boundary and F1-2. Staring from F1-1, every second fracture in model-s2 has a shorter length or height than fractures with the same initial height in model-s1, which indicates that this reservoir cannot be fractured very well if fracture spacing equal to 30m. Shape (propagation paths) of fractures are shown in Fig.26. It is found that fractures propagate vertically to x-direction which is the direction of minimum horizontal stress. Fracture propagation paths parallel to each other. Comparing Fig.26(b) with Fig.26(a), fractured zones induced by F2-1 to F2-4 in model-s2 are smaller than fracture zone induced by F2 in model-s1. Similarity, fractured zones 17 / 24

ACCEPTED MANUSCRIPT induced by F3-1, F3-2, and F3-4 in model-s2 are smaller than fracture zone induced by F3 in model-s1. Different colors in Fig.26 represent the magnitude of apertures. Most of the apertures are smaller than 0.0035m. Fig.27 and Fig.28 illustrate the minimum stress field and seepage field for t=80min. In Fig.27(a), compared with fracture F2 and F3, fracture F1 disturbs a wider stress field. It appears that perforations

RI PT

with larger initial length will form a wider stress disturbance region. In Fig.27(b), fracture spacing of model-s2 is reduced to 30m. Disturbance regions formed by two neighboring zones between fractures are under higher compressive stress. Therefore, the extension of some fractures, such as F2-1, F2-3, and F3-1, is delayed by the high compression field induced by their left and right fractures. Comparing Fig.27 to Fig.28, seepage disturbance regions are much smaller than stress disturbance region. The

SC

permeability of this reservoir is low and fracturing fluid only flow to several elements near the fracture. Therefore, the stress shadow effect is induced by fracture propagation rather than seepage.

M AN U

From Fig.27 and Fig.28, it is clearly shown that propagation paths of fractures in model-s2 alternately present as undulate, which is due to the inhomogeneity of the model and stress shadow effect. During hydraulic fracturing, fractures start propagating after the injection pressure reaching their breakdown pressure. Because of the inhomogeneity of the model, fracture propagation paths are asymmetry. When the left and right fractures of any fracture prefer propagating in the same direction, they will hinder the propagation of the middle fracture to this direction. So, propagating in the opposite

TE D

direction is easier for the middle fracture. As every second fracture influences the fracture between them, fracture paths and zones between them become waves. Fig.29 shows the elastic modulus of elements in model-s1 and model-s2 when t=80min. It is found that the barrier between the oil layer and the water layer stops the fracture very well and avoids

EP

fractures connecting the water layer. This certifies that the mass flow rate and position of perforations are reasonable. When fracture spacing is 100m, fracture spacing is wide enough to avoid influence by stress shadow effect. Besides of stress showdown effect, the choice of fracture spacing in field

AC C

fracturing operation is determined by many other factors, such as drainage radius and fracturing method. Finally, the spacing of 100m and 60m were chosen for fractures near the heal and toe of horizontal well separately.

These simulations have attempted to investigated to understand how fracture spacing affects

fracture propagation in a field-scale model. As the fracturing treatment method used in field fracturing operation is open hole ball drop method, different initial fracture height should be considered in the simulation. As expected, we found that fracture spacing and initial fracture height influence the volume of the final fractured zone. Fractures with larger initial height will form a longer fracture and a wider stress disturbance region. Fracture spacing should be large enough to ensure the propagation of fractures. This study has contributed advancement to choose a reasonable hydraulic fracture spacing in 18 / 24

ACCEPTED MANUSCRIPT field operations.

5.

Conclusions In

this

paper,

A

geotechnical

model

is

proposed

to

study

the

coupled

hydraulic-mechanical-damage processes in geological media. In this model, elastic damage mechanics describe the damage process of elements. Heterogeneity leads to non-linearity and causes progressive

RI PT

failure behavior observed in models. The maximum tensile stress/strain criterion and the Mohr-Coulomb criterion are utilized as the damage criterions. Priority is given to maximum tensile strain criterion in the determination of element state. Darcy’s Law and unsteady seepage flow equation describe fluid flow in elements, where compression of fluid and rock are considered. The apertures of

SC

hydraulic fractures are determined based on Sneddon and Elliott’s model (Sneddon and Elliot, 1946).

The coupled simulation process is fulfilled by Biot’s consolidation theory describing how pore-pressure affects stress field and the stress/aperture/mesh-dependent permeability which has three

M AN U

different formulations depending on the state of elements and represents the influence of stress to seepage field. The proposed HMD model is implemented in the calculation module of RFPA-Petrol. A reasonably good agreement is found between the analytically and numerically results for thick-wall cylinder model (Zhan et al., 2009).

An example of double hydraulic fractures propagation from two discal parallel perforations was conducted to demonstrate the capability of RFPA-Petrol in simulating the progressive damage of

TE D

laboratory-scale brittle rocks with hydraulic-mechanical processes involved. This example was selected because it was related to a common theme. The initiation pressures, fracture propagation paths, pore-pressure field and injection pressure variation tendency obtained by the simulation show good agreement with the laboratory test results reported in Rahman et al.’s studies (Rahman et al., 2002).

EP

Moreover, the simulation successfully reproduced the initiation, growth, and deflection of fractures in the model. Fracture propagation paths are asymmetrical due to the heterogeneity of rock and stress

AC C

shadow effects near the perforations. Moreover, the primary fracture’s curvature is smaller than the secondary fracture’s curvature. This asymmetrical effect has been successfully simulated by RFPA-Petrol, which has not been possible by analytical formulas. The RFPA-Petrol was used to investigate the three parallel hydraulic fractures propagation from a

horizontal well. The fracture spacing has a significant effect on the length of the continuous fractured zone and the direction of minimum principal stress. It is easier for the middle fracture to form a dominating hydraulic fracture with a larger fracture spacing. The horizontal stress ratio also influences the propagation of the middle fracture when the fracture spacing is large enough. The middle fracture will merge into other fractures if the ratio is low and will form a separate fracture when it is higher than 1. A horizontal well drilled by SHENGLI oilfield is simulated to discuss how fracture spacing and 19 / 24

ACCEPTED MANUSCRIPT initial fracture height influence fracturing effect. It is found that larger fracture spacing and initial fracture height induce a larger average volume of the final fractured zone. Simulation results can help to judge whether the fracture spacing is large enough to form efficient fractured zones. HMD-based RFPA-Petrol is a reliable technique for determining hydraulic fracture spacing in field operations. Despite the success in applying HMD based RFPA-Petrol to hydraulic fracturing problems, one of

RI PT

the drawbacks with the current RFPA-Petrol simulator is that the influence of anisotropy on the evolution of deformation and permeability is not considered. Moreover, fracture proppant is not taken into account. These limitations are the subjects of the forthcoming work. Although these issues require further clarification, the numerical results showed that RFPA-Petrol can simulate the entire failure process including fracture initiation, fracture growth, fracture deflection of fractured geological media

6.

SC

in a more natural manner.

Acknowledgments

M AN U

The study was jointly supported by the grants from National Science and Technology Major Project (Grant No. 2017ZX05072), National Natural Science Foundation of China (Grant Nos. 51761135102 and 51479024). The authors are grateful for these supports.

References

TE D

Adachi, A., Siebrits, E., Peirce, A. and Desroches, J., 2007. Computer Simulation of Hydraulic Fractures. International Journal of Rock Mechanics and Mining Sciences, 44(5): 739-757. http://doi.org/10.1016/j.ijrmms.2006.11.006. Belytschko, T. and Black, T., 1999. Elastic Crack Growth in Finite Elements with Minimal Remeshing. International Journal for Numerical Methods in Engineering, 45(5): 601-620. http://doi.org/10.1002/(SICI)1097-0207(19990620)45:5<601::AID-NME598>3.0.CO;2-S. Biot, M.A., 1941. General Theory of Three‐Dimensional Consolidation. Journal of Applied Physics,

AC C

EP

12(2): 155-164. Bryant, E.C., Hwang, J. and Sharma, M.M., 2015. Arbitrary Fracture Propagation in Heterogeneous Poroelastic Formations Using a Finite Volume-Based Cohesive Zone Model, SPE Hydraulic Fracturing Technology Conference. Society of Petroleum Engineers, The Woodlands, Texas, USA, pp. 17. Chen, C., Wan, J., Zhan, H. and Shen, Z., 2004. Physical and Numerical Simulation of Seepage-Pipe Coupling Model. Hydrogeology and Engineering Geology, 31(1): 1-8. Crawford, B.R. and Yale, D.P., 2002. Constitutive Modeling of Deformation and Permeability: Relationships between Critical State and Micromechanics, SPE/ISRM Rock Mechanics Conference. Society of Petroleum Engineers, Irving, Texas, pp. 10. Cundall, B.P.A., 1971. A computer model for simulating progressive large scale movements in blocky rock systems, Proc. Sympo. Int. Soc. Rock Mech. Damjanac, B. and Cundall, P., 2016. Application of Distinct Element Methods to Simulation of Hydraulic Fracturing in Naturally Fractured Reservoirs. Computers and Geotechnics, 71: 283-294. http://doi.org/10.1016/j.compgeo.2015.06.007.

20 / 24

ACCEPTED MANUSCRIPT

Freeze, R.A. and Cherry, J.A., 1980. Groundwater., 13(3): 455–458.

RI PT

Deng, S., Li, H., Ma, G., Huang, H. and Li, X., 2014. Simulation of Shale-Proppant Interaction in Hydraulic Fracturing by the Discrete Element Method. International Journal of Rock Mechanics and Mining Sciences, 70(9): 219-228. http://doi.org/10.1016/j.ijrmms.2014.04.011. Economides and Michael, J., 2007. Modern fracturing. ET Publishing. Fader, S.W., 1974. Ground water in the Kansas River valley, Junction City to Kansas City, Kansas. State Geological Survey, University of Kansas. Ferronato, M., Castelletto, N. and Gambolati, G., 2010. A Fully Coupled 3-D Mixed Finite Element Model of Biot Consolidation. Journal of Computational Physics, 229(12): 4813-4830. http://doi.org/10.1016/j.jcp.2010.03.018.

AC C

EP

TE D

M AN U

SC

Hazeghian, M. and Soroush, A., 2016. DEM-aided Study of Shear Band Formation in Dip-Slip Faulting through Granular Soils. Computers and Geotechnics, 71: 221-236. http://doi.org/10.1016/j.compgeo.2015.10.002. He, J., Liu, Q., Ma, G. and Zeng, B., 2016. An Improved Numerical Manifold Method Incorporating Hybrid Crack Element for Crack Propagation Simulation. International Journal of Fracture, 199(1): 21-38. http://doi.org/10.1007/s10704-016-0084-z. He, Q., Suorineni, F.T., Ma, T. and Oh, J., 2017. Effect of Discontinuity Stress Shadows On Hydraulic Fracture Re-Orientation. International Journal of Rock Mechanics and Mining Sciences, 91: 179-194. http://doi.org/10.1016/j.ijrmms.2016.11.021. Hu, M., Rutqvist, J. and Wang, Y., 2016. A Practical Model for Fluid Flow in Discrete-Fracture Porous Media by Using the Numerical Manifold Method. Advances in Water Resources, 97: 38-51. http://doi.org/10.1016/j.advwatres.2016.09.001. Huang, G., 2012. The Study On the Technologies for Large Scale Fracturing in G892 and F154 Blocks, China University of Petroleum (East China). Hubbert, M.K. and Willis, D.G.W., 1972. Mechanics of Hydraulic Fracturing. Mem. - Am. Assoc. Pet. Geol.; (United States), 18(6): 153-163. Khoei, A.R., Hosseini, N. and Mohammadnejad, T., 2016. Numerical Modeling of Two-Phase Fluid Flow in Deformable Fractured Porous Media Using the Extended Finite Element Method and an Equivalent Continuum Model. Advances in Water Resources, 94: 510-528. http://doi.org/10.1016/j.advwatres.2016.02.017. Klimczak, C., Schultz, R.A., Parashar, R. and Reeves, D.M., 2010. Cubic Law with Aperture-Length Correlation: Implications for Network Scale Fluid Flow. Hydrogeology Journal, 18(4): 851-862. Lecampion, B., 2009. An Extended Finite Element Method for Hydraulic Fracture Problems. Communications in Numerical Methods in Engineering, 25(2): 121-133. http://doi.org/10.1002/cnm.1111. Li, G., Allison, D. and Soliman, M.Y., 2011. Geomechanical Study of the Multistage Fracturing Process for Horizontal Wells, 45th U.S. Rock Mechanics / Geomechanics Symposium. American Rock Mechanics Association, San Francisco, California, pp. 6. Li, G., Tang, C., Li, L. and Li, H., 2016. An Unconditionally Stable Explicit and Precise Multiple Timescale Finite Element Modeling Scheme for the Fully Coupled Hydro-Mechanical Analysis of Saturated Poroelastic Media. Computers and Geotechnics, 71: 69-81. http://doi.org/10.1016/j.compgeo.2015.09.003. Li, L. et al., 2016. The Behaviour of Fracture Growth in Sedimentary Rocks: A Numerical Study Based on Hydraulic Fracturing Processes. Energies, 9(3): 169. http://doi.org/10.3390/en9030169. Li, L.C. et al., 2012. Numerical Simulation of 3D Hydraulic Fracturing Based on an Improved Flow-Stress-Damage Model and a Parallel FEM Technique. Rock Mechanics and Rock Engineering, 45(5SI): 801-818. http://doi.org/10.1007/s00603-012-0252-z.

21 / 24

ACCEPTED MANUSCRIPT

AC C

EP

TE D

M AN U

SC

RI PT

Li, L.C., Tang, C.A., Wang, S.Y. and Yu, J., 2013. A Coupled Thermo-Hydrologic-Mechanical Damage Model and Associated Application in a Stability Analysis On a Rock Pillar. Tunnelling and Underground Space Technology, 34(1): 38-53. http://doi.org/10.1016/j.tust.2012.10.003. Li, Y., Jiang, Y., Zhao, J., Liu, C. and Zhang, L., 2015. Extended Finite Element Method for Analysis of Multi-Scale Flow in Fractured Shale Gas Reservoirs. Environmental Earth Sciences, 73(10): 6035-6045. http://doi.org/10.1007/s12665-015-4367-x. Li, Z. et al., 2017. Numerical Investigation On the Propagation Behavior of Hydraulic Fractures in Shale Reservoir Based On the DIP Technique. Journal of Petroleum Science and Engineering, 154: 302-314. http://doi.org/10.1016/j.petrol.2017.04.034. Li, Z. et al., 2018. A Numerical Investigation On the Effects of Rock Brittleness On the Hydraulic Fractures in the Shale Reservoir. Journal of Natural Gas Science and Engineering, 50: 22-32. http://doi.org/10.1016/j.jngse.2017.09.013. Liang, Z.Z., Tang, C.A., Li, H.X., Xu, T. and Zhang, Y.B., 2004. Numerical Simulation of 3-D Failure Process in Heterogeneous Rocks. International Journal of Rock Mechanics and Mining Sciences, 41: 323-328. http://doi.org/10.1016/j.ijrmms.2004.03.061. Liu, Y.H., Chen, J.S., Gao, S. and Peng, S.U., 2012. Drilling Fluid Technology for the Slim-Hole Horizontal Well of Fan154-ping. Inner Mongolia Petrochemical Industry. Louis, C., 1975. Rock Hydraulics. International Journal of Rock Mechanics & Mining Science & Geomechanics Abstracts, 12(4): 59-59. MA, G., AN, X. and HE, L., 2011. THE NUMERICAL MANIFOLD METHOD: A REVIEW. International Journal of Computational Methods, 07(01): 1000204-. Ma, G.W., An, X.M., Zhang, H.H. and Li, L.X., 2009. Modeling Complex Crack Problems Using the Numerical Manifold Method. International Journal of Fracture, 156(1): 21-35. http://doi.org/10.1007/s10704-009-9342-7. Marina, S., Derek, I., Mohamed, P., Yong, S. and Imo-Imo, E.K., 2015. Simulation of the Hydraulic Fracturing Process of Fractured Rocks by the Discrete Element Method. Environmental Earth Sciences, 73(12): 8451-8469. http://doi.org/10.1007/s12665-014-4005-z. Niu, G., 2017. Simulation-Based Research On Prediction for Tight Sandstone Reservoir's Mechanical Parameters, Dalian University of Technology, Dalian. Oron, A.P. and Berkowitz, B., 1998. Flow in Rock Fractures: The Local Cubic Law Assumption Reexamined. Water Resources Research, 34(11): 2811-2825. http://doi.org/10.1029/98WR02285. Pan, L., Zhang, S., Cheng, L., Lu, Z. and Liu, K., 2014. A Numerical Simulation of the Inter-Cluster Interference in Multi-Cluster Staged Fracking for Horizontal Wells. Natural Gas Industry, 34(1): 74-9. http://doi.org/10.3787/j.issn.1000-0976.2014.01.011. Peirce, A. and Bunger, A., 2015. Interference Fracturing: Nonuniform Distributions of Perforation Clusters that Promote Simultaneous Growth of Multiple Hydraulic Fractures. Spe Journal, 20(2). http://doi.org/10.2118/172500-PA. Prochazka, P.P., 2004. Application of Discrete Element Methods to Fracture Mechanics of Rock Bursts. Engineering Fracture Mechanics, 71(4-6): 601-618. http://doi.org/10.1016/S0013-7944(03)00029-8. Rahman, M.M., Hossain, M.M., Crosby, D.G., Rahman, M.K. and Rahman, S.S., 2002. Analytical, Numerical and Experimental Investigations of Transverse Fracture Propagation From Horizontal Wells. Journal of Petroleum Science and Engineering, 35(PII S0920-4105(02)00236-X3-4): 127-150. http://doi.org/10.1016/S0920-4105(02)00236-X. Rankine, W.J.M., 2013. Manual of applied mechanics. C.griffen Co. Ren, Y.Z., Han, X.Z. and Zhang, Y.T., 2011. Slim Hole Horizontal Well Drilling Technology with Multi-Packer in Fan 154-P1 Well. Science Technology & Engineering.

22 / 24

ACCEPTED MANUSCRIPT Rios, A.M. et al., 2013. Stress Shadow Evaluations for Chicontepec–Evaluating New Completion

AC C

EP

TE D

M AN U

SC

RI PT

Options, 47th US Rock Mechanics/Geomechanics Symposium. American Rock Mechanics Association, San Francisco, California. Sepehri, J., Soliman, M.Y. and Morse, S.M., 2015. Application of Extended Finite Element Method to Simulate Hydraulic Fracture Propagation From OrientedPerforations, SPE Hydraulic Fracturing Technology Conference. Society of Petroleum Engineers, The Woodlands, Texas, USA, pp. 15. Sharafisafa, M. and Nazem, M., 2014. Application of the Distinct Element Method and the Extended Finite Element Method in Modelling Cracks and Coalescence in Brittle Materials. Computational Materials Science, 91(2): 102-121. http://doi.org/10.1016/j.commatsci.2014.04.006. Shi, G., 1992. Discontinuous Deformation Analysis: A New Numerical Model for the Statics and Dynamics of Deformable Block Structures. Engineering Computations, 9(2): 157-168. http://doi.org/10.1108/eb023855. Shimizu, H., Murata, S. and Ishida, T., 2011. The Distinct Element Analysis for Hydraulic Fracturing in Hard Rock Considering Fluid Viscosity and Particle Size Distribution. International Journal of Rock Mechanics and Mining Sciences, 48(5): 712-727. http://doi.org/10.1016/j.ijrmms.2011.04.013. Sisavath, S., Al-Yaaruby, A., Pain, C.C. and Zimmerman, R.W., 2003. A Simple Model for Deviations from the Cubic Law for a Fracture Undergoing Dilation or Closure. Pure and Applied Geophysics, 160(5): 1009-1022. http://doi.org/10.1007/PL00012558. Sneddon, I.N. and Elliot, H.A., 1946. The Opening of a Griffith Crack Under Internal Pressure. Quarterly of Applied Mathematics, 4(3): 262-267. Sneddon, I.N. and Elliott, H.A., 1946. The opening of a Griffith crack under internal pressure. Sukumar, N., Chopp, D.L. and Moran, B., 2003. Extended Finite Element Method and Fast Marching Method for Three-Dimensional Fatigue Crack Propagation. Engineering Fracture Mechanics, 70(PII S0013-7944(02)00032-21): 29-48. http://doi.org/10.1016/S0013-7944(02)00032-2. Tan, Y., Yang, D. and Sheng, Y., 2009. Discrete Element Method (DEM) Modeling of Fracture and Damage in the Machining Process of Polycrystalline SiC. Journal of the European Ceramic Society, 29(6): 1029-1037. http://doi.org/10.1016/j.jeurceramsoc.2008.07.060. Tang, C., 1997. Numerical Simulation of Progressive Rock Failure and Associated Seismicity. International Journal of Rock Mechanics & Mining Sciences, 34(2): 249-261. Tang, C.A., Tham, L.G., Lee, P., Yang, T.H. and Li, L.C., 2002. Coupled Analysis of Flow, Stress and Damage (FSD) in Rock Failure. International Journal of Rock Mechanics and Mining Sciences, 39(PII S1356-1609(02)00023-04): 477-489. http://doi.org/10.1016/S1365-1609(02)00023-0. Tang, C.A., Wang, W.T., Fu, Y.F. and Xu, X.H., 1998. A New Approach to Numerical Method of Modelling Geological Processes and Rock Engineering Problems - Continuum to Discontinuum and Linearity to Nonlinearity. Engineering Geology, 49(3-4): 207-214. http://doi.org/10.1016/S0013-7952(97)00051-3. Unconventional Project Department of Shengli Oilfield, 2013, Practice and Understanding of Tight Oil Unconventional Production Technique, https://wenku.baidu.com/view/7a22f1039b6648d7c0c7465f.html, 04/17/2018 Wang, H., 2016. Numerical Investigation of Fracture Spacing and Sequencing Effects On Multiple Hydraulic Fracture Interference and Coalescence in Brittle and Ductile Reservoir Rocks. Engineering Fracture Mechanics, 157: 107-124. http://doi.org/10.1016/j.engfracmech.2016.02.025. Wang, P., Mao, X., Lin, J. and Du, C., 2009. Study of the Borehole Hydraulic Fracturing and the Principle of Gas Seepage in the Coal Seam. Procedia Earth and Planetary Science, 1(1): 1561 1573. http://doi.org/10.1016/j.proeps.2009.09.241.

23 / 24

ACCEPTED MANUSCRIPT Witherspoon, P.A., Wang, J.S., Iwai, K. and Gale, J.E., 1980. Validity of Cubic Law for Fluid Flow in a Deformable Rock Fracture. Water Resources Research, 16(6): 1016–1024.

AC C

EP

TE D

M AN U

SC

RI PT

Wu, R., Kresse, O., Weng, X., Cohen, C. and Gu, H., 2012. Modeling of Interaction of Hydraulic Fractures in Complex Fracture Networks, SPE Hydraulic Fracturing Technology Conference. Society of Petroleum Engineers, The Woodlands, Texas, USA, pp. 14. Wu, Z., Fan, L., Liu, Q. and Ma, G., 2017. Micro-Mechanical Modeling of the Macro-Mechanical Response and Fracture Behavior of Rock Using the Numerical Manifold Method. Engineering Geology, 225(SI): 49-60. http://doi.org/10.1016/j.enggeo.2016.08.018. Yang, Y., Tang, X., Zheng, H., Liu, Q. and He, L., 2016. Three-Dimensional Fracture Propagation with Numerical Manifold Method. Engineering Analysis with Boundary Elements, 72: 65-77. http://doi.org/10.1016/j.enganabound.2016.08.008. Yuan, S.C. and Harrison, J.P., 2005. Development of a Hydro-Mechanical Local Degradation Approach and its Application to Modelling Fluid Flow During Progressive Fracturing of Heterogeneous Rocks. International Journal of Rock Mechanics and Mining Sciences, 42(7-8): 961-984. http://doi.org/10.1016/j.ijrmms.2005.05.005. Zeeb, C. and Konietzky, H., 2015. Simulating the Hydraulic Stimulation of Multiple Fractures in an Anisotropic Stress Field Applying the Discrete Element Method. Energy Procedia, 76: 264-72. http://doi.org/10.1016/j.egypro.2015.07.859. Zhan, M., Sheng, J. and Cen, J., 2009. Experimental and Analytical Study of Hydraulic Fracturing of Cylinder Sample. Advances in Water Resources and Hydraulic Engineering: 252-256. Zhang, G., Li, X. and Li, H., 2015. Simulation of Hydraulic Fracture Utilizing Numerical Manifold Method. Science China-Technological Sciences, 58(9SI): 1542-1557. http://doi.org/10.1007/s11431-015-5901-5. Zhang, Q., Lin, S., Xie, Z. and Su, H., 2016. Fractured Porous Medium Flow Analysis Using Numerical Manifold Method with Independent Covers. Journal of Hydrology, 542: 790-808. http://doi.org/10.1016/j.jhydrol.2016.09.054. Zheng, A. and Luo, X., 2015. Hydro-Mechanical Modeling of Impermeable Discontinuity in Rock by Extended Finite Element Method. Journal of Central South University, 22(11): 4337-4346. http://doi.org/10.1007/s11771-015-2982-z. Zheng, H., Liu, F. and Li, C., 2014. The MLS-based Numerical Manifold Method with Applications to Crack Analysis. International Journal of Fracture, 190(1-2): 147-166. http://doi.org/10.1007/s10704-014-9980-2. Zhou, L. and Hou, M.Z., 2013. A New Numerical 3D-model for Simulation of Hydraulic Fracturing in Consideration of Hydro-Mechanical Coupling Effects. International Journal of Rock Mechanics and Mining Sciences, 60(2): 370-380. http://doi.org/10.1016/j.ijrmms.2013.01.006. Zhou, W., Yan, G. and Yang, R., 1998. Elasto Brittle Damage Model for Rockmass Based On Field Tests in Laxiwa Arch Dam Site. Chinese Jounal of Geotechnical Engineering, 5. Zhu, W.C., Liu, J., Yang, T.H., Sheng, J.C. and Elsworth, D., 2006. Effects of Local Rock Heterogeneities on the Hydromechanics of Fractured Rocks Using a Digital-Image-Based Technique. International Journal of Rock Mechanics and Mining Sciences, 43(8): 1182-1199. http://doi.org/10.1016/j.ijrmms.2006.03.009.

24 / 24

ACCEPTED MANUSCRIPT Table 1 Mechanical parameters of the thick-wall cylinder model value

Poisson ratio

0.25

pore-pressure coefficient

1

homogeneity index

1000

AC C

EP

TE D

M AN U

SC

RI PT

parameter

ACCEPTED MANUSCRIPT Table 2 Seepage and boundary input parameters in RFPA-Petrol symbol

value

unit

volume

V

0.4*0.4*0.4

m3

the radius of the perforation region

r

0.02

m

permeability

K

8e-20

m2

Poisson ratio

ν

0.2

-

storativity

Ss

1.38e-6

m-1

flux

Q

1.10e-8

m3·s-1

fracture fluid density

ρ

897

kg·m-3

fracture fluid viscosity

µ

0.005

Pa·s

maximum in-situ horizontal stress

σx, σH

6

Mpa

vertical in-situ stress

σy, σv

7

minimum in-situ horizontal stress

σz, σh

SC 5

Parameters were selected to be similar to the laboratory test DC-3(Rahman et al., 2002)

AC C

EP

TE D

M AN U

1)

RI PT

parameter

Mpa

Mpa

ACCEPTED MANUSCRIPT Table 3 Macroscopic, microscopic and RFPA-Petrol input physical parameters symbol

value

unit

macroscopic elastic modulus

Em

7000

MPa

macroscopic strength

pm

35

MPa

homogeneity index

m

4

-

numerical input elastic modulus

Ei

8300

MPa

numerical input strength

Pi

130

Mpa

numerical output elastic modulus

Eo

6910

MPa

numerical output strength

Po

34

MPa

AC C

EP

TE D

M AN U

SC

RI PT

parameter

ACCEPTED MANUSCRIPT Table 4 Min principal stress and pore-pressure of elements at key times obtained in RFPA-Petrol t

min principal stress

hydraulic pressure

P

(s)

(MPa)

(MPa)

a

40

b

480

c

AC C

EP

TE D

M AN U

SC

RI PT

K

880

ACCEPTED MANUSCRIPT t

min principal stress

hydraulic pressure

P

(s)

(MPa)

(MPa)

d

1160

e

1360

f

AC C

EP

TE D

M AN U

SC

RI PT

K

1560

ACCEPTED MANUSCRIPT t

min principal stress

hydraulic pressure

P

(s)

(MPa)

(MPa)

g

2000

h

2480

AC C

EP

TE D

M AN U

SC

RI PT

K

ACCEPTED MANUSCRIPT

dz (mm)

dz/r

1

16

0.8

2

28

1.4

3

40

2

4

52

2.6

5

64

3.2

6

76

3.8

7

88

4.4

AC C

EP

TE D

M AN U

SC

model No.

RI PT

Table 5 Fracture spacings of numerical model 1~7

ACCEPTED MANUSCRIPT Table 6 The maximum horizontal stress conditions of numerical model 8~15 σH - σh, σx -σz

σH / σh, σx /σz

(MPa)

(MPa)

-

8

3

-2

0.6

9

4

-1

0.8

10

5

0

1

11

6

1

1.2

12

8

3

13

10

5

14

12

7

15

15

10

RI PT

σH, σx

1.6 2

2.4 3

AC C

EP

TE D

M AN U

SC

model No.

ACCEPTED MANUSCRIPT Table 7 Mesh size, initial conditions, and boundary conditions parameter

symbol

unit

value 1)

volume

V

m3

300×100×68

element number

nele

-

2040000

node number

nn

2097669 -1

m ·min

fracture fluid density

ρ

-3

kg·m

fracture fluid viscosity

µ

Pa·s

coupling coefficient

α

-

minimum in-situ horizontal stress

σx, σh

Mpa

maximum in-situ stress

σy, σH

Mpa

vertical in-situ horizontal stress

σz, σv

1)

4 1000

RI PT

Q

SC

flux

3

Mpa

0.001 1.8 39

45

62

Boundary conditions are gained from field-measured files (Ren et al., 2011; Huang, 2012; Liu et al., 2012;

AC C

EP

TE D

M AN U

Unconventional Project Department of Shengli Oilfield, 2013).

ACCEPTED MANUSCRIPT Table 8 Material parameters of different layers symbol

unit

oil layer1)

barrier layer2)

macroscopic elastic modulus

Em

GPa

12

16

uniaxial compressive strength

pm

MPa

60

90

homogeneity index

m

-

10

10

compression-tension ratio

-

-

10

10

friction angle

φ

°

30

30

Poisson ratio

ν

-

0.11

0.11

3e-15

1e-17

14.9

0.1

8.6e-7

6.0e-7

0.1

0.1

permeability

k

m

porosity

n

% -1

storativity

Ss

m

residual strength coefficient

λ

-

ultimate tensile strain coefficient

ηt

-

ultimate compressive strain coefficient

ηc

-

1)

5

5

200

200

Parameters of oil layers are gained from experiment tests of 54 core samples drilled from well Fan155-1(Huang, 2012; Liu et al.,

M AN U

2012; Ren et al., 2011).

Parameters of barrier layers are gained from experiment tests and numerical tests of core samples drilled from well Fan144(Niu,

EP

TE D

2017).

AC C

2)

SC

2

RI PT

parameter

ACCEPTED MANUSCRIPT Table 9 Height and length of fractures at t=80min coordinate x (m)

fracture height (m)

fracture length (m)

F1

50

18

21

F1-1

15

19

19

F1-2

45

18

25

F2

150

15

13

F2-1

75

15

11

F2-2

105

15

F2-3

135

10

F2-4

165

14

F3

250

2

F3-1

195

2

F3-2

225

2

F3-3

255

2

F3-4

285

2

RI PT

fracture No.

12 10 14 5

AC C

EP

TE D

M AN U

SC

1 3 3 5

ACCEPTED MANUSCRIPT

List of symbols Symbol

Unit

Property

a

m

mesh size

A

2

m

b

-

unit area of the fracture coupling coefficient -1

cM

MPa

df

m

equivalent diameter

d0

m

fracture aperture

dz

m

fracture spacing

D

-

damage variable

E

MPa

elastic modulus

Ee

MPa

elastic modulus of elastic element

Ed

MPa

modulus of damaged element

Ei

MPa

numerical input elastic modulus

Em

MPa

macroscopic elastic modulus

Eo

MPa

numerical output elastic modulus

f

MPa

body force

SC

M AN U

-1

RI PT

vertical uniaxial compressibility

f ( x)

MPa

g

m/s2

gravity

G

MPa

shear modulus

hf

m

fracture height

H

m

water head

Jf

-

hydraulic gradient of fluid

lf

m

m

-

n

-

p

MPa

pf

MPa

the pressure inside the fracture

MPa

numerical input strength

MPa

macroscopic strength

MPa

numerical output strength

P

MPa

fluid pressure

Pa

MPa

internal pore-pressure of the thick-wall cylinder model

Pb

MPa

external pore-pressure of the thick-wall cylinder model

Pnet

MPa

net pressure

q

3

m /s

volume sources

q0

3

m /s

unit discharge

r

m

distance from the center of the thick-wall cylinder model

ra

m

inside radius of the thick-wall cylinder model

rb

m

outside radius of the thick-wall cylinder model

pm

TE D

AC C

po

fracture length

homogeneity index porosity

pore-pressure

EP

pi

distribution density

-1

S

m

t

s

time

w

m

fracture aperture

storativity

ACCEPTED MANUSCRIPT Symbol

Unit

x

MPa

Property material parameters of the element (such as Young’s modulus and strength) obtained from experiments

x0

MPa

α

-

the mean value of elements put into the numerical simulation program coupling coefficient

β

MPa

ε

-

strain tensor

ε

-

equivalent principal strain

ε1

-

maximum principal strain

ε2

-

middle principal strain

ε3

-

minimum principal strain

ε c0

-

critical strain under peak compressive stress

ε cu

-

ultimate compression strain

ε t0

-

threshold strain

ε tu

-

ultimate tensile strain

η

-

κe

m

initial permeability

permeability of elastic element

permeability of damage element

m

2

κd

m

2

κf

m2

λ

MPa

Lamé constant

µ

Pa·s

fluid dynamic viscosity

ν

-

permeability of failure element

Poisson ratio

ρ

kg/m

σ

MPa

σ1

MPa

σ3

MPa

σc

MPa

σ c0

MPa

σ rc

TE D

fracturing fluid density stress tensor

maximum principal stress minimum principal stress fracture closure pressure uniaxial compressive strength

MPa

residual tensile stress

MPa

residual compressive stress

MPa

uniaxial tensile strength

AC C

σ t0

3

EP

σ rt

σ ii / 3

MPa

φ

-

SC

ultimate tensile stress factor 2

RI PT

fracturing fluid compressibility

M AN U

κ0

-1

average total stress friction angle

M AN U

SC

RI PT

ACCEPTED MANUSCRIPT

AC C

EP

TE D

Figure 1 Elastic damage constitutive law of an element under uniaxial stress

1 / 29

M AN U

SC

RI PT

ACCEPTED MANUSCRIPT

AC C

EP

TE D

Figure 2 Damage variable of an element (the x-axis is equivalent principal strain for tensile failure mode and maximum principal strain for shear failure mode)

2 / 29

RI PT

ACCEPTED MANUSCRIPT

AC C

EP

TE D

M AN U

SC

Figure 3 Schematics showing the elastic elements, damaged elements and failure elements in RFPA

3 / 29

RI PT

ACCEPTED MANUSCRIPT

a) an elastic element

b) a damaged element

AC C

EP

TE D

M AN U

SC

Figure 4 Schematic illustration of permeability of an element due to tensile failure

4 / 29

AC C

EP

TE D

M AN U

SC

RI PT

ACCEPTED MANUSCRIPT

Figure 5 The numerical analysis process and computation schematic of RFPA-Petrol

5 / 29

SC

RI PT

ACCEPTED MANUSCRIPT

AC C

EP

TE D

M AN U

Figure 6 Cross section of the thick-wall cylinder validation model

6 / 29

RI PT

ACCEPTED MANUSCRIPT

AC C

EP

TE D

M AN U

SC

Figure 7 Distribution of stresses in thick-wall cylinder numerical model

7 / 29

RI PT

ACCEPTED MANUSCRIPT

AC C

EP

TE D

M AN U

SC

Figure 8 Stress (radial and tangential stresses at θ= 0) comparison between numerical results and analytical results.

8 / 29

M AN U

SC

RI PT

ACCEPTED MANUSCRIPT

AC C

EP

TE D

Figure 9 Elastic modulus distribution in numerical uniaxial compression model

9 / 29

RI PT

ACCEPTED MANUSCRIPT

AC C

EP

TE D

M AN U

SC

Figure 10 Strain-stress curve of the numerical uniaxial compression model

10 / 29

(b) Fracture test configuration for Rahman et al.’s laboratory test

EP

TE D

M AN U

(a) Configuration for Rahman et al.’s laboratory test

SC

RI PT

ACCEPTED MANUSCRIPT

(c) Cross-section of the simulated model built in RFPA-Petrol

(d) Configuration for simplified simulation model

AC C

Figure 11 Graphical representation of the laboratory model and simulated model for a single horizontal well with two perforations (Fig.11(a) and Fig.11(b) are redrawn based on (Rahman et al., 2002) )

11 / 29

M AN U

SC

RI PT

ACCEPTED MANUSCRIPT

AC C

EP

TE D

(a) Injection pressure and rate record during laboratory multiple fracturing test, DC-3 (results from (Rahman et al., 2002))

(b) Injection pressure and rate variations during multiple fracturing simulation obtained from simulated analysis

Figure 12 Injection pressure response and rate variations during multiple fracturing (Fig.12(a) is redrawn based on (Rahman et al., 2002))

12 / 29

(b) t=1160s

(c) t=1360s

M AN U

SC

(a) t=880s

RI PT

ACCEPTED MANUSCRIPT

(d) t=1560s

(e) t=2000s

(f) t=2480s

AC C

EP

TE D

Figure 13 Fracture aperture of elements at different time obtained from RFPA-Petrol

13 / 29

M AN U

SC

RI PT

ACCEPTED MANUSCRIPT

AC C

EP

TE D

Figure 14 Propagation paths of multiple radial fractures obtained from laboratory test DC-3 (Rahman et al., 2002) and RFPA-Petrol analysis

14 / 29

RI PT

ACCEPTED MANUSCRIPT

(b) HMD simulation paths (this paper)

M AN U

SC

(a)Rahman et al.’s laboratory result

(c) Wang’s simulation paths

(d) Bryant et al.’s simulation paths

(e) Wu et al.’s simulation paths

AC C

EP

TE D

Figure 15 Comparison of numerical fracture propagation paths for two initially parallel fractures (Fig.15(a)(c)(d)(e) are redrawn based on (Rahman et al., 2002), (Wang, 2016), (Bryant et al., 2015), (Wu et al., 2012) separately.)

15 / 29

RI PT

ACCEPTED MANUSCRIPT

AC C

EP

TE D

M AN U

SC

Figure 16 The partial section of the numerical model with three parallel initial fractures

16 / 29

RI PT

ACCEPTED MANUSCRIPT

AC C

EP

TE D

M AN U

SC

Figure 17 Maximum lengths of the continuous fractured zone in the x and the y direction of the middle fracture in models with different fracture spacing at t=2500s

17 / 29

RI PT

ACCEPTED MANUSCRIPT

AC C

EP

TE D

M AN U

SC

Figure 18 The breakdown pressure (or peak pressure) of the middle fracture in models with different fracture spacing

18 / 29

(a) model1: dz=16mm

(b) model4: dz=52mm

RI PT

ACCEPTED MANUSCRIPT

(c) model7: dz=88mm

AC C

EP

TE D

M AN U

SC

Figure 19 Min principal stress fields of models with fracture spacing 16mm, 52mm, and 88mm

19 / 29

RI PT

ACCEPTED MANUSCRIPT

AC C

EP

TE D

M AN U

SC

Figure 20 Maximum lengths of the continuous fractured zone in the x and the y direction of the middle fracture in models with different horizontal stress ratios at t=2500s

20 / 29

RI PT

ACCEPTED MANUSCRIPT

AC C

EP

TE D

M AN U

SC

Figure 21 The breakdown pressure (or peak pressure) of the middle fracture in models with different horizontal stress ratios

21 / 29

(a) model 9: σH/σh = 0.8

(b) model 11: σH/σh = 1.2

RI PT

ACCEPTED MANUSCRIPT

(c) model 13: σH/σh = 2

AC C

EP

TE D

M AN U

SC

Figure 22 Fracture aperture of elements in different models obtained from numerical analysis

22 / 29

M AN U

SC

(a) reservoir profile of Fan154-P1

RI PT

ACCEPTED MANUSCRIPT

(b) conceptual model of part of Fan154-P1

AC C

EP

TE D

Figure 23 Reservoir profile and conceptual model (Fig.23(a) is redrawn based on (Ren et al., 2011))

23 / 29

TE D

M AN U

SC

(a) model-s1

RI PT

ACCEPTED MANUSCRIPT

(a) model-s2

AC C

EP

Figure 24 Elastic modulus map of models: (a) fracture spacing is 100m; (b) fracture spacing is 30m

24 / 29

M AN U

SC

RI PT

ACCEPTED MANUSCRIPT

AC C

EP

TE D

Figure 25 Injection pressure response and injection flux variations during fracturing ( field measured data obtained from (Unconventional Project Department of Shengli Oilfield, 2013))

25 / 29

M AN U

SC

(a) model-s1

RI PT

ACCEPTED MANUSCRIPT

(b) model-s2

AC C

EP

TE D

Figure 26 Fracture aperture map of models for t=80min

26 / 29

TE D

M AN U

SC

(a) model-s1

RI PT

ACCEPTED MANUSCRIPT

(b) model-s2

AC C

EP

Figure 27 Minimum principal stress field of models at t=80min

27 / 29

TE D

M AN U

(a) model-s1

SC

RI PT

ACCEPTED MANUSCRIPT

(b) model-s2

AC C

EP

Figure 28 Seepage field of models at t=80min

28 / 29

M AN U

SC

(a) model-s1

RI PT

ACCEPTED MANUSCRIPT

(b) model-s2

AC C

EP

TE D

Figure 29 Elastic modulus map of models at t=80min

29 / 29

ACCEPTED MANUSCRIPT Highlights Coupled hydraulic-mechanical-damage (HMD) geotechnical model for simulating the entire failure process in geological media Simulation of asynchronous initiation, asymmetric growth and non-equal deflection of

RI PT

fractures under hydraulic fracturing Investigation of how fracture spacing and horizontal stress ratio influence propagation and reorientation of three hydraulic fractures from a horizontal well

AC C

EP

TE D

M AN U

SC

Certification of HMD based simulator’s ability on field-scale models