Seepage–stress coupled analysis on anisotropic characteristics of the fractured rock mass around roadway

Tunnelling and Underground Space Technology 43 (2014) 11–19

Contents lists available at ScienceDirect

Tunnelling and Underground Space Technology journal homepage: www.elsevier.com/locate/tust

Seepage–stress coupled analysis on anisotropic characteristics of the fractured rock mass around roadway T.H. Yang a,b, P. Jia a,b,⇑, W.H. Shi a,b, P.T. Wang a,b, H.L. Liu a,b, Q.L. Yu a,b a b

Key Laboratory of Ministry of Education on Safe Mining of Deep Metal Mines, Northeastern University, Shenyang 110819, China School of Resources & Civil Engineering, Northeastern University, Shenyang 110819, China

a r t i c l e

i n f o

Article history: Received 23 May 2013 Received in revised form 14 February 2014 Accepted 5 March 2014

Keywords: Fractured rock mass Anisotropy Discrete fracture network (DFN) model Seepage–stress coupled analysis Numerical simulation

a b s t r a c t Anisotropic properties of the fractured rock masses are investigated considering the coupled effect of the seepage and stress. The equivalent permeability and damage tensor of the fractured rock mass are initially examined using a series of Discrete-Fracture-Network (DFN) models with varied size and orientations from the geological investigation data of the sandstone roadway on the floor of 12# coal seam in Fangezhuang Coal Mine. A seepage–stress cross-coupling anisotropic model considering the coupled effect of the seepage and stress is described and applied to analyze the influence of the principal orientations of the joint sets on the anisotropic properties of the rock mass. It appears that the anisotropic properties of the rock mass have a great influence on the stress distribution, hydraulic conductivity coefficient and damage zone. The model may contribute to a more reasonable explanation on the dominant effect of the joint sets on deformation and failure of rock mass. Ó 2014 Elsevier Ltd. All rights reserved.

1. Introduction Discontinuities such as joints have a significant and complex effect on seepage and elastic properties of a rock mass, thus they can directly deteriorate the stability of roadways in coal mines. In situ discontinuities can be classified into five levels according to the structural features (Sun, 1991, 1997, 1998), in which class-IV and class-V are the most significant structural planes for rock mass stabilities. Numerical models based on the traditional rock mass classification methods like RMR, Q or GSI (Hoek, 2000) have been applied to assess the rock mass stabilities (Benardos and Kaliampakos, 2005), however, the distributions of the joint sets are usually neglected in these methods. This simplification may result in numerical inaccuracies especially when the rock mass properties are anisotropic. Thus, detailed characterization on rock mass discontinuities and appropriate involvement of these discontinuities in a numerical model may be the key to accurately simulate the mechanical behavior of the rock mass. When the equivalent continuum theory is used for the analysis of fractured rock mass, anisotropic properties must be considered due to the shear failure along these initial discontinuities (Gerrard, 1982; Fossum, 1985; Min and Jing, 2003). Anisotropic characteristics generally originate from the mineral foliation in metamorphic ⇑ Corresponding author. Tel.: +86 24 83687705. E-mail address: [email protected] (P. Jia). http://dx.doi.org/10.1016/j.tust.2014.03.005 0886-7798/Ó 2014 Elsevier Ltd. All rights reserved.

rocks, stratification in sedimentary rocks, and discontinuities in rock mass (Cho et al., 2012). Extensive efforts have been made to investigate the mechanical behavior of the transversely anisotropic rock material. Formulae used for characterizing the seepage tensor were presented to simulate different fracture networks by Snow (1968) by simulating different fracture networks; the influencing factors of the equivalent seepage tensor for a fractured rock mass were investigated by Witherspoon et al. (1980). Oda (1985) statistically processed the crack orientation data and analyzed the anisotropic properties of the fractured granite using stereographic projection. Based on the superposition principle of liquid dissipation energy, Zhou et al. (2008) proposed an analytical model to determine the permeability tensor of the fractured rock mass. The discrete fracture network (DFN) models were also wildly used in examining the equivalent permeability tensor (Baghbanan and Jing, 2007). Liu et al. (2009) and Sun et al. (2011) suggested that the anisotropic property should be applied to the permeability tests of the fractured rock mass. Liu et al. (2009) provided a multi-point-statistics method to predict the characteristics of the flow behavior in a fractured model. Sun and Zhao (2010) numerically studied the effects of the anisotropic permeability Sun et al. (2011) extended a similar research from 2 dimensional to 3 dimensional and they indicated that the anisotropic permeability has a significant impact on water pressure or water quantities. Yu et al. (2013) presented an analysis method for the anisotropic seepage field in the fractured rock mass using FLAC3D code, which

12

T.H. Yang et al. / Tunnelling and Underground Space Technology 43 (2014) 11–19

effectively simulated the seepage field of a water-sealed underground oil storage cavern. Many attempts had been made to determine the anisotropic properties of a rock mass. Mechanical properties such as elastic modulus and strength were investigated based on standard laboratory tests on transverse isotropic rock specimens (Nunes, 2002) or mica gneiss (Hakala et al., 2007; Gonzaga et al., 2008). The results showed some correlation between the mechanical properties and geological characteristics of the rock mass (Wen et al., 2000), which have provided some important references for further study of the rock mass anisotropic properties (Exadaktylos, 2001; Chen et al., 2008). For an engineering project, the deformation, strength and seepage characteristics of the fractured rock mass are critical parameters for design, construction and regular maintenance (Yang, 2009). However, the scale of a fractured rock specimen in the laboratory tests is insufficient to represent the equivalent continuum of the field rock mass; it is therefore necessary to develop a method that can be used to postulate the rock mass properties. The aim of this work is to develop a mathematical model to study the anisotropy of seepage and elasticity of a fractured rock mass. The DFN models are generated using Monte Carlo method, where the input parameters are derived from mapping on the engineering project related digital photogrammetry model using ShapMetrix3D. The permeability tensors and elasticity tensors of the rock mass are calculated based on the discrete medium seepage method and anisotropic damage theory respectively. By using the COMSOL Multiphysics numerical method, an anisotropy seepage model considering the influence of the discontinuity orientations, is built and applied to the seepage and stress filed analysis of the roadway on the floor of 12# coal seam in Fangezhuang Coal Mine. The influences of the principal direction of seepage or elasticity on stress, seepage and damage zone are discussed.

2. Constitute model of the anisotropic rock mass Determination of the constitute model of an anisotropic rock mass is very important for calculating the corresponding mechanical properties. Once the constitutive model is defined, a more efficient continuum numerical analyses based on the equivalent properties can be adopted. To describe the constitutive model, three premises should be confirmed initially (Wang et al., 2013): (i) Anisotropy of the rock mass is predominantly controlled by structural planes. (ii) Rock mass can be treated as an elastic, homogeneous and anisotropic object. (iii) Seepage tensor and damage tensor can be defined in the Representative Element Volume (REV) of the rock mass. 2.1. Stress

e11

3

2

3

2

1 E1



m231

E3 6 6   6 7 6 6 e22 7 ¼ ½S½r ¼ 6  m21 þ m231 4 5 E1 E3 6 4 e12 0

  m2  mE121 þ E313 1 E1

 0

m231 E3

0 0 2ð1þm12 Þ E1

3

2 3 7 r11 76 7 76 74 r22 7 5 7 5

r12

where [S] is the flexibility matrix of the material under the plane strain condition. The inverse matrix [S]1 (or [E]) is expressed by:

E21 ðE3 m21 þE1 m231 Þ

D E21 ðE1 m231 E3 Þ

D

0

0 0 E1 2ð1þm12 Þ

3

2 3 7 e11 76 7 74 e22 5 5

e12

ð2Þ where D is expressed as follow:

D ¼ E1 E3 þ E1 E3 v 221 þ E21 v 231 þ E21 v 231 þ 2E21 v 21 v 231

ð3Þ

The Young’s moduli E1 and E3 and Poisson’s ratio v31 and v12 can be determined by uniaxial compressive or tensile tests in the 1(or 2) and 3 directions. When pore water pressure is considered, the stress equilibrium equation is denoted by Eq. (4) (Terzaghi et al., 1996).

rij ¼ Dijkl ekl  aij Pdij

ð4Þ

where rij denotes the total stress tensors; Dijkl denotes the elastic tensors of the solid phase; ekl denotes the strain tensors; aij denotes a positive constant; P denotes the hydraulic pressure and dij denotes the Kronecher delta function. The strain tensors ekl are expressed by

ekl ¼ ðU k;l þ U l;k Þet =2

ð5Þ

where Uk,l and Ul,k are the displacements; ev is the volume strain. 2.2. Seepage Based on Biot’s theory (Biot, 1941), the relationship between the hydraulic pressure P and the hydraulic conductivity Kij is expressed by:

K ij r2 P ¼ 0

ð6Þ

where Kij is the hydraulic conductivity, P is the hydraulic pressure. The equation of hydraulic conductivity is expressed by Eq. (7) for a plane strain model. The coordinate transformation form is expressed by Eq. (8)

K ij ¼



K 11 K 21

K 12 K 22

 ð7Þ

8 K 0 þK 0 K 0 K 0 > > K ¼ 11 2 22 þ 11 2 22 cos 2h > < 11 K 0 K 0

K 12 ¼ K 21 ¼  11 2 22 sin 2h > > > : K ¼ K 011 þK 022  K 011 K 022 cos 2h 22 2 2

ð8Þ

where K 011 is the hydraulic conductivity coefficient along the strike of joint planes, K 022 is the hydraulic conductivity coefficient perpendicular to the strike of joint planes and h is measured from the x axis to the optimal direction of permeability. The ground water permeability is controlled by the normal stress on joint planes. The stress is coupled to seepage according to Eq. (9) by Louis (1974).

ð9Þ

where Kf is the current groundwater hydraulic conductivity; K0 is the initial hydraulic conductivity; r is the normal stress on the joint plane, and a is the coupling parameter that reflects the influence of stress. The principal hydraulic conductivities are recast by Eq. (10):

K 011 ¼ K 0110 ear2 ð1Þ

E2 ðE3 E1 m2 Þ

1 31 r11 D 6 6 E21 ðE3 m21 þE1 m2 Þ 6 7 31 4 r22 5 ¼ ½E½e ¼ 6 4 D r12 0

K f ¼ K 0 ear

For orthotropic material in plane problems, the elastic constitutive model in plane strain problems excluding the factor of pore water pressure is given by:

2

2

K 022 ¼ K 0220 ear1

) ð10Þ

where K 0110 and K 0220 are the initial hydraulic conductivities in the principal permeability directions; r2 and r1 are the normal stresses in the elastic principal directions. LU decomposition is used to calculate the equations with a convergence precision of 106.

T.H. Yang et al. / Tunnelling and Underground Space Technology 43 (2014) 11–19

13

2.3. Coordinate transformation Generally, joint planes are inclined at an angle from the principal stress direction, as shown in Fig. 1. The angle h is the inclination of joints measured from the x-axis to the 1-axis. Considering the influence of the hydraulic pressure, Eq. (4) is expressed by

8 9 > < rx > = > :

ry sxy

8 > < ex ¼ ½T r 1 ½Q ½T e  ey > > :c ;

8 9 >

=  aij P dij > > : > ; ; 0 xy 9 > =

ð11Þ

where [Tr]1 is the reverse matrix of stress coordinates transformation, [Te] is the strain coordinates transformation matrix, following:

2 ½T r 1

6 ¼4

2

cos2 h

sin h

2

cos2 h

sin h

2 cos h sin h

3

7 2 sin h cos h 5

ð12Þ

Fig. 2. Image captured in a sandstone roadway.

2

sin h cos h  sin h cos h cos2 h  sin h 2 6 ½T e  ¼ 4

2

cos2 h

sin h

2

cos2 h

sin h

2 cos h sin h

3

7  sin h cos h 5

ð13Þ

2

2 sin h cos h 2 sin h cos h cos2 h  sin h

Joint set: 2#

Joint set: 1#

3. Case study The described constitutive model is used to analyze the seepage and stress distribution of the roadway on the floor of 12# coal seam in Fangezhuang Coal Mine. The roadway is 530 m underground and the in situ water pressure is 1–2 MPa. The anisotropic properties of the seepage and stress field are discussed below. 3.1. Mechanical parameters Fig. 3. Distribution map of joint fissures.

Using a 3D contact-free ShapeMetrix3D system (Austria Startup company, 2008), the exposed discontinuities characterized with joint trace length, joint orientation and joint spacing in the field were mapped and shown in Figs. 2 and 3 respectively. The DFN model was then simulated using the Monte Carlo method (Fig. 4). A 16 m  16 m DFN model generated from the statistical data of the mapped joints is shown in Fig. 5. Five concentric statistical windows with varied size (3 m, 5 m, 8 m, 9 m and 10 m) are used and rotated every 15° anticlockwise to quantify the REV size. The permeability tensors were captured using the method proposed by Yang (2009). The hydraulic conductivity coefficients in different statistical windows are shown in Fig. 6. The window size is scaled from 3 m to 10 m with a different step length until the difference of the equivalent parameters in different directions becomes

1

2 y K22'

Fig. 4. The fracture network of the 3D roadway.

x ' 11

K

O

Joint plane

Fig. 1. Schematic diagram showing directions of joint planes, principal stress and permeability (K 011 is the coefficient of hydraulic conductivity along the strike of joint planes and K 022 is the coefficient of hydraulic conductivity perpendicular to the strike of joint planes. The angle 1Ox, i.e. h, describes the direction of joint planes from x-axis.).

minimal. As the size of the DFN model increases, the hydraulic conductivity coefficient decreases gradually. When the size increases to 10 m, the change of the principal values of hydraulic conductivity coefficient becomes very small compared with that of the 9 m window, the difference being 5.79%. The scale of the REV for the current discontinuity distribution is thus determined as 10 m  10 m. Moreover, the anisotropy of the seepage varies for a different discontinuity inclination. The principal direction angle

14

T.H. Yang et al. / Tunnelling and Underground Space Technology 43 (2014) 11–19

15m

0 345

3.5

330

15

10m

30

3

315

45

2.5 2

300

5m 60

1.5 1

285

75

0.5

270

90

0

255

105

240

120 225

135 210

150 195

165 180

Fig. 7. Damage tensor (101).

also 10 m and the initial damage tensor of REV can be obtained. The principal damage variable Dmin and Dmax for the fractured sample is 0.13 and 0.35 respectively and the principal direction angle h of the damage tensor is approximate 135°, perpendicular to the direction of the principal permeability. According to the principal of energy equivalence by Sidoroff (1981), the flexibility matrix for the fractured rock sample is denoted by

Fig. 5. Generated fracture network (16 m  16 m).

0 345

2

3m 15

330

30 1.5

315 300

5m 45 60

1

285

75

0.5

270

90

0

255

8m

S0ij ¼ ð1  Di Þ1 Si;j ð1  Dj Þ1

9m

where Sij is the flexibility matrix for intact rock; Di and Dj are the principal damage variables in i and j direction respectively. The related Young’s modulus and Poisson’s ratio are listed in Table 1 according to the principle by Zhang (Zhang, 2006).

10

ð15Þ

105

3.2. Numerical model 240

120 225

In order to simulate the anisotropic stress and seepage field of the fractured rock mass after excavation, a roadway model is built. The size of the model is 18 m  15 m with a u-shaped roadway of 3 m  3 m, as shown in Fig. 8.

135 210

150 165

195 180

Fig. 6. Permeability tensor (106 m/s).

is about 45°. The maximum and minimum hydraulic conductivity coefficients are 1.42 (106 m/s) and 6.67 (107 m/s) respectively, the maximum to minimum ratio being 2.13. The damage tensors are defined by Eq. (14) (Kawamoto et al., 1988) according to the geometric information of a fractured sample:

Dij ¼

N lX aðkÞ ðnðkÞ  nðkÞ Þ ði; j ¼ 1; 2; 3Þ V k¼1

ð14Þ

where N is the number of joints; l is the minimum spacing between joints; V is the volume of the rock mass; n(k) is the normal vector of the kth joint; and a(k) is the trace length of the kth joint for 2 dimensions. Similarly, the scale effect of the damage variables is studied for different statistical windows. The size of the fracture network is chosen as 5 m, 10 m and 15 m respectively. The principal damage variable for each statistic window is calculated every 15° anticlockwise. The damage tensors for these fracture networks are showed in Fig. 7. The damage variable tends to be same when the size of the DFN model increases from 10 m to 15 m, thus the REV size is

(1) Stress conditions: the left and right boundaries of the model are fixed in x direction and the bottom boundary is fixed in all three directions. A pressure of 11 MPa is applied on the top of the model to simulate the 523 m underground in situ load. (2) Permeability conditions: a hydrostatic pressure of 2 MPa (pw) is applied on the bottom boundary. Non-flow conditions are set on the other three boundaries. No initial water pressure is applied on the roadway boundary. 4. Results and discussion All the formulae described above are used by COMSOL Multiphysics, a PDE-based multiphysics modeling environment (COMSOL A.B., 2005). Fig. 9 shows the water pressure nephogram and

Table 1 Mechanical parameters of rock mass. Elastic modulus (MPa)

Poisson’s ratio

E1 = 2352 E2 = 1826 E3 = 2119

v12 = 0.30 v13 = 0.27 v23 = 0.33

15

T.H. Yang et al. / Tunnelling and Underground Space Technology 43 (2014) 11–19

The nephogram of the principal hydraulic conductivity coefficients K11 and K22 are plotted. In Fig. 10, the maximum value mainly concentrates at roof and floor of the roadway, as the inclination of joints is 45° in this simulation. The permeability factors at the two side walls are therefore lower compared with that at the floor and crown of the roadway, due to the effect of the normal stress on these joints. As a comparison, the distribution of K22 in Fig. 11 shows an opposite distribution from that in Fig. 10. To quantify the influence of the stress on hydraulic properties, a comparison between two simulations is presented here. The principal hydraulic conductivity coefficients K11 and K22 along A–A0 (showed in Fig. 11) line are calculated by a seepage–stress coupled model and a decoupled model respectively, as in Figs. 12 and 13. Results show that the values of K11 and K22 decrease when stress is considered during simulation.

16 14 12 10 8 3m

3m

6

I-I

4 2 0 0

2

4

6

8

10

12

14

16

18

Fig. 8. The numerical model and stress boundary conditions.

Maximum 8.5×107 Pa/m

7

10

8

5. Further discussion To further study the effect of the rock mass anisotropy on seepage field, stress field and evolution of damage, six models with joint inclination ranging from 0° to 150° are simulated. 5.1. Seepage

7 6

Fig. 14 shows the fluid pressure distribution and the flow direction (in red) in these models. Results suggest that the anisotropy of the rock mass has an obvious influence on the seepage field of the

5 4

×10-7

Maximum 7.64×10-7 m/s 3

7.5

2

7.0

1

6.5

Minimum 1.48×10 Pa/m

-6

Maximum 1.53×10 m/s

A’

A

Fig. 9. The pressure gradient nephogram.

6.0 5.5

×10-6

5.0

1.5 1.45

4.5 -7

1.4

Minimum 4.11×10 m/s Fig. 11. Nephogram of conductivity coefficient K22.

1.35

1.25 1.2 1.15 -6

Minimum 1.13×10 m/s Fig. 10. Nephogram of conductivity coefficient K11.

flow directions (in red1). The seepage field is apparently asymmetric due to the influence of the inclined joints. Fluid vectors distribute mainly along the principal direction of permeability.

Hydraulic condutivity coefficient K11 (m/s)

1.3 1.45 1.4 1.35 1.3 1.25

Decoupled model

1.2

Coupled model

1.15 1.1 1.05 1

0

2

4

6

8

X coordinate (m) 1 For interpretation of color in Figs. 9 and 14, the reader is referred to the web version of this article.

Fig. 12. The conductivity coefficient K11 along A–A0 line.

10

T.H. Yang et al. / Tunnelling and Underground Space Technology 43 (2014) 11–19

surrounding rock mass. Asymmetries of the fluid pressure as well as flow vectors are observed when the seepage principal axis is not consistent with the co-ordinate axis. Fig. 15 shows the X and Y flow velocities along I-I section. In Fig. 15(a), the X flow velocity changes from positive to negative at the X coordinate of 9 m along I-I line for models with joint inclination 0° and 90°, while for other models the corresponding X coordinate values differ with each other. Similarly, the Y flow velocity presents anisotropic characteristics in different models.

6.8 6.6 6.4 6.2 6 Decoupled model

5.8

Coupled model 5.6

5.2. Stress

5.4 5.2

0

2

4

6

8

Fig. 16 shows the shear stress nephograms of the models with different joint orientations. The shear stress distribution changes relative to the principal direction. The maximum shear stress in different models changes from 1.80 MPa to 4.56 MPa when the angle is 30° and 150°, respectively.

10

X coordinate (m) Fig. 13. The conductivity coefficient K22 along A–A0 line.

Maximum 7.65×105 Pa/m

Maximum 9.76×105 Pa/m

θ = 0°

θ = 30°

Maximum 1.09×106 Pa/m

θ= 60°

Minimum 695.13 Pa/m

Minimum 268.59Pa/m

Minimum 45.87Pa/m

Maximum 9.61×105 Pa/m

6

5

Maximum 1.08×10 Pa/m

Maximum 9.89×10 Pa/m

θ = 90°

θ = 120°

Minimum 63.76Pa/m

θ =150°

Minimum 694.71Pa/m

Minimum 277.27Pa/m

Fig. 14. The pressure gradient of different scenarios.

1.5

× 10 −6

6

X coordinate (m)

X flow velocity (m/s)

1 0.5 0 -0.5 7.5 -1 -1.5 -2 -2.5

8

8.5

9

9.5

0° 30° 60° 90° 120° 150°

10

10.5

Y flow velocity (m/s)

Hydraulic condutivity coefficient K22 (m/s)

16

0° 30° 60° 90° 120° 150°

5 4 3 2 1 0 7.5

-3

×10-6

8

8.5

9

9.5

X coordinate (m)

(a)

(b)

Fig. 15. Seepage velocities alone the horizontal line of 1 m below the roadway floor.

10

10.5

17

T.H. Yang et al. / Tunnelling and Underground Space Technology 43 (2014) 11–19

Maximum 2.02×106 Pa

Maximum 1.80×106 Pa

Maximum 3.33×106 Pa

=30°

=60°

Minimum -3.33×106 Pa

Minimum -4.52×106 Pa

Minimum -3.27×106 Pa

Maximum 2.71×106 Pa

Maximum 3.11×106 Pa

Maximum 4.56×106 Pa

=0°

=90°

=150°

=120°

Minimum -2.87×106 Pa

Minimum -1.76×106 Pa

Minimum -2.08×106 Pa Fig. 16. The shear stress nephogram of different scenarios.

0

Y displacement (mm)

-10

7.5

8

8.5

9

9.5

10

10.5

-20

discussed by other researchers (Cho et al., 2012; Gonzaga et al., 2008). In this section, the Hoffman anisotropic strength criterion is used to assess the damage zone in this numerical model as shown in Eq. (16).

r21 Xt Xc

-30 -40 0° 30° 60° 90° 120° 150°

-50 -60 -70

X coordinate (m) Fig. 17. The Y displacements on roadway floor.

Table 2 The strength parameters of rock mass. Tensile strength (MPa)

Compressive strength (MPa)

Shear strength (MPa)

Xt = 0.252 Yt = 0.526

Xc = 12.83 Yc = 23.12

S = 10.12



r1 r2 Xt Xc

þ

r22 YtYc

þ

Xc  Xt Y  Yt s2 r1 þ c r2 þ 122 ¼ 1 Xt Xc YtYc S

where the parameters Xt and Yt represent the tensile strength parallel and perpendicular to joint planes, respectively. Xc and Yc represent the compressive strength parallel and perpendicular to joint planes respectively. S is the shear strength of the rock mass along joint planes. r1 represents the normal stress along the principle direction of elasticity and r2 represents the normal stress perpendicular to the principle direction of elasticity. s12 represents the shear stress. Table 2 lists the mechanical parameters of the fractured rock mass by reducing the strength of rock samples (Liu et al., 2010; Yang et al., 2009). Damage zones of different models are shown in Fig. 18, in which the damage direction is perpendicular to the seepage principal direction.

The Y displacements along the floor of the roadway are showed in Fig. 17, where asymmetry is also captured by comparing simulation results of the six models, implying that the subsidence of the roadway floor is influenced by the anisotropy of the seepage. 5.3. Damage zone The anisotropy of the strength parameters such as tensile and compressive strength of the fractured rock mass has been

ð16Þ

Fig. 18. Damaged zones of different scenarios.

18

T.H. Yang et al. / Tunnelling and Underground Space Technology 43 (2014) 11–19

anisotropic property of seepage and stress in jointed rock masses. A further validation of the proposed model in site-specific rock engineering needs to be performed. Acknowledgements The work presented in this paper is financially supported by the National Basic Research Program of China (No. 2013CB227902), the Natural Science Foundation of China (Grant Nos. 51174045, 51034001, 41172265, 50934006 and 51304036), the basic scientific research funds of Ministry of Education of China (Nos. N090101001, N120601002), and the project supported by the Specialized Research Fund for the Doctoral Program of Higher Education of China (Grant No. 20120042120053). The authors are grateful for these supports. Moreover, the authors would like to thank Professor Jian Zhao and the reviewers for the valuable comments and helpful suggestions during the paper review process. References

Fig. 19. Comparison between the damage zone and failure mode of surrounding rock. Top is observed by Vietor et al. (2010) and bottom by Li (2013).

It is observed that the shape and size of the damage zones are significantly influenced by the orientation of the joint plane. Take model with joint inclination 90° as an example, the tensile strength in the direction perpendicular to the joint planes is much lower than that in model with joint inclination 0°. Damage zones mainly concentrate at roof and bottom of the roadway in the form of roof falling and floor heave. Similar failure mode is also noted in other models. The simulation results are in good agreement with results given by Vietor et al. (2010) and Li (2013). The potential failure mode of the anisotropic fractured rock mass around roadway is simulated realistically and is feasible to be applied in a rock engineering project. In addition, a correlation exists between the shear stress shown in Fig. 15 and the damage zone shown in Fig. 19, the mechanism however needs a further study. 6. Conclusions A seepage–stress cross-coupling anisotropic model to investigate the influence of joint orientation on the anisotropic characteristics of the stress and seepage is described in the present paper. The proposed model is numerically applied in the seepage and stress distribution analysis of a roadway in 12# coal seam in Fangezhuang Coal Mine. Numerical results showed that the seepage and stress filed and the shape and size of damage zones are greatly influenced by the anisotropy of the rock mass controlled by preexisting discontinuities. The seepage field is apparently asymmetric due to the influence of the inclined joints and the flow vectors distribute mainly along the principal direction of permeability. The maximum value of hydraulic conductivity coefficient mainly concentrates at roof and floor of the roadway at the inclination of joints of 45°. The permeability factors at the two side walls are therefore lower compared with that at the floor and crown of the roadway, due to the effect of the normal stress on these joints. A correlation is captured between the shear stress and damage zone and the simulation results are in good agreement with the in situ engineering practice. The work presented in this paper indicates that the proposed anisotropic model based on equivalent continuum mechanics is reasonable in stability analysis of the fractured rock mass. It should be noted that, the work reported in this paper is an initial effort on the influence of joint orientation on the

Austria Startup Company, 2008. ShapeMetriX3D System Manuals. Europe and China Co. LTD., Shengyang. Baghbanan, A., Jing, L., 2007. Hydraulic properties of fractured rock masses with correlated fracture length and aperture. Int. J. Rock Mech. Min. Sci. 44, 704–719. Benardos, A.G., Kaliampakos, D.C., 2005. Hydrocarbon storage in unlined rock caverns in Greek limestone. Tunn. Undergr. Space Technol. 20, 175–182. Biot, M.A., 1941. General theory of three-dimensional consolidation. J. Appl. Phys. 12, 155–164. Chen, W.Z., Yang, J.P., Zou, X.D., 2008. Research on macro mechanical parameters of fractured rock masses. Chinese J. Rock Mech. Eng. 27 (8), 1569–1575. Cho, J.W., Kim, H., Jeon, S., Min, K.B., 2012. Deformation and strength anisotropy of Asan gneiss, Boryeong shale, and Yeoncheon schist. Int. J. Rock Mech. Min. Sci. 50, 158–169. COMSOL, A.B., 2005. COMSOL Multiphysics version 3.2. Stockholm. Exadaktylos, G.E., 2001. On the constraints and relations of elastic constants of transversely isotropic geomaterials. Int. J. Rock Mech. Min. Sci. 38, 941–956. Fossum, A.F., 1985. Effective elastic properties for a randomly fractured rock mass. Int. J. Rock Mech. Min. Sci. Geomechanics. Abstract 22 (6) 467–470. Gerrard, C.M., 1982. Equivalent elastic moduli of a rock mass consisting of orthorhombic layers. Int. J. Rock Mech. Min. Sci. Geomechanics. Abstract 19 (1), 9–14. Gonzaga, G.G., Leite, M.H., Corthesy, R., 2008. Determination of anisotropic deformability parameters from a single standard rock specimen. Int. J. Rock Mech. Min. Sci. 45, 1420–1438. Hakala, M., Kuula, H., Hudson, J.A., 2007. Estimating the transversely isotropic elastic intact rock properties for in situ stress measurement data reduction: a case study of the Olkiluoto mica gneiss, Finland. Int. J. Rock Mech. Min. Sci. 44, 14–46. Hoek, E., 2000. Practice Rock Engineering. A.A. Balkema, Rotterdam. Kawamoto, T., Ichikawa, Y., Kyoya, T., 1988. Deformation and fracturing behaviour of discontinuous rock mass and damage mechanics theory. Int. J. Numer. Anal. Meth. Geomech. 12 (1), 1–30. Li, X.L., 2013. TIMODAZ: a successful international cooperation project to investigate the thermal impact on the EDZ around a radioactive waste disposal in clay host rocks. J. Rock Mech. Geotech. Eng. 5, 231–242. Liu, H.L., Yang, T.H., Yu, Q.L., Chen, S.K., Wei, C.H., 2010. Numerical analysis on the process of water inrush from the floor of seam 12 in Fangezhuang coal mine. Coal Geol. Explor. 38 (3), 27–31. Liu, X.Y., Zhang, C.Y., Liu, Q.S., Birkholze, J., 2009. Multiple-point statistical prediction on fracture networks at Yucca Mountain. Environ. Geol. 57, 1361– 1370. Louis, C., 1974. Rock hydraulics. Rock Mechanics, L. Muller, ed., Springer-Verlag, New York, pp. 287–299. Min, K.B., Jing, L., 2003. Numerical determination of the equivalent elastic compliance tensor for fractured rock masses using the distinct element method. Int. J. Rock Mech. Min. Sci. 40 (6), 795–816. Nunes, A.L., 2002. A new method for the determination of transverse isotropic orientation and associated elastic parameters for intact rock. Int. J. Rock Mech. Min Sci. 39, 257–273. Oda, M., 1985. Permeability tensor for discontinuous rock masses. Geotechnique 35 (4), 483–495. Sidoroff, F., 1981. Description of anisotropic damage application to elasticity. Proceedings of IUTAM Colloquium on Physical Nonlinearities in structural analysis. pp. 237–244. Snow, D.T., 1968. Rock fracture spacing, openings and porosities. Soil Mech. Found. Div. Proc. ASCE 94 (SMI), 73–91. Sun, G.Z., 1991. The progress of rock mechanics-rock mass structure mechanics. Chinese J. Rock Mech. Eng. 10 (2), 112–116. Sun, G.Z., 1998. Rock Mass Structure Mechanics. Science Press, Beijing, pp. 23–57.

T.H. Yang et al. / Tunnelling and Underground Space Technology 43 (2014) 11–19 Sun, J.P., Zhao, Z.Y., 2010. Effects of anisotropic permeability of fractured rock masses on underground oil storage caverns. Tunn. Undergr. Space Technol. 25, 629–637. Sun, J.P., Zhao, Z.Y., Zhang, Y., 2011. Determination of three dimensional hydraulic conductivities using a combined analytical/neural network model. Tunn. Undergr. Space Technol. 26 (2), 310–319. Sun, Y.K., 1997. Rock mass structural mechanism – a new advance in rock engineering. Geol. Mech. 5 (4), 292–294. Terzaghi, K., Peck, R.B., Mesri, G., 1996. Soil Mechanics in Engineering Practice. Wiley-Interscience. Vietor, T, Li, X.L., Chen, G.J., Verstricht, J, Fisch, H., Fierz, T., 2010. Small scale in situ tests: bore-hole experiments at HADES and Mont Terri rock laboratories. In: Deliverable 8, TIMODAZ project. Wang, P.T., Yang, T.H., Xu, T., Yu, Q.L., Liu, H.L., 2013. A model of anisotropic property of seepage and stress for jointed rock mass. J. Appl. Math., vol. 2013, Article ID 420536, 19 pages, 2013. DOI: http://dx.doi.org/10.1155/ 2013/420536. Wen, F., Maohong, Y.U., Tonglu, L.I., Jian-Bin, P., 2000. Failure pattern and numerical simulation of landslide stability of stratified rock. Chinese J. Rock Mech. Eng. 19, 983–986.

19

Witherspoon, P.A., Wang, J.S.Y., Iwai, K., Gale, J.E., 1980. Validity of cubic law for fluid-flow in a deformable rock fracture. Water Resour. Res. 16 (6), 1016–1024. Yang, J.P., 2009. Study of macro mechanical parameters of fractured rock mass. Wuhan institute of the Chinese Academy of Sciences Geotechnical Engineering, Wuhan, pp. 57–74. Yang, T.H., Yu, Q.L., Chen, S.K., 2009. Rock mass structure digital recognition and hydro-mechanical parameters characterization of sandstone in Fangezhuang coal mine. Chinese J. Rock Mech. Eng. 28 (12), 2482–2488. Yu, C., Deng, S.C., Li, H.B., Li, J.C., Xia, X., 2013. The anisotropic seepage analysis of water-sealed underground oil storage caverns. Tunn. Undergr. Space Technol. 38, 26–37. Zhang, G.K., 2006. Study on equivalent orthotropic mechanical parameters and yield criterion of jointed rock mass and its engineering application. Hohai University, Nanjing, pp. 20–48. Zhou, C.B., Sharmam, R.S, Chen, Y.F., Rong, G., 2008. Flow-stress coupled permeability tensor for fractured rock masses. Int. J. Numer. Anal. Method Geomechanics 32, 1289–1309.