SciELO - Scientific Electronic Library Online

vol.16 número6Higher-order explicit time integration methods for numerical analyses of structural dynamics índice de autoresíndice de assuntospesquisa de artigos
Home Pagelista alfabética de periódicos  

Serviços Personalizados




Links relacionados


Latin American Journal of Solids and Structures

versão impressa ISSN 1679-7817versão On-line ISSN 1679-7825

Lat. Am. j. solids struct. vol.16 no.6 Rio de Janeiro  2019  Epub 29-Jul-2019 


Analytical solution of circular tunnel in elastic-viscoplastic rock mass

aState Key Laboratory for Geomechanics and Deep Underground Engineering, School of Mechanics & Civil Engineering, China University of Mining & Technology, Xuzhou 221116, PR China. Email:,,,,


This paper deals with a deep circular tunnel excavated in infinite homogeneous and isotropic elasto-viscoplastic rock mass subjected to a hydrostatic initial stresses. The tunnel is divided into the initial plastic, viscoplastic, and elastic zones. By combining the generalized Bingham model with Mohr-Coulomb yield criteria, analytical solutions of the circular tunnel are derived with the consideration of non-associated flow rule and elasto-viscoplasticity. The initial plastic zone is defined as the instantaneous change of rock mass excavation. Based on the initial plastic zone, the stresses in viscoplastic zone are transferred to the deep surrounding rock with time due to the initial earth stresses. The results are compared with the elastic-brittle solution at static conditions, and the solutions of this paper are validated. Moreover, the presented results shows that the stress and displacement of the surrounding rock varies with time, and the solutions can be obtained at different instants of time.

Keywords: Circular tunnel; Surrounding rock; Elasto-viscoplasticity; Analytic solution

1. Introduction

Many coal mines and metal mines have reached a mining depth of more than one kilometer. Owing to the increase in mining depth, dangers to the safety of the mines significantly increases, even with minor negligence [1]. Increased initial earth stress causes the surrounding rock to be in a state of rupture. The expansion of the surrounding rock rupture range is temporally dependent. This leads to a nonlinear continuous deformation and a difficulty in maintaining the stability of surrounding rock. For soft rock mass, the creep displacement is significant with a deformation ratio of tens of millimeters per day. If the supporting intensity is weak, the deformation ratio can reach tens of centimeters per day [2]. The stress and displacement distribution of the rock surrounding the roadway rock influence the stability and safety of roadway support system. Therefore, it is of great significance to study the rheological properties of surrounding rock in deep roadway [3]. Typical rheological failures of the roadway are shown in Fig. 1.

Fig. 1 Rheological failures of laneway [4

The axisymmetric deformation analysis of the rock around a circular tunnel is a classic problem. The ground responses of elasto-perfectly-plastic, elasto-brittle-plastic and strain-softening have been studied. Park and Kim [5] derived analytical solutions of the stress and deformation of the surrounding rock using an elasto-brittle model based on different elastic deformation laws in the plastic zone. Sharan [6-7] derived the closed-form expressions of stress and deformation in the elasto-brittle-plastic Hoek-Brown rock mass for circular tunnels. Jiang et al. [8] studied the stress and deformation solution of circular tunnels in the elasto-plastic-brittle rock mass on the account that the deep underground rock mass possesses a higher plastic deformation capacity. Stille et al. [9] combined the Mohr-Coulomb (M-C) criterion with the elasto-brittle-plastic model, and considered the dilatancy characteristics of the surrounding rock after yielding. Brown et al. [10] proposed the distribution law of stress and strain fields of surrounding rock, and considered the material characteristics of plastic softening and volume expansion.

Strictly speaking, the deformation and stress responses of the roadway obtained by elastic and elastoplastic mechanics were instantaneous in the past. In fact, the deformation and stress which were measured and maintained by the engineering support were caused by rheology. Therefore, solving the practical problems of rock mechanics cannot be separated from the analysis of rheology. In consideration of the elasto-viscoplastic constitutive theory, Perzyna proposed an elasto-viscoplastic model in 1966 [11]. Subsequently, Perzyna’s general theory of elasto-viscoplasticity was applied to metal alloys by Lemaitre and Chaboche [12], and was extended by Pellet et al. [13] to the study of the time-dependent behavior of rocks. Boidy et al. [14] and Bonini et al. [15] combined numerical simulations with Lemaitre’s viscoplastic model. In turn, Pellet [16] combined the Burger and the elastoplastic models, and calibrated the constitutive parameters using back-analyses of convergence data.

Barla et al. [17], Debernardi et al. [18], and D. Sterpi et al. [19], have conducted corresponding studies of deep tunnels using numerical simulations at squeezing conditions. Roatesi [20] obtained the numerical solution using the elasto-viscoplastic Cristescu constitutive law, which was used to perform an analysis of the evolution of the damage zones around a circular lined cavity in a viscoplastic material subjected to a nonisotropic far stress field. Many scholars have conducted corresponding research on the analytical solution of surrounding rock regarding rheology. Specifically, Sulem [21], Ladanyi [22], Gill [23], and Wang [24] et al., analyzed the stress and strain of the surrounding rock of the roadway based on the viscoelastic model, but neglected the plastic mechanical properties. This led to the overestimation of the support pressure in the support design.

In recent times, research related to the viscoplastic analysis of the surrounding rock was less extensive. Fritz [25] used the viscoplastic constitutive relation of the surrounding rock in a manner similar to the form of the linear elastic Hooke's law, and the elastic-viscoplastic analysis was carried out by ensuring that the stress and strain satisfy the basic equations of the problem and the stress intensity condition. This study was based on the generalized Bingham model in accordance to the correspondence principle that has been extensively adopted in plastic and viscoelastic mechanics. Additionally, its three-dimensional elastic-viscoplastic constitutive equation is established by rigorous mechanical analyses, and the stress and deformation solutions of elastic-viscoplastic surrounding rock are obtained in circular roadways. Combined with engineering examples, the variation law of stress and displacement of the surrounding rock was analyzed at different times given its significance for roadway engineering design and construction.

2. Problem definition

2.1 Mechanical model

The purpose of this paper is to obtain the distributions of stress and displacement of surrounding rock around a circular tunnel in elastic-viscoplastic rock mass. The characteristics of brittle rock mass are that brittle failure occurs and steps into the post-failure region once the stresses exceed the elastic limit.

When the tunnel is located in brittle rock mass, the deformation can be divided into two parts after stress redistribution, namely, the elastic region and the fractured region as shown in Fig. 2. The emerging conditions of the two parts can be described in accordance to two conditions: a) when the stress of the surrounding rock does not reach the plastic condition, the deformation is considered to be in the elastic state, and b) when the stress reaches and exceed the plastic yield condition (that is, when no deformation occurs and when the stress decreases at the same time), the deformation is considered to be in the fractured region. For convenience of analysis, the stress-strain curve of the actual rock is replaced by the fold line in Fig. 2.

Fig. 2 Stress-strain curves for elastic-brittle rocks [26

The plane strain problem requires the deformation of long roadway to be completed instantaneously. In fact, the deformation of the surrounding rock is gradually carried out. In this paper, the surrounding rock is assumed in an elastic state at the beginning. Stress concentration of surrounding rock leads to plastic yielding and fracture of rock and causes viscoplastic deformation with the excavation. Finally, when the stress distribution of the surrounding rock reaches equilibrium again, the viscoplastic deformation no longer changes with time.

Fig.3 shows a circular tunnel with radius a in a homogeneous infinite isotropic rock mass subjected to a hydrostatic stress p 0 in the cross section of the tunnel. A support pressure p i is applied on the inner surface. The surrounding rock consists of three parts: the initial plastic region, the viscoplastic region and the elastic region. The initial plastic zone is generated due to excavation of rock mass, and it is transiently formed. The viscoplastic zone continues to develop with time on the basis of the above. There is no stress disturbance beyond the viscoplastic zone, it is in an elastic region.

Fig. 3 Elasto-viscoplastic analysis model of a circular opening 

Assuming the rock mass is uniform, isotropic and the tunnel is a problem of plane strain, the equations are expressed in polar coordinates, the radial stress, radial strain, tangential stress and tangential strain are expressed as σr , εr , σθ and εθ respectively. The displacement is represented by u, variable superscript “c” and “e” represents fractured zone and elastic zone, respectively.

2.2 Governing equations

The differential equation of equilibrium for axisymmetric problems in each region can be expressed as (not considering the body force of rock mass)

σrr=σθσrr (1)

The geometric equation of this axisymmetric problem can be denoted as

{εr=urεθ=ur (2)

In this study, the M-C yield criterion was considered when the plastic yield and rupture of the surrounding rock occurred. The strength theory is suitable for homogeneous isotropic rock, and the calculated results were also generally consistent with the actual situation. The M-C criterion cannot only describe the failure characteristics of brittle materials, but also describes the damage characteristics of the friction materials. It can be uniformly represented as:

F=σθNσrS=0 (3)

where N and S are material parameters. For the intact and failure rock mass, N and S can be expressed by the cohesion strength and the friction angle as follows:

{Np=(1+sinφp)/(1sinφp)Sp=2cpcosφp/(1sinφp) (3a)

{Nc=(1+sinφc)/(1sinφc)Sc=2cccosφc/(1sinφc) (3b)

where c is the cohesive strength and φ is the frictional angle. The subscript “c” and “p” denote the fractured region and the plastic region, respectively.

The complexity of the elasto-plastic problem lies in the nonlinear relationship between the stress and the strain. According to the theory of plasticity, plastic strain depends on the plastic potential, corresponding M-C criteria, the plastic potential [27, 28] here is:

Φ=σθασr (4)

where α=(1+sinψc)/(1sinψc) is material parameter, ψc is the dilation angle of the fractured region.

The plastic strain can be obtained from the potential theory, as shown in the following equation:

εrp=λΦσr,εθp=λΦσθ (5)

On the basis of the plastic potential (4), combined the above results with the elastic strain, the total strain of rock mass induced by excavation can be deduced as

εr=[(1u)(σrp0)μ(σθp0)]/(2G)λαεθ=[(1u)(σθp0)μ(σrp0)]/(2G)+λ} (6)

Where G is the shear modulus, μ is poisson ratio, λ is plastic strain. When λ=0 , equation (6) is the constitutive equation of the elastic region.

It is assumed that after the stress of the surrounding rock exceeds the yield strength σs , viscoplastic deformation begins to appear. And this constitutive model, which describes this kind of phenomenon, is characterized by generalized Bingham mechanical model with a combination of a Hookean (spring), Newtonian (dashpot) and St.Venant (slider) elements as shown in Fig. 4.

Fig. 4 Generalized Bingham mechanical model 

One dimensional constitutive equation of generalized Bingham model is succinct, as shown in the following equation:

ε˙=σ˙E+σσsη (7)

where dot “ ∙ “ denotes the derivative of time, E is elastic modulus, σs is yield strength of plastic elements, η is viscosity coefficient of Newtonian element.

By the correspondence principle, the three-dimensional constitutive equation of the generalized Bingham model can be denoted as:

e˙ij=s˙ij2G+sij2ηε˙kk=σ˙kk3K+σkk2η'} (8)

where sij is deviatoric stress tensor, eij is deviatoric strain tensor, σij is stress tensor, εij is strain tensor, both stress and strain tensor can be decomposed into deviatoric stress and hydrostatic pressure, namely, sij=σijσkkδij/3 , eij=εijεkkδij/3 , σkk and εkk indicate the first invariants of stress and strain. δij is Kroneker symbol, when i=j , δij=1 , when ij , δij=0 . G is shear modulus, and K is bulk modulus, η is viscosity coefficient of deviatoric stress and deviatoric strain, η' is viscosity coefficient of hydrostatic pressure ( σkk ) and volume strain ( εkk ).

It is generally considered that the volume deformation of the rock mass is elastic. Here, in order to describe the general situation, it is considered that the volume deformation of the surrounding rock and the deviatoric strain tensor both present the elastic-viscoplasticity.

For circular tunnel, the Eq. (8) is unfolded in a cylindrical coordinate system:

ε˙r=12G[σ˙r(132G9K)Θ˙]+12η[σr(13η3η)Θ˙]ε˙θ=12G[σ˙θ(132G9K)Θ˙]+12η[σθ(13η3η)Θ˙]ε˙z=12G[σ˙z(132G9K)Θ˙]+12η[σz13(1ηη)Θ˙]} (9)

where Θ=σr+σθ+σz is the volume stress.

By substituting plane strain condition ( εz=ε˙z=0 ) into the third equation of Eq. (9), axial stress can be written as:

σz=η'η2η'+η(σr+σθ)σ˙z=3K2G6K+2G(σ˙r+σ˙θ)} (10)

Under the condition of plane strain, axial stress is zero, therefore, the relationship between the viscosity coefficient η and η' in the Eq. (10) can be obtained:

η'=3K2Gη (11)

By substituting Eq. (11) into Eq. (9), adopting the form of Eq. (6), it is considered that the elastic-viscoplastic deformation is produced after the elastic deformation occurs, (3K2G)/(6K+2G) is replaced by the Poisson ratio μ , the three-dimensional constitutive equation of generalized Bingham body under the condition of plane strain is obtained:

ε˙r=12G[(1μ)σ˜˙rμσ˜˙θ]+12η[(1μ)σ˜rμσ˜θ]ε˙θ=12G[(1μ)σ˜˙θμσ˜˙r]+12η[(1μ)σ˜θμσ˜r]} (12a)

σ˜r=σrσrcσ˜θ=σθσθc} (12b)

where σ˜r and σ˜θ indicates that viscoplastic flow happens after the stress exceeds the plastic yield limit.

2.3 Boundary conditions and initial conditions

The stress boundary condition, initial condition and contact conditions can be expressed as follows:

r=a,σr=pir=Rc,σr(Rc,t)=σr(Rc+),u(Rc,t)=u(Rc+)r,σr=p0t=0,σr=σr0e,σθ=σθ0e,u=u0e} (13)

where Rc and Rc+ represents the left and right boundary of r= Rc , respectively.

3 Deformation analysis of broken rock in the final state

In the final state of the surrounding rock, the fractured region and the elastic region are formed, and the stress components and the displacement can be expressed as follows [8, 26].

3.1 Cracked zone (arR c )

Eqs. (1), (3), (3b) and (13) will be used to get the stress of the cracked zone:

σrc=pic(r/a)Nc1Dcσθc=Ncpic(r/a)Nc1Dc} (14)

where Dc=Cccotϕc , pic=pi+Dc .

By using the geometric Eq. (2), the constitutive Eq. (6) and the deformation coordination equation, the corresponding expressions of strain and displacement can be obtained:

εrc=[NcB1(r/a)Nc1B2+Aαc(Rc/r)αc+1]/(2G)εθc=[B1(r/a)Nc1B2A(Rc/r)αc+1]/(2G)uc=rεθc} (15)

where B1=[(1μ)(Ncαc+1)/(Nc+αc)μ](pi+Cccotφc) , B2=(12μ)(p0+Cccotφc) , A is an integral constant.

3.2 Elastic zone (Rc ≤ r <∞)

The stress, strain and displacement of the elastic zone are easily obtained:

σre=p0+(pcp0)Rc2r2σθe=p0(pcp0)Rc2r2} (16)

εre=(pcp0)Rc2r2/(2G)εθe=(p0pc)Rc2r2/(2G)ue=rεθe} (17)

where pc is radial stress of the interface between the elastic zone and the cracked zone, Rc is the radius of the interface, they can be represented as [29]:

pc=(pi+Cccotφc)(Rc/a)Nc1Cccotφc (18)

Rc=a{p0(1sinφp)Cpcosφp+Cccotφcpi+Cccotφc}1/(Nc1) (19)

Using the condition of continuous displacement at the junction of the cracked zone and the elastic zone, i.e. the second equation of Eq. (13), we can get the solution for A:

A=(1μ)(Nc+1)(αc+1)Nc+αc(pi+Cccotφc)(Rca)Nc12(1μ)(p0+Cccotφc) (20)

By substituting Eq. (20) into Eq. (17), the expression of the strain and displacement in the cracked zone can be obtained. If the rock has not been broken, the above results can be returned to the famous Kastner solution, therefore, the aforementioned analysis is a supplement of Kastner solution.

4. Elastic-viscoplastic solution

Elastic deformation instantaneously occurs at the moment of tunnel excavation, and the surrounding rock near the boundary of the roadway may generate a viscoplastic area with radius Rt due to the release of stress, and beyond this radius is the elastic zone.

4.1 Transient elastic stress and deformation of surrounding rock

When the surrounding rock is in the elastic state, its stress component and displacement can be expressed as:

σr0e=p0+(pip0)a2r2σθ0e=p0(pip0)a2r2u0e=12G(p0pi)a2r} (21)

where variable superscript “e” represents elastic state, the subscript “0” represents the initial moment (t=0).

Through analyzing rheological test data [30], elastoplastic deformation occurs, but the magnitude of the elastic deformation is far greater than the plastic deformation, the plastic deformation is ignored here, only elastic deformation exists. Hence, at the moment of t=0+ , this problem can be described by Eq. (21).

4.2 Initial plastic zone stress and relative displacement of surrounding rock

Due to the excavation of rock mass, the stress of the surrounding rock is redistributed, we assume that an initial plastic zone produces at the moment of excavation (i.e. t=0+ ).

When the stress of the surrounding rock reaches the plastic yield condition, it enters the critical plastic state. By substituting the expression of elastic stress Eq. (21) into yield condition Eq. (3) and Eq. (3b), initial critical plastic zone radius is:

ρ0=a(1+Np)(pip0)(1Np)p0Sp (22)

Combined with the basic equations (1), (2) and (11a), we can get an equation:

rr[σ˜˙r+σ˜˙θ+β(σ˜r+σ˜θ)]=0 (23)

where β=G/η .

Eq. (23) is partial differential equation, its solution can be expressed as:

σ˜r+σ˜θ=f(t)+eβtg(r) (24)

where, f(t) and g(r) is functions of t and r respectively. It is an undetermined function and for simplicity its integral and other equations can be used in the original form, without affecting the result of the problem.

By using the first and fourth equation of the boundary and initial conditions Eq. (13), and the equilibrium Eq. (1), the radial stress and the circumferential stress can be expressed as:

σr=σrc(1eβt)+σr0eeβt+(1a2r2)f(t)σθ=σθc(1eβt)+σθ0eeβt+(1+a2r2)f(t)} (25)

By substituting Eq. (25) into the second equation of Eq. (12a), and integrate Eq. (12a) on t, the displacement u can be obtained:

u=r/(2G){[(1μ)σθ0eμσr0e]eβt+[(1μ)σθcμσrc](1eβt)+(12μ+a2r2)f(t)}+r/(2η){[(1μ)(σθ0eσθc)μ(σr0eσrc)](1eβt)/β}+r/(2η)(12μ+a2r2)f(t)dt+g1(r) (26)

The displacement contact condition is replaced by Eq. (26), that is, the fourth equation of Eq. (13):

g1(r)=r(2μ1)p0/(2G) (27)

Combined Eq. (27) with Eq. (26), the displacement of the initial plastic zone is shown as:

u=r/(2G){[(1μ)σθ0eμσr0e]+(12μ+a2r2)f(t)}+r/(2η)(12μ+a2r2)f(t)dt+r(2μ1)p0/(2G) (28)

The above equation, σr0e and σθ0e are both equivalent to the values under r=a conditions.

4.3 Stress and relative displacement of viscoplastic zone of surrounding rock

Suppose total stress exceeds yield strength, outside the scope of the initial plastic zone (i.e. ρ0rRt ). This hypothesis affects the evolution of stress and strain with time. When the cracked zone extends to the radius Rt at the moment of t, the stress and displacement are described as similar to the elastoplastic states. Hence, the stress in viscoplastic zone is similar to the form of Eq. (14), the displacement is similar to the form of Eq. (15), in these equations, a and pi is replaced by ρ0 and P(t) respectively. The forms of stress, strain and displacement after transformation are shown as:

σrvp=(P(t)+Cccotφc)(r/ρ0)Nc1Cccotφcσθvp=Nc(P(t)+Cccotφc)(r/ρ0)Nc1Cccotφc} (29)

uvp=r[B1vp(r/ρ0)Nc1B2Avp(Rt/r)αc+1]/(2G) (30)

where, B1vp , Avp and Rt , as shown in Eq. (31). In these equations, the radius and the stress are also replaced accordingly.

B1vp=[(1μ)(Ncαc+1)/(Nc+αc)μ](P(t)+Cccotφc)Avp=(1μ)(Nc+1)(αc+1)Nc+αc(P(t)+Cccotφc)(Rta)Nc12(1μ)(p0+Cccotφc)Rt=ρ0{p0(1sinφp)Cpcosφp+CccotφcP(t)+Cccotφc}1/(Nc1)} (31)

In order to separate the stress, strain, and displacement of the cracked zone from the corresponding variables of the viscoplastic zone, these variables are labeled with the superscript “vp”. In the expression of Rt , the radial stress is expressed in P(t) , and is defined in the form of the Eq. (25), as shown in Eq. (32).

P(t)=σrc(1eβt)+σr0eeβt+(1a2ρ02)f(t) (32)

In essence, the deformation of the surrounding rock refers to the process of constant stress adjustment. During the process of surrounding rock deformation, the stress in the viscoplastic zone changes with time. In addition, when time t , the stress in the viscoplastic zone is consistent with the stress obtained by the elastoplastic solution at the position of the corresponding radius. Certainly, the initial plastic zone will not increase indefinitely, and the development of viscoplasticity will eventually become plastic. In the case of partial brittle rock mass, when the plastic deformation occurs, brittle failure occurs, and enters into the post peak plastic behaviors. In the case of the rock mass material, the post peak residual strength cannot be ignored, and the radius of the final state plastic zone is the radius of the cracked zone.

The functions f(t) in Eq. (28) can be determined by the equal displacement at the position of r=ρ0 , i.e. Eq. (28) and Eq. (30) at the boundary r=ρ0 are equal, an integral equation about f(t) is obtained:

[(1μ)(Ncαc+1)/(Nc+αc)μ](P(t)+Cccotφc)(1μ)[(Nc+1)(αc+1)Nc+αc(P(t)+Cccotφc)(Rtρ0)Nc+αc2(p0+Cccotφc)(Rtρ0)αc+1](12μ)Cccotφc+μσr0e(1μ)σθ0e(12μ+a2ρ02)f(t)(12μ+a2ρ02)βf(t)dt=0 (33)

With the continuous development of time, the viscosity of the surrounding rock gradually disappeared, finally, the viscoplastic zone no longer changes, the boundary between the viscoplastic zone and the elastic zone becomes the junction of the cracked zone and the elastic zone, the Bingham elasto-viscoplastic model degenerates into an elastoplastic model.

4.4 Calculation procedure

By inputting mechanical parameters of surrounding rock, f(t) and f(t) can be calculated by MATLAB software. According to the derived analytical solutions, the radius of the viscoplastic region at different time can be obtained. When the radius of the viscoplastic region no longer changes, it indicates that the surrounding rock has reached equilibrium again. A flow chart which summarizes the implementation of the proposed method is shown in Fig. 5.

Fig. 5 Flowchart of the implementation 

5. Example analysis

In order to illustrate the correctness of the derived equation, an engineering example is discussed. Tang Kou coal mine was excavated in crosscut, which is located in Jining, Shandong, China. After the construction of chambers and tunnels, large deformation and failure occurred. The stress measurement of the original rock was carried out at the depth of 1028 m in tunnel, the vertical and horizontal stress components of the 2# measuring point is 23.77 and 23.04 MPa, respectively. The stress of the original rock is basically isobaric in two directions, herein, the original crustal stress is taken as 23.4 MPa. The physical and mechanical indexes of rock mass as shown in Table 1.

Table 1 Physical and mechanical parameters of rock mass 

Rock mass parameters Value
Radius, a (m) 1
Initial ground stress p 0 (MPa) 23.4
Modulus of elasticity E (MPa) 5700
Viscosity coefficient η ( Days·MPa ) 6000
Poisson ratio μ (-) 0.3
Internal friction angle φ p (Deg.) 35.82
Residual internal friction angle φ c (Deg.) 22.99
Cohesive force c p (MPa) 2.51
Residual cohesive force c c (MPa) 1.43
Dilatancy angle ψ c (deg.) 13

5.1 The solution of function f(t)

Eq. (33) is a transcendental equation of the function f(t) , an analytical expression of f(t) can not be obtained directly. By imposing the proper boundary conditions and using MATLAB, the approximate solutions of f(t) and f(t) can be obtained. The accuracy of the approximate solutions is related to the selected time step, for the requirement of setting precision, the stickiness of the specific problem should be taken into consideration, and the reasonable time step should be selected.

By substituting the rock mass parameters into the Eq. (33), analytical solutions of f(t) and f(t) can be obtained, as shown in Figs. 6 and 7. Fig. 6 denotes the curve of f(t) with time, with (Et)/η as abscissa and f(t)/p0 as ordinate; Fig. 7 shows the curve of f(t) with time, with (Et)/η as abscissa and (Ef(t))/(ηp0) as ordinate.

Fig. 6 Analytic solution of function f(t) under different support force 

Fig. 7 Analytic solution of function f(t) under different support force 

It is obvious that f(t) and f(t) vary with different support force. When the support force becomes larger, the initial value of f(t) becomes larger, and all of them increase first and then decrease, but eventually it tends to 0. The time required for f(t) to achieve balance is shortened when the support force is increased.

5.2 Distributions of stresses

When the support force is 0 MPa, the radius of the initial plastic zone is 1.2197 m and the final plastic zone is 2.5236 m. There is a viscoplastic zone between the above two zones. The calculation time is taken as 0, 10, 50 and 500 days, respectively, the corresponding stress states are States I, II, III and IV. The distributions of radial and tangential stresses of the above four states are shown in Fig. 8. The state I is fully elastic, plastic deformation doesn’t occur. The radius of the state II is 1.7381 m, the radius of the state III is 2.3303 m, and the radius of the state IV is 2.5236 m. Thus it can be seen, the development of the viscoplastic zone is a process of increasing radius. When the radius no longer changes, it means that the viscoplastic zone no longer develops, and the surrounding rock reaches a state of balance. The difference of the radial stress and circumferential stress at the radius of the initial plastic zone decreases with time.

Fig. 8 Stress variation of surrounding rock under different times 

5.3 Comparison of elasto-viscoplasticity and elastic-brittle displacement solution

This paper mainly studies the variation of stress and displacement of surrounding rock with time. Here the cases of states II, III and IV are discussed especially. Fig. 9 shows the distributions of viscoplastic displacement corresponding to States II, III and IV, respectively. Obviously, the displacement of above four states has similar distribution forms. With the increase in radius, the displacement gradually decreases. The displacement is continuous at the interface of the initial plastic zone and the viscoplastic zone. The maximum displacement occurs on the free face. And the displacements of the three states are 0.019 m, 0.04 m and 0.049 m separately in the free face of the tunnel. The three points (i.e. A, B and C, shown in Fig. 9) represent the Rt values of the latter three states respectively. The radius of the viscoplastic zone increases obviously. The elasto-viscoplastic displacement solution at the State III is in good agreement with that of the elastic-brittleness [26].

Fig.9 Elasto-viscoplastic and elastic-brittle displacement solutions 

The three-dimensional elasto-viscoplastic model established in this study is suitable for homogeneous isotropic rock. The results based on the elasto-viscoplastic model are in agreement with the static analytical solution of the elasto-brittle state. The proposed analytical solutions consider the progressive evolution process in plastic region, which is an improvement for the static analytical solution.

6. Conclusions

In this study, the elasto-viscoplastic analytical solutions of stress and deformation around a circular tunnel were derived according to the generalized Bingham model using the Mohr-Coulomb criterion and the non-associated flow rule.

The rock is assumed to be homogeneous and isotropic. When excavating tunnel, rock mass behaves elasticity. Then, when t=0+, the radius of the initial plastic zone reaches ρ0 . Up to the actual radius Rt of the plastic zone (i.e. for ρ0rRt ), expressions (31) are relevant and for rRt equations (16) and (17) are relevant, whereby pc , Rc and the index c have to be replaced by P(t) , Rt and t respectively. The unknown function f(t) (cf. Fig. 6 and 7) can be evaluated for an arbitrary time form (33). With the development of time, the radius of the viscoplastic zone increases, and the elasto-viscoplasticity will eventually develop into an elastic brittleness state, the two final radii is consistent.

The emergence of viscoplastic zone causes redistribution of stress in the surrounding rock, that is, the peak stress is transferred to the deep surrounding rock, and rheological analyses can thus better describe the changing process of developed stress and strain in the surrounding rock of deep tunnel. This method can effectively predict the failure process of the roadway and provides a reference for engineering support.


This study was supported by the Fundamental Research Funds for the Central Universities of China (No. 2015XKMS015) and the support is gratefully acknowledged.


[1] He M. C., Xie H. P., Peng S. P., Jiang Y. D. (2005). Study on rock mechanics in deep mining engineering. Chinese Journal of Rock Mechanics and Engineering 24: 2803-2813 (in Chinese). [ Links ]

[2] Zhang Q., Li C., Min M., Jiang B. S., Yu L. Y. (2017). Elastoplastic Analysis of Circular Openings in Elasto-Brittle-Plastic Rock Mass Based on Logarithmic Strain. Mathematical Problems in Engineering 2017: 1-9. [ Links ]

[3] Sun J. (1999). Rheology and engineering application of geotechnical materials. Beijing: China Construction Industry Press. [ Links ]

[4] Yang S. Q., Chen M., Jing H. W., Chen K. F., Meng B. (2017). A case study on large deformation failure mechanism of deep soft rock roadway in Xin'An coal mine, China. Engineering Geology 217: 89-101. [ Links ]

[5] Park K. H., Kim Y. J. (2006). Analytical solution for a circular opening in an elastic-brittle-plastic rock. International Journal of Rock Mechanics and Mining Sciences 43: 616-622. [ Links ]

[6] Sharan S. K. (2003). Elastic-brittle-plastic analysis of circular openings in Hoek-Brown media. International Journal of Rock Mechanics and Mining Sciences 40: 817-824. [ Links ]

[7] Sharan S. K. (2008). Analytical solutions for stresses and displacements around a circular opening in a generalized Hoek-Brown rock. International Journal of Rock Mechanics and Mining Sciences 40: 78-85. [ Links ]

[8] Jiang B. S., Zhang Q., He Y. N., Han L. J. (2007). Elastoplastic analysis of cracked surrounding rocks in deep circular openings. Chinese Journal of Rock Mechanics and Engineering 26: 982-986 (in Chinese). [ Links ]

[9] Stille H., Holmberg M., Nord G. (1989). Support of weak rock with grouted bolts and shotcrete. International Journal of Rock Mechanics & Mining Sciences & Geomechanics Abstracts 26: 99-113. [ Links ]

[10] Brown E. T., Bray J. W., Ladanyi B., Hoek E. (1983). Ground response curves for rock tunnels. Journal of Geotechnical Engineering 109: 15-39. [ Links ]

[11] Perzyna P. (1966). Fundamental problems in viscoplasticity. Advances in applied mechanics 9: 244-368. [ Links ]

[12] Lemaitre, J., Chaboche, J. L. (1985). Me’canique desmate’riaux solides. Dunod, Paris. [ Links ]

[13] Pellet F., Hajdu A., Deleruyelle F., Besnus F. (2005). A viscoplastic constitutive model including anisotropic damage for the time dependent mechanical behaviour of rock. International Journal for Numerical and Analytical Methods in Geomechanics 29: 941-970. [ Links ]

[14] Boidy E., Bouvard A., Pellet F. (2002). Back analysis of time-dependent behaviour of a test gallery in claystone. Tunnelling and Underground Space Technology 17: 415-424. [ Links ]

[15] Bonini M., Debernardi D., Barla M., Barla G. (2009). The Mechanical Behaviour of Clay Shales and Implications on the Design of Tunnels. Rock Mechanics & Rock Engineering 42: 361-388. [ Links ]

[16] Pellet F. (2009). Contact between a Tunnel Lining and a Damage-Susceptible Viscoplastic Medium. Computer Modeling in Engineering & Sciences 52: 279-295. [ Links ]

[17] Barla G., Debernardi D., Sterpi D. (2012). Time-Dependent Modeling of Tunnels in Squeezing Conditions. International Journal of Geomechanics 12: 697-710. [ Links ]

[18] Debernardi D., Barla G. (2009). New Viscoplastic Model for Design Analysis of Tunnels in Squeezing Conditions. Rock Mechanics & Rock Engineering 42: 259-288. [ Links ]

[19] Sterpi D., Gioda G. (2009). Visco-Plastic Behaviour around Advancing Tunnels in Squeezing Rock. Rock Mechanics & Rock Engineering 42: 319-339. [ Links ]

[20] Roatesi S. (2013). Analytical approach of the tunnel face effect on the rock-support interface in viscoelastic rock mass. American Institute of Physics: 1354-1358. [ Links ]

[21] Sulem J.,Panet M.,Guenot A. (1987). An analytical solution for time-dependent displacements in a circular tunnel. International Journal of Rock Mechanics & Mining Sciences & Geomechanics Abstracts 24: 155-164. [ Links ]

[22] Ladanyi B. (1980). Direct determination of ground pressure on tunnel lining in a non-linear viscoelastic rock underground rock engineering. Proceedings of the Thirteenth Canadian Rock Mechanics Symposium, Montreal: Canadian Institute of Mining and Metallurgy: 126-132. [ Links ]

[23] Gill D. E., Ladanyi B. (1987). Time-dependent ground response curves for tunnel lining design. Proceedings of 6th International Congress on Rock Mechanics, Montreal: Taylor & Francis 2: 917-921. [ Links ]

[24] Wang H. N., Nie G. H. (2010). Analytical expressions for stress and displacement fields in viscoelastic axisymmetric plane problem involving time-dependent boundary regions. Acta Mechanica 210: 315-330. [ Links ]

[25] Fritz P. (2010). An analytical solution for axisymmetric tunnel problems in elasto-viscoplastic media. International Journal for Numerical & Analytical Methods in Geomechanics 8: 325-342. [ Links ]

[26] Zhang Q., Jiang B. S., Wu X. S., Zhang H. Q., Han L. J. (2012). Elasto-plastic coupling analysis of circular openings in elasto-brittle-plastic rock mass. Theoretical & Applied Fracture Mechanics 60: 60-67. [ Links ]

[27] Zhang Q., Jiang B. S., Wang S. L., Ge X. R., Zhang H. Q. (2012). Elasto-plastic analysis of a circular opening in strain-softening rock mass. International Journal of Rock Mechanics & Mining Sciences 50: 38-46. [ Links ]

[28] Zhang Q., Zhang C. H., Jiang B. S., Li N., Wang Y. C. (2018). Elastoplastic Coupling Solution of Circular Openings in Strain-Softening Rock Mass Considering Pressure-Dependent Effect. International Journal of Geomechanics 18: 04017132. [ Links ]

[29] Lv A. Z., Xu G. S., Sun F., Sun W. Q. (2010). Elasto-plastic analysis of a circular tunnel including the effect of the axial in situ stress. International Journal of Rock Mechanics and Mining Sciences 47: 50-59. [ Links ]

[30] Yang S. Q., Xu P., Ranjith P. G. (2015). Damage model of coal under creep and triaxial compression. International Journal of Rock Mechanics & Mining Sciences 80: 337-345. [ Links ]

Received: June 21, 2019; Revised: June 24, 2019; Accepted: June 29, 2019; Published: July 02, 2019

* Corresponding author

Author's Contributions:

Conceptualization, X Wang and BS Jiang; Methodology, X Wang; Funding acquisition, MM Lu; Formal analysis, X Wang and Q Zhang; Software, X Wang, MM Lu and M Chen; Writing-original draft, X Wang, BS Jiang and Q Zhang.


Marcílio Alves

Creative Commons License This is an open-access article distributed under the terms of the Creative Commons Attribution License