An elasto-viscoplastic model for soft rock around tunnels considering overconsolidation and structure effects

Computers and Geotechnics 50 (2013) 6–16

Contents lists available at SciVerse ScienceDirect

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

An elasto-viscoplastic model for soft rock around tunnels considering overconsolidation and structure effects Hehua Zhu a, Bin Ye a,⇑, Yongchang Cai a, Feng Zhang b a b

Department of Geotechnical Engineering, College of Civil Engineering, Tongji University, Siping Road 1239, Shanghai 200092, China Department of Civil Engineering, Nagoya Institute of Technology, Showa-ku, Gokiso-cho, Nagoya 4668555, Japan

a r t i c l e

i n f o

Article history: Received 28 October 2012 Received in revised form 13 December 2012 Accepted 13 December 2012

Keywords: Elasto-viscoplasticity Constitutive model Soft rock Structure Tunnel

a b s t r a c t When evaluating the long-term stability of existing tunnels, the creep behavior of soft rock around the tunnel should be properly considered. It is also important to understand the failure mechanism of soft rock when designing the mitigation and remediation of the failure around a tunnel. In this paper, an elasto-viscoplastic model is first modified so that the overconsolidation effect and the structure effect of soft rock can be considered. Then, the performance of the modified model is confirmed with drained triaxial compression tests and creep tests on a manmade rock produced with gypsum and diatom clay. Based on the modified model, finite element analyses are conducted to simulate the model tests of an existing tunnel constructed within manmade rock. Two kinds of model tests are simulated: one is loading failure test and the other creep failure test. The good agreement between the numerical results and the test data validates the performance of the modified constitutive model and the applicability of the corresponding FEM for evaluating the creep failure behavior of an existed tunnel constructed in soft rock. Ó 2012 Elsevier Ltd. All rights reserved.

1. Introduction Soft rock is linked to various geotechnical engineering problems, such as the instability of slopes, foundations and tunnels. In general, the mechanical behavior of soft rock is elasto-plastic, shear-dilatant and strain softening. Sometimes it may also show time-dependent behavior such as creep and strain-rate dependency. In geotechnical engineering, ground improvement with cement mixture is often used to enhance the strength of soft ground, in which the mixture usually behaves like an unconsolidated soft rock with high porosity. When evaluating the long-term stability of existing tunnels, some of which were constructed in the soft rock, it is important to properly consider the above-mentioned mechanical behavior of the ground surrounding the tunnel. To design the mitigation and remediation of cracked tunnels constructed in soft rock with high risk of creep failure, the ground materials also need to be properly considered. Many elasto-plastic or elasto-viscoplastic constitutive models, such as those by Adachi and Oka [1], Ohnaka et al. [2], Nicolae [3], Zhou and Zhu [4], Weng et al. [5], have been proposed to describe the mechanical behavior of soft rock. Although these models can successfully predict the deformation features of some soft rocks at different conditions, they have an obvious limitation in ⇑ Corresponding author. Tel./fax: +86 21 65982384. E-mail address: [email protected] (B. Ye). 0266-352X/$ - see front matter Ó 2012 Elsevier Ltd. All rights reserved. http://dx.doi.org/10.1016/j.compgeo.2012.12.004

that they do not fully considerate the structure effect of soft rock. The structure of a soft rock, as mentioned here, has two meanings. First it is just a geologic term for a soft rock called random fabric structure formed in its sedimentation process, sometimes also called bounding. Second it means that under the same loading, the void ratio of a structured material can be maintained at a higher level than that of a non-structured material because of the extra strength from the structure or bounding of the geomaterials. The existing constitutive models for soft rock mostly take into account the first structure effect, but neglect the second one. The second structure effect can have a great influence on the behavior of soft rock after loading. For example, when a load exceeds the strength of a structured soft rock, strain-softening behavior accompanied by volume contraction, not expansion (dilation), might occur because of the high porosity. Different volume changes of surrounding soils might cause totally different forces acting on tunnel linings [6]. Therefore, it is necessary to consider the second structure effect when evaluating the stability of tunnels excavated in soft rock. Numerical analysis is the main tool for evaluating the stability of tunnels. Zhu et al. [7] carried out a numerical simulation of the deformation and failure of the rock masses surrounding tunnels with finite element method (FEM). Bizjak and Petkovsek [8] calculated the propagation of a tunnel destressed zone and stressfield around a tunnel with finite difference method (FDM). Sterpi and Cividini [9] made a finite element simulation to investigate

7

H. Zhu et al. / Computers and Geotechnics 50 (2013) 6–16

the stability of shallow tunnels in strain softening rock. Lee et al. [10] studied the tunnel stability and the arching effect of the surrounding soft ground using FEM. Funatsu et al. [11] investigated the effects of ground support and reinforcement on the stability of tunnels using distinct element method (DEM). Jia and Tang [12] carried out a FEM analysis on the stability of a tunnel in jointed rock mass. Genis [13] conducted a three-dimensional stability analysis of a tunnel portal using FEM. Maghous et al. [14] numerically investigated rock deformation in bolt-supported tunnels by a simplified homogenization approach. In these analyses, although different numerical techniques were adopted, the structure effect that might cause great volume contraction of surrounding ground was not considered. Some other prior studies have discussed the shear contraction behavior of surrounding ground (e.g. Hellmich and Mang [6], Rutqvist et al. [15]), but these studies mainly dealt with the normally compressed soils whose contraction effect is not so significant as the structured soil. The main purpose of this study is to propose an effective method for evaluating the long-term stability of tunnels constructed in structured soft rock. To do so, first an elastoviscoplastic constitutive model proposed by Zhang et al. [16] is modified in order to properly consider both the effect of consolidation or density and the structure effect. It is known that the concept of subloading proposed by Hashiguchi and Ueno [17] or the concept of boundary surface proposed by Dafalias and Popov [18] can be used to describe the overconsolidation effect. The concept of superloading proposed by Asaoka et al [19] is known to be able to describe the structured soil properly. Therefore, it is natural to introduce this concept to the original model for the soft rock with structure effect. The original model was proposed in the framework of a transformed stress tensor named tij stress tensor [20] that can properly consider the influence of the intermediate principal stress. The performance of the modified model is then confirmed with drained triaxial compression tests and triaxial creep tests on a manmade rock produced from the mixture of gypsum and diatom clay. Second, based on the modified elasto-viscoplastic model, numerical simulations by finite element method (FEM) are conducted to analyze the long-term stability of an existing tunnel in soft rock due to the creep behavior of the ground. The calculations were executed on two types of model tunnel tests, the loading failure test and the creep failure test [21]. The ground in the model tests is manmade rock produced from a mixture of gypsum and diatom clay, whose mechanical behavior is very similar to that of a soft sedimentary rock with high porosity. By comparing the finite element analysis results with the corresponding model test data, the effectiveness of the proposed numerical method is verified.

2. Derivation of the modified constitutive model Zhang et al. [16] proposed an elasto-viscoplastic model (hereafter called as the original model) that can properly take into consideration of the influence of intermediate principal stress by introducing the tij concept [20]. The model, however, cannot consider the structural efffect that may lead to a volumatric shrinkege during post-peak shearing process due to the collapse of the strucutre, which may occur in soft rocks or manmade rock produced from the mixture of gypsum and diatom clay with high porosity. In order to describe the structural effect, the original model is modified in this paper by introducing the superloading concept [19] so that it can describe not only the strain-softning behavior accompanied by positive dilatancy but also the strainsoftning behavior accompanied by negative dilatancy or shear shrinkege.

Fig. 1. Yielding function for normally consolidated soil in tij space after Nakai and Hinokio [22].

Fig. 2. Stress-dilatancy relation in tij space after Nakai and Hinokio [22].

2.1. Description of model As shown in Fig. 1, the yield surface of the proposed model is defined in a tS–tN stress space, in which tN and tS are respectively the modified mean effective stress and the deviatoric stress in the tij stress space. Detailed description of the modified stress tensor tij can be referred to Appendix I. Fig. 2 shows the relationship between the stress ratio and incremental plastic strain ratio. Based on the normality rule, the following relation can be obtained: p dtN  dep SMP þ dt S  dcSMP ¼ 0

ð1Þ

p where, the directions of dep SMP and dcSMP are the plastic volumetric and shear strains defined in the tij stress space and coincide with those of dtN and dtS respectively, because coaxialty between stress and plastic strain increment are assumed. Substituting Eq. (1) into the relation shown in Fig. 2, the following can be derived

Y¼

dep dtS M b  X b SMP ¼ p ¼  dt N dcSMP X b1

ð2Þ

where b is a material parameter that controls the shape of the stress-dilatancy relation shown in Fig. 2. Integrating Eq. (2), the following yield function for normally consolidated soil can be obtained:

f ¼ ln

tN tN tN1 þ 1ðXÞ ¼ ln þ 1ðXÞ  ln ¼0 t N1 tN0 t N0

1ðXÞ ¼

 b 1 X ; b M

X¼

tS tN

ð3Þ

ð4Þ

where, tN and X  tS/tN are respectively the mean stress and the stress ratio based on tij-concept, and tN1 determines the size of the yield surface (the value of tN at tS = 0). tN0, a reference mean stress, is a reference value and is taken as 98 kPa in this paper. From the equation shown in Fig. 2, the M in Eq. (4) can be easily derived using the principal stress ratio XCS  (tS/tN)CS and the incremental p plastic strain ratio Y CS  ðdep SMP =dcSMP ÞCS at critical state:

 1=b M ¼ X bCS þ X b1 CS Y CS

ð5Þ

8

H. Zhu et al. / Computers and Geotechnics 50 (2013) 6–16

The ratios XCS and YCS can be calculated by RCS, the principal stress ratio at critical state in triaxial compression test, as below [23]:

X CS ¼

pffiffiffi   2 pffiffiffiffiffiffiffi 1 RCS  pffiffiffiffiffiffiffi ; 3 RCS

pffiffiffiffiffiffiffi 1  RCS Y CS ¼ pffiffiffi pffiffiffiffiffiffi 2 ð Rcs þ 0:5Þ

RCS ¼ ðr1 =r3 Þcritical state

ð6Þ ð7Þ

In this model, since the associated flow rule is adopted, the yielding function is the same as the plastic potential. Eq. (3) just gives a normal yielding surface that can only be used for normally consolidated soil. For overconsolidated and structured soils, however, further modification is needed. In order to consider the effect of density (or overconsolidation) and structure of a soft rock, the concepts of subloading [17] and of superloading [19] are introduced into the proposed model. Fig. 3 shows the interrelation between the normal yielding surface, subloading surface and superloading surface. All these yielding surfaces are supposed to be similar to each other [19]. As shown in Fig. 3b, under a prescribed mean stress tN1, there are three different void ratios. The lowest one represents the current state (overconsolidated state) and the middle one the normally consolidated state. The uppermost one represents the structural state whose physical meaning is that the void ratio at this state can be maintained at a higher level than that at the normally consolidated state under the same stress level because of the extra strength from the structure or bounding of the geomaterials. From the figure, it is easy to calculate the void ratio difference between the structural and the current states q and the void ratio difference between the structural and the normally consolidated  espectively as states q

q ¼ ðk  jÞ ln

tN ¼ ðk  jÞ ln R; t N1S

¼ ðk  jÞ ln R

q ¼ ðk  jÞ ln

t N1e tN1S



ept ¼

kj t N1e tN1e tN1e ept ln ¼ C p ln ) ln ¼ ; 1 þ e0 tN0 t N0 t N0 Cp

Cp ¼

kj 1 þ e0

where

fs ¼ ln

tN tN tN1 þ 1ðXÞ ¼ ln þ 1ðXÞ  ln ¼0 t N1 tN0 t N0

ð11Þ

tN t N1e tN1 tN1e þ 1ðXÞ  ln  ln þ ln ¼0 t N0 tN0 t N1S t N1S

ð12Þ

or

fs ¼ ln

As shown in Fig. 3a, tN1, tN1e and tN1S represent the size of the subloading surface, the normal yielding surface and the superloading surface (the cross sections of the surfaces with the axis tS = 0) respectively. By substituting Eqs. (9) and (10) into Eq. (12), the subloading surface can be rewritten as

fs ¼ fr 

1 p ðe þ C p ln R  C p ln R Þ ¼ 0 Cp t

represent respectively the overconsolidation state variable and the structure state variable, both of which vary from zero to one. R = 0 means that the soil is at a heavily overconsolidated state while R = 1 means that the soil is at a normally consolidated state. R = 0 means that the soil is at a very high structured state while R = 1 means that the soil is at the state where the structure is destroyed completely.

ð13Þ

where

fr ¼ ln

 b tN tN 1 X þ 1ðXÞ ¼ ln þ t N0 t N0 b M

ð14Þ

The consistency equation of the yielding function can then be expressed as

ð15Þ

Meanwhile, the stress rate can be calculated by the following equation

r_ ij ¼ Eijkl ½e_ kl  e_ pkl  ð9Þ

ð10Þ

Similar to the reasoning of Eq. (3), it is easy to obtain the subloading surface in the following equation,

_  R_  =R Þ ¼ 0 f_ s ¼ f_ r  ðe_ pt =C p þ R=R ð8Þ

tN 1 tN1e R¼ ¼ ð0 6 R 6 1Þ ð0 6 R 6 1Þ and R ¼ tN1S OCR t N1S

From Fig. 3a and based on the definition of normal yielding surface for normally consolidated clay, the plastic volumetric strain epv , accumulated from the initial state (tN0, 0) to the normally consolidated state (tN, tS) or (tN1, 0), can be calculated by the following equation

ð16Þ

In the proposed model, the associated flow rule is defined in the tij stress space and the plastic strain rate tensor is calculated by

e_ pij ¼ K

@fs @tij

ð17Þ

From Eqs. (16) and (17), it is easy to obtain the following relation,

@fr @f @f @f @f f_ r ¼ r_ ij ¼ r Eijkl ðe_ kl  e_ pkl Þ ¼ r Eijkl e_ kl  r Eijkl s K @ rij @ rij @ rij @ rij @t kl

Fig. 3. Yield surfaces using subloading and superloading concepts and their corresponding e-lntN relations.

ð18Þ

9

H. Zhu et al. / Computers and Geotechnics 50 (2013) 6–16

In order to give a complete expression for K, it is necessary to define properly the evolution equations for the overconsolidation state variable R and the structure state variable R. The evolution equation for R can be defined by the following expression

1 K ln R  hðtÞR R_ ¼ mR RC n lnð1þt=t1 Þ C p tN

ð19Þ

In the above equation, parameter mR is used to control the losing rate of overconsolidation while the parameter Cn, the same as the original model, is used to describe the strain rate dependency of the soft rock. As to the function for time effect, the same form as adopted in the original model is used:

hðtÞ ¼ e_ 0t ð1 þ t=t 1 Þa

ð20Þ

where e_ 0v is the volumetric strain rate at time t = 0, in which the time t = 0 represents the moment when shearing begins. t1 is a unit time used to standardize the time. It should keep in mind that a creep stress is usually not applied abruptly to a sample and the loading process for the creep stress should be taken into consideration. The parameter a is a time-dependency parameter mainly used to control the gradient of strain rate vs. time in logarithmic axes. As to the evolution equation for the structure, the idea proposed by Asaoka et al. [19] is adopted and its detailed expression is given below

1 K R_  ¼ mR R2 ð1  R Þ C p tN

ð21Þ

where, the parameter mR is used to control the decaying rate of structure. By substituting Eqs. (18), (19), and (21) into Eq. (15), it is easy to calculate the positive variable K as,



 p @fr hðtÞ @fr @fs h Eijkl e_ kl þ Eijkl þ Cp @ rij @ rij @t kl C p   @fr hðtÞ D Eijkl e_ kl þ ¼ @ rij Cp

K¼

ð22Þ

where p

D¼

@fr @fs h Eijkl þ @ rij @tkl C p

ð23Þ

and hp is calculated by p

h ¼

  @fs ln R Cn lnð1þt=t1 Þ R  mR þ mR R ð1  R Þ tN R @t ii

ð24Þ

The loading criteria are given in the same way as in the original model: 8_ > < f r > 0 hardening ke_ pij k > 0 if K > 0 and f_ r < 0 softening ke_ pij k ¼ 0 if K 6 0 elastic > :_ f r ¼ 0 creep

ð25Þ Different from any other constitutive models, the loading criteria can not only judge the strain hardening and strain softening, but also the creep state which is regarded as a loading process. The rate of viscoplastic strain tensor can then be calculated by @f

e_ pij ¼ K

r Emnkl e_ kl þ hðtÞ=C p @fs @fs @ rmn ¼ @t ij D @tij

ð26Þ

r_ ij ¼ ðEijkl  Epijkl Þe_ kl  AEijqr where Epijkl ¼

@fs @tqr

r Eijqr Emnkl @@f rmn

D

@fs @tqr

ð28Þ

;

A¼

hðtÞ D  Cp

ð29Þ

Eq. (28) is used to calculate the rate of stress tensor from the rate of strain tensor. 2.2. Determaination of the parameters involved in the proposed model Ten parameters are involved in the modified model. Compared with the original model [16], an additional parameter mR is added to control the decaying rate of the structure. The methods for determining the two parameters Cn and a are almost the same as those in the proposed model. Parameter Cn is determined by the difference of peak strengths of two rock samples under drained triaxial compression tests with different constant shear strain rates. Parameter a is the gradient of strain rate vs. time in logarithmic axes and can be uniquely determined with drained triaxial creep test. The parameter b controls the shape of the stress-dilatancy relation shown in Fig. 2. If b takes a value equal to 1, the plastic potential will be the same as that of the Cam-clay model. Anyway, it is very easy to determine the parameter b based on the test data of stress-dilatancy relation. The compression index k and the swelling index j can easily be determined by isotropic consolidation test of soft rock. RCS, the stress ratio at failure state, can also be easily determined with triaxial compression test, and can be calculated by equation 0 0 RCS = (1 + sin /0 )/(1  sin / ), where / is an effective internal frictional angle. The overconsolidation parameter mR, which controls the losing rate of the overconsolidation, can be determined from the peak strengths of soft rock specimens under triaxial compression test by curve fitting, in which at least two samples, with the same physical conditions but different initial overconsolidation, should be tested under the same loading condition. The structure parameter mR , used to control the decaying rate of structure, is determined from the dilatancy relation of soft rock specimens under triaxial compression tests by curve fitting, in which at least two samples, with the same physical conditions but different initial structure state, should be tested under the same loading condition. It should be pointed out that, unlike the overconsolidation that can be rather easily identified, the structure state is rather difficult to be identified. Therefore, the curve-fitting of the structure parameter mR is the most time-consuming task in determining all the ten parameters involved in the proposed model. 2.3. Performance of the proposed model The validity of the modification is verified by the conventional triaxial compression and creep tests on the manmade soft rock under drained condition. The manmade soft rock is produced from gypsum and diatom clay and its composition and material properties are listed in Tables 1 and 2. Conventional triaxial compression tests were conducted under different confining stresses and shear

As to the rate of stress tensor, it can be evaluated by Table 1 Material weight ratio of manmade rock.

r_ ij ¼ Eijkl ðe_ kl  e_ pkl Þ ¼ Eijkl e_ kl  Eijqr e_ pqr ¼ Eijkl e_ kl 

r Eijqr Emnkl @@f rmn

D

@fs @t qr

e_ kl



Eijqr hðtÞ=C p @fs D @tqr

ð27Þ

Gypsum

Diatom clay

Water

Retardation material (%)

1.0

0.3

1.0

0.4

10

H. Zhu et al. / Computers and Geotechnics 50 (2013) 6–16

Table 2 Physical properties of manmade rock. Unit weight (g/cm3)

Water content (%)

Uniaxial strength (MPa)

1.45–1.51

67.9–73.5

1.68–1.83

Table 3 Material parameters of manmade rock. Passion’s ratio v Reference void ratio N(p0 = 98 kPa)

0.15 1.05 10.9 (0.1 MPa) 6.90 (0.3 MPa) 3.69 (l.0 MPa) 0.084 0.0051 1.5 0.40 2.1 0.65 0.050

Stress ratio RCS( = r1/r3)critical Compression index k Swelling index j Parameter of potential shape b Overconsolidation parameter mR Structure parameter mR Secondary consolidation parameter a Time dependency parameter Cn

Table 4 Initial values of state parameters for manmade rock.

Initial R0 Initial R0

r3 = 0.1 MPa

r3 = 0.3 MPa

r3 = 1.0 MPa

0.042 0.25

0.025 0.25

0.035 0.25

strain rates at the constant-strain-rate loading condition. The confining stresses used in the tests are r3 = 0.1, 0.3 and 1.0 MPa and the loading rate are 1.0, 0.1 and 0.01%/min respectively. The creep stresses used in the creep tests are 80% of the peak stresses of the specimen at the same confining stress as in triaxial compression test. Table 3 lists the material parameters of the manmade

:1.0 MPa : 0.3 MPa : 0.1 MPa

Strain rate 0.1%/min Solid line: Theory Dash line: Test

4

10

Volumetric strain (MPa)

Stress difference (MPa)

5

rock and Table 4 shows the initial values of the state parameters for the manmade rock. Fig. 4 shows the comparison of theoretical and test results of the stress-strain-dilatancy relation under different confining stresses in triaxial compression tests. The test results revealed the fact that if the confining stress increased from 0.1 MPa to 1.0 MPa, the stress-strain relation changed from brittle to ductile and the dilatancy changed from negative to positive, showing a typical behavior of overconsolidated geomaterial. It is known from Eq. (9) that the overconsolidation ratio OCR is estimated to be 40.0 under the condition of r3 = 0.3 MPa. Nevertheless, the dilatancy still shows negative value, that is, shear contractive. The reason why such behavior happened can only be explained with the fact that the structure can sustain a relatively higher stress at the same porosity. When strain softening occurred, the porosity collapsed, leading to a contraction instead of dilation that can be observed in overconsolidated materials without the structure effect. The theoretical prediction can satisfactorily simulate this behavior. Fig. 5 shows the comparison of theoretical and test results with different loading rate in triaxial compression tests. In the test, the strain rate dependency is quite evident. When the loading is faster, the peak strength gets higher and the dilatancy changes from negative to positive. Although the theoretical prediction can qualitatively describe this behavior, accuracy is not so good, implying that further modification is needed. Fig. 6 shows the comparison of theoretical and test results in triaxial creep tests, in which the creep stresses under two different confining stresses were set to be 80% of each peak strength obtained in the corresponding triaxial compression tests. The test results were reproduced quite well. It is confirmed that the strength and the dilatant behavior of the specimens under drained and undrained conditions are almost the same and that there is no water observed along the shear band of the specimen at residual state, which means that pore water has

3 2 1 0

0

-10 0

5

10

15

0

5

Axial strain (%)

10

15

Axial strain (%)

Fig. 4. Comparison of theoretical and test results with different confining stress in triaxial compression tests.

10

Confirming stress = 0.1MPa Solid line: Theory Dash line: Test

4 3

:1.0 %/min.

2

: 0.1 %/min. : 0.01 %/min.

Volumetric strain (MPa)

Stress difference (MPa)

5

1 0

0

-10 0

5

10

Axial strain (%)

15

0

5

10

Axial strain (%)

Fig. 5. Comparison of theoretical and test results with different loading rate in triaxial compression tests.

15

11

H. Zhu et al. / Computers and Geotechnics 50 (2013) 6–16

Photo 1. Failure pattern of tunnel surface in loading failure test. Fig. 6. Comparison of theoretical and test results in triaxial creep tests (creep stress = 80% of peak strength).

500 mm

Loading frame Oil jacks

500mm

Model tunnel

106 mm

500 mm

Loading Plate

Flat jack

Loading frame 500mm Fig. 9. Finite element mesh. Fig. 7. Loading device used in model tests.

3. Numerical simulation on model tests of existing Tunnel

Loading vertical stress (MPa)

In the work by Tasaka et al. [21], model tests of an existing tunnel in manmade rock as discussed in Section 2.3 were conducted to investigate the long-term stability of tunnels. The main advantage of using manmade soft rock as the model ground is that the mechanical properties of the material are almost the same for all specimens while the natural soft sedimentary rocks have mechanical properties which are often widely scattering. 1.5

Loading rate 0.025MPa/min Vertical

1.0 Loading failure test

0.5 0.0

Horizontal 0.3MPa 0

10

20

30

40

50

60

The device for the model test is shown in Fig. 7. The model ground is 500 mm high, 500 mm wide and 145 mm thick with a circle tunnel of diameter 106 mm located at the center. The loads were applied with three jacks in the vertical direction and two flat jacks in the horizontal direction on two side surfaces so that it is possible to apply vertical and horizontal stresses uniformly on the model ground. Strain gauges along three different directions were coated on the surface of the tunnel and on the front vertical planes. Displacement transducers and photo devices were also used to measure the displacements of the model tunnel in three directions [21]. Two types of loading tests were conducted. One is the loading failure test, in which the load was applied continuously at a

Loading vertical stress (MPa)

little influence on the mechanical behavior of the manmade soft rock.

1.5 Loading rate 0.025MPa/min

Vertical

1.0

1.1MPa Creep failure test

0.5 0.0

Horizontal 0.3MPa 0

10

20

30

Time (min)

Time (min) Fig. 8. Loading paths of model tests.

40

50

60

12

H. Zhu et al. / Computers and Geotechnics 50 (2013) 6–16

pffiffiffiffi pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi Fig. 10. Calculated distribution of deviatoric stress ( J2 ¼ 0:5Sij Sij ) at different loading stages.

pffiffiffiffi pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi Fig. 11. Calculated distribution of deviatoric strain ( I2 ¼ 0:5eij eij ) at different loading stages.

constant rate of 0.025 MPa/min until the model failed. The failed vertical load is rv = 1.4 MPa under a confining stress of rh = 0.3 MPa. The other test is the creep failure test in which the vertical load was first applied to a prescribed value rv = 1.1 MPa under a confining stress rh = 0.3 MPa and then the load was kept for about 80 min untill the ground around the tunnel surface failed. Fig. 8 shows the loading paths of the types of two tests. The tests inidcate that the failure of the model tunnel always starts along the tunnel surface in both the loading failure test and the creep failure test (see Photo 1). Because the cracks are evenly distributed along the thickness direction, it can be concluded that the geometry and loading condition are in plane-strain condition.

To simulate the model tests, a FEM code ‘SOFT’ [24,25]; was used, in which the mechanical behavior of the model ground was described by the modified model presented in previous section. Considering the symmetric geometry and loading conditions, only half of the model ground was considered. Although the boundary value problem (BVP) of the model tests are three-dimensional (3D), the mechanical behavior of the model tunnel is found to be mainly in plane-strain condition and thus, for simplicity, the analyses were conducted in two-dimensional (2D) condition. Fig. 9 shows the finite element mesh used in the calculation. The boundary conditions in the analysis are the same as those in the Fig. 10 shows the calculated distribution of shear stress pffiffiffiffitests. pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ( J 2 ¼ 0:5Sij Sij , Sij: deviatoric stress tensor) at different loading

H. Zhu et al. / Computers and Geotechnics 50 (2013) 6–16

Fig. 12. Test and calculated displacements around tunnel periphery in loading failure test.

Fig. 13. Test and calculated displacements around tunnel periphery in creep failure test.

13

stages. It can be seen that during the loading failure test, the deviatoric stress increased to a very large value in a wide area. While in the creep failure test, the stress was almost unchanged even at the creep failure stage. On the other hand, if we take a look at the calpffiffiffiffi pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi culated distribution of shear strain ( I2 ¼ 0:5eij eij , eij: deviatoric stain tensor) in the creep failure test shown in Fig. 11, it can be pffiffiffiffi seen that the shear strain I2 keeps increasing, showing the typical creep behavior of a ground. It is clear from Fig. 11 that, after the creep stress was kept for 80 min, a concentrated large strain occurred at the tunnel shoulder area. This trend of strain concentration was also observed by Hellmich and Mang [6]. Their simulation results showed that a large plastic zone extended upwards from the tunnel shoulder area along an angle of 45° from the horizontal axis. Figs. 12 and 13 show the comparison between the calculated and test results of displacements versus time around tunnel periphery at different points respectively for the loading failure test and the creep failure test. The numerical results are in good agreement with the test results, especially for the creep failure test, indicating that the finite element analysis based on a suitable constitutive law can appropriately describe the mechanical behavior of tunnels in soft rock. Fig. 14 shows the comparison between the calculated and measured periphery strains at the locations 5 mm away from the tunnel surface. It can be seen that in the loading test, the calculated and measured strains are in good agreement, while in the creep failure test, the calculated creep strains are much smaller than the measured ones. The main reason for this discrepancy might be due to the fact that discontinued shear bands or cracks occurred at localized areas in the test, while the calculation is still kept in continuous form. Overall, however, the tendency of the development of creep strains observed in the test can be simulated quite well. The same conclusion can be drawn for the periphery creep strains at points 20 mm and 40 mm away from the tunnel surface, as shown in Figs. 15 and 16.

Fig. 14. Test and calculated periphery strains at the place 5 mm away from tunnel surface in creep failure test.

Fig. 15. Test and calculated periphery strains at the place 20 mm away from tunnel surface in creep failure test.

14

H. Zhu et al. / Computers and Geotechnics 50 (2013) 6–16

Fig. 16. Test and calculated periphery strains at the place 40 mm away from tunnel surface in creep failure test.

4. Conclusions In this paper, a viscoplastic model formerly proposed by Zhang et al. [16] was modified by introduing the concept of structure to describe not only the strain softening and time-dependent behavior of soft rock but also the volumatric shrinkege due to the collapse of the strucutre during shearing. Based on the modified model, 2D finite element analyses were conduced to simulate the model tunneltests, whose ground was made of the abovementioned manmade rock. Two types of the model tests, loading failure test and creep failure test, were simulated. The folllowing conclusions can be drawn based on the work in this research: Compared with the original model [16], only one parameter, the structural state parameter mR , is added to control the decaying rate of the structure. Ten parameters are involved in the modified model, among which five parameters, Poisson’s ratio m, reference void ratio N, shear stress ratio at critical state RCS (=r1/r3))critical, compression index k and swelling index j, are the same as the Cam-clay model. These parameters can be determined simply with triaxial compression and creep tests, with the only exception that the structure parameter mR has to be determined using the curvefitting method in which the structure state is rather difficult to be identified. The validity of the modified model is verified by the conventional triaxial compression and creep tests on the manmade soft rock under drained condition. The comparisons between the element test results and the theoretical predictions indicate that the modified model can not only simulate the triaxial compression tests but also the triaxial creep tests uniquely with a rather satisfactory accuracy. The model tunnel tests, whose ground is made of the same manmade rock, aimed to investigate the long-term stability of an existing tunnel, were simulated uisng two-dimensional finite element analysis based on the modified elasto-viscoplastic model. Because the physical and mechanical properties of the manmade soft

rock have very smaller scattering, the model test results are suitable to verify the performance of the simulations. The results indicate that the finite element analysis based on the modified constitutive model can appropriately simulate the mechanical behavior of the model tunnel under the failure loading test. As for the creep failure test, the long-term stability of the model tunnel can also be simulated qualitatively well. Acknowledgements This work was supported by the National Natural Science Foundation of China (Nos. 41130751 and 41002094), National Basic Research Program of China (973 Program: 2011CB013800), the Program for Changjiang Scholars and Innovation Research Team in Universities (PCSIRT, IRT1029). The authors wish to thank Professor Lianyang Zhang at University of Arizona for providing valuable advice. Appendix I. Definition of tij [20] As shown in Fig. A-1, the Spatial Mobilized Plane (SMP) defined by Nakai and Matsuoka can be expressed in a principal stress space (rI, rII, rIII) as,

rI r1

rII r2

rIII r3

pffiffiffiffiffiffi þ pffiffiffiffiffiffi þ pffiffiffiffiffiffi ¼ 1

ðA:I:1Þ

A normal to the plane (a1, a2, a3) can be evaluated by the following equation,

sffiffiffiffiffiffiffiffiffi I3 ai ¼ ði ¼ 1; 2; 3Þ I 2 ri

ðA:I:2Þ

where I1, I2 and I3 are respectively the first, second and third invariants of the effective stress tensor and can be expressed by the following using three principal stresses:

Fig. A-1. Explanation of SMP [18].

15

H. Zhu et al. / Computers and Geotechnics 50 (2013) 6–16

Fig. A-2. Explanation of stress and strain in principal-value space of stress tensor tij [18].

I1 ¼ r1 þ r2 þ r3

Appendix II. Deductions of differential

I 2 ¼ r1 r2 þ r2 r3 þ r3 r1

ðA:I:3Þ @f @ rij

I 3 ¼ r1 r2 r3 In tij clay and sand models, a symmetric tensor tij is expressed by a product of the stress tensor rik and a tensor akj as:

t ij ¼ rik akj

ðA:I:4Þ

where tensor akj can be evaluated by Cayley-Hamilton’s theory as

sffiffiffiffi I3 1 r ; aij ¼ I2 ij

r ik r kj ¼ rij ;

r ij

¼ ðrik þ I2 dik ÞðI1 rkj þ I3 dkj Þ1

ðA:I:5Þ

Because rij is a one-half power function of rij, rij is a symmetric tensor and has the same principal directions as rij and its principal valpffiffiffiffiffiffi pffiffiffiffiffiffi pffiffiffiffiffiffi ues are ( r1 ; r1 ; r1 ). Therefore, aij is also a symmetric tensor and has the same principal directions as rij and has the principal values of (a1, a2, a3). As a result, tij is also a symmetric tensor and has the same principal directions as aij and rij. The normal component tn and the tangential component ts of the principal-value vector of tij can be determined by (see Fig. A-2a),

t N ¼ t 1 a1 þ t 2 a2 þ t3 a3 ¼ tij aij  rSMP ¼ r1 a21 þ r2 a22 þ r3 a23 ¼ 3

I3 I2

ðA:I:6Þ

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ðt21 þ t22 þ t 23 Þ  t 2N ¼ t ij t ij  ðtij aij Þ2  sSMP qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ¼ ðr1  r2 Þ2 a21 a22 þ ðr2  r3 Þ2 a22 a23 þ ðr3  r1 Þ2 a23 a21 qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi I1 I2 I3  9I23 ¼ I2

›f ›rij

Can be written as

^ mn @f @f @ r ¼ ^ @ rij @ rmn @ rij

where @r^@fmn can be obtained by just substituting rij, the stress tensor in normal stress space, and tij, the stress tensor in tij stress space, ^ ij tensor and ^t ij tensor in the expression of @ r@f , as dewith relevant r mn scribed in detail by Nakai and Hinokio [25]. ^ mn @r ^ mn In order to calculate @ rij , the transformation stress tensor r expressed with Eq. (A-II-1) can be rewritten in a suffix form as,

r^ mn ¼ ða þ c  2bÞmkl rlk mmn þ crmn þ ðb  cÞðmmk rkn þ rmk mkn Þ

^ mn @r @ðmkl rlk mmn Þ @ rmn ¼ ða þ c  2bÞ þc @ rij @ rij @ rij  @ðmmk rkn Þ @ðrmk mkn Þ ða þ c  2bÞ þ ðb  cÞ þ @ rij @ rij   @mkl @m  rlk mmn þ mkl dli dkj mmn þ m0kl rlk mn þ cdmi dnj @ rij @ rij  @mmk @m þ ðb  cÞ rkn þ mmk dki dnj þ dmi dkj mkn þ rmk kn @ rij @ rij ðA:II:3Þ Being the orientation tensor of the sedimentation plane, mij is a constant with respect to the stress tensor rij, and it is valid that,

ðA:I:7Þ

It is assumed that the directions of plastic principal strain increments coincide with those of principal axes of tij, and the plastic strain increment can also be given by a normal component (deP SMP ) and a tangential component (dcp SMP ) of the principal plastic strain increment vector on the SMP (see Fig. A-2b),

de

p 1 a1

¼ de

p 2 a2

þ de

p 3 a3

þ de

p ij aij

¼ de

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi p2 p2 p2 ðdep2 1 þ de2 þ de3 Þ  deSMP qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ¼ depij depij  ðdepij aij Þ2

ðA:II:2Þ

Therefore,

tS ¼

P SMP

ðA:II:1Þ

ðA:I:8Þ

dcp SMP ¼

@mkl @mmn @mmk @mkn ¼ ¼ ¼ ¼0 @ rij @ rij @ rij @ rij

ðA:II:4Þ

Consequently, @@r^rmn can be expressed as, ij

^ mn @r ¼ ða þ c  2bÞmkl dli dkj mmn þ cdmi dnj þ ðb  cÞ @ rij  ðmmk dki dnj þ dmi dkj mkn Þ ¼ ða þ c  2bÞmji mmn þ cdmi dnj þ ðb  cÞðmmi dnj þ dmi mjn Þ

ðA:II:5Þ

References

ðA:I:9Þ

[1] Adachi T, Oka F. An elasto-plastic constitutive model for soft rock with strain softening. Int J Numer Anal Meth Geomech 1995;19(4):233–47.

16

H. Zhu et al. / Computers and Geotechnics 50 (2013) 6–16

[2] Ohnaka M, Akatsu M, Mochizuki H, Odedra A, Tagashira F, Yamamoto Y. A constitutive law for the shear failure of rock under lithospheric conditions. Tectonophysic 1997;277(1–3):1–27. [3] Nicolae M. Non-associated elasto-viscoplastic models for rock salt. Int J Eng Sci 1999;37(3):269–97. [4] Zhou CY, Zhu FX. An associated elastic–viscoplastic constitutive model for sandstone involving shear-induced volumetric deformation. Int J Rock Mech Min 2010;47(8):385–95. [5] Weng MC, Tsai LS, Hsieh YM, Jeng FS. An associated elastic–viscoplastic constitutive model for sand stone involving shear-induced volumetric deformation. Int J Rock Mech Min 2010;47(8):1263–73. [6] Hellmich C, Mang A. Influence of the dilation of soil and shotcrete on the load bearing behavior of NATM-tunnels. Felsbau 1999;17(1):35–43. [7] Zhu W, Li S, Li S, Chen W, Lee CF. Systematic numerical simulation of rock tunnel stability considering different rock conditions and construction effects. Tunn Undergr Sp Tech 2003;18(5):531–6. [8] Bizjak KF, Petkovsek B. Displacement analysis of tunnel support in soft rock around a shallow highway tunnel at Golovec. Eng Geol 2004;75(1):89–106. [9] Sterpi D, Cividini A. A physical and numerical investigation on the stability of shallow tunnels in strain softening media. Rock Mech Rock Eng 2004;37(4):277–98. [10] Lee CJ, Wu BR, Chen HT, Chiang KH. Tunnel stability and arching effects during tunneling in soft clayey soil. Tunn Undergr Sp Tech 2006;21(2):119–32. [11] Funatsu T, Hoshino T, Sawae H, Shimizu N. Numerical analysis to better understand the mechanism of the effects of ground supports and reinforcements on the stability of tunnels using the distinct element method. Tunn Undergr Sp Tech 2008;23(5):561–73. [12] Jia P, Tang CA. Numerical study on failure mechanism of tunnel in jointed rock mass. Tunn Undergr Sp Tech 2008;23(5):500–7. [13] Genis M. Assessment of the dynamic stability of the portals of the Dorukhan tunnel using numerical analysis. Int J Rock Mech Min 2010;47(8):1231–41. [14] Maghous S, Bernaud D, Couto E. Three-dimensional numerical simulation of rock deformation in bolt-supported tunnels: a homogenization approach. Tunn Undergr Sp Tech 2012;31:68–79. [15] Rutqvist J, Borgesson L, Chijimatsu M, Hernelind J, Jing L, Kobayashi A, et al. Modeling of damage, permeability changes and pressure responses during

[16]

[17]

[18] [19]

[20] [21]

[22]

[23] [24]

[25]

excavation of the TSX tunnel in granitic rock at URL, Canada. Environ Geol 2009;57(6):1263–74. Zhang F, Yashima A, Nakai T, Ye GL, Aung H. An elasto-viscoplastic model for soft sedimentary rock based on tij concept and subloading yield surface. Soils Found 2005;45(1):65–73. Hashiguchi K, Ueno M. Elastoplastic constitutive laws of granular material. In: Constitutive Equations of Soils. Murayama S, Schofield AN, editors. Proceedings of the 9th international conference of, soil mechanics and foundation engineering; 1977, p. 73–82. Dafalias YF, Popov EP. Boundary surface plasticity. II: Application to isotropic cohesive soils. J Eng Mech Div, ASCE 1976;112:1263–91. Asaoka A, Nakano M, Noda T. Super loading yield surface concept for the saturated structured soils. In: Proceedings of the fourth European conference on numerical methods in geotechnical engineering-NUMGE98; 1998, p. 232– 42. Nakai T, Mihara Y. A new mechanical quantity for soils and its application to elastoplastic constitutive models. Soils Found 1984;24(2):82–94. Tasaka Y, Azuma H, Ohmori T, Sekine Y, Kameya H, Zhang F. Model test on long-term deformation behavior of tunnel and its numerical simulation, Part1: Model test of tunnel considering creep behavior of soft rocks. In: Proceedings of annual national congress of JGS; 2007. vol. 42(775). p. 1543–44 (in Japanese). Nakai T, Hinokio M. A simple elastoplastic model for normally and over consolidated soils with unified material parameters. Soil Found 2004;44(2):53–70. Nakai T, Matsuoka H. A generalized elastoplastic constitutive model for clay in a three-dimensional stresses. Soils Found 1986;26(3):81–9. Ye GL, Zhang F, Naito K, Aung H, Yashima A. Test on soft sedimentary rock under different loading paths and its interpretation. Soils Found 2007;47(5):897–909. Sekine Y, Zhang F, Tasaka Y, Kurose H, Ohmori T. Model tests and numerical analysis on the evaluation of long-term stability of existing tunnel. In: Proceedings of the 17th international conference of, soil mechanics and foundation engineering; 2009, p. 1848–54.