A pore-scale numerical model for non-Darcy fluid flow through rough-walled fractures

Computers and Geotechnics 87 (2017) 139–148

Contents lists available at ScienceDirect

Computers and Geotechnics journal homepage: www.elsevier.com/locate/compgeo

Research Paper

A pore-scale numerical model for non-Darcy fluid flow through roughwalled fractures Wei Zhang a,b, Beibing Dai a,b, Zhen Liu a,c, Cuiying Zhou a,b,⇑ a

Research Center for Geotechnical Engineering & Information Technology, Sun Yat-sen University, Guangzhou 510275, China School of Engineering, Sun Yat-sen University, Guangzhou 510275, China c School of Marine Sciences, Sun Yat-sen University, Guangzhou 510275, China b

a r t i c l e

i n f o

Article history: Received 18 October 2016 Received in revised form 3 January 2017 Accepted 13 February 2017

Keywords: Non-Darcy flow Pore-scale Seepage Fracture Numerical Model

a b s t r a c t This paper aims at developing a pore-scale numerical model for non-Darcy fluid flow through roughwalled fractures. A simple general relationship between the local hydraulic conductivity and the flow velocity is proposed. A new governing equation for non-Darcy fluid flow through rough-walled fractures is then derived by introducing this relationship into the Reynolds equation. Based on the non-linear finite element method, a self-developed code is used to simulate the non-Darcy fluid flow through fractures. It is found that the macroscopic results obtained by the numerical simulation agree well with the experimental results. Furthermore, some interesting experimental observations can be reproduced. Ó 2017 Elsevier Ltd. All rights reserved.

1. Introduction The natural rock mass in geotechnical engineering usually consists of rock matrix and fractures at different scales. In many situations, the permeability of rock matrix is negligible in comparison with that of fractures, and thus the hydraulic behaviour of rock mass is controlled by fractures. In many fields related to rock mechanics, such as hydraulic engineering, slope engineering, mineral resource exploitation, underground oil storage, geothermal development, nuclear waste disposal and so on, the fluid flow in rock fractures plays an important role. According to statistical results, 30–40% of hydropower project dam failures, more than 90% of rock slope failures, and around 60% of mine accidents were related to groundwater. Therefore, the research on the fluid flow through rock fractures is an important fundamental issue. Traditionally, the flow through fracture was treated as linear laminar flow. Individual fractures were idealized as parallel plates in order to obtain a tractable mathematical description, namely, the cubic law [1]. However, it has been widely accepted that the surfaces of most fractures are rough-walled and the fracture aperture is inhomogeneous, and thus the macroscopic flow behaviour

⇑ Corresponding author at: Research Center for Geotechnical Engineering & Information Technology, Sun Yat-sen University, Guangzhou 510275, China. E-mail address: [email protected] (C. Zhou). http://dx.doi.org/10.1016/j.compgeo.2017.02.011 0266-352X/Ó 2017 Elsevier Ltd. All rights reserved.

in a single fracture can be determined by the pore-scale void geometrical properties of the fracture. The local cubic law, which assumes that the cubic law is applicable to each local void space to estimate the local transmissivity value, has been widely examined for calculating the macroscopic discharge through fractures with inhomogeneous apertures [2–6]. So far, many experimental studies on the fluid flow in rough-walled fracture have been carried out. The phenomena of tortuous flow path and preferential flow channels have been observed. The concepts of aperture, roughness, area contact ratio, tortuosity, fractal dimension have been proposed to describe the fracture geometrical properties. The effect of the fracture roughness on the linear flow in roughwalled fractures has been well studied. However, more studies show that the non-Darcy flow is a significant existence. With the increase of flow rate, deviation from linear Darcy’s law may occur [7–12]. Recently, many experimental studies on non-Darcy flow through rough-walled fractures have been carried out. For example, Zhang and Nemcik [13] investigated the fluid flow regimes and non-Darcy flow characteristics through rough-walled fractures, by conducting water flow tests through both mated and non-mated sandstone fractures. Javadi et al. [14] experimentally investigated the role of shear processes on the variation of critical Reynolds number and non-Darcy flow through rough-walled fractures. Chen et al. [15] experimentally evaluated the Forchheimer equation coefficients for non-Darcy flow through

140

W. Zhang et al. / Computers and Geotechnics 87 (2017) 139–148

rough-walled fractures. Singh et al. [16] carried out water flow tests through single fractured granite and quantified the sources of non-linear relationship between discharge and base pressure. Zhang et al. [17] carried out triaxial flow experiments to investigate the macroscopic flow behaviour through rough-walled fractures, and also examined the microscopic flow behaviour by introducing an advanced micro-fluidic technology. Although experimental investigations play an important and irreplaceable role in the study of the fluid flow through roughwalled fractures, there exist some deficiencies due to the limitation of experimental techniques. Besides that the cost is higher and the observed data is less, the measurement of pore-scale physical quantity and the reveal of pore-scale mechanisms are insufficient. Therefore, pore-scale numerical simulation has become an important supplementary means. Although some pore-scale numerical investigations have been carried out to improve the understanding of the linear laminar fluid flow through rough-walled fractures (e.g. [18–26]), only a few studies focused on the non-Darcy flow in fractures. Javadi et al. [27] generated the computational domain of an artificial three-dimensional fracture and performed both laminar and turbulent flow simulations. Trykozko et al. [28] studied the flow in hydraulically fractured formations and paid special attention to the conductivity reduction due to non-Darcy effects at high flow rates. Both of these studies applied the technique of direct numerical simulations (DNS), i.e. direct numerical solution of the Navier-Stokes (N-S) equations describing the flows at pore-scale. Theoretically, the non-Darcy flow in rough-walled fractures can be fully described by the N-S equation. However, DNS is generally rather time-consuming, and thus it is necessary to develop a less rigorous but more computationally efficient numerical model. This is the main reason why the Reynolds equation has been widely used as an alternative approximation of linear flow in fractures (e.g. [6,29–31]), although the DNS results have shown that the parabolic vertical profiles of flow velocity assumed in the Reynolds equation may not be tenable [32]. This study proposes a pore-scale numerical model for nonDarcy fluid flow through rough-walled fractures. In the model, the inhomogeneous aperture distribution of fractures is generated based on the random field theory. Starting with the N-S equations, a simple general relationship between the local hydraulic conductivity and the flow velocity is proposed. A new governing equation for non-Darcy fluid flow through rough-walled fractures is then derived by introducing this relationship into the Reynolds equation. By resorting to the non-linear finite element method for the solution of the governing equation, a self-developed code is used to simulate the non-Darcy fluid flow through rough-walled fractures. Eventually, the results from the pore-scale numerical simulation are compared with the experiments, and the macroscopic and pore-scale non-Darcy fluid flow through rough-walled fractures are investigated numerically.

2. Geometrical model of rough-walled fractures In order to simulate the non-Darcy flow through rough-walled fractures, it is essential to generate the homogeneous aperture field representing a rough-walled fracture. Many laboratory observations suggest that the void aperture field between two roughwalled fracture surfaces follows a log-normal or normal distribution (e.g. [23,33,34]). This type of aperture field distribution has been used in most numerical studies to generate fracture void geometry (e.g. [18–20,23]). In this study, in order to generate the homogeneous aperture field, the two rough surfaces of a fracture are generated separately, and then they are brought together with both tangent and normal position offsets. Thus the void with variable shape and aperture, as

well as the regions of contact between the surfaces, are formed. This generation process is similar to those of natural rock fractures and the sample preparation procedures in recent experiments (e.g. [15–17]). In order to generate the rough fracture surface, the height of a rough fracture surface is assumed to be random variables characterized statistically by normal distribution. The covariance matrix decomposition method [35] is adopted to generate a random field. The isotropic autocovariance function is adopted as[36]:


 4h AðhÞ ¼ r2 exp  l

ð1Þ

where, r is the standard deviation, l is the spatial correlation length, and h is the distance between the two locations of the fracture plane. Fig. 1 illustrates two numerical samples of the fracture aperture field generated by the present method. For each sample, the sample size is 50 mm  100 mm and it is divided into 50  100 square elements with 51  101 nodes. The two surfaces of a fracture are assumed to possess the same geometrical shape, and the spatial correlation length is 5 mm. The tangent and normal position offsets between the two fracture surfaces are 0 mm and 0:5 mm, respectively. The only difference of the two samples lies in that the standard deviation r of sample 1 is 0:25 mm and that of sample 2 is 0:5 mm. From Fig. 1 it can be seen that sample 1 is smoother than sample 2. 3. Pore-scale numerical model for non-Darcy flow 3.1. Theoretical background For fluid flow through rough rock fractures, the motion of fluid obeys the basic principles of fluid dynamics, i.e. the mass, momentum and energy conservations. The well-known N-S equations can accurately describe the travel of the fluid in rock fractures. For incompressible Newtonian viscous fluid, such as water, the N-S equations in the steady-state can be given as follows:

lr2 u  qðu  rÞu ¼ rp

ð2Þ

where q is the fluid density, p is the hydraulic pressure, l is the fluid viscosity, and u is the flow velocity vector. The first term

lr2 u represents the viscous forces; the second term qðu  rÞu is the advective acceleration term, and represents the inertial forces; and the third term rp represents the pressure gradient. The presence of the advective acceleration term generally causes the equations to be non-linear, and consequently makes the computation time-consuming. The requirement for convenient engineering prediction has stimulated the development of simple governing laws. For example, the linear Darcy’s law has been widely utilized to describe the laminar flow. By assuming that a single rock fracture consists of parallel planar plates, the cubic law [1] can be derived from the N-S equations: rp ¼

12lQ we3

ð3Þ

where Q is the volumetric flow rate, w is the fracture width and e is the hydraulic aperture. The assumption of parallel planar model is overly-simplified, because natural rock fractures have rough-walled walls and spatially variable apertures, as well as regions where the two opposing fracture walls are in contact with each other. Considering the effect of the aperture variation, the Reynolds equation has been widely examined for calculating the flow rate through rough-walled fractures [2–6]. The Reynolds equation is derived from the N-S equa-

141

(a)

(b)

[mm]

100

[mm]

100

2.5

1.4 75

75

1.2

2

Frequency [−]

W. Zhang et al. / Computers and Geotechnics 87 (2017) 139–148

0.1 0.05 0

0.6 0.4

25

50

25 X [mm]

50

0.5 1 1.5 Aperture [mm]

(d) 1

25

0.5

0.2 0 0

0

1.5

0 0

25

Frequency [−]

0.8

Y [mm]

Y [mm]

1.0 50

(c)

0.15

0.4 0.3 0.2 0.1 0 0 0.5 1 1.5 2 2.5 3

50

Aperture [mm]

X [mm]

Fig. 1. (a) Aperture distribution of Sample 1; (b) aperture distribution of Sample 2; (c) aperture histogram of Sample 1; (d) aperture histogram of Sample 2.

tion by assuming that inertial forces in the fluid are negligibly small compared with the viscous and pressure forces, and the changes in aperture occur gradually. The Reynolds equation can be written as



 @ @h @ @h þ þ Qs ¼ 0 k k @x @x @y @y

ð4Þ

where h is the hydraulic head, Q s is the source/sink term, and k is the hydraulic conductivity and can be calculated as follows:

k¼

qge2 12l

ð5Þ

According to the Reynolds equation, the fluid flow through rough fractures can be described by merely assuming that the cubic law holds locally at each point in the fracture, and then invoking the principle of mass conservation.

Fig. 2. Side view of a cross-section of a rough-walled rock fracture.

and the order of magnitude of the inertial terms is

mag½qðu  rÞu 

3.2. Generalized Reynolds equation For non-Darcy flow through rough fractures, the inertial forces in the fluid cannot be neglected, and thus the traditional Reynolds equation is not applicable. However, according to the principle of mass conservation, a generalized Reynolds equation can be obtained as:



 @ @h @ @h þ þ Qs ¼ 0 kg ðv Þ k g ðv Þ @x @x @y @y

ð6Þ

in which kg ðv Þ is the generalized hydraulic conductivity. Due to the presence of non-linear inertial term, for a local point in the fracture, kg ðv Þ is not constant, but is a function of the flow velocity. 3.3. Relationship between the local hydraulic conductivity and the flow velocity The relationship between the local hydraulic conductivity and the flow velocity is non-linear. It is impossible to derive this nonlinear relationship accurately from the N-S equations. However, it is available to estimate the orders of magnitude of the three terms in Eq. (2). Let U be a characteristic magnitude of the velocity, which could be thought of as the average velocity along the thickness of fracture. Define a characteristic length K in the x direction, which may be the wavelength of the aperture variations (as shown in Fig. 2). Hence the order of magnitude of the viscous terms lr2 u can be estimated to be [5]

mag½lr2 u 

lU e2

ð7Þ

qU 2

ð8Þ

K

In order to estimate the proportion of the inertial forces, the reduced Reynolds number Re can be defined to be the ratio between the order of magnitude of the inertial terms and that of the viscous terms, as follows:

Re ¼

mag½qðu  rÞu mag½lr u 2

¼

qUe2 lK

ð9Þ

It can be seen that the reduced Reynolds number Re is the product of the traditional Reynolds number Re ¼ qUe=l and the geometrical parameter e=K. When the flow velocity is low enough, Re is far less than 1. The inertial forces are negligible, and the flow through fractures is laminar flow. In this case, the linear Darcys law and Reynolds equation are applicable. With the increase of the flow velocity, Re gradually increases, indicating that the proportion of inertial forces gradually increases and the inertial forces become not negligible. Thus a gradual transition from laminar flow to fully developed turbulent flow occurs. For the fully developed turbulent flow, Re is far greater than 1. In this case, which is opposite to the above one, the viscous forces become negligible. According to Eqs. (2) and (8), the pressure gradient is in proportion to the square of flow velocity, i.e.

rp ¼ kv 2

ð10Þ

where k is a constant coefficient. From Eqs. (6) and (10), the relationship between the local hydraulic conductivity and the flow velocity for the fully developed turbulent flow can be calculated as: 0

k ¼

qg kv

ð11Þ

142

W. Zhang et al. / Computers and Geotechnics 87 (2017) 139–148

By using the critical Reynolds number Rec to distinguish the laminar flow and turbulent flow and ignoring the transitional flow regime, a general Darcian-like relationship between the local hydraulic conductivity and the flow velocity can be obtained as:

( k g ðv Þ ¼

k ¼ q12gel ; Re < Rec 2

0

qg

k ¼ kv ;

ð12Þ

Re P Rec

From Eq. (12), it can be seen that when the Reynolds number Re is smaller than the critical Reynolds number Rec , the local hydraulic conductivity k is constant, indicating that the fluid flow is linear. When Re is greater than Rec , the local hydraulic conductivity gradually decreases with the increase of the velocity, and the fluid flow is non-Darcian. As the critical Reynolds number Rec is independent of the local hydraulic aperture, it can be seen as a constant parameter for varied hydraulic apertures. According to the continuity condition of Eq. (12), when Rec is known, the coefficient k can be calculated as:

k¼

12q eRec

ð13Þ

Thus, we have derived a simple generalized relationship between the local hydraulic conductivity and the flow velocity. Compared with the classic local cubic law, there is only one more independent parameter, the critical Reynolds number Rec . The critical Reynolds number Rec , as a pore-scale parameter, cannot be directly measured through experiments. However, it can be determined through a careful parameter calibration procedure based on the macroscopic experimental results. Fig. 3 shows an example for this relationship. The critical Reynolds number Rec is 50, the den3

sity is 103 kg=m , the gravitational acceleration g is 9:8 m=s2 and the fluid viscosity is 103 Pa  s. 3.4. Numerical implementation The generalized Reynolds Eq. (6) can be solved numerically by using a standard non-linear finite-element approach. The finiteelement based discrete formulation is:

½Kðv Þfhg ¼ fFg: where

Kðv Þ ¼ Z F¼

C

ð14Þ

Z X

BT DBdX

ð15Þ

q/dC

ð16Þ

Hydraulic conductivity [m/s]

0.035 0.1mm 0.15mm 0.2mm

0.03 0.025 0.02 0.015 0.01 0.005 0 0

0.2

0.4

0.6

0.8

1

Velocity of flow [m/s] Fig. 3. An example for the relationship between the local hydraulic conductivity and the flow velocity with different local hydraulic aperture.

Bij ¼ /j;i ;

i ¼ 1; 2;

j ¼ 1; 2; 3; 4

ð17Þ

where / is the shape function, and D is the non-linear hydraulic conductivity matrix and can be expressed as:



D¼

kg ðv Þ



0

ð18Þ

k g ðv Þ

0

Newton’s method is adopted for the solution of the non-linear discrete formulation (14). The iteration formulation is given as:

fWk g ¼ ½Kðv k Þfh g  fFg k

kþ1

fh

1

k

g ¼ fh g  ½K kt  fWk g

ð19Þ

where the superscripts k and k þ 1 represent the iteration steps. The iteration convergence condition is given as k

e¼

kh  h

kþ1

kþ1

kh

k1

k1

<2

ð20Þ

where 2 is the convergence precision and is taken to be 0:01 in this study. 4. Model verification and numerical study Recently, Zhang and Nemcik [13] carried out water flow experiments through both mated and non-mated sandstone fractures in triaxial cell with confining stress varying from 1:0 MPa to 3:5 MPa. They found that linear Darcy’s law holds for water flow through mated fracture samples, due to high flow resistance caused by the small aperture and high tortuosity of the flow pathway, while non-Darcy flow occurs for non-mated fracture due to the enlarged aperture. In terms of non-mated fracture, totally 4 samples with different rock joint profiles were prepared. According to the experimental results, the relationship between pressure gradient and volumetric flow rate of Sample 4 is very close to that of Sample 1. Therefore, in this study, the flow experiments of non-mated fracture Samples 1–3 under the confining stress of 1:0 MPa were chosen for numerical simulation. The macroscopic relationships between the pressure gradient and volumetric flow rate obtained by the flow experiments are used to verify the present numerical model. Furthermore, the pore-scale fracture flow behaviour is studied numerically. Note that under the present experimental technology, many pore-scale physical quantities cannot be measured. 4.1. Experimental setup The procedures of non-mated fracture sample preparation are illustrated in Fig. 4. Intact cylindrical samples were cored from a fine grain sandstone block in the laboratory. The cored samples were then split into halves using two types of splitting wedges (straight and sinusoidal-shaped), to achieve different rock joint profiles. In order to obtain ultimate non-mated fracture sample, the two fracture halves were displaced by 2 mm along the flow direction, and then the extended parts were cut off and the sample ends were polished. After the intact sample was split into two halves, the fracture surface roughness was measured using a non-contact 3D laser scanner in fine scanning mode. Based on scanned data, the peak asperity height n was estimated by averaging over the nine equally spaced profiles parallel to the flow direction. And then the Joint Roughness Coefficient (JRC) was estimated according to the empirical equation of Barton and de Quadros [37]:

JRC ¼ 400

n L

ð21Þ

143

W. Zhang et al. / Computers and Geotechnics 87 (2017) 139–148

Flow out

Splitting

Mis-matching

Cutting and Polishing

Flow in Fig. 4. Procedures of non-mated fracture sample preparation.

where L is the fracture length. Finally, the water flow experiments were conducted in triaxial cell. During the flow experiments, pressurized water from the tank enters the fracture at the bottom of the sample, then exits through the top cap. When the flow was stable at a given pressure, the flow rate was measured using a precise electrical balance scale.

4.2. Numerical simulation

q ¼ 1:0  103 kg=m3 and viscosity l ¼ 1:03 Pa  s are used. The

In order to reproduce the flow experiments numerically, it is essential that the aperture distribution of the numerical sample is similar to that of the experimental sample. In this study, the standard deviation r of random field, which is related to the roughness of the fracture surface, is determined from the joint roughness coefficients (JRC) on the basis of the following empirical equation [38],

JRC ¼ 32:2 þ 32:47 ln

r

Dx

Table 1 provides the geometrical parameters of the three numerical samples. The values of tangent offset Os for the three samples are all 2:0 mm, which is consistent with the sample preparation process in the experiments. The values of normal offset On , which were not given in the literature [13] due to the difficulty of measurement, are derived through a careful parameter calibration procedure. In the numerical simulation, the water density

ð22Þ

where Dx is the element size discretizing the random field. It should be mentioned that the JRC value is defined in 1-D, i.e. along the measurement line; however, the fracture aperture distribution is in 2-D. In the experiments, for each fracture, nine equally spaced measurement lines parallel to the flow direction are used. The peak asperity height n was estimated by averaging over the nine measurement lines, and then the JRC value was estimated according to an empirical equation. Therefore, the 2-D aperture distribution characteristics of fracture samples had been reflected. It is regretful that no more statistical information about the spatial distribution of fracture apertures can be found in the literature [13]. However, after the random field corresponding to the fracture surface was generated numerically, the peak asperity height n of the generated surface can be estimated in the same way as adopted in the laboratory experiments. That is to say, it can be calculated by averaging over the nine equally spaced lines parallel to the flow direction. Due to the randomness of the fracture surface generation process, we can generate a fracture surface repeatedly until the peak asperity height n is consistent with that estimated based on the laser scanned data in the experiment. After the fracture surface was generated, the fracture aperture distribution was generated by bringing the two surfaces with the same geometrical shape together, with both tangent and normal position offsets from each other. For each numerical sample, the two surfaces are assumed to possess the same geometrical shape. Since the fractures in the experiments were generated by splitting the intact cylindrical samples, this assumption is rational.

value of critical Reynolds number Rec is calibrated to be 50. The boundary conditions are given as follows (Fig. 5): the constant water head boundary is applied to the top of fracture plane, the constant flow rate boundary is applied to the bottom of fracture plane, and the impermeable boundaries are applied to both sides of fracture plane. 4.3. Macroscopic fracture flow behaviour Fig. 6 shows the numerical results of the relationship between the pressure gradient and volumetric flow rate for the three samples, together with the corresponding experimental results which were obtained by the regression analysis with Forchheimer equation. Forchheimer equation has been widely used to macroscopically characterize the non-Darcy flow process, which can be written as:

dp 2 ¼ aQ þ bQ dx

ð23Þ

where a and b are model coefficients, representing pressure drop components caused by linear and non-linear effects, respectively. For all the three samples, the numerical results agree well with the experimental results, indicating that the numerical model can reasonably predict the macroscopic non-Darcy flow behaviour through the rough-walled fractures. In the relationship between the local hydraulic conductivity and the flow velocity, there is only one more independent parameter, the critical Reynolds number Rec , as compared with the classic local cubic law. To investigate the effect of Rec on the non-Darcy flow, the values of Rec are varied in the simulation of the flow through Sample 1. Fig. 7 shows the relationship between the pressure gradient and volumetric flow rate at different Rec . It is clear that the deviation degree from linear flow gradually increases with the decrease of Rec . The reason is that the local flow tends to be turbulent when the Rec is small. Thus, the overall flow resistance is higher and this increases the non-linear degree. When the value of Rec is adequately large (e.g. Rec ¼ 500), the local turbulent flow hardly occurs, and thus the overall fracture flow is linear. As men-

144

W. Zhang et al. / Computers and Geotechnics 87 (2017) 139–148

Table 1 Computation parameters of the three samples employed in the numerical simulation. Sample No.

Fracture width W (mm)

Fracture length L (mm)

JRC

1 2 3

50.0 50.0 50.0

100.0 100.0 100.0

5.5 7.0 9.2

Standard deviation 0.440 0.460 0.492

Impermeab ble boundary

100 0 mm

Impermeable boundary

Normal offset On (mm)

2.0 2.0 2.0

0.268 0.330 0.330

4.4.2. Distribution of Reynolds number Fig. 9 shows the distribution of Reynolds number Re. It can be seen that the region of Re < 50 occupies the majority of the fracture. However, the region of Re > 50 can be found locally. In this region, the fracture flow is turbulent, and the local cubic law cannot be applied to estimate the fracture flow rate. Therefore, the present numerical model is more suitable than the local cubic law. It can be concluded that due to the inhomogeneous flow distribution, the locally developed turbulent flow contributes to the macroscopic non-linearity of the fracture flow.

y

ow direction Macroscopic flo

50 mm

0

Tangent offset Ot (mm)

order to analyze the pore-scale flow behaviour, the water flow through Sample 1 with the hydraulic gradient dp=dx ¼ 6 MPa=m is chosen as a representative. Fig. 8(a) shows the distribution of flow rate. Obvious inhomogeneous flow distribution and preferential flow channels can be observed. Fig. 8(b) shows the distribution of flow velocity. It is confirmed that the flow distribution is inhomogeneous. Because the relationship between the local hydraulic conductivity and the flow velocity is based on the local cubic law, the inhomogeneous aperture distribution leads to a much more inhomogeneous flow distribution. In Fig. 8(b), an enlarged view of flow velocity distribution is also shown. Note that an obvious water flow around the contact region can be observed. It is clear that the present numerical model can well capture the pore-scale flow behaviour, such as the inhomogeneous flow distribution and the preferential flow channels.

Constant water head boundary

Model of a fracture with spatially variable apertures

r (mm)

x Constant flow rate boundary

Fig. 5. Boundary conditions in the pore-scale numerical simulation of flow through a rock fracture with spatially variable apertures.

tioned above, when Rec ¼ 50, the non-linear degree of the numerical results fit the experimental results well for all the three samples, indicating that the local critical Reynolds number Rec might physically be close to 50.

4.4.3. Tracing of flow The flow paths from several points are traced using the velocity simulated by the present model. Fig. 10 shows these flow paths. In this figure, it can be observed that each flow path is tortuous. The lengths of flow paths 1–4 are respectively 125:8 mm; 131:1 mm; 126:9 mm and 131:1 mm. It is clear that the actual flow path is 25.8–31.1% longer than the straight line. In general, the hydraulic gradient is defined by the ratio of the

4.4. Pore-scale fracture flow behaviour 4.4.1. Distribution of flow rate and flow velocity The pore-scale flow behaviour is generally similar for all the three samples with different hydraulic gradients. Therefore, in

7 6

dp/dx [MPa/m]

5 4

Experimental result of Sample 1 Experimental result of Sample 2 Experimental result of Sample 3 Numerical result of Sample 1 Numerical result of Sample 2 Numerical result of Sample 3

3 2 1 0 0

5

10

15

20 Q [ml/s]

25

30

35

40

Fig. 6. Comparison between the numerical results and the experimental results of the relationship between the pressure gradient and volumetric flow rate.

145

W. Zhang et al. / Computers and Geotechnics 87 (2017) 139–148

(a)

(b)

50

6 Re =5

Re =50

c

40

c

Re =10

5

c

Re =50 c

30

20

10

0 10

c

Re =250 dp/dx [MPa/m]

dp/dx [MPa/m]

Rec=25

Re =100

4

c

Rec=500

3 2 1

15

20

25

30

0 10

15

20

Q [ml/s]

25

30

Q [ml/s]

Fig. 7. Relationship between the pressure gradient and volumetric flow rate with different Rec .

Fig. 8. (a) Distribution of flow rate; (b) distribution of flow velocity on the aperture contour map and enlarged view.

100

100

2238 2000 1000

75

75

200

50

100

Y [mm]

Y [mm]

500

Flow path 1 Flow path 2 Flow path 3 Flow path 4

50

50 25

20

25

10 0

0

25

50

0

X [mm] Fig. 9. Distribution of Reynolds number.

0

0

25

50

X [mm] Fig. 10. Several representative flow paths.

hydraulic head difference and the specimen size along the flow direction. Based on these results, the hydraulic gradient is actually smaller than that calculated using the traditional method. On the other hand, the travel time of the flow paths is not constant. The travel time of flow paths 1–4 is respectively 0:168 s; 0:148 s; 0:260 s and 0:324 s. This finding agrees with the results of the tracer experiments conducted by Sawada and Sato [39]. Due to the inhomogeneous aperture distribution, the pore-scale flow behaviour through rough-walled fractures is quite

complex, and the locally developed turbulent flow further increases this complexity. 4.5. Evolution process with the increase of the flow rate Given the limitation of the current experiment condition, the pressure gradient utilized in water flow experiments cannot be

W. Zhang et al. / Computers and Geotechnics 87 (2017) 139–148

(a)

(b)

4

10

1

1 Ratio of turbulent elements Value of E

Regression result Numerical result 0.8

Value of E

dp/dx [MPa/m]

2

0.8

10

0

10

−2

0.6

0.6

0.4

0.4

0.2

0.2

10

0

0 −2 10

2

10

10

0

2

10

Q [ml/s]

Ratio of turbulent elements

146

0

10

Q [ml/s]

Fig. 11. (a) Relationship between the pressure gradient and volumetric flow rate obtained by the numerical simulation and the regression result based on the existing experimental data; (b) changing process of the ratio of turbulent elements and the value of E with the increase of the flow rate.

(b)

75

75

75

50

50 25

0

25

0

50

Y [mm]

75

Y [mm]

100

0

50 25

0

X [mm]

25

0

50

50 25

0

X [mm]

(e)

25

0

50

0

X [mm]

(f)

(g)

75

75

75

75

25 0

25

0

25

0

50

Y [mm]

100

Y [mm]

100

50

50 25

0

X [mm]

25

0

50

50 25

0

X [mm]

(i)

25

0

50

0

X [mm]

(j)

(k)

75

75

75

0

25

0

25

X [mm]

50

0

Y [mm]

75

Y [mm]

100

Y [mm]

100

25

50 25

0

25

50

X [mm]

0

50

(l)

100

50

25

X [mm]

100

50

50

(h)

100

50

25

X [mm]

100

Y [mm]

Y [mm]

(d)

100

25

Y [mm]

(c)

100

Y [mm]

Y [mm]

(a) 100

50 25

0

25

X [mm]

50

0

0

25

50

X [mm]

Fig. 12. Detailed evolution process of the turbulent elements with the increase of the flow rate: (a) Q ¼ 0:1 ml/s; (b) Q ¼ 0:2 ml/s; (c) Q ¼ 0:5 ml/s; (d) Q ¼ 1 ml/s; (e) Q ¼ 2 ml/s; (f) Q ¼ 5 ml/s; (g) Q ¼ 10 ml/s; (h) Q ¼ 20 ml/s; (i) Q ¼ 50 ml/s; (j) Q ¼ 100 ml/s; (k) Q ¼ 200 ml/s; (l) Q ¼ 500 ml/s.

very high. In this study, the numerical model is used to further study the non-Darcy fracture flow under high pressure gradient. Once again, Sample 1 is chosen in the numerical simulation. The

maximum pressure gradient is approximately 1400 MPa=m, and the corresponding flow rate is approximately 500 ml/s. The maximum water pressure in the numerical simulation is more than

W. Zhang et al. / Computers and Geotechnics 87 (2017) 139–148

200 times of that in the laboratory experiments. In terms of the laboratory experiments, because of distinct hydro-mechanical effect and potential leakage, it is impossible to conduct the water flow experiment under such high water pressure. Fig. 11(a) shows the relationship between the pressure gradient and volumetric flow rate obtained by the numerical simulation. The regression result based on the existing experimental data is also shown in this figure. It can be seen that the numerical results agree well with the regression result, indicating that both the numerical model and the regression model using Forchheimer equation can reasonably predict the non-Darcy fracture flow under high pressure gradient. It has been known that the effect of non-linear deviation becomes more significant under higher flow rate. To quantitatively estimate non-linear degree in fluid flow through rock fractures, a factor E is introduced to estimate the non-linear degree [40]:

E¼

bQ

2

aQ þ bQ

2

ð24Þ

where a and b are coefficients a and b in Forchheimer equation. According to the regression analysis based on the experimental data [13], a is taken to be 1:2  1010 Pa  s  m4 , and b is taken to be 6:7  1015 Pa  s2  m7 . To study the pore-scale flow regime evolution process, the concept of turbulent element is introduced. The turbulent elements are defined as the elements possessing a Reynolds number Re greater than the critical Reynolds number Rec . Fig. 11(b) shows the variation process of the ratio of turbulent elements and the value of E with the increase of the flow rate. It is clear that with the increase of the flow rate, the ratio of turbulent elements and the value of E increase synchronously, indicating that more turbulent elements appear and the non-linear degree gradually increases. Therefore, the increase of turbulent elements might be the pore-scale mechanism of the increase of non-linear deviation. A detailed evolution process of the turbulent elements with the increase of the flow rate is shown in Fig. 12. This figure clearly shows that the turbulent elements gradually occur with the increase of the flow rate. It should be mentioned that even the flow rate is rather low (for example, Fig. 12(c), Q ¼ 0:5 ml/s), turbulent flow still occurs in a few elements, indicating that although the volumetric flow rate is rather low and the macroscopic fracture flow is generally considered to be linear, local turbulent flow would occur due to the inhomogeneous flow distribution. 5. Conclusions In this paper, a pore-scale numerical model for non-Darcy fluid flow through rough-walled fractures is presented. Based on the assumption that the fracture aperture field follows a normal distribution, the fracture geometrical model is generated. Starting with the N-S equations, a simple general relationship between the local hydraulic conductivity and the flow velocity is proposed. A new governing equation for non-Darcy fluid flow through roughwalled fractures is then derived by introducing this relationship into the Reynolds equation. By Replacing the Darcy’s law in twodimensional finite element method for linear seepage analysis, the non-linear discrete equations are thus obtained. Eventually, Newton’s method is used to solve the non-linear discrete equations. The present numerical model is validated by comparing the numerical results with the results obtained from a recent experimental study. In the present model, there is only one more independent parameter, the critical Reynolds number Rec , when compared with the classic local cubic law. Numerical results show that, with the decrease of Rec , the deviation degree from linear flow gradually increases. By incorporating the experimental results into

147

the numerical model, it can be deduced that the local critical Reynolds number Rec might physically be close to 50. As the aperture distribution is inhomogeneous, the flow distribution is also inhomogeneous. Because the relationship between the local hydraulic conductivity and the flow velocity is based on the local cubic law, the flow distribution is more inhomogeneous than the aperture distribution. Thus, an obvious phenomenon of the preferential flow channels can be observed in the numerical simulations. Furthermore, the flow paths are generally tortuous, and thus the hydraulic gradient is actually smaller than that calculated with the traditional method. The numerical model is used to further study the non-Darcy fracture flow under high pressure gradient, which is generally difficult to be studied experimentally. It is found that both the numerical model and the regression model using Forchheimer equation can reasonably predict the non-Darcy fracture flow under high pressure gradient. Due to the inhomogeneous flow distribution, the turbulent flow occurs locally, even when the volumetric flow rate is rather small and the macroscopic flow is generally considered to be linear. In the local turbulent flow region, the local cubic law does not work and the present model is more suitable. With the increase of the flow rate, the local turbulent flow regions gradually become large and the flow resistances increase correspondingly. This leads to the gradual increase of non-linear deviation degree. Acknowledgements The research is supported by the Natural Science Foundation of China (NSFC) through Grant Nos. 41530638, 41372302 and High level talent project in Guangdong Province through Grant No. 20143900042010003. References [1] Witherspoon PA, Wang JSY, Iwai K, et al. Validity of cubic law for fluid flow in a deformable rock fracture. Water Resour Res 1980;16(6):1016–24. [2] Walsh JB. Effect of pore pressure and confining pressure on fracture permeability. Int J Rock Mech Min Sci 1981;18(5):429–35. [3] Brown SR. Fluid flow through rock joints: the effect of surface roughness. J Geophys Res–Atmos 1987;92(B2):1337–47. [4] Renshaw CE. On the relationship between mechanical and hydraulic apertures in rough-walled fractures. J Geophys Res–Atmos 1995;100(B12):24629–36. [5] Zimmerman RW, Bodvarsson GS. Hydraulic conductivity of rock fractures. Transp Porous Media 1996;23(1):1–30. [6] Brush DJ, Thomson NR. Fluid flow in synthetic rough-walled fractures: NavierStokes, Stokes, and local cubic law simulations. Water Resour Res 2003;39 (4):1037–41. [7] Lee SH, Lee KK, Yeo IW. Assessment of the validity of Stokes and Reynolds equations for fluid flow through a rough-walled fracture with flow imaging. Geophys Res Lett 2014;41(13):4578–85. [8] Quinn PM, Cherry JA, Parker BL. Quantification of non-Darcian flow observed during packer testing in fractured sedimentary rock. Water Resour Res 2011;47(9):178–87. [9] Cherubini C, Giasi CI, Pastore N. Bench scale laboratory tests to analyze nonlinear flow in fractured media. Hydrol Earth Syst Sci 2012;16(8):2511–22. [10] QIAN JZ, Wang M, Zhang Y, et al. Experimental study of the transition from non-Darcian to Darcy behavior for flow through a single fracture. J Hydrodyn 2015;27(5):679–88. [11] Kolditz O. Non-linear flow in fractured rock. Int J Numer Methods Heat Fluid Flow 2001;11(6):547–75. [12] Nowamooz A, Radilla G, Fourar M. Non-Darcian two-phase flow in a transparent replica of a rough-walled rock fracture. Water Resour Res 2009;45(7):4542–8. [13] Zhang Z, Nemcik J. Fluid flow regimes and nonlinear flow characteristics in deformable rock fractures. J Hydrol 2013;477(1):139–51. [14] Javadi M, Sharifzadeh M, Shahriar K, et al. Critical Reynolds number for nonlinear flow through rough-walled fractures: the role of shear processes. Water Resour Res 2014;50(2):1789–804. [15] Chen YF, Zhou JQ, Hu SH, et al. Evaluation of Forchheimer equation coefficients for non-Darcy flow in deformable rough-walled fractures. J Hydrol 2015;529:993–1006. [16] Singh KK, Singh DN, Ranjith PG. Laboratory simulation of flow through single fractured granite. Rock Mech Rock Eng 2015;48(3):987–1000.

148

W. Zhang et al. / Computers and Geotechnics 87 (2017) 139–148

[17] Zhang Z, Nemcik J, Ma S. Micro- and macro-behaviour of fluid flow through rock fractures: an experimental study. Hydrogeol J 2013;21(8):1717–29. [18] Auradou H, Drazer G, Boschan A, et al. Flow channeling in a single fracture induced by shear displacement. Geothermics 2006;35:576–88. [19] Jeong WC, Bruel D, Cho YS. Numerical experiments of flow and transport in a variable-aperture fracture subject to effective normal stresses. J Hydraul Res 2006;44(2):259–68. [20] Watanabe N, Hirano N, Tsuchiya N. Determination of aperture structure and fluid flow in a rock fracture by high-resolution numerical modeling on the basis of a flow-through experiment under confining pressure. Water Resour Res 2008;44(6):W06412. [21] Petchsingto T, Karpyn ZT. Deterministic modeling of fluid flow through a CTscanned fracture using computational fluid dynamics. Energy Sources Part A 2009;31(11):897–905. [22] Koyama T, Li B, Jiang Y, et al. Numerical modelling of fluid flow tests in a rock fracture with a special algorithm for contact areas. Comput Geotech 2009;36:291–303. [23] Crandall D, Bromhal G, Karpyn ZT. Numerical simulations examining the relationship between wall-roughness and fluid flow in rock fractures. Int J Rock Mech Min Sci 2010;47(5):784–96. [24] Pan PZ, Feng XT, Xu DP, et al. Modelling fluid flow through a single fracture with different contacts using cellular automata. Comput Geotech 2011;38 (8):959–69. [25] Liu R, Jiang Y, Li B, et al. A fractal model for characterizing fluid flow in fractured rock masses based on randomly distributed rock fracture networks. Comput Geotech 2015;65:45–55. [26] Pyrak-Nolte LJ, Nolte DD. Approaching a universal scaling relationship between fracture stiffness and fluid flow. Nat Commun 2016. http://dx.doi. org/10.1038/ncomms10663. [27] Javadi M, Sharifzadeh M, Shahriar K. A new geometrical model for non-linear fluid flow through rough fractures. J Hydrol 2010;389:18–30. [28] Trykozko A, Peszynska M, Dohnalik M. Modeling non-Darcy flows in realistic pore-scale proppant geometries. Comput Geotech 2016;71:352–60.

[29] Zimmerman RW, Kumar S, Bodvarsson GS. Lubrication theory analysis of the permeability of rough-walled fractures. Int J Rock Mech Min Sci 1991;28 (4):325–31. [30] Yeo IW, Ge S. Applicable range of the Reynolds equation for fluid flow in a rock fracture. Geosci J 2005;9(4):347–52. [31] Koyama T, Neretnieks I, Jing L. A numerical study on differences in using Navier-Stokes and Reynolds equations for modeling the fluid flow and particle transport in single rock fractures with shear. Int J Rock Mech Min Sci 2008;45 (7):1082–101. [32] Koyama T, Tsukahara T, Jing L, et al. Numerical simulation of laboratory coupled shear–flow tests for rock fractures: a comparative study on differences in using Reynolds and Navier-Stokes equations to simulate the fluid flow. In: Thermo-hydromechanical and chemical coupling in geomaterials and applications, proceedings of the third international symposium GeoProc; 2008b. p. 525–32. [33] Hakami E, Larsson E. Aperture measurements and flow experiments on a single natural fracture. Int J Rock Mech Min Sci 1996;33:395–404. [34] Keller AA. High resolution cat imaging of fractures in consolidated materials. Int J Rock Mech Min Sci 1997;34(3–4):358. [35] Fenton GA. Error evaluation of three random-field generators. J Eng Mech 1994;120(12):2478–97. [36] Liu L, Neretnieks I. Analysis of fluid flow and solute transport through a single fracture with variable apertures intersecting a canister: comparison between fractal and gaussian fractures. Phys Chem Earth 2006;31(10–14):634–9. [37] Barton N, de Quadros EF. Joint aperture and roughness in the prediction of flow and groutability of rock masses. Int J Rock Mech Min Sci 1997;34(3–4). 252– 252. [38] Tse R, Cruden DM. Estimating joint roughness coefficients. Int J Rock Mech Min Sci 1979;16(5):303–7. [39] Sawada A, Sato H. A study of hydraulic properties in a single fracture with inplane heterogeneity. Nucl Eng Technol 2010;42(1):9–16. [40] Zeng Z, Grigg R. A criterion for non-Darcy flow in porous media. Transp Porous Media 2006;63(1):57–69.