Analysis model for the excavation damage zone in surrounding rock mass of circular tunnel

Tunnelling and Underground Space Technology 35 (2013) 78–88

Contents lists available at SciVerse ScienceDirect

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

Analysis model for the excavation damage zone in surrounding rock mass of circular tunnel Haiqing Yang ⇑, Da Huang, Xiuming Yang, Xiaoping Zhou School of Civil Engineering, Chongqing University, Chongqing 400045, China Key Laboratory of New Technology for Construction of Cities in Mountain Area, Ministry of Education, Chongqing 400045, China

a r t i c l e

i n f o

Article history: Received 17 April 2012 Received in revised form 17 September 2012 Accepted 26 December 2012 Available online 24 January 2013 Keywords: Excavation damage zone Rock mass Tunnel Microcracks Equivalent crack method

a b s t r a c t The time dependent viscoelastic model is extended so that excavation lining process of tunnel can be taken into account, with special attention on the analysis of the excavation damage zone around a circular tunnel under hydrostatic condition. Combining the sliding crack model and equivalent crack method, an analytical model for the evolution of fractured zone in viscoelastic surrounding rock mass is proposed. The whole construction process of tunnel is divided into three stages, in which the stress and displacement are considered separately. Then, the time dependent evolution of EDZ is implemented in a short computation algorithm, in which both extent and direction of fractured zones are determined. It is key important for the settlement of rock bolt. Through sensitivity analysis of parameters, the dependence of EDZ on time, excavation velocity, structure of heterogeneous rock mass, lining supporting system, initial geostress and the radius of tunnel are determined. Because of the time dependence of rock mass, the failure and deformation behavior continue increasing, even if the excavation process finished. In addition, the displacement after lining is set up is far more than the one at the moment of excavation just finished. Ó 2012 Elsevier Ltd. All rights reserved.

1. Introduction Excavation of an underground tunnel will alter the mechanical properties of the surrounding rock mass around the opening. Knowledge of the direction and extent of the excavation damage zone (EDZ) is important for the design and construction of deep underground tunnel. EDZ, which is also referred to as the ‘‘excavation disturbed zone’’, is the zone around an opening where local failure occurs in rock mass. The localized failure zone results from stress redistribution which induced by excavation and timedependent behavior of rock mass. EDZ could also form a permeable pathway of groundwater flow which would reduce the safety (concern) of the tunnel. It is therefore important to characterize the EDZ for the design and construction of underground rock engineering. It is deemed that there are mainly three aspects, which are excavation process, time dependent behavior of rock mass and the evolution of fractured zones in rock mass, may contribute to formation of the EDZ around a tunnel. To our general knowledge (Shen and Barton, 1997; Mitaim and Detournay, 2004), the formation of EDZ is sensitive to microcracks initiation and propagation.

⇑ Corresponding author at: School of Civil Engineering, Chongqing University, Chongqing 400045, China. Tel.: +86 23 65121982; fax: +86 23 65123511. E-mail address: [email protected] (H. Yang). 0886-7798/$ - see front matter Ó 2012 Elsevier Ltd. All rights reserved. http://dx.doi.org/10.1016/j.tust.2012.12.006

When a tunnel is excavated in a crystalline rock at great depth, the surrounding rock mass commonly break as a result of initiation and growth of microcracks. In the past, EDZ was studied by various means, particularly in the rock mechanic experiment field. For example, Meglis et al. (2001) measured ultrasonic wave velocities within a Min-By tunnel at Canada’s Underground Research Laboratory (URL) during excavation period. Martino and Chandler (2004) studied EDZ at the URL and measured the properties of the damage zone surrounding the opening. Some researchers tried to derive the theoretical model to analyze the damage evolution process. Lots of analytical solutions have been proposed to estimate the change of stress and deformation in the ground due to the excavation of a tunnel and to predict the stress and displacement distribution in the surrounding rock mass around tunnel. Perhaps the best known analytical solution is that formulated by Schwartz and Einstein (1978, 1980), which is based on the relative stiffness between the support and the surrounding ground. An elastic plane strain solution for stresses and displacements around a lined circular tunnel in an isotropic geomaterial due to uniform ground loads and internal pressure is obtained by Li and Wang (2008). Bobet (2009) deduced the complete analytical solutions for a shallow tunnel in saturated ground. Fahimifar et al. (2010) presented a closed-form solution in order to determine tunnel wall displacement and ground pressure imposed on tunnel support.

H. Yang et al. / Tunnelling and Underground Space Technology 35 (2013) 78–88

Unfortunately, we have to admit that in spite of their great importance, theoretical model cannot be effectively used to determine the failure behavior of EDZ, which might end up with tunnel collapse. In order to overcome this difficulty, some numerical models have been developed by numerous authors. Hommand-Etienne et al. (1998), for example, proposed a model for numerically analyzing the damage zone around an underground opening excavated in brittle rocks, and found that the highly damage zones are formed in association with the zones of high compressive stress. Zhu and Bruhns (2008) simulated excavation damaged zone around a circular opening under hydrostatic condition by using RFPA2D. In view of this, it can be concluded that: (1) for theoretical model, although the solution may be exact enough in theory, it is difficult to determine the extent and direction of EDZ and (2) while, as far as numerical models, the effect of excavation process and time dependent behavior of rock mass on development of EDZ are not easily taken into account by most of preexisting numerical models. In order to overcome this difficulty, theoretical and numerical methods are combined in this paper. Based on the stress and displacement distribution in surrounding rock mass around tunnel, construction effects of different excavation stage are discussed. Both radial and axial excavation velocity of tunnel is considered. Subsequently, the microcracking process due to local stress concentrations around pre-existing cracks are believed to be responsible for the evolution of EDZ. Through a short computation algorithm, the evolution process of EDZ is visualized. Finally, parameters study is carried on.

2. Analytical model for the EDZ in surrounding rock mass of circular tunnel 2.1. Definition of the problem In this paper, the problem of excavating a circular tunnel in a viscoelastic rock mass is considered. The following assumptions are made in this regard: (1) The initial state of far field stress is hydrostatic. So the expression of stress distribution will be simple and axisymmetric. Without this condition, the expression will be too complex to solve, and the anisotropic behavior should be taken into account. (2) The gravity effect at the top of the tunnel wall is neglected, because the tunnel is buried at great depth. The damage evolution is mainly affected by geostress. (3) The surrounding rock mass are supposed to be heterogeneous and incompressible. In order to simplify the viscoelastic constitutive model, the rock mass is supposed to be incompressible. Thus, the bulk modulus is supposed to be infinite, which can simplify the constitutive model. Before the tunnel is excavated, there exist many microcracks in the rock mass. For simplicity, faults are not considered in this model. If there are macroscopic faults, it is not easy to determine the stress distribution around tunnel. (4) The lining structure in this model mainly refers to the shotcrete. Moreover, the lining is approximately regarded as an elastic ring, and the inner radius is r1. (5) The excavation of tunnel is simplified as the enlarge process of tunnel radius. That is to say, at the stage of excavation, the radius increases evenly or unevenly as time increasing. (6) It is also supposed that the effects of microcracks on the macroscopic stress and displacement distribution can be ignored. For theoretical discussing, I think this assumption is acceptable in order to simplify the theoretical model.

79

Moreover, many scholars had adopted similar assumption in their papers, such as the works of Mitaim and Detournay (2004) and Mortazavi and Molladavoodi (2012). The sliding crack model adopted in this paper is only used for assessing the extent of excavation damage zone (EDZ). Although present model is not an exact solution, it is helpful for the assessment of excavation disturbed or damage zone. This is key important for the construction of tunnel, especially for rock ground. Fig. 1 depicts a circular tunnel of radius r0 excavated in an infinite isotropic rock mass. Regarding the above assumptions, the stress and deformation distribution around tunnel is symmetrical with respect to the axis of the tunnel. The hydrostatic far field stress is r0, and the mechanical response of rock is dependent on time t, which is start at the beginning of tunnel being excavated. Both polar and cartesian coordinates are originated at the center of tunnel section. In Fig. 1, t1 represents the finish of excavation, while the lining structure start to work from t = t2. In this model, the total construction process is divided into three stages. Firstly, when 0 < t < t1, the radius of tunnel grows from 0 to r0. Secondly, when t1 < t < t2, the radius of tunnel is fixed at r0, and no lining has been set up. Usually, this stage is also called delay phase between excavation and reinforcement process. It is obvious that the contact stresses on the tunnel surface at these two stages are equal to zero. Finally, when t > t2, the lining structure already act on the supporting system. Then, the contact stress p(t) between lining and surrounding rock mass, which varies with time, start to take effects on the deformation and failure behaviors of surrounding rock mass. Furthermore, the excavation velocity both in section and perpendicular to section (axis direction of tunnel) are taken into account. In the first stage, the radius of tunnel is time-dependent. In order to take the excavation velocity in section into account, the radius of tunnel is written as

 RðtÞ ¼

aðtÞ; ð0 6 t < t1 Þ r0 ;

ðt P t1 Þ

ð1Þ

where a(t), which represents the enlarge process of the radius of tunnel, can be arbitrary. It is up to the construction manner. For simplicity, the excavation velocity in section is supposed to be even, so a(t) = r0  t/t1 is adopted. t1 represents the finish of excavation. In fact, every tunnel has a length along axial direction. So, a tunnel cannot be formed instantaneously. Because of the tunnel face, the stress and displacement in tunnel face is dependent on the distance from the tunnel face. Indeed, the analysis of stresses and displacements near the tunnel face is a three-dimensional problem. However, some experience formula can be used to simplify it into 2D plane-strain problem. It is reported in the works of Liu and Du (2004) that there is exponential relationship between the initial geostress and the axial excavation velocity. It is also said that because of the axial excavation, the initial geostress should be modified. Based on the results of numerical simulation, the following equation is proposed in the works of Liu and Du (2004).

p0 ðtÞ ¼ r0 ð1  0:7emt Þ

ð2Þ

where, the spatial effective coefficients is m = 1.575V/Rs. V represents axial excavation velocity. Rs depends on the excavation plan. If the section of tunnel is formed as a whole, we have Rs = r0. Otherwise, if the section of tunnel is formed by using bench method, Rs is limited in the range from 0 and r0. It is clear that by using Eq. (2) we can transform the 3D problem into plane-strain one. Moreover, the axial excavation velocity is taken into account in the experience equation.

80

H. Yang et al. / Tunnelling and Underground Space Technology 35 (2013) 78–88

Fig. 1. Circular tunnel in a viscoelastic rock mass with initial distribution of microcracks.

Since the rock mass is supposed to be viscoelastic, the constitutive model should be time-dependent. In this paper, Boltzmann viscoelastic model is adopted. In fact, other constitutive models such as Kelvin and Maxwell model are also fit for this discussion. A schematic representation of Boltzmann viscoelastic model is shown in Fig. 2. Since the rock mass is assumed to be incompressible, the bulk modulus will be infinite. The constitutive is well known as follows.

GðtÞ ¼

G1 þG2 G22 G1 G2 e g t þ G2 þ G1 G1 þ G2

ð3aÞ

the stress components, which are functions of radius r only, do not depend on the angle h. This implies that the hoop and radial stresses are principal stresses and there is no shear stress. For viscoelastic solution of circular tunnel under hydrostatic condition, a lot of expressions have been obtained. Some of them are obtained according to the correspondence principle. Others are derived by the Laplace transformation method. In this paper, a series of explicit solution, which is published in the works of Wang and Zhong (2012), is adopted. The detail expressions, which are derived by Laplace transformation method, are as follows. When 0 6 t < t1 , we have

ð3bÞ

rr ¼ r0 ð1  0:7emt Þ½1  R2 ðtÞ=r2 

ð4Þ

where G1 and G2 are shear modulus, while g is viscosity of Boltzmann model. These parameters can be easily obtained indoor creep test.

rh ¼ r0 ð1  0:7emt Þ½1 þ R2 ðtÞ=r2 

ð5Þ

2.2. Macroscopic stress field of surrounding rock mass at different stage

ur ðr; tÞ ¼

KðtÞ ¼ 1

According to the construction process of tunnel, the stress field is calculated separately. Since the far field stress is assumed to be hydrostatic, the stress distribution is axisymmetrical. Furthermore,

The radial displacement can be expressed as follows

ð6Þ

When t1 6 t < t2 , the stress field is written as

rr ¼ r0 ð1  0:7emt Þ½1  r20 =r2 

ð7Þ

rh ¼ r0 ð1  0:7emt Þ½1 þ r20 =r2 

ð8Þ

ur ðr; tÞ ¼

Fig. 2. Boltzmann viscoelastic model.

1 2 R ðtÞr0 ð1  0:7emt Þ 2G2 r Z G1 1 Gg1 t t 2 þ R ðsÞr0 ð1  0:7ems Þe g s ds e 2gr 0

1 2 r r0 ð1  0:7emt Þ 2G2 r 0 Z G1 1 Gg1 t t 2 þ r 0 r0 ð1  0:7ems Þe g s ds e 2gr 0

ð9Þ

H. Yang et al. / Tunnelling and Underground Space Technology 35 (2013) 78–88

wing cracks will growth from the tip of initial microcracks, in an angle with respect to the initial microcracks. The wing crack is also assumed to be open. As mentioned above, the disturbance of macroscopic stress field, caused by microcracks, is ignored. Therefore, the shear and normal stresses on the initial crack plane can be determined. They can be calculated by

When t P t2 , the stress field is determined by

rr ¼ r0 ð1  0:7emt Þ þ

r20 r2

(

r0 ð1  0:7emt Þ

0:7r0 r 0 G1  mg 0:7r0 r 0 G1 emt  emt2 2AG2 G1  kg  mg 2AG2 G1  kg  G 0:7r0 r 0 mg2 1 k emt2 e k g ðtt2 Þ 2AG2 ðG1  kg  mgÞðG1  kgÞ    ) G B G1 r0 Gg1 t2 k g1 ðtt 2 Þ  e e 1 2A G1  kg þ

rh ¼ r0 ð1  0:7emt Þ þ

r20 r2

ð10Þ

r0 ð1  0:7emt Þ

G1

r0 where A ¼ 2G þ 2G1 c 2 G1

r 0 r 21

r 20 r21

rn ¼ rxx sin2 a  2rxy sin a cos a þ ryy cos2 a

ð15bÞ

r0s ¼ ðrxx  ryy Þ sin b cos b þ rxy ðcos2 b  sin2 bÞ

ð16aÞ

r0n ¼ rxx sin2 b  2rxy sin b cos b þ ryy cos2 b

ð16bÞ

ð11Þ 3. Evolution of excavation damage zone in the rock mass

pðtÞ ¼

 p0 ðsÞe g s ds

ð15aÞ

where b is the angle between wing crack and horizontal line. It will be determined in the following section.

In this stage, the lining system starts to take effects on the stress and displacement distribution around tunnel. It is obvious that the contact stress between lining and rock mass is time-dependent. Moreover, the contact stress is determined by the relative displacement between lining and the surrounding rock mass. On the basis of the compatibility conditions of deformation, the contact stress between lining and rock mass is obtained by Wang and Zhong (2012).

0:7r0 r 0 G1  mg 0:7r0 r 0 G1 emt  emt2 2AG2 G1  kg  mg 2AG2 G1  kg  G 0:7r0 r 0 mg2 1 k emt2 e k g ðtt2 Þ 2AG2 ðG1  kg  mgÞðG1  kgÞ     G1 B G1 r0 Gg1 t2 þ e e k g ðtt2 Þ  1 2A G1  kg Z r2 1 Gg1 t t 2 ur ðr; tÞ ¼ 0 ½pðtÞ  p0 ðtÞ þ e R ðsÞ½pðsÞ 2G2 r 2gr 0

rs ¼ ðrxx  ryy Þ sin a cos a þ rxy ðcos2 a  sin2 aÞ

where a is the dip angle of microcracks. Similarly, the shear and normal stresses on the wing crack plane. They can be calculated by

(

0:7r0 r 0 G1  mg 0:7r0 r 0 G1 emt  emt2 þ 2AG2 G1  kg  mg 2AG2 G1  kg  G 0:7r0 r 0 mg2 1 k emt2 e k g ðtt2 Þ 2AG2 ðG1  kg  mgÞðG1  kgÞ    ) G B G1 r0 Gg1 t2 k g1 ðtt 2 Þ þ e e 1 2A G1  kg

81

ð12Þ

In order to model the evolution of excavation damage zone in the surrounding rock mass, the equivalent crack method is combined into the sliding crack model here. The evolution process is implemented in a numerical iteration. The basic principle is summarized as follows. Firstly, the SIF at the tip of preexisting microcracks is obtained by using sliding crack model. Meanwhile, the length and dip angle of wing crack is determined, as shown in Fig. 3a. Subsequently, the crack tips of wing crack are connected with a straight line, with a length of 2an+1 as depicted in Fig. 3b. Then, this new line is defined as a virtual straight preexisting microcrack. This is the core process of equivalent model. So, through this equivalence, a kinked crack problem is simplified as a series of straight crack problems. The detailed of the calculation process is detailed as follows. Although interaction of cracks is ignored here, this model also shows good agreement with engineering practices. If the interaction of cracks is taken into account, the problem would become

ð13Þ þ 1þK lc c

G1

r30

r0 Ec Ec ; k ¼  2A g; Gc ¼ 2ð1þlc Þ; K c ¼ 12lc ; G1 R t1 ð1  0:7ems ÞR2 ðsÞe g s ds; Gc rep0

r20 r 21

0 eð g mÞt1 þ g1r0 B ¼  Gr01 e g t1 þ G0:7r 1 mg

resents the shear modulus of lining, while Kc is the bulk strain modulus of lining structure. Eqs. (4)–(11) can be change into x–y coordinate system. The elastic stresses in the x–y coordinate are

rxx ¼ rr cos2 h þ rh sin2 h

ð14aÞ

rxy ¼ ðrr  rh Þ sin h cos h

ð14bÞ

ryy ¼ rr sin2 h þ rh cos2 h

ð14cÞ

Of all the possible sources of microcracks evolution, the mechanism of frictional sliding on pre-existing cracks is generally accepted as the primary stress concentrator that produces tensile wing crack growth after the excavation of tunnel. In this model, the failure of rock mass is induced by the propagation of microcracks. In order to analyze the growth rules, the stress distribution on the crack plane should be determined previously. The sliding crack model (Paliwal and Ramesh, 2008; Zhou and Yang, 2007), which is widely used in fracture mechanics, is adopted here. The initial microcracks are supposed to be straight. After excavation,

(a) schematic of kinked crack (solid line)

(b) Equivalent crack model for the iterative process Fig. 3. Computational model for the evolution of microcracks, Francois and Dascalu (2010).

82

H. Yang et al. / Tunnelling and Underground Space Technology 35 (2013) 78–88

3.1. Stress intensity factor and fracture criterion

70 r(t=3days) r(t=7days) r(t=20days) (t=3days) (t=7days) (t=20days)

Stress (MPa)

60 50 40 30 20 10 0

3

5

7

9

11

13

15

17

r (m) Fig. 4. Radial and hoop stress distribution around the tunnel perimeter at different time.

We now evaluate the stress intensity factor (SIF) at the crack tips. Owing to the excavation of tunnel, the stress distribution will disturbed. The change of stress field results in tractions acting on the preexisting microcracks and the wing-crack surfaces. These traction vectors have two components which are normal and parallel to the respective plane. trs represents tractions on the plane corresponding to the closed sliding crack, while, trw represents tractions on the plane corresponding to the straight open wing crack Hence, at the closed sliding crack location, we have

( t rs

ðt rs Þk ¼ rs ðt rs Þ? ¼ rn

ð17Þ

Radial displacement (m)

Similarly, at the straight open wing-crack location,

(

0.12

t rw

t=3days t=7days t=20days

0.1 0.08

ðt rw Þk ¼ r0s ðt rw Þ? ¼ r0n

ð18Þ

where the superscripts ‘‘k’’ and ‘‘\’’ indicates the directions parallel and perpendicular to the surfaces, respectively. We can write an expression for the tractions acting on the sliding flaw surfaces, as

0.06 0.04 0.02

(

0 3

5

7

9

11

13

15

17

r (m) Fig. 5. Radial displacement distribution around the tunnel perimeter at different time.

too complex to be solved. It is also not indeed needed in engineering application.

ts

ðt s Þ? ¼ ðt rs Þ? ðt s Þk ¼ sc  lðt rs Þ?

ð19Þ

where, l is the coefficient of dry friction for the closed cracks and sc is the cohesion.

KI ¼

pffiffiffiffiffi 2sðseff Þk cos a pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi  ðseff Þ? pl  pðl þ l Þ

Fig. 6. Evolution of fractured zones in the surrounding rock mass.

ð20aÞ

83

H. Yang et al. / Tunnelling and Underground Space Technology 35 (2013) 78–88 70

60

60 r(t1=1day) r(t1=3day) r(t1=5day) (t1=1day) (t1=3day) (t1=5day)

40 30 20

Stress (MPa)

Stress (MPa)

50

30

r(t2=15days)

20

r(t2=24days)

10 -10

0

5

10

15

20

25

30

35

r(t2=6days)

40

(t2=6days)

0

10 0

50

40

0

5

10

15

45

20

25

30

35

40

45

40

45

t (days)

t (days) Fig. 12. Stress versus time for different delay time.

0.3

Radial displacement (m)

Radial displacement (m)

Fig. 7. Stress versus time for different excavation velocity (r = 1.1r0).

0.25 0.2 0.15 0.1

t1=1day t1=3day t1=5day

0.05 0

0

5

10

15

20

25

30

35

40

45

0.16 0.14 0.12 0.1 t2=6days t2=15days t2=24days

0.08 0.06 0.04 0.02 0

0

5

10

15

t (days)

70

35

50 40 30 20

Fig. 13. Radial displacement versus time for different delay time.

60

Stress (MPa)

Stress (MPa)

30

70

r(V=0.5m/day) r(V=1.5m/day) r(V=3.0m/day) (V=0.5m/day) (V=1.5m/day) (V=3.0m/day)

60

10 0

50 40 30 20 10 0

-10 0

5

10

15

20

25

30

35

40

-10

45

0

5

10

15

t (days)

Radial displacement (m)

0.16 0.12 0.08

V=0.5m/day V=1.5m/day V=3.0m/day

0.04 0 5

10

15

20

25

30

35

40

45

0.06 0.05 0.04 0.03 0.02 0.01 1

1.5

40

45

0.12 0.1 t=12days t=20days t=30days

0.08 0.06

0

0.2

0.4

0.6

0.8

1

1.2

2

2.5

3

Fig. 15. Radial displacement versus thickness of lining for different times.

ðseff Þk ¼ sc  lðtrs Þ?  ðt rs Þk

ð20bÞ

ðseff Þ? ¼ r0n

ð20cÞ

Substituting Eqs. (17) and (18) into Eq. (20), the SIF at crack tip can be determined as follows

t=3days t=7days t=9days

0.5

35

r1-r0(m)

Fig. 10. Radial displacement versus time for different axial excavation velocity.

0

30

0.14

t (days)

0

25

Fig. 14. Stress versus time for different thickness of lining.

0.2

0

20

t (days)

Fig. 9. Stress versus time for different axial excavation velocity (r = r0).

Radial displacement (m)

25

t (days)

Fig. 8. Radial displacement versus time for different excavation velocity.

Radial displacement (m)

20

pffiffiffiffiffi 2sseff cos a K I ¼ pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi þ r0n pl pðl þ 0:27sÞ

ð21Þ

pffiffiffiffiffi 2sseff sin a K II ¼ pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi  r0s pl pðl þ 0:27sÞ

ð22Þ

3.5

V (m/day) Fig. 11. Radial displacement versus axial excavation velocity for different time.

84

H. Yang et al. / Tunnelling and Underground Space Technology 35 (2013) 78–88

where seff = sc  lrn  rs. The following fracture criterion is used to distinct the growth of crack tip

K I ¼ K ICC

ð23Þ

anþ1 ¼

sinð/max Þ dl sinðdhÞ

ð27Þ

where an and an+1 are the length of the equivalent straight crack ate steps n and n + 1.

where KICC is fracture toughness of intact viscoelastic rock. 3.3. Computation algorithm

3.2. Determine the growth trajectory of cracks According to the research of Schutte and Bruhns (2002), this kinking angle can be expressed as the following function: 2

/ ¼ sgnðK II Þ½0:70966k31  0:097725 sin ð3:9174k1 Þ  13:1588 tanhð0:15199k1 Þ

ð24Þ

where sgn is the signum function and k1 is a mode mixity factor that combines the stress intensity factors of mode I, KI, and mode II, KII, of the straight crack

k1 ¼

jK II j K I þ jK II j

ð25Þ

The equivalent crack is determined as follows for the crack orientation after each step of calculation by means of geometrical relationships.

sinð/max Þ þ cosð/max Þ dl

tanðdhÞ ¼ an

ð26Þ

For the updated crack length:

Due to the dependency of the failure behavior on the stress distribution in the surrounding rock mass, the problem is highly nonlinear. In this paper, the computation process is carried on step by step. For the whole process, the problem is solved by an iterative procedure as follows: Step 1. Parameters initialization (t = 0, n = 1). Step 2. According the construction stage of tunnel, calculate the macroscopic stress and displacement field. Step 3. Discrete the zone of surrounding rock mass into lots of small equivalent stress cell. In each position of cell, the stress is set to be equal to the stress at center. Step 4. In each cell, one initial microcrack is defined at the center. Then, update the length and dip angle of cracks. Step 5. Determination of the stress intensity factors. Step 6. Update of the crack length. Step 7. Update the growth trajectory of crack. Step 8. If the SIF meets the growth criterion, set s = an, go back to step 4.

Fig. 16. Evolution of fractured zones in the surrounding rock mass for different dip angle.

H. Yang et al. / Tunnelling and Underground Space Technology 35 (2013) 78–88

Step 9. Otherwise, set t = t + Dt, goback to step 1, keep on calculating till the time step finished. Then, enter post processing subroutine. where n is the iteration step, t is the time step number, Dt is the size of the time step. 4. Numerical results The damage evolution model described in Section 2 has been implemented in a simple MATLAB code. First, an investigation devoted to illustrate the main features of the proposed model is presented, then, sensing analysis of parameters are performed. The following parameters are considered: r0 = 3 m, r1 = 2.85 m, r0 = 30 MPa, g = 20 GPa d, G1 = 1 GPa, G2 = 2 GPa, Ec = 24 GPa, lc = 0.2, t1 = 5 d, t2 = 8 d, V = 1.5 m/d, s = 2 mm, sc = 0.3 MPa, a = 25°, E = 25 GPa, KICC = 0.5 MPa m1/2. In present model, the excavation is simplified as two-dimensional model. In order to take the spatial effects of excavation into account, an approximate equation, which can transfer the spatial effects to plane model, is adopted here. Under hydrostatic condition, the stress and displacement distribution is axisymmetric. The stress and radial displacement distribution is shown as follows. The radial and hoop stress distribution around the tunnel perimeter for different time is presented in Fig. 4. It is shown in Fig. 4 that the disturbance zone of stress distribution is extended to the distance which is nearly two times of the radius.

85

Furthermore, after lining system is set up, the difference between radial and hoop stress in the surrounding rock mass decreased as the time increasing. It can be concluded that the hoop stress will increase after excavation Fig. 5 shows the radial displacement distribution around the tunnel perimeter at different time. It is indicated in Fig. 5 that radial displacement is less than 0.12 m, which is 4% of the radius of tunnel. This value shows good agreement with engineering practices. In addition, the disturbed region of radial displacement is limited in 3 m from the tunnel boundary. The evolution process of fractured zones in the surrounding rock mass is presented in Fig. 6. Owing to the heterogeneous of rock mass, the microdefects will propagate as the rock mass experience unloading process. As shown in Fig. 6, at the moment of the excavation just finished, the fractured zone has not connected into a whole. However, the fractured zones grow as the time increasing, despite the lining system is set up. It is also revealed from Fig. 6e that the fractured zone will localized in a small extent. That means some cracks coalescence in local part, while other cracks, although have grown in a period, will stop propagating at some length. Obviously, from the crack growth pattern, it is indicated that the localized macro cracks can be controlled by rock bolt. 4.1. Influence of excavation velocity on the stress and displacement distribution in surrounding rock mass In present model, the excavation time in the tunnel section is defined as t1. If r0 is constant, a longer t1 represents a slower

Fig. 17. Evolution of fractured zones in the surrounding rock mass for different initial length of cracks.

86

H. Yang et al. / Tunnelling and Underground Space Technology 35 (2013) 78–88

excavation velocity. Owing to the creep behavior of surrounding rock mass, radial and hoop stress will vary with time. Fig. 7 depicts radial and hoop stress at r = 1.1r0 versus time for different excavation velocity. It is seen from Fig. 7 that the difference between radial and hoop stress increases as the excavation velocity increasing, especially in the stage before lining has been set up. As we all know, the difference between radial and hoop stress plays key role in the failure extent of rock mass. In view of the mechanism of stress redistribution in rock mass, it can be concluded that a faster excavation velocity leads to a fiercer disturbance on the stress distribution of surrounding rock mass. The dependence of radial displacement on time for different excavation velocity is depicted in Fig. 8. As illustrated in Fig. 8, after excavation, radial displacement grows as the time increasing. However, 40 days after excavation, the radial displacement gradually tends to be constant. Meanwhile, the final value of radial displacement remarkably increases with an increase of excavation velocity. It is well known that the stress distribution around tunnel is also affected by the axial excavation process, especially in the section which is close to the tunnel face. In present model, the axial excavation velocity is represented by V. Fig. 9 shows radial and hoop stress at r = r0 versus time for different axial excavation velocity. It is inferred from Fig. 9 that when the axial excavation velocity is faster than 1.5 m/day, the difference between radial and hoop stress increases with an increase of the excavation velocity, especially in the stage after lining is been set up. In order to control the stress disturbance, the axial excavation velocity should be limited.

The relationship between radial displacement and time for different axial excavation velocity V is presented in Fig. 10. It is obvious that the radial displacement increases with an increase of axial excavation velocity V. Moreover, as shown in Fig. 10, we can see that the radial displacement increases rapidly before ling is set up, while the growth rate of radial displacement is slowdown after the supporting system working. The final value of displacement increases as the axial excavation velocity V increasing. Fig. 11 depicts the radial displacement versus axial excavation velocity at different time. As shown in Fig. 11, it can be seen that the radial displacement grow as the axial excavation velocity increasing. However, when the axial excavation velocity is faster than 1.5 m/day, the growth rate of radial displacement is slowdown. It can be concluded from Figs. 7–11 that the excavation velocity takes great effects on the stress and displacement. A reasonable excavation velocity is less than 1.5 m/day. In engineering practice, this value is also recommended by many codes and engineers. Thus, this model can be used to optimize the excavation velocity in engineering practice.

4.2. Influence of lining on the stress and displacement distribution in surrounding rock mass In engineering practice, lining system will not be set up immediately after the excavation finished. The lining is always delayed for a short period, especially in the NATM (New Austrian Tunneling Method) construction. In this section, we will discuss the dependence of stress and radial displacement on the delay time of lining

Fig. 18. Evolution of fractured zones in the surrounding rock mass for different initial geostress.

H. Yang et al. / Tunnelling and Underground Space Technology 35 (2013) 78–88

construction. Moreover, the relationship between the thickness of lining structure and stress distribution is investigated here. Fig. 12 shows radial and hoop stress at r = r0 versus time for different delay time. It is inferred from Fig. 12 that the difference between radial and hoop stress grow with the delay time increases. That is to say if we delay a long time the loads bear by lining system would remarkably increase. This is unsafe for lining system. So, the delay time should be controlled in a reasonable extent. The dependence of radial displacement on time for different delay time is depicted in Fig. 13. As illustrated in Fig. 13, after excavation, radial displacement grows as the time increasing. Furthermore, the final radial displacement increases with an increase of delay time. If the lining system is delayed too long, the deformation of rock mass will sharply increase, which may cause the unstability of surrounding rock mass. Fig. 14 depicts the radial and hoop stress versus time for different thickness of lining. As shown in Fig. 14, it can be seen that the difference between radial and hoop stress reduced as the thickness of lining increasing. But, as the thickness of lining increasing, the costs of supporting system will grow evidently. The relationship between radial displacement and the thickness of lining for different time is presented in Fig. 15. In Fig. 15, (r1  r0) represents the thickness of lining system. It is obvious that the radial displacement decreases with an increase of the thickness of lining system. Moreover, as far as this example, the reasonable thickness would be less than 0.3 m. In engineering practice, it is very important to select the time of setting up the lining system. In present model, we can choose the time through the stress and displacement distribution, because the time-dependent behavior of rock is taken into account.

87

4.3. Influence of rock structure on the failure behaviors of surrounding rock mass In this section, we will discuss the influence of rock structure on the failure behaviors of surrounding rock mass, such as dip angle and initial length of cracks. The evolution process of fractured zones in the surrounding rock mass is presented in Fig. 16. As shown in Fig. 16, the location of fractured zones varies with the dip angle. It is one of the innovations of present model, because the fractured zones are visualized dynamically by using a simple method. As far as classical elastoplastic theoretical model, only the extent of plastic zone can be determined, but the growth direction of cracks cannot be obtained. However, it is revealed in Fig. 16 that the fractured zones will localized in some location, instead of homogeneous distribution. This is very important for the construction of rock bolt. Fig. 17 shows the evolution process of fractured zones in the surrounding rock mass for different initial length of cracks. As shown in Fig. 17, the extent of fractured zones increases with the initial length of cracks increasing. Moreover, it is revealed in Fig. 17 that the fractured zones will localized in some location, which is a little distance from the tunnel profile instead of on the surface of tunnel. There are many microcracks in the rock mass. The dip angle and initial length of cracks have significant influence on the damage evolution of rock mass. If we have detail geologic files, we can choose a reasonable routine which the extent of damage evolution is the minimum. Moreover, because the rock structure may vary in the excavation period, we can predict the damage evolution according the dip angle and initial length of cracks.

Fig. 19. Evolution of fractured zones in the surrounding rock mass for different radius of tunnel.

88

H. Yang et al. / Tunnelling and Underground Space Technology 35 (2013) 78–88

4.4. Influence of initial geostress and radius of tunnel on the failure behaviors of surrounding rock mass In this section, we will discuss the influence of geostress and radius of tunnel on the failure behaviors of surrounding rock mass. Fig. 18 depicts the evolution process of fractured zones in the surrounding rock mass for different initial geostress. It can be seen from Fig. 18 that the extent of fractured zones grows with an increase of geostress. Fig. 19 presents the evolution process of fractured zones in the surrounding rock mass for different radius of tunnel. It can be concluded from Fig. 19 that the extent of fractured zones increase as the radius of tunnel increasing. In addition, it is seen that the extent of fractured zones is nearly equal to the radius of tunnel. The initial geostress and radius of tunnel take great effects on the failure behaviors of surrounding rock mass. Before excavation, we should test the geostress, otherwise it may cause serious accidents. 5. Conclusions Although the EDZ problem has been studied in many literatures, the information is mainly limited to the locations where the surrounding rock mass is disturbed. Elasto-plastic theory is widely accepted in the previous study. However, as far as the crackweakened rock mass sensitive to the rock structure is considered, (which is sensitive to the rock structure,) the growth process of cracks should be taken into account. In this paper, an analysis model based on damage evolution mechanism is developed. In order to consider the time dependence of displacement, the viscoelastic rock mass is supposed to be Boltzamnn type. Subsequently, the stress and displacement distribution in the surrounding rock mass around tunnel is determined by using Laplace transformation method. Combining the sliding crack model and equivalent crack method, the timedependent evolution of cracks is visualized. In the present study, the whole construction process, from excavation to lining, is taken into account. Meanwhile, the failure and deformation behavior at these stages are obtained respectively. It is obvious from the failure process that not only the extent but also the direction of fractured zones is determined. This is key important for the construction of rock bolt. Through sensitivity analysis of parameters, we have: (1) Both excavation velocity in and perpendicular to the tunnel section take effects on the displacement of surrounding rock mass. The final displacement increases with an increase of excavation velocity. (2) Because of the time dependence of rock mass, the failure and deformation behavior continue increasing, even if the excavation process finished. In addition, the displacement after lining is set up is far more than the one at the moment of excavation just finished. The fractured zones characterized as localization in some part of the surrounding rock mass. (3) This model can be used to optimize the construction velocity, thickness of lining and time for setting up the lining system.

Even if there are still many assumptions adopted in present model, the analysis model shows good agreement with engineering practices. In addition, this model cannot consider the fault at the present stage, and the harmonization of crack growth is simplified. Therefore, there still need some improvements. In the future work, in order to couple the micro and macroscopic problem, the two-scale method and numerical simulation will be adopted. Acknowledgements This work is supported by the National Natural Science Foundation of China (Nos. 40902078, 41172243), the Fundamental Research Funds for the Central Universities and the visiting scholar foundation of key laboratory of new technology for construction of cities in mountain area in Chongqing University. References Bobet, A., 2009. Elastic solution for deep tunnels application to excavation damage zone and rockbolt support. Rock Mech. Rock Eng. 42, 147–174. Fahimifar, A., Tehrani, F.M., Hedayat, A., Vakilzadeh, A., 2010. Analytical solution for the excavation of circular tunnels in a visco-elastic Burger’s material under hydrostatic stress field. Tun. Undergr. Sp. Tech. 25, 297–304. Francois, B., Dascalu, C., 2010. A two-scale time-dependent damage model based on non-planar growth of micro-cracks. J. Mech. Phys. Solids. 58, 1928–1946. Hommand-Etienne, F., Hoxha, D., Shao, J.F., 1998. A continuum damage constitutive law for brittle rocks. Comput. Geotech. 22 (2), 135–151. Li, S.C., Wang, M.B., 2008. Elastic analysis of stress–displacement field for a lined circular tunnel at great depth due to ground loads and internal pressure. Tun. Undergr. Sp. Tech. 23, 609–617. Liu, B.G., Du, X.D., 2004. Visco-elastical analysis on interaction between supporting structure and surrounding rocks of circle tunnel. Chin. J. Rock Mech. Eng. 23 (4), 561–564 (in Chinese). Martino, J.B., Chandler, N.A., 2004. Excavation-induced damage studies at the underground research laboratory. Int. J. Rock Mech. Min. Sci. 41, 1413– 1426. Meglis, L.L., Chow, T., Young, R.P., 2001. Assessing microcrack damage around a tunnel at the underground research laboratory using ultrasonic velocity tomography. In: Thirty-Eighth Rock Mechanics Symposium, DC Rocks, July 7– 10. Washington DC, USA, pp. 919–925. Mitaim, S., Detournay, E., 2004. Damage around a cylindrical opening in a brittle rock mass. Int. J. Rock Mech. Min. Sci. 41, 1447–1457. Mortazavi, A., Molladavoodi, H., 2012. A numerical investigation of brittle rock damage model in deep underground openings. Eng. Fract. Mech. 90, 101– 120. Paliwal, B., Ramesh, K.R., 2008. An interacting micro-crack damage model for failure of brittle materials under compression. J. Mech. Phys. Solids 56, 896–923. Schutte, H., Bruhns, 2002. On a geometrically nonlinear damage model based on a multiplicative decomposition of the deformation gradient and the propagation of microcracks. J. Mech. Phys. Solids 50, 827–853. Schwartz, C.W., Einstein, H.H., 1978. Improvement of ground-support performance by full consideration of ground displacements. Transp. Res. Rec. 684, 14–20. Schwartz, C.W., Einstein, H.H., 1980. Simplified analysis for groud-structure interaction in tunneling. In: Proceedings of the 21st US Symposium on Rock Mechanics. University of Missouri-Rolla, pp. 787–796. Shen, B., Barton, N., 1997. The disturbed zone around tunnels in jointed rock masses. Int. J. Rock Mech. Min. Sci. 34 (1), 1117–1125. Wang, H.N., Zhong, Z., 2012. Analytical analysis of supporting circular tunnel construction in viscoelastic rock mass. Chin. J. Appl. Mech. 29 (1), 1–8 (in Chinese). Zhou, X.P., Yang, H.Q., 2007. Micromechanical modeling of dynamic compressive responses of mesoscopic heterogenous brittle rock. Theor. Appl. Frac. Mech. 48, 1–20. Zhu, W.C., Bruhns, O.T., 2008. Simulating excavation damaged zone around a circular opening under hydromechancial conditions. Int. J. Rock Mech. Min. Sci. 45, 815–830.