Equivalent analysis of orthogonal viscoelastic jointed rock via an adaptive algorithm in time domain

Finite Elements in Analysis and Design 46 (2010) 875–888

Contents lists available at ScienceDirect

Finite Elements in Analysis and Design journal homepage: www.elsevier.com/locate/finel

Equivalent analysis of orthogonal viscoelastic jointed rock via an adaptive algorithm in time domain Yi Ren, Haitian Yang n State Key Laboratory of Structural Analysis for Industrial Equipment, Department of Engineering Mechanics, Dalian University of Technology, Dalian 116024, China

a r t i c l e in fo

abstract

Article history: Received 29 July 2009 Received in revised form 1 June 2010 Accepted 2 June 2010 Available online 17 June 2010

This paper presents a recursive equivalent constitutive model that is nested in an FE (finite element) based recursive adaptive process to estimate equivalent behaviors of orthogonal viscoelastic jointed rock. An adaptive algorithm in the time domain and a simple equivalence assumption are combined to convert an inhomogeneous time dependent problem into a series of recursive homogeneous time independent problems, which can be solved either by finite element method (FEM) or other well developed numerical schemes. Ten equivalent analysis based numerical solutions are compared with the FE-based heterogeneous solutions via ANSYS. In the viewpoint of balance between computing accuracy and cost, the results achieved by the proposed approach are fairly good. Crown Copyright & 2010 Published by Elsevier B.V. All rights reserved.

Keywords: Jointed rock Viscoelasticity Equivalent analysis Time-domain Adaptive algorithm

1. Introduction Equivalent analysis is an economic way in the term of computational expense for the numerical simulation of highly jointed rock mass [1,2]. This kind of approach regards jointed rock as an equivalent macroscopic continuum with homogeneous constitutive models, based on which the numerical simulation of jointed rock will become more convenient and simple [3]. Since the numerical simulation of a viscoelastic problem is computationally much more expensive than an elastic problem in general, the equivalent analysis of viscoelastic jointed rock seems more meaningful than that of elastic cases, for which there have been a variety of related works [1,3–18]. Some simple equivalent models for the viscoelastic jointed rock were summarized by Chen [19,20] who recently presented an equivalent rheology model, under some simple assumptions for the jointed rock masses, reinforced by fully grouted bolt and shotcrete lining, and exhibited a comparison of Load–Displacement curves obtained by the equivalent analysis and a FE-based heterogeneous solution. But he did not give sufficient computing details, such as the comparison of computing efficiency. The equivalent analysis of viscoelastic jointed rock can also be investigated by micromechanical modeling approaches which have been successfully used to describe the homogeneous constitutive behavior/effective response of elastic, linear/nonlinear viscoelastic, and viscopalstic heterogeneous medium

n

Corresponding author. Tel.: +86 4114708394; fax: + 86 4114708393. E-mail address: [email protected] (H. Yang).

[21–37], such as the double scale asymptotic expansion technique [34], translated fields and self-consistent techniques [25], and the affine homogenization approach [28], etc. A prediction of equivalent viscoelastic properties of unidirectional jointed rock is given by Liu who transformed the local constitutive laws into fictitious linear thermo-elastic relations so that an asymptotic scheme can be applied, and the temporal response is computed after a numerical inversion of Laplace transform [38,39]. The crucial step of such a process lies in the ability to compute the Laplace inverse transform that generally depends on the numerical procedures [40]. The major advantage of micromechanical modeling approaches is rigorous mathematically and physically, but such approaches seem a bit complex and uneasy for those who are not profound of micromechanics. No matter based on a simple physical assumption or a micromechanical frame, and no matter discretized by an FE or other numerical schemes in the spatial domain, a homogenized viscoelastic constitutive model (either in the incremental or in the full value forms) usually needs to be embedded in a step by step based numerical process in the time domain to evaluate the equivalent response of the field. Thereby, a discretized temporal algorithm with sufficient accuracy and fitted for the homogenized viscoelastic constitutive model is required. However, there seems no direct report on this issue. The objective of this paper is to present a simple recursive equivalent constitutive model for the orthogonal viscoelastic jointed rock, and an adaptive temporal algorithm to estimate the equivalent behaviors of the homogenized viscoelastic field. Thereby a simple equivalence assumption and a piecewised expansion based algorithm in the time domain are combined each other, and the local

0168-874X/$ - see front matter Crown Copyright & 2010 Published by Elsevier B.V. All rights reserved. doi:10.1016/j.finel.2010.06.001

876

Y. Ren, H. Yang / Finite Elements in Analysis and Design 46 (2010) 875–888

truncation order at a discretized time interval can be controlled by the increase of powers of the expansion adaptively. Regarding to the different combinations of constitutive relationships for rock and joint, different relative width and orientations of joints, different geometry and different step sizes, ten equivalent analysis based numerical solutions are provided, and compared with an FE-based heterogeneous solutions via ANSYS. In the viewpoint of balance between computing accuracy and expense, the results achieved by the proposed approach are fairly good.

2. Recursive constitutive relationship It is assumed that both joint and rock are viscoelastic, the constitutive relationship of 2D problem is specified by [41] 8 9 Zt  L   L < L  L   L  dJL ðtt1 Þ = tZ0 dt1 e ðtÞ ¼ K s ðt1 Þ J ð0Þ s ðtÞ þ : ; dðtt1 Þ 0þ

ð1Þ where T

fsL ðtÞg ¼ fsLx , sLy , tL g feL ðtÞg ¼ feLx , eLy , gL g L

ð2Þ

T

ð3Þ

For the plane stress problem L L ¼ K22 ¼ 1; K11

L K33 ¼ 2ð1 þ nL Þ;

L K12 ¼ nL

L L ¼ K22 ¼ 1ðnL Þ2 ; K11

L K33 ¼ 2ð1 þ nL Þ;

L K12 ¼ nL ð1 þ nL Þ

EL1 ,

EL2

X

feLm gsm

ð14Þ

where m refers to the power order. n oT fsLm g ¼ sLxm sLym tLm

ð15Þ

n feLm g ¼ eLxm

ð16Þ

eLym gLm

oT

fsLm g and feLm g are the expansion coefficients of {sL(t)} and {eL(t)}. The conversion relationship of differentials from t to s is 2

d 1 d ¼ ; dt T ds

2

d 1 d ¼ 2 dt 2 T ds2

ð17Þ

Substituting Eqs. (13) and (14) into Eq. (9) and equating the power of the two sides of the equation give L fsLm g ¼ ½DL feLm g þ fCm g,

where 0 B ½DL  ¼ @

DL11

DL12 DL22

symm

m ¼ 1, 2, 3, :::

0

ð18Þ

1

0 C A DL33

ð19Þ

For the plane stress problem ðnL Þ2 ; s2 L

DL11 ¼ s1 L þ  L Cm ¼

DL22 ¼

(

nL wL2 sL2

þ wL1



1 ; s2 L

wL2 sL2



DL33 ¼ wL2 sL3

1 ; s3 L

DL12 ¼

nL s2 L

ð20Þ

)T

qL1 pL1

ð21Þ

ð22Þ

sL2 ¼ ð1ðnL Þ2 Þ

pL1 qL1

ð23Þ

ð7Þ sL3 ¼ 2ð1þ nL Þ ð8Þ

EL1

ð13Þ

ð6Þ

where

ZL1

fsLm gsm

m¼0

sL1 ¼

For three-parameter solid model [41]  1 1 L J L ðtÞ ¼ L þ L 1et=t1 E2 E1

L

feL ðtÞg ¼

ð5Þ

For the plane strain problem

X

m¼0

L

{s (t)} and {e (t)} represent stress and strain vectors, respectively, t denotes time, superscript L¼r or L¼ j, ‘r’ and ‘j’ refer to ‘rock’ and ‘joint’, respectively. 0 L 1 L K11 K12 0 B L 0 C K22 ½K L  ¼ @ ð4Þ A L symm K33

tL1 ¼

fsL ðtÞg ¼

wL1 ¼

L 1

nL T 1 m pL1

pL1 qL1

sLyðm1Þ þ

ð24Þ T qL0 L T e  sL m pL1 xðm1Þ pL1 m xðm1Þ

ð25Þ

and Z are constitutive parameters of joint or rock. Eq. (1) can be written in a differential form    dfeL ðtÞg  L   L  dfsL ðtÞg ¼ K ðt 40Þ ð9Þ qL0 eL ðtÞ þ qL1 s ðtÞ þpL1 dt dt

wL2 ¼

TqL ð1ðnL Þ2 ÞT 1 L nL T qL0 L syðm1Þ  exðm1Þ  L 0 eLyðm1Þ L L m m q1 q1 q1 m

ð26Þ



wL3 ¼

2ð1 þ nL ÞT 1 L T qL0 L t  g m qL1 m1 m qL1 m1

ð27Þ

n,

  1   e ðtÞ ¼ L K L sL ðtÞ E2 L

ðt ¼ 0Þ

ð10Þ

For the plane strain problem, EL1 , EL2 and nL in Eqs. (20)–(27) L EL1 EL2 need to be replaced by 1-ðnL Þ2 , 1-ðnL Þ2 and 1n nL , respectively.

where pL1 ¼

ZL1 EL1 þ EL2

;

qL0 ¼

EL1 EL2 ; EL1 þ EL2

qL1 ¼

EL2 ZL1 EL1 þ EL2

ð11Þ

We divide time domain into a number of intervals, where the initial points and sizes of time intervals are defined by t0, t1, t2, y, tk, y and T1, T2, y, Tk, y, respectively. At a discretized time interval, {sL(t)} and {eL(t)} are expanded in the term of s defined by s¼

ttk1 , Tk

s A ½0, 1,

k ¼ 1, 2, 3, . . .

ð12Þ

At the first time interval fsL0 g ¼ ½DL feL0 g At the (k+ 1)th time interval 8 P < feL0 gðk þ 1Þ ¼ i ¼ 0 feLi gk P , k ¼ 1, 2, 3, . . . : fsL0 gðk þ 1Þ ¼ i ¼ 0 fsLi gk

ð28Þ

ð29Þ

where subscript k+1 and k refer to (k+ 1)th and kth time intervals, respectively.

Y. Ren, H. Yang / Finite Elements in Analysis and Design 46 (2010) 875–888

3. Equivalent assumption and recursive equivalent constitutive equation 3.1. Equivalent assumptions The rock intersected by the orthogonal joints is presented in Fig. 1, where base cell is small enough. b and a denote relative width and dip angle of joint, respectively. In local coordinate system x y, the stress and strain in the I, II, III domains are represented by n fsI ðtÞg ¼ sIx ðtÞ

sIy ðtÞ tI ðtÞ

n feI ðtÞg ¼ eIx ðtÞ

oT

n feII ðtÞg ¼ eIIx ðtÞ n

ð30Þ

eIy ðtÞ gI ðtÞ

n fsII ðtÞg ¼ sIIx ðtÞ

III

oT

sIIy ðtÞ tII ðtÞ eIIy ðtÞ gII ðtÞ

III x ðtÞ

III y ðtÞ

fs ðtÞg ¼ s

s

ð31Þ oT

ð32Þ

oT

III

t ðtÞ

ð33Þ oT

ð34Þ

877

It is assumed that 8 sx ðtÞ ¼ bsIIx ðtÞ þ ð1bÞsIx ðtÞ > > > > I > s ðtÞ ¼ bsIII > y ðtÞ þ ð1bÞsy ðtÞ > < y I III sx ðtÞ ¼ sx ðtÞ > > > sIy ðtÞ ¼ sIIy ðtÞ > > > > : tðtÞ ¼ tI ðtÞ ¼ tII ðtÞ ¼ tIII ðtÞ

ð38Þ

8 > ex ðtÞ ¼ eIIx ðtÞ > > > > > e ðtÞ ¼ eIII > y ðtÞ > y > > < ex ðtÞ ¼ beIII ðtÞ þð1bÞeI ðtÞ x x II I > > ey ðtÞ ¼ bey ðtÞ þð1bÞey ðtÞ > > > > > gðtÞ ¼ a1 bgII ðtÞ þ a2 ð1bÞgI ðtÞ > > > : gII ðtÞ ¼ gIII ðtÞ

ð39Þ

where a1, a2 are parameters which need to be determined by experiments or other means. 3.2. Equivalent recursive constitutive equation We expand {s(t)}, {e(t)}, {si(t)} and {ei(t)} in the form X fsðtÞg ¼ fsm gsm

ð40Þ

m¼0

n feIII ðtÞg ¼ eIII x ðtÞ

oT

III eIII y ðtÞ g ðtÞ

ð35Þ

X

feðtÞg ¼

fem gsm

ð41Þ

m¼0

The equivalent stress and strain are represented by  T fsðtÞg ¼ sx ðtÞ sy ðtÞ tðtÞ

ð36Þ

 feðtÞg ¼ ex ðtÞ

ð37Þ

ey ðtÞ gðtÞ

T

fsi ðtÞg ¼

X

fsi m gsm

ð42Þ

m¼0

fei ðtÞg ¼

X

fei m gsm

ð43Þ

m¼0

where i¼I, II, III, fsm g ¼ fsxm

dx )dx

n fsi m g ¼ sixm

y 1−β)d

IV

n fei m g ¼ eixm

(

III

dy

II

ð44Þ

eym gm gT

fem g ¼ fexm

βdy

(1−β

βdx

sym tm gT

ð45Þ

siym tim eiym gim

oT

ð46Þ

oT

ð47Þ

i m}

i m}

and {e are the expansion coefficients of {sm}, {em}, {s {s(t)}, {e(t)}, {si(t)} and {ei(t)}, respectively. Substituting Eqs. (18), (40) and (41) into Eqs. (38) and (39) and equating the power of the two sides of the equations yield

I

α

fsm g ¼ ½DH fem g þ fCm g,

m ¼ 1, 2, 3, . . .

ð48Þ

Rock

Eq. (48) stands for a recursive constitutive equation in the local coordinate system x  ywhere 0 H 1 D11 DH 0 12 B H H 0 C D22 ½D  ¼ @ ð49Þ A symm DH 33

Joint

fcm g ¼ fcm1

Y y

Base cell

cm2

cm3 gT

ð50Þ

It is assumed that nr ¼ nj ¼ n For the plane stress problem

x

DH 11 ¼

α

2

X Fig. 1. Schematic illustration of the orthogonal viscoelastic jointed rock.

2

2

2

3

2 2 2 r j2 2 r2 j ðn2 þ 1Þc bEj3 2 þ ½ðn þ 2Þcb þ ðn c þ n c Þb þ c E2 E2 þ ð1n Þðb þ cbÞE2 E2 2

2

j r r2 2 2 2 ð1n2 Þc Ej2 2 þ 2ð1n ÞcbE2 E2 þ ð1n Þ b E2

ð51Þ

878

Y. Ren, H. Yang / Finite Elements in Analysis and Design 46 (2010) 875–888

h DH 12

i





2 r2 j 3 3 b2 n3 c þ ðnc þ c2 nÞb þ c2 n Er2 Ej2 2 þ ðn þ nÞcn þ n b E2 E2

¼

ð1n

2 Þc2 Ej2 þ 2ð1 2 ÞcbEj Er þ ð1 2 Þ2 b2 Er2 2 2 2 2

n

n

ð52Þ

 Cm1   Cm1 Cm1  b ¼ b1 , b2 , . . ., bCm1 12

ð61Þ

 Cm2   Cm2 Cm2  b ¼ b1 , b2 , . . ., bCm2 12

ð62Þ 2

2

j r r j 2 r r 2 r 2 j bCm1 ¼ T cEj2 1 2 E2 E1 Z1 ðcbn E2 þ cbE2 2b n E2 cbn E2 þ c E2 Þ

DH 22 ¼ 2

2

2

2

ð63Þ

3

2 2 2 r j2 2 r2 j ðn2 þ1Þc bEj3 2 þ ½ðn þ 2Þcb þ ðn c þ n c Þb þ c E2 E2 þ ð1n Þðb þ cbÞE2 E2 2

2

j r r2 2 2 2 ð1n2 Þc Ej2 2 þ2ð1n ÞcbE2 E2 þ ð1n Þ b E2

2

ð53Þ

2

DH 33 ¼

h

2

¼ T cEj2 Zj1 ðn1Þðn þ 1ÞðEr2 þ Er1 Þðc Ej2 þ cbEr2 b n2 Er2 Þ bCm1 2 2

2

bCm1 ¼ T cnEj2 Er2 Er1 Zj1 ðb n2 Er2 2cbEr2 þ b Er2 þ cbEj2 c Ej2 Þ 3

Er2 Ej2

2ð1þ nÞ a1 bEr2 þ a2 cEj2

i

2

2

¼ T cnEj2 Zj1 ðn1Þðn þ 1ÞðEr2 þEr1 Þð2b n2 Er2 2c Ej2 bCm1 4 2

þ cbEj2 þ b Ej2 3cbEr2 Þ 2

2

2

j2 r2 r j 2 ¼ T bEj1 Ej2 Zr1 ðc Er2 bCm1 5 2 þ b E2 þ2bcE2 E2 þ n c E2 2 2 j r 2 2 j r E2 E2 b 2 Er2 2  c E2 E2 Þ

þ bcn

n

n

ð67Þ

h i T ba1 Er2 Zr1 ðEj1 þ Ej2 Þ þ ca2 Ej2 Zj1 ðEr2 þEr1 Þ mZr1 Zj1 ða1 bEr2 þ a2 cEj2 Þ

Y 2275 KN/m

tm1

a1 bTEr2 Ej1 Ej2 gIIm1 a2 cTEj2 Er1 Er2 gIm1 þ 2mZj1 ð1 þ nÞða1 bEr2 þ a2 cEj2 Þ 2mZr1 ð1 þ nÞða1 bEr2 þ a2 cEj2 Þ

974.8KN/m

ð56Þ

50 m

þ

ð66Þ

ð55Þ

fbCm2 g dfCompm1 g J

Cm3 ¼ 

Y 1

5m 2

h i 2 j r 2 2 2 2 r2 J ¼ mZr1 Zj1 ð1n2 Þc Ej2 2 þ2ð1n ÞcbE2 E2 þ ð1n Þ b E2

ð58Þ

c ¼ 1b

ð59Þ

 n Compm1 ¼ eIxðm1Þ

50 m

sIxðm1Þ eIyðm1Þ sIyðm1Þ eIIxðm1Þ sIIxðm1Þ

eIIyðm1Þ sIIyðm1Þ eIII sIII eIII sIII xðm1Þ xðm1Þ yðm1Þ yðm1Þ

oT

ð60Þ

X

45º 45º

3

50 m

where

1167.5KN/m

ð57Þ



ð65Þ

ð54Þ

fbCm1 g dfCompm1 g Cm1 ¼ J

Cm2 ¼

ð64Þ

X

4 5

Fig. 2. Schematic illustration of the structures concerned with Cases 1–5: (a) the structures concerned with Cases 1–5 and (b) the locations of sample points.

Table 1 Cases list. Case no.

Constitutive models Rock

1 2 3 4 5 6 7 8 9 10

Elastic Viscoelastic Viscoelastic Viscoelastic Viscoelastic Viscoelastic Viscoelastic Viscoelastic Viscoelastic Viscoelastic

Joint

Viscoelastic elastic Viscoelastic Viscoelastic Viscoelastic Viscoelastic Viscoelastic Viscoelastic Viscoelastic Viscoelastic

Geometric characteristics a

b

0.2 0.2 0.2 0.1 0.05 0.2 0.2 0.2 0.2 0.2

a

b is relative width of joint (Fig. 1).

b

a is dip angle of joint (Fig. 1).

a

Total DOFs

Shape of cross section of tunnels

Types of tunnel

(1)b

Equivalent model

Heterogeneous model

0 0 0 0 0 30 60 30 30 30

Circular Circular Circular Circular Circular Circular Circular Non-circular Circular Circular/non-circular

Single tunnel Single tunnel Single tunnel Single tunnel Single tunnel Single tunnel Single tunnel Single tunnel Twin tunnel Twin tunnel with cross passage and service tunnel

675 675 675 675 675 1118 1118 850 1396 3040

173,586 173,586 173,586 308,062 312,290 316,360 319,821 320,028 333,946 359,563

Table 2 Constitutive parameters. Material

Rock Joint

Elastic properties

Viscoelastic properties

Elastic modules, EL2 ðpaÞ

Poisson’s ratios, nL

EL1 ðpaÞ

EL2 ðpaÞ

ZL1 ðpaÞ

nL

3.5  1010 2  108

0.3 0.3

1  108 1  109

3.5  1010 2  108

9  109 1  1010

0.3 0.3

Y. Ren, H. Yang / Finite Elements in Analysis and Design 46 (2010) 875–888

Joint (magenta):

Joint (magenta): visco88

visco88/plane82

1.25 m

0.0625 m

1.25 m

0.125 m

1.25 m

0.25 m 8-nodes isoparametric element

Rock (cyan): visco88

Rock (cyan): visco88

Rock (cyan): *visco88/plane82

879

Joint (magenta): visco88

Fig. 3. FE meshes of different models for Cases 1–5: (a) the equivalent model (Cases 1–5); (b) the heterogeneous model (Cases 1–3); (c) the heterogeneous model (Case 4) and (d) the heterogeneous model (Case 5). *Note: visco88 and plane82 are ANSYS 2-D 8-Node viscoelastic and structural solid element.

18

12 10 8 Equ

6

CPU Time

Het

4

CPU Time

=

355 . 698 (s)

= 38522 . 294 (s)

2 200

400 600 Time (day)

800

Step=5(day)

Case 1

(10 m)

Step=1(day) Step=0.1(day)

-2

1

Equ,y

0 -2 -3 -4 0

200

400 600 Time (day)

800

1000

8 6

CPU Time

4

CPU Time

bCm1 ¼ T bZr1 ðn1Þðn þ 1ÞðEj2 þ Ej1 ÞðcEj2 þ nbEr2 þ bEr2 ÞðcEj2 nbEr2 þ bEr2 Þ 6

ð68Þ 2 j r 2 j2 2 r2 r j bEj1 Ej2 r1 ðb2 2 Er2 2 3bcE2 E2 þ c E2 E2 2c E2 b E2 Þ

Z

200

12 10 8 6 4 2 0 -2 -4 -6 -8 -10 -12 -14

Equ Het

=

356 . 670 (s)

= 42903 . 516 (s)

400 600 Time (day)

800

Step=5(day) Step=1(day) Step=0.1(day)

0

Fig. 4. Error(t) curve and the comparison of uðtÞHet,y of Case 1: (a) the Error(t) curve 1 and comparison of CPU time and (b) the comparisons of uðtÞHet,y with different step 1 sizes (equivalent model).

¼ Tn

10

0

u(t)1

Equ,y

-2

(10 m) u(t)1

3

bCm1 7

12

1000

4

-1

14

2 0

2

Case 2

16

14

Error(t) (%)

Error(t) (%)

18

Case 1

16

200

1000

Case 2

400 600 Time (day)

800

1000

Fig. 5. Error(t) curve and the comparison of uðtÞHet,y of Case 2: (a) the Error(t) curve 1 and comparison of CPU time and (b) the comparisons of uðtÞHet,y with different step 1 sizes (equivalent model).

2

2

2

j r j2 bCm1 ¼ T nbZr1 ðn1Þðn þ 1ÞðEj2 þEj1 Þð2b n2 Er2 8 2 þ c E2 E2 3c E2 2

5bcEr2 Ej2 2b Er2 2 Þ

ð70Þ

n

ð69Þ

¼ T nbEj1 Ej2 Er2 Zr1 ðbn2 Ej2 þ bEr2 bn2 Er2 þ cEj2 Þ bCm1 9

ð71Þ

880

Y. Ren, H. Yang / Finite Elements in Analysis and Design 46 (2010) 875–888 j j j r r 2 j r 2 r bCm1 10 ¼ T cbE2 Z1 ðn1Þðn þ1ÞðE2 þ E1 Þðbn E2 þ bE2 bn E2 þ cE2 Þ

18

Error(t) (%)

ð72Þ

Case 3

16 14

j j r r j 2 r 2 j r bCm1 11 ¼ T cnbE1 E2 E2 Z1 ðbn E2 bn E2 bE2 cE2 Þ

12

ð73Þ

j j j r r 2 r 2 j r bCm1 12 ¼ 2T cnbE2 Z1 ðn1Þðn þ 1ÞðE2 þ E1 Þðbn E2 bn E2 bE2 cE2 Þ

10 8

ð74Þ

6

CPU Time

4

CPU Time

Equ Het

=

357. 683 ( s)

= 50975. 921 ( s)

2 0

200

400

600

800

1000

Time (day)

2

2

2

bCm2 ¼ T cnEj2 Er2 Er1 Zj1 ðb Er2 2cbEr2 þ b n2 Er2 c Ej2 þ cbEj2 Þ 1

ð75Þ

bCm2 ¼ T ncbEj2 Zj1 ðn1Þðn þ 1ÞðEr2 þ Er1 ÞðbEr2 þ cEj2 cEr2 Þ 2

ð76Þ

2

2

¼ T cEj2 Er2 Er1 Zj1 ðcbEr2 þ cEr2 bn2 2Er2 b n2 Ej2 cbn2 þEj2 c Þ bCm2 3 ð77Þ

3 Step=5(day) Step=1(day) Step=0.1(day)

Case 3

¼ T cEj2 j1 ð 1Þð þ 1ÞðEr2 þ Er1 Þð2cb 2 Er2 2cb 2 Ej2 2 2 3 2 b Er2 þ c Ej2 þ cbEr2 Þ

Z n

n

n

n

n

¼ T ncbEr2 Ej1 Zr1 ðbn2 Er2 bn2 Ej1 bEr2 cEj2 Þ bCm2 5

ð79Þ

-1

bCm2 ¼0 6

ð80Þ

-2

bCm2 ¼ T cbEj1 Ej2 Er2 Zr1 ðbn2 Ej2 þ bEr2 bn2 Er2 þ cEj2 Þ 7

ð81Þ

¼ T cbEr2 Zr1 ðn1Þðn þ1ÞðEj2 þ Ej1 Þðbn2 Ej2 þ bEr2 bn2 Er2 þ Ej2 Þ bCm2 8

ð82Þ

Equ,y

-2

(10 m)

ð78Þ

1

u(t)1

2

bCm2 4

0

-3

0

200

400 600 Time (day)

800

1000

2

Fig. 6. Error(t) curve and the comparison of uðtÞHet,y of Case 3: (a) the Error(t) curve 1 and comparison of CPU time and (b) the comparisons of uðtÞHet,y with different step 1 sizes (equivalent model).

2

2

2

2 r2 r j r2 bCm2 ¼ T nbEj1 Ej2 Zr1 ðc Ej2 Er2 2c Ej2 9 2 þ b n E2 3bcE2 E2 b E2 Þ

ð83Þ j r j j j r r bCm2 10 ¼ T cnbE2 Z1 ðn1Þðn þ 1ÞðE2 þ E1 ÞðcE2 bE2 cE2 Þ 2

2

ð84Þ

2

2

j j r j r j r j2 j2 2 2 r2 2 2 bCm2 11 ¼ T bE1 E2 Z1 ðcn bE2 E2 b n E2 n c E2 E2 þ c E2 þ n c E2 2

r j þ b Er2 2 þ2bcE2 E2 Þ

18

18

Case 4

16

10 8 6

CPU Time

4

CPU Time

Equ Het

=

357 . 683 (s)

Error(t) (%)

Error(t) (%)

14

12

= 142417 . 062 (s)

0

200

400 600 Time (day)

800

8 CPU Time

6

1000

CPU Time 0

-2

Het

=

357 . 683 (s)

= 150336 . 406 (s)

200

400 600 Time (day)

Step=5(day) Step=1(day) Step=0.1(day)

(10 m)

-2

1

Equ,y

2

-1

Equ

800

1000

3

Case 4

0

u(t) 1

-2

(10 m) u(t)1

Equ,y

2

10

2

3 Step=5(day) Step=1(day) Step=0.1(day)

12

4

2

0

Case 5

16

14

1

ð85Þ

-1

case 5

-2

-3 0

200

400

600

800

1000

Time (day) Fig. 7. Error(t) curve and the comparison of uðtÞHet,y of Case 4: (a) the Error(t) curve 1 and comparison of CPU time and (b) the comparisons of uðtÞHet,y with different step 1 sizes (equivalent model).

-3 0

200

400 600 Time (day)

800

1000

Fig. 8. Error(t) curve and the comparison of uðtÞHet,y of Case 5: (a) the Error(t) curve 1 and comparison of CPU time and (b) the comparisons of uðtÞHet,y with different step 1 sizes (equivalent model).

Y. Ren, H. Yang / Finite Elements in Analysis and Design 46 (2010) 875–888 2

2

2

j j j r r 2 j2 2 2 r2 bCm2 12 ¼ T bZ1 ðn1Þðn þ 1ÞðE2 þE1 Þð2c n E2 2n c E2 E2 b n E2

þ 2cn

2

2 j2 r j bEr2 Ej2 þ b2 Er2 2 þ c E2 þ 2bcE2 E2 Þ

The relationship between strain and displacement is ð86Þ

feðtÞg ¼ AT fuðtÞg

ForL the plane strain problem, EL1 , EL2 , n need to be replaced by EL1 E2 , , n , respectively. 2 1-n 1-n2 1n At the first time interval   fs0 g ¼ DH fe0 g ð87Þ (

At the (k+1)th time intervals P fe0 gðk þ 1Þ ¼ i ¼ 0 fei gk P , k ¼ 1, 2, 3, . . . fs0 gðk þ 1Þ ¼ i ¼ 0 fsi gk

where subscript k+ 1 and k refer to (k+ 1)th and kth time intervals, respectively.

The boundary condition is specified by fngfsðtÞg ¼ fpðtÞg, ~ fuðtÞg ¼ fuðtÞg,

4. FE based recursive equivalent field analysis

fng ¼

2275 KN/m

ð92Þ

x, y A Gu

nx

0

ny

ð93Þ

0

ny

nx

!

974.8K N/m

ð89Þ

50 m

x, yA Gs

where

It is assumed that the equilibrium equation in the equivalent field is

974.8 KN/m

ð90Þ

where {f(t)} represents a vector of body force, {u(t)} designates a vector of displacement 0 1 @ @ @ B @x @z @y C B C B @ @ @ C B C ð91Þ A¼B @y @z @x C B C B C @ @ @ @ A @z @y @x

ð88Þ

AfsðtÞg þff ðtÞg ¼ 0

881

ð94Þ

X 1 2

8

Y

45.0° 7

50 m 50 m

45.0° 1167.5 KN/m

1167.5 KN/m

X

50 m

6

3 Y

4 5

5m

Fig. 9. Schematic illustration of the structures concerned with Cases 6 and 7: (a) the structures concerned with Cases 6 and 7 and (b) the locations of sample points.

Rock (cyan): visco88

Joint (magenta): visco88

Joint (magenta): visco88

Rock (cyan): visco88

60° 30º

1.4 0.2

m .29

0

8-nodes isoparametric element

4m

9m

4m

1.4

Fig. 10. FE meshes of different models for Cases 6 and 7: (a) the equivalent model (Cases 6 and 7); (b) the heterogeneous model (Case 6) and (c) the heterogeneous model (Case 7).

882

Y. Ren, H. Yang / Finite Elements in Analysis and Design 46 (2010) 875–888

~ nx, ny denote the outside normals along Gs, {p(t)} and fuðtÞg are the prescribed functions, G ¼ Gu + Gs represents the whole boundary of the problem. According to the virtual work principle, we have [42] Z Z Z dfegT fsðtÞgdv ¼ dfugT ff ðtÞgdv þ dfugT fpðtÞgdG ð95Þ V

V

dfuðtÞg ¼ ½Ndfue ðtÞg

ð97Þ

where [N] represents a matrix of shape functions, fue ðtÞg and dfue ðtÞg represent nodal vector of {u(t)} and d{u(t)} at element level. Substituting Eq. (96) into Eq. (90) gives feðtÞg ¼ AT ½Nfue ðtÞg ¼ ½Bfue ðtÞg

m¼0

X

fuem gsm

ð100Þ

m¼0

X

fpðtÞg ¼

dfue gT

Z

e

Ve

e

Ve

e

½Bfsm gdv ¼

Ve

e

ð102Þ

½Kfum g ¼ fFm g,

Equ

CPU Time

6

½K ¼

XZ

Het

CPU Time

m ¼ 0, 1, 2, . . .

ð106Þ

fFm g ¼

½BT ½DH ½Bdv

ð107Þ

8P R R < e V ½NT ff0 gdv þ G ½NT fp0 gdG; e e P R P R P R :  e V ½BfCm gdv þ e V ½NT ffm gdv þ e G ½NT fpm gdG, e e e

m¼0 m ¼ 1, 2, 3, . . .

ð108Þ At the first time interval fu0 g is provided by the solution of Eq. (106) with m¼0 At the (k+ 1)th time interval interval fu0 g is given by X fui gk , k ¼ 1, 2, 3, . . . ð109Þ fu0 gk þ 1 ¼ The increase of m will gradually adapt computing accuracy to a criteria specified by ! m 1 X uki sk rx ð110Þ abs umi sm = k¼0

s¼1

Case 7

16 12 10 8

951. 310 (s)

6

= 155332. 760 (s)

4

=

2 0

200

400 600 Time (day)

800

1000

CPU Time CPU Time 0

200

Equ Het

=

951. 310 (s)

= 157447. 878 (s)

400 600 Time (day)

800

1000

3 Case 6

-1

-2

(10 m)

1

Equ,y

0

Step=5(day) Step=1(day) Step=0.1(day)

2

u(t) 1

Step=5(day) Step=1(day) Step=0.1(day)

2 -2

½NT fpm gdG

Ge

e

Ve

e

3 (10 m)

Z

where

Error(t) (%)

Error(t) (%)

8 4

Equ,y

dfue gT

14

10

u(t) 1

Ve

X

Substituting Eqs. (48), (87) and (101) into Eq. (105) and considering the arbitrariness of dfue gT yield

Case 6

12

Case 7

0

-2

-2 -3

½NT ffm gdv þ

18

14

-1

Z

ð105Þ

Ge

16

1

dfue gT

e

18

2

X

i¼0

ð101Þ

Substituting Eq. (40) into Eq. (95) yields XZ XZ XZ dfegT fsm gdv ¼ dfugT ffm gdv þ dfugT fpm gdG

ð104Þ

Substituting Eq. (98) into Eq. (102) gives X

Substituting Eqs. (100) and (41) into Eq. (98) gives fem g ¼ ½Bfuem g

fpm gsm

m¼0

ð98Þ

Let fuðtÞg represents the general nodal vector of {u(t)}, and expand both fuðtÞg and fue ðtÞg in the form X fum gsm ð99Þ fuðtÞg ¼

ð103Þ

m¼0

G

where d{e} and d{u} stand for virtual strain and displacements, respectively. In the implementation of FEM, the {u(t)} and d{u(t)} are approximated by   fuðtÞg ¼ ½N ue ðtÞ ð96Þ

fue ðtÞg ¼

where {fm} and {pm} are coefficient vectors of X ff ðtÞg ¼ ffm gsm

0

200

400 600 Time (day)

800

1000

Fig. 11. Error(t) curve and the comparison of uðtÞHet,y of Case 6: (a) the Error(t) 1 curve and comparison of CPU time and (b) the comparisons of uðtÞHet,y with 1 different step sizes (equivalent model).

-3

0

200

400 600 Time (day)

800

1000

Fig. 12. Error(t) curve and the comparison of uðtÞHet,y of Case 7: (a) the Error(t) 1 curve and comparison of CPU time and (b) the comparisons of uðtÞHet,y with 1 different step sizes (equivalent model).

Y. Ren, H. Yang / Finite Elements in Analysis and Design 46 (2010) 875–888

where x is a prescribed error tolerance, umi stands for the i-th component of fum g and i ¼1, 2, 3, y Eq. (110) is checked when each fum g is obtained, if it is satisfied continually for three times, the computation at the current time interval will stop, and then step into the next one. If the global coordinate system X Y is not identical with the local coordinate system x  y (Fig. 1), a transformation of stress and strain is required. We assume that the equivalent stress and strain in the X  Y system are defined by fsuðtÞg ¼ sux ðtÞ

euy ðtÞ guðtÞ

fsuðtÞg ¼ ½TfsðtÞg

ð113Þ

feuðtÞg ¼ ½TT feðtÞg

ð114Þ

By expanding {s (t)} and {e (t)} in the form X fsuðtÞg ¼ fsum gsm 0

0

ð115Þ

m¼0

feuðtÞg ¼

X

feum gsm

ð116Þ

m¼0

ð111Þ

and substituting Eqs. (115), (40) and (116), (41) into Eqs. (113) and (114) yields

T

ð112Þ

974.8 KN/m

 feuðtÞg ¼ eux ðtÞ

suy ðtÞ tuðtÞ

T

We have [43]

50 m

2275 KN/m

fsum g ¼ ½Tfsm g

ð117Þ

974.8 KN/m



883

Y 1

Y

8

2

50 m

50 m

5m

45º

7

3 X

5m

X 1167.5 KN/m

1167.5 KN/m

50 m

45º

6

4

5

Fig. 13. Schematic illustration of structure concerned with Case 8: (a) the structure concerned with Case 8 and (b) the locations of sample points.

Joint (magenta): visco88 Rock (cyan): visco88

9m

0.2 8-nodes isoparametric element

4m

1.4

Fig. 14. FE meshes of different models for Case 8: (a) the equivalent model and (b) the heterogeneous model.

884

Y. Ren, H. Yang / Finite Elements in Analysis and Design 46 (2010) 875–888

feum g ¼ ½TT fem g

ð118Þ (

Substituting Eqs. (117) and (118) into Eq. (48) gives fsum g ¼ ½DuH feum g þ fCum g,

m ¼ 1, 2, 3, . . .

ð119Þ

Eq. (119) stands for a recursive constitutive equation in the global coordinate system X Ywhere ½DuH  ¼ ½T1 ½DH ½T T

ð120Þ

fCum g ¼ ½T1 fCm g

ð121Þ

0 B ½T ¼ @

sin2 a

sin2 a

cos2 a

sin a cos y

For FE formulation in the global coordinate system X Y, [DH] and {Cm} in Eq. (106) needs to be replaced by [D0 H] and {C0 m}, respectively. Eight-node isoparametric elements are used in this paper [42].

sin a cos a

1

2 sin a cos a C 2 sin a cos a A

ð122Þ

2

cos2 asin a

At the first time interval fs00 g ¼ ½DuH feu0 g

ð123Þ

18 Case 8

16 14 Error(t) (%)

ð124Þ

5. Numerical verification

cos2 a

This section presents ten 2-D numerical simulations of fictitious tunnels in the orthogonal viscoelastic jointed rock. The results obtained by the proposed equivalent model are compared with those given by an FE based heterogeneous analysis via ANSYS, in terms of computing accuracy and cost. The comparisons are carried out on a PC with a 2.4 GHz AMD Athlon processor and 2 GB of memory. ‘CPU time’ referred in this section does not include the expense on an FE pre-processing. All the details relevant to constitutive models, and computing parameters, etc. are listed in Tables 1 and 2, for all the numerical tests, x ¼0.001. The relative deviation between two models are estimated by ErrorðtÞ ¼ , , !

 

  

 

 Int P



Equ, x y y

y x

uðtÞiHet, x þ uðtÞEqu, uðtÞHet, uðtÞHet, uðtÞHet,



uðtÞi i i i i

12 10

i¼1

8 6

CPU Time

4

CPU Time

2

At the (k+ 1)th time intervals P feu0 gðk þ 1Þ ¼ i ¼ 0 feui gk P fsu0 gðk þ 1Þ ¼ i ¼ 0 fsui gk

0

200

Equ Het

400

 100%

Int

655 . 469 (s)

=

ð125Þ

= 160791 . 028 (s) 600

800

where u(t)i is displacement of sample point i, Int stands for the number of sample points, the superscripts Equ and Het refer to equivalent and heterogeneous models, respectively, and the superscripts x and y refer to x and y directions, respectively.

1000

Time (day) 3 Step=5(day) Step=1(day) Step=0.1(day)

(10 m)

-2

1

Equ,y

0

u(t) 1

2

-1

Case 8

5.1. Different combinations constitutive models

-2 -3

0

200

400 600 Time (day)

800

1000

974.8 KN/m

2275 KN/m

50 m

Y

100 m

1167.5 KN/m

100 m

Y

1 8

45.0° 3 45.0° 6 R5m

9

4 5

45.0° 11 45.0°

15 14 R5m

15 m

10

16

2

7

X

50 m

1167.5 KN/m

974.8 KN/m

Fig. 15. Error(t) curve and the comparison of uðtÞHet,y of Case 8: (a) the Error(t) 1 curve and comparison of CPU time and (b) the comparisons of uðtÞHet,y with 1 different step sizes (equivalent model).

Cases 1–3 consider three kinds of different combinations of constitutive models of joint and rock for the structure shown in Fig. 2a, where only a half of the structure is outlined because the problem is symmetric. The locations of five sample points are labeled in Fig. 2b, and descriptions of FE meshes for two models (equivalent and heterogeneous models) are exhibited in Fig. 3a and b. Figs. 4a, 5a and 6a give three t Error(t) curves and CPU time for two models. For the Cases 1–3, Error(0) are 4.16%, 4.47%, and 4.66%, respectively, and Error(1000) are 13.77%, 13.60%, and 13.81%, respectively. The ratios of CPU time are about 108, 120, and 143.

X

12 13

15 m

Fig. 16. Schematic illustration of structure concerned with Case 9: (a) the structure concerned with Case 9 and (b) the locations of sample points.

Y. Ren, H. Yang / Finite Elements in Analysis and Design 46 (2010) 875–888

885

Joint (magenta): visco88

Rock (cyan): visco88

8m

2.8

8m

0.5

8-nodes isoparametric element

Fig. 17. FE meshes of different models for Case 9: (a) the equivalent model and (b) the heterogeneous model.

Figs. 4b, 5b, and 6b give three comparisons of the solutions obtained via different step sizes for an equivalent model.

18 Case 9

16 14

Cases 4 and 5 consider the impact of different b on the solution for the structure of Case 3. Fig. 3 gives descriptions of FE meshes for two models. Figs. 7a and 8a give two t  Error(t) curves and CPU time for two models. For the Cases 4 and 5, Error(0) are 4.97% and 4.91%, respectively, and Error(1000) are 13.66% and 13.60%, respectively. However, the ratios of CPU time are about 398 and 420, respectively. Figs. 7b and 8b give two comparisons of the solutions obtained via different step sizes for an equivalent model.

Error(t) (%)

5.2. Different relative width of joint (b)

12 10 8 6

CPU Time

4 2

CPU Time 0

Cases 8–10 consider more complicated geometry, as shown in Figs. 13a, 16a and 19a. The locations of sample points are labeled in Figs. 13b, 16b and 19b. Figs. 14, 17 and 20 show descriptions of FE meshes for two models. Figs. 15a, 18a and 21a give three t  Error(t) curves and CPU time for two models. For Cases 8–10, Error(0) are 4.67%, 5.52% and 6.23%, respectively, and Error(1000) are 13.57%, 14.20% and

(10 m)

Equ,y

-2

5.4. Different geometry of tunnel

400

=

1408 . 094 (s)

= 172825 . 086 (s) 600

800

1000

3

u(t) 1

Cases 6 and 7 consider the impact of different a on the solution, as shown in Fig. 9a, where a whole structure is outlined because the problem is not symmetric, which is different with Cases 1–5. The locations of eight sample points are labeled in Fig. 9b. Fig. 10 shows descriptions of FE meshes for two models. Figs. 11a and 12a give two t  Error(t) curves and CPU time for two models. For Cases 6 and 7, Error(0) are 5.01% and 5.22%, respectively, and Error(1000) are 13.90% and 13.91%, respectively. The ratios of CPU time are about 163 and 165. Figs. 11b and 12b give two comparisons of the solutions obtained via different step sizes for an equivalent model.

Het

Time (day) Step=5(day) Step=1(day) Step=0.1(day)

2 5.3. Different dip angle of joint (a)

200

Equ

1

Case 9

0 -1 -2 -3 -4

0

200

400

600

800

1000

Time (day) Fig. 18. Error(t) curve and the comparison of uðtÞHet,y of Case 9: (a) the Error(t) 1 curve and comparison of CPU time and (b) the comparisons of uðtÞHet,y with 1 different step sizes (equivalent model).

14.91%, respectively. The ratios of CPU time are about 245, 123 and 32. (Figs. 16–20) Figs. 15b, 18b and 21b give three comparisons of the solutions obtained via different step sizes for an equivalent model. From Cases 1 to 10, the comparison of an equivalent analysis based solutions obtained via different step sizes is only given at one sample point for each Case, actually similar comparison can be observed at other sample points, but not given here due to the capacity restriction of a paper.

Y. Ren, H. Yang / Finite Elements in Analysis and Design 46 (2010) 875–888

974.8 KN/m

50 m

2275 KN/m

X

amin ¼

amax ¼

Dr33



Dj33

1167.5 KN/m

125 m

125 m

1  



ð126Þ

1     1b=Dr33

ð127Þ

b=Dj33  1b=Dr33



b=Dj33

when a1 is fixed

a2 min ¼

! Dj 1 a1 b Dj33 a1 b  ¼ 33  r j D33 D 1b Dr33 1b 33

Main tunnel

Cross passage 1

X

18

2

R 4m 3

16

15

R5m

Y

50 m

1167.5 KN/m

974.8 KN/m

1. In comparison with the FE based heterogeneous analysis via ANSYS, the proposed equivalent model can provide a very fast numerical approximation with a reasonable computing accuracy. 2. The computing accuracy of proposed approach can be maintained via an adaptive process when the step size varies, because the criteria controlling the truncation error at a discretized time interval does not vary (Figs. 4b, 5b, 6b, 7b, 8b, 11b, 12b, 15b, 18b and 21b). 3. The heterogeneous model incorporates detailed joint geometries and the corresponding FE pre-processing is not straightforward by an ANSYS GUI (graphical user interface), and becomes more tough and time-consuming when the density of joints is relatively high. Even by using a ANSYS Parametric Design Language (APDL) based FE pre-processor we created, the ratio of CPU Time of the FE pre-processing between the heterogeneous and equivalent models still ranges 3273–7234 (Table 3). On the other hand, when a and b of a structure change, the heterogeneous model needs a new mesh generation, but the equivalent model does not. 4. The major reason inducing the computing inaccuracy of the equivalent analysis mainly attributes to (1) The adopted equivalent assumption that is simple and easy to implement, is not rigorous enough both physically and mathematically.

(2) The density of joints in the computing domain, which is closely concerned with the accuracy of an equivalent analysis. The higher the density, the more accurate the results (Fig. 22). As the increase of density of joints in the computing domain, we expect a better comparison with the FE-based heterogeneous analysis that unfortunately cannot work on the PC when number of joints is more than 80  40 because of the increase of unknowns in addition to the unbearable FE pre-processing and computing expense. (3) The selection of a1 and a2, which is an important factor relevant to the results. when a1 ¼ a2 ¼1.0, 1.5, 1.9, 2.0, 2.1, 2.2, 2.3 Error(0) are 13.50%, 7.33%, 5.74%, 4.16%, 6.09%, 7.42%, 8.13%, respectively (Case 1). In the all above presented numerical tests, we take a1 ¼ a2 ¼2. r Physically there is Dr33 rDH 33 rD33 in general, so that we can estimate the range of a1 and a2 for certain cases.When a1 ¼ a2

17

ð128Þ

Main tunnel

Service tunnel 7

4 5

6

11

14 13 12

8

10 3.5 m

6. Computing remarks and discussions

4m

886

Y 9

R 5m

4m

4m 20 m

20 m

Fig. 19. Schematic illustration of structure concerned with Case 10: (a) the structure concerned with Case 10 and (b) the locations of sample points.

Joint (magenta): visco88

Rock (cyan): visco88

0m

3.6

2m

0.7

8-nodes isoparametric element

Fig. 20. FE meshes of different models for Case 10: (a) the equivalent model and (b) the heterogeneous model.

Y. Ren, H. Yang / Finite Elements in Analysis and Design 46 (2010) 875–888

It seems a bit difficult to select the best pair of a1 and a2 within amin, amax at the present stage. Further detailed work related with

18 Case 10

16

this issue is required with more physical/mathematic and experimental supports.

14 Error(t) (%)

887

12 10

7. Conclusions

8

Equ

6

CPU Time

4

CPU Time

2

Het

0

200

=

5826. 172 (s)

The major merits of this paper include

= 189318. 566 (s)

400 600 Time (day)

800

1. By utilizing a piecewised expanding technique at a discretized time interval, a recursive equivalent constitutive model is presented for an orthogonal viscoelastic jointed rock. 2. Based on the recursive equivalent constitutive model and FEM, an adaptive numerical model is developed to simulate equivalent viscoelastic field, consequently a reasonable, convenient and economic numerical simulation can be carried out. Sufficient computing accuracy can be maintained by the proposed adaptive algorithm for the different choices of step size, but it may not be the case for other temporal algorithms. [44] 3. Ten numerical tests are presented to verify the proposed model, in the viewpoint of balance between computing accuracy and cost, the results are fairly good, and may be considered as a primary estimation of a field problem.

1000

3 Step=5(day) Step=1(day) Step=0.1(day)

(10 m)

2 -2

1

Case 10

Equ,y

0

u(t) 1

-1 -2 -3 -4

0

200

400 600 Time (day)

800

1000

Fig. 21. Error(t) curve and the comparison of uðtÞHet,y of Case10: (a) the Error(t) 1 curve and comparison of CPU time and (b) the comparison of uðtÞHet,y with 1 different step sizes (equivalent model). Table 3 Comparison of CPU time of pre-processing for two models. Case no.

Equivalent model (s)

Heterogeneous model (s)

1 2 3 4 5 6 7 8 9 10

1.579 1.579 1.579 1.579 1.579 2.281 2.281 1.828 2.984 4.188

5367.938 5367.938 5367.938 10,389.516 11,422.032 11,993.134 12,107.354 12,335.760 12,906.897 13,706.439

An equivalent assumption adopted in this paper is relatively simple, and leads to a relatively simple mathematical derivation and an economic computing expense with reasonable numerical result. However, it is definitely not rigorous enough mathematically and physically, and results in computing inaccuracy. An asymptotic theory [45,46] based equivalent model of visocelastic jointed rock, which is rigorous physically and mathematically, is currently being investigated by authors. Some primary results indicate that model seems much more expensive computationally than the proposed model in this paper.

Acknowledgements The research leading to this paper is funded by NSF (10421002), NSF (10472019), NSF (10772035), NSF (10721062), NKBRSF [2005CB321704, 2010CB832703], and the fund of disciplined leaders of young and middle age faculty in colleges of Liaoning Province.

40

References

Error (1000) (%)

35 30 25 20 15 10

20×10 30×15 40×20 50×25 60×30 70×35 80×40 Number of joints

Fig. 22. The variation of Error(1000) with the increase of joints in the computing domain (Case 1).

a2 max ¼

1 Dj33

!



a1 b Dj33

Dj33 a1 b ¼ 1 1b 1b

Similar estimation can be obtained when a2 is fixed.

ð129Þ

[1] S. Maghous, P. de Buhan, L. Dormieux, Non-linear global elastic behaviour of a periodically jointed material, Mech. Res. Commun. 29 (2002) 45–51. [2] H. Yoshida, H. Horii, Micromechanics-based continuum model for a jointed rock mass and excavation analyses of a large-scale cavern, Int. J. Rock Mech. Min. Sci. 41 (2004) 119–145. [3] T.G. Sitharam, J. Sridevi, N. Shimizu, Practical equivalent continuum characterization of jointed rock masses, Int. J. Rock Mech. Min. Sci. 38 (2001) 437–448. [4] C. Esposito, S. Martino, G. Scarascia, Mountain slope deformations along thrust fronts in jointed limestone : an equivalent continuum modelling approach, Geomorphology 90 (2007) 55–72. [5] C.y. Ku, J.s. Lin, J.c. Chern, Modeling of jointed rock masses using a combined equivalent-continuum and discrete approach, Int. J. Rock Mech. Min. Sci. 41 (2004) 2004. [6] J.s. Lin, C.y. Ku, Two-scale modeling of jointed rock masses, Int. J. Rock Mech. Min. Sci. 43 (2006) 426–436. [7] W.G. Pariseau, S. Puri, S.C. Schmelter, A new model for effects of impersistent joint sets on rock slope stability, Int. J. Rock Mech. Min. Sci. 45 (2008) 122–131. [8] M. Cai, A constitutive model for highly jointed rock masses, Mech. Mater. 13 (1992) 217–246. [9] E.P. Chen, Constitutive model for jointed rock mass with orthogonal sets of joints, J. Appl. Mech. 56 (1989) 25–32. [10] C.M. Gerrard, Elastic models of rock masses having one, two and three sets of joints, Int. J. Rock Mech. Min. Sci. Geomech. Abstr. 19 (1982) 15–23.

888

Y. Ren, H. Yang / Finite Elements in Analysis and Design 46 (2010) 875–888

[11] T.H. Huang, C.S. Chang, Z.Y. Yang, Elastic-moduli for fractured rock mass, Rock Mech. Rock Eng. 28 (1995) 135–144. [12] B. Singh, Continuum characterization of jointed rock masses: Part i—the constitutive equations, Int. J. Rock Mech. Min. Sci. Geomech. Abstr. 10 (1973) 311–335. [13] B. Singh, Continuum characterization of jointed rock masses: Part ii—significance of low shear modulus, Int. J. Rock Mech. Min. Sci. Geomech. Abstr. 10 (1973) 337–349. [14] J. Sridevi, T.G. Sitharam, Analysis of strength and moduli of jointed rocks, Geotech. Geol. Eng. 18 (2000) 3–21. [15] H.T. Yang, R. Wu, G. Wang, X. Li, Finite element simulating computation for joint-rock based on composite constitutive models, Chin. J. Geotech. Eng. 18 (1996) 69–76. [16] O.C. Zienkiewicz, D.W. Kelly, P. Bettess, The coupling of the finite element method and boundary solution procedures, Int. J. Numer. Methods Eng. 11 (1976) 355–375. [17] W. Zhang, H.X. Zhang, Elastic models of jointed rock masses, Chin. J. Geotech. Eng. 9 (1987) 33–44. [18] D.P. Adhikary, A.V. Dyskin, Continuum model of layered rock masses with non-associative joint plasticity, Int. J. Numer. Anal. Methods Geomech. 22 (1998) 245–261. [19] S.h. Chen, P. Egger, Three dimensional elasto-viscoplastic finite element analysis of reinforced rock masses and its application, Int. J. Numer. Anal. Methods Geomech. 78 (1999) 61–78. [20] S.h. Chen, C.H. Fu, S. Isam, Finite element analysis of jointed rock masses reinforced by fully grouted bolts and shotcrete lining, Int. J. Rock Mech. Min. Sci. 46 (2009) 19–30. [21] J. Aboudi, Micromechanical characterization of the non-linear viscoelastic behavior of resin matrix composites, Compos. Sci. Technol. 38 (1990) 371–386. [22] J. Aboudi, Micromechanically established constitutive equations for multiphase materials with viscoelastic–viscoplastic phases, Mech. Time-Dependent Mater. 9 (2005) 121–145. [23] R. Haj-ali, J. Aboudi, International journal of solids and structures nonlinear micromechanical formulation of the high fidelity generalized method of cells, Int. J. Solids Struct. 46 (2009) 2577–2592. [24] R.M. Haj-ali, A.H. Muliana, Micromechanical constitutive framework for the nonlinear viscoelastic behavior of pultruded composite materials, Int. J. Solids Struct. 40 (2003) 1037–1057. [25] C. Mareau, V. Favier, M. Berveiller, Micromechanical modeling coupling timeindependent and time-dependent behaviors for heterogeneous materials, Int. J. Solids Struct. 46 (2009) 223–237. [26] A. Muliana, R. Haj–Ali, A micromechanical model for the nonlinear viscoelastic behavior of laminated composites, in: Proc 16th ASCE Environ. Eng. Conf. (2003) 1–6. [27] A.H. Muliana, R.M. Haj-ali, Nested nonlinear viscoelastic and micromechanical models for the analysis of pultruded composite materials and structures, Mech. Mater. 36 (2004) 1087–1110.

[28] O. Pierard, J. Llorca, J. Segurado, I. Doghri, Micromechanics of particlereinforced elasto-viscoplastic composites: finite element simulations versus affine homogenization, Int. J. Plast. 23 (2007) 1041–1060. [29] S. Maghous, G.J. Creus, Periodic homogenization in thermoviscoelasticity: case of multilayered media with ageing, Int. J. Solids Struct. 40 (2003) 851–870. [30] P.M. Suquet, E. Sanchez-Palencia, A. Zaoui, Homogenization Techniques for Composite Media, Springer, Berlin/Heidelberg, 1987. [31] W. Chen, J. Fish, A dispersive model for wave propagation in periodic heterogeneous media based on homogenization with multiple spatial and temporal scales, J. Appl. Mech. 68 (2001) 153–161. [32] S. De, J. Fish, M.S. Shephard, S.K. Kumar, Multiscale modeling of polymer rheology, Phys. Rev. E. 74 (2006) 1–4. [33] H. Waisman, J. Fish, A heterogeneous space–time full approximation storage multilevel method for molecular dynamics simulations, Int. J. Numer. Methods Eng. 73 (2008) 407–426. [34] Q. Yu, J. Fish, Multiscale asymptotic homogenization for multiphysics problems with multiple spatial and temporal scales: a coupled thermoviscoelastic example problem, Int. J. Solids Struct. 39 (2002) 6429–6452. [35] J. Aboudi, The generalized method of cells and high-fidelity generalized method of cells micromechanical models—a review, Mech. Adv. Mater. Struct. 11 (2004) 329–366. [36] L.C. Brinson, W.G. Knauss, Thermorheologically complex behavior multiphase viscoelastic materials, J. Mech. Phys. Solids 39 (1991) 859–880. [37] J. Goering, Analysis of thermoviscoelastic behavior of unidirectional fiber composites, Compos. Sci. Technol. 29 (1987) 103–131. [38] S.t. Liu, C.y. Chang, H.t. Yang, Prediction of visco-elastic property of unidirectionally jointed rock, Chin. J. Rock. Mech. Eng. 22 (2003) 582–588. [39] S.T. Liu, K.Z. Chen, X.A. Feng, Prediction of viscoelastic property of layered materials, Int. J. Solids Struct. 41 (2004) 3675–3688. [40] V.F.P. Dutra, S. Maghous, A. Campos, A.R. Pacheco, A micromechanical approach to elastic and viscoelastic properties of fiber reinforced concrete, Cem. Concr. Res. 40 (2010) 460–472. [41] R.M. Christensen, Theory of Viscoelasticity, second ed., Academic Press, New York, 1982. [42] O.C. Zienkiewicz, R.L. Taylor, The Finite Element Method, fifth ed., Butterworth-Heinemann, Boston, 2000. [43] R.M. Jones, Mechanics of Composite Materials, second ed., Taylor & Francis, Philadelphia, 1999. [44] H.T. Yang, A precise algorithm in the time domain to solve the problem of heat transfer, Numer. Heat Tr. B—Fund 35 (1999) 243–249. [45] B. Hassani, E. Hinton, A review of homogenization and topology optimization i—homogenization theory for media with periodic structure, Comput. Struct. 69 (1998) 707–717. [46] B. Hassani, E. Hinton, A review of homogenization and topology optimization ii—analytical and numerical solution of homogenization equations, Comput. Struct. 69 (1998) 719–738.