Analytical solution of circular tunnel in elastic-viscoplastic rock mass

Abstract

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 [¹]. 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 [²]. 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 [³]. Typical rheological failures of the roadway are shown in Fig. 1.

Fig. 1
Rheological failures of laneway [⁴]

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 [⁵] 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 [⁶-⁷] derived the closed-form expressions of stress and deformation in the elasto-brittle-plastic Hoek-Brown rock mass for circular tunnels. Jiang et al. [⁸] 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. [⁹] 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. [¹⁰] 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 [¹¹]. Subsequently, Perzyna’s general theory of elasto-viscoplasticity was applied to metal alloys by Lemaitre and Chaboche [¹²], and was extended by Pellet et al. [¹³] to the study of the time-dependent behavior of rocks. Boidy et al. [¹⁴] and Bonini et al. [¹⁵] combined numerical simulations with Lemaitre’s viscoplastic model. In turn, Pellet [¹⁶] combined the Burger and the elastoplastic models, and calibrated the constitutive parameters using back-analyses of convergence data.

Barla et al. [¹⁷], Debernardi et al. [¹⁸], and D. Sterpi et al. [¹⁹], have conducted corresponding studies of deep tunnels using numerical simulations at squeezing conditions. Roatesi [²⁰] 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 [²¹], Ladanyi [²²], Gill [²³], and Wang [²⁴] 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 [²⁵] 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 [²⁶]

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 ₀ 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 ${\sigma }_r$, ${\epsilon }_r$, ${\sigma }_{\theta }$ and ${\epsilon }_{\theta }$ 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)

$\frac{\partial {\sigma }_r}{\partial r}=\frac{{\sigma }_{\theta }-{\sigma }_r}{r}$ (1)

The geometric equation of this axisymmetric problem can be denoted as

$\{\begin{array}{l}{\epsilon }_r=\frac{\partial u}{\partial r}\\ {\epsilon }_{\theta }=\frac{u}{r}\end{array}$ (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={\sigma }_{\theta }-N{\sigma }_r-S=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:

$\{\begin{array}{l}{N}_p=(1+\sin {\phi }_p)/(1-\sin {\phi }_p)\\ S_p=2c_p\cos {\phi }_p/(1-\sin {\phi }_p)\end{array}$ (3a)

$\{\begin{array}{l}{N}_c=(1+\sin {\phi }_c)/(1-\sin {\phi }_c)\\ S_c=2c_c\cos {\phi }_c/(1-\sin {\phi }_c)\end{array}$ (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 [²⁷, ²⁸] here is:

$\Phi ={\sigma }_{\theta }-\alpha {\sigma }_r$ (4)

where $\alpha =\left(1+\sin {\psi }_c\right)/\left(1-\sin {\psi }_c\right)$ is material parameter, ${\psi }_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:

$\begin{array}{cc}{\epsilon }_r^p=\lambda \frac{\partial \Phi }{\partial {\sigma }_r},& {\epsilon }_{\theta }^p=\lambda \frac{\partial \Phi }{\partial {\sigma }_{\theta }}\end{array}$ (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

$\begin{array}{l}{\epsilon }_r=[(1-u)({\sigma }_r-p_0)-\mu ({\sigma }_{\theta }-p_0)]/(2G)-\lambda \alpha \\ {\epsilon }_{\theta }=[(1-u)({\sigma }_{\theta }-p_0)-\mu ({\sigma }_r-p_0)]/(2G)+\lambda \end{array}\}$ (6)

Where $G$ is the shear modulus, $\mu$ is poisson ratio, $\lambda$ is plastic strain. When $\lambda =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${\sigma }_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:

$\dot{\epsilon }=\frac{\dot{\sigma }}{E}+\frac{\sigma -{\sigma }_s}{\eta }$ (7)

where dot “ ^{∙ “} denotes the derivative of time, $E$ is elastic modulus, ${\sigma }_s$ is yield strength of plastic elements, $\eta$ is viscosity coefficient of Newtonian element.

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

$\begin{array}{l}{\dot{e}}_{ij}=\frac{{\dot{s}}_{ij}}{2G}+\frac{s_{ij}}{2\eta }\hfill \\ {\dot{\epsilon }}_{kk}=\frac{{\dot{\sigma }}_{kk}}{3K}+\frac{{\sigma }_{kk}}{2\eta \text{'}}\hfill \end{array}\}$ (8)

where $s_{ij}$ is deviatoric stress tensor, $e_{ij}$ is deviatoric strain tensor, ${\sigma }_{ij}$ is stress tensor, ${\epsilon }_{ij}$ is strain tensor, both stress and strain tensor can be decomposed into deviatoric stress and hydrostatic pressure, namely, $s_{ij}={\sigma }_{ij}-{\sigma }_{kk}{\delta }_{ij}/3$, $e_{ij}={\epsilon }_{ij}-{\epsilon }_{kk}{\delta }_{ij}/3$, ${\sigma }_{kk}$ and ${\epsilon }_{kk}$ indicate the first invariants of stress and strain. ${\delta }_{ij}$ is Kroneker symbol, when $i=j$, ${\delta }_{ij}=1$, when $i\ne j$, ${\delta }_{ij}=0$. G is shear modulus, and K is bulk modulus, $\eta$ is viscosity coefficient of deviatoric stress and deviatoric strain, $\eta \text{'}$ is viscosity coefficient of hydrostatic pressure (${\sigma }_{kk}$) and volume strain (${\epsilon }_{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:

$\begin{array}{l}{\dot{\epsilon }}_r=\frac{1}{2G}[{\dot{\sigma }}_r-(\frac{1}{3}-\frac{2G}{9K})\dot{\Theta }]+\frac{1}{2\eta }[{\sigma }_r-(\frac{1}{3}-\frac{\eta }{3{\eta }^{\prime }})\dot{\Theta }]\hfill \\ {\dot{\epsilon }}_{\theta }=\frac{1}{2G}[{\dot{\sigma }}_{\theta }-(\frac{1}{3}-\frac{2G}{9K})\dot{\Theta }]+\frac{1}{2\eta }[{\sigma }_{\theta }-(\frac{1}{3}-\frac{\eta }{3{\eta }^{\prime }})\dot{\Theta }]\hfill \\ {\dot{\epsilon }}_z=\frac{1}{2G}[{\dot{\sigma }}_z-(\frac{1}{3}-\frac{2G}{9K})\dot{\Theta }]+\frac{1}{2\eta }[{\sigma }_z-\frac{1}{3}(1-\frac{\eta }{{\eta }^{\prime }})\dot{\Theta }]\hfill \end{array}\}$ (9)

where $\Theta ={\sigma }_r+{\sigma }_{\theta }+{\sigma }_z$ is the volume stress.

By substituting plane strain condition (${\epsilon }_z={\dot{\epsilon }}_z=0$) into the third equation of Eq. (9), axial stress can be written as:

$\begin{array}{l}{\sigma }_z=\frac{\eta \text{'}-\eta }{2\eta \text{'}+\eta }({\sigma }_r+{\sigma }_{\theta })\hfill \\ {\dot{\sigma }}_z=\frac{3K-2G}{6K+2G}({\dot{\sigma }}_r+{\dot{\sigma }}_{\theta })\hfill \end{array}\}$ (10)

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

$\eta \text{'}=\frac{3K}{2G}\eta$ (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, $(3K-2G)/(6K+2G)$ is replaced by the Poisson ratio $\mu$, the three-dimensional constitutive equation of generalized Bingham body under the condition of plane strain is obtained:

$\begin{array}{c}{\dot{\epsilon }}_r=\frac{1}{2G}[(1-\mu ){\dot{\tilde{\sigma }}}_r-\mu {\dot{\tilde{\sigma }}}_{\theta }]+\frac{1}{2\eta }[(1-\mu ){\tilde{\sigma }}_r-\mu {\tilde{\sigma }}_{\theta }]\\ {\dot{\epsilon }}_{\theta }=\frac{1}{2G}[(1-\mu ){\dot{\tilde{\sigma }}}_{\theta }-\mu {\dot{\tilde{\sigma }}}_r]+\frac{1}{2\eta }[(1-\mu ){\tilde{\sigma }}_{\theta }-\mu {\tilde{\sigma }}_r]\end{array}\}$ (12a)

$\begin{array}{l}{\tilde{\sigma }}_r={\sigma }_r-{\sigma }_r^c\hfill \\ {\tilde{\sigma }}_{\theta }={\sigma }_{\theta }-{\sigma }_{\theta }^c\hfill \end{array}\}$ (12b)

where ${\tilde{\sigma }}_r$ and ${\tilde{\sigma }}_{\theta }$ 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:

$\begin{array}{l}r=a,{\sigma }_r=p_i\\ r=R_c,{\sigma }_r(R_c{}^{-},t)={\sigma }_r(R_c{}^{+}),u(R_c{}^{-},t)=u(R_c{}^{+})\\ r\to \infty ,{\sigma }_r=p_0\\ t=\mathrm{0,}{\sigma }_r={\sigma }_{r0}^e,{\sigma }_{\theta }={\sigma }_{\theta 0}^e,u=u_0^e\end{array}\}$ (13)

where $R_c^{-}$and $R_c^{+}$ represents the left and right boundary of r=$R_c$, 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 [⁸, ²⁶].

3.1 Cracked zone (a≤r≤R _c )

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

$\begin{array}{l}{\sigma }_r^c=p_i^c{(}^r-D_c\\ {\sigma }_{\theta }^c=N_c{p}_i^c{(}^r-D_c\end{array}\}$ (14)

where $D_c=C_c\cot {\varphi }_c$, $p_i^c=p_i+D_c$.

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:

$\begin{array}{l}{\epsilon }_r^c=[N_c{B}_1{(}^r-B_2+A{\alpha }_c{(}^{R_c}]/(2G)\\ {\epsilon }_{\theta }^c=[B_1{(}^r-B_2-A{(}^{R_c}]/(2G)\\ u_c=r{\epsilon }_{\theta }^c\end{array}\}$ (15)

where $B_1=[(1-\mu )(N_c{\alpha }_c+1)/(N_c+{\alpha }_c)-\mu ](p_i+C_c\cot {\phi }_c)$,$B_2=(1-2\mu )(p_0+C_c\cot {\phi }_c)$, A is an integral constant.

3.2 Elastic zone (Rc ≤ r <∞)

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

$\begin{array}{l}{\sigma }_r^e=p_0+(p_c-p_0)R_c^2{r}^{-2}\\ {\sigma }_{\theta }^e=p_0-(p_c-p_0)R_c^2{r}^{-2}\end{array}\}$ (16)

$\begin{array}{l}{\epsilon }_r^e=(p_c-p_0)R_c^2{r}^{-2}/(2G)\hfill \\ {\epsilon }_{\theta }^e=(p_0-p_c)R_c^2{r}^{-2}/(2G)\hfill \\ u^e=r{\epsilon }_{\theta }^e\hfill \end{array}\}$ (17)

where $p_c$ is radial stress of the interface between the elastic zone and the cracked zone, $R_c$ is the radius of the interface, they can be represented as [²⁹]:

$p_c=(p_i+C_c\cot {\phi }_c){(R_c/a)}^{N_c-1}-C_c\cot {\phi }_c$ (18)

$R_c=a{\left\{\frac{p_0(1-\sin {\phi }_p)-C_p\cos {\phi }_p+C_c\cot {\phi }_c}{p_i+C_c\cot {\phi }_c}\right\}}^{1/(N_c-1)}$ (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-\mu )\frac{(N_c+1)({\alpha }_c+1)}{N_c+{\alpha }_c}(p_i+C_c\cot {\phi }_c){\left(\frac{R_c}{a}\right)}^{N_c-1}-2(1-\mu )(p_0+C_c\cot {\phi }_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 $R_t$ 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:

$\begin{array}{c}{\sigma }_{r0}^e=p_0+(p_i-p_0)\frac{a^2}{r^2}\\ {\sigma }_{\theta 0}^e=p_0-(p_i-p_0)\frac{a^2}{r^2}\\ u_0^e=\frac{1}{2G}(p_0-p_i)\frac{a^2}{r}\end{array}\}$ (21)

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

Through analyzing rheological test data [³⁰], 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:

${\rho }_0=a\sqrt{\frac{(1+N_p)(p_i-p_0)}{(1-N_p)p_0-S_p}}$ (22)

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

$r\frac{\partial }{\partial r}\left[{\dot{\tilde{\sigma }}}_r+{\dot{\tilde{\sigma }}}_{\theta }+\beta ({\tilde{\sigma }}_r+{\tilde{\sigma }}_{\theta })\right]=0$ (23)

where $\beta =G/\eta$.

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

${\tilde{\sigma }}_r+{\tilde{\sigma }}_{\theta }=f(t)+e^{-\beta t}g(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:

$\begin{array}{c}{\sigma }_r={\sigma }_r^c(1-e^{-\beta t})+{\sigma }_{r0}^e{e}^{-\beta t}+(1-a^2{r}^{-2})f(t)\\ {\sigma }_{\theta }={\sigma }_{\theta }^c(1-e^{-\beta t})+{\sigma }_{\theta 0}^e{e}^{-\beta t}+(1+a^2{r}^{-2})f(t)\end{array}\}$ (25)

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

$\begin{array}{l}u=r/(2G)\{\left[(1-\mu ){\sigma }_{\theta 0}^e-\mu {\sigma }_{r0}^e\right]e^{-\beta t}+\left[(1-\mu ){\sigma }_{\theta }^c-\mu {\sigma }_r^c\right](1-e^{-\beta t})\\ +(1-2\mu +a^2{r}^{-2})f(t)\}+r/(2\eta )\left\{\left[(1-\mu )({\sigma }_{\theta 0}^e-{\sigma }_{\theta }^c)-\mu ({\sigma }_{r0}^e-{\sigma }_r^c)\right](1-e^{-\beta t})/\beta \right\}\\ +r/(2\eta )(1-2\mu +a^2{r}^{-2}){\displaystyle \int f(t)}dt+g_1(r)\end{array}$ (26)

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

$g_1(r)=r(2\mu -1)p_0/(2G)$ (27)

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

$\begin{array}{l}u=r/(2G)\{\left[(1-\mu ){\sigma }_{\theta 0}^e-\mu {\sigma }_{r0}^e\right]+(1-2\mu +a^2{r}^{-2})f(t)\}\\ +r/(2\eta )(1-2\mu +a^2{r}^{-2}){\displaystyle \int f(t)}dt+r(2\mu -1)p_0/(2G)\end{array}$ (28)

The above equation, ${\sigma }_{r0}^e$ and ${\sigma }_{\theta 0}^e$ 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. ${\rho }_0\le r\le R_t$). This hypothesis affects the evolution of stress and strain with time. When the cracked zone extends to the radius $R_t$ 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 $p_i$ is replaced by ${\rho }_0$ and $P(t)$ respectively. The forms of stress, strain and displacement after transformation are shown as:

$\begin{array}{l}{\sigma }_r^{vp}=(P(t)+C_c\cot {\phi }_c){(r/{\rho }_0)}^{N_c-1}-C_c\cot {\phi }_c\hfill \\ {\sigma }_{\theta }^{vp}=N_c(P(t)+C_c\cot {\phi }_c){(r/{\rho }_0)}^{N_c-1}-C_c\cot {\phi }_c\hfill \end{array}\}$ (29)

$u^{vp}=r\left[B_1^{vp}{(r/{\rho }_0)}^{N_c-1}-B_2-A^{vp}{(R_t/r)}^{{\alpha }_c+1}\right]/(2G)$ (30)

where, $B_1^{vp}$, $A^{vp}$ and $R_t$, as shown in Eq. (31). In these equations, the radius and the stress are also replaced accordingly.

$\begin{array}{l}{B}_1^{vp}=[(1-\mu )(N_c{\alpha }_c+1)/(N_c+{\alpha }_c)-\mu ](P(t)+C_c\cot {\phi }_c)\\ A^{vp}=(1-\mu )\frac{(N_c+1)({\alpha }_c+1)}{N_c+{\alpha }_c}(P(t)+C_c\cot {\phi }_c){\left(\frac{R_t}{a}\right)}^{N_c-1}-2(1-\mu )(p_0+C_c\cot {\phi }_c)\\ R_t={\rho }_0{\left\{\frac{p_0(1-\sin {\phi }_p)-C_p\cos {\phi }_p+C_c\cot {\phi }_c}{P(t)+C_c\cot {\phi }_c}\right\}}^{1/(N_c-1)}\end{array}\}$ (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 $R_t$, 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)={\sigma }_r^c(1-e^{-\beta t})+{\sigma }_{r0}^e{e}^{-\beta t}+(1-a^2{\rho }_0{}^{-2})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\to \infty$, 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={\rho }_0$, i.e. Eq. (28) and Eq. (30) at the boundary $r={\rho }_0$ are equal, an integral equation about$f(t)$is obtained:

$\begin{array}{l}\left[\left(1-\mu \right)\left(N_c{\alpha }_c+1\right)/\left(N_c+{\alpha }_c\right)-\mu \right]\left(P(t)+C_c\cot {\phi }_c\right)\\ -(1-\mu )[\frac{\left(N_c+1\right)\left({\alpha }_c+1\right)}{N_c+{\alpha }_c}(P(t)+C_c\cot {\phi }_c){\left(\frac{R_t}{{\rho }_0}\right)}^{N_c+{\alpha }_c}-2(p_0+C_c\cot {\phi }_c){\left(\frac{R_t}{{\rho }_0}\right)}^{{\alpha }_c+1}]\\ -\left(1-2\mu \right)C_c\cot {\phi }_c+\mu {\sigma }_{r0}^e-\left(1-\mu \right){\sigma }_{\theta 0}^e-\left(1-2\mu +a^2{\rho }_0^{-2}\right)f\left(t\right)\\ -\left(1-2\mu +a^2{\rho }_0^{-2}\right)\beta {\displaystyle \int f\left(t\right)dt}=0\end{array}$ (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 ${\displaystyle \int 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 ${\displaystyle \int 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 ${\displaystyle \int f(t)}$ can be obtained, as shown in Figs. 6 and 7. Fig. 6 denotes the curve of $f(t)$ with time, with$(Et)/\eta$as abscissa and$f(t)/p_0$as ordinate; Fig. 7 shows the curve of ${\displaystyle \int f(t)}$with time, with$(Et)/\eta$as abscissa and$(E{\displaystyle \int f(t)})/(\eta p_0)$as ordinate.

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

Fig. 7
Analytic solution of function${\displaystyle \int f(t)}$under different support force

It is obvious that$f(t)$and${\displaystyle \int 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${\displaystyle \int 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 $R_t$ 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 [²⁶].

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 ${\rho }_0$. Up to the actual radius $R_t$ of the plastic zone (i.e. for ${\rho }_0\le r\le R_t$), expressions (31) are relevant and for $r\ge R_t$ equations (16) and (17) are relevant, whereby $p_c$, $R_c$ and the index c have to be replaced by $P(t)$, $R_t$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.

Acknowledgments

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

References

• ^[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).
• ^[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.
• ^[3]
  Sun J. (1999). Rheology and engineering application of geotechnical materials. Beijing: China Construction Industry Press.
• ^[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.
• ^[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.
• ^[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.
• ^[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.
• ^[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).
• ^[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.
• ^[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.
• ^[11]
  Perzyna P. (1966). Fundamental problems in viscoplasticity. Advances in applied mechanics 9: 244-368.
• ^[12]
  Lemaitre, J., Chaboche, J. L. (1985). Me’canique desmate’riaux solides. Dunod, Paris.
• ^[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.
• ^[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.
• ^[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.
• ^[16]
  Pellet F. (2009). Contact between a Tunnel Lining and a Damage-Susceptible Viscoplastic Medium. Computer Modeling in Engineering & Sciences 52: 279-295.
• ^[17]
  Barla G., Debernardi D., Sterpi D. (2012). Time-Dependent Modeling of Tunnels in Squeezing Conditions. International Journal of Geomechanics 12: 697-710.
• ^[18]
  Debernardi D., Barla G. (2009). New Viscoplastic Model for Design Analysis of Tunnels in Squeezing Conditions. Rock Mechanics & Rock Engineering 42: 259-288.
• ^[19]
  Sterpi D., Gioda G. (2009). Visco-Plastic Behaviour around Advancing Tunnels in Squeezing Rock. Rock Mechanics & Rock Engineering 42: 319-339.
• ^[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.
• ^[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.
• ^[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.
• ^[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.
• ^[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.
• ^[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.
• ^[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.
• ^[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.
• ^[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.
• ^[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.
• ^[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.

Edited by

Editor:

Marcílio Alves

Publication Dates

• Publication in this collection
  29 July 2019
• Date of issue
  2019

History

• Received
  21 June 2019
• Reviewed
  24 June 2019
• Accepted
  29 June 2019
• Published
  02 July 2019