Geomech. Geophys. Geo-energ. Geo-resour. (2016) 2:343–363 DOI 10.1007/s40948-016-0040-4

ORIGINAL ARTICLE

Numerical simulation of fluid flow through deformable natural fracture network Peijie Yin . Gao-Feng Zhao

Received: 8 September 2015 / Accepted: 1 October 2016 / Published online: 15 October 2016 Ó Springer International Publishing Switzerland 2016

Abstract In this paper, fluid flow through natural fracture network is studied using computational fluid dynamics. To investigate the influence of fracture roughness, normal deformation and shear deformation on the fracture transmissivity/permeability, numerical tests of fluid flow through 3D rock fracture are conducted using the lattice Boltzmann method (LBM) in a middle size cluster. An empirical equation was obtained from the numerical results. Following this, natural fracture networks are built for fluid dynamics simulation of fluid flow through rock fracture network. It is found that the pipe network model enriched with the derived empirical equation can produce similar results compared with the LBM simulation which further confirms the empirical equation’s applicability. Finally, influences of fracture length, fracture density, and deformation of the fracture network on the fluid flow are studied preliminarily from coupling LBM with the discrete fracture network model and discrete element model.

P. Yin School of Highway, Chang’an University, Xi’an 710064, China G.-F. Zhao (&) State Key Laboratory of Hydraulic Engineering Simulation and Safety, School of Civil Engineering, Tianjin University, Tianjin 300072, China e-mail: [email protected]

Keywords Fluid flow  Fracture roughness  Deformation  Fracture network  Pipe network model  Lattice Boltzmann method

1 Introduction The flow behavior in fracture networks has been a research focus over the past half century. The discrete fracture network model (DFN) model has become the most widely used method since the work by Long et al. (1982). A DFN typically combines deterministic and stochastic discrete fractures, which presents the same geological statistics properties as observations, such as fracture density, distribution of location, orientation, size and hydraulic aperture. Extensive work was conducted to investigate the fluid flow behaviors in rock fractures using DFN model (e.g. Long and Witherspoon 1985; Chen et al. 1999; Dershowitz et al. 2004; Kim et al. 2007; Parker 2007; Liu et al. 2014). However, in DFN, the fluid flow is calculated based on the cubic law under the assumption that the fractures are plate surfaces. In practice, the fractures are rough with variety of profiles and aperture distribution which can be changed dynamically under normal and shear deformation. An accurate prediction of hydraulic behavior in fracture network requires a clear understanding of fluid flow though single fracture under these coupled conditions.

123

344

Geomech. Geophys. Geo-energ. Geo-resour. (2016) 2:343–363

hydraulic property of fracture to the stress rather than deformation of the fracture. The stress-permeability relationship is complex, which is influenced by many factors, such as stress condition and fracture profiles as well as the deformation. The mechanism of the flow behavior behind these experiments is not clearly understood because the geometry within the fracture is not easy to be controlled and obtained. Paralleling with the experimental study, extensive theoretical analysis was conducted (e.g. Zimmerman and Bodvarsson 1996). Theoretically, the flow of incompressible Newtonian viscous fluid is governed by the Navier–Stokes equation (Batchelor 1967). However, the Navier–Stokes equation cannot be solved in closed form when dealing with realistic fracture with rough surfaces. Alternatively, numerical approaches provided the opportunity to obtain the solution of fluid flow though rough surfaces under complex boundary conditions. There are varieties of traditional methods developed for fluid simulations, which are based on discretized partial differential equations, such as finite differences

Comprehensive work has been conducted to investigate the flow behavior in single fracture including experimental investigation, theoretical analysis and numerical simulation. The early work on fluid flow in single fracture was conducted experimentally by Lomize (1951). The cubic law was found essentially valid for laminar flow in rock joints based on the assumption of parallel flat surface. However, fracture walls contain irregularities which reduce fluid flow and lead to a local channeling effect of preferential flow. A large number of laboratory studies were carried out, the validation of cubic law was discussed and different empirical correction of the cubic law were proposed (e.g. Iwai 1976; Witherspoon et al. 1980; Neuzil and Tracy 1981; Tsang 1984; Barton et al. 1985; Brown 1987; Barton and Quadros 1997). Recently, influence of deformation on fluid flow in single fracture receives more attention. For example, Koyama et al. (2008) conducted the coupled shearflow tests for rock fractures. Indraratna et al. (2014) investigate the fluid flow through deformable rough rock joints. However, most of the works relate the Fig. 1 The D2Q9 model and D3Q15 model. a D3Q15 model, b D2Q9 model

5 10(14)

z

2 11(7)

5

6

y

y 2

x

1

3

1

3D

x 2D

8(12)

13(9)

7

8

6

4

(a)

(b)

Fig. 2 Bounce back scheme

f5

f2

f2

f7

f5

f6

Streaming f3

f1

f3

f7

f8

f1

f6

f4 f8 f4 wall

123

wall

Geomech. Geophys. Geo-energ. Geo-resour. (2016) 2:343–363

345 1.2

Values

Density (q)

1.0

Reynolds number (Re)

1.0

Resolution (N, l.u.)

10–100

dx

1/N

dt

*dx2

(e.g. Ames 1977; Morton and Mayers 1994), finite volumes (Bryan 1969) or finite element (e.g. Zienkiewicz and Taylor 1991). For example, Brown (1989) used the finite difference method to calculate the volume flow rate and electric current in simulated fractures composed of rough surfaces generated with a fractal algorithm. Rasouli and Hosseinian (2011) used the FEM based software (FLUENT) to develop a correlation to estimate the hydraulic parameters through channel of combined JRC profiles under different minimum closures. Indraratna et al. (2014) adopted the finite-volume method to solve the flow problem in deformable rough rock joints, where the three-dimensional Navier–Stokes equation was converted to an equivalent 2D flow model by considering the hydraulic aperture distribution. However, most of the traditional methods present the drawbacks such as long computation time, poor convergence and numerical instabilities, and the difficulties in dealing with complex boundaries (Wolf-Gladrow 2000). Alternatively, the ‘‘bottom up’’ approaches, such as lattice Boltzmann method (LBM) developed in the past two decades received more popularity in characterizing the flow problems. Originating from the 60.00% Resolution = 100

Maximum Relative error

50.00%

Resolution = 50 40.00%

Resolution = 20

30.00% 20.00% 10.00% 0.00% 0

0.5

1

1.5

2

Parameter k

Fig. 3 Impact of parameter choice of k

2.5

3

3.5

1 0.8 0.6

Analytical solution N = 100

0.4

N = 80 N = 60

0.2

N = 40

0 -0.5

-0.4

-0.3

-0.2

-0.1

0

0.1

0.2

0.3

0.4

0.5

y (distance from center line)

(a) 14.00%

Maximum Relative Errors

Parameters

Velosity (dimensionless)

Table 1 Parameters used in the simulation of fluid flow in single fracture

12.00% 10.00% 8.00% 6.00% 4.00% 2.00% 0.00% 0 5 10 15 20 25 30 35 40 45 50 55 60 65 70 75 80 85 90 95 100

Resolution

(b) Fig. 4 The simulation results based on the relationship dt = dx2. a Velocity profiles at different resolution, b maximum relative errors at different resolution and parameters

kinetic theory, the LBM has the appealing features of programming simplicity, intrinsic parallelism, and straightforward resolution of complex solid boundaries and multiple fluid species (e.g. Succi 2001; Higuera and Jime´nez 1989; Inamuro et al. 1995; He et al. 2010; Guo et al. 2002; Latt et al. 2008; Yan et al. 2011). For example, Eker and Akin (2006) presented studies of flow through two dimensional synthetically created fracture apertures using LBM. The permeability of fracture is found to be related to the mean aperture, fractal dimension and anisotropy factor of the synthetic fracture. However, the aforementioned numerical work is limited or simplified to 2D, the geometry description of the fracture is not accurate and the deformation cannot be involved properly. According to the literature review, there is still no comprehensive study on fluid flow in single fracture considering both the roughness and deformation. Moreover, the direct investigation of flow in fracture networks with roughness is rarely reported as well.

123

346

Therefore, it is important to explore the mechanism of fluid flow in the natural fracture and fracture network. This paper is structured as follows. Firstly, the fluid flow behavior in rough fracture is investigated considering the fracture’s deformation. The fracture is characterized by the mathematical model proposed by Brown (1995) and the fluid flow is simulated through LBM (Succi et al. 1995; Chen et al. 1998). The accurate of LBM for study of fluid flow through rock fracture is firstly verified through the comparison with the Poiseuille’s Law. After that, numbers of fluid flow simulations on realistic synthetic 3D rock fractures are conducted. A two parameters equation is developed based on the simulation results to predict the flow in rough fracture. The proposed equation can be used to characterize the fluid flow in single fracture involving both roughness and deformation. Then, LBM is used to investigate the fluid flow in natural fracture network, in which the roughness effect is directly incorporated. A good agreement is obtained between LBM simulation and the modified pipe network model using the derived empirical

Fig. 5 Fracture generator (Ogilvie et al. 2006)

123

Geomech. Geophys. Geo-energ. Geo-resour. (2016) 2:343–363

equation. Finally, the fluid flow behavior of stochastic DFN under deformation is preliminarily studied using LBM.

2 Single phase incompressible LBGK model 2.1 Basic concept In the incompressible LBGK model (Guo et al. 2000), the evolution equation of the density distribution function is expressed as fi ðx þ cei Dx; t þ DtÞ ¼ fi ðx; tÞ þ Xi ðfi ðx; tÞÞ ði ¼ 0; 1; . . .; MÞ

ð1Þ

where c ¼ Dx=Dt. Dx, ei and Dt are the lattice grid spacing, discrete velocity direction and time step, respectively. There are two commonly used lattice models for 2D and 3D problems (as illustrated in Fig. 1).

Geomech. Geophys. Geo-energ. Geo-resour. (2016) 2:343–363 Table 2 Parameters used to generate fracture for different fractal dimension, standard deviation and aperture Fractal dimension

SD (mm)

Mean aperture (mm)

1.0; 1.2; 1.4; 1.6; 1.8; 2.0

1; 1.5; 2; 2.5

2; 2.5; 3; 3.5; 4

347

Xi ðfi ðx; tÞÞ is the collision operator given by 1 Xi ¼  ðfi  fieq Þ s

ð2Þ

where s is the dimensionless relaxation time, and fieq is the equilibrium distribution function defined as

Fig. 6 Effect of fractal dimension on fracture profile (SD = 2 mm). a Fractal dimension = 1.2, b fractal dimension = 1.4, c fractal dimension = 1.6, d fractal dimension = 1.8

123

348

Geomech. Geophys. Geo-energ. Geo-resour. (2016) 2:343–363

Fig. 6 continued

fieq ðxÞ ¼ xi qðxÞ½1 þ si ðuÞ

ð3Þ

q¼

in which, xi is the weight index, and si ðuÞ ¼ 3

ei  u 9 ðei  uÞ2 3 u2 þ  2 c2 2 c4 2c

fi

ð5Þ

i¼1

ð4Þ

where u is the macroscopic velocity. The macroscopic density and velocity can be obtained as

123

M X

PM u¼

i¼1 fi ei

q

ð6Þ

Detailed explanation of the incompressible LBGK model can be found in the work by Guo et al. (2000), and

Geomech. Geophys. Geo-energ. Geo-resour. (2016) 2:343–363

349

Fig. 7 Effect of standard deviation on fracture profile (fractal dimension = 1.6). a SD = 1.0 mm, b SD = 1.5 mm, c SD = 2.0 mm, d SD = 2.5 mm

the corresponding incompressible Navier-Stoke equations were derived through multi-scaling expansion as, ru¼0

ð7Þ

ou þ r  ðuuÞ ¼ rp þ mr2 u ot

ð8Þ

pffiffiffi where p ¼ c2s q is the pressure, cs ¼ c= 3 is the sound speed, and m ¼ ð2s  1Þc2 Dt=6 is the kinetic viscosity.

2.2 Boundary conditions There are bunch of boundary conditions have been implemented in LBM. In this work, the boundary conditions are classified in two groups: the boundary condition at the open end (inlet and outlet) and the boundary condition at the solid interface. At the inlet and outlet, pressure boundary is applied to produce the pressure gradient. The no-slip boundary is used at the

123

350

Geomech. Geophys. Geo-energ. Geo-resour. (2016) 2:343–363

Fig. 7 continued

solid surface (wall), which is implemented through the bounce-back scheme. The so-called bounce-back means that when a fluid particle reaches solid (wall) nodes, the particle will scatter back to the fluid along with its coming direction as shown in Fig. 2. Both of the boundary conditions are implemented according to the work by Zou and He (1996).

123

2.3 Palabos The numerical simulations are conducted using the Palabos library (http://www.palabos.org/), which is a framework for general-purpose computational fluid dynamics (CFD) with a kernel based on LBM (Jonas Latt 2008). Its programming interface is

Geomech. Geophys. Geo-energ. Geo-resour. (2016) 2:343–363

351

Fig. 8 Geometry for the fracture and velocity distributions at slices from LBM simulation XZ Plane

YZ Plane

Flow direction (X)

Outlet 200 l.u.

Inlet

500 l.u. 500 l.u.

Fig. 9 Dependence of coefficients on fractal dimension and standard deviation of fracture surfaces. a Relation between a and fracture geometry, b relation between b and fracture geometry

0 -0.1 -0.2 SDT = 1 mm

-0.3

SDT = 1.5 mm

-0.4

SDT = 2.0 mm

-0.5

SDT = 2.5 mm -0.6 -0.7 -0.8 1

1.2

1.4

1.6

1.8

2

Fractal Dimension

(a) 0.3 0.2 0.1

SDT = 1 mm

0

SDT = 1.5 mm

-0.1 SDT = 2.0 mm

-0.2

SDT = 2.5 mm

-0.3 -0.4 -0.5 -0.6 1

1.2

1.4

1.6

1.8

2

Fractal Dimension

(b)

123

352

Geomech. Geophys. Geo-energ. Geo-resour. (2016) 2:343–363

straightforward, which makes it possible to set up fluid flow simulations with relative ease. Meanwhile, the LBM has the features of intrinsic parallelism. Programs written with Palabos can be automatically parallelized and the parallelization is performed with the message-passing paradigm of the MPI library. The Leonardi (http://leonardi.unsw.wikispaces.net/), a middle size cluster, is used to implement the parallel computation of the following numerical simulations.

3 Fluid flow in single fracture 3.1 Validation Firstly, the fluid flow between two parallel platesis simulated so as to verify the accuracy of Palabos. In LBM, it is necessary to convert the physical system to a discrete system so that the LBM simulation can be conducted. The flow domain is set as dimensionless system with length 2 and width 1. The pressure boundary condition is set at the left and right open side so as to produce unit velocity at the center, and the noslip boundary condition is set at top and bottom. And the parameters used in the simulation are presented in Table 1. The simulation is considered as converged when the ratio between the standard deviation and average of the velocity is less than 1e-6. The simulations are compared with the Poiseuille’s Law,  1  u ¼ rP  ðh=2Þ2  y2 ð9Þ 2qm where rP is the pressure gradient, q is the density, m is the viscosity, h is the aperture and y is the distance from the center line.

Table 3 Parameters in the fluid flow under shear displacement Physical size (cm)

595

Resolution

500 9 500

Fractal dimension

2.0

SD (mm) Initial aperture (mm)

2.0 5

Normal displacement (mm)

0; 1; 2

Shear displacement (mm)

0–10

123

Fig. 10 Assumptions made to implement shear displacement

It is necessary to mention that, discrete variables dx and dt are important parameters which have impact on the accuracy and the stability of a simulation, which should follow the relationship dt*dx2 (Jonas Latt 2008). Once the resolution changes, the time interval should change accordingly. However, the exact relationship has not been reported. In this study, the microscopic parameters are firstly calibrated and the dt = jdx2 is used to obtain a reasonable value for the choice of dt. The influence of j on simulation accuracy is investigated at different resolutions and relative error is summarized in Fig. 3. It is found that, when j = 1, stable results are obtained. Based on the calibrated relationship, dt = dx2, the resolution effect on the simulation accuracy is investigated. The velocity profile at different resolution, which is compared with the analytical solution (Fig. 4a), and relative error is presented in Fig. 4b.

Geomech. Geophys. Geo-energ. Geo-resour. (2016) 2:343–363

353

Fig. 11 Geometry of the fracture and the aperture change induced by displacement. a Assumption 1, b assumption 2

3 mm

2 cm 5 cm

(a) Shear displacement 2 mm

4 mm

6 mm

8 mm

(b) The flow rate along the pressure gradient can be obtained through the integration of the analytical solution, Z h=2  i dP 1 h 2 2 Qx ¼  ðh=2Þ  y  dy dl 2qm h=2 dP 1 h3   ð10Þ ¼ dl 2qm 12 Therefore, the transmissivity is described as the cubic law, T¼

h3 12

ð11Þ

(c)

3.2 Fluid flow through 3D synthetic fracture The parallel plate model can only be considered a qualitative description of flow through real fractures. Real fracture surfaces are not smooth parallel plates but are rough and contact each other at discrete points (Brown 1995). There are a number of parameters proposed to characterize the fracture roughness, such as Z2 (Myers 1962), joint roughness coefficient (Barton 1973), and fractal dimension (Xie et al. 1998). However, it is not always possible to characterize the fracture roughness by single parameter

123

354

Geomech. Geophys. Geo-energ. Geo-resour. (2016) 2:343–363 2 Normal displacement = 2 mm 1.5

Normal displacement = 1 mm

Aperture change (mm)

Normal displacement = 0 mm 1

0.5

0

-0.5 0

2

4

6

8

10

12

Shear displacement (mm)

(a) 0

Aperture change (mm)

-0.05

Normal displacement = 2 mm

-0.1

Normal displacement = 1 mm

-0.15

Normal displacement = 0 mm

-0.2 -0.25 -0.3 -0.35 -0.4 0

2

4

6

8

10

12

Shear displacement (mm)

(b) Fig. 12 Relationship between the aperture changes and shear displacement under two assumptions. a Assumption 1, b assumption 2 14

Transmissivity (mm3)

12

Assumption 1 Assumption 2

10

Equation 5.15 8 6 4 2 0 0

1

2

3

4

5

6

7

Mean aperture (mm)

In this part, a mathematical model developed by Brown (1995) is used to characterize the fracture roughness. In the mathematical model, the roughwalled fractures are dominated by three main parameters: the fracture dimension, the standard deviation of the surface profile, and a length scale describing the degree of mismatch between the two fracture surfaces. The software SynFrac (Ogilvie et al. 2006) is used to generate the synthetic fracture, which shares the same geometrical statistics as the natural fractures. The GUI of SynFrac is shown in Fig. 5 and the Brown model is selected among three modules. In this study, the fracture is assumed as two parallel surfaces without considering the mismatch effect. Therefore, there are three parameters to characterize the flow behavior in rough fracture, which are the fractal dimension, the standard deviation and the aperture. In order to investigate influence of roughness and develop a generic model to characterize the flow behavior in fracture, the physical size is set as 5 cm with resolution of 500 9 500, the mismatch length is set as 5 mm for reference only, and the other parameters are presented in Table 2. Examples of fracture profiles with different fractal dimension and standard deviation are presented in Figs. 6 and 7. A total number of 120 cases are generated based on the parameters in Table 2 which are imported to LBM for fluid dynamics simulation. The geometry of the 3D fracture and velocity profiles are presented in Fig. 8. The flow rate is calculated through the integration of the velocity at every lattice in the flow domain and the transmissivity of fractures under different profiles is derived. Based on the simulation results, a two coefficients model (Eq. 12) is developed to describe the fluid flow in a natural fracture. It is found from Fig. 9 that, the transmissivity decreases with the increase of fractal dimension and standard deviation according to the value of a and the proposed equation. However, the relationship between the coefficient b and fractal dimension is only clear when the standard deviation is larger than 1.5 mm.

Fig. 13 Relationship between the transmissivity and mean aperture

T¼

because the roughness of fracture surface in rock depends on the sample size or scale of observation (Bandis et al. 1981; Brown and Scholz 1985).

where T is the transmissivity, w is the width (equal to 1 in 2D), hm is the mean aperture, a and b are the coefficients that characterize the fracture roughness.

123

wð1 þ a  ebhm Þh3m 12

ð12Þ

Geomech. Geophys. Geo-energ. Geo-resour. (2016) 2:343–363

355

Fig. 14 Process of implement the fracture network in pipe network model and LBM simulation. a Random fracture network, b block cutting by DC, c pipe network model in Matlab, d geometry used in LBM simulation

3.3 Effect of displacement It is notable that the permeability of single fracture depends on fracture geometry and the stress condition (e.g. Pyrak-Nolte and Morris 2000; Koyama et al. 2008; Indraratna et al. 2014). Even with extensive studies on stress-flow coupled problem, there is still no existing law or principle to characterize the flow behavior in single fracture under specific fracture deformation. Therefore, it is necessary to investigate the flow behavior involves both the roughness and displacement. In this study, the fluid flow in single fracture is characterized by the mean aperture, hm , which is expressed as, hm ¼ h0  Dhn þ Dhs

ð13Þ

where h0 is the initial mechanical aperture, Dhn is normal displacement and Dhs is the shear displacement induced aperture change.

The effect of normal displacement is simple to change the mean aperture which has been investigated in the previous part. The following equation is suggested to predict of transmissivity of fracture under normal displacement,   w 1 þ a  ebðh0 Dhn Þ ðh0  Dhn Þ3 T¼ ð14Þ 12 The effect of shear displacement on fluid flow in rough fracture is studied as follows. The geometry of fracture is set as follows in Table 3. In reality, the two rough surfaces may overlap with each other at some points under shear displacement. However, the degree of overlap is not clear unless the stress-deformation relationship is known. Therefore, two assumptions are made to consider the fracture’s deformation. In assumption 1, there is only one contact point between the two surfaces, which means that one of the surfaces may skim over the opposite surface when they get contacted as shown in

123

Alplitude (cm)

356

Geomech. Geophys. Geo-energ. Geo-resour. (2016) 2:343–363 5 0 -5

0

10

20

30

40

50

Fracture length (cm)

(a) un

us

j Roughness distribution

i

Aperture

(b)

(c) Fig. 15 Generation of fracture network with roughness. a The geometry of individual fracture roughness, b implementation of roughness at individual fracture, c The geometry of fracture network Table 4 The coefficients of rough fracture Maximum amplitude, cm

Coefficients a

0

b 0

0

0.05

-0.1003

1.7314

0.10 0.15

-0.2448 -0.4632

0.1609 -0.7938

0.20

-0.5610

-0.4510

0.25

-0.7409

-0.6770

0.30

-0.9258

-0.5951

Fig. 10b. In assumption 2, the two surfaces will overlapand generate contact areas (damaged) as shown in Fig. 10c. Based on these two assumptions,

123

Fig. 16 Velocity distribution of fluid flow in fracture network under different maximum amplitude through LBM simulation. a Maximum amplitude = 0.05 cm, b maximum amplitude = 0.10 cm, c maximum amplitude = 0.15 cm, d maximum amplitude = 0.20 cm, e maximum amplitude = 0.25 cm, f maximum amplitude = 0.30 cm

one example of the generated fracture is demonstrated in Fig. 11. To visualize the shear displacement induced aperture change, the relationship between the mean aperture change and shear displacement under different normal displacements is summarized in Fig. 12. The flow simulations are conducted in the fracture under different shear displacement and normal displacement. The relationship between the transmissivity and mean aperture is shown in Fig. 13, which matches the proposed Eq. (15) well.

Geomech. Geophys. Geo-energ. Geo-resour. (2016) 2:343–363

357

Transmissivity (cm3)

0.00400 0.00350

Modified PNM

0.00300

LBM

Hpipe

 h3pipe 1 ¼ 12llpipe 1

1 1

ð17Þ

where hpipe is the aperture of the pipe, l is the dynamic viscosity of the fluid, and lpipe is the length of the pipe. A global flux pressure matrix is assembled in the way that the flow into a bubble equals the flow out of the bubble, so that the flux pressure relationship can be given as

0.00250 0.00200 0.00150 0.00100 0.00050 0.00000 0

0.05

0.1

0.15

0.2

0.25

0.3

0.35

Maximum amplitude (cm)

Fig. 17 Simulation results of fluid flow in fracture network involving roughness

  w 1 þ a  ebðh0 Dhn Dhs Þ ðh0  Dhn  Dhs Þ3 T¼ 12

Q ¼ Hglobal P

ð18Þ

where Q is the vector of bubble fluid flux, Hglobal is the global flux pressure matrix, and P is the vector of bubble pressure. Similar with classical FEM simulation, the global system matrix can be assembled from each single pipe (fracture) using Eq. (18). Fluid transmissivity of the fracture network can be obtained from solving the linear system equation together with the prescribed boundary conditions.

ð15Þ 4.2 Pipe model generation

4 Fluid flow through fracture network 4.1 Pipe network model (PNM)

In order to calculate the fluid flow in fracture network through the pipe network model, it needs to identify the location of the bubbles or joints and the links between them. To implement this, we generate the randomly sized polygonal blocks (Fig. 14a) using UDEC. The fractures of blocks are exported to the block cutting code of DDA to obtain the number index of the (bubbles) joints as well as the information of the links (pipes) (Fig. 14b). After that, information is imported to the Matlab program to generate the pipe network model (Fig. 14c). In order to verify the accuracy of LBM in simulating the fluid flow in fracture network, the corresponding LBM model is generated as shown in Fig. 14d.

In the pipe network model, the fluid flow in fracture network is represented by a discrete network made up from bubble and pipes, where the bubbles are the intersection points of the fractures, and the pipes are the links between two bubbles. The fluid flow variables are defined for each bubble as the fluid pressure p and fluid flux q. For each pipe, the fluid flux and pressure constitutive relationship can be obtained from Darcy’s law,



qi p ¼ Hpipe i ð16Þ qj pj

4.3 Fluid flow in natural fracture network

where qi and pi are the fluid flux and pressure of the ith bubble, Hpipe is the flux pressure matrix given by

As illustrated in the pipe network model, the flow in individual fracture of the fracture network model is

Table 5 Parameters of DFN models Block size (cm)

20.0 9 20.0

Density (/cm2)

0.4; 0.6; 0.8

Fracture length (cm)

Aperture (cm)

Mean

SD

Mean

SD

4; 5; 6; 7; 8

20 % of the mean value

0.3

0

123

358

Geomech. Geophys. Geo-energ. Geo-resour. (2016) 2:343–363

Fig. 18 DFN models with different fracture densities and lengths. a Density = 0.4, length = 4; b density = 0.4, length = 6; c density = 0.4, length = 8; d density = 0.6, length = 4; e density = 0.6, length = 6; f density = 0.6, length = 8

Fig. 19 Velocity distributions in DFN from LBM simulations. a Density = 0.6, length = 4; b density = 0.6, length = 6; c density = 0.6, length = 8; d density = 0.8, length = 4; e density = 0.8, length = 6; f density = 0.8, length = 8

characterized by the cubic law. However, the fracture is rough rather than a plate in most circumstances, and the influence of roughness on fluid flow is significant as explained in the previous part. Accordingly,

123

roughness effect on hydraulic behavior of fracture network receives increasing attention. However, most of studies require the empirical equations that relate the hydraulic aperture to mechanical aperture. For

Geomech. Geophys. Geo-energ. Geo-resour. (2016) 2:343–363 0.14 Density = 0.4

0.12

FLow rate (cm3/s)

Density = 0.6 0.1

Density = 0.8

0.08 0.06 0.04 0.02 0 4

5

6

7

8

9

Mean fracture length (cm)

Fig. 20 Dependence of flow rate on mean fracture length and density 0.14

FLow rate (cm3/s)

0.12 0.1 0.08 0.06 0.04 0.02 0

1

2

3

4

5

6

7

LD (/cm)

Fig. 21 Relationship between LD and flow rate 0.05 0.045

FLow rate (cm3/s)

0.04 0.035 0.03 0.025 0.02 0.015

LD = 1.5 cm

0.01

LD = 2.0 cm

0.005

LD = 2.5 cm

0 0.1

0.3

0.5

0.7

0.9

Fracture density (/cm2)

Fig. 22 Relationship between fracture density and flow rate under different LD

example, Zhao et al. (2014) adopted the dimensionless parameter Z2 proposed by Myers (1962) to obtain the hydraulic aperture from mechanical aperture.

359

However, the directly investigation on roughness effect on fluid flow through fracture network and the validation of empirical equation is still not reported. In order to investigate the roughness effect, the individual fracture between two intersection points is simplified as uniform distributed rough surface with constant aperture (Fig. 15a). The rough surface at individual fracture is implemented through the geometrical method illustrated in Fig. 15b. There are several parameters that characterize the geometry of fracture between intersection i and j, which includes the direction tensors (un , us ), the roughness distribution along us and aperture. The 2D fracture network with roughness is generated by creating the fracture between connected intersections (Fig. 15c). To obtain the two coefficients that characterize the roughness effect, the fluid flow through single fracture is firstly simulated using LBM. The length of fracture is 5 cm, the maximum amplitude ranges from 0.20 to 0.40 and the aperture ranges from 0.20 to 0.40. The parameters used in the simulation are the same as Table 1. The coefficients under different maximum amplitude are derived based on the simulation results which are shown in Table 4. The fluid flow through the fracture network is numerically implemented by using both the modified pipe network model and LBM. In the modified pipe network model, the transmissivity of pipes is calculated based on the proposed equation. In all the simulation, the topology is the same as the one shown in Fig. 14d with size of 20 cm 9 20 cm, the resolution in the LBM is set as 100 l.u./cm, and the aperture is set as 0.2 cm for all simulation. It is necessary to mention that, the velocity distribution in the single fracture can be simulated directly in the LBM (Fig. 16), whereas only the flow rate in the pipe network can be obtained for its macroscopic description. The simulation results from both models are summarized Fig. 17. Meanwhile, it is found that the simulation result from LBM is close to the computed transmissivity from pipe network model under cubic law assumption when no roughness is presented in the fracture network (amplitude = 0), which shows the ability of LBM in dealing with the fracture network flow problems. Furthermore, the proposed two coefficient equation can be employed directly into the pipe network model to represent the roughness effect, which can produce reasonable results through the comparison to the LBM.

123

360

Geomech. Geophys. Geo-energ. Geo-resour. (2016) 2:343–363

Fig. 23 Implementation of ‘‘indirect’’ hydromechanical coupling. a The boundary condition, b the flow region in LBM Rotation angle

Flow region

(a) 5 Fluid flow in discrete fracture network model (DFN) under deformation In the conventional methods that involve the DFN, it is always necessary to implement the block cutting before the flow simulation, the detection of the topology of fracture network as well as the isolated or dead-ends is a complex process. In contrast, the block cutting is not a problem for LBM, the isolated and dead-ends of fractures are naturally detected because they have no contribution to the flow at the micro scale. Meanwhile, it is always difficult or impossible to characterize the flow under large deformation, such as sliding, block rotation, aperture opening or closing, which may change the topology of the fracture network. For example, the existing flow channels may disappear during aperture closing, and the new flow channel may generate, which are not easy to be detected directly by using conventional method. In this section, the interconnectivity of fracture network and the anisotropic flow due to deformation is investigated by using LBM. 5.1 Fluid flow in DFN According to the work conducted by Long et al. (1982), the interconnection is dominated by the fracture density and fracture extent or size. To verify this, the fluid flow in the stochastic DFN is numerically simulated by using LBM. The geometry parameters of DFN are presented in Table 5. The location of the fracture and the orientation is uniformly distributed, and the aperture is set as 0.3 cm.

123

(b) The generated DFN model is imported LBM with a resolution of 100 l.u./cm so as to obtain relatively accurate results. Examples of the DFN models with different density and fracture length are shown in Fig. 18, where the black is the solid matrix, and the white is the flow channel. The pressure boundaries are set at the left inlet and right outlet. The velocity distribution is shown in Fig. 19, it is clear that there is no flow in the isolated fractures and dead-ends of fractures. The flow rate of DFN model under different mean length and density in Table 5 is calculated. All the simulation results from LBM are shown in Fig. 20. It is clear that, the larger density and longer length, the higher flow rate or degree of interconnection of the fracture network. In this part, we adopt the lengthdensity parameter LD introduced by Long et al.(1982) to summarize the results, which is shown in Fig. 21. The fracture network with larger LD tends to have higher flow rate and degree of interconnection. However, the fracture network have the same LD may behave differently from each other as indicated in the red rectangle in Fig. 21. Therefore, it is necessary to explore the mechanism behind this. We keep the LD constant as 1.5, 2.0 and 2.5, the fracture density varies from 0.1 to 0.8. The simulation results are shown in Fig. 22. It is notable that, under smaller LD, e.g. at LD = 1.5 cm, the fracture network with longer fracture and lower density will have higher flow rate and degree of interconnection. However, at larger LD, e.g. with LD = 2.5 cm, the flow rate increases with the fracture density to certain point and then decrease.

Geomech. Geophys. Geo-energ. Geo-resour. (2016) 2:343–363

0

30

0

30

0

30

0

30

361

(a)

(b)

(c)

(d)

60

90

60

90

60

90

60

90

Fig. 24 Conductivity of flow region with different orientation and axial strain. a Axial strain = 0, b axial strain =0.0002, c axial strain = 0.0006, d axial strain = 0.0008

5.2 Indirect coupling between LBM and DEM Hydro-mechanical coupling is an important consideration in fractured rock mass for rock mechanics and hydrogeology applications. A large number of attempts are conducted to investigate the influence of stresses on permeability of fracture network. The key factors that affect the flow behavior in fracture network include

opening, closure, sliding and dilation. For example, Min et al. (2004) conducted the simulation on fluid flow through DFN under stress condition. It is found that, the stress ratio is the main reason for dilation and the permeability of fracture network increases with the stress ratio under certain circumstances. In this part, an indirect hydro-mechanical coupling is conducted by using DEM (discrete element model)

123

362

Geomech. Geophys. Geo-energ. Geo-resour. (2016) 2:343–363

Fig. 25 Transmissivity of flow region with different orientation and axial strain

0.0008 0.00075

Transmissivity (cm3)

Axial Strain 0.0007

0 0.20%

0.00065

0.40% 0.0006

0.60% 0.80%

0.00055

1.00% 0.0005 0.00045 0.0004 0

30

60

90

120

150

180

210

Flow direction (Degree)

and LBM. The ‘‘indirect’’ coupling means that, the applied stress won’t produce the change in fluid pressure, and vice versa. Extensive work has been conducted in the study on stress-deformation relationship of the fracture network, which is not our key consideration. This work mainly focuses on the ability and advantage of LBM in characterizing the anisotropic flow behavior of rock masses under uniaxial deformation. To implement this, the stochastic DFN model is firstly generated in rock mass as shown in Fig. 23a. The displacement boundaries are set at the top and bottom with constant velocity of 0.01 cm/s. Detailed configuration of the mechanical model in UDEC can be found in the work by Kazerani and Zhao (2010). The geometry of fracture network is updated under certain axial strains. In order to characterize the anisotropic behavior, the conductivity in fracture network is measured in different directions. Figure 23b shows the flow region chosen with angle of h with respect to the horizontal plane. The LBM simulations are conducted on the deformed DNF models produced by UDEC, and the conductivity is measured relative to the orientation of the flow region. In each of the flow model, the aperture closing and opening due to mechanical displacement is naturally detected, which has shown the advantage of LBM in dealing with the hydro-mechanical coupled problems. The simulation results are presented in Fig. 24 and the velocity distribution of different flow direction and axial strain can be clear observed. The transmissivity

123

of flow regions under different axial strains are calculated, which is presented in Fig. 25. It is clear that, at the initial state without any displacement, the fracture network behaves as an isotropic media. The anisotropic flow behavior is observed when the vertical displacement is applied on the rock mass. The transmissivity at the horizontal direction decreases with axial strain because of the decrease of aperture due to compression. On contrary, the transmissivity increase as the axial strain increase which is induced by the aperture opening.

6 Concluding remarks The fluid flow through single rough fracture is extensively simulated in 3D through the LBM approach, which shows the ability of LBM in dealing with complex geometries. The two coefficients equation is proposed to characterize the flow behavior by considering both the roughness and displacements. The 2D fluid flows through fracture network are simulated by using LBM and the direct investigation of roughness effect on flow in fracture network is conducted. It is found that, the LBM and modified pipe network model could effectively take into account the roughness effect on fluid flow in fracture networks. Meanwhile, the interconnectivity of fracture network is investigated by LBM, and it is found that interconnectivity is dominant by the length density parameter and the correlation between the DFN parameter and

Geomech. Geophys. Geo-energ. Geo-resour. (2016) 2:343–363

transmissivity is analyzed. The ‘indirect’ hydromechanical coupled problem can be properly simulated by combing the distinct element method and LBM. The difficulty of flow prediction under large deformation could easily coped by LBM and the anisotropic flow behavior due to deformation is readily captured by the hydro-mechanical coupled analysis.

References Bandis S, Lumsden AC, Barton NR (1981) Experimental studies of scale effects on the shear behaviour of rock joints. Int J Rock Mech Min Sci Geomech Abstr 18(1):1–21 Barton N (1973) Review of a new shear-strength criterion for rock joints. Eng Geol 7:287–332 Barton NR, Bandis S, Bakhtar K (1985) Strength, deformation and conductivity coupling of rock joints. Int J Rock Mech Min Sci Geomech Abstr 22(3):121–140 Batchelor GK (1967) An introduction to fluid dynamics. Cambridge University Press, New York, pp 147–150 Brown SR (1987) Fluid flow through rock joints: the effect of surface roughness. J Geophys Res 92(B2):1337–1347 Brown SR (1989) Transport of fluid and electric current through a single fracture. J Geophys Res 94(B7):9429–9438 Brown SR (1995) Simple mathematical model of a rough fracture. J Geophys Res Solid Earth 100(B4):5941–5952 Brown SR, Scholz CH (1985) Broad bandwidth study of the topography of natural rock surfaces. J Geophys Res 90(12):12575–12582 Chen M, Bai M, Roegiers JC (1999) Permeability tensors of anisotropic fracture networks. Math Geol 31(4):335–373 Dershowitz WS, La Pointe PR, Doe TW (2004) Advances in discrete fracture network modeling. Proceedings. US EPA/ NGWA fractured rock conference, Portland, pp 882–894 Eker E, Akin S (2006) Lattice Boltzmann simulation of fluid flow in synthetic fractures. Transp Porous Media 65: 363–384 Guo Z, Shi B, Wang N (2000) Lattice BGK model for incompressible Navier–Stokes equation. J Comput Phys 165:288–306 He YL, Tao YJ, Yang LZ (2010) Experimental research on hydraulic behaviors in a single joint with various values of JRC. Chin J Rock Mech Eng 29(Supp. 1):3235–3241 Iwai K (1976) Fundamental studies of fluid flow through a single fracture. Ph.D. dissertation, University of California, Berkeley Indraratna B, Premadasa W, Brown ET, Gens A, Heitor A (2014) Shear strength of rock joints influenced by compacted infill. Int J Rock Mech Min Sci 70:296–307 Kazerani T, Zhao J (2010) Micromechanical parameters in bonded particle method for modelling of brittle material failure. Int J Numer Anal Meth Geomech 34(18): 1877–1895

363 Kim HM, Ryu DW, Tanaka T, Ando K (2007) Hydro-geological descriptive model using discrete fracture network (DFN) in a site characterization for subsea tunnels construction. Chin J Rock Mech Eng 26(11):2217–2225 Koyama T, Li B, Jiang Y, Jing L (2008) Coupled shear-flow tests for rock fractures with visualization of the fluid flow and their numerical simulations. Int J Geotech Eng 2(3):215–227 La Pointe PR, Hudson JA (1985) Characterization and interpretation of rock mass joint patterns. Geol Soc Am Spec Pap 199:1–38 Liu R, Jiang Y, Li B, Wang X, Xu B (2014) Numerical calculation of directivity of equivalent permeability of fractured rock masses network. Rock Soil Mech 35(8):2394–2400 Lomize GM (1951) Flow in fractured rocks. Gosenergoizdat, Moscow (in Russian) Long JCS, Witherspoon PA (1985) The relationship of the degree of interconnection to permeability in fracture networks. J Geophys Res 90(B4):3087–3098 Long JCS, Remer JS, Wilson CR, Witherspoon PA (1982) Porous media equivalents for networks of discontinuous fractures. Water Resour Res 18(3):645–658 Min KB, Rutqvist J, Tsang CF, Jing L (2004) Stress-dependent permeability of fractured rock masses: a numerical study. Int J Rock Mech Min Sci 41(7):1191–1210 Myers NO (1962) Characterization of surface roughness. Wear 5:182–189 Neuzil CE, Tracy JV (1981) Flow through fractures. Water Resour Res 1(3):191–199 Ogilvie SR, Isakov E, Glover PWJ (2006) Fluid flow through rough fractures. II: a new matching model for rough rock fractures. Earth Planet Sci Lett 241:454–465 Parker BL (2007) Investigating contaminated sites on fractured rock using the DFN approach. In: Proceedings of 2007 U.S. EPA/NGWA fractured rock conference: state of the science and measuring success in remediation, September 24–26. National Ground Water Association, Portland, Maine, Westerville Rasouli V, Hosseinian A (2011) Correlations developed for estimation of hydraulic parameters of rough fractures through the simulation of JRC flow channels. Rock Mech Rock Eng 44(4):447–461 Tsang YW (1984) The effect of tortuosity on fluid flow through a single fracture. Water Resour Res 20(9):1209–1215 Witherspoon PA, Wang JSY, Iwai K, Gale JE (1980) Validity of cubic law for fluid flow in a deformable rock fracture. Water Resour Res 16(6):1016–1024 Xie H, Wang JA, Stein E (1998) Direct fractal measurement and multifractal properties of fracture surface. Phys Lett A 242(1–2):41–50 Zhao ZH, Li B, Jiang YJ (2014) Effects of fracture surface roughness on macroscopic fluid flow and solute transport in fracture networks. Rock Mech Rock Eng 47(6):2279–2286 Zimmerman RW, Bodvarsson GS (1996) Hydraulic conductivity of rock fractures. Transp Porous Media 23(1):1–30 Zou Q, He X (1996) On pressure and velocity flow boundary conditions for the lattice Boltzmann BGK model. Phys Fluids 9:1591–1598

123