Elastic-viscoplastic solution of a circular tunnel in Hoek-Brown rock mass

ABSTRACT

A closed-form analytical solution to rheological phenomena is obtained based on Hoek-Brown (HB) yield criterion and the generalized Bingham model. Once the circular tunnel is excavated, an initial plastic region is formed. When the total stresses fulfil the yield condition, the initial plastic region steps into the viscoplastic region, in which the displacement develops with time. The one-dimensional constitutive equation of the generalized Bingham model is transformed into a three-dimensional expression in polar coordinates by combining the plane strain condition. With the help of equilibrium differential equation, geometric equation and boundary condition, the stress and displacement in the initial plastic region can be expressed. At the moment of t, the stresses and displacements are described similar to the solutions of the elasto-plastic state. The stresses, strains and displacements after transformation can be expressed in the viscoplastic region. The proposed solution was verified by the closed-form solution of circular tunnels in elastic-brittle-plastic rock mass.

KEYWORDS:
Circular tunnel; Brittle-plastic; Viscoplastic; Analytical; Closed-form solution

Graphical Abstract

List of Symbols
• a  Excavation radius
• σ_r , σ _θ  Radial and tangential stresses
• ε_r , ε _θ  Radial and tangential strains
• ${\sigma }_r^{vp}$, ${\sigma }_{\vartheta }^{vp}$  Radial and tangential viscoplastic stresses
• ${\epsilon }_r^{vp}$, ${\epsilon }_{\vartheta }^{vp}$  Radial and tangential viscoplastic strains
• u  Radial displacement
• $u^{vp}$  Radial viscoplastic displacement
• E  Young’s modulus
• G  Shear modulus
• σ _c  Uniaxial compression strength
• μ  Poisson’s ratio
• ψ  Dilatancy angle
• c  Cohesive strength
• φ  Friction angle
• m, s  The first and section strength parameters
• p _i  Internal supporting pressure
• p ^c  Radial stress at elasto-plastic interface
• ρ ₀  The radius of the initial plastic region
• R _t  The radius of the viscoplastic region
•  Derivative of variable to time

1. Introduction

Because of the urgent need of mineral resources, the mining dramatically deep into the earth, many coal and metal mines have been mined to a depth of more than 1 km. In order to facilitate the passage, deep buried tunnel construction has been common. However, deep buried roadway is located in the harsh environment of high confining pressure, high temperature and high permeability, which is easy to produce engineering geological disasters such as rheology, rock burst and large deformation of surrounding rock. Especially in the soft rock roadway, rock lithology metamorphism and surrounding rock failure occur frequently, the construction process is full of challenges [¹]. This will seriously affect the stability of roadway surrounding rock, cause safety hazards to construction and traffic personnel, and cause huge economic losses to the roadway repeatedly renovated and maintained. Correctly predicting the deformation, stress and failure mode of surrounding rock can provide reliable reference for the stability evaluation of underground engineering. The expansion of surrounding rock fracture range is time dependent. This leads to the nonlinear continuous deformation of surrounding rock, which makes it more difficult to maintain the stability of surrounding rock. For soft rock mass, the creep displacement is significant, and the deformation of soft rock mass can reach tens of millimeters every day, and when the support strength is weak, the deformation can reach tens of centimeters every day. In the 1960s, Lecian’s concept of aging was first introduced into the analysis of stress and strain around tunnels [²]. An increasing number of researchers and engineers are paying attention to the time-dependent characteristics of the surrounding rock in underground engineering. Rheological research provides a reference for support design in rock foundations, slopes, and underground engineering. The deformation of these projects is time-dependent, and the deformation of the surrounding rock increases with time. After a certain period of time, large cumulative deformation results in damage to the tunnel support structure, and rock rheology is closely related to the long-term stability of geotechnical engineering. In deep underground engineering, due to high stress, the rheologic phenomena of the rock mass are more obvious, and research on the mechanical state of the rheological rock mass is of great value to engineering stability and safety design. Therefore, it is necessary to derive analytical solutions considering rheological factors.

Axisymmetric analysis of a circular tunnel is a classical problem. Many studies in existing research reports have been performed on the analytical solutions after excavation, and many engineering problems have been solved. In 1983, Brown et al. [³] proposed a method for analyzing the distribution of the stress and displacement in the plastic region into surrounding rock by using the Hoek-Brown (HB) criterion. Sharan, S. K. [⁴-⁵] presented stress and displacement field analysis and analytical solutions of the rock surrounding tunnels under the HB criterion. Carranza et al. [⁶-⁷] presented a closed multidimensional elastoplastic solution of the rock surrounding tunnels under the HB criterion and elasto-plastic solutions of the rock surrounding tunnels based on the Mohr-Coulomb (MC) and HB criteria, respectively. Zhang et al. [⁸] dealt with elasto-plastic coupling solutions for the prediction of displacements around circular openings in elasto-brittle-plastic rock masses. Wang et al. [⁹] divided the strain softening region into k infinitesimal rings, and this method considered the effect of parameter degradation in the strain softening process. Jiang et al. [¹⁰] divided the surrounding rock into the elastic region, the plastic region, and the plastic region, and obtained the stress and deformation solutions of each region. Wu et al. [¹¹] established a three-dimensional twin tunnels scale model utilizing the discrete element method (DEM) with PFC3D, and investigated the displacement (in horizontal and vertical directions) and deformation of the first tunnel lining in four different cases, in which the clear distance of twin tunnels is 5, 10, 15 and 20 m during the second tunnel construction process. These above solutions are instantaneous without considering the rheological factors. However, in engineering practice, the rheological effect is often an important factor for the deformation and stability evaluation of underground engineering.

In addition to elastic-plastic research, many scholars have carried out extensive research on the rheological properties of the surrounding rocks. H. N. Wang et al. [¹²] derived the analytical solution of stress and displacement of lining tunnel by viscoelastic method. P. Nomikos et al. [¹³] derived the analytical formulas of stress and displacement of an axisymmetric tunnel with the Burgers model under elastic lining support, and analyzed the influence of viscoelastic parameters. A. Fahimifar et al. [¹⁴] derived the analytical solution of tunnel wall deformation under hydrostatic pressure and predicted the time-dependent displacement of the tunnel wall, assuming that the rock mass is an incompressible Burgers viscoelastic material. The above analysis does not take into account the plastic deformation of the rock mass and cannot fully explain the behavior of underground engineering. D. Sterpi and G. Gioda [¹⁵] used rheological models that consider elastic-plastic, viscoelastic and viscoplastic properties to simulate three-stage creep behavior of rocks in squeezing rock of deep-buried tunnels. Panet [¹⁶] and J. Sulem [¹⁷] analyzed the aging characteristics of a circular tunnel and considered the plastic mechanical properties of the surrounding rock, overcoming the disadvantage that the viscoelastic model cannot calculate the influence of strain softening and plastic dilatancy of the surrounding rock on rheological deformation. M. E. Gurtin and E. Sternberg [¹⁸] presented the application of Stieltjes integration to the treatment of discontinuous boundary conditions in viscoelastic problems. G. A. C. Graham [¹⁹] studied the viscoelastic solution of a thick-walled cylinder in temperature and stress fields. Nevertheless, in the solution of the above theories, the researchers all assumed that the rheological properties of the surrounding rocks follow Kelvin's viscoelastic equation and then introduced the MC yield criterion. P. Fritz [²⁰] presented an analytical solution of a modified St. Venant slider, and the viscoplastic constitutive equations were obtained by a similar method. Based on the research of these scholars, the evolution law of stress and displacement with radius and time is analyzed by using generalized Bingham model and strict theory of solid mechanics.

When analyzing the rheological behavior of the rock mass in deep tunnel engineering, the rheological constitutive model and mechanical calculation model are two key factors affecting the results. In analytical research, the rheological constitutive model has been widely used. Empirical models, component models, internal variable models without yield surface, and rheological constitutive models [²¹] based on damage mechanisms are widely used in the literature. The generalized Bingham model used in this paper is one of the component models. In many practical cases, especially in jointed rock masses, the nonlinear HB empirical yield criterion is more suitable. It can better reflect the inherent characteristics and nonlinear failure characteristics of the surrounding rock in deep tunnels, as well as the influence of the number of structural planes, rock strength, and stress state on the strength of the surrounding rocks [²²].

By combining the HB yield criterion and generalized Bingham elasto-viscoplastic constitutive model, a three-dimensional constitutive equation can be established. The analytical solutions of the elasto-viscoplastic stress and displacement are obtained by MATLAB software. The results are compared with those of the MC criterion under equivalent parameters, the analytical solutions in this paper are compared with the previous results. Finally, the influences of associated parameters, such as the support force and dilation angle, are systematically investigated.

2. Problem description

Although the rock mass is actually discontinuous medium, Fazelabdolabadi et al. [²³] develops a Bayesian framework to quantify the absolute permeability of water in a porous structure from the geometry and clustering parameters of its underlying pore-throat network. Herein, A circular tunnel with a radius a is excavated in a homogeneous infinite isotropic rock mass subjected to a hydrostatic stress (p ₀), the influence of porosity and permeability is not considered. The inner surface of the opening is subjected to a support pressure p _i, and the inner support pressure p _i equals the hydrostatic stress p ₀ before the excavation. The elasto-viscoplastic circular tunnel is shown in Figure 1. The analytical model of the surrounding rock consists of three concentric annuli, i.e., the elastic region, the viscoplastic region, and the initial plastic region.

Fig. 1
Physical model of a circular opening.

To derive the analytical solution, the following assumptions are necessary: the rock mass is isotropic, uniform, and homogeneous; the plane strain condition is considered and the problem illustrated in Figure 1 can be considered to be axisymmetric. For plane strain problems, it can be assumed that the radial stress σ _r is the minor principal stress, and the tangential stress σ _θ is the major principal stress. Additionally, the strains are assumed to be small. The notation u represents the radial displacement, and the superscripts “c” and “e” represent plastic region and elastic region, respectively. The displacements and stresses are functions of the distance r to the center of the cross section and the time t when gravity is ignored.

The support pressure p _i gradually decreases, and the deformation of the rock mass converges, when the roadway is excavated. When p _i is lower than the elasto-plastic support pressure p ^c, the plastic region begins to appear outside the excavation face. After that, the material properties of the plastic region reach the residual value. Therefore, the macro brittle failure of the rock is the result of crack initiation, propagation, and connection in the medium. The elasto-brittle-plastic model is shown in Figure 2, where AD represents the ideal plastic model, ACE represents the actual stress-strain curve of the rock, and ABE represents the assumed brittle-plastic model, the assumed brittle plastic model is used to analyze the stress and deformation of surrounding rock.

Fig. 2
Mechanical behavior model for elastic-btittle-plastic material.

For the problem of circular openings, the tangential stress (σ _θ) and radial stress (σ _r) are the major and minor stresses in circular tunnels, respectively. In this paper, two of the most commonly-used yield criteria are considered, the linear MC criterion

${\sigma }_1=N{\sigma }_3+S$ (1)

and the nonlinear HB criterion

${\sigma }_1={\sigma }_3+\sqrt{m{\sigma }_s{\sigma }_3+s{\sigma }_s^2}$ (2)

where σ ₃ is the minor principal stress, σ ₁ is the major principal stress, $N=(1+\sin \phi )/(1-\sin \phi )$, $S=2c\cos \phi /(1-\sin \phi )$, $c$ and φ are the cohesion and friction angle of the rock mass, respectively, m and s are material parameters that depend on the properties of the rock, and σ _s is the uniaxial compressive strength of the intact rock mass.

The yield criteria of Eqs. (1) and (2) for both intact and residual rock masses can be expressed as, where Eq. (3) and Eq. (5) represent MC and HB intact rock mass, respectively; Eq. (4) and Eq. (6) represent MC and HB residual rock mass, respectively.

${\sigma }_{\theta }=N{\sigma }_r+Sfor\mathrm{int}actrockmass$ (3)

${\sigma }_{\theta }=N_c{\sigma }_r+S_{c}forresidualrockmass$ (4)

${\sigma }_{\theta }={\sigma }_r+\sqrt{m{\sigma }_s{\sigma }_r+s{\sigma }_s^2}for\mathrm{int}actrockmass$ (5)

${\sigma }_{\theta }={\sigma }_r+\sqrt{m_c{\sigma }_s{\sigma }_r+s_c{\sigma }_s^2}forresidualrockmass$ (6)

where the $N_c$, $S_c$, $m_c$ and $s_c$ are residual values of the rock properties, ${\sigma }_s$ is the uniaxial compressive strength and $N_c=(1+\sin {\phi }_c)/(1-\sin {\phi }_c)$, and $S_c=2c_c\cos {\phi }_c/(1-\sin {\phi }_c)$, in which ${\phi }_c$ and $c_c$ are the residual friction angle and cohesion, respectively.

For axisymmetric problems, the stresses satisfy the differential equation of equilibrium in each region (without considering the body force of the surrounding rock), as shown in Eq. (7)

$r\partial {\sigma }_r/\partial r={\sigma }_{\theta }-{\sigma }_r$ (7)

The geometric equation of this axisymmetric problem can be denoted as

$\begin{array}{l}{\epsilon }_r=\partial u/\partial r\\ {\epsilon }_{\theta }=u/r\end{array}\}$ (8)

By the geometric equation of Eq. (8), the compatibility equation can be shown as Eq. (9)

$\frac{\partial {\epsilon }_{\theta }}{\partial r}-\frac{{\epsilon }_r-{\epsilon }_{\theta }}{r}=0$ (9)

The difficulty of elastic-plastic problems is that the plastic deformation is nonlinear. By the theory of plastic, plastic strain is determined by the plastic potential, and rock mass materials generally obey the MC plastic potential function, which can be shown as Eq. (10)

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

where $\psi$ is the dilatancy angle of the rock mass, and $\alpha =(1+\sin \psi )/(1-\sin \psi )$.

The classical potential theory can be shown in the following equations, the radial and tangential plastic strains can be obtained:

$\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}$ (11)

Combined the integration results of Eq. (11) with the elastic strain, and considering the influence of roadway excavation, the total strain of the rock mass can be expressed as

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

where $G=E/\left(2\left(1+\mu \right)\right)$ is the shear modulus, $\mu$ is Poisson’s ratio, and $\lambda$ is the plastic strain multiplier. When $\lambda =0$, Eq. (12) represents the constitutive equation of the elastic rock mass.

The initial condition, contact condition and stress boundary condition of the present problem can be expressed as Eq. (13)

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

where $p_i$ is the support force, $R_c^{-}$ and $R_c^{+}$ represent the inner and outer boundaries, respectively, at $r=R_c$. When the radius tends to infinity, the radial stress equals the original in situ stress. Parameters ${\sigma }_{r0}^e$ and ${\sigma }_{\theta 0}^e$ are the radial stress and tangential stress, respectively, of rock mass in the elastic region at $t=0$.

3. Surrounding rock model

In the quiescent state, the Bingham viscoplastic model can keep nonzero deviatoric stresses when the stress is less than the yield stress, and no plastic flow appears. Once the stress reaches the critical value, plastic flow occurs [²⁴]. The generalized Bingham body is composed of Hookean, Newtonian, and St. Venant elements, so it is slightly different from the Bingham model. Elastic deformation occurs before applying stress is less than the yield stress, and the material manifests elasto-viscoplastic deformation when the yield stress is reached. Figure 3 shows the generalized Bingham model.

Fig. 3
Calculated mechanical model.

The strain of the generalized Bingham model is a function of time, the one-dimensional constitutive equation can expressed as Eq. (14)

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

where dot “∙” on the top of σ and ε denotes the derivative of time, $E$ is the elastic modulus, ${\sigma }_s$ is the yield strength of the slider, and $\eta$ is the viscosity coefficient of the dashpot.

With the aid of the correspondence law, convert E to 2G, convert η to 2η, convert ε to e _ij , convert σ to s _ij . It is generally considered that the volume deformation of rock mass is elastic, for the purpose of general description, the volume deformation is considered as elastic viscoplastic deformation as the deviator strain tensor. According to the above conversion method, one dimensional constitutive equation can be transformed into three-dimensional constitutive equation, so the three-dimensional equation 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}\}$ (15)

where $e_{ij}$ is the deviatoric strain tensor; $s_{ij}$ is the deviatoric stress tensor; ${\sigma }_{ij}$ is the stress tensor, ${\epsilon }_{ij}$ is the strain tensor, and both the stress and the strain tensors can be decomposed into the deviatoric stress and hydrostatic pressure, i.e., $s_{ij}={\sigma }_{ij}-{\sigma }_{kk}{\delta }_{ij}/3$ and $e_{ij}={\epsilon }_{ij}-{\epsilon }_{kk}{\delta }_{ij}/3$; ${\sigma }_{kk}$ and ${\epsilon }_{kk}$ represent the first invariants of the stress and strain, respectively; ${\delta }_{ij}$ is the Kronecker symbol, when $i=j$, ${\delta }_{ij}=1$, when $i\ne j$, ${\delta }_{ij}=0$; G is the shear modulus; K is the bulk modulus; $\eta$ is the viscosity coefficient of the deviatoric stress and deviatoric strain; ${\eta }^{\prime }$ is the viscosity coefficient of the volume strain (${\epsilon }_{kk}$) and hydrostatic pressure (${\sigma }_{kk}$).

In a circular opening, Eq. (15) can be unfolded in a polar coordinate system as Eq. (16)

$\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}\}$ (16)

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

By substituting the plane strain condition (${\epsilon }_z={\dot{\epsilon }}_z=0$) into the third equation of Eq. (16), the axial stress can be reformulated 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}\}$ (17)

The coefficients of the first equation and the second equation of Eq. (17) are equal. Therefore, η’ can be represented by η as follows

$\eta \text{'}=\frac{3K}{2G}\eta$ (18)

By combining Eq. (18) with Eq. (16), then using the equilibrium differential equation Eq. (7), $(3K-2G)/(6K+2G)$ can be replaced by Poisson’s ratio $\mu$, by simplification, the three-dimensional constitutive equation of this model in a polar coordinate system can be finally obtained as Eq. (19)

$\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}\}$ (19)

where ${\tilde{\sigma }}_r={\sigma }_r-{\sigma }_r^c$, ${\tilde{\sigma }}_{\theta }={\sigma }_{\theta }-{\sigma }_{\theta }^c$, and ${\sigma }_r^c$ and ${\sigma }_{\theta }^c$ are the radial stress and tangential stress, respectively, of the rock mass in the plastic region. Parameters ${\tilde{\sigma }}_r$ and ${\tilde{\sigma }}_{\theta }$ indicate that rheological deformation occurs when the stress exceeds the plastic yield limit.

4. Analytical solutions

4.1. Elasto-plastic solutions under the HB yield criterion

The elastic-plastic solutions in the surrounding rock under the MC criterion have been studied, and the results will not be listed here. In the HB rock mass, the plastic region and the elastic region are formed, and the analytical solution of the stress and displacement can be expressed:

(1) Plastic region ($a\le r\le R_c$)

By combining Eqs. (1), (5), (7) and (13), the stresses of the plastic region can be obtained as Eq. (20)

$\begin{array}{l}{\sigma }_r^c=A{\ln }^2\left(r/a\right)+B\ln \left(r/a\right)+p_i\\ {\sigma }_{\theta }^c=A{\ln }^2\left(r/a\right)+\left(2A+B\right)\ln \left(r/a\right)+B+p_i\end{array}\}$ (20)

where $A=m_c{\sigma }_s/4$, $B=\sqrt{m_c{\sigma }_s{p}_i+s_c{\sigma }_s^2}$.

By utilizing the geometric Eq. (8), the compatibility equation for deformation Eq. (9) and the total strain Eq. (12), the strains and displacements can be obtained as Eq. (21)

$\begin{array}{l}{\epsilon }_r^c=\left[d_1{\ln }^2\left(r/a\right)+\left(2d_1+d_2\right)\ln \left(r/a\right)+d_2+d_3-C\alpha r^{-1-\alpha }\right]/\left(2G\right)\hfill \\ {\epsilon }_{\theta }^c=\left(d_1{\ln }^2\left(r/a\right)+d_2\ln \left(r/a\right)+d_3+C{r}^{-\alpha -1}\right)/\left(2G\right)\hfill \\ u^c=r{\epsilon }_{\theta }^c\hfill \end{array}\}$ (21)

where $d_1=(1-2\mu )A$, $d_2=(1-2\mu )B+2dA$, $d_3=(1-2\mu )\left(p_i-p_0\right)+dB-2dA/\left(\alpha +1\right)$, $d=\left(1-\mu \right)\left(\alpha -1\right)/\left(\alpha +1\right)$, and $C$ is an integration constant.

(2) Elastic region ($R_c\le r<\infty$)

The stress, strain, and displacement of the elastic region can be easily obtained as Eqs. (22) and (23)

$\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}\}$ (22)

$\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}\}$ (23)

where $p_{}^c$ is the radial stress of elastic-plastic interface and $R_c$ is the radius of the interface. Because the radial stress and displacement at the elastic-plastic interface are continuous, the integration constant C, the radius, and the stress at the interface can be obtained as Eq. (24)

$R_c=a\exp \left[Q-\frac{2}{m_c{\sigma }_s}{\left(m_c{\sigma }_s{p}_i+s_c{\sigma }_s^2\right)}^{1/2}\right]$ (24)

where $Y=\frac{1}{2}{\left[{\left(\frac{m}{4}\right)}^2+m\frac{p_0}{{\sigma }_s}+s\right]}^{1/2}-\frac{m}{8}$, $Q=\frac{2}{m_c{\sigma }_s}{\left(m_c{\sigma }_s{p}_0+s_c{\sigma }_s^2-m_c{\sigma }_s^{2}Y\right)}^{1/2}$.

$C=R_c^{\alpha +1}\left[p_0-p^c-d_1{\ln }^2\frac{R_c}{a}-d_2\ln \frac{R_c}{a}-d_3\right]$ (25)

$p_{}^c={\sigma }_{r=R_c}=p_0-Y{\sigma }_s$ (26)

By combining Eq. (25) with Eq. (21), the solution of the displacement can be easily expressed as an explicit form. For elastic brittle plastic mechanical model, if the peak mechanical parameters are equal to the residual mechanical parameters, the above analytical solution can be converted to the Kastner solution. Therefore, the elasto-plastic analytical solution is a further elaboration of the Kastner method.

4.2. Elasto-viscoplastic solutions of the HB surrounding rock

Because only the cross section of the tunnel is studied, sequential excavation is not considered. It is assumed that the long tunnel is excavated at one time here. In the context of the generalized Bingham calculation model, when the tunnel is excavated, elastic deformation occurs instantaneously. Because the sudden stress release on the excavation surface, the stress redistributed in the surrounding rock need to achieve a new equilibrium state, and a viscoplastic region with radius R _t generates near the boundary of the tunnel.

4.2.1. Elastic solution of the rock mass at t=0

Because rheological test data indicate that elastic deformation is much larger than plastic deformation [²⁵]. Consequently, the mechanical behavior of rock mass can be assumed to be elastic at $t=0$. Namely, once the tunnel is excavated, only elastic deformation occurs. Hence, when $t=0$, this solutions can be represented by Eq. (22) and Eq. (23), $p_{}^c$ is replaced by $p_i$, and $R_c^{}$ by the excavation radius a.

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

where the superscript “e” and the subscript “0” represent the solutions of the elastic state and the time $t=0$, respectively.

4.2.2. The solutions of the initial plastic region ($a\le r\le {\rho }_0$)

The Newtonian element begins to play its role in the region with time when the stress reaches the yield condition. At $t=0^{+}$, the initial plastic region occurs. By substituting the elastic stress of Eq. (27) into the yield criteria of Eqs. (2) and (5), the radius of the initial plastic region can be expressed as follows

$4{\left(p_i-p_0\right)}^2{(a/{\rho }_0^{})}^4-m{\sigma }_s(p_i-p_0){(a/{\rho }_0^{})}^2=m{\sigma }_s{p}_0+s{\sigma }_s^2$ (28)

Eq. (28) is one variable equation of four times, and it can be calculated by the mathematical method.

Differentiation with respect to time in the geometric relation leads to

$r\frac{\partial {\dot{\epsilon }}_{\theta }}{\partial r}={\dot{\epsilon }}_r-{\dot{\epsilon }}_{\theta }$ (29)

Because the equilibrium equation of Eq. (7) should be fulfilled all the time, it leads to

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

Combining Eqs. (19), (29) and (30) leads to the differential equation as follows

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

where $\beta =G/\eta$.

Integrating with respect to r and t, Eq. (31) can be solved as Eq. (32)

${\tilde{\sigma }}_r+{\tilde{\sigma }}_{\theta }=f(t)+e^{-\beta t}g(r)$ (32)

where $f(t)$ is the function of the time t, and $g(r)$ is the function of the radius r, yet to be determined.

By means of the first and fourth equations of Eq. (13), and utilizing Eq. (30), the radial and tangential 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}\}$ (33)

It should be noted that the equations of Eq. (33) are valid within the initial plastic region, and it was assumed that the surrounding rock at radius r is always in a plastic state.

By combining Eq. (33) with the second equation of Eq. (19) and by integration of the tangential strain ${\dot{\epsilon }}_{\theta }$ with respect to t in Eq. (18), the displacement u of the region can be obtained as Eq. (34)

$\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}$ (34)

The unknown function $g_1(r)$ can be determined by the fourth equation in the initial condition of Eq. (13), namely,

$t=\mathrm{0,}u=u_0^e$

where the displacement $u_0^e$ at $t=0$ is defined in Eq. (27), the function g ₁(r) can be obtained as

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

Combining Eqs. (34) and (35), the displacement of the initial plastic region can be derived as follows

$\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)(fora\le r\le {\rho }_0)\end{array}$ (36)

4.2.3. The analytical solutions of the viscoplastic region

Outside the initial plastic region (i.e., ${\rho }_0^{}\le r\le R_t$), as a further assumption, the condition is introduced that the total stresses may fulfil the yield condition, but they may never exceed it [²⁰]. This hypothesis influences the development of stresses and displacements with time, but not the final state. At the moment of t, when the radius of the viscoplastic region extends to Rt, the stresses and displacements are described similar to the solutions of the elasto-plastic state. Namely, the stresses in the viscoplastic region are defined similarly to Eq. (19), and the displacements are similar to Eq. (21) and the radius Rt to Eq. (24). In these new equations, a and pi are replaced by ρ0 and P(t), respectively. The stresses, strains and displacements after transformation can be expressed as Eqs. (37), (38), (39) and (40)

$\begin{array}{l}{\sigma }_r^{vp}=A_{}^{vp}{\ln }^2(r/{\rho }_0^{})+B_{}^{vp}\ln (r/{\rho }_0^{})+P_{}(t)\hfill \\ {\sigma }_{\theta }^{vp}=A_{}^{vp}{\ln }^2(r/{\rho }_0^{})+(2A_{}^{vp}+B_{}^{vp})\ln (r/{\rho }_0^{})+B_{}^{vp}+P_{}(t)\hfill \end{array}\}$ (37)

$\begin{array}{l}\begin{array}{l}{\epsilon }_r^{vp}=[d_1^{vp}\ln \left(r/{\rho }_0^{}\right)\left(2+\ln \left(r/{\rho }_0^{}\right)\right)+d_2^{vp}\left(1+\ln \left(r/{\rho }_0^{}\right)\right)+d_3^{vp}\\ +r\left(\ln \left(r/{\rho }_0^{}\right)\left(1-2\mu \right)+d\right)m_c{\sigma }_s\left(\partial P\left(t\right)/\partial r\right)/\left(2\sqrt{m_c{\sigma }_{s}P\left(t\right)+s_c{\sigma }_s^2}\right)\\ +r\left(1-2\mu \right)\left(\partial P\left(t\right)/\partial r\right)-\alpha C^{vp}{r}^{(-{\alpha }_c-1)}]/(2G)\end{array}\hfill \\ {\epsilon }_{\theta }^{vp}=\left[d_1^{vp}{\ln }^2(r/{\rho }_0^{})+d_2^{vp}\ln (r/{\rho }_0^{})+d_3^{vp}+C^{vp}{r}^{-\alpha -1}\right]/(2G)\hfill \end{array}\}$ (38)

$u_{}^{vp}=\frac{r}{2G}\left[d_1^{vp}{\ln }^2(r/{\rho }_0^{})+d_2^{vp}\ln (r/{\rho }_0^{})+d_3^{vp}+C^{vp}{r}^{-\alpha -1}\right]$ (39)

$R_t={\rho }_0\exp \left[Q-\frac{2}{m_c{\sigma }_s}{\left(m_c{\sigma }_{s}P(t)+s_c{\sigma }_s^2\right)}^{1/2}\right]$ (40)

where $A^{vp}=m_c{\sigma }_s/4$, $B_{}^{vp}=\sqrt{m_c{\sigma }_{s}P(t)+s_c{\sigma }_s^2}$, $d_1^{vp}=(1-2\mu )A^{vp}$, $d_2^{vp}=(1-2\mu )B^{vp}+2d{A}^{vp}$, $d_3^{vp}=(1-2\mu )\left(P\left(t\right)-p_0\right)+d{B}^{vp}-2d{A}^{vp}/\left(\alpha +1\right)$, $P(t)={\sigma }_r^c(1-e^{-\beta t})+{\sigma }_{r0}^e{e}^{-\beta t}+[1-{\left(a/{\rho }_0^{}\right)}^2]f(t)$ and $C^{vp}=R_t^{\alpha +1}\left[p_0-p^c-d_1^{vp}{\ln }^2\frac{R_t}{{\rho }_0}-d_2^{vp}\ln \frac{R_t}{{\rho }_0}-d_3^{vp}\right]$.

Because the condition of displacement continuity at the position of $r={\rho }_0^{}$, the function $f(t)$ in Eq. (33) and Eq. (36) can be determined, i.e., Eq. (36) is equal to Eq. (39) at the radius $r={\rho }_0^{}$, and an equation about the function $f(t)$ can be expressed as

$\begin{array}{l}(1-\mu ){\sigma }_{\theta 0}^e-\mu {\sigma }_{r0}^e+(2\mu -1)p_0-d_3^{vp}-C^{vp}{(}_{{\rho }_0^{}}^{)}\\ +[1-2\mu +{(}^a]f(t)+\beta [1-2\mu +{(}^a]{\displaystyle \int f(t)dt}=0\end{array}$ (41)

The boundary condition of $f(t)$ is as Eq. (42)

$t=0:f\left(t\right)={\displaystyle \int f\left(t\right)dt=0}$ (42)

Eq. (41) is a transcendental equation, which can be solved numerically by a simple computer program.

In fact, the deformation of rock mass is the process that the stress is adjusted continuously to reach balance after the tunnel excavation. In the viscoplastic region, the stress and deformation of the rock mass changes with time. When the time $t\to \infty$, the solutions of this paper and those gotten by the elasto-plastic method are consistent at the corresponding position, and at this time, even if the time increases, the viscoplastic region will not change again. Eventually, the initial plastic region and the viscoplastic region becomes the plastic region, and the solution of the calculated model corresponds to the foregoing section 4.1.

4.3 Calculation methodology

By inputting geometric and 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 Figure 4.

Fig. 4
Flowchart of calculation methodology

5. Example verification

In order to analyze the influence of dilatancy angle on stress and displacement in the plastic region, dilation angles of 0º and 30º are selected. The parameters of the HB criterion are converted into equivalent MC parameters. By doing so, the solutions of the two yield criteria can be compared under similar geological conditions. The specific conversion method [²⁶] is as follows: by equating the areas covered by the two envelopes of the HB and MC criteria, Eq. (43) and Eq. (44) are formed, which may be solved for the equivalent MC rock mass parameters when the HB constants are known [²⁷]. This results in the following equations for the angle of friction φ and cohesive strength c:

$\phi ={\sin }^{-1}\left[\frac{3Jm{\left(s+m{\sigma }_{3n}\right)}^{J-1}}{\left(1+J\right)\left(2+J\right)+3Jm{\left(s+m{\sigma }_{3n}\right)}^{J-1}}\right]$ (43)

$c=\frac{{\sigma }_s\left[\left(1+2J\right)s+\left(1-J\right)m{\sigma }_{3n}\right]{\left(s+m{\sigma }_{3n}\right)}^{J-1}}{\left(1+J\right)\left(2+J\right)\sqrt{1+\left[6Jm{\left(s+m{\sigma }_{3n}\right)}^{J-1}\right]/\left[\left(1+J\right)\left(2+J\right)\right]}}$ (44)

Eq. (43) and Eq. (44) are derived from the generalized HB criterion. In this paper, J represents empirical parameters, where J =0.5. ${\sigma }_{3n}={\sigma }_{3\max }/{\sigma }_s$, where ${\sigma }_s$ is the uniaxial compressive strength and ${\sigma }_{3\max }$ is the maximum confining pressure. In tunnel engineering, the relationship between ${\sigma }_{3\max }$ and the rock mass strength ${\sigma }_{cm}$ is as follows

$\frac{{\sigma }_{3\max }}{{\sigma }_{cm}}=0.47{\left[\frac{{\sigma }_{cm}}{\gamma H}\right]}^{-0.94}$ (45)

${\sigma }_{cm}={\sigma }_s\frac{\left[m+4s-J(m-8s)\right]{(m/4+s)}^{J-1}}{2(1+J)(2+J)}$ (46)

where, $\gamma$ is the bulk density of the rock mass and H is the buried depth of the rock mass.

To validate the analytical solutions, the typical input data used by Carranza-Torres [⁷] were used. After equivalent conversion, the mechanical parameters of the HB and MC criteria are listed in Table 1.

Table 1
Geometric and mechanical parameters

Scheme	Value
Radius of tunnel, a (m)	1
Initial stress, p 0 (MPa)	15
Young’s modulus, E (MPa)	5700
Poisson’s ratio, μ (-)	0.3
Viscosity coefficient, η (Days·MPa)	6000
Uniaxial compressive strength, σ s (MPa)	30
Dilation angle, 𝜓 (Deg.)	0, 30
MC rock mass
c (MPa)	1.4184
φ (Deg.)	31.64
cc (MPa)	1.1214
φc (Deg.)	26.00
HB rock mass
m (-)	1.7
s (-)	0.0039
mc (-)	0.85
sc (-)	0.0019

5.1. The solutions of functions f(t) and $\int f(t)dt$

The f(t) function in implicit equation (41) is solved by MATLAB software. Figure 5 (a) and Figure 5 (c) show the curves of f(t)/p ₀ with the support force, dilation angle and time under the HB and MC criteria. Figure 5 (b) and Figure 5 (d) show the curves of $E\int f(t)dt/(\eta p_0)$ with the support force, dilation angle and time under the HB and MC criteria. f(t) show a trend of first increasing and then decreasing under the two yield criteria; when $t\to \infty$, the value of f(t) eventually tends to zero. $\int f(t)dt$ always increase with time. In MC rock mass, when the support force increases from 0MPa to 0.75MPa, the peak value of f(t) decreases from 2.984 to 2.243, and the decrease percentage is 24.83%. However, when the supporting force increases from 0 MPa to 0.75 MPa in HB rock mass, the peak value of f(t) decreased from 3.458 to 2.286, the percentage of decrease was 33.89%. It can be seen that the f(t) function of HB rock mass is more affected by the support force. In any rock mass, the value of $\int f(t)dt$ decreases with the increase of support force. In MC and HB rock mass, when the support force is 0 MPa, the peak value of f(t) increases from 2.984 and 3.458 to 4.552 and 4.899, respectively, with the increase percentage of 52.55% and 41.67%. It can be seen that the f(t) function of MC rock mass is more affected by the dilatancy angle.

Fig. 5
The solutions of f(t)/p ₀ and $E\int f(t)dt/(\eta p_0)$ with time under different p _i and ψ in the HB and MC rock masses: (a) f(t)/p ₀ in HB rock mass, (b) $E\int f(t)dt/(\eta p_0)$ in HB rock mass, (c) f(t)/p ₀ in MC rock mass, (d) $E\int f(t)dt/(\eta p_0)$ in MC rock mass.

5.2. Analytical solutions of stress and displacement

Figure 6 shows the variations in radial stresses and tangential stresses with radii at different times under the HB criterion. As shown in Figure 6, with increasing time, under the HB yield criterion, At the outer radius of the initial plastic region, i.e. r=ρ ₀, the difference of tangential stress becomes smaller and smaller. Eventually, the difference becomes zero. The stresses are compared at varying times under the HB criterion. When r/a=1, the value of the tangential stress at t=2, t=6 and t=60 days are 19.94 MPa, 11.43 MPa and 1.69 MPa, respectively. At the free face of roadway, compared with that at t=2 days, the tangential stress at t= 60 days decreases by 91.52%. The radial and tangential stresses in the plastic region at any other time are always smaller than at t=0 under the HB criterion.

Under the HB criterion, the radius of the initial plastic region and the final plastic region are 1.2996 m and 2.4168 m, respectively, when there is no supporting force at the free face of the roadway. The viscoplastic region is the result of the plastic region expansion outwards with time. When t=2, t=6, t=60 days, successively, continuous deformation occurs in the viscoplastic region. When $t\to \infty$, the radius of the viscoplastic region remains unchanged, meaning that the rock mass reaches the equilibrium state. By numerical calculation, the critical pressure p ^c is 6.12 MPa. A. B and C represent the positions of the radius of the viscoplastic region in three time states, respectively. With the increase of time, the radius of viscoplastic region expands outward, and the difference of circumferential stress values on both sides of ρ ₀ is smaller and smaller. When t=60 days, the circumferential stress in viscoplastic region and initial plastic region is continuous.

Fig. 6
The changes in the radial and tangential stresses with time and radius under the HB criterion.

The dilation angle is an important factor affecting the displacement. Figure 7 shows the variation in the radial displacement of the tunnel with time and dimensionless r/a under the MC criterion and HB criterion when the dilation angle is 0º and 30º. Obviously, the radial displacement increases gradually with increasing time. Meanwhile, the displacement changes rapidly at the value of r/a from 1 to 3 and tends to be flat after r/a=3 when ψ=0º, which indicates that large and rapid deformation occurs near the excavation face of the tunnel. For the condition of ψ=30º, the displacement varies rapidly when r/a increases from 1 to 2. Compared with the displacement of ψ=0º, the displacement of the excavation face increases 3.44 times when t=60 days. The radial displacement increases with increasing dilatancy angle. The radial displacement obtained in HB rock mass is larger than that in MC rock mass.

Fig. 7
The radial displacement of the HB and MC criteria with radius changes at different times: (a) ψ=0º, and (b) ψ=30º

5.3. Displacement Verification of the HB yield Criterion

To prove the correctness of the results obtained in this paper, the verification results of the HB rock mass are visualized in Figure 8. After calculation, the state at t = 300 days is selected. Because the viscoplastic region does not change at this time, the final state becomes an elastic-plastic rock mass. Some examples for the HB rock mass are given (Brown [³], Sharan [⁴], Zhang [⁸], Lee and Pietruszczak [²⁸], and Park [²⁹]). Herein, these previous results are utilized to verify this new approach. The mechanical parameters of the HB rock mass are a=5m, p _i =5MPa, p ₀ =30MPa, E=5.5GPa, μ=0.25, σ _s =30MPa, m=1.7, s=0.0039, m_c =1.0, s_c =0, and ψ=30º. As Figure 8 indicates, the radial displacements calculated in this paper are fully consistent with the solution of Lee & Pietruszczk, Zhang, and Park when ψ=30º, while the displacement solutions of Sharan overestimates and Brown et al.’s approximation underestimates the displacements. The calculated results of the radial displacement at t = 300 days are compared with the solutions of Zhang [⁸], and the two curves overlap. In this paper, the radial displacement on the excavation surface is 19.868 mm for ψ = 30º, which is approximately 2.8 times of 7.194 mm for ψ = 0º. And the analytical solution in this paper is compared with the result obtained by Fritz, the results are basically the same.

Fig. 8
Comparison of the radial displacement solutions for the Hoek-Brown rock mass (^a denotes ψ=0º)

5.4. Sensitivity analysis

5.4.1 Influence of support force

The dimensionless final plastic radii (R _c /a) with different supporting forces (p _i ) are calculated and compared with the analytical results of Zhang [⁸], as shown in Table 2. As seen from the table, the radii of the plastic region calculated in this paper are consistent with Zhang [⁸]. With the increase in the supporting force, the dimensionless ρ ₀/a and R _c /a both decrease. This law shows that applying the support force can effectively reduce the extension of the plastic region into the surrounding rock of the tunnel.

Figure 9 shows the initial and final plastic radii with different supporting forces (p _i ) for the HB and MC rock masses. The relationship between the radius of the initial plastic region and the supporting force is linear, and the R-square is as high as 99.9%. Under the two yield criteria, the slope of the straight line is the same. However, the radii of the initial plastic region under the HB yield criterion are larger under the same supporting force. The radii of the final plastic region also decrease with increasing support force. The fitting curve in Figure 9 (b) shows an approximate linear relationship under the MC yield criterion; the R-square is 99.4%. Under the HB yield criterion, the relationship between the support force and radii of the final plastic region is nonlinear. When the support force is less than 0.4MPa, the radii of the final plastic region under the HB yield criterion are larger than that under the MC yield criterion. When the support force is greater than 0.4MPa, the opposite is true.

Table 2
Computed results of the dimensionless initial and final plastic radii with different p _i (data in parentheses are solutions by the method of Zhang [⁸])

Rock mass	pi (MPa)	ρ0/a	Rc /a
MC rock	0	1.2856	2.2622 (2.2622)
	0.15	1.2791	2.1724 (2.1724)
	0.3	1.2726	2.0913 (2.0913)
	0.75	1.2530	1.8879 (1.8879)
HB rock	0	1.2996	2.4168 (2.4168)
	0.15	1.2930	2.2266 (2.2266)
	0.3	1.2865	2.1066 (2.1066)
	0.75	1.2667	1.8720 (1.8720)

Fig. 9
Comparison of the initial (a) and final (b) plastic radii with different supporting forces (p _i) for the HB and MC rock masses.

5.4.2 Influence of Young’s modulus

In order to analyze the influence of Young's modulus on the radius of viscoplastic region and the radial stress at the outer radius of initial plastic region, Young's modulus is taken as 2, 5, 10, 20 and 40GPa respectively, the creep time is taken as 4, 10 and 60 days respectively. As shown in Figure 10 (a), in MC rock mass, the radius of viscoplastic region increases with the increase of Young's modulus. Under three different time states, the radius of viscoplastic region increases from 1.573m, 1.677m and 2.095m to 2.178m, 2.246m and 2.262m respectively with the increase of Young's modulus from 2GPa to 40GPa, and the percentages of increase are 38.46%, 33.93% and 7.97% respectively. Similarly, as shown in Figure 10 (b), in HB rock mass, the radius of viscoplastic region increases with the increase of Young's modulus. Under three different time states, the radius of viscoplastic region increases from 1.455m, 1.630m and 2.170m to 2.219m, 2.381m and 2.417m respectively with the increase of Young's modulus from 2GPa to 40GPa, and the percentages of increase are 52.51%, 46.07% and 11.38% respectively. The location of the elastic-plastic interface, that is, the radius of the viscoplastic region, increases with the increase of time and Young's modulus. The larger the young's modulus is, the faster the viscoplastic region develops and the earlier the rock mass reaches equilibrium. The radius of the viscoplastic region in HB rock mass under the same Young's modulus, the viscoplastic radius of HB rock mass develops faster than that of MC rock mass.

Fig. 10
Influence of Young's modulus on the radius of viscoplastic region for MC (a) and HB (b) rock mass

In MC rock mass, when Young's modulus increases from 2GPa to 40GPa, the radial stress at the outer radius of the initial plastic region decreases with time. Under three different time states, the radial stress decreases from 3.704MPa, 3.131MPa and 1.539MPa to 1.312MPa, 1.412MPa and 1.104MPa respectively, and the percentage of decrease was 64.58%, 54.90% and 28.27%, respectively, as shown in Figure 11 (a). But in the HB rock mass, the radial stress decreases from 4.780MPa, 3.602MPa and 1.356MPa to 1.223MPa, 0.851MPa and 0.780MPa respectively, and the percentage of decrease was 74.41%, 76.37% and 42.48%, respectively, as shown in Figure 11 (b). Through comparative analyses, the radial stress at the outer radius of the initial plastic region in MC and HB rock mass decreases with the increase of time and Young's modulus, but the decrease percentage of HB rock mass is larger than that of MC rock mass.

Fig. 11
Effect of Young's modulus on radial stress at the outer radius of initial plastic region for MC (a) and HB (b) rock mass

6. Conclusions

The paper presented elasto-viscoplastic solutions for the stress and displacement of a circular roadway in an elasto-brittle-plastic rock mass with the generalized Bingham model and HB failure criterion. Rock rheology is the result of the continuous accumulation of internal damage, which eventually leads to rock burst. The analytical solution is time dependent, it can obtain the stress and displacement at any time in the process of surrounding rock deformation, this is the biggest difference from the static solution. The HB yield criterion can be used to predict the failure of surrounding rock to a greater extent and provide more reliable guidance for engineering design.

According to the mechanical parameters of HB rock mass, those of MC are deduced. The analytical results obtained under the two yield criteria are compared. When ${\rho }_0^{}\le r\le R_t$, the rock mass exhibits viscoplastic properties, and the radius increases with time. When R _t =R _c , the rock mass behaves brittle-plastic. The unknown function f(t) and its integral (cf. Figure 5) can be calculated for an arbitrary time form through Eq. (41).

The radial and tangential stresses in the plastic region at the moment which $t>0$ are always smaller than at t=0 under the HB criterion. Comparing the results of this study with previous research results in the literature, the displacement under the HB rock mass shows very good agreement with the solution of Lee & Pietruszczk, Zhang and Park, which indicates the correctness of the new approach. In this paper, the influence of the dilation angle and the support force is studied. When the dilation angle becomes 30º, the displacement increases noticeably, and the displacement under the HB yield criterion is larger than that under the MC yield criterion at the same time. With the increase in the supporting force, the dimensionless ρ ₀/a and R _c /a both decrease. The radius of viscoplastic region increases with the increase of Young's modulus, and the viscoplastic radius in HB rock mass develops faster. With Young's modulus increases, the radial stress at the outer radius of the initial plastic region decreases with time, but the decrease percentage in HB rock mass is larger.

Acknowledgments

The authors would like to express appreciation of the support provided by the Fundamental Research Funds for the Central Universities (Grant number: 2018ZZCX04). The authors gratefully appreciate the editors and reviewers for their valuable comments, which have greatly improved this paper.

References

• ^[1]
  Sadjadi F, Khalkhali A B. “Geotechnical Challenges of Tehran Metro Line 7 (South Northern Route),” Civil Engineering Journal, vol. 4, no. 5, pp. 1117, 2018. DOI: https://doi.org/10.28991/cej-0309161
  » https://doi.org/10.28991/cej-0309161
• ^[2]
  J. Jaak, and K. Daemen, “Tunnel support loading caused by rock failure,” International Journal of Rock Mechanics and Mining Sciences & Geomechanics Abstracts, vol. 12, no. 12, pp. 174, 1975. DOI: https://doi.org/10.1016/0148-9062(76)91601-6
  » https://doi.org/10.1016/0148-9062(76)91601-6
• ^[3]
  E. T. Brown, M. Asce, J. W. Bray, B. Ladanyi, F. Asce, and E. Hoek, “Ground response curves for rock tunnels,” Journal of Geotechnical Engineering, vol. 109, no. 1, pp. 15-39, 1983. DOI: https://doi.org/10.1016/0148-9062(83)91081-1
  » https://doi.org/10.1016/0148-9062(83)91081-1
• ^[4]
  S. K. Sharan, “Elastic-brittle-plastic analysis of circular openings in Hoek-Brown media,” International Journal of Rock Mechanics and Mining Sciences, vol. 40, no. 6, pp. 817-824, 2003. DOI: https://doi.org/10.1016/s1365-1609(03)00040-6
  » https://doi.org/10.1016/s1365-1609(03)00040-6
• ^[5]
  S. K. Sharan, “Exact and approximate solutions for displacements around circular openings in elastic-brittleplastic Hoek-Brown rock,” International Journal of Rock Mechanics and Mining Sciences, vol. 42, no. 4, pp. 542-549, 2005. DOI: https://doi.org/10.1016/j.ijrmms.2005.03.019
  » https://doi.org/10.1016/j.ijrmms.2005.03.019
• ^[6]
  C. Carranza-Torres, and C. Fairhurst, “The elastoplastic response of underground excavations in rock masses that satisfy the Hoek-Brown failure criterion,” International Journal of Rock Mechanics and Mining Sciences, vol. 36, no. 6, pp. 777-809, 1999. DOI: https://doi.org/10.1016/s0148-9062(99)00047-9
  » https://doi.org/10.1016/s0148-9062(99)00047-9
• ^[7]
  C. Carranza-Torres, “Dimensionless graphical representation of the exact elasto-plastic solution of a circular tunnel in a Mohr-Coulomb material subject to uniform far-field stresses,” Rock Mechanics and Rock Engineering, vol. 36, no. 3, pp. 237-253, 2003. DOI: https://doi.org/10.1007/s00603-002-0048-7
  » https://doi.org/10.1007/s00603-002-0048-7
• ^[8]
  Q. Zhang, B. S. Jiang, X. S. Wu, H. Q. Zhang, and L. J. Han, “Elasto-plastic coupling analysis of circular openings in elasto-brittle-plastic rock mass,” Theoretical and Applied Fracture Mechanics, vol. 60, no. 1, pp. 60-67, 2012. DOI: https://doi.org/10.1016/j.tafmec.2012.06.008
  » https://doi.org/10.1016/j.tafmec.2012.06.008
• ^[9]
  S. L Wang, X. T. Yin, H. Tang, and X. R. Ge, “A new approach for analyzing circular tunnel in strain-softening rock masses,” International Journal of Rock Mechanics & Mining Sciences, vol. 47, no. 1, pp. 170-178, 2010. DOI: https://doi.org/10.1016/j.ijrmms.2009.02.011
  » https://doi.org/10.1016/j.ijrmms.2009.02.011
• ^[10]
  B. S. Jiang, L. Yang, and L. P. Shi, “Stress analysis of cracked surrounding rock based on Hoek-Brown criterion,” Chinese Journal of Solid Mechanics, vol. 32, no. S1, pp. 300-305, 2011. DOI: https://doi.org/10.19636/j.cnki.cjsm42-1250/o3.2011.s1.051
  » https://doi.org/10.19636/j.cnki.cjsm42-1250/o3.2011.s1.051
• ^[11]
  Wu L, Zhang X, Zhang Z, et al. “Displacement and Deformation of the First Tunnel Lining During the Second Tunnel Construction,” Civil Engineering Journal, vol. 5, no. 2, pp. 332, 2019. DOI: https://doi.org/10.28991/cej-2019-03091248
  » https://doi.org/10.28991/cej-2019-03091248
• ^[12]
  H. N. Wang, S. Utili, and M. J. Jiang, “An analytical approach for the sequential excavation of axisymmetric lined tunnels in viscoelastic rock,” International Journal of Rock Mechanics and Mining Sciences, vol. 68, pp. 85-106, 2014. DOI: https://doi.org/10.1016/j.ijrmms.2014.02.002
  » https://doi.org/10.1016/j.ijrmms.2014.02.002
• ^[13]
  P. Nomikos, R. Rahmannejad, and A. Sofianos, “Supported axisymmetric tunnels within linear viscoelastic Burgers rocks,” Rock Mechanics and Rock Engineering, vol. 44, no. 5, pp. 553-564, 2011. DOI: https://doi.org/10.1007/s00603-011-0159-0
  » https://doi.org/10.1007/s00603-011-0159-0
• ^[14]
  A. Fahimifar, F. M. Tehrani, A. Hedayat, and A. Vakilzadeh, “Analytical solution for the excavation of circular tunnels in a visco-elastic Burger’s material under hydrostatic stress field,” Tunnelling and Underground Space Technology, vol. 25, no. 4, pp. 297-304, 2010. DOI: https://doi.org/10.1016/j.tust.2010.01.002
  » https://doi.org/10.1016/j.tust.2010.01.002
• ^[15]
  D. Sterpi, and G. Gioda, “Visco-plastic behaviour around advancing tunnels in squeezing rock,” Rock Mechanics and Rock Engineering, vol. 42, no. 2, pp. 319-339, 2009. DOI: https://doi.org/10.1007/s00603-007-0137-8
  » https://doi.org/10.1007/s00603-007-0137-8
• ^[16]
  M. Panet, “Time-dependent deformations in underground works,” International Society for Rock Mechanics, Proceeding of the fourth international congress on rock mechanics, Rotterdam, pp. 279-289, 1978. DOI: https://doi.org/10.1016/0148-9062(81)90337-5
  » https://doi.org/10.1016/0148-9062(81)90337-5
• ^[17]
  J. Sulem, M. Panet, and A. Guenot, “An analytical solution for time-dependent displacements in a circular tunnel,” International Journal of Rock Mechanics & Mining Sciences & Geomechanics Abstracts, vol. 24, no. 3, pp. 155-164, 1987. DOI: https://doi.org/10.1016/0148-9062(87)92647-7
  » https://doi.org/10.1016/0148-9062(87)92647-7
• ^[18]
  M. E. Gurtin, and E. Sternberg, “On the linear theory of viscoelasticity,” Archive for Rational Mechanics and Analysis, vol. 11, no. 1, pp. 291-356, 1962. DOI: https://doi.org/10.1007/BF00253942
  » https://doi.org/10.1007/BF00253942
• ^[19]
  G. A. C. Graham, “The solution of mixed boundary value problems that involve time-dependent boundary regions for viscoelastic materials with one relaxation function,” Acta Mechanica, vol. 8, no. 3-4, pp. 188-204, 1969. DOI: https://doi.org/10.1007/bf01182260
  » https://doi.org/10.1007/bf01182260
• ^[20]
  P. Fritz, “An analytical solution for axisymmetric tunnel problems in elasto-viscoplastic media,” International Journal for Numerical & Analytical Methods in Geomechanics, vol. 8, no. 4, pp. 325-342, 1984. DOI: https://doi.org/10.1002/nag.1610080403
  » https://doi.org/10.1002/nag.1610080403
• ^[21]
  Wang Z, Chen X, Xue X, et al. “Mechanical Parameter Inversion in Sandstone Diversion Tunnel and Stability Analysis during Operation Period,” Civil Engineering Journal, vol. 5, no. 9, pp. 1917-1928, 2019. DOI: https://doi.org/10.28991/cej-2019-03091382
  » https://doi.org/10.28991/cej-2019-03091382
• ^[22]
  M. H. Zhao, Y. Lei, and B. H. Ma, “Determination of ultimate bearing capacity of rock-socketed pile based on HB strength criterion,” Shuili Xuebao/Journal of Hydraulic Engineering, vol. 42, no. 9, pp. 1058-1064+1074, 2011. DOI: https://doi.org/10.13243/j.cnki.slxb.2011.09.016
  » https://doi.org/10.13243/j.cnki.slxb.2011.09.016
• ^[23]
  Babak F., Mohammad H. G., “Towards Bayesian Quantification of Permeability in Micro-scale Porous Structures - The Database of Micro Networks,” HighTech and Innovation Journal, vol. 1, no. 4, pp. 148-160, 2020. DOI: https://doi.org/10.28991/hij-2020-01-04-02
  » https://doi.org/10.28991/hij-2020-01-04-02
• ^[24]
  B. Amadei, and W. Z. Savage, “An analytical solution for transient flow of Bingham viscoplastic materials in rock fractures,” International Journal of Rock Mechanics and Mining Sciences, vol. 38, no. 2, pp. 285-296, 2001. DOI: https://doi.org/10.1016/s1365-1609(00)00080-0
  » https://doi.org/10.1016/s1365-1609(00)00080-0
• ^[25]
  S. Q. Yang, P. Xu, and P. G. Ranjith, “Damage model of coal under creep and triaxial compression,” International Journal of Rock Mechanics & Mining Sciences, vol. 80, pp. 337-345, 2015. DOI: https://doi.org/10.1016/j.ijrmms.2015.10.006
  » https://doi.org/10.1016/j.ijrmms.2015.10.006
• ^[26]
  E. Hoek, C. Carranza-Torres, and B. Corkum, “Hoek-Brown failure criterion-2002 edition,” Proceedings of the North American Rock Mechanics Society Meeting, Toronto, 2002.
• ^[27]
  A. I. Sofianos, and P. P. Nomikos, “Equivalent Mohr-Coulomb and generalized Hoek-Brown strength parameters for supported axisymmetric tunnels in plastic or brittle rock,” International Journal of Rock Mechanics & Mining Sciences, vol. 43, no. 5, pp. 683-704, 2006. DOI: https://doi.org/10.1016/j.ijrmms.2005.11.006
  » https://doi.org/10.1016/j.ijrmms.2005.11.006
• ^[28]
  Y. K. Lee, and S. Pietruszczak, “A new numerical procedure for elasto-plastic analysis of a circular opening excavated in a strain-softening rock mass,” Tunnelling and Underground Space Technology, vol. 23, no. 5, pp. 588-599, 2008. DOI: https://doi.org/10.1016/j.tust.2007.11.002
  » https://doi.org/10.1016/j.tust.2007.11.002
• ^[29]
  Kyung-Ho Park, and Yong-Jin Kim, “Analytical solution for a circular opening in an elastic-brittle-plastic rock,” International Journal of Rock Mechanics & Mining Sciences, vol. 43, no. 4, pp. 616-622, 2006. DOI: https://doi.org/10.1016/j.ijrmms.2005.11.004
  » https://doi.org/10.1016/j.ijrmms.2005.11.004

Edited by

Editor:

Marcílio Alves.

Publication Dates

• Publication in this collection
  08 Mar 2021
• Date of issue
  2021

History

• Received
  10 Nov 2020
• Reviewed
  24 Dec 2020
• Accepted
  28 Dec 2020
• Published
  07 Jan 2021