An analytical solution for circular tunnels excavated in rock masses exhibiting viscous elastic-plastic behavior

An analytical solution for circular tunnels excavated in rock masses exhibiting viscous elastic-plastic behavior

International Journal of Rock Mechanics & Mining Sciences 124 (2019) 104128 Contents lists available at ScienceDirect International Journal of Rock ...

1MB Sizes 0 Downloads 252 Views

International Journal of Rock Mechanics & Mining Sciences 124 (2019) 104128

Contents lists available at ScienceDirect

International Journal of Rock Mechanics and Mining Sciences journal homepage:

An analytical solution for circular tunnels excavated in rock masses exhibiting viscous elastic-plastic behavior Ali Reza Kargar School of Mining Engineering, College of Engineering, University of Tehran, Iran



Keywords: Circular tunnels Viscous elastic-plastic model Time-dependent deformation

Investigating time-dependent behavior of tunnels driven in rocks has been the topic of interest for many years. The majority of researches were conducted with the assumption that the behavior of the surrounding rock masses is viscoelastic. In this paper, an analytical solution is proposed for circular tunnels under hydrostatic stress field, whose surrounding rock mass is isotropic, homogenous, obeying viscous elastic-plastic constitutive model. It is assumed that the rate of excavation is infinitely large. Therefore, the solution is proposed for unlined tunnels when the advance of tunnel face is halted, or after lining installation when the tunnel face is far from the studied section. The results were compared with those predicted by finite difference numerical model which exhibited a good agreement for both lined and unlined tunnels, except for large times where numerical model results slightly exceed the values predicted by the proposed solutions for lined tunnels.

1. Introduction Once a tunnel is excavated in a rock mass, two different factors could affect the excavation. First the advances of the tunnel face which is modeled by the distance between the current section and face of the tunnel. Another factor is the time-dependent behavior of surrounding rock mass which is occurred after excavation process. The advance of tunnel face is usually modeled by many researchers through applying elastic-plastic stress-strain analysis along with convergence-confinement concept. Among those, Duncan-Fama pre­ sented a solution for stress and deformation field around a circular hole in a Mohr Coulomb medium.1 Later, Brown et al. introduced a solution by considering Hoek-Brown yield function as the failure criterion for rock mass.2 Carranza-Torres and Fairhurst also proposed a closed-form solution by applying dimensionless method to their derivation, and managed to predict stress and deformation field, along with radius of plastic zone more accurately.3 Sharan developed a new solution using elastic-brittle plastic constitutive model mass.4,5 Finally, Alonso et al. took a milestone and presented a solution for ground response curve of rock masses with strain softening behavior.6 In addition to the above mentioned researches, many other investigations have been carried out considering time-dependent behavior of rock masses. Gnrik and Johnson determined deformation field analytically around a circular shaft with and without lining support, considering viscoelastic behavior for surrounding rock mass.7 Sakurai, by introducing an

equivalent initial stress, developed an analytical solution and obtained pressure acting on the tunnel supports vs. time and tunnel advance­ ment.8 Ladanyi and Grill investigated the effect of long term rock deformation on lining pressure through analytical methods.9 They employed both linear and non-linear viscoelastic Maxwell models into their calculations. Goodman used correspondence principle in deter­ mining deformation around a circular hole in order to apply back analysis of data from field tests to find viscoelastic constants of the rock masses.10 Pan and Dong proposed a solution for time-dependent tunnel convergence by employing tunnel advancement effect and support installation in a viscoelastic medium.11,12 Fahimifar et al. suggested a closed form solution for excavation of circular tunnels in viscoelastic burgers’ rock masses, in the case of both lined and unlined tunnels.13 Moreover, Nikomos et al. derived a solution for tunnels in viscoelastic materials by considering support installation time.14 As mentioned the majority of researchers considered rock masses as viscoelastic media, and rarely any researchers applied an analytical viscous elastic-plastic analysis to obtain tunnel responses to time-dependent effect due to excavation. Only Cristescu, and Sulem et al. proposed solutions to obtain deformation around a circular tunnel in rock masses with elastoplastic time-dependent behavior, and managed to predict the lining support reaction.15–17 However, their solutions suffered from a more compre­ hensive rheological model. In addition to above researches, there were other scientists who implemented numerical methods to investigate time-dependent effect on tunnel convergence.18–21

E-mail addresses: [email protected], [email protected] Received 28 February 2019; Received in revised form 4 October 2019; Accepted 12 October 2019 Available online 25 October 2019 1365-1609/© 2019 Elsevier Ltd. All rights reserved.

A.R. Kargar

International Journal of Rock Mechanics and Mining Sciences 124 (2019) 104128

Fig. 1. Rheological models: a) Kelvin, b) Maxwell, c) Burgers, and d) CVISC bodies.

support pressure on the tunnel wall.23 In this case, the support pressure could be introduced as follows:

The aim of this paper is to propose an analytical solution for stress and deformation field around lined and unlined circular tunnels driven in rock masses with Burger-creep viscous elastic-plastic behavior. Sub­ sequently, the pressure acting on lining support system is determined by considering time-dependent interaction between the lining support and surrounding rock mass, and tacking into account the unloading effect.

pi ¼ ð1


λe Þσ 0

where σ 0 is the initial hydrostatic stress, and λe , which is equal to the ratio of the current to ultimate radial displacements of the tunnel wall, is a parameter depending on the distance between the studied section and tunnel face. There are different formulations representing the parameter λe . However, in this study, due to elastoplastic behavior of the rock mass, the equation introduced by Panet and Guenot is used as follows24: 12 3 2 0

2. General consideration In order to investigate the problem of excavating a circular tunnel in a rock mass with viscous elastic-plastic behavior, several assumptions are made as follows:

6 λe ðxÞ ¼ 0:28 þ 0:7241

1. Homogenous, isotropic rock mass material following viscous elasticplastic constitutive model, 2. Viscos elastic-plastic model is a series combination of a Burgers’ body and a slider, known as CVISC constitutive model, 3. Mohr-Coulomb as the yield criterion, 4. Perfectly plastic behavior for slider body in CVISC constitutive model, 5. Extremely large rate of excavation, which let us to neglect the viscous behavior of rock mass during excavation process, 6. No elapsed time before or after support installation.

B @

0:84 C 7 A5 0:84 þ x=r



As seen in Eq. (2), λe varies from 0, related to a section in front of tunnel face, to 1, for a far section. In order to deal with this problem analytically, considering viscous elastic-plastic behavior for rock mass, viscoelastic and plastic behaviors, as well as the viscous elastic-plastic model used in this paper, known as CVISC constitutive model, will be briefly introduced in the following sections. 2.1. Viscoelastic behavior

The time-dependent behavior is investigated for unlined and lined tunnels only when the tunnel excavation is stopped, or after lining support is installed and tunnel excavation proceeds instantly. Due to extremely short periods of excavation process and support installation compared to long-term time studied for determining the deformations, no viscous behavior is considered within these periods. Owing to these assumptions, the time effect is employed to the solution either when the studied section is far from the tunnel face for lined tunnels (more than four times of the tunnel radius), or when the tunnel advance is stopped in case the tunnel has no support. Therefore, the instant elastic-plastic behavior only dominates the model at initial time due to the excava­ tion process. The solution for a circular hole in an infinite plate with linear elastic behavior was first introduced by Kirch.22 Assuming the plane strain condition with hydrostatic stress field the solution is reduced to that for a circular tunnel as an axisymmetric problem. Since the excavation of a tunnel is a three dimensional problem, the phenomenon could be investigated using a plain strain framework through applying a fictitious

Viscoelastic behavior is a combination of viscous, related to fluids, and elastic, associated with solids behaviors, which are represented by dashpot and spring, respectively, in rheological models. The two simplest viscoelastic models are Kelvin and Maxwell bodies which are shown in Fig. 1. The constitutive equations dominating these models are as follows for Kelvin and Maxwell bodies, respectively:

σ ¼ Gk γ þ ηk γ’






σ ¼ γ’ ηm

where Gk and ηk are Kelvin shear modulus and viscosity, Gm and ηm are shear modulus and viscosity of Maxwell model, and notation (’) denotes derivative with respect to time, as well. As shown in Fig. 1(c), the Bur­ gers’ body is a series combination of Maxwell and Kelvin’s models. Therefore, its constitutive relation could be expressed as the following equation: 2

A.R. Kargar

International Journal of Rock Mechanics and Mining Sciences 124 (2019) 104128

2.2. Plastic behavior Plastic behavior which is concerned with irreversible deformation of solids is presented by a slider in rheological models. Equations pre­ dicting the threshold stress level, which is the onset of plastic behavior, and the induced deformations are the primary concerns in the theory of plasticity. Generally, the yield function for a material could be expressed by the following equation26: � (12) f σ ij ; kij ¼ 0 where kij is hardening tensor, which controls the shape of yield function during plastic flow. The plastic strain rate is also obtained through applying gðσij ; kij Þ, called plastic potential function, according to the following equation26: � ∂g σ ij ; kij ε_ ij ¼ λ_ (13)


where λ_ is a scalar multiplier which could change during plastic flow. The notation (:) here means ∂τ∂ where τ could be fictitious time. If the functions f and g are the same it is called associated flow rule, otherwise it is named as non-associated flow rule in literature. Although there exist diverse plastic models -perfectly plastic, hardening and softening plastic behavior-the focus of this study is based on perfectly plastic behavior.

Fig. 2. A circular tunnel in an infinite medium.

ηk γ} þ Gk γ’ ¼

ηk Gm

τ} þ 1 þ

� Gk ηk ’ Gk þ τ þ τ G m ηm ηm


Defining the linear differential operator with respect to time, Dt � ∂∂t, Eq. (5) could be written as a continuation of the following general equation form10,25: n X

pi Dit γ ¼


m X

2.3. CVISC model CVISC model introduced by Itasca is a Burger-creep viscous elasticplastic model which is shown in Fig. 1(d). It consists of a series combi­ nation of a Burgers’ body and a slider, modelling viscoelastic and plastic behavior, respectively. Consequently, the total deformation could be expressed as follows:


qi Dit τ

i¼0 i

where Dit � ∂∂ti . For Burgers’ model n ¼ m ¼ 2. Eq. (6) could generally be solved by the help of linear differential operators as follows: Q γ¼ τ 2P

dεij ¼ dεpij þ dεve ij


By employing Eqs. (11) and (13), Eq. (14) could be rewritten as follows: � ∂g σ ij dsij dσ kk dεij ¼ dλ þ þ δij (15) ∂σij 2G 3K

where n X

2P ¼

pi Dit ; Q ¼


m X


qi Dit


where the first term in the right side of the equation is plastic strain emerged from slider body, and the next two terms are viscoelastic strains produced by Burgers’ model. In the following sections, the stress and deformation field around a circular tunnel assuming viscous elasticperfectly plastic behavior for surrounding rock mass is determined based on CVISC constitutive model, and subsequently, rock-support interaction is investigated with respect to time for the lining support installed.

Assuming τ as a constant variable as τ0 (as it is in creep problems), the deformation is obtained as follows: 0 1 γ¼

1 1 1 B ¼ @1 τ0 ; ​ ​ ​ GðtÞ GðtÞ Gk





1 t þ Gm ηm


In three dimensional problems, generally time-dependent de­ formations are assumed only due to deviatoric stress components, and hydrostatic stress component produces only elastic dilation.10,25 Therefore, by decomposition of stress and strain tensors, σ ij and εij , as follows:

εij ¼ δij


σij ¼ δij


σ kk 3

þ eij

3. Solution In this section, at first the solution for a circular tunnel in an elasticperfectly plastic media is presented. Subsequently, the solution is developed for a circular tunnel excavated in a rock mass with viscous elastic-perfectly plastic behavior. Finally, support-ground interaction is studied, and two different analytical methods are presented to deter­ mine stress and deformation fields in rock mass. Moreover, the pressure acting on the lining is predicted as well.


þ Sij

where eij and Sij are components of deviatoric strain and stress tensors, respectively, and δij is Kronecker delta, the constitutive stress-strain relation could be expressed as the following equations: Qsij ¼ 2Peij σii ¼ 3K εii In other words, the operator

in viscoelastic equations.

3.1. Elastic-perfectly plastic constitutive model

(11) P Q


Fig. 2 shows a circular tunnel excavated in a homogenous isotropic rock mass with elastic-perfectly plastic behavior, under initial hydro­ static stress p0 . The opening surface is subjected to an internal pressure pi . As pi decreases gradually, radial displacement increases and a plastic

plays the role of the shear modulus G


A.R. Kargar

International Journal of Rock Mechanics and Mining Sciences 124 (2019) 104128

zone develops around the tunnel. As mentioned above, the linear MohrCoulomb failure criterion dominates in the plastic region.

equation could be found for the radial displacement: u_ _ du_ þ β ¼ hðrÞ r dr


σ 1 ¼ σ c þ Kp σ 3 where σ c is the uniaxial compressive strength, and Kp ¼

1þsin ϕ 1 sin ϕ,


in which

ϕ denotes the friction angle of the rock mass. This yield criterion could be expressed with respect to radial and circumferential stresses within plastic zone as follows:

By integrating Eq. (26) from the onset of plastic deformation to the current fictitious time at each point in plastic zone, the fictitious time derivatives would be removed, and following equation would be resulted:

where σθ and σr are circumferential and radial stresses, respectively. Differential equation of equilibrium for this axisymmetric problem could also be expressed as the following equation: dσ r σ r σθ þ ¼0 dr r

du u þ β ¼ hðrÞ dr r

By substituting σ θ from Eq. (17) into Eq. (18) and solving the dif­ ferential equation, using the boundary condition of σ r ¼ pi at r ¼ ri , the radial and circumferential stresses are determined as follows: � �Kp 1 σc r σr ¼ þ D1 ri Kp 1 (19) � � �Kp 1 � Kp r σθ ¼ σc 1 þ Kp D 1 ri Kp 1 Kp 1.


hðrÞ ¼

σre þ D2

ue ¼



where D2 ¼ Kp



, and σre is the radial stress at elastoplastic interface

2p0 σc Kp þ 1


_ pr _ pθ

ε_ r ¼ ε þ ε ε_ θ ¼ ε þ ε


βσ r


p0 Þg


σ re Þ

H1 ¼

D 1 ri � βKp þ 1 Kp þ β

H2 ¼

� D2 ν ðβ þ 1Þ Kp þ 1 þ νð2p0 βþ1

νðβ þ 1Þ Kp þ 1


σc Þ

(32) � β D2 Kp σc þ D2 βþ1




p0 2G



As mentioned above, in order to predict time-dependent behavior of surrounding rock mass, it is assumed that tunnel is excavated instantly, and then its advance is stopped. To determine deformation field around the circular opening considering CVISC constitutive model for surrounding rock mass, cor­ respondence principle is used in this study. As stated in section 2.3, the induced deformation is divided into plastic and viscoelastic parts, whose viscoelastic part could be easily computed by applying the equivalent shear and bulk modulus. This concept could be incorporated into the derivation of equations by defining shear modulus and poisson’s ratio with respect to linear differential operators P and Q 10, 25:



ψ where β ¼ 1þsin 1 sin ψ , and ψ is the dilation angle; and applying Eq. (13), the

P 3K 2G G¼ ; ​ ​ ​ ​ ν ¼ Q 6K þ 2G

following equation could be found with respect to plastic strain rates:

ε_ pr þ β_εpθ ¼ 0

νÞðσ θ

3.2. Viscous elastic-perfectly plastic constitutive model (CVISC)

As mentioned earlier the notation (.) here means ∂τ∂ which τ is ficti­ tious time that could be the radius of plastic zone or the internal pressure on the opening surface. Note that this fictitious time is different from the long term time which is of interest in the next section. Assuming the plastic potential function as the following equation: g ¼ σθ

re ðp0 2G

uðrÞ ¼

where in axisymmetric problems u_ du_ εr ¼ ; ​ ​ ​ εθ ¼ dr r


p0 Þ þ ðβ

The stress and deformation field in elastic region are given according to Lame’s solution: �r �2 σr ¼ p0 ðp0 σ re Þ e r (34) �r �2 σθ ¼ p0 þ ðp0 σre Þ e r

It is obvious that as long as pi < σ re , a plastic zone develops around the opening, otherwise there exists only elastic region in the vicinity of the tunnel. In order to determine deformation field, it is assumed that non-associated flow rule dominates in this study. Generally, in incre­ mental theory of plasticity for infinitesimal deformations, the total radial and circumferential strain rates in plastic region could be decomposed into elastic and plastic parts as follows: _ er _ eθ

βνÞðσ r

where H1 and H2 are calculated as follows:


(r ¼ re ) which is obtained as follows:

σ re ¼


which results in the following equation for deformation field in plastic region � �Kp �� � � �Kp �r �β � 1 r re e 2Gue þ H1 (31) uðrÞ ¼ H1 þ H2 r þ H2 re 2G ri r ri

The radius of plastic region re could also be

�K 1 1

1 fð1 2G

where σr and σ θ are substituted from Eq. (19). Eq. (28) may be solved by using the following boundary condition at elastoplastic interface:

determined based on the continuity of the radial stress according to the following equation: re ¼ ri


Using Hooke’s law in the theory of elasticity for plane strain condi­ tion, hðrÞ could be found as follows:


where D1 ¼ pi þ


_ ¼ ε_ e þ β_εe hðrÞ r θ


σ θ ¼ σ c þ Kp σ r




Considering Burger’s body and the creep condition, the equivalent shear modulus GðtÞ could be expressed as Eq. (9), and poisson’s ratio νðtÞ is obtained based on Eq. (36). Substituting GðtÞ and νðtÞ into the derived

By applying Eqs. (22), (23) and (25) the following differential 4

A.R. Kargar

International Journal of Rock Mechanics and Mining Sciences 124 (2019) 104128

Fig. 3. A schematic representation of ground reaction curve for different rheological times.

Fig. 6. Tunnel wall displacement, predicted by Eq. (37) and finite difference numerical model (FLAC), vs. different rheological times when the tunnel sec­ tion is far (more than four times of the tunnel radius from the tunnel face).

uðr; t; pi Þ ¼

p0 σre �re ðpi Þ�2 r r 2GðtÞ


In case u ¼ uðri ; t; pi Þ, the time-dependent ground reaction curve (GRC) may be obtained with respect to internal pressure and tunnel wall displacement. As observed in Fig. 3, by increasing time, GRC is theo­ retically expected to shift to the right. In fact, in creep problems, it is assumed that when time elapses only deformation field changes, while the stress field remains the same. The topic of interest is when a support is installed, since the creep assumption is no more valid due to the fact that stress field changes during time because of the variation of internal support pressure owing to interaction between tunnel wall and support system. This problem which includes an elastic unloading is investigated in the next part by installing lining support inside the tunnel.

Fig. 4. Time-dependent ground reaction curve and lining support response curve (SRC).

3.3. Lining support installation Let lining support be installed at the distance x0 ¼ xðt0 Þ from the tunnel face. As assumed, the tunnel is excavated instantly and there is no elapsed time before or after support installation. Consequently, it is possible to let t0 � 0. The pressure acting on lining is subjected to the tunnel closure by the following equation: � � u uin ps ¼ K s (39) ri ri where uin ¼ uðx0 Þ is the tunnel closure before support installation, which is equal to λe ðx0 Þumax , and umax ¼ uðr ¼ ri ; t ¼ t0 ; pi ¼ 0Þ according to Eq. (37). Ks is also the lining support stiffness obtained as follows27: � Ec r2ex r2in � � Ks ¼ (40) ð1 þ νc Þ ð1 2νc Þr2ex þ r2in

Fig. 5. Ground reaction curve for different rheological times predicted by Eq. (37).

Ec and νc , in Eq. (40), are concrete Young’s modulus and Poisson’s ratio, rex and rin are also outer and inner lining radius, respectively (it is assumed that the lining behaves elastically). After installing the support, initial equilibrium of rock-support system is achieved when pi ¼ ps . In other words, u0 , the opening displacement subjected to this equilibrium position is obtained based on Eq. (37) as follows: � � �� u0 uin u0 ¼ u ri ; t0 ; Ks (41) ri ri

equations in section 3.1, and defining the relevant parameters as func­ tion of t and pi , the displacement around the opening is determined as follows in plastic region: � �Kp � �r ðp Þ�β � 1 r e i uðr; t; pi Þ ¼ H1 ðt; pi Þ þ H2 ðtÞr 2GðtÞ ri r �K � �� re ðpi Þ p 2GðtÞue ðt; pi Þ þ H1 ðt; pi Þ þ H2 ðtÞre ðpi Þ (37) ri

This equation could be solved numerically, and u0 could be obtained. As Fig. 4 depicts, from that point on, time-dependent interaction be­ tween lining and surrounding rock mass is started. As observed in Fig. 4, by each time increment Δt, there would be a displacement increment Δu

Furthermore, in elastic zone, radial displacement is as the following equation:


International Journal of Rock Mechanics and Mining Sciences 124 (2019) 104128

A.R. Kargar

Fig. 7. Radial displacement, predicted by Eqs. (37) and (38), and FLAC finite difference software, vs. scaled radial distance for different rheological times.

and corresponding support pressure increment Δps . Δu and Δps could be computed in an incremental method such as the following procedure: tnþ1 ¼ tn þ Δt Δui;nþ1 ¼ u ri ; tnþ1 ; ps;n Δps;nþ1 ¼ Ks

u ri ; tn ; ps;n

Δui;nþ1 ri


ui;nþ1 ¼ ui;n þ Δui;nþ1 ps;nþ1 ¼ ps;n þ Δps;nþ1 where n denotes the number of steps, and subscript i is subjected to the related parameter at the boundary of the tunnel. This trend may also be implemented in order to determine radial displacement and, radial and circumferential stresses at any point in the surrounding rock mass such as following method: Fig. 8. Tunnel wall displacement with respect to number of time steps in Eq. (42) for different rheological times.

tnþ1 ¼ tn þ Δt Δunþ1 ¼ u r; tnþ1 ; ps;n

Δui;nþ1 ¼ u ri ; tnþ1 ; ps;n Δps;nþ1 ¼ Ks

u r; tn ; ps;n �

u ri ; tn ; ps;n

Δui;nþ1 ri

unþ1 ¼ un þ Δunþ1 ps;nþ1 ¼ ps;n þ Δps;nþ1 8 � �Kp 1 σc r > > ps;nþ1 < σre ; ​ r < re þ D > 1 > > ri Kp 1 > > > > < �r �2 e σr ¼ ps;nþ1 < σ re ; ​ r � re > p0 ðp0 σre Þ r > > > > > > > ��ri �2 > : p0 p0 ps;nþ1 ps;nþ1 � σre r 8 � � � Kp 1 � Kp r > > σc 1 ps;nþ1 < σre ; ​ r < re þ Kp D1 > > > r K 1 p i > > > > < �r �2 e σθ ¼ p0 þ ðp0 σre Þ ps;nþ1 < σre ; ​ r � re > > r > > > > > > ��ri �2 > : p0 þ p0 ps;nþ1 ps;nþ1 � σre r

Fig. 9. Tunnel wall displacement, predicted by Eqs. (42) and (47), and FLAC finite difference code, vs. time.


Alternatively, this problem may also be solved by differentiating the 6

A.R. Kargar

International Journal of Rock Mechanics and Mining Sciences 124 (2019) 104128

Fig. 10. Radial displacement, predicted by a) Eq. (43), b) Eq. (48), and FLAC finite difference code, vs. scaled radial distance.

function u ¼ uðr; t; ps Þ as follows: du ¼

∂u ∂u ∂u *dps dr þ dt þ ∂r ∂t ∂ps

By considering dps ¼ Ks duri i , Eq. (46) could be written as follows: (44)

dui ¼ dt 1

where the notation (*) denotes convolution operator employed due to variation of pressure dps , which is defined as follows: Z tþdt ∂u ∂u * dps ¼ ðt þ dt τÞdps ðτÞ (45) ∂ ps ∂ps t

∂ui ∂ui *dp dt þ ∂t ∂ps s

jt ¼ 0


which is associated with the initial condition of ui ¼ u0 at t ¼ t0 . Regarding the fact that the right side of Eq. (47) are also function with respect to t and ps , where ps itself is equal to Kris ðui uin Þ, a nonlinear

equation is resulted whose solution may be the induced displacement along the opening wall considering time effects. Furthermore, radial displacement at any point in the rock mass could be found by applying Eqs. (44) and (45) as follows:

when r ¼ ri , the above equation is reduced to the following: dui ¼

∂ui ∂t Ks ∂ui ri ∂ps


du ∂u ∂u Ks dui ¼ þ jt ¼ 0 dt ∂t ∂ps ri dt

where the subscript (i) denotes the amount of studied variable for r ¼ ri . 7


A.R. Kargar

International Journal of Rock Mechanics and Mining Sciences 124 (2019) 104128

Fig. 11. Radial stress, predicted by a) Eq. (43), b) Eq. (48), and FLAC finite difference code, vs. scaled radial distance.

It is important to know that the calculation of relative derivatives of ui and u presented in Eqs. (47) and (48), respectively, are based on Eq. (37) or Eq. (38) depending on whether a plastic zone exists at the studied

point or not. Therefore, plastic region as follows:

� �K �r ðp Þ�β � re ðps Þ p e s 2GðtÞue ðt; ps Þ þ H1 ðt; ps Þ r ri �� �Kp � � 1 r ∂ ∂ re ðps Þ β ​ ​ ​ ​ ​ ​ ​ ​ ​ ​ þ H2 ðtÞre ðps Þ�g þ ½ 2G’ ðtÞue ðt; ps Þ H1 ðt; ps Þ þ r H2 ðtÞ ∂t ∂t 2GðtÞ ri r �K � �� ∂ re ðps Þ p ∂ ∂ ​ ​ ​ ​ ​ ​ ​ ​ ​ 2GðtÞ ue ðt; ps Þ þ H1 ðt; ps Þ þ re ðps Þ H2 ðtÞ ∂t ∂t ∂t ri �

∂u ∂t

is determined in viscous elastic-perfectly

� � Kp

∂u G’ ðtÞ r ¼ H1 ðt; ps Þ ri ∂t 2G2 ðtÞ

þ H2 ðtÞr



A.R. Kargar

International Journal of Rock Mechanics and Mining Sciences 124 (2019) 104128

in which: �∂ ∂ D 1 ri νðtÞ H ðt; ps Þ ¼ ðβ þ 1Þ Kp þ 1 ∂t 1 ∂t Kp þ β


�∂ � � ∂ νðtÞ H ðtÞ ¼ D2 Kp þ 1 þ 2p0 σ c ∂t 2 ∂t ∂ G’ ðtÞ ðp0 σre Þre ðps Þ ue ðt; ps Þ ¼ ∂t 2G2 ðtÞ

0 G ðtÞ ¼

∂ui ∂ps

in Eq. (47) are consequently determined by

4. Discussion

1 In this section the proposed solution for unlined circular tunnels is compared with finite difference numerical model (FLAC code) through an example. Thereafter, the results of the proposed stepwise procedure introduced by Eqs. (42) and (43), and numerical solution of Eq. (47) along with Eq. (48) for lined tunnels are compared with those obtained by FLAC finite difference code.28. The opening was supposed as a cir­

Be Gηkk t 1C C B þ CG2 ðtÞ B @ ηk ηm A


substituting r ¼ ri into Eqs. (49), (51), (53) and (54). Solving Eq. (48) in conjunction with Eq. (47), could give radial displacement at any point in the surrounding rock. Radial and circum­ ferential stresses may also be obtained similarly through applying Eq. (19) or (34) by replacing pi ¼ Kris ðui uin Þ.


∂ 18K νðtÞ ¼ G’ ðtÞ ∂t ð6K þ 2GðtÞÞ2

∂ui ∂t

In addition, ∂∂pus could be written as:

�� �Kp r ∂ β �re ðps Þ�β 1 ∂ H ðt; ps Þ r ðp Þ½ 2GðtÞue ðt; ps Þ ∂ps 1 ∂ps e s ri r r � � � �K re ðps Þ p re ðps Þ�β ∂ þ H2 ðtÞre ðps Þ ½ 2GðtÞ ue ðt; ps Þ ​ ​ ​ ​ ​ ​ ​ ​ ​ ​ ​ ​ ​ þ H1 ðt; ps Þ ∂ps ri r � � �� �K �K 1 re ðps Þ p ∂ Kp re ðps Þ p ∂ ∂ ​ ​ ​ ​ ​ ​ ​ ​ ​ ​ ​ ​ ​ ​ ​ ​ ​ ​ ​ ​ ​ þ H1 ðt; ps Þ ​ þ H1 ðt; ps Þ re ðps Þ þ H2 ðtÞ re ðps Þ ri ∂ps ∂ps ∂ps ri ri

∂u 1 ¼ ∂ps 2GðtÞ


νðtÞðβ þ 1Þ Kp þ 1

∂ p0 σ re ∂ ue ðt; ps Þ ¼ re ðps Þ ∂ps 2GðtÞ ∂ps


∂ ri 1 σre þ D2 r ðp Þ ¼ ∂ps e s Kp 1 D 1 D1 On the other hand,

region as:

∂u ¼ ∂t

G’ ðtÞ ðp0 2G2 ðtÞ

∂u ∂t

σ re Þr


cular tunnel with the radius of 3 m, under 7 MPa initial hydrostatic stress field. The data for viscoelastic properties of surrounding rock mass was adopted from Zhang et al. study on rock salt.28 Therefore, the co­ efficients for Burgers body were selected as Gk ¼ 6:73�104 ​ MPa, ηk ¼ 1:99�104 ​ MPa*day, Gm ¼ 1:99�104 ​ MPa, ηm ¼ 1:39�105 ​ MPa*day, and K ¼ 4:31�104 ​ MPa. Furthermore, the plastic properties of rock mass were assumed as c ¼ 0:8 ​ MPa, ϕ ¼ 23o and ψ ¼ 0. The concrete lining was supposed to have 30 ​ cm thickness, with elastic modulus and Poisson’s ratio of 1:65�104 ​ MPa and 0:2, respectively.


∂ ri � H ðt; ps Þ ¼ βKp þ 1 ∂ps 1 Kp þ β

∂u ∂ps

�r ðp Þ�2 e i r

∂u p0 σre �re ðps Þ� ∂ ¼ r ðp Þ ∂ ps e s r ∂ ps GðtÞ


�K 1 1 p

in Eq. (48) is calculated in viscoelastic

4.1. Comparison of the proposed solution for unlined circular tunnels with FLAC finite difference results


The ground reaction curves obtained according to Eqs. (37) and (38) for the aforementioned data set are depicted for different rheological times in Fig. 5. Assuming the studied section is a far section from the tunnel face, Fig. 6 shows tunnel wall displacement predicted from Eq. (37) and FLAC finite difference code. Moreover, Fig. 7 illustrates radial displacement vs. radial distance from the tunnel center, determined by Eqs. (37), (38), and FLAC. As observed a close agreement is evident between results from Eqs. (37), (38), and FLAC.


4.2. Comparison of the results predicted by proposed solution and numerical finite difference model (FLAC) for lined circular tunnels In this section, the lining is supposed to be installed 3 m from the tunnel face. Note that the solution of Eqs. (47) and (48) are employed into analysis through applying finite difference numerical method. The number of steps is chosen as 100 steps for Eqs. (42) and (43), since after which the difference in radial displacement for each time period is negligible according to Fig. 8. The radial displacement of the tunnel wall predicted by Eqs. (42), (47), and FLAC finite difference code is depicted vs. time in Fig. 9. Radial displacement and stress, according to Eqs. (43) and (48), are also illustrated with respect to scaled radial

Fig. 12. Thrust force, predicted by Eqs. (42) and (47), and FLAC finite differ­ ence code, vs. time. 9

A.R. Kargar

International Journal of Rock Mechanics and Mining Sciences 124 (2019) 104128

distance from the tunnel center in Figs. 10 and 11, respectively. It is observed from Fig. 9 that the radial displacement predicted by Eqs. (42) and (47) are close to each other in such a manner that they make an upper and a lower bound, and the curve resulted from FLAC software passes between these two bounds, except for large time values where this curve slightly exceeds. This trend could be occurred since the lining support reaction in the analytical methods and FLAC finite difference model is marginally different to tunnel wall displacement. In fact, lining reaction in the analytical model is included through lining stiffness, which is supposed here to be calculated according to Eq. (39), proposed by Brady and Brown based on Lame’s solution, while FLAC software computes lining reaction numerically by using liner elements. Since there is no momentum and shear force in lining due to hy­ drostatic stress field, the axial force induced in lining support may be determined as follows29: N¼

rin þ rex ps 2

References 1 Hoek E. Practical Rock Engineering. 2000. 2 Brown ET, Bray JW, Ladanyi B, Hoek E. Ground response curves for rock tunnels. J. Geotech. Eng. 1983;109(1):15–39. 3 Carranza-Torres C, Fairhurst C. The elasto-plastic response of underground excavations in rock masses that satisfy the Hoek–Brown failure criterion. Int J Rock Mech Min Sci. 1999;36(6):777–809. 4 Sharan S. Elastic–brittle–plastic analysis of circular openings in Hoek–Brown media. Int J Rock Mech Min Sci. 2003;40(6):817–824. 5 Sharan S. Exact and approximate solutions for displacements around circular openings in elastic–brittle–plastic Hoek–Brown rock. Int J Rock Mech Min Sci. 2005; 42(4):542–549. 6 Alonso E, Alejano L, Varas F, Fdez-Manin G, Carranza-Torres C. Ground response curves for rock masses exhibiting strain-softening behaviour. Int J Numer Anal Methods Geomech. 2003;27(13):1153–1185. 7 Gnirk PF, Johnson RE. The deformational behavior of a circular mine shaft situated in a viscoelastic medium under hydrostatic stress. In: The 6th US Symposium on Rock Mechanics (USRMS). American Rock Mechanics Association; 1964. 8 Sakurai S. Approximate time-dependent analysis of tunnel support structure considering progress of tunnel face. Int J Numer Anal Methods Geomech. 1978;2(2): 159–175. 9 Ladanyi B, Gill D. Design of tunnel linings in a creeping rock. Int J Min Geol Eng. 1988;6(2):113–126. 10 Goodman RE. Introduction to Rock Mechanics. New York: Wiley; 1989. 11 Pan Y-W, Dong J-J. Time-dependent tunnel convergence—II. Advance rate and tunnel-support interaction. International journal of rock mechanics and mining sciences & geomechanics abstracts. 1991;vol. 28:477–488. Elsevier. 12 Pan Y-W, Dong J-J. Time-dependent tunnel convergence—I. Formulation of the model. Int J Rock Mech Min Sci Geomech Abstr. 1991;28:469–475. Elsevier. 13 Fahimifar A, Tehrani FM, Hedayat A, Vakilzadeh A. Analytical solution for the excavation of circular tunnels in a visco-elastic Burger’s material under hydrostatic stress field. Tunn Undergr Space Technol. 2010;25(4):297–304. 14 Nomikos P, Rahmannejad R, Sofianos A. Supported axisymmetric tunnels within linear viscoelastic Burgers rocks. Rock Mech Rock Eng. 2011;44(5):553–564. 15 Sulem J, Panet M, Guenot A. An analytical solution for time-dependent displacements in a circular tunnel. Int J Rock Mech Min Sci Geomech Abstr. 1987;24: 155–164. Elsevier. 16 Cristescu N. Viscoplastic creep of rocks around a lined tunnel. Int J Plast. 1988;4(4): 393–412. 17 Critescu N. Viscoplastic creep of rocks around horizontal tunnels. Int. J. Rock Mech. 1985;22(6). 18 Ghaboussi J, Gioda G. On the time-dependent effects in advancing tunnels. Int J Numer Anal Methods Geomech. 1977;1(3):249–269. 19 Gioda G. A finite element solution of non-linear creep problems in rocks. Int J Rock Mech Min Sci Geomech Abstr. 1981;18:35–46. Elsevier. 20 Paraskevopoulou C, Diederichs M. Analysis of time-dependent deformation in tunnels using the Convergence-Confinement Method. Tunn Undergr Space Technol. 2018;71:62–80. 21 Peila D, Oreste P, Rabajoli G, Trabucco E. The pretunnel method, a new Italian technology for full-face tunnel excavation: a numerical approach to design. Tunn Undergr Space Technol. 1995;10(3):367–374. 22 Kirsch G. Theory of elasticity and application in strength of materials. Zeitschrift des Vereins Deutscher Ingenieure. 1898;42(29):797–807. � 23 PMaG P. Contribution au problEme de l’�etude du soutenement d’un tunnel derriere de front de taillevol. 1. Denver: ISRM.; 1974:1163–1168. 24 Penet M, Guenot A. Analysis of convergence behind the face of a tunnel. Int. Symposium Tunneling. 1982;82:197–204. 25 Mase GT, Mase GE. Continuum Mechanics for Engineers. CRC press; 1999. 26 Khan AS, Huang S. Continuum Theory of Plasticity. John Wiley & Sons; 1995. 27 Brady BH, Brown ET. Rock Mechanics: For Underground Mining. Springer Science & Business Media; 2013. 28 Zhang H, Wang Z, Zheng Y, Duan P, Ding S. Study on tri-axial creep experiment and constitutive relation of different rock salt. Saf Sci. 2012;50(4):801–805. 29 Carranza-Torres C, Labuz J. Class notes on underground excavations in rock. Topic. 2006;8:1–6.


Where rin and rex denote internal and external radii of the lining structure, respectively. Thrust forces are calculated based on Eqs. (42) and (47), and compared with FLAC software results in Fig. 12, which show a good agreement. 5. Conclusion An analytical solution is proposed for lined and unlined circular tunnels in rock masses with viscous elastic-perfectly plastic behavior which could simply predict stress and displacement fields in the sur­ rounding rock mass. Stress field variation resulted from rock-support interaction, and due unloading effect were taken into account in deri­ vation of formulas for lined tunnels. The results were compared with those of FLAC numerical software which coincided closely for unlined circular tunnels. Similarly, for lined circular tunnels, a close agreement was evident between results predicted by analytical solutions and FLAC finite difference software except for large times, where the values pre­ dicted by FLAC finite difference software slightly exceed those of analytical solutions. This trend could be due to the difference between methods implemented by analytical solutions and FLAC numerical software in order to calculate lining support stiffness against tunnel wall displacement. The former ones’ calculation uses the magnitude of lining stiffness proposed by Brady and Brown based on Lame’s solution, while the later one is based on numerical computation through applying liner elements. Therefore, since there exist diverse alternatives for computing lining support stiffness, it is expected that both analytical solution and FLAC numerical modelling converge into one for large times in case the magnitude of lining stiffness selected in analytical solution equals to that resulted from FLAC numerical modelling. Acknowledgement I would like to acknowledge Dr. Vahidoddin Fattahpour for his main role in language editing of this paper.