Analytical solution of deep tunnels in a strain-hardening elasto-plastic rock mass

Abstract

Excavation of tunnels produces a redistribuition of stresses and induces deformations in the rock mass around the tunnel’s cross section. In the case of elasto-plastic behavior of rock mass, plastic zones may appear. It is important to quantify the influence of this zone on the overall response of the tunnel. In this paper, we deduce a fully analytical solution in terms of displacements and stresses around a circular deep tunnel. The aim here is not to replace a 3D numerical calculation. This kind of analytical calculation are only useful to have a good understanding of the tunnel behavior in the preliminary phases of the project. For example, to perform parametric studies useful to choosing good parameters to introduce in a 3D numerical calculation. A homogeneous and isotropic rock mass is considered. For elasto-plastic behavior, the Tresca’s constitutive model with associate flow rule and Mohr-Coulomb’s constitutive model with non-associate flow rule are considered. For both, the idealized stress-strain curve presents a linear istropic hardening law. A geostatic-hydrostatic state of initial stresses and infinitesimal strains is assumed. The analytical solutions are compared with the FEM solutions demonstrating excellent agreement.

Keywords:
Deep tunnels; analytical solution; constitutive model; elastoplasticity

Graphical Abstract

1 INTRODUCTION

During the analysis of a deep tunnel it is essential to understand the behavior of the rock mass after the excavation. The excavation induces a disturbance zone around the tunnel that depends on the properties of the rock mass, the excavation method, in situ stress field, the geometry of the cross section of the tunnel and the type and stiffness of the support and the timing of its installation. This plastic zone is defined as a region in the rock mass where its stiffness parameters have changed due to the process of excavation and support of the tunnel. One of the tools commonly used to understand the interdependence of these factors is analytical solutions.

The analytical solutions consist of the use of equations deduced from the assumptions of continuous mechanics. However, for the closed-form solution, it is necessary to adopt some simplifying hypotheses, such as circular section, infinite homogeneous medium, axissimetry boundary conditions and geostatic-hydrostatic stress field. In any case, these solutions have the potential of providing a quick overview of the problem and are very useful to identify the relationship between the most important variables. The aim here is not to replace a three-dimensional numerical calculation. This kind of analytical calculation are only useful to have a good understanding of the tunnel behavior in the preliminary phases of the project. In addition, as they are quick calculations, they allow parametric studies that can be useful when choosing the parameters of a three-dimensional numerical calculation that takes more time.

The first analytical solution employed in deep tunnels derives from theoretical development to determine the field of stresses and deformations around openings in elastic media. Lamé (1866) proposed the first solution (in plane state of deformations) for cylindrical openings in linear elastic medium submitted to an initial state of hydrostatic geostatic stresses. Afterwards, Kirsch (1898) proposed a solution in elastic media generalizing the initial stresses field. In the following century, along with the development of the phenomenological theory of plasticity, several analytical solutions were proposed. In general, considering different yield criteria (Mohr-Coulomb, Hoek-Brown, Tresca for example), stress-strain behavior (elasto-plastic perfect, hardening/softening or brittle) and strategies for considering volumetric reduction during plastic deformation (such as using a non-associated flow rule). Studies such as Fenner (1938), Morrison and Coates (1955), Salencon (1969), Peck (1969), Daemen and Fairhurst (1970), Panet (1976), Berest et al. (1978), Detournay and Fairhurst (1987), Corbetta (1990) are exemples of this contribution. Brown et al. (1983) present a good review of these solutions until 1980s and presents its own analytical solution.

In recent analytical solutions, part of the researchers is studying the influence of a non-hidrostatic initial stress field considering the stress along the axis of the tunnel the major principal. Studies such as Reed (1988), Pan and Brown (1996), Lu et al. (2010), Zhou and Li (2011), Wang et al. (2012), including finite strain, such Vrakas and Anagnostou (2014) are examples of this contribution.

This article presents two analytical solution: the first considers a Tresca’s constitutive model with linear hardening behavior with associated flow rule. The second case, more complex, considers a Mohr-Coulomb’s constitutive model with linear hardening behavior with non-associated flow rule, in wich post-peak dilatancy occurs at a constant rate with major principal strain. In both cases, an initial geostatic-hydrostatic stress field and infinitesimal deformation is adopted.

2 DEFINITION AND HYPOTHESIS OF THE PROBLEM

The Figure 1 shows the geometry and boundary conditions of the problem. A circular deep tunnel of radius $R_i$ in a homogeneous and isotropic rock mass with an elastic-plastic hardening behavior. A plane state of deformation is assumed, i.e., the section is far from the three-dimensional effects of the excavation face, toghether with axial symmetry boundary conditions for uniform loading $P_{\infty }$. $P_i$ is the pressure on the wall of the tunnel that simulates the decompression caused by the advance of the excavation and placement of the lining. For example, it is used for generate the characteristic curve of rock mass in the Convergence-Confinement method (Panet et al., 2001).

Figure 1
Geometry and boundary conditions of the problem

When $P_i$ decreases from $P_{\infty }$ and the yield criterion is achived, a cylindrical plastic zone begins to develop around the tunnel. In an elastic-plastic rock mass with hardening, this zone can have two internal domains: a hardening elastic-plastic zone and a perfect plastic zone. In the absence of lining, when $P_i=0$, depending on the properties of rock mass, the advancement of this plastic zone can lead to high deformations, indicating the possible failure of the rock mass.

3 TRESCA CRITERION

3.1 Equations of the problem

Considering the axial symmetry condition of the problem in Figure 1, the resulting stress state at a distance $r$ is defined by the radial stress ${\sigma }_{rr}(r)$ and circumferential stress${\sigma }_{\theta \theta }(r)$, which are the minor and major principal stresses, ${\sigma }_3(r)$ and ${\sigma }_1(r)$ respectively, when the section of interest is far from the excavation face and with initial state of geostatic-hydrostatic stresses. The resulting displacement is defined by the radial displacement $u(r)$.

The strain tensor $\underset{\_}{\underset{\_}{\epsilon }}=\nabla \underset{\_}{u}$ in polar coordinates is

$\underset{\_}{\underset{\_}{\epsilon }}=\left(\begin{array}{ccc}\frac{\partial u}{\partial r}& 0& 0\\ 0& \frac{u}{r}& 0\\ 0& 0& 0\end{array}\right)$ (1)

and the stress tensor is written as

$\underset{\_}{\underset{\_}{\sigma }}=\left(\begin{array}{ccc}{\sigma }_{rr}& 0& 0\\ 0& {\sigma }_{\theta \theta }& 0\\ 0& 0& {\sigma }_{zz}\end{array}\right)$ (2)

Applying the condition $div(\underset{\_}{\underset{\_}{\sigma }})=0$, i. e., the problem is statically admissible, have the following differential equation:

${\sigma }_{\theta \theta }={\sigma }_{rr}+r\frac{\partial {\sigma }_{rr}}{\partial r}$ (3)

Considering the boundary conditions of Figure 1 which are

${\sigma }_{rr}(r=R_i)={\sigma }_i=-P_i$ (4)

$\underset{\_}{\underset{\_}{{\sigma }_0}}=\left(\begin{array}{ccc}{\sigma }_{\infty }& 0& 0\\ 0& {\sigma }_{\infty }& 0\\ 0& 0& {\sigma }_{\infty }\end{array}\right)where{\sigma }_{\infty }=-P_{\infty }$ (5)

In linear elasticity behaviour, where $E{\underset{\_}{\underset{\_}{\epsilon }}}^e=\left(1+\nu \right)(\underset{\_}{\underset{\_}{\sigma }}-\underset{\_}{\underset{\_}{{\sigma }_0}})-\nu tr(\underset{\_}{\underset{\_}{\sigma }}-\underset{\_}{\underset{\_}{{\sigma }_0}})\underset{\_}{\underset{\_}{1}}$ (Hooke’s law), the solution of this boundary value problem is

$\frac{u(r)}{r}=\frac{(1+\nu )}{E}(P_{\infty }-P_i){\left(\frac{R_i}{r}\right)}^2$ (6)

${\sigma }_{rr}(r)=(P_{\infty }-P_i){\left(\frac{R_i}{r}\right)}^2-P_{\infty }$ (7)

${\sigma }_{\theta \theta }(r)=-(P_{\infty }-P_i){\left(\frac{R_i}{r}\right)}^2-P_{\infty }$ (8)

${\sigma }_{zz}(r)=-P_{\infty }$ (9)

In infinitesimal plasticity the strain deformation $\underset{\_}{\underset{\_}{\epsilon }}$ is the sum of the tensor of elastic deformations ${\underset{\_}{\underset{\_}{\epsilon }}}^e$ and the tensor of the plastic deformation ${\underset{\_}{\underset{\_}{\epsilon }}}^p$, i.e.,

$\underset{\_}{\underset{\_}{\epsilon }}={\underset{\_}{\underset{\_}{\epsilon }}}^e+{\underset{\_}{\underset{\_}{\epsilon }}}^p$ (10)

This deformation ${\underset{\_}{\underset{\_}{\epsilon }}}^p$ will occur when the stresses reach the Tresca yield criterion written as

$F(\underset{\_}{\underset{\_}{\sigma }})={\sigma }_{rr}-{\sigma }_{\theta \theta }-2C(\alpha )$ (11)

It should be noted that the out-of-plane stress ${\sigma }_{zz}(r)$ is taken as the intermediate principal stress ${\sigma }_2(r)$. This assumption is valid with geostatic-hydrostatic in situ stress conditions and sections far from the excavation face. In this solution the plastic deformation was assumed to be associated, i.e.,$G=F$, and given by the following expression

${\underset{\_}{\underset{\_}{\epsilon }}}^p={\epsilon }^p\frac{\partial G}{\partial \underset{\_}{\underset{\_}{\sigma }}}={\epsilon }^p\frac{\partial F}{\partial \underset{\_}{\underset{\_}{\sigma }}}=\left(\begin{array}{ccc}{\epsilon }^p& 0& 0\\ 0& -{\epsilon }^p& 0\\ 0& 0& 0\end{array}\right)$ (12)

Here we are interested in the particular dependence of the cohesion $C(\alpha )$, shown in Figure 2, which introduces the linear hardening law.

Figure 2
Variation of the cohesion: A -linear hardening behavior and B - perfect plasticity behavior

And therefore,

$C(\alpha )=\{\begin{array}{c}{C}_0+C^{\prime }\alpha ,0\le \alpha \le {\epsilon }_0^p\\ C_1,\alpha \ge {\epsilon }_0^p\end{array}$ (13)

where

$C\text{'}=\frac{C_1-C_0}{{\epsilon }_0^p}$ (14)

and the internal variable $\alpha$ is the plastic strain magnitude ${\epsilon }^p$.

Introducing (12) in (10) and replacing ${\underset{\_}{\underset{\_}{\epsilon }}}^e$ in Hook's law gives

$E\left(\frac{\partial u}{\partial r}-{\epsilon }^p\right)={\sigma }_{rr}-\nu ({\sigma }_{\theta \theta }+{\sigma }_{zz})+(2\nu -1){\sigma }_{\infty }$ (15)

$E\left(\frac{u}{r}+{\epsilon }^p\right)={\sigma }_{\theta \theta }-\nu ({\sigma }_{rr}+{\sigma }_{zz})+(2\nu -1){\sigma }_{\infty }$ (16)

Making the addition of (15) and (16):

$E\frac{1}{r}\frac{\partial u}{\partial r}=(1-2\nu )(\nu +1)({\sigma }_{rr}+{\sigma }_{\theta \theta }-2{\sigma }_{\infty })$ (17)

So, using the (17) in (3) we obtain two equations valid in all the domain: (18) and (19) that we will use in the following deductions.

$E\frac{u}{r}-(1-2\nu )(\nu +1)({\sigma }_{rr}-{\sigma }_{\infty })=\frac{A}{r^2}$ (18)

$2{\epsilon }^p=r\frac{\partial }{\partial r}\left[\frac{(\nu +1)}{E}{\sigma }_{rr}+\frac{u}{r}\right]$ (19)

Replacing (18) in (19) the last one becomes:

${\epsilon }^p=\frac{1}{E^{\prime }}r\frac{\partial {\sigma }_{rr}}{\partial r}-\frac{A^{\prime }}{r}$ (20)

where

$A^{\prime }=\frac{A}{E}$ (21)

$E^{\prime }=\frac{E}{1-{\nu }^2}$ (22)

Considering a monotonic decreasing function in $P_i$ starting from $P_{\infty }$, we have the different behavior zones in the rock mass, show in Figure 3. A first phase, when $P_i$ remains close enough to the $P_{\infty }$, the entire rock mass remains in elasticity, having only zone 1. However, when $P_i$ is below a limit value, necessary for the stress state to reach the yield criterion at $r=R_i$, the plastic zone 2 with linear hardening behavior begins. When $P_i$ decreases further, a third zone, of perfect plasticity, may appear.

Figure 3
Elasto-plastic and elastic zones around the deep tunnel with hardening behaviour

3.2 Solution in the elastic zone 1 $(r\ge y)$

The solution in this zone is classic (Mandel, 1966). We define the radius $r=y$ those where the elastic zone begins. The conditions ${\epsilon }^p=0$ and $F=0$ in $r=y$ allows obtaining the value of the constant $A\text{'}$ using (18) and (19):

$A^{\prime }=-2C_0\frac{y^2}{E^{\prime }}$ (23)

Using the condition ${\epsilon }^p=0$ in the elastic domain and the boundary condition in $r=\infty$ we obtain the complete solution in this zone:

$\frac{u(r)}{r}=-\frac{C_0(1+\nu )}{E}{\left(\frac{y}{r}\right)}^2$ (24)

${\sigma }_{rr}(r)={\sigma }_{\infty }+C_0{\left(\frac{y}{r}\right)}^2$ (25)

${\sigma }_{\theta \theta }(r)={\sigma }_{\infty }-C_0{\left(\frac{y}{r}\right)}^2$ (26)

${\sigma }_{zz}(r)={\sigma }_{\infty }$ (27)

3.3 Solution in the plastic zone 3 $(R_i\le r\le x)$

We define the radius $r=x$ as the limit between the plastics zones 3 and 2. In this perfect plastic zone, the yield criterion is given by

$F(\underset{\_}{\underset{\_}{\sigma }})={\sigma }_{rr}-{\sigma }_{\theta \theta }-2C_1$ (28)

As far as the yield criterion is zero in this zone, the equilibrium equation (3) becomes

$\frac{1}{r}\frac{\partial {\sigma }_{rr}}{\partial r}=-2C_1$ (29)

The integration of the (19) together with the boundary condition of (4) leads to the following equation of the radial stress:

${\sigma }_{rr}(r)=2C_1\ln \left(\frac{R_i}{r}\right)+{\sigma }_i$ (30)

Replacing (30) into (18) we obtain the expression of the radial displacement given by

$\frac{u(r)}{r}=\frac{(1+\nu )}{E}\left\{(1-2\nu )\left[2C_1\ln \left(\frac{R_i}{r}\right)+{\sigma }_i-{\sigma }_{\infty }\right]-2C_0(1-\nu ){\left(\frac{y}{r}\right)}^2\right\}$ (31)

Afterwards, using the equilibrium equation and the plane state condition of deformations, the rest of the solution is given by

${\sigma }_{\theta \theta }(r)=2C_1\ln \left(\frac{R_i}{r}\right)+{\sigma }_i-2C_1$ (32)

${\sigma }_{zz}(r)=\nu ({\sigma }_{rr}+{\sigma }_{\theta \theta })+(1-2\nu ){\sigma }_{\infty }$ (33)

3.4 Solution in the plastic zone 2 $(x\le r\le y)$

In this zone of plasticity with hardening behaviour the yield criterion is written as

$F(\underset{\_}{\underset{\_}{\sigma }})={\sigma }_{rr}-{\sigma }_{\theta \theta }-2C(\alpha )$ (34)

As far as this criterion is null in this zone, the integration of the equilibrium equation with the condition in $r=x$ leads to

${\sigma }_{rr}(r)={\sigma }_x+\frac{2C_0}{1+2\frac{C^{\prime }}{E^{\prime }}}\left\{\frac{C^{\prime }}{E^{\prime }}\left[\frac{1}{r^2}-\frac{1}{x^2}\right]y^2+\ln \left(\frac{x}{r}\right)\right\}$ (35)

where ${\sigma }_x$ is the value of the radial stress in $r=x$. Due to the continuity of the radial stress between zone 3 and 2, equation (30) provides

${\sigma }_x={\sigma }_{rr}(r=x)=2C_1\ln \left(\frac{R_i}{x}\right)+{\sigma }_i$ (36)

In the same way that for the zone 3, the expression of the displacement is obtained using (18), which gives:

$\frac{u(r)}{r}=\frac{(1-2\nu )(\nu +1)}{E}\left\{{\sigma }_x+\frac{2C_0}{1+2\frac{C^{\prime }}{E^{\prime }}}\left[\frac{C^{\prime }}{E^{\prime }}\left(\frac{1}{r^2}-\frac{1}{x^2}\right)y^2+\ln \left(\frac{x}{r}\right)\right]-{\sigma }_{\infty }\right\}-\frac{2C_0}{E^{\prime }}{\left(\frac{y}{r}\right)}^2$ (37)

Using the equilibrium equation and the plane state condition of deformations, the rest of the solution is given by:

${\sigma }_{\theta \theta }(r)={\sigma }_x-\frac{2C_0}{1+2\frac{C^{\prime }}{E^{\prime }}}\left[\frac{C^{\prime }}{E^{\prime }}\left(\frac{1}{r^2}+\frac{1}{x^2}\right)y^2-\ln \left(\frac{x}{r}\right)+1\right]$ (38)

${\sigma }_{zz}(r)=\nu ({\sigma }_{rr}+{\sigma }_{\theta \theta })+(1-2\nu ){\sigma }_{\infty }$ (39)

3.5 Calculation of the plastic radius $x$and$y$

In order to finish the resolution of the problem we must calculate the plastic radius $x$ and $y$.In $r=x$ we have

${\epsilon }^p(r=x)={\epsilon }_0^{p}andC(\alpha )=C_1$ (40)

Substituting (40) in (20) we obtain

$x^2=\frac{2C_0}{E^{\prime }{\epsilon }_0^p+2C_1}{y}^2$ (41)

The value of $y$ is finally obtained writing the continuity of the radial stress in $r=y$, i.e.,

${\sigma }_{\infty }+C_0={\sigma }_x+\frac{2C_0}{1+2\frac{C^{\prime }}{E^{\prime }}}\left[\frac{C^{\prime }}{E^{\prime }}\left(\frac{1}{y^2}-\frac{1}{x^2}\right)y^2+\ln \left(\frac{x}{r}\right)\right]$ (42)

Together with (36) and (41) the value of $y$ is obtain in an explicit way:

$y^2=\frac{R_i}{{\omega }_0}{e}^{\frac{({\omega }_1+{\omega }_2+{\omega }_3)}{C_1}}$ (43)

with the following constants:

${\omega }_0=\frac{2C_0}{E^{\prime }{\epsilon }_0^p+2C_1}$ (44)

${\omega }_1={\sigma }_i-{\sigma }_{\infty }-C_0$ (45)

${\omega }_2=\frac{C_0}{1+2\frac{C^{\prime }}{E^{\prime }}}\ln ({\omega }_0)$ (46)

${\omega }_3=\frac{\frac{C^{\prime }}{E^{\prime }}\left[2(C_0-C_1)-{\epsilon }_0^p\right]}{1+2\frac{C^{\prime }}{E^{\prime }}}$ (47)

3.6 Comparison between the numerical and analytical solution

The numerical solution was obtained through the elastoplasticity algorithm implemented in the GEOMEC 91 finite element code. The local elastoplasticity integration algorithm of the stresses and internal variables and the global equilibrium iteration algorithm can be obtained in Bernaud (1991). The finite element mesh with 100 isoparametric elements with 9 nodes per element, together with their boundary conditions is shown in Figure 4.

Figure 4
Finite element mesh and boundary conditions

The characteristics of the following example are those of the Boom clay (Bernaud, 1991) obtained from samples at 240m deep: $E=1430MPa$, $\nu =0.4$, $C_0=0.21MPa$, $C_1=0.56MPa$, ${\epsilon }_0^p(\%)=2.4$, $P_{\infty }=4.5MPa$ and $P_i=4.5MPa\to 0$. Figure 5 illustrates the convergence curve, $U_i=-u(r=R_i)/R_i$, for a given application. The point marked in Figure 5 corresponds in FEM numerical calculation and we can see that numerical and analytical solutions are identical.

Figure 5
Convergence curve $P_{i}x{U}_i$for analytical and numerical solutions for the Tresca criterion

The Figure 6 shows the plastic radius for the internal pressure range.

Figure 6
Plastic radius $x$ and $y$ versus the internal pressure $P_i$ for the Tresca criterion

For a given pressure $P_i=2.5MPa$ the Figure 7 illustrates the stress $\underset{\_}{\underset{\_}{\sigma }}$ versus the radius $r$ for analytical and numerical solutions.

Figure 7
Stress versus the radius for analytical and numerical solutions for the Tresca criterion

As before, we can show the excellent agreement between the analytical and numerical solution. In Figure 7, through ${\sigma }_{\theta \theta }(r)$, the three zones that delimit the elastoplastic regime can be clearly seen.

4 MOHR-COULOMB CRITERION

4.1 Equations of the problem

In this case, the same strain and stresses tensor given by (1) and (2), equilibrium equation (3), the boundary conditions (4) and (5) and the superposition of the deformations elastic and plastic (10) are used. But the yield function is written as

$F(\underset{\_}{\underset{\_}{\sigma }})={\sigma }_{rr}-{\sigma }_{\theta \theta }+(K_p-1)({\sigma }_{rr}-H(\alpha ))$ (48)

where

$K_p=t{g}^2\left(\frac{\pi }{4}+\frac{\phi }{2}\right)$ (49)

$H(\alpha )=C(\alpha )\cot g(\phi )$ (50)

and $\phi$ is the friction angle. In this model the hardening behaviour was adopted only in the cohesive parameter $H$. Assuming a non-associated flow rule, i.e., $G\ne F$, where

$G(\underset{\_}{\underset{\_}{\sigma }})={\sigma }_{rr}-{\sigma }_{\theta \theta }+(\beta -1)({\sigma }_{rr}-H(\alpha ))$ (51)

and

$\beta =\frac{1+\sin ({\phi }_{\beta })}{1-\sin ({\phi }_{\beta })}$ (52)

obtained for plastic deformation

${\underset{\_}{\underset{\_}{\epsilon }}}^p={\epsilon }^p\frac{\partial G}{\partial \underset{\_}{\underset{\_}{\sigma }}}=\left(\begin{array}{ccc}\beta {\epsilon }^p& 0& 0\\ 0& -{\epsilon }^p& 0\\ 0& 0& 0\end{array}\right)$ (53)

where ${\phi }_{\beta }$ is the material dilatation angle. The internal variable associated with the isotropic hardening behavior is given by:

$\alpha =(1+\beta ){\epsilon }^p$ (54)

Here we are interested in the particular dependence of the $H(\alpha )$ in (55), shown in Figure 8, which introduces the linear hardening law.

Figure 8
Variation of the function$H(\alpha )$: A - linear hardening behaviour and B - perfect plasticity behavior

And therefore,

$H(\alpha )=\{\begin{array}{c}{H}_0+H^{\prime }{\epsilon }^p,0\le \alpha \le {\alpha }_0\\ H_1,\alpha >{\alpha }_0\end{array}$ (55)

where

$H^{\prime }=\frac{H_1-H_0}{{\epsilon }_0^p}$ (56)

${\alpha }_0=(1+\beta ){\epsilon }_0^p$ (57)

The equilibrium equation (3) and the Hooke’s law give the following relations (that are valid in all zones of behavior)

${\epsilon }^p=\frac{r}{1+\beta }\frac{\partial }{\partial r}\left[\frac{u}{r}+\frac{1+\nu }{E}{\sigma }_{rr}\right]$ (58)

$\left[\beta -\nu (1+\beta )\right]r\frac{\partial {\sigma }_{rr}}{\partial r}+(1+\beta )(1-2\nu )({\sigma }_{rr}-{\sigma }_{\infty })=\frac{E}{r^{\beta }(1+\nu )}\frac{\partial }{\partial r}\left(r^{\beta }u\right)$ (59)

Considering a monotonic decreasing function in $P_i$ starting from $P_{\infty }$, we have the different behavior zones in the rock mass. A first phase, when $P_i$ remains close enough to the $P_{\infty }$, the entire rock mass remains in elasticity, having only zone 1. However, when $P_i$ is below a limit value, necessary for the stress state to reach the yield criterion at $r=R_i$, the plastic zone 2 with linear hardening behavior begins. When $P_i$ decreases further, a third zone, of perfect plasticity, may appear.

In the following, we describe the deductions corresponding to the two cases of plasticity. For the first, we have only the plastic with linear hardening zone. The second is more complex because exist two plastics zones in the rock mass: the zone with perfect plasticity and the zone of plasticity with hardening behaviour.

4.2 First case: only one plastic zone

In this case exists only the plastic with hardening zone corresponding to the zone 2 illustrates in the Figure 9.

Figure 9
Plastic zone 2 and elastic zone 1 around the tunnel

4.2.1 Solution in the elastic zone 1 $(r\ge y)$

The solution in this zone is obtained writing the conditions:

${\epsilon }^p=0forr\ge y$ (60)

${\epsilon }^p=0andF=0whenr=y$ (61)

$\underset{\_}{\underset{\_}{{\sigma }_0}}=\left(\begin{array}{ccc}{\sigma }_{\infty }& 0& 0\\ 0& {\sigma }_{\infty }& 0\\ 0& 0& {\sigma }_{\infty }\end{array}\right)where{\sigma }_{\infty }=-P_{\infty }$ (62)

These conditions applied to the (58) and (59) gives the following solution:

${\sigma }_{rr}(r)={\sigma }_{\infty }-\frac{(K_p-1)}{(K_p+1)}({\sigma }_{\infty }-H_0){\left(\frac{y}{r}\right)}^2$ (63)

$\frac{u(r)}{r}=\frac{(1+\nu )}{E}\frac{(K_p-1)}{(K_p+1)}({\sigma }_{\infty }-H_0){\left(\frac{y}{r}\right)}^2$ (64)

4.2.2 Solution in the plastic zone 2 $(R_i\le r\le y)$

Considering the continuity of ${\sigma }_{rr}$ and $u$ in $r=y$:

${\sigma }_y={\sigma }_{rr}(r=y)={\sigma }_{\infty }-\frac{(K_p-1)}{(K_p+1)}({\sigma }_{\infty }-H_0)$ (65)

$u_y=u(r=y)=\frac{(1+\nu )}{E}\frac{(K_p-1)}{(K_p+1)}({\sigma }_{\infty }-H_0)y$ (66)

the displacements are obtained by integration of (59) between $x$ and $y$. Expressed as a function of ${\sigma }_{rr}$ we have

$\frac{u(r)}{r}=a{\sigma }_{rr}+b{\left(\frac{y}{r}\right)}^{\beta +1}+\frac{c}{r^{(\beta +1)}}{\displaystyle {\int }_y^r{r}^{\beta }{\sigma }_{rr}dr+d}$ (67)

with

$a=\frac{(1+\nu )}{E}\left[\beta -\nu (1+\beta )\right]$ (68)

$b=\frac{(1-{\nu }^2)}{E}\left[\frac{2(K_p-\beta ){\sigma }_{\infty }-(\beta +1)(K_p-1)H_0}{(K_p+1)}\right]$ (69)

$c=\frac{(1-{\nu }^2)}{E}(1-{\beta }^2)$ (70)

$d=\frac{(2\nu -1)(1+\nu )}{E}{\sigma }_{\infty }$ (71)

Due to the nullity of the yield criterion in this zone, the equilibrium equation writes

$r\frac{\partial {\sigma }_{rr}}{\partial r}=(K_p-1)({\sigma }_{rr}-H_0-H^{\prime }{\epsilon }^p)$ (72)

Due to the fundamental (58) we can eliminate ${\epsilon }^p$. So we obtain

$r\frac{\partial {\sigma }_{rr}}{\partial r}\left[1+\frac{(K_p-1)}{(1+\beta )}\frac{(1+\nu )}{E}{H}^{\prime }\right]-(K_p-1)({\sigma }_{rr}-H_0)=-\frac{(K_p-1)}{(1+\beta )}{H}^{\prime }r\frac{\partial }{\partial r}\left(\frac{u}{r}\right)$ (73)

Equation where the second member can be written as a function of ${\sigma }_{rr}$, given in the (63). So, an integral-differential equation in ${\sigma }_{rr}$ is expressed as

$a_1{r}^{\beta }\frac{\partial {\sigma }_{rr}}{\partial r}+a_2{r}^{(\beta -1)}{\sigma }_{rr}+\frac{a_3}{r^2}+\frac{a_4}{r^2}{\displaystyle {\int }_y^r{r}^{\beta }{\sigma }_{rr}}+a_5{r}^{\beta -1}=0$ (74)

$a_1=1+\frac{(1-{\nu }^2)}{E}(K_p-1)H^{\prime }$ (75)

$a_2=(K_p-1)\left[H^{\prime }\frac{(1-\beta )(1-{\nu }^2)}{E}-1\right]$ (76)

$a_3=(1-K_p)\frac{(1-{\nu }^2)}{E}{H}^{\prime }\left\{\frac{2(K_p-\beta ){\sigma }_{\infty }-(\beta +1)(K_p-1)H_0}{(K_p+1)}\right\}{y}^{(\beta +1)}$ (77)

$a_4=H^{\prime }\frac{(1-{\nu }^2)}{E}(K_p-1)(1-{\beta }^2)$ (78)

$a_5=\left(K_p-1\right)H_0$ (79)

The method of resolution (74) uses technics purely mathematics. Therefore, we expose in the following directly the solution obtained:

${\sigma }_{rr}(r)=k_1{c}_1{r}^{\left[k_1-(\beta +1)\right]}+k_2{c}_2{r}^{\left[k_2-(\beta +1)\right]}-(\beta +1)c_3$ (80)

with

$c_3=\frac{H_{0}E}{2H^{\prime }(1-{\beta }^2)(1-{\nu }^2)-E(\beta +1)}$ (81)

$k_1=\frac{-d_1+\sqrt{d_1^2-d_2}}{2}$ (82)

$k_2=\frac{-d_1-\sqrt{d_1^2-d_2}}{2}$ (83)

$d_1=\frac{-2\beta (K_p-1)(1-{\nu }^2)H^{\prime }-E(K_p+\beta )}{E+(1-{\nu }^2)(K_p-1)H^{\prime }}$ (84)

$d_2=\frac{4(K_p-1)(1-{\nu }^2)(1-{\beta }^2)H^{\prime }}{E+(1-{\nu }^2)(K_p-1)H^{\prime }}<0$ (85)

The constants $c_1$ and $c_2$ are obtained taken into account the boundary and continuity conditions of stresses, given by:

${\sigma }_{rr}(r=R_i)={\sigma }_i$ (86)

${\sigma }_{rr}(r=y)={\sigma }_{y}with{\sigma }_y={\sigma }_{\infty }-\frac{(K_p-1)}{(K_p+1)}({\sigma }_{\infty }-H_0)$ (87)

Therefore:

$c_1=\frac{{\sigma }_i-k_2{c}_2{R}_i^{\left[k_2-(\beta +1)\right]}+(\beta +1)c_3}{k_1{R}_i^{\left[k_1-(\beta +1)\right]}}$ (88)

$c_2=\frac{(\beta +1)c_3+{\sigma }_y-\left[{\sigma }_i+(\beta +1)c_3\right]{\left(\frac{y}{R_i}\right)}^{\left[k_1-(\beta +1)\right]}}{k_2{y}^{\left[k_2-(\beta +1)\right]}-k_2{R}_i{}^{(k_2-k_1)}{y}^{\left[k_1-(\beta +1)\right]}}$ (89)

Finally, the displacement is a function only of the variable $r$

$\frac{u(r)}{r}=\overline{A}{r}^{\left[k_1-(\beta +1)\right]}+\overline{B}{r}^{\left[k_2-(\beta +1)\right]}+\overline{C}{\left(\frac{y}{R_i}\right)}^{(\beta +1)}+\overline{D}{\left(\frac{y}{r}\right)}^{(\beta +1)}+\overline{E}$ (90)

with the follow constants:

$\overline{A}=\frac{c_1}{E}\left\{(1+\nu )\left[\beta -\nu (1+\beta )\right]k_1+(1-{\nu }^2)(1-{\beta }^2)\right\}$ (91)

$\overline{B}=\frac{c_2}{E}\left\{(1+\nu )\left[\beta -\nu (1+\beta )\right]k_2+(1-{\nu }^2)(1-{\beta }^2)\right\}$ (92)

$\overline{C}=\frac{(1-{\nu }^2)}{E}\left\{\frac{2(K_p-\beta ){\sigma }_{\infty }-(1+\beta )(K_p-1)H_0}{(K_p+1)}\right\}$ (93)

$\overline{D}=-\frac{(1-{\nu }^2)(1-{\beta }^2)}{E}\left\{\frac{2(K_p-\beta ){\sigma }_{\infty }-(1+\beta )(K_p-1)H_0}{(1-{\beta }^2)(K_p+1)}\right\}$ (94)

$\overline{E}=-\frac{(1+\nu )}{E}\left\{\left\{(\beta +1)\left[\beta -\nu (\beta +1)\right]+(1-\nu )(1-{\beta }^2)\right\}{C}_3-(2\nu -1){\sigma }_{\infty }\right\}$ (95)

4.2.3 Calculation of the plastic radius $y$

In order to complete the solution of the problem, we must calculate the plastic radius $y$. This is performed writing the nullity of ${\epsilon }^p$ and of the yield criterion in $r=y$ in the (58):

${r\frac{\partial }{\partial r}\left(\frac{u}{r}\right)|}_{r=y}+\frac{(1+\nu )}{E}(K_p-1)({\sigma }_y-H_0)=0$ (96)

The ${r\frac{\partial }{\partial r}\left(\frac{u}{r}\right)|}_{r=y}$ deduction can be performed using the (67) and (80) and obtaining an implicit equation of $y$:

$e_1{y}^{\left[k_2-(\beta +1)\right]}+e_2{y}^{(k_2-k_1)}+e_3=0$ (97)

with

$e_1=\frac{(1-{\nu }^2){(\beta +1)}^2(1-\beta )}{E{R}_i^{\left[k_1-(\beta +1)\right]}}\left(1-\frac{k_2}{k_1}\right)\left[{\sigma }_i+\frac{H_{0}E(\beta +1)}{2H^{\prime }(1-{\beta }^2)(1-{\nu }^2)-E(\beta +1)}\right]$ (98)

$\begin{array}{l}{e}_2=\frac{{\sigma }_y}{E}\left\{(1-{\nu }^2)(\beta +1)\left[k_2(K_p-\beta )-(1-{\beta }^2)\right]\right\}+\frac{4{\sigma }_{\infty }(1-{\nu }^2)({\beta }^2+\beta )k_2}{E(K_p+1)}+\\ H_0\left\{\frac{k_2(1+\nu )(K_p-1)}{E(K_p+1)}\left[2(1-\nu ){(\beta +1)}^2-(K_p+1)\right]+\frac{(1-{\nu }^2)(1-{\beta }^2)(k_2-1)}{2H^{\prime }(1-\beta )(1-{\nu }^2)-E}\right\}\end{array}$ (99)

$\begin{array}{l}{e}_3=k_2{R}_i^{(k_2-k_1)}\{\frac{{\sigma }_y}{E{k}_1}(1-{\nu }^2)(\beta +1)\left[1-{\beta }^2-k_1(K_p-\beta )\right]-\frac{4{\sigma }_{\infty }(1-{\nu }^2)({\beta }^2+\beta )}{E(K_p+1)}+\\ H_0\left[\frac{(1+\nu )(K_p-1)\left[K_p+1-2(1-\nu ){(\beta +1)}^2\right]}{E(K_p+1)}\right]+\frac{(1-{\nu }^2)(1-{\beta }^2)(\beta +1-k_1)}{k_1\left[2H^{\prime }(1-\beta )(1-{\nu }^2)-E\right]}\}\end{array}$ (100)

This solution is only valid as long as the stress values ${\sigma }_i$ maintain the condtion $\alpha \le {\alpha }_0$.

4.3 Second case: two plastics zones

In this case, the second zone (perfect plasticity) appears in $r=R_i$ (see Figure 10).

Figure 10
Plastic zones (2 and 3) and elastic zone 1 around the tunnel

4.3.1 Solution in the plastic zone 3 $(R_i\le r\le x)$

In this zone the yield criterion writes:

$F={\sigma }_{rr}-{\sigma }_{\theta \theta }+(K_p-1)({\sigma }_{rr}-H_1)$ (101)

Due to the nullity of $F$ in this zone, the equilibrium equation is:

$r\frac{\partial {\sigma }_{rr}}{\partial r}-(K_p-1){\sigma }_{rr}=(1-K_p)H_1$ (102)

The resolution of the differential equation above with the limit condition ${\sigma }_{rr}(r=R_i)={\sigma }_i$ gives the classical solution:

${\sigma }_{rr}(r)=({\sigma }_i-H_1){\left(\frac{r}{R_i}\right)}^{(K_p-1)}+H_1$ (103)

The calculation of the displacements is performed by the integration of the fundamental (59) between the radius $r=R_i$ and $r=x$:

$\begin{array}{l}\frac{u(r)}{r}=\frac{(1+\nu )}{E}\left\{\frac{(K_p-1)\left[\beta -\nu (1+\beta )\right]+(1+\beta )(1-2\nu )}{(\beta +K_p)}\right\}\frac{({\sigma }_i-H_1)}{R_i^{(K_p-1)}}\left[r^{(K_p-1)}-\frac{x^{(\beta +K_p)}}{r^{(\beta +1)}}\right]+\\ \frac{(1+\beta )}{E}(1-\nu -2{\nu }^2)(H_1-{\sigma }_{\infty })\left[1-{\left(\frac{x}{r}\right)}^{(\beta +1)}\right]+\frac{x^{\beta }}{r^{(\beta +1)}}{u}_x\end{array}$ (104)

Where $u_x$ representes the displacement in the point $r=x$. Due to the continuity of $u_x$ the displacements can be calculate by the equation of $u$ in the second plastic zone.

4.3.2 Solution in the plastic zone 2 $(x\le r\le y)$

In this zone, the stresses and displacements are written as the same way that in the first calculation, (67) and (80). The difference concerns the limit condition ${\sigma }_{rr}(R_i)={\sigma }_i$ of the first calculation that is replaced by the condition ${\sigma }_{rr}(x)={\sigma }_x$. So, the radial stress in this zone is written with the (80) but with others values of $C_1$ and $C_2$ that are:

${\sigma }_{rr}(r)=k_1{c}_1{r}^{\left[k_1-(\beta +1)\right]}+k_2{c}_2{r}^{\left[k_2-(\beta +1)\right]}-(\beta +1)c_3$ (105)

$c_1=\frac{{\sigma }_x-k_2{x}^{\left[k_2-(\beta +1)\right]}{c}_2+(\beta +1)c_3}{k_1{x}^{\left[k_1-(\beta +1)\right]}}$ (106)

$c_2=\frac{(\beta +1)c_3+{\sigma }_y-\left[{\sigma }_x+(\beta +1)c_3\right]{\left(\frac{y}{x}\right)}^{\left[k_1-(\beta +1)\right]}}{k_2{y}^{\left[k_2-(\beta +1)\right]}-k_2{x}^{(k_2-k_1)}{y}^{\left[k_1-(\beta +1)\right]}}$ (107)

All the others constants are the same as of the first calculation. We remember:

${\sigma }_x=({\sigma }_i-H_1){\left(\frac{x}{R_i}\right)}^{(K_p+1)}+H_1$ (108)

${\sigma }_y={\sigma }_{\infty }-\frac{(K_p-1)}{(K_p+1)}({\sigma }_{\infty }-H_0)$ (109)

The displacements in this zone are given by the (90) with the constants of (91) to (95). The constants $c_1$ and $c_2$ must be calculated with the (106) and (107). We remember this equation:

$\frac{u(r)}{r}=\overline{A}{r}^{\left[k_1-(\beta +1)\right]}+\overline{B}{r}^{\left[k_2-(\beta +1)\right]}+\overline{C}{\left(\frac{y}{r}\right)}^{(\beta +1)}+\overline{D}{\left(\frac{y}{r}\right)}^{(\beta +1)}+\overline{E}$ (110)

4.3.3 Calculation of the plastic radius $x$ and $y$

As before, we need to calculate the plastic radius $x$ and $y$ to complete the solution of the problem. With this goal we write the conditions about ${\epsilon }^p$ in $x$ and $y$:

${\epsilon }^p=\{\begin{array}{c}\mathrm{0,}r=y\\ {\epsilon }_0^p,r=x\end{array}andH(\alpha )=\{\begin{array}{c}{H}_0,r=y\\ H_1,r=x\end{array}$ (111)

Considering the (58) the first condition is

$y{\frac{\partial }{\partial r}\left(\frac{u}{r}\right)|}_{r=y}+\frac{(1+\nu )}{E}(K_p-1)({\sigma }_y-H_0)=0$ (112)

The derivation in function of $r$ for the (67) written in $r=y$ and the value of ${\sigma }_y$ allow to write the above equation as a function of $\gamma$ and $x$:

$\begin{array}{l}{q}_1+\frac{(1-{\nu }^2)}{E}{(}^{\beta }(1-\beta ){\gamma }^{\left[k_1-(\beta +1)\right]}(q_0{x}^{(K_p-1)}+q_2)\left[\frac{k_1-k_2}{k_1{k}_2(1-{\gamma }^{(k_1-k_2)})}\right]+\\ \frac{(1-{\nu }^2)}{E}{(}^{\beta }(1-\beta )q_3\left[\frac{1}{k_1({\gamma }^{(k_2-k_1)}-1)}-\frac{1}{k_2({\gamma }^{(k_2-k_1)}-1)}\right]=0\end{array}$ (113)

where $\gamma =y/x$ and the constants:

$q_0=\frac{{\sigma }_i-H_1}{R_i^{(K_p-1)}}$ (114)

$\begin{array}{l}{q}_1=H_0\left\{\frac{(1+\nu )(K_p-1)\left\{(1-\nu ){(1+\beta )}^2+2(1-K_p)\left[\beta -\nu (1+\beta )\right]\right\}}{E(K_p+1)}+\frac{(1-{\nu }^2)(1-{\beta }^2)}{2H^{\prime }(1-\beta )(1-{\nu }^2)-E}\right\}+\\ \frac{{\sigma }_{\infty }}{E}\left\{\frac{2(1+\nu )\left[\nu (1-{\beta }^2)+{\beta }^2\right]+(1-\nu )(1-{\beta }^2)}{(K_p+1)}-2(\beta +1)(1-{\nu }^2)(K_p-\beta )\right\}\end{array}$ (115)

$q_2=H_1+\frac{H_{0}E}{2H^{\prime }(1-\beta )(1-{\nu }^2)-E}$ (116)

$q_3=H_0\left[\frac{E}{2H^{\prime }(1-\beta )(1-{\nu }^2)-E}+\frac{K_p-1}{K_p+1}\right]+\frac{2{\sigma }_{\infty }}{(K_p+1)}$ (117)

The second condition (${\epsilon }^p={\epsilon }_0^p$ in $r=x$) gives:

${\epsilon }_0^p=x{\frac{\partial }{\partial r}\left(\frac{u}{r}\right)|}_{r=x}+\frac{(1+\nu )}{E}(K_p-1)({\sigma }_i-H_1){\left(\frac{x}{R_i}\right)}^{(K_p-1)}$ (118)

The derivation of the (67) is written in $r=x$ and we define $\gamma =y/x$ in the (118). Finally, we can express $x$ as a function of $\gamma$.

$x^{(K_p-1)}=\frac{q_4({\gamma }^{(k_2-k_1)}-1)+q_5({\gamma }^{(k_2-k_1)}-1){\gamma }^{(\beta +1)}-q_6{\gamma }^{(\beta +1-k_1)}+q_7}{q_8({\gamma }^{(k_2-k_1)}-1)-q_9}$ (119)

with the following constants:

$q_4=(1+\beta ){\epsilon }_0^p-(1-{\nu }^2)(1-{\beta }^2)\left(1-\frac{\beta +1}{k_1}\right)\left[\frac{H_1}{E}+\frac{H_0}{2H^{\prime }(1-\beta )(1-{\nu }^2)-E}\right]$ (120)

$q_5=\frac{2(1-{\nu }^2)(1+\beta )}{E}\left[\frac{2(K_p-\beta ){\sigma }_{\infty }-(1+\beta )(K_p-1)H_0}{(K_p+1)}\right]$ (121)

$q_6=\frac{(1-{\nu }^2){(1+\beta )}^2(1-\beta )}{E}\left(\frac{k_2}{k_1}-1\right)\left\{H_0\left[\frac{E}{2H^{\prime }(1-\beta )(1-{\nu }^2)-E}+\frac{(K_p-1)}{(K_p+1)}\right]+\frac{2{\sigma }_{\infty }}{(K_p+1)}\right\}$ (122)

$q_7=\frac{(1-{\nu }^2){(1+\beta )}^2(1-\beta )}{E{k}_1{k}_2}(k_2-k_1)\left[H_1+\frac{H_{0}E}{2H^{\prime }(1-\beta )(1-{\nu }^2)-E}\right]$ (123)

$q_8=\frac{{\sigma }_i-H_1}{R_i^{(K_p-1)}}\left[\frac{(1-{\nu }^2)(1-{\beta }^2)}{E}\left(1-\frac{\beta +1}{k_1}\right)+\frac{(1+\nu )}{E}(K_p-1)(1-\nu )(1+\beta )\right]$ (124)

$q_9=\frac{(1-{\nu }^2){(1+\beta )}^2(1-\beta )}{E}\left(\frac{k_2}{k_1}-1\right)\frac{{\sigma }_i-H_1}{R_i^{(K_p-1)}}$ (125)

So, the (113) and (119) are a system of two implicit equations for $x$ and $y$ that must be solve numerically. First we replace the expression of $y$ in function of $\gamma$ given by the (119) in the (113). With a numerical calculation (calculation of the zero of a function) we find the value of $\gamma =y/x$ for a given ${\sigma }_i$. This value introduced in the (119) gives the value of $x$ and consequently the value of $y$.

4.4 Comparison between the numerical and analytical solution

The numerical solution uses the same program as the the previous solution (GEOMEC 91) whose mesh and boundary conditions are shown in Figure 4. The characteristics of the following example are those of the Boom clay (Bernaud, 1991) obtained from samples at 240m deep: $E=1430MPa$, $\nu =0.4$, $\varphi =4^o$, ${\varphi }_{\beta }={2.07}^o$, $C_0=0.21MPa$, $C_1=0.56MPa$, ${\epsilon }_0^p(\%)=2.4$, $P_{\infty }=4.5MPa$ and $P_i=4.5MPa\to 0$. Figure 11 illustrates the convergence curve for a given application.

Figure 11
Convergence curve $P_{i}x{U}_i$for analytical solution for the Mohr-Coulomb criterion

The Figure 12 shows the plastic radius for spectrum of internal pressures.

Figure 12
Plastic radius $x$ and $y$ versus the internal pressure $P_i$ for the Mohr-Coulomb criterion

For a given pressure $P_i=1.0MPa$ (with $\nu =0.5$) Figure 13 illustrates the stress versus the radius r for analytical and numerical solutions.

Figure 13
Stress versus the radius for analytical and numerical solutions for the Mohr-Coulomb criterion

5 ANALYTICAL AND NUMERICAL CONVERGENCES - COMPARISON

Table 1 summarizes the result of analytical and numerical convergences $U_i=-u(r=R_i)/R_i$ of six different cases. Two considering the Tresca model (T) and five considering the Mohr-Coulomb model (MC). The characteristics of rock mass correspond to that of Boom’s clay (Bernaud, 1991).

Table 1
Characteristics mechanics and Convergences Ui (%) analytical and numerical

CASE	$F(\underset{\_}{\underset{\_}{\sigma }})$	$E$ (MPa)	$\nu$	$\phi$ ${}_{(°)}$	${\phi }_{\beta }$ $(°)$	$C_0$ (MPa)	$C_1$ (MPa)	${\epsilon }_0^p$ (%)	$P_i$ (MPa)	$P_{\infty }$ (MPa)	$U_i(\%)$ Analytical	$U_i(\%)$ Numerical
1	T	1430	0.4	-	-	0.21	0.56	2.4	2.5	4.5	5.87	5.91
2	T	1430	0.5	-	-	0.56	0.56	2.4	1.5	4.5	4.54	4.47
3	MC	1430	0.4	10	4	0.21	0.56	2.4	0.5	4.5	5.11	5.22
4	MC	1430	0.5	10	4	0.56	0.56	2.4	0.0	4.5	4.53	4.61
5	MC	1430	0.4	4	2.07	0.21	0.56	2.4	1.0	4.5	10.81	11.03
6	MC	1430	0.5	4	2.07	0.21	0.56	2.4	1.0	4.5	10.23	10.43

• The maximal difference between analytical and numerical results is 2% indicate the good agreement of the results.

6 CONCLUSIONS

In this paper, simple closed-form solutions of the displacements and stresses for an elastic-plastic rock mass of the circular deep tunnel with Tresca’s and Mohr-Coulomb’s linear hardening model was presented. Three zones appear around the tunnel: an elastic zone, a plastic hardening zone and a perfect plastic zone. The analytical solution for these three zones is shown. The analytical results were compared with the FEM results and a perfectly agreement was obtained.

References

• Berest, P., Minh, D. N., & Panet, M. (1978). Contribution a l’étude de la stabilité d’une cavité souterraine dans un milieu avec radoucissement. Revue Française de Géotechnique, (4), 63-72
• Bernaud, D. (1991). Tunnels profonds dans les milieux viscoplastiques: approches expérimentale et numérique (Doctoral dissertation, Ecole Nationale des Ponts et Chaussées).
• Brown, E. T., Bray, J. W., Ladanyi, B., & Hoek, E. (1983). Ground response curves for rock tunnels. Journal of Geotechnical Engineering, 109(1), 15-39.
• Corbetta, F. (1990). Calculs Analytiques et Numériques de Tunnels Profonds (Doctoral dissertation, Ecole Nationale Supérieure des Mines de Paris)
• Daemen, J. J. K., & Fairhurst, C. (1970). Influence of failed rock properties on tunnel stability. In The 12th US Symposium on Rock Mechanics (USRMS). American Rock Mechanics Association.
• Detournay, E., & Fairhurst, C. (1987). Two-dimensional elastoplastic analysis of a long, cylindrical cavity under non-hydrostatic loading. In International Journal of Rock Mechanics and Mining Sciences & Geomechanics Abstracts (Vol. 24, No. 4, pp. 197-211). Pergamon.
• Fenner, R. (1938). Untersuchungen zur erkenntnis des gebirgsdrucks. Glückauf, 74, 681-695 and 705-715.
• Kirsch, C. (1898). “Die theorie der elastizitat und die bedurfnisse der festigkeitslehre.” Zeitschrift des Vereines Deutscher Ingenieure, 42, 797-807.
• Lamé, G. (1866). Leçons sur la théorie mathématique de l’elasticité des corps solides par G. Lamé. Gauthier-Villars.
• Lu, A. Z., Xu, G. S., Sun, F., & Sun, W. Q. (2010). Elasto-plastic analysis of a circular tunnel including the effect of the axial in situ stress. International Journal of Rock Mechanics and Mining Sciences, 47(1), 50-59.
• Mandel, J. (1966). Cours de mécanique des milieux continus. Cours de mécanique des milieu continus. Gauthier-Villars.
• Morrison, R. G. K., & Coates, D. F. (1955). Soil mechanics applied to rock failure in mines. The Canadian Mining and Metallurgical Bulletin, 48(523), 701-711.
• Pan, X. D., & Brown, E. T. (1996). Influence of axial stress and dilatancy on rock tunnel stability. Journal of Geotechnical Engineering, 122(2), 139-146.
• Panet, M. (1976). Analyse de la stabilité d'un tunnel creusé dans un massif rocheux en tenant compte du comportement après la rupture. Rock mechanics, 8(4), 209-223.
• Panet, M., Givet, P. D. C. O., Guilloux, A., Duc, J. L. D. G. N. M., Piraud, J., & Wong, H. T. S. D. H. (2001). The convergence-confinement method. AFTES-recommendations des Groupes de Travait.
• Peck, R. B. (1969). Deep excavations and tunneling in soft ground. Proc. 7th ICSMFE, 1969, 225-290.
• Reed, M. B. (1988). The influence of out‐of‐plane stress on a plane strain problem in rock mechanics. International Journal for Numerical and Analytical Methods in Geomechanics, 12(2), 173-181.
• Salencon, J. (1969). Contraction quasi-statique d’une cavite a symetrie spherique ou cylindrique dans un milieu elastoplastique. In Annales des ponts et chaussées (Vol. 4, pp. 231-236).
• Vrakas, A., & Anagnostou, G. (2014). A finite strain closed‐form solution for the elastoplastic ground response curve in tunnelling. International Journal for Numerical and Analytical Methods in Geomechanics, 38(11), 1131-1148.
• Wang, S., Wu, Z., Guo, M., & Ge, X. (2012). Theoretical solutions of a circular tunnel with the influence of axial in situ stress in elastic-brittle-plastic rock. Tunnelling and Underground Space Technology, 30, 155-168.
• Zhou, X. P., & Li, J. L. (2011). Hoek-Brown criterion applied to circular tunnel using elastoplasticity and in situ axial stress. Theoretical and applied fracture mechanics, 56(2), 95-103.

Edited by

Editor:

Pablo Andrés Muñoz Rojas.

Publication Dates

• Publication in this collection
  18 Sept 2020
• Date of issue
  2020

History

• Received
  12 Mar 2020
• Reviewed
  30 June 2020
• Accepted
  31 July 2020
• Published
  05 Aug 2020