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 com putational 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 latt ice Boltzmann method (LBM)
in a middle size cluster. An empirical equatio n 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 applicabil-
ity. 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.
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 defor-
mation. An accurate prediction of hydraulic behavior
in fracture network requires a clear under standing of
fluid flow though single fracture under these coupled
conditions.
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: gaofeng.zhao@tju.edu.cn
123
Geomech. Geophys. Geo-energ. Geo-resour. (2016) 2:343–363
DOI 10.1007/s40948-016-0040-4
Comprehensive work has been conducted to inves-
tigate 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 shear-
flow tests for rock fractures. Indraratna et al. (
2014)
investigate the fluid flow through deformable rough
rock joints. However, most of the works relate the
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 Bod-
varsson
1996). Theoretically, the flow of incompress-
ible 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 pro-
vided 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
(a) (b)
1
11(7)
5
2
6
1
5
2
6
3
7
4
8
10(14)
13(9) 8(12)
y
x
2D
z
x
3D
y
Fig. 1 The D2Q9 model
and D3Q15 model. a D3Q15
model, b D2Q9 model
f
5f
6
f
3
f
7
f
4
f
8
f
2
f
1
f
5
f
6
f
3
f
7
f
4
f
8
f
2
f
1
Streaming
wall
wall
Fig. 2 Bounce back
scheme
344 Geomech. Geophys. Geo-energ. Geo-resour. (2016) 2:343–363
123
(e.g. Ames 1977; Morton and Mayers 1994), finite
volumes (Bryan 1969) or finite element (e.g. Zienkie-
wicz 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 con-
verted 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 numer-
ical 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 charac-
terizing the flow problems. Originating from the
kinetic theory, the LBM has the appealing features
of programming simplicity, intrinsic parallelism, and
straightforward resolution of complex solid bound-
aries 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 permeabil-
ity 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.
Table 1 Parameters used in the simulation of fluid flow in
single fracture
Parameters Values
Density (q) 1.0
Reynolds number (Re) 1.0
Resolution (N, l.u.) 10–100
dx 1/N
dt
*dx
2
0.00%
10.00%
20.00%
30.00%
40.00%
50.00%
60.00%
0 0.5 1 1.5 2 2.5 3 3.5
Resolution = 100
Resolution = 50
Resolution = 20
Parameter k
Maximum Relative error
Fig. 3 Impact of parameter choice of k
(a)
(b)
0
0.2
0.4
0.6
0.8
1
1.2
-0.5 -0.4 -0.3 -0.2 -0.1 0 0.1 0.2 0.3 0.4 0.5
Analytical
solution
N = 100
N = 80
N = 60
N = 40
y (distance from center line)
Velosity (dimensionless)
0.00%
2.00%
4.00%
6.00%
8.00%
10.00%
12.00%
14.00%
0 5 10 15 20 25 30 35 40 45 50 55 60 65 70 75 80 85 90 95 100
Resolution
Maximum Relative Errors
Fig. 4 The simulation results based on the relationship
dt = dx
2
. a Velocity profiles at different resolution, b maximum
relative errors at different resolution and parameters
Geomech. Geophys. Geo-energ. Geo-resour. (2016) 2:343–363 345
123
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 simula-
tion 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
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
f
i
ðx þ ce
i
Dx; t þ DtÞ¼f
i
ðx; tÞþX
i
ðf
i
ðx; tÞÞ
ði ¼ 0; 1; ...; MÞð1Þ
where c ¼ Dx=Dt. Dx, e
i
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).
Fig. 5 Fracture generator (Ogilvie et al. 2006)
346 Geomech. Geophys. Geo-energ. Geo-resour. (2016) 2:343–363
123
X
i
f
i
ðx; tÞðÞis the collision operator given by
X
i
¼
1
s
f
i
f
eq
i
ðÞ ð2Þ
where s is the dimensionless relaxation time, and f
eq
i
is
the equilibrium distribution function defined as
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
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
Geomech. Geophys. Geo-energ. Geo-resour. (2016) 2:343–363 347
123