The 3rd International Conference on ″Computational Mechanics and Virtual Engineering″″ COMEC 2009 29 – 30 OCTOBER 2009, Brasov, Romania

NUMERICAL STUDY OF THE BEHAVIOUR FOR ELASTICVISCOPLASTIC ROCK AROUND CIRCULAR GALLERIES Iuliana Paraschiv-Munteanu1 1

Faculty of Mathematics and Computer Science, University of Bucharest, Romania [email protected]

Abstract : The variation of stress during creep convergence of a horizontal circular galleries excavated in rock salt is studied. Examples are given for rock salt by N. Cristescu ([1], [2]). A non-associated elasto-viscoplastic constitutive equation is used to describe both compressibility and/or dilatancy during transient and steady-state creep, as well as evolutive damage possibly leading to failure. An in-house FEM numerical method and iterative method is used for this purpose ([4], [6]). The variation in time of radial convergence of the galleries walls and of the stress state will be illustrated by several figures. The significance of these variations for long-term stability is discussed. Numerical results are obtained using MATLAB ([8], [9]). Keywords: elastic-viscoplastic model, rock mechanics, numerical methods.

1. INTRODUCTION Study of stress distribution during creep of the rock surronding a circular horizontal tunnel is a very important problem, mainly for mining engineering. At big depths, an opening excavated in rock can close completely after time intervals which are of the order of several tens of years. For the design of underground cavities one must be able to predict quite accurately not only the stress and strain distribution around them, but also the apperance and possibly slow spreading in time of a microcraked domain produced just by the excavation. Since the microcraking is related to dilatancy, the irreversible volumetric changes, either dilatancy or compressibility have started to be studied too. If the stress is in the dilatancy domain damage by microcraking can develop steadily in time, ultimately leading to a major underground failure. That is why it is important to study the stress variation during creep when microcraks are also developing. In this paper we study the distribution of stresses, deformations and displacement around a circular cylindrical gallery. We tried to determine a numerical solution without using hypothesis that stress state remains constant in time as in case of simplified solution. For numerical solution we used the scheme proposed by Paraschiv-Munteanu ([3], [4], [5]) using finit element method for spatial integration and a complet implicit method for integration in time. In most cases we observed that in proximity of underground opening the stress becomes relaxed relatively to the moment of excavation. However, for short period of time the creep solution and the numerical solution are very close. The elasto-viscoplastic model that we consider in this paper is described by equation

ε& =

σ& R

W I (t ) ∂F 1  R ∂S  1 + − (σ ) + k S (σ ) ,  σ& 1 + k T 1 − 2G  3K 2G  H (σ ) ∂σ ∂σ

(1)

where K and G are elastic moduls, k T and k S are viscosity constants, H is the plasticity function, F is the viscoplastic potential for transitory creep, S is the viscoplastic potential for stationary creep. The constitutive equation (1) can describe the following mechanical properties exhibited by most rocks: transitory and stationary creep, work-hardening during transient creep, volumetric compressibility and/or dilatancy, as well as short-term failure. All these properties are incorporated into the constitutive equation via the procedure used to determine the constitutive functions (Cristescu [1], [2], Cristescu and Hunche [3]).

2. PROBLEM FORMULATION The stress distribution just after excavation is obtained by exact elastic solution. Let a the initial radius of the cavity and m ∈ N , m ≥ m0 ≥ 5 , number of radius which defined the limits of the domain, Ω = [a, ma ] × [0,2π ) . We assume that in all horizontal directions the primary stresses are the same, σ h , and the depth is sufficient great to consider that

495

σ v , the vertical initial (primary) stress, is not variable in the domain Ω ( σ h corresponds for the axis of the tunnel). The conditions for r → ∞ (in case of infinitely domain) has been considered on the external boundary of the domain Ω .

{

}

Proposition 1. If on the walls of the cavity, Γ = (a, θ ) θ ∈ [0, 2π ) , a pressure p is acting (due to various reasons 1 and which may be constant or variable): S (a, θ ) = p, σ S (a, θ ) = 0, ∀θ ∈ [0, 2π ) σ rr (2) rθ and on the external boundary of the domain Ω , Γ2 = { (ma, θ ) θ ∈ [0,2π ) } , we have:

1 1  S σ rr ( ma,θ ) = 2 (σ h +σ v ) + 2 (σ h − σ v )cos 2θ  1 1  S σ θθ (ma, θ ) = (σ h +σ v ) − (σ h − σ v )cos 2θ , 2 2  1  S σ rθ ( ma,θ ) = − 2 (σ h − σ v )sin 2θ  then the stress state just after excavation is: C  6C 4D  σ~rrS (r ,θ ) = 2 A1 + 21 +  − 2 A2 − 42 − 2 2  cos 2θ r r r   C  6C  S (r ,θ ) = 2 A1 − 21 +  2 A2 + 12 B2 r 2 + 42  cos 2θ σ~θθ r r   6C 2D   σ~rSθ (r ,θ ) =  2 A2 + 6 B2 r 2 − 42 − 2 2  sin 2θ r r   

(3)

(4)

 4 D2   cos 2θ  , 2 r     where A1 , C1 , A2 , B2 , C 2 , D2 are constants : 

σ~ zzS (r ,θ ) = ν 4 A1 + 12 B2 r 2 −

A1 = A2 =

m2 4(m 2 − 1)

(σ h + σ v ) −

m 2 ( m 4 + m 2 + 4) 6

4(m − 1) m 4 a 4 ( m 2 + 1)

1 2(m 2 − 1)

(σ h − σ v ) ,

B2=

m2a2  1   p − 2 (σ h + σ v ) , 2 ( m − 1)  

m2 2

2a ( m 6 − 1)

(σ h − σ v ) ,

m2a2

(σ h − σ v ) . 4( m 6 − 1) 2(m 2 − 1) Proof. The components of stress are obtained from equilibrum equation using the Airy function and the constants result from conditions (2) and (3). ■ C 2= −

(σ h − σ v ) ,

p , C1 =

D2=

From (4) it is easy to obtain: Proposition 2. The components of deformation corresponding stress state (5) are:  C1  6C 2 4 D2 1 +ν   2 ε~rr ( r ,θ ) = 2(1 − 2ν ) A1 + 2 + − 2 A2 − 12νB2 r − 4 − 2 (1 − ν ) cos 2θ  E  r r r   

  C1  6C 2 4 D2  2 2(1 − 2ν ) A1 − 2 + 2 A2 + 12(1 −ν ) B2 r + 4 + 2 ν  cos 2θ  , r r r     6 4 C D 1 + ν   ε~rθ ( r ,θ ) =  2 A2 + 6 B2 r 2 − 42 − 2 2  sin 2θ , E  r r  and the components of the displacement are:  2C 4 D2 C1  1 +ν   u~r (r ,θ ) = + − 2 A2 r − 4νB2 r 3 + 32 + (1 −ν ) cos 2θ  2(1 − 2ν ) A1r − E  r  r r  

ε~θθ (r ,θ ) =

1 +ν E

1 +ν u~θ ( r ,θ ) = 2E

4C 2 4 D2   3 4 A2 r + 4(3 − 2ν ) B2 r + 3 + r ( 2ν − 1)  sin 2θ . r  

(5)

(6)

Observation. It is easy to observe that in the case of infinitely domain, when m → ∞ , the stress, deformation and displacement components are the same like in papers of Cristescu and Paraschiv ([3], [4]), because, when m → ∞ , we have:

496

A1 →

1 (σ h + σ v ) , C1 → a 2  p − 1 (σ h + σ v ) , 4 2  

4 2 1 (σ h − σ v ) , B2 → 0 , C 2 → − a (σ h − σ v ) , D 2 → a (σ h − σ v ) . 4 4 2 So, this result proves in one way that taking m ≥ m0 ≥ 5 it is acceptable for moving the infinitely condition on the

A2 → −

boundary r = ma . Elastic solution is used as initial data for the integration in long time intervals using finite elements methods. The general formulation of the problem of determination the stress distribution around a circular horizontal tunnel in elastoviscoplastic rock, like a cvasistatic problem, is:

(

)

find the displacement function u r , uθ : R + × Ω → R 2 , the stress function σ : R + × Ω → S 3 and the irreversible stress work function W I : R + × Ω → R such that: R+ ×Ω Div σ R (t , (r , θ )) = 0 in (7)

σ& R = 2Gε& + ( 3K − 2G ) ε& 1 + kT 1 −

W I (t ) H(σ )

∂F   ( 2G − 3K ) ∂F 1 + 2G +  σ 3 ∂ ∂ σ  

∂S   ( 2G − 3K ) ∂S kS  1 + 2G ∂σ ∂σ  3  W I (t ) ∂F W& I = k T 1 − ⋅σ H (σ ) ∂σ

in

in

R+ × Ω

R+ × Ω

(9)

σ R (t , (a, θ )) = p − σ P (θ ) rr  rr , ∀ t > 0 , θ ∈ [0,2π ) σ rRθ (t , (a, θ )) = 0  u (t , ( ma, θ )) = 0 , ∀ t > 0 , θ ∈ [0,2π ) σ S (0, (r ,θ )) = σ P (θ ) + σ~ (r ,θ )

(10)

(u r , uθ )(0, (r ,θ )) = (u~r , u~θ )(r ,θ ) , ∀ (r ,θ ) ∈ Ω W (0, (r , θ )) = H (σ I

P

(8)

(11)

(θ ))

where σ~ and (u~r , u~θ ) are the stress and, respectively, the displacement corresponding for the moment of excavation 1 1 and σ rrP (θ ) = (σ h + σ v ) + (σ h − σ v ) cos 2θ . 2 2

3. THE NUMERICAL APPROACH For the problem (7)-(11) we determine a numerical solution based on some results presented by Rosca and Sofonea [7] using a complete implicit method for integration in time (see Paraschiv-Munteanu [5]). If (u , σ R , W I ) , where u = (u r , uθ ) , is the solution of the problem (7)-(11) then we determine:

u = u − u~ , σ = σ R − σ R , such that u (t , ma, θ ) = 0 , ∀ t > 0 , θ ∈ [0,2π ) Div σ R (t , (r , θ )) = 0

in

(12)

R+ × Ω

(13)

σ (t , (a, θ )) n = o , ∀t > 0 , θ ∈ [0,2π ) . So, we have to solve the problem: find the displacement function u : R + × Ω → R 2 , the stress function σ : R + × Ω → S 3 and the irreversible stress work function W I : R + × Ω → R such that:

497

σ& R = 2Gε& (u& ) + (3 K − 2G ) ε& (u& )1 + kT 1 −

W I (t ) H (σ + σ~ + σ P )

∂F  ( 2G − 3K ) ∂F  (σ + σ~ + σ P ) 1 + 2G (σ + σ~ + σ P )  +  σ 3 ∂ ∂ σ  

∂S  ( 2G − 3K ) ∂S  (σ + σ~ + σ P ) kS  (σ + σ~ + σ P ) 1 + 2G σ ∂ ∂ σ 3   W& I = kT 1 −

in

W I (t ) ∂F (σ + σ~ + σ P ) ⋅ (σ + σ~ + σ P ) H (σ + σ~ + σ P ) ∂σ

in

R+ × Ω

(14)

R+ × Ω

(15)

σ (0, (r ,θ )) = 0 u (0, (r ,θ )) = 0 , ∀ (r ,θ ) ∈ Ω W (0, (r , θ )) = H (σ I

P

(16)

(θ )) .

In order to determine a numerical approach of the solution of the problem (14)-(16) we consider an interval [0, T ], T > 0 . Let us note

{

V1 = v = (v1 , v2 ,0) vi ∈ L2 (Ω ) , vi = vi ( r ,θ ) , vi (ma,θ ) = 0 , i = 1,2

[

V 2 = σ ∈ L2 (Ω ) 

]

3

}

(17)

σ = σ (r , θ ), Div σ = 0 , σ (a, θ )n = 0}

(18)

From (13) result that the solution (u , σ , W I ) of the problem (14)-(16) has the properties u ∈V 1 and σ ∈ V2 . t be the step time and Let M ∈ N, M > 2 , ∆t = M t 0 = 0 , t n+1 = t n + ∆ t , n = 0,K, M − 1 . Let us consider Vh ⊂ V1 a finite-dimensional subspace constructed using the finite element method. We determine

( )

 u n , σ n +1 , W I  h h 

n +1  h 



approach of the solution (u , σ , W I ) on the moment t n .

Let B = {ϕ1 ,K, ϕ I } ⊂ Vh be a base in Vh , dim Vh = I . Taking u h0 = 0 we determine u hn+1 ∈ Vh , n > 0 , such that:

u hn+1 =

I

∑α

n+1 j ϕj

,

j =1

where the constants α nj +1 , j = 1,K, I are the solution of a linear system. For the stress approach and irreversible stress work approach we consider σ 0j = 0 and

(W )

I 0 j

( )

= H (σ P ) and we determine σ jn +1 and W I

n +1 j ,

n ≥ 0 , using the

following implicite scheme:

σ hn +1 = σ hn + 2G[ε (u hn +! − ε (u hn )] + (3K − 2G) [ε (u hn +! ) − ε (u hn )] 1 + kT 1 −

W I (t ) H (σ hn+1 + σ~ + σ P )

∂F  ( 2G − 3K ) ∂F  (σ hn +1 + σ~ + σ P ) 1 + 2G (σ hn+1 + σ~ + σ P ) +  3 ∂ σ ∂ σ  

∂S  ( 2G − 3K ) ∂S  kS  (σ hn+1 + σ~ + σ P ) 1 + 2G (σ hn +1 + σ~ + σ P ) ∂ 3 ∂ σ σ   and, respectively W& I = kT 1 −

W I (t ) ∂F (σ hn +1 + σ~ + σ P ) ⋅ (σ hn +1 + σ~ + σ P ) . n +1 P ~ H (σ h + σ + σ ) ∂σ

,

(19)

(20)

For numerical solution we use the scheme proposed by Paraschiv-Munteanu ([5]) using finite elements methods for spatial integration and a complet implicit method for integration in time. For short period of time the creep solution and the numerical solution are very close. Thus deformation by creep and stress variation can simultaneously be described. The similar results for deep boreholes are obtained in papers [5] and [6].

498

REFERENCES [1] Cristescu N.: Constitutive equation of rock salt and mining applications, In: Seventh Int. Symp. on Salt, April 6-9, 1992, Kyoto, Japan, Elsevier Science, Amsterdam, vol. I, 1992, 105-115. [2] Cristescu N.: A general constitutive equation for transient and stationary creep for rock salt, Int. J. Rock Mech. Min. Sci., vol. 30, 1993, 125-140. [3] Cristescu N., Hunsche, U.: Time Effects in Rock Mechanics, Chichester-New York-Weinheim-BrisbaneSingapore-Toronto,1998. [4] Paraschiv I., Cristescu N.: Deformability response of rock salt around circular mining excavations, Rev. Roum. Sci. Techn. – Mec. Appl., tome 38, 3, 1998, 257-276. [5] Paraschiv-Munteanu I.: Metode numerice in geomecanica, Teza de doctorat, Universitatea din Bucuresti, 1997. [6] Paraschiv-Munteanu I., Cristescu N.D.: Stress relaxation during creep of rocks around deep boreholes, Int. J. of Eng. Sci., vol. 39, 7, 2001, 737-754. [7] Rosca I., Sofonea M.: Error estimates of an iterative method for a quasistatic elasto-viscoplastic problem, Applications of Mathematics, vol, 39, 6, 1994, 401-414. [8] Paraschiv-Munteanu I., Massier D.: Probleme de mecanica rezolvate in MATLAB, Editura Universitatii Bucuresti, ISBN 978-973-737-518-6, 2008, 206 pag. [9] ***MATLAB, User Guide, 2001.

499