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

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 geostatichydrostatic state of initial stresses and infinitesimal strains is assumed. The analytical solutions are compared with the FEM solutions demonstrating excellent agreement.


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 andLi (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.

DEFINITION AND HYPOTHESIS OF THE PROBLEM
The Figure 1 shows the geometry and boundary conditions of the problem. A circular deep tunnel of radius i R 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 ∞ . i P 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).  When i P decreases from P ∞ 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 0 i P = , 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.

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 ( ) rr r σ and circumferential stress ( ) r θθ σ , which are the minor and major principal stresses, 3 ( ) r σ and 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 u ε = ∇ in polar coordinates is and the stress tensor is written as Applying the condition div( ) 0 σ = , i. e., the problem is statically admissible, have the following differential equation: Considering the boundary conditions of Figure 1 which are ( ) In infinitesimal plasticity the strain deformation ε is the sum of the tensor of elastic deformations e ε and the tensor of the plastic deformation p ε , i.e., e p ε ε ε = + This deformation p ε will occur when the stresses reach the Tresca yield criterion written as It should be noted that the out-of-plane stress ( ) zz r σ is taken as the intermediate principal stress 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 Here we are interested in the particular dependence of the cohesion ( ) C α , shown in Figure 2, which introduces the linear hardening law. , 0 ( ) , 14) and the internal variable α is the plastic strain magnitude p ε .
Introducing (12) in (10) and replacing e ε in Hook's law gives Making the addition of (15) and (16): 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.
Replacing (18) in (19) the last one becomes: Considering a monotonic decreasing function in i P starting from P ∞ , we have the different behavior zones in the rock mass, show in Figure 3. A first phase, when i P remains close enough to the P ∞ , the entire rock mass remains in elasticity, having only zone 1. However, when i P is below a limit value, necessary for the stress state to reach the yield criterion at i r R = , the plastic zone 2 with linear hardening behavior begins. When i P decreases further, a third zone, of perfect plasticity, may appear. The solution in this zone is classic (Mandel, 1966). We define the radius r y = those where the elastic zone begins.
The conditions 0 p ε = and 0 F = in r y = allows obtaining the value of the constant ' A using (18) and (19): Using the condition 0 p ε = in the elastic domain and the boundary condition in r = ∞ we obtain the complete solution in this zone:

Solution in the plastic zone 3 ( )
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 As far as the yield criterion is zero in this zone, the equilibrium equation (3) becomes The integration of the (19) together with the boundary condition of (4) leads to the following equation of the radial stress: Replacing (30) into (18) we obtain the expression of the radial displacement given by Afterwards, using the equilibrium equation and the plane state condition of deformations, the rest of the solution is given by

Solution in the plastic zone 2 ( ) x r y ≤ ≤
In this zone of plasticity with hardening behaviour the yield criterion is written as As far as this criterion is null in this zone, the integration of the equilibrium equation with the condition in r x = leads to where 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 In the same way that for the zone 3, the expression of the displacement is obtained using (18), which gives: Using the equilibrium equation and the plane state condition of deformations, the rest of the solution is given by:

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 1 0 Substituting (40) in (20) we obtain The value of y is finally obtained writing the continuity of the radial stress in r y = , i.e., Together with (36) and (41) the value of y is obtain in an explicit way: with the following constants:

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. The characteristics of the following example are those of the Boom clay (Bernaud, 1991) obtained from samples at 240m deep: Figure 5 illustrates the convergence curve,  The Figure 6 shows the plastic radius for the internal pressure range.  As before, we can show the excellent agreement between the analytical and numerical solution. In Figure 7, through ( ) r θθ σ , the three zones that delimit the elastoplastic regime can be clearly seen.

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 and ϕ 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 F ≠ , where where β ϕ is the material dilatation angle. The internal variable associated with the isotropic hardening behavior is given by: ( Here we are interested in the particular dependence of the ( ) H α in (55), shown in Figure 8, which introduces the linear hardening law.
The equilibrium equation ( Considering a monotonic decreasing function in i P starting from P ∞ , we have the different behavior zones in the rock mass. A first phase, when i P remains close enough to the P ∞ , the entire rock mass remains in elasticity, having only zone 1. However, when i P is below a limit value, necessary for the stress state to reach the yield criterion at i r R = , the plastic zone 2 with linear hardening behavior begins. When i P 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.

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.
These conditions applied to the (58) and (59) gives the following solution: Due to the nullity of the yield criterion in this zone, the equilibrium equation writes Due to the fundamental (58) we can eliminate p ε . So we obtain Equation where the second member can be written as a function of rr σ , given in the (63). So, an integral-differential equation in rr σ is expressed as  The method of resolution (74) uses technics purely mathematics. Therefore, we expose in the following directly the solution obtained: The constants 1 c and 2 c are obtained taken into account the boundary and continuity conditions of stresses, given by:

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 p ε and of the yield criterion in r y = in the (58): The r y u r r r = ∂     ∂   deduction can be performed using the (67) and (80) This solution is only valid as long as the stress values i σ maintain the condtion 0 α α ≤ .

Second case: two plastics zones
In this case, the second zone (perfect plasticity) appears in i r R = (see Figure 10). In this zone the yield criterion writes: Due to the nullity of F in this zone, the equilibrium equation is: The resolution of the differential equation above with the limit condition ( ) gives the classical solution: [ ] All the others constants are the same as of the first calculation. We remember: The displacements in this zone are given by the (90) with the constants of (91) to (95). The constants 1 c and 2 c must be calculated with the (106) and (107). We remember this equation:

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 p ε in x and y : The derivation in function of r for the (67) written in r y = and the value of y σ allow to write the above equation as a function of γ and x : where / y x γ = and the constants: 2 2 1 2 9 ( 1) 1 (1 )(1 ) (1 ) 1 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 γ given by the (119) in the (113). With a numerical calculation (calculation of the zero of a function) we find the value of / y x γ = for a given i σ . This value introduced in the (119) gives the value of x and consequently the value of y .

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: → . Figure 11 illustrates the convergence curve for a given application. The Figure 12 shows the plastic radius for spectrum of internal pressures.   (Bernaud, 1991).   The maximal difference between analytical and numerical results is 2% indicate the good agreement of the results.

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.