Rock Mech Rock Eng DOI 10.1007/s00603-015-0816-9

ORIGINAL PAPER

Combined Finite-Discrete Element Method for Simulation of Hydraulic Fracturing Chengzeng Yan1 • Hong Zheng2 • Guanhua Sun2 • Xiurun Ge2

Received: 16 August 2014 / Accepted: 10 August 2015  Springer-Verlag Wien 2015

Abstract Hydraulic fracturing is widely used in the exploitation of unconventional gas (such as shale gas).Thus, the study of hydraulic fracturing is of particular importance for petroleum industry. The combined finitediscrete element method (FDEM) proposed by Munjiza is an innovative numerical technique to capture progressive damage and failure processes in rock. However, it cannot model the fracturing process of rock driven by hydraulic pressure. In this study, we present a coupled hydro-mechanical model based on FDEM for the simulation of hydraulic fracturing in complex fracture geometries, where an algorithm for updating hydraulic fracture network is proposed. The algorithm can carry out connectivity searches for arbitrarily complex fracture networks. Then, we develop a new combined finite-discrete element method numerical code (Y-flow) for the simulation of hydraulic fracturing. Finally, several verification examples are given, and the simulation results agree well with the analytical or experimental results, indicating that the newly developed numerical code can capture hydraulic fracturing process correctly and effectively. Keywords FDEM  Y-flow  Hydro-mechanical coupling  Fracture network connectivity  Seepage  Hydraulic fracturing

& Hong Zheng [email protected] 1

Key Laboratory of Urban Security and Disaster Engineering, Ministry of Education, Beijing University of Technology, Beijing 100124, China

2

State Key Laboratory of Geomechanics and Geotechnical Engineering, Institute of Rock and Soil Mechanics, Chinese Academy of Sciences, Wuhan 430071, China

Abbreviations M Lumped mass diagonal matrices C Damping diagonal matrices F Nodal force vector X Vector of nodal coordinates Dp Pressure differential between two adjacent nodes p1 Pressure of Node 1 p2 Pressure of Node 2 qw Mass density of fluid g Gravitational acceleration y1 Vertical coordinate of Node 1 y2 Vertical coordinate of Node 2 q Flow rate l Dynamic coefficient of viscosity a Average aperture of crack element a0 Initial value of aperture ares Minimum value of aperture amax Maximum value of aperture un Normal displacement of crack element L Length of crack element fs Function of saturation fn Total fluid pressure on edge of triangular element s Saturation of current time step s0 Saturation of previous time step Q Total flow rate q21 Flow rate from Node 1 to 2 q23 Flow rate from Node 3 to 2 q24 Flow rate from Node 4 to 2 p Pressure of the node at the current time step p0 Pressure of the node at the previous time step Kw Bulk modulus of the fluid ki Permeability factor of the ith crack element Dt Time step Dtf Critical time step V Volume of the node at the current time step

123

C. Yan et al.

V0 DV Vm V2 Vj21 Vj23 Vj24 h h1 h2 x B P P0 T t l f

Volume of the node at the previous time step Volume change Average volume Volume of Node 2 Volume of crack element j21 Volume of crack element j23 Volume of crack element j24 Height of free surface Water level on the left Water level on the right X-axis coordinate component Width of the dam Hydraulic pressure at the location with a distance to the impermeable boundary of (l - x) Hydraulic at the left end Non-dimensional time Time Length of single crack lx l

1 Introduction During unconventional gas exploitation (i.e., exploitation of shale gas), hydraulic fracturing is normally used to create fractures in the rock in order to extract gas resources. This particular process involves deformation and damage of the rock masses as well as the fluid flow in fractures, which is a typical hydro-mechanical coupling problem. However, the capability of existing numerical methods for dealing with these complex problems is very limited. Therefore, the study of hydraulic fracturing is of importance for petroleum industry. An international research project, DECOVALEX, has been conducted for many years, focussing on the problems of rock mechanics under multi-field coupling (Jing et al. 1995; Hudson et al. 2001; Tsang et al. 2005; Nowak et al. 2011). The project is divided into two parts: experimental study and numerical simulations. The numerical simulations are mainly based on continuum mechanics (Tsang et al. 2005), which is not suitable for heterogeneous and anisotropic materials. The conventional distinct element method (DEM) is also used in that project (Blum et al. 2005; Mas Ivars 2006). Although DEM is capable of modelling the joint fissures and the commercial software products such as UDEC and 3DEC can deal with the problems of multi-field coupling, it cannot simulate crack propagation directly. The particle flow methods (PFC2D and PFC3D) are very suitable for simulation of crack propagation, which provide a potential

123

way for tackling hydraulic fracturing (Al-Busaidi et al. 2005; Han and Cundall 2011, 2013; Shimizu et al. 2011). A theoretical model named FSD has also been commonly used to simulate hydraulic fracturing (Tang et al. 2002; Yang et al. 2004). For example, Wang et al. (2009) introduced a numerical model (F-RFPA2D) that can consider the coupling effect of seepage and deformation, to investigate the mechanism of crack initiation and propagation around a 2D cylindrical cavity in heterogeneous stiff soils during hydraulic fracturing. Li et al. (2012) applied a three-dimensional (3D) finite element model that simulates the coupling of seepage and deformation based on an improved flow-stress-damage model. In order to model the transition from continuous to discontinuous behaviour in fracture and fragmentation processes, Munjiza (1995, 2004) proposed the combined finite-discrete element method (FDEM) and developed the FDEM programme, known as Y-Code (Munjiza et al. 2011). The Y-code blends FEM techniques with DEM concepts, which is capable of simulating the complete process from the initiation and propagation of micro cracks to the macroscopic failure of rock. Later on, Mahabadi et al. (2010) developed a graphical user interface for Y-code. Furthermore, several algorithmic developments have been implemented to specifically address a broad range of rock mechanics problems (Mahabadi et al. 2012). Latham et al. (2013) investigated the influence of in situ stresses on the flow in fractured rock mass using FDEM and the Complex Systems Modelling Platform (CSMP??), where FDEM was merely used to model the deformation of the rock mass. Lisjak et al. (2014) illustrated the new developments of FDEM in modelling layered materials. Grasselli et al. (2014) presented the preliminary results obtained using FDEM to study the interaction between fluid driven fractures and rock mass discontinuities. Despite the great potential of Y-Code, hydro-mechanical coupling has not been directly incorporated in Y-Code. The purpose of this paper is to develop a hydromechanical coupling procedure for hydraulic fracture modelling based on FDEM, which is hereafter referred to as Y-flow. The paper is organised as follows. In Sect. 2, the fundamental principles of FDEM are briefly explained. In particular, the unique features of FDEM are highlighted. In Sect. 3, a two-dimensional coupled hydro-mechanical model is established for simulation of hydraulic fracturing with FDEM. In Sect. 4, an updating algorithm of hydraulic fracture network is proposed; In Sect. 5, we give the entire implementation procedure of hydro-mechanical coupling in FDEM. Finally, in Sect. 6, four verification examples are presented, with the comparisons with the analytical or experimental results.

Y-Flow: New Combined Finite-Discrete Element Numerical Code for Simulation of Hydraulic… Fig. 1 Intact solid representation in FDEM

Joint element

Intact solid

Triangular element

2 Fundamental Principles of FDEM In FDEM, each intact body is discreted with a mesh consisting of 3-node triangular elements (Lisjak and Grasselli 2014). In order to simulate the fracture of intact body, joint elements with initially no thickness are inserted at the common edges of the triangular elements, as shown in Fig. 1. The elastic deformation of the bulk material is captured by the constant-strain, linear-elastic triangular elements with impenetrability between elements enforced by a penalty function method. Fracture nucleation within the continuum is simulated by the breakage of the joint elements. Since fractures can nucleate along the boundaries of the triangular elements, the crack initiation location and propagation direction is traced in a more nature way compared with continuum methods (in which we have to deploy cracks in the model). The nodal coordinates are updated at each simulation time step, according to the governing equations for a FDEM system M€ x þ Cx_ ¼ FðxÞ

4 1 3 2

(a) newly generated fracture

4

1

3

ð1Þ 2

where M and C are the lumped mass and damping diagonal matrices of the system; x is the vector of nodal coordinates; and F is the nodal force vector, which includes contributions from the external loads, contact forces, deformation forces of triangular elements, and bonding forces of the joint elements. An explicit second-order finite-difference integration scheme is employed to solve the Eq. (1); the stiffness matrix does not need to be generated in the solving process. Therefore, the required memory is small even in simulating complex or large scale engineering problems. Unlike conventional DEM, the NBS algorithm and potential contact force are used to deal with contact detection and interaction in FDEM, respectively (Munjiza and Andrews 1998, 2000). The NBS contact detection algorithm is of linear complexity in time. The potential contact force avoids difficulties in tackling corner–corner

(b) Fig. 2 a Initial rock fracture network, b newly generated fracture connected with Crack no. 1

n Δa pold

Qold

aold

anew

Qnew

pnew

n

(a)

(b)

Fig. 3 Hydro-mechanical coupling in discontinuum

123

C. Yan et al. Fig. 4 Fluid flow in the hydraulic fracture network

Vj24/2 Vj21/2

contacts that exist in conventional DEM and the discontinuous deformation analysis, abbreviated by DDA (Zheng and Jiang 2009). The full explanation of FDEM can be found in Munjiza (2004); the NBS contact detection algorithm was introduced in Munjiza and Andrews (1998) and the potential contact forces in Munjiza and Andrews (2000).

a amax

a0

3 Hydro-Mechanical Coupling Model 3.1 Basic Assumptions The fractured joint elements are called crack elements in this study. The hydro-mechanical coupling model proposed in this study is based on the following hypotheses. (1) The triangular elements are impermeable and seepage flow occurs only in the crack elements. (2) Only those crack elements connecting with the hydraulic fracture network are involved in the seepage calculation. By a hydraulic fracture network we mean a collection of crack elements that connect with the hydraulic boundary. (3) The flow is assumed to be laminar, viscous, and incompressible, and the fluid pressure is always positive. That is to say, the fractal characteristics of the crack walls are ignored, because, otherwise, more complication will be involved, see Fan and Zheng (2013). As seen in Fig. 2a, the dashed line denotes the hydraulic boundary. Crack no. 1 is isolated and is not connected to the hydraulic boundary or the existing hydraulic fracture network. Thus, it is not a part of the hydraulic fracture network and will not been considered in the coming seepage calculation. However, when some newly generated fractures intersect with Crack no. 1 as seen in Fig. 2b, it becomes a part of the hydraulic fracture network and will participate in the coming seepage calculation.

123

Vj23/2

ares compressive

un

Fig. 5 Relationship between hydraulic aperture and the normal displacement of crack element

Therefore, we need an algorithm to dynamically determine which crack elements join a hydraulic fracture network and accordingly participate in seepage calculation. In Sect. 4, we will propose such an update algorithm to deal with this problem. 3.2 The Basic Idea of the Hydro-Mechanical Coupling The fluid in a crack applies fluid pressure on the two walls of the crack. The crack widens or narrows under the joint action of the fluid flow and the external force, which in turn affects the fluid flow in the crack. As shown in Fig. 3a, the fluid flow in the crack applies fluid pressure on the two walls of the crack. Under the fluid pressure, the aperture of the crack will increase (Fig. 3b). According to the cubic law, the unit width flow rate of fluid flowing in the crack is in cubic relationship with aperture of the crack, and the crack opening affects the flow rate of the fluid and the magnitude of the pressure.

Y-Flow: New Combined Finite-Discrete Element Numerical Code for Simulation of Hydraulic…

3.3 Fluid Flow Algorithm The fluid flow modelling is based on the hydraulic fracture network consisting of crack elements shown in Fig. 4, in which solid nodes denote nodes without fluid flowage and white nodes denote nodes with fluid flowage. The path connecting any two adjacent white nodes is a flow path as a part of the hydraulic fracture network, for example, the path between Nodes 4 and 2.

The procedure to solve the flow rate between two adjacent nodes is explained as follows. The pressure difference between the two nodes is first calculated. Take Node 2 as an example, the total pressure difference between Node 2 and Node 1 can be given by Dp ¼ p2  p1 þ qw gðy2  y1 Þ;

where p1 is the fluid pressure of Node 1, p2 is the fluid pressure of Node 2, and y1 and y2 are the vertical coordinates of Nodes 1 and 2, respectively. Using the cubic law, the flow rate from Nodes 1–2 can be calculated as: q¼

1

fn

p1

2 p2

fn

Fig. 6 Fulid pressure acting on triangular element boundary

ð2Þ

1 3 Dp a 12l L

ð3Þ

where l is the viscosity of the fluid, L is the length of the crack element. a is the average aperture of the crack element, its dependence on the average normal displacement of the crack element is illustrated in Fig. 5 (Itasca 2005). The introduction of the negative sign is to follow the convention: flow-out is negative and flow-in is positive. However, according to Eq. (3), when the pressures at both Node 1 and Node 2 are zero, the flow between the two nodes still occurs due to the effect of gravity. This is true only when both the nodes are fully saturated. The valid scenario is that the flow rate entering into Node 2 should decrease as the saturation of Node 1 reduces; when the saturation of Node 1 is zero, the flow rate entering from

Fig. 7 Schematic diagram of the algorithm for updating hydraulic fracturing network

(a)

(c)

(b)

(d)

123

C. Yan et al.

Node 1 into Node 2 is also zero. Therefore, the flow rates calculated by Eq. (3) are multiplied with a coefficient fs that is related to the saturation of Node 1, leading to q21 ¼ fs

1 3 Dp a ; 12l L

Node Forces caused by triangular and joint element deformation

ð4Þ Are new joint elements

where fs is a function of saturation, reading fs ¼ s2 ð3  2sÞ:

broken ?

ð5Þ

The empirical function has the following characteristics: fs = 0 if s = 0; fs = 1 if s = 1; and 0 B fs B 1 for any s[[0, 1]. In other words, the actual flow rate from Node 1 into Node 2 will decrease when the Node 1 is unsaturated. As Node 2 is connected with Nodes 3 and 4, the flow rate from Nodes 3 and 4 to Node 2, q23, q24 can be computed similarly. Therefore, the total flow rate of Node 2 can be calculated as: Q ¼ q21 þ q23 þ q24 :

Y Update hydraulic fracture network

Seepage analysis

ð6Þ

The fluid pressure of Node 2 can be calculated as (Itasca 2005) Dt DV p ¼ p0 þ Kw Q  Kw ; V Vm

Compute node forces caused by fluid pressure

ð7Þ Contact detection and node force

where p0 is the pressure of Node 2 at the previous time step, Kw is the bulk modulus of the fluid, Q is the total flow rate of Node 2, Dt is the time step, DV = V - V0, and

caused by contact force

0Þ Vm ¼ ðVþV 2 . Here, V and V0 are the volumes of Node 2 at current time step and previous time step, respectively. The volume of a node in the hydraulic fracture network is defined as half of the total volume of all crack elements connected to the node. Taking Node 2 in Fig. 4 as an instance, the volume of Node 2 is defined as

V2 ¼

Vj21 Vj23 Vj24 þ þ : 2 2 2

Dt DV  ; V Vm

Move to the next time step

Fig. 8 Flow chart of the calculation process in Y-Flow

ð9Þ

where s0 is the saturation of Node 2 at the previous time step. It should be indicated that if s \ 1, p is set zero and Eq. (9) is used instead of Eq. (7). If s [ 1, we set s = 1 and Eq. (7) is applied. To assure the numerical stability of the present explicit fluid flow algorithm, the time step needs to be limited to   V P Dtf ¼ min ; ð10Þ Kw i ki

123

Update node velocity and coordinates

ð8Þ

If p [ 0 is derived from Eq. (7), the change in volume of Node 2 leads to the change of pressure at the node. Otherwise, if p \ 0, the fluid in the volume of Node 2 is insufficient to fill that volume. In this case, we set p = 0 and update the saturation of Node 2 as s ¼ s0 þ Q

N

h

B=8m

h1=4m h2=1m x

Fig. 9 Dam model, hydraulic boundary and schematic diagram of free surface (Itasca 2005)

where V is the volume of the node, ki is the permeability factor of the ith crack element connected to the node. The minimum value of Dt over all nodes is used in the algorithm.

Y-Flow: New Combined Finite-Discrete Element Numerical Code for Simulation of Hydraulic… Fig. 10 Unstructured mesh for the dam in Fig. 9, where all inner triangular elements boundaries are specified joints

distributed loading. As shown in Fig. 6, assume p1 and p2 are the fluid pressures of node 1 and node 2, respectively. The total fluid pressure on this edge can be calculate as

Hydraulic pressure (Pa) 0

8000

16000

24000

32000

40000

fn ¼

ðp1 þ p2 Þ L: 2

ð11Þ

Then, the fluid pressure is converted to nodal forces that are applied at the nodes of the relevant triangular elements. Consequently, fluid pressure accumulated in each node is acting on triangular elements.

4 Algorithm for Updating Hydraulic Fracture Network

Fig. 11 Hydraulic pressure distribution in the dam

3.4 Fluid Pressure According to the fluid flow algorithm above, the fluid pressures of all nodes are calculated. The fluid pressure is acting on a triangular element boundary as a linearly

Fig. 12 Comparison of free surface between simulation result and analytical solution

The hydraulic fracture network essentially is a connected graph consisting of crack elements and their nodes. When a new crack element is created, the existing hydraulic fracture network needs to be updated. It should be noted that

4 Analytical Solution Numerical Solution

h(m)

3

2

1

0

1

2

3

4

5

6

7

8

x(m)

123

C. Yan et al.

the initial hydraulic network only consists of the nodes on the hydraulic boundary. The following algorithm is used when updating the hydraulic network. We denote N to be the set of the newly generated crack elements. The following operation is conducted. Select a crack element from N and mark it as J. Here, one of the nodes of J is required to belong to the existing hydraulic fracture network. We connect J and the other node into the hydraulic fracture network and then delete J from N. Repeat the above loop until the abovementioned J no longer exists in N. At this moment, the crack elements in N are isolated and cannot be connected to the existing hydraulic fracture network. In other words, they are not considered in the seepage calculation. Figure 7 demonstrates the abovementioned updating process, where the solid lines represent the crack elements in the existing hydraulic fracture network, and the dashed lines constitute the set N of newly generated crack elements.

The 5 newly generated crack elements are marked with dashed lines in Fig. 7a, constituting the set N. Wherein the difference between the newly generated crack elements 3 and 5 is that crack element 3 is connected to the existing hydraulic fracture network, and water will flow into it, and hence it will participate in the coming seepage calculation. On the contrary, crack element 5 is an isolated crack and not connected to the existing hydraulic fracture network. Thus, there is no water to enter crack element 5 and it will not participate in the coming seepage calculation.

J4

J6 J3

J1

J5

J2

J7

Fig. 15 Hydraulic fracturing model with complex pre-crack

P0 x

Table 1 Parameters for example 3 [Rock and Joint element parameters from Lisjak et al. (2013)]

l

Parameters

Value

Rock

Fig. 13 Rock band containing a single crack and mesh model

Bulk density, q(kg/m3)

2300

Young’s modulus, E(GPa)

3

Poisson’s ratio, v

0.29

Joint element 1

T=0.1 Numerical Solution T=0.1 Analytical Solution T=0.3 Numerical Solution T=0.3 Analytical Solution T=0.5 Numerical Solution T=0.5 Analytical Solution

0.8

Tensile strength, ft(MPa)

3

Internal cohesion, c(MPa) Friction angle of intact material, / ()

15 35

Friction angle of fractures, / ()

35

Mode I fracture energy release rate, GI(J/m2)

0.6

P/P0

2

0.4

2.0

Mode II fracture energy release rate, GII(J/m )

10

Normal contact penalty, pn(GPa)

30

Tangential contact penalty, pt(GPa/m)

3

Fracture penalty, pf(GPa)

15

Fluid flow parameters

0.2

0

0.2

0.4

0.6

0.8

x(m)

Fig. 14 Hydraulic pressure distributions in crack at different time

123

1

Bulk modulus of the fluid, Kw(GPa)

2.2

Viscosity of the fluid, l(Pa.s)

0.001

Density of water (kg/m3)

1000

Initial value of aperture, a0 (m)

5e–5

Minimum value of aperture, ares(m)

2.5e–5

Y-Flow: New Combined Finite-Discrete Element Numerical Code for Simulation of Hydraulic…

Hydraulic pressure (Pa) 2e5

1e6

1.8e6

2.6e6

4

3.4e6

J4

6 J6

3

J3 J1 1

step 1000

J2

J5

2

5

step 4000

J7

7

(a) 6

x 10

step 6000

step 8000

hydraulic pressure (Pa)

3.5 3

point 1 point 2

2.5

point 3

point 4 point 5 point 6 point 7

2 1.5 1 0.5 0

0.5

1

1.5

2

step

2.5

3 x 10 4

(b) Fig. 17 a Position of the feature points, b variation of hydraulic pressure with time at feature points step 10000

step 13000

Figure 7b–d show the sequence of the crack elements in set N to be connected to the hydraulic fracture network. Finally, set N only contains crack elements 2 and 5 as shown in Fig. 7d, which are isolated from the hydraulic fracture network.

step 17000

step 25000

step 19000

step 28000

Fig. 16 Crack propagation and hydraulic pressure distribution under hydraulic fracturing

5 Implementation Procedure of HydroMechanical Coupling in FDEM Based on the above procedure, we developed a new finitediscrete element method (FDEM) code—Y-flow, which implements the simulation of hydro-mechanical coupling in the original Y-code. As seen in Fig. 8, a typical time step in Y-Flow is stated as follows. Firstly, nodal forces are calculated based on triangular and joint element deformation at the current time step. Secondly, update the hydraulic fracture network when new crack elements are generated. Thirdly, the fluid pressure of nodes is

123

C. Yan et al.

computed and converted into linearly distributed tractions acting on triangular element boundaries. Fourthly, the nodal forces caused by contact force are calculated. Finally, the nodal coordinates and velocity are updated according to Eq. (1).

6 Illustrative Examples 6.1 Example 1: Validation of Pure Fracture Seepage with Free Surface Illustrated in Fig. 9 is a seepage flow through a homogeneous dam. The dam model is 8 m in width and 4 m in height. The water level on the left and right are h1 = 4 m and h2 = 1 m, respectively. The bottom is impermeable. According to Dupuit’s formula (Harr 1962), the free surface is given by rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi h2  h22 hð xÞ ¼ h21  1 x; ð12Þ B Fig. 18 A plaster specimen with pre-set cracks: a geometry of sample, b cross section of Mode I, c cross section of Mode II [Source: Jiao et al. (2015)]

123

with B = 8 m the width of the dam traverse section. Since the free surface given by Dupuit’s formula is independent of fluid parameters, we use this example to verify the correctness of the fluid flow algorithm regardless of the fluid parameter inputted in Y-flow. We replace the porous media in the dam by the joint media as shown in Fig. 10. The dam model is represented by the assembly of triangular elements bonded by joint elements to simulate homogenous isotropic material, where the impermeable material is composed of 5119 triangles with an average edge size of 0.16 m. As only the seepage flow analysis is validated, it is assumed that the crack aperture is constant a = 10-4 m, and the rock deformation is ignored. The hydraulic parameters include: the bulk modulus of water Kw = 2.2 GPa, the viscosity coefficient of water l = 0.001 Pa s, the density of water qw = 1000 kg/m3, and acceleration of gravity g = 10 m/s2. A constant integration time step of 7 9 10-6 s is used.

Y-Flow: New Combined Finite-Discrete Element Numerical Code for Simulation of Hydraulic…

As shown in Fig. 11, the hydraulic pressure is large at the left side and small at the right side, large at the bottom and small at the top. In addition, as shown in Fig. 12, the simulated free surface is very close to that evaluated by Dupuit’s formula (10). 6.2 Example 2: Validation of Pure Time-Dependent Seepage in a Single Crack Shown in Fig. 13 is a single crack in the rock sample with the dimension of 1 m 9 0.1 m. The right end is an impermeable boundary and a hydraulic pressure P0 is suddenly added to the left side. We study the variation of hydraulic pressure with time at different locations in the crack. The analytical solution of this example is (Carslaw and Jaeger 1959) 1 P 4X ¼ 1þ P0 p n¼0 " !#   ð2n þ 1Þp ð1Þnþ1 ð2nþ1Þ2 ðT=4Þp2 f   e cos ; 2 2n þ 1 ð13Þ where f ¼

lx l ,

T is non-dimensional time calculated as

T ¼ Kw

ða2 =ð12lÞÞt ; l2

ð14Þ

P is the hydraulic pressure at the location with a distance to the impermeable boundary of (l - x), and P0 is the hydraulic pressure at the left end. The hydraulic parameters include: the bulk modulus of water Kw = 2.2 GPa, the aperture of crack element a = 3910-5 m, the crack length l = 1 m, p0 = 9.5 MPa, and l = 1 9 10-3 Pa s. The mesh is shown in Fig. 13, where 80 identical isosceles right triangular elements are used. A constant integration time step of 2 9 10-9 s is taken. Hydraulic pressure distributions at different time are shown in Fig. 14, together with the analytical solution. As shown in Fig. 14, the hydraulic pressure increases gradually from left to right, the pressure at the same location increases with time. The simulated results agree well with the analytical solution. 6.3 Example 3: Validation of the Updating Algorithm of Hydraulic Fracture network In order to verify whether the updating algorithm of hydraulic fracture network is able to automatically search out all the crack elements in connection with the borehole,

Fig. 19 Experimental results of pre-set crack sample: a Mode I test, b Mode II test [Source: Jiao et al. (2015)]

123

C. Yan et al.

(a)

(b)

Fig. 20 Computational model and mesh of two specimen with pre-cracks under hydraulic fracturing: a mesh of Mode I

we designed an example as show in Fig. 15, where a square sample with a side length of 4 m, a central borehole with radius of 0.08 m. The six pre-set cracks are connected to the borehole boundary. In addition, there are also other preset fractures in the sample (a black line in the Fig. 15 represents a pre-set crack with an index of J1 to J7), which are far from the borehole. The model is meshed into 9916 triangular elements. Along the borehole boundary—the hydraulic boundary condition, a constant fluid pressure 3.5 MPa is loaded. A constant integration time step of 1 9 10-7 s is taken. Table 1 reports all mechanical and hydraulic parameters. Since the rock is impervious and water only flows in the crack elements, water in the borehole will first flow into the six pre-set cracks connected with the borehole boundary, then the pre-set cracks began to propagate under the hydraulic pressure, and may be intersected with other fractures far from the borehole. Afterward, water will flow into some of the fractures, and produce hydraulic pressure in these fractures. The entire fracturing process takes the following two features: (1) before the pre-set cracks extend, hydraulic pressure exists only in the pre-set cracks, which constitutes the initial hydraulic network; (2) as the hydraulic fracture network propagates and intersects with some cracks far from the borehole, hydraulic pressure in the cracks is observed; but no hydraulic pressure exists in the cracks not in connection with the hydraulic fracture network. We will examine whether the simulation results can reproduce the two features. The simulation result of crack propagation and hydraulic pressure distribution are shown in Fig. 16.

123

Table 2 Parameters for example 4 Parameters

Value

Rock Bulk density, q(kg/m3)

1.1e3

Young’s modulus, E(GPa)

4.0

Poisson’s ratio, v

0.3

Joint element Tensile strength, ft(MPa)

0.3

Internal cohesion, c(MPa)

1.0

Friction angle of intact material, / ()

32

Friction angle of fractures, / ()

32

Mode I fracture energy release rate, GI(J/m2)

1.0

Mode II fracture energy release rate, GII(J/m2) Normal contact penalty, pn(GPa)

10 4.0

Tangential contact penalty, pt(GPa/m)

0.4

Fracture penalty, pf(GPa)

4.0

Fluid flow parameters Bulk modulus of the fluid, Kw(GPa)

2.2

viscosity of the fluid, l (Pa.s)

0.001

Density of water (kg/m3)

1000

Initial value of aperture, a0 (m)

5e–5

Minimum value of aperture, ares(m)

2.5e–5

As shown in Fig. 16, water first enters the six pre-set cracks till step 1000. At step 4000, the six pre-set cracks have propagated to some distance but have not been in connection with other pre-set cracks far from the borehole. And no hydraulic pressure distribution is observed in these cracks, reproducing the first feature. As the hydraulic fracture network, namely the

Y-Flow: New Combined Finite-Discrete Element Numerical Code for Simulation of Hydraulic… Fig. 21 Displacement contours of Mode I

step 1000

step 10000

step 30000

step 40000

step 5000

step 20000

step 35000

step 46000

123

C. Yan et al.

step 5000

step 40000

step 30000

step 46000

Fig. 22 Velocity vector diagram of Model I

network of those cracks in connection with the borehole, propagates further, it begins to intersect at step 6000 with some cracks far from the borehole and the hydraulic pressure emerges inside these cracks at step 8000 and later steps. This reproduces the second feature. Figure 17b displays the change in hydraulic pressure of some feature points marked in Fig. 17a. Initially, the hydraulic pressure of all feature points is zero. Then the hydraulic pressure of point 1 first increases. When the pre-set crack J2 extends to point 2, hydraulic pressure in the feature point 2 is observed, which corresponds to step 4000 in Fig. 16. The pre-set crack J2 continues to grow and eventually intersects with the isolated crack J3. Because J3 and J4 are connected, point 3, located on crack J4, begins to have water pressure, as shown in Fig. 16 at step 6000 and step 8000. Then the upper end of crack J4 begins to extend under the action of hydraulic pressure. When crack J4 extends to point 4, point 4

123

begins to have water pressure as shown in Fig. 16 at step 19000. With the further propagation of the hydraulic fracture network, the newly generated cracks are connected to the isolated fracture J5, and point 5 on J5 begins to have hydraulic pressure, as shown in Fig. 16 at step 25000. However, the pressure of points 6 and 7 are always zero because cracks J6 and J7, on which points 6 and 7 are located, respectively, have never connected to the hydraulic fracture network, as shown in Fig. 16 at step 28000. According to the above simulation results, the algorithm for updating hydraulic fracture networks captures exactly the main features in the fracturing process. 6.4 Example 4: validation of hydraulic fracturing of a plaster sample with pre-set cracks Jiao et al. (2015) designed an experiment to verify the validity of his code DDARF to simulate hydraulic fracturing. A cube plaster sample with a hole in the centre,

Y-Flow: New Combined Finite-Discrete Element Numerical Code for Simulation of Hydraulic… Fig. 23 The maximum principal stress contours of Mode I (positive denote tension)

step 1000

step 5000

step 10000

step 20000

step 30000

step 35000

step 40000

step 46000

123

C. Yan et al. Fig. 24 Hydraulic pressure distribution of model I

123

step 1000

step 5000

step 10000

step 20000

step 30000

step 35000

step 40000

step 46000

Y-Flow: New Combined Finite-Discrete Element Numerical Code for Simulation of Hydraulic… Fig. 25 Displacement contours of Mode II

step 1000

step 10000

step 30000

step 40000

step 5000

step 20000

step 35000

step 46000

123

C. Yan et al. Fig. 26 Velocity vector diagram of Model II

step 5000

step 30000

step 40000

step 46000

which is fractured by injecting high-pressure water into the borehole, and no other loads are applied. The geometry of the sample is shown in Fig. 18a. They made two groups of samples with pre-set cracks, with inclination angles of 0 (Model I) and 30 (Mode II), respectively. As shown in Fig. 18b, c, the cross-section is perpendicular to the borehole, two cracks are symmetrically set along the diameter direction of the borehole. The length and depth of the preset crack is 20 and 90 mm, respectively. Figure 19 shows the failure patterns of the two groups of tested samples. As can be seen, under the action of highpressure water, the sample extends strictly from the crack tips and along the direction of original pre-set cracks, until the complete failure. Strictly speaking, this example does not fall into the categories of plane problems. Since both the sample and the borehole are quite thick, however, the plane strain assumption should not bring about a too much deviation. As a result, we use this experiment to verify the effectiveness of Y-flow to simulate hydraulic fracturing by building two numerical models: Model I (the inclination angle of crack is 0) and Model II (the inclination angle of crack is 30). As shown in Fig. 20a, these two models are meshed into 10148 and 9540 triangular elements,

123

respectively. Now we are increasing water pressure in the borehole with 5 Pa in each step. A constant integration time step of 1 9 10-7 s is taken. Mechanical parameters and hydraulic parameters are listed in Table 2. The macromechanical parameters are from the literature (Jiao et al. 2015) and the micro-mechanical parameters are calibrated according to the literature (Tatone and Grasselli 2015). Hydraulic fracturing simulation results of these two models are shown in Figs. 21, 22, 23, 24, 25, 26, 27, 28. For the simulation result of Model I, Fig. 21 is the displacement contour of fracturing process. Under the action of high-pressure water, the two pre-set cracks initiate from their tips and then propagate, as shown in Fig. 21 at step 35000. The crack path propagates along the direction of pre-set crack strictly; the specimen is eventually cut into two symmetrical blocks, as shown in Fig. 21 at step 46000. Figure 22 is a velocity vector diagram, from which we can see that the sample was separated in the opposite direction under the effect of high-pressure water. Figure 23 is the evolution of maximum principal stress (tension is positive). Initially, step 1000 say, the tensile stress concentration only appears on the left and right of borehole. Afterwards, step 5000 say, water enters the pre-set cracks, tensile stresses are concentrated on the pre-crack tips.

Y-Flow: New Combined Finite-Discrete Element Numerical Code for Simulation of Hydraulic… Fig. 27 The maximum principal stress contours of Mode II (positive denote tension)

step 1000

step 5000

step 10000

step 20000

step 30000

step 40000

step 35000

step 46000

123

C. Yan et al. Fig. 28 Hydraulic pressure distribution of model II

step 1000

step 5000

step 10000

step 20000

step 30000

step 40000

123

step 35000

step 46000

Y-Flow: New Combined Finite-Discrete Element Numerical Code for Simulation of Hydraulic…

Subsequently, step 35000 say, the pre-set cracks begin to extend, and the tensile stress concentration zones moves with the crack tip. This suggests that the crack propagation of the sample is controlled by the tensile stress. Figure 24 is hydraulic pressure distribution with hydraulic fracturing. Water enters the pre-set crack firstly (step 1000, say); then the pre-set cracks are filled with water (step 5000). With the increase of hydraulic pressure, the pre-set cracks begin to extend (step 35000), and water flows into the newly generated cracks. The sample is cut into two pieces finally. The farther the feature points are from the borehole, the smaller hydraulic pressure of feature points will be (step 46000 in Fig. 24). The fracturing process of Model II is similar as shown in Figs. 25, 26, 27, 28, the propagation is also along the direction of the pre-set cracks, and eventually the sample is cutting into two symmetrical blocks. Through comparison of Figs. 21 and 19a, Figs. 25 and 19b, we can see the simulation results agree well with the experimental results.

7 Conclusions Through the numerical tests including those in this study, the developed code Y-flow is able to simulate the process of hydraulic fracture as well as to conduct the hydraulic analysis of fissured rock networks. This greatly extends the application scope of FDEM. Meanwhile, the proposed algorithm for updating the hydraulic fracture network can accurately update most complicated hydraulic fracture network during the process of hydraulic fracturing, which can be easily extended to the 3-dimensional situations. Since cracking only extends along the solid element boundaries, the simulation results are unavoidably influenced by the size and mesh configuration, although such influences would diminish with the mesh refinement. Acknowledgments This work has been supported by the National Natural Science Foundation of China under the Grant number 11202223; and the Natural Basic Research Program of China under the Grant numbers: 2011CB013505 and 2014CB047100. Special thanks to the anonymous reviewer and editor-in-chief who have given their time and expertise to improve this article.

References Al-Busaidi A, Hazzard JF, Young RP (2005) Distinct element modeling of hydraulically fractured Lac du Bonnet granite. J Geoph Res-Solid Earth 110:BO6302 Blum P, Mackay R, Riley MS, Knight JL (2005) Performance assessment of a nuclear waste repository: upscaling coupled

hydro-mechanical properties for far-field transport analysis. Int J Rock Mech Min Sci 42(5):781–792 Carslaw HS, Jaeger JC (1959) Conduction of heat in solids, 2nd edn. Clarendon Press, Oxford Fan H, Zheng H (2013) MRT-LBM-based numerical simulation of seepage flow through fractal fracture networks. Sci China Technol Sci 52(12):3115–3122 Grasselli G, Lisjak A, Mahabadi OK, Tatone BSA (2014) Influence of pre-existing discontinuities and bedding planes on hydraulic fracturing initiation. Euro J Environ Civil Eng 19(5):580–597 Han Y, Cundall PA (2011) Lattice Boltzmann modeling of pore-scale fluid flow through idealized porous media. Int J Numer Meth Fluids 67(11):1720–1734 Han Y, Cundall PA (2013) LBM–DEM modeling of fluid–solid interaction in porous media. Int J Numer Anal Meth Geomech 37(10):1391–1407 Harr ME (1962) Groundwater and seepage. McGraw-Hill Book Co., Inc., New York Hudson JA, Stephansson O, Andersson J, Tsang CF, Jing L (2001) Coupled T–H–M issues relating to radioactive waste repository design and performance. Int J Rock Mech Min Sci 38(1):143–161 Itasca (2005) UDEC version 4.0 user’s manuals, Itasca Consulting Group Inc, Minnesota, USA Jiao Y-Y, Zhang H-Q, Zhang X-L, Li H-B, Jiang Q-H (2015) A twodimensional coupled hydromechanical discontinuum model for simulating rock hydraulic fracturing. Int J Numer Anal Meth Geomech 39(5):457–481 Jing L, Tsang CF, Stephansson O (1995) DECOVALEX—An international co-operative research project on mathematical models of coupled THM processes for safety analysis of radioactive waste repositories. Int J Rock Mech Mining Sci Geomech Abst 32(5):389–398 Latham J-P, Xiang J, Belayneh M, Nick HM, Tsang C-F, Blunt MJ (2013) Modelling stress-dependent permeability in fractured rock including effects of propagating and bending fractures. Int J Rock Mech Min Sci 57:100–112 Li LC, Tang CA, Li G, Wang SY, Liang ZZ, Zhang YB (2012) Numerical simulation of 3d hydraulic fracturing based on an improved flow-stress-damage model and a parallel FEM technique. Rock Mech Rock Eng 45(5):801–818 Lisjak A, Grasselli G (2014) A review of discrete modeling techniques for fracturing processes in discontinuous rock masses. J Rock Mech Geotech Eng 6(4):301–314 Lisjak A, Liu Q, Zhao Q, Mahabadi OK, Grasselli G (2013) Numerical simulation of acoustic emission in brittle rocks by two-dimensional finite-discrete element analysis. Geophys J Int 195(1):423–443 Lisjak A, Tatone BA, Grasselli G, Vietor T (2014) Numerical modelling of the anisotropic mechanical behaviour of opalinus clay at the laboratory-scale using FEM/DEM. Rock Mech Rock Eng 47(1):187–206 Mahabadi OK, Grasselli G, Munjiza A (2010) Y-GUI: a graphical user interface and pre-processor for the combined finite-discrete element code, Y2D, incorporating material heterogeneity. Comput Geosci 36(2):241–252 Mahabadi O, Lisjak A, Munjiza A, Grasselli G (2012) Y-Geo: new combined finite-discrete element numerical code for geomechanical applications. Int J Geomech 12(6):676–688 Mas Ivars D (2006) Water inflow into excavations in fractured rock— a three-dimensional hydro-mechanical numerical study. Int J Rock Mech Min Sci 43(5):705–725 Munjiza A (2004) The combined finite-discrete element method. Wiley, Chichester Munjiza A, Andrews KRF (1998) NBS contact detection algorithm for bodies of similar size. Int J Numer Meth Eng 43(1):131–149

123

C. Yan et al. Munjiza A, Andrews KRF (2000) Penalty function method for combined finite-discrete element systems comprising large number of separate bodies. Int J Numer Meth Eng 49(11):1377–1396 Munjiza A, Owen DRJ, Bicanic N (1995) A combined finite-discrete element method in transient dynamics of fracturing solids. Eng Comput 12(2):145–174 Munjiza A, Knight E, Rougier E (2011) Computational mechanics of discontinua. Wiley Nowak T, Kunz H, Dixon D, Wang W, Go¨rke U-J, Kolditz O (2011) Coupled 3-D thermo-hydro-mechanical analysis of geotechnological in situ tests. Int J Rock Mech Min Sci 48(1):1–15 Shimizu H, Murata S, Ishida T (2011) The distinct element analysis for hydraulic fracturing in hard rock considering fluid viscosity and particle size distribution. Int J Rock Mech Min Sci 48(5):712–727 Tang CA, Tham LG, Lee PKK, Yang TH, Li LC (2002) Coupled analysis of flow, stress and damage (FSD) in rock failure. Int J Rock Mech Min Sci 39(4):477–489

123

Tatone BSA, Grasselli G (2015) A calibration procedure for twodimensional laboratory-scale hybrid finite–discrete element simulations. Int J Rock Mech Min Sci 75:56–72 Tsang C-F, Jing L, Stephansson O, Kautsky F (2005) The DECOVALEX III project: a summary of activities and lessons learned. Int J Rock Mech Min Sci 42(5):593–610 Wang SY, Sun L, Au ASK, Yang TH, Tang CA (2009) 2D-numerical analysis of hydraulic fracturing in heterogeneous geo-materials. Constr Build Mater 23(6):2196–2206 Yang TH, Tham LG, Tang CA, Liang ZZ, Tsui Y (2004) Influence of heterogeneity of mechanical properties on hydraulic fracturing in permeable rocks. Rock Mech Rock Eng 37(4):251–275 Zheng H, Jiang W (2009) Discontinuous deformation analysis based on complementary theory. Sci China Ser E: Technol Sci 52(9):2547–2554