Elastic-plastic analysis of circular tunnel based on unified strength theory and considering multiple plastic zones

Abstract

Considering that the internal friction angle of the surrounding rock is not a constant but a function of the mean normal stress, and the unified strength theory (UST) is used as the plastic condition of the tunnel surrounding rock, the plastic zone of the circular tunnel surrounding rock is divided into multiple annular regions. For different plastic zones, there are different yield conditions represented by different stress functions, the elastic-plastic analysis of circular tunnels is performed by using the piecewise linearization of nonlinear yield function, stress and displacement in the elastic-plastic zone and the radius of the plastic zone are obtained by combining the equilibrium equations. It is shown that the radius of the plastic zone increases, the radial stresses in the elastic-plastic zone and the circumferential stresses in the plastic zone become smaller and the circumferential stresses in the elastic zone becomes larger compared with that of only one plastic zone. By using single factor analysis method, the calculated values of the tunnel displacements and plastic radii by UST were compared with those obtained by the Mohr-Coulomb (M-C) criterion, Drucker-Prager (D-P) criterion, and Zienkiewicz-Pande (Z-P) criterion. The analysis shows that under the same conditions, the solution of UST is smaller than M-C criterion and D-P criterion, therefore, the selection of rock strength criterion has great influence on the calculation of rock mechanics and engineering, the proper application of UST will guarantee the safety of engineering practice and have more practical value.

Introduction

After tunnel excavation, the stresses in the surrounding rock are redistributed, and in order to guarantee the stability of the tunnel, an elastic-plastic analysis must be performed to provide a basis for the tunnel support design. At present, the elastic-plastic analysis of the tunnel surrounding rock adopts M-C criterion^1–3 or H-B criterion,^4–6 etc. However, neither the M-C criterion nor the H-B criterion considers the effect of the intermediate principal stress. In fact, the rock mass strength is often observed to be dependent on the intermediate principal stress, and this effect is associated with the rock mass properties and stress states. ⁷ The UST, which considers all the stress components and their effects on the material yield and failure, can quantify the effect of the intermediate principal stress of various materials with different stress states and take full advantage of the material strength. ⁸

In routine engineering application, there are many examples of elastic-plastic analysis of surrounding rock based on the assumption that the internal friction angle of the rock is constant and to be governed by the UST. Tu et al. ⁹ investigated the radial subgrade modulus in an elastic-brittle-plastic rock mass, and modified analytical solutions were derived that are compatible with a generalized nonlinear unified strength criterion. Based on the unified strength theory and a non-associated linear flow rule, Zeng et al. ¹⁰ given an elastoplastic unified solution for displacements around a deep circular tunnel which taken the true elastic strains in the plastic zone into account. Sun et al. ¹¹ present the analytical solutions for the responses of tunnels excavated in rock masses exhibiting strain-softening behavior by introducing the unified strength theory to analyze the tunnel response. Zou et al. ¹² discussed the effect of the initial elastic displacement of the surrounding rock, stiffness of the support structure, and the coefficient b of the intermediate principal stress on the plastic zone. Hu and Yu ¹³ applied the unified strength theory to explore the strength properties of rocks under triaxial compressive loading and the elastic-plastic analysis of circular tunnel surrounds under hydrostatic pressure conditions. Zhang et al. ¹⁴ proposed the displacements accounting for three different definitions for elastic strains and different Young’s modulus in the plastic zone.

However, Matin and Chandler ¹⁵ proved that in the process of rock mass failure, the internal friction angle and cohesion were not constant values, and the cohesion strength and friction strength components could not be exerted simultaneously. Kennedy and Lindberg ¹⁶ showed the effect of yielding according to a nonlinear Mohr-Coulomb-type function in which the angle of internal friction decreases with increasing mean normal stress. Florence and Schwer ¹⁷ extended the linear Mohr–Coulomb criterion solution. Park et al. ¹⁸ discussed four different combinations of dilatancy angle and softening parameters to investigate the effects of elastic strain increments and variable dilatancy in the plastic region. Yang et al. ¹⁹ divided the plastic zone into multiple plastic zones based on the D-P criterion, and obtained a more accurate solution.

In this paper, the UST is used as the yield condition, and the elastic-plastic analysis of circular tunnels is performed by using the piecewise linearization of nonlinear yield function, the effects of variable internal friction angle within plastic region is considered. When the plastic zone is divided into seven circles, the solution of UST exhibits convergence and reaches stabilization, and this solution is compared with a plastic zone in which the internal friction angle of the rock remains unchanged. Lastly the UST solution is compared with M-C criterion solution, D-P criterion solution, and Z-P criterion solution, the accuracy and practical application of the proposed solution are illustrated by solving some examples.

Four yield criterion

Unified strength theory

If the compressive stresses are specified as positive and the tensile stresses are negative, the UST that which can appropriately reflect the effect of the intermediate principal stress can be written as

$$\begin{array}{c}\mathrm{when}{\sigma }_2\le \frac{\alpha {\sigma }_1+{\sigma }_3}{1+\alpha },F=\alpha {\sigma }_1-\frac{b{\sigma }_1+{\sigma }_3}{1+b}={\sigma }_t\hfill \\ \mathrm{when}{\sigma }_2\ge \frac{\alpha {\sigma }_1+{\sigma }_3}{1+\alpha },F\prime =\frac{\alpha }{1+b}({\sigma }_1+b{\sigma }_2)-{\sigma }_3={\sigma }_t\hfill \end{array}$$

(1)

where b is the UST parameter that reflects the influence of the intermediate principal stress and also a parameter of failure criterion, $\alpha$ is the tensile-compression ratio of the uniaxial strength of the material, ${\sigma }_t$ is the uniaxial tensile stress of the material. The UST can be written in terms of the cohesion $c$ and the internal friction angle $\phi$ as follows:

$$\begin{array}{c}\mathrm{when}{\sigma }_2\le \frac{{\sigma }_1+{\sigma }_3}{2}-\frac{{\sigma }_1-{\sigma }_3}{2}\sin \phi ,\hfill \\ F={\sigma }_1(1-\sin \phi )-\frac{b{\sigma }_2+{\sigma }_3}{1+b}(1+\sin \phi )=2c\cos \phi \hfill \\ \mathrm{and}\mathrm{when}{\sigma }_2\ge \frac{{\sigma }_1+{\sigma }_3}{2}-\frac{{\sigma }_1-{\sigma }_3}{2}\sin \phi ,\hfill \\ F\prime =\frac{{\sigma }_1+b{\sigma }_2}{1+b}(1-\sin \phi )-{\sigma }_3(1+\sin \phi )=2c\cos \phi \hfill \end{array}$$

(2)

For the strength theory, b has a range of [0, 1]. The UST degenerates into M-C criterion with b = 0, which defines the lower bound of the UST; the UST degenerates into the twin shear stress strength theory with b = 1, which defines the upper bound of the UST. When 0 < b < 1, a series of failure criteria between the two extremes can be obtained to describe the strength behaviors for various materials.

Z-P yield criterion

Zienkiewicz et al. have modified and generalized the M-C criterion. The Z-P condition used a quadratic yield curve on the meridian plane, the general expression of the hyperbolic yield function is given by

$$F=\beta p^2+{\alpha }_{1}p-k+{\left[\frac{q}{g\left({\theta }_{\sigma }\right)}\right]}^2=0$$

(3)

Where $p$ is the average value of the three principal stresses, $q$ is the generalized shear stress, $g\left({\theta }_{\sigma }\right)$ is the yield curve shape function of the $\pi$ plane, ${\alpha }_1$ and $\beta$ are coefficients, $k$ is the yield parameter.

where

$$\begin{array}{c}\tan \overline{\phi }=\frac{6\sin \phi }{3-\sin \phi },\overline{c}=\frac{6c\cos \phi }{3-\sin \phi }\hfill \\ \beta =-{\tan }^2\overline{\phi },{\alpha }_1=-2\overline{c}{\tan }^2\overline{\phi },k={\overline{c}}^2-a^2{\tan }^2\overline{\phi }\hfill \\ g\left({\theta }_{\sigma }\right)=\frac{2K}{(1+K)-(K-1)\sin 3{\theta }_{\sigma }}\hfill \\ (\sin 3{\theta }_{\sigma }=-\frac{3\sqrt{3}}{2}\frac{J_3}{\sqrt{J_2^3}},K=\frac{3-\sin \phi }{3+\sin \phi })\hfill \\ p=\frac{1}{3}({\sigma }_1+{\sigma }_2+{\sigma }_3)\hfill \\ q=\frac{1}{\sqrt{2}}{[{({\sigma }_1-{\sigma }_2)}^2+{({\sigma }_2-{\sigma }_3)}^2+{({\sigma }_1-{\sigma }_3)}^2]}^{\frac{1}{2}}\hfill \end{array}$$

(4)

Substituting ${\alpha }_1$ , $\beta$ , and $k$ into equation gives

$${\left[\frac{q}{g\left({\theta }_{\sigma }\right)}\right]}^2+a^2{\tan }^2\overline{\phi }={(p\tan \overline{\phi }+\overline{c})}^2$$

(5)

The Z-P criterion yield curve can be taken arbitrarily close to the M-C straight line with a = 0, equation can be rewritten as

$$\frac{q}{g\left({\theta }_{\sigma }\right)}=p\tan \overline{\phi }+\overline{c}$$

(6)

In engineering, the intermediate principal stress coefficient $n$ can be written in terms of the three principal stresses as follows

$$n=\frac{{\sigma }_2-{\sigma }_3}{{\sigma }_1-{\sigma }_3}$$

(7)

Equation can be written as ²⁰ ${\sigma }_2=n{\sigma }_1+(1-n){\sigma }_3$ , the $p$ , $q$ , $J_2$ , and $J_3$ is defined by ${\sigma }_1$ and ${\sigma }_3$ .

$$\begin{array}{c}p=\frac{1}{3}[(n+1){\sigma }_1+(2-n){\sigma }_3]\hfill \\ q=\sqrt{n^2-n+1}({\sigma }_1-{\sigma }_3)\hfill \\ J_2=\frac{1}{3}(n^2-n+1){({\sigma }_1-{\sigma }_3)}^2\hfill \\ J_3=\frac{1}{27}(2n^3-3n^2-3n+2){({\sigma }_1-{\sigma }_3)}^3\hfill \end{array}$$

(8)

Where $m=\frac{\sqrt{n^2-n+1}}{g\left({\theta }_{\sigma }\right)}$ , substituting equation (8) into equation (6), the Z-P criterion can be written in terms of $m,n,{\sigma }_1,{\sigma }_3,\overline{\phi }$ , and $\overline{c}$ as follows

$$[\frac{\tan \overline{\phi }}{3}(n+1)-m]{\sigma }_1+[\frac{\tan \overline{\phi }}{3}(2-n)+m]{\sigma }_3+\overline{c}=0$$

(9)

M-C yield criterion

In 1773, Cullen proposed that the failure of rock was mainly shear failure, the shear strength criterion in the plane is given by

$$\left|\tau \right|=c+\sigma \tan \phi$$

(10)

Where $\tau$ is the shear stress (shear strength) on the shear plane, $\sigma$ is the normal stress on the shear plane, $c$ is the cohesion, $\phi$ is the internal friction angle. $\sigma$ and $\tau$ can be written in terms of the first and third principal stresses as follows

$$\begin{array}{c}\sigma =\frac{1}{2}({\sigma }_1+{\sigma }_3)+\frac{1}{2}({\sigma }_1-{\sigma }_3)\cos 2\theta \hfill \\ \tau =\frac{1}{2}({\sigma }_1-{\sigma }_3)\sin 2\theta \hfill \end{array}$$

(11)

Substituting equation (11) into equation (10), the M-C criterion can be expressed by the principal stress as

$${\sigma }_1=\frac{1+\sin \phi }{1-\sin \phi }{\sigma }_3+\frac{2c\cos \phi }{1-\sin \phi }$$

(12)

D-P yield criterion

The D-P criterion with consideration of the influence of intermediate principal stress and hydrostatic pressure, can be written as

$$F=\alpha I_1+\sqrt{J_2}-k=0$$

(13)

Where $\alpha$ and $k$ are experimental constants, which are related to the internal friction angle $\phi$ and cohesion $c$ . Different stress conditions correspond to different values of $\alpha$ and $k$ . $\alpha$ and $k$ with consideration of the axial strain equal to zero, can be expressed as

$$\begin{array}{c}\alpha =\frac{\sin \phi }{\sqrt{3}\sqrt{3+{\sin }^2\phi }}\hfill \\ k=\frac{\sqrt{3}c\cos \phi }{\sqrt{3+{\sin }^2\phi }}\hfill \end{array}$$

(14)

Where $I_1$ is the first invariant of the stress tensor and $J_2$ is the second invariant of the deviatoric stress tensor, which can be written in terms of the three principal stresses ${\sigma }_1,{\sigma }_2$ , and ${\sigma }_3$ as follows

$$\begin{array}{c}{I}_1={\sigma }_1+{\sigma }_2+{\sigma }_3\hfill \\ J_2=\frac{1}{6}[{({\sigma }_1-{\sigma }_2)}^2+{({\sigma }_2-{\sigma }_3)}^2+{({\sigma }_1-{\sigma }_3)}^2]\hfill \end{array}$$

(15)

The relationship between the three principal stresses and the intermediate principal stress coefficient $n$ , appeared in equation (7), can be modified as ${\sigma }_2=n{\sigma }_1+(1-n){\sigma }_3$ , $I_1$ and $J_2$ are then expressed explicitly as

$$\begin{array}{c}{I}_1=\frac{3}{2}({\sigma }_1+{\sigma }_3)+(n-\frac{1}{2})({\sigma }_1-{\sigma }_3)\hfill \\ J_2=\frac{1}{3}(n^2-n+1){({\sigma }_1-{\sigma }_3)}^2=m^2{({\sigma }_1-{\sigma }_3)}^2\hfill \end{array}$$

(16)

Where $m=\sqrt{\frac{1}{3}(n^2-n+1)}$ , substituting equation (16) into equation (13), the D-P criterion can be written in terms of ${\sigma }_1,{\sigma }_3,\alpha ,k,m$ , and $n$ as follows

$$(m-n\alpha -\alpha ){\sigma }_1-(m-n\alpha +2\alpha ){\sigma }_3-k=0$$

(17)

Elastic-plastic analysis of tunnel surrounding rock based on UST

Basic assumptions and stress state analysis

At first the following basic assumptions are made. (1) The surrounding rock is homogeneous, isotropic, linear elastic, and no creep or viscous behavior. (2) The tunnel section is circular, in the infinite length of the tunnel the nature of the surrounding rock is consistent. (3) When the deep burial is greater than or equal to 20 times the radius of the tunnel, ignoring the self-weight of the rock within the influence of the tunnel, the error with the original problem does not exceed 5%, so the horizontal stress of primary rock can be simplified to be homogeneous. (4) The lateral pressure coefficient $\lambda =1$ , the tunnel is subjected to a uniform stress field which the stress of primary rock is equal. In this way, the original problem constitutes a plane strain circular hole problem which both the load and structure are axisymmetric (as Figure 1).

Figure 1.

Mechanics model of axisymmetric tunnel.

Where $R_k$ is the radius of the plastic zone, $r_a$ is the tunnel radius, $p_1$ is the support resistance, $p_0$ is the stress of primary rock; the principal stress ${\sigma }_{\theta }$ , the radial stress ${\sigma }_r$ , and axial stress ${\sigma }_z$ are the major, minor, and intermediate stresses (i.e. ${\sigma }_1={\sigma }_{\theta }$ , ${\sigma }_2={\sigma }_z$ , and ${\sigma }_3={\sigma }_r$ ) respectively, and are defined by the following relationships: ${\sigma }_z=\frac{1}{2}({\sigma }_r+{\sigma }_{\theta })$ , the condition (i.e. ${\sigma }_2>\frac{1}{2}({\sigma }_1+{\sigma }_3)-\frac{1}{2}({\sigma }_1-{\sigma }_3)\sin \phi$ ) is satisfied, substituting ${\sigma }_{\theta }$ and ${\sigma }_r$ into equation (2), the UST can be expressed by ${\sigma }_r$ , ${\sigma }_{\theta }$ , and ${\sigma }_z$ as

$$\frac{{\sigma }_{\theta }+b{\sigma }_z}{1+b}(1-\sin \phi )-{\sigma }_r(1+\sin \phi )=2c\cos \phi$$

(18)

Plastic zone stress

The stresses in the plastic zone satisfy the static equilibrium equation and also satisfy the plastic yield criterion. Notably, it is an efficient way to calculate the elastic-plastic solution of tunnel excavation by using piecewise linear yield function to represent the nonlinear yield criterion. Based on the influence of the normal stress on the internal friction angle of the rock, by dividing the plastic region with a number of thin annular rings. Figure 2 show the plastic zone is divided into $k$ concentric annular ring finite elements, and it is assumed that the radial stress ${\sigma }_{\mathrm{rp}}$ decreases uniformly along each ring, from the elastic-plastic interface ${\sigma }_{\mathrm{rp}}={\sigma }_R$ , at $r=R_k$ ; to the tunnel wall ${\sigma }_{\mathrm{rp}}=p_1$ , $r=r_a$ .

Figure 2.

Dividing plastic zone into k circles.

The equilibrium equation in the plastic region $R_1(r_a\le r\le R_1)$ is given by

$$r\frac{d{\sigma }_{\mathrm{rp}}}{\mathrm{dr}}+{\sigma }_{\mathrm{rp}}-{\sigma }_{\theta p}=0$$

(19)

Where r is the distance from any point on the surrounding rock of the tunnel section to the center of the tunnel circle. ${\sigma }_{\mathrm{rp}}$ and ${\sigma }_{\theta p}$ are the radial stresses and circumferential stresses of the plastic zone. $A_1$ and $B_1$ are defined by the following relationships:

$$\begin{array}{c}{A}_1=\frac{4(1+b)\sin {\phi }_1}{(2+b)(1-\sin {\phi }_1)}\hfill \\ B_1=\frac{4c(1+b)\sin {\phi }_1}{(2+b)(1-\sin {\phi }_1)}\hfill \end{array}$$

(20)

Equation may be modified as

$${\sigma }_{\mathrm{rp}}-{\sigma }_{\theta p}=-A_1{\sigma }_{\mathrm{rp}}-B_1$$

(21)

Substituting equation into equation yields

$$r\frac{d{\sigma }_{\mathrm{rp}}}{\mathrm{dr}}-A_1{\sigma }_{\mathrm{rp}}-B_1=0$$

(22)

The general solution can be obtained by solving equation

$${\sigma }_{\mathrm{rp}}=C_1{r}^{A_1}-\frac{B_1}{A_1}$$

(23)

Where $C_1$ is the integration constant, and it is determined by the boundary conditions (i.e. ${\sigma }_{\mathrm{rp}}{|}_{r=r_a}=p_1$ ), where $p_1$ is the support resistance, the integration constant $C_1$ can be expressed as

$$C_1=(p_1+\frac{B_1}{A_1}){r_a}^{-A_1}$$

(24)

Substituting equation into equations and (21), the stresses in the plastic zone $R_1(r_a\le r\le R_1)$ can be expressed as

$$\begin{array}{c}{\sigma }_{\mathrm{rp}}=(p_1+\frac{B_1}{A_1}){\left(\frac{r}{r_a}\right)}^{A_1}-\frac{B_1}{A_1}\hfill \\ {\sigma }_{\theta p}=(1+A_1)(p_1+\frac{B_1}{A_1}){\left(\frac{r}{r_a}\right)}^{A_1}-\frac{B_1}{A_1}\hfill \end{array}$$

(25)

Similarly, substituting the yield condition of plastic region $R_i$ into the equilibrium equation yields

$$r\frac{d{\sigma }_{\mathrm{rp}}}{\mathrm{dr}}-A_i{\sigma }_{\mathrm{rp}}-B_i=0$$

(26)

where

$$\begin{array}{c}{A}_i=\frac{4(1+b)\sin {\phi }_i}{(2+b)(1-\sin {\phi }_i)}\hfill \\ B_i=\frac{4c(1+b)\sin {\phi }_i}{(2+b)(1-\sin {\phi }_i)}\hfill \end{array}$$

(27)

Considering the boundary condition ${\sigma }_{\mathrm{rp}}{|}_{r=R_{i-1}}=p_i(2\le i\le k)$ , $p_i$ is the radial stresses at the interface between the $i$ th ring and the $i-1$ th ring, integrating equation (26), the stresses in the plastic region $R_i$ can be expressed as

$$\begin{array}{c}{\sigma }_{\mathrm{rp}}=(p_i+\frac{B_i}{A_i}){\left(\frac{r}{R_{i-1}}\right)}^{A_i}-\frac{B_i}{A_i}\hfill \\ {\sigma }_{\theta p}=(1+A_i)(p_i+\frac{B_i}{A_i}){\left(\frac{r}{R_{i-1}}\right)}^{A_i}-\frac{B_i}{A_i}\hfill \end{array}$$

(28)

For the plastic region $R_k$ , substituting the yield condition into the equilibrium equation yields

$$r\frac{d{\sigma }_{\mathrm{rp}}}{\mathrm{dr}}-A_k{\sigma }_{\mathrm{rp}}-B_k=0$$

(29)

Integrating equation (29), the stresses in the plastic region $R_k$ can be expressed as

$$\begin{array}{c}{\sigma }_{\mathrm{rp}}=C_k{r}^{A_k}-\frac{B_k}{A_k}\hfill \\ {\sigma }_{\theta p}=(1+A_k)C_k{r}^{A_k}-\frac{B_k}{A_k}\hfill \end{array}$$

(30)

the stresses in the elastic zone can be written as

$$\begin{array}{c}{\sigma }_{\mathrm{re}}=p_0+\frac{D}{r^2}\hfill \\ {\sigma }_{\theta e}=p_0-\frac{D}{r^2}\hfill \end{array}$$

(31)

By utilizing the boundary conditions ${\sigma }_{\mathrm{re}}={\sigma }_{\mathrm{rp}}$ and ${\sigma }_{\theta e}={\sigma }_{\theta p}$ at $r=R_k$ , combining equations (30) and (31), $C_k$ is then expressed explicitly as

$$C_k=\frac{2p_0{A}_k+2B_k}{(2+A_k)A_k{R}_k^{A_k}}$$

(32)

Substituting equation (32) into equation (30), the radial stress ${\sigma }_{\mathrm{rp}}$ and circumferential stress ${\sigma }_{\theta p}$ at the elastic-plastic interface, in which $r=R_k$ , is given by

$$\begin{array}{c}{\sigma }_R={\sigma }_{\mathrm{rp}}{|}_{r=R_k}=\frac{2p_0-B_k}{2+A_k}\hfill \\ {\sigma }_{\theta p}{|}_{r=R_k}=\frac{2(1+A_k)p_0+B_k}{2+A_k}\hfill \end{array}$$

(33)

Radius of plastic zone

The boundary radius $R_1$ , $R_2$ , …, $R_{k-1}$ , $R_k$ of each region can be determined by considering the continuity of the radial stresses at the interface between the $i$ th ring and the $i-1$ th ring (i.e. at $r=R_{i-1},{\sigma }_{\mathrm{rp}}=p_i$ ), utilizing equation by considering the radial stresses continuous conditions, the plastic radius $R_k$ is obtained as

$$\begin{array}{c}{R}_1=r_a{\left[\frac{p_2{A}_1+B_1}{p_1{A}_1+B_1}\right]}^{\frac{1}{A_1}}\hfill \\ R_2=R_1{\left[\frac{p_3{A}_2+B_2}{p_2{A}_2+B_2}\right]}^{\frac{1}{A_2}}\hfill \\ R_k=R_{k-1}{\left[\frac{{\sigma }_R{A}_k+B_k}{p_k{A}_k+B_k}\right]}^{\frac{1}{A_k}}\hfill \end{array}$$

(34)

Stress in the elastic zone

For the elastic zone of the tunnel surrounding rock, it can be considered as a thick-walled cylinder with the outer boundary approaches infinity and the inner boundary as $R_k$ . The internal boundary conditions are the same as that of the external boundary conditions. There is a stress of primary rock $p_0$ acting on the external surface, and there is a constraint reaction ${\sigma }_R$ acting on the inner interface. Applying the theory of axisymmetric plane strain circular hole problem, the stresses in the elastic zone can be expressed as

$$\begin{array}{c}{\sigma }_{\mathrm{re}}=p_0(1-\frac{{R_k}^2}{r^2})+{\sigma }_R\frac{{R_k}^2}{r^2}\hfill \\ {\sigma }_{\theta e}=p_0(1+\frac{{R_k}^2}{r^2})-{\sigma }_R\frac{{R_k}^2}{r^2}\hfill \end{array}$$

(35)

Displacement of elastic-plastic zone

In the plastic region $R_k$ , the compatibility equation and strain–displacement relations are given by

$$\frac{d{\epsilon }_{\theta }}{\mathrm{dr}}+\frac{{\epsilon }_{\theta }-{\epsilon }_r}{r}=0$$

(36)

$${\epsilon }_r=\frac{\mathrm{du}}{\mathrm{dr}},{\epsilon }_{\theta }=\frac{u}{r}$$

(37)

The law of plastic flow can be expressed as

$$d{\epsilon }_{\mathrm{ij}}^p=d\lambda \frac{\partial G}{\partial {\sigma }_{\mathrm{ij}}}$$

(38)

Where $d\lambda \ge 0$ is the proportional coefficient, $G$ is the plastic potential function.

If $G\not\equiv F$ , it is called the associated flow law, and if $G\not\equiv F$ , it is called the non-associated flow rule. ²¹

$$G=(1+A_k^{\left(p\right)}){\sigma }_{\mathrm{rp}}-{\sigma }_{\theta p}+B_k^{\left(p\right)}=0$$

(39)

$$F=(1+A_k){\sigma }_{\mathrm{rp}}-{\sigma }_{\theta p}+B_k=0$$

(40)

Where $A_k^{\left(p\right)}=\frac{4(1+b)\sin {\varphi }_k}{(2+b)(1-\sin {\varphi }_k)}$ , $B_k^{\left(p\right)}=\frac{4c(1+b)\sin {\varphi }_k}{(2+b)(1-\sin {\varphi }_k)}$ , ${\varphi }_k$ is the dilatancy angle of rock in region $k$ , which can be obtained from triaxial compression experiment.

Substitute equation (39) into equation (38), the plastic strain increments are obtained as

$$d{\epsilon }_{\theta }^p=-d\lambda ,d{\epsilon }_2^p=0,d{\epsilon }_r^p=(1+A_k^{\left(p\right)})d\lambda$$

(41)

Integrate equation yields

$${\epsilon }_{\theta }^p=-\lambda ,{\epsilon }_2^p=0,{\epsilon }_r^p=(1+A_k^{\left(p\right)})\lambda =-(1+A_k^{\left(p\right)}){\epsilon }_{\theta }^p$$

(42)

Plastic volumetric strain is given by

$${\epsilon }_v^p={\epsilon }_2^p+{\epsilon }_r^p+{\epsilon }_{\theta }^p=A_k^{\left(p\right)}\lambda$$

(43)

Considering the plane strain condition ${\epsilon }_2={\epsilon }_2^e+{\epsilon }_2^p=0$ (i.e. ${\epsilon }_2^e=0$ ), ${\sigma }_2$ can be given as

$${\sigma }_2=v({\sigma }_r+{\sigma }_{\theta })$$

(44)

the constitutive equation is given by

$$\begin{array}{c}{\epsilon }_r^e=\frac{1-v^2}{E}({\sigma }_r-\frac{v}{1-v}{\sigma }_{\theta })\hfill \\ {\epsilon }_{\theta }^e=\frac{1-v^2}{E}({\sigma }_{\theta }-\frac{v}{1-v}{\sigma }_r)\hfill \end{array}$$

(45)

the total strains can be decomposed into elastic and plastic parts as follows:

$$\begin{array}{c}{\epsilon }_{\theta }={\epsilon }_{\theta }^e+{\epsilon }_{\theta }^p\hfill \\ {\epsilon }_r={\epsilon }_r^e+{\epsilon }_r^p={\epsilon }_r^e-(1+A_k^{\left(p\right)}){\epsilon }_{\theta }^p\hfill \end{array}$$

(46)

Substituting equation (46) into equation (36), the differential equation for ${\epsilon }_{\theta }^p$ , can be written as

$$r\frac{d{\epsilon }_{\theta }^p}{\mathrm{dr}}+(2+A_k^{\left(p\right)}){\epsilon }_{\theta }^p={\epsilon }_r^e-{\epsilon }_{\theta }^e-r\frac{d{\epsilon }_{\theta }^e}{\mathrm{dr}}$$

(47)

Substituting the radial stresses and circumferential stresses in the plastic region $R_k$ (i.e. equation (30)) into equation (45), the elastic strains in the $k$ region can be written as

$${\epsilon }_r^e=\frac{1-v^2}{E}[\frac{(2v-1)B_k}{(1-v)A_k}+\left(\frac{1-2v-v{A}_k}{1-v}\right)(p_k+\frac{B_k}{A_k}){\left(\frac{r}{R_{k-1}}\right)}^{A_k}]$$

(48)

$${\epsilon }_{\theta }^e=\frac{1-v^2}{E}[\frac{(2v-1)B_k}{(1-v)A_k}+(1+A_k-\frac{v}{1-v})(p_k+\frac{B_k}{A_k}){\left(\frac{r}{R_{k-1}}\right)}^{A_k}]$$

(49)

Solving equation (47), ${\epsilon }_{\theta }^p$ is then expressed explicitly as

$${\epsilon }_{\theta }^p=\frac{(1-v^2)(-2A_k-A_k^2)}{E(A_k^{\left(p\right)}+A_k+2)}(p_k+\frac{B_k}{A_k}){\left(\frac{r}{R_{k-1}}\right)}^{A_k}+C_k{r}^{-(2+A_k^{\left(p\right)})}$$

(50)

Where $C_k$ is the integral constant which can be determined by using the elastic-plastic boundary condition ${\epsilon }_{\theta }^p{|}_{r=R_k}=0$ .

$$C_k=\frac{(1-v^2)(2A_k+A_k^2)}{E(A_k^{\left(p\right)}+A_k+2)}(p_k+\frac{B_k}{A_k}){\left(\frac{R_k}{R_{k-1}}\right)}^{A_k}{R_k}^{2+A_k^{\left(p\right)}}$$

(51)

Considering the boundary condition ${\epsilon }_{\theta }^k={\epsilon }_{\theta }^{k-1}$ at $r=R_{k-1}$ , ${\epsilon }_{\theta }^k$ and ${\epsilon }_{\theta }^{k-1}$ are the circumferential strains at the interface between in the $k$ th ring and the $k-1$ th ring, respectively. $C_{k-1}$ can be written as

$$\begin{array}{c}{C}_{k-1}={\epsilon }_{\theta }^k{R}_{k-1}^{2+A_{k-1}^{\left(p\right)}}+\frac{(1-v^2)(A_{k-1}^2+2A_{k-1})}{E(2+A_{k-1}+A_{k-1}^{\left(p\right)})}\hfill \\ (p_{k-1}+\frac{B_{k-1}}{A_{k-1}}){\left(\frac{R_{k-1}}{R_{k-2}}\right)}^{A_{k-1}}{R}_{k-1}^{2+A_{k-1}^{\left(p\right)}}\hfill \end{array}$$

(52)

Substituting equations (49)–(51) into equation (37), the displacement of plastic zone $R_k$ can be obtained by

$$u=r{\epsilon }_{\theta }=r({\epsilon }_{\theta }^e+{\epsilon }_{\theta }^p)$$

(53)

where

$$\begin{array}{c}{u}_k=r\frac{1-v^2}{E}\hfill \\ \left\{\begin{array}{c}(\frac{-A_k^2-2A_k}{A_k+A_k^{\left(p\right)}+2}+A_k+\frac{1-2v}{1-v})(p_k+\frac{B_k}{A_k}){\left(\frac{r}{R_{k-1}}\right)}^{A_k}+\frac{(2v-1)B_k}{(1-v)A_k}\hfill \\ \multicolumn{1}{c}{+\frac{A_k^2+2A_k}{A_k+A_k^{\left(p\right)}+2}(p_k+\frac{B_k}{A_k}){\left(\frac{R_k}{R_{k-1}}\right)}^{A_k}{\left(\frac{R_k}{r}\right)}^{2+A_k^{\left(p\right)}}}\end{array}\right\}\hfill \end{array}$$

(54)

Similarly, the strains and displacement in other regions can be obtained.

In the elastic region, substituting equation (35) into equation (45), the elastic strains can be written as

$${\epsilon }_{\theta }^e=\frac{1-v^2}{E}[(p_0-{\sigma }_R)(1+\frac{v}{1-v})\frac{R_k^2}{r^2}+(1-\frac{v}{1-v})p_0]$$

(55)

where

$${\epsilon }_{\theta }^p={\epsilon }_r^p={\epsilon }_2^p=0$$

(56)

The displacement of the elastic region (i.e. $u_e={\epsilon }_{\theta }r=({\epsilon }_{\theta }^p+{\epsilon }_{\theta }^e)r={\epsilon }_{\theta }^{e}r$ ), can be written as

$$u_e=r\frac{1-v^2}{E}[(p_0-{\sigma }_R)(1+\frac{v}{1-v})\frac{R_k^2}{r^2}+(1-\frac{v}{1-v})p_0]$$

(57)

Example and analysis

An ideal circular tunnel arranged horizontally is used as the calculation model, the radius of the tunnel $r_a$ is 3 m, the physical and mechanical parameters of the surrounding rock of $p_0=25$  MPa, $p_1=1$  MPa, $c=2$  MPa, $\phi =30°$ , $\nu =0.249$ , $G=5$  GPa, and $E=12.49$  GPa are adopted for the analysis. The intermediate principal stress coefficient $n$ tends to be 0.5 after the surrounding rock enters the plastic state. ²² At the same time, the linear approximation of the D-P criterion can be obtained by taking $b$ as 0.5 in the UST. Therefore, from the viewpoint of objectivity and accuracy, $n=0.5$ , $b=0.5$ .

UST solution for the elastic-plastic zone of the tunnel surrounding rock

According to the effect of normal stress on internal friction angle of rock, the plastic zone is divided into multiple circles.

When $p_1\le {\sigma }_{\mathrm{rp}}<0.4{\sigma }_R$ , it is ${\phi }_1=21°$ .

When $0.4{\sigma }_R\le {\sigma }_{\mathrm{rp}}<0.5{\sigma }_R,$ it is ${\phi }_2=22.5°$ .

When $0.5{\sigma }_R\le {\sigma }_{\mathrm{rp}}<0.6{\sigma }_R,$ it is ${\phi }_3=24°$ .

When $0.6{\sigma }_R\le {\sigma }_{\mathrm{rp}}<0.7{\sigma }_R,$ it is ${\phi }_4=25.5°$ .

When $0.7{\sigma }_R\le {\sigma }_{\mathrm{rp}}<0.8{\sigma }_R,$ it is ${\phi }_5=27°$ .

When $0.8{\sigma }_R\le {\sigma }_{\mathrm{rp}}<0.9{\sigma }_R,$ it is ${\phi }_6=28.5°$ .

When $0.9{\sigma }_R\le {\sigma }_{\mathrm{rp}}<{\sigma }_R,$ it is ${\phi }_7=30°$ .

When the surrounding rock is divided into multiple plastic zones and only one plastic zone (with an internal friction angle of $\phi =30°$ ), the distribution of the radial stress ${\sigma }_r$ , the circumferential stress ${\sigma }_{\theta }$ , the radial strains ${\epsilon }_r$ , and the circumferential strains ${\epsilon }_{\theta }$ are shown in Figures 3 and 4.

Figure 3.

Stresses distribution of surrounding rock divided into one plastic zone and multiple plastic zones.

Figure 4.

Strains distribution of surrounding rock divided into one plastic zone and multiple plastic zones.

From Figure 3, the circumferential stress reaches its maximum value at the elastic-plastic interface, whereas the radial stress increases monotonically along the radial direction. By dividing the plastic zone into multiple annular rings, the radius of the plastic zone increases. In the plastic zone, the radial and circumferential stresses decrease; in the elastic zone, the circumferential stress increases and the radial stress decrease, therefore, the effect of variable internal friction angle on the tunnel responses is remarkable. The stresses distribution curves almost coincide when the plastic zone is divided into six, seven, and eight circles. It can be seen that when the plastic zone is divided into seven circles, the stresses and radius of the plastic zone all converge to a solution of the problem.

From Figure 4, in comparison, the calculated values of the circumferential strain obtained by the proposed method is much larger than that obtained by the traditional method of calculation (the internal friction angle is constant). This once again proves the importance of considering the effect of the variable internal friction angle. So viewing from the accuracy and engineering safety, multiple plastic zones can be substitute for one plastic zone and applied to the practical engineering design with the piecewise linearization of nonlinear yield function method.

Comparison of UST with D-P criterion, Z-P criterion, and M-C criterion

In order to better analyze the radius and displacement of the plastic zone, the UST solution is compared with the D-P criterion solution, Z-P criterion solution, and M-C criterion solution. The effects of cohesion $c$ , internal friction angle $\phi$ , primary rock stress $p_0$ , and support resistance $p_1$ on the radius and displacement of the plastic zone based on the four different yield criterion are considered separately by using the single factor analysis. Referring to the relevant literature, the values of cohesion $c$ is range from $0.5$ to $7$  MPa, the values of internal friction angle $\phi$ is range from 15° to 60°, the values of primary rock stress $p_0$ is range from $10$ to $50$  MPa, and the values of support resistance $p_1$ is range from 0 to $2$  MPa.

Figure 5 shows the stresses distribution curves of the tunnel surrounding rock under the four yield criterion, from which it can be seen that the circumferential stresses in the elastic zone given by M-C criterion and D-P criterion are nearly equal, which are greater that of Z-P criterion and the UST, and the radial stresses curve of M-C criterion almost coincided with that of D-P criterion, which are smaller that of Z-P criterion and the UST. In comparison, the positions where the peak values are located are quite different, the circumferential stress at the elastic-plastic interface based on the M-C criterion is the largest and that D-P criterion is the smallest close to the UST, the radius of the plastic zone based on the D-P criterion is the largest and that Z-P criterion is the smallest.

Figure 5.

Stress distribution based on different yield criterion.

Figures 6 to 9 show that the radius and displacement of the plastic zone decrease with the increase of the cohesion $c$ and the internal friction angle $\phi$ based on four different yield criterions, and the change trends are similar. Under the same conditions, the calculated values of the tunnel displacements and plastic radii by UST were lower than those obtained by the M-C criterion and D-P criterion, and they had a similar variation trend when the cohesion $c\ge 4$  MPa and the internal friction angle $\phi \ge 35°$ , and the solutions between different criterions has large difference when $c<4$  MPa.

Figure 6.

Relationship between cohesion and plastic zone radius.

Figure 7.

Relationship between internal friction angle and plastic zone radius.

Figure 8.

Relationship between cohesion and plastic zone displacement.

Figure 9.

Relationship between internal friction angle and plastic zone displacement.

It can be seen from Figures 10 and 11 that the radius and displacement of the plastic zone increase with the increase of the stress $p_0$ of primary rock. The growth rate of the D-P criterion solution is the largest, followed by the M-C criterion solution, the UST solution is smaller, and the growth rate of the Z-P criterion solution is the smallest. There was little difference between the four yield criterion when the stress of primary rock $p_0$ is small, but with the increase of $p_0$ , there was a significant difference in the four yield criterion.

Figure 10.

Relationship between stress of primary rock and plastic zone radius.

Figure 11.

Relationship between stress of primary rock and plastic zone displacement.

It is clear from Figures 12 and 13 that the calculated values of the radius and displacement by the four yield criterion are approximately linear distributing along the support resistance, they all decrease with the increase of the support resistance, and had a similar variation trend. Under the same support resistance condition, the solution of D-P criterion is the largest, the solution of M-C criterion is the second, the solution of UST is smaller, and the solution of Z-P criterion is the smallest.

Figure 12.

Relationship between support resistance and plastic zone radius.

Figure 13.

Relationship between support resistance and plastic zone displacement.

It can be seen that $R_k$ and $u_p$ under the four criteria have similar trends with $c$ , $\phi$ , $p_0$ , and $p_1$ , respectively. In comparison, D-P criterion has the largest solution and Z-P criterion has the smallest solution. The M-C criterion completely ignores the influence of the intermediate principal stress on the deformation and failure of surrounding rock, while the UST considers both the influence of the intermediate principal stress and the influence of hydrostatic pressure. For elastoplastic analysis of tunnels, the radius of plastic zone based on UST is relatively small and the distribution is relatively uniform. Therefore, the application of UST criterion can improve the applicability and accuracy of the project. The conclusion obtained by comparison can be used for reference in engineering. It is of great significance to academics and engineering.

Conclusions

The surrounding rock suffers various degrees of loosening caused by tunnel excavation, and the internal friction angle is no longer a constant but a function of the mean normal stress. In this paper, based on the unified strength theory, the plastic zone of the surrounding rock is divided into $k$ circles, the elastic-plastic analysis of circular tunnels is performed by using the piecewise linearization of nonlinear yield function, and the following main conclusions can be drawn:

• (1) Compared with only one plastic zone, after dividing the plastic zone of the surrounding rock into seven annular rings with different internal friction angles, the solution of UST exhibits convergence and reaches stabilization. The variable internal friction angle showed a significant effect on the tunnel responses: the radius of the plastic zone increases, the radial stress of the elastic-plastic zone and the circumferential stress of the plastic zone decrease, and the circumferential stress of the elastic zone increases.

• (2) Both the radius and displacement of the plastic zone decrease with the increase of cohesion, internal friction angle, and support resistance of the rock mass, while increase with the increase of stress of primary rock. Under the same conditions, the solution of the D-P criterion is the largest, the solution of the M-C criterion is the second largest, the solution of the UST is smaller, and the solution of the Z-P criterion is the smallest.

• (3) The calculated values of the tunnel displacements and plastic radii by UST were lower than those obtained by the M-C criterion, indicating that the M-C criterion over-estimated the tunnel responses. The effect of intermediate principal stress helps prevent the surrounding rock from yielding and resisting tunnel displacement. The selection of rock strength criterion has great influence on the calculation of rock mechanics and engineering. The selection of strength criterion that conforms to the actual conditions of rock is the basis of the correct analysis of rock engineering mechanics.

• (4) Although the solution proposed in this paper is limited in scope, it appears to be useful for the preliminary design of circular rock tunnels to provide a basis for the tunnel support design. Generally, correctly considering the effects of variable internal friction angle within plastic region is shown to be important in elastic-plastic analysis of circular tunnel. By comparison and considering the effect of the variable internal friction angle, the UST can comprehensively reflect the various strength characteristics of rock, which is closer to the engineering practice.