J. González-Co, F. Vrs, F.G. Bstnte, L.R. Alejno,?

aDepartment of Natural Resources and Environmental Engineering, School of Mines, University of Vigo, 36.310 Vigo, Spain

bDepartment of Applied Mathematics, University of Vigo, 36.310 Vigo, Spain

1. Introduction

The convergence-confinement method (CCM) is a technique that quantifies the interplay between tunnel and installed support in terms of strains and stresses. The method relies on three components: the ground reaction curve (GRC), which describes the relationship between diminishing internal pressure and deformation on the tunnel spring-line; the longitudinal deformation profile,which relates the deformation on the spring-line to distance to the face; and the support characteristic curve, which represents the stress-strain relationship for the support system. These together enable the designer to estimate the performance of the support system. In geomechanical practice, the method is usually applied, for design purposes, to tunnels where stresses are expected to surpass rock mass strength, that is, deep tunnels excavated in average to poor quality rock masses.

The first CCM developments were analytical solutions for materials with simple behaviors such as elasto-perfectly plastic(EPP) rock masses (Duncan-Fama, 1993; Panet, 1993; Carranza-Torres and Fairhurst, 1999). More complex analytical models were subsequently developed to simulate elasto-brittle (EB) materials(Carranza-Torres, 1998, 2004). It is usually not possible to find the analytical solutions for the problem of a tunnel excavated in more complex materials, such as strain-softening (SS) materials;therefore, the GRC has to be obtained using a numerical approach.There are a number of approaches to obtain the GRC for SS rock masses (Alonso et al., 2003; Guan et al., 2007; Lee and Pietruszczak,2008; Park et al., 2008; Wang et al., 2010; Zhang et al., 2012). For a homogeneous material, the most efficient solution is numerical approximation of the self-similar solution, which requires a single initial value associated with an ordinary differential equation system to be solved. Note that it has been shown that the response of rock will differ depending on the selected behavior model (Alejano et al., 2009a).

The curves can also be obtained for any type of material using 2D or 3D numerical models, but the process is typically inefficient and very time-consuming (usually requiring several hours).Moreover, observed under particular circumstances are mesh dependency and bifurcation and localization phenomena that produce non-symmetric solutions, which significantly increasecomputation times. If there is bifurcation, the symmetric solution is no longer stable. For practical purposes, however, symmetric solutions derived from self-similar approaches are still relevant to a homogeneous material (Varas et al., 2005). This indicates that even if this problem may slightly affect results in 2D models, the variations used to be not very relevant at the scale of the problem. When dealing with 1D models, as in our case, localization cannot take place. So it is expected that symmetric solutions for non-homogeneous materials, as shown below, behave in the same way.

Under particular circumstances, a tunnel of radius a or an underground excavation may be surrounded by an aureole or ring-like zone of radius Rcmade of a material whose behavior differs from that of the rock mass (Fig. 1), where an external radius, Rf, large enough can be selected. In such a way, the elastic conditions are encountered at that distance of the excavation. The GRC of this case can be solved as a series of problems where an internal pressure, pi, is applied to the boundary of the tunnel and this pressure is reduced from σ0, the field stress, to a null value.

This may happen when there is an aureole of blast-damaged rock around a tunnel or when a regular pattern of reinforcement (rock bolts or cables) is installed so that the bolted rock behaves like a stronger material. Other situations reflecting different types of reinforcement, such as grouted zones or forepole umbrellas, might also be considered under this approach. In these cases, the rock mass is not homogeneous and solutions are not selfsimilar, so it is impossible to use the traditional rapid approaches to obtain the GRC. It would nevertheless be convenient to have a rapid technique available for obtaining the GRC of the excavation when the two different materials appear concentrically around the tunnel. We therefore sought to develop an efficient numerical approach for obtaining a symmetrical solution for such cases.

Fig. 1. Unloading scheme for an excavation of radius a, in a non-homogeneous rock mass under a hydrostatical stress σ0. r is the radial coordinate and τ is the evolution parameter.

The remainder of the article is laid out as follows. Section 2 states the problem and describes the mathematical background, the rock mass behavior model, and the characterization approach adopted to estimate all the parameters. Section 3 presents the numerical solution to the problem. Section 4 describes the technique for estimating GRCs in tunnels surrounded by an aureole of blast-damaged material (along with some examples) and Section 5 describes an approach for reinforced tunnels. The article closes with a discussion of our conclusions, and also an appendix that verifies the approach for simple cases.

2. Statement of the problem

2.1. The unloading problem

The unloading problem for a particular excavation section Ω can be formulated as the search for a stress field σ, depending both on the spatial coordinates, r, and on an additional parameter, τ,which measures the evolution of the unloading process. Note that this additional parameter will also serve as the natural parameter in the formulation of elasto-plastic material behavior in terms of the rate-independent, incremental plasticity theory.

Neglecting inertial terms, the solution to the unloading problem over section Ω leads to a continuum of equilibrium problems,associated with the unloading conditions on the excavation wall Γ and suitable asymptotic conditions:

where→n is a unitary vector normal to the excavation wall; the dot represents derivatives with respect to the unloading evolution parameter τ, which, in line with the previous comment, can also be considered to be the plasticity evolution parameter.

In order to effectively compute the GRC, we now consider a semi-discretization of the unloading problem. In particular, we assume that the stress field associated with a pressure pkion the wall excavation has already been computed and we consider a further decrease in the pressure from pkito pk+1i. This semidiscretization requires solving an equilibrium problem:

with a boundary condition associated with the unloading on the excavation wall:

where

In the semi-discretized problem stated above, incremental stresses Δk+1σ will be given by integration of the SS material behavior, in terms of the incremental plasticity theory. More precisely,Δk+1σ will be obtained through:

starting from stresses (and internal variables related to the plastic state) at the k-th step. It should be remarked that the load path appearing in this integral (connecting the plastic states corresponding to τkand τk+1) must, in fact, be computed as part of the solution.

In order to initialize computation of the GRC, suitable initial conditions corresponding to the stress field before excavation must be described. These initial conditions are dependent on the spatial variable in non-homogeneous rock masses.

Let us consider, for example, an excavation in an otherwise homogeneous rock mass where a cylindrical region (corresponding to the tunnel and a surrounding ring-like region) is composed of a different material. If this rock mass (composed of two homogeneous materials) is subjected to hydrostatic stress, a radial distribution of stresses (including a jump in the circumferential stress at the contact radius) will appear (see Fig. 2). From this profile, we obtain the initial stresses (including the initial internal pressure on the excavation wall) to compute the unloading curve.

In some cases, the corresponding stress profile (for an elastic solution) can be obtained analytically. For more complex cases,computation times. If there is bifurcation, the symmetric solution is no longer stable. For practical purposes, however, symmetric solutions derived from self-similar approaches are still relevant to a homogeneous material (Varas et al., 2005). This indicates that even if this problem may slightly affect results in 2D models, the variations used to be not very relevant at the scale of the problem. When dealing with 1D models, as in our case, localization cannot take place. So it is expected that symmetric solutions for non-homogeneous materials, as shown below, behave in the same way.

Under particular circumstances, a tunnel of radius a or an underground excavation may be surrounded by an aureole or ring-like zone of radius Rcmade of a material whose behavior differs from that of the rock mass (Fig. 1), where an external radius, Rf, large enough can be selected. In such a way, the elastic conditions are encountered at that distance of the excavation. The GRC of this case can be solved as a series of problems where an internal pressure, pi, is applied to the boundary of the tunnel and this pressure is reduced from σ0, the field stress, to a null value.

This may happen when there is an aureole of blast-damaged rock around a tunnel or when a regular pattern of reinforcement (rock bolts or cables) is installed so that the bolted rock behaves like a stronger material. Other situations reflecting different types of reinforcement, such as grouted zones or forepole umbrellas, might also be considered under this approach. In these cases, the rock mass is not homogeneous and solutions are not selfsimilar, so it is impossible to use the traditional rapid approaches to obtain the GRC. It would nevertheless be convenient to have a rapid technique available for obtaining the GRC of the excavation when the two different materials appear concentrically around the tunnel. We therefore sought to develop an efficient numerical approach for obtaining a symmetrical solution for such cases.

Fig. 1. Unloading scheme for an excavation of radius a, in a non-homogeneous rock mass under a hydrostatical stress σ0. r is the radial coordinate and τ is the evolution parameter.

The remainder of the article is laid out as follows. Section 2 states the problem and describes the mathematical background, the rock mass behavior model, and the characterization approach adopted to estimate all the parameters. Section 3 presents the numerical solution to the problem. Section 4 describes the technique for estimating GRCs in tunnels surrounded by an aureole of blast-damaged material (along with some examples) and Section 5 describes an approach for reinforced tunnels. The article closes with a discussion of our conclusions, and also an appendix that verifies the approach for simple cases.

2. Statement of the problem

2.1. The unloading problem

The unloading problem for a particular excavation section Ω can be formulated as the search for a stress field σ, depending both on the spatial coordinates, r, and on an additional parameter, τ,which measures the evolution of the unloading process. Note that this additional parameter will also serve as the natural parameter in the formulation of elasto-plastic material behavior in terms of the rate-independent, incremental plasticity theory.

Neglecting inertial terms, the solution to the unloading problem over section Ω leads to a continuum of equilibrium problems,associated with the unloading conditions on the excavation wall Γ and suitable asymptotic conditions:

where→n is a unitary vector normal to the excavation wall; the dot represents derivatives with respect to the unloading evolution parameter τ, which, in line with the previous comment, can also be considered to be the plasticity evolution parameter.

In order to effectively compute the GRC, we now consider a semi-discretization of the unloading problem. In particular, we assume that the stress field associated with a pressure pkion the wall excavation has already been computed and we consider a further decrease in the pressure from pkito pk+1i. This semidiscretization requires solving an equilibrium problem:

with a boundary condition associated with the unloading on the excavation wall:

where

In the semi-discretized problem stated above, incremental stresses Δk+1σ will be given by integration of the SS material behavior, in terms of the incremental plasticity theory. More precisely,Δk+1σ will be obtained through:

starting from stresses (and internal variables related to the plastic state) at the k-th step. It should be remarked that the load path appearing in this integral (connecting the plastic states corresponding to τkand τk+1) must, in fact, be computed as part of the solution.

In order to initialize computation of the GRC, suitable initial conditions corresponding to the stress field before excavation must be described. These initial conditions are dependent on the spatial variable in non-homogeneous rock masses.

Let us consider, for example, an excavation in an otherwise homogeneous rock mass where a cylindrical region (corresponding to the tunnel and a surrounding ring-like region) is composed of a different material. If this rock mass (composed of two homogeneous materials) is subjected to hydrostatic stress, a radial distribution of stresses (including a jump in the circumferential stress at the contact radius) will appear (see Fig. 2). From this profile, we obtain the initial stresses (including the initial internal pressure on the excavation wall) to compute the unloading curve.

In some cases, the corresponding stress profile (for an elastic solution) can be obtained analytically. For more complex cases,

Fig. 2. Example of radial and tangential initial stress (σrand σθ) distributions(remark the jump in σθ) and initial displacement (ur) for the case of a tunnel (a = 2 m)surrounded by an annulus (Rc= 3 m) of a softer material.

In some cases, the corresponding stress profile (for an elastic solution) can be obtained analytically. For more complex cases,where plasticity arises before excavation begins or where variations in material properties occur, an analytical solution may not be possible, so a numerical approach will be required.

An asymptotic condition associated with the field stress is used in the statement of the unloading problem. Instead of using a large computational domain to ensure this condition, asymptotic behavior can be imposed by matching the elastic solution at an artificial external boundary (in the sequel it will be taken as r = Rf) chosen beyond the reach of the plastic zone. This will reduce the computational load.

Due to the axisymmetric geometry, the problem can be written in terms of radial and circumferential stresses as functions of the radial coordinate:

where X refers to the cohesion and the friction. Fig. 4 shows the peak and residual failure criteria and their impact on the stress-strain behavior of the rock mass, together with the typical functions of c and φ as defined by Eq. (11). This approach can be implemented

with the unloading condition on the excavation wall (at r = a):

Matching with the elastic solution at the external radius r = Rfis then provided by

In the methodology developed in Section 3, the material parameters are allowed to vary arbitrarily in the radial coordinate(corresponding to the axial symmetry assumed throughout this paper) and a simpler configuration (corresponding to a homogeneous rock mass except for the presence of a ring-like region of damaged or reinforced material) will be used in the examples.

In order to compute the unloading curve, once the problem has been semi-discretized, a set of equilibrium problems must be solved (Eqs. (2) and (3)). This can be viewed as a continuation method in which internal pressure plays a role of the continuation parameter. Although natural, this approach can be quite inefficient.In practice, very large displacements at the excavation walls are to be expected for a rock mass with strong SS behavior. As a consequence, the corresponding GRC (obtained by plotting the pointafter solving each equilibrium problem) can exhibit a very shallow step in the final stretch of the curve.Consequently,the step in internal pressure variation must be very small in order to ensure convergence in the equilibrium problem solution, as a good initial iteration will otherwise be difficult.

As an alternative, the incremental displacement at the wall can be used as the control parameter in order to compute the unloading curve. Even better is to use an arc-length continuation technique,where the distance in the ur(a)-piplane between two consecutive steps is used as a control parameter.

Fig. 3. Rock mass behavior models.

2.2. Rock mass behavior models

Rock mass behavior models can be classified according to postfailure behavior in EPP, EB and SS materials. SS can be considered as the most general case, since EPP and EB are simply particular types of SS: EPP is a type of SS in which the peak and residual failure criteria coincide, and EB is a particular case of SS in which stress jumps from the peak strength to the residual strength without strain relaxation during failure (Fig. 3). Since SS behavior covers every type of rock mass, it is the focus of the present study.

2.2.1. Strain-softening behavior

SS behavior can be derived from the incremental, rateindependent plasticity theory (Hill, 1950), described by means of a failure criterion, f, and a plastic potential, g, both of which depend not only on the stress tensor σijbut also on the so-called softening(or plasticity) parameter η.

For the sake of simplicity, we have chosen the Mohr-Coulomb failure criterion, although the method can be easily adapted to any other failure criterion and plastic potential. Accordingly, the general forms of the failure criterion and the plastic potential in this case are given by

where

where c is the cohesion, φ is the friction and ψ is the dilation angle.

The proposal by Alonso et al. (2003) is followed to define the evolving Mohr-Coulomb failure criterion. The cohesion c and the friction φ are characterized by means of bilinear functions of the plastic parameter η, in such a way that, starting from their peak values (cpand φp), they linearly decrease to their residual values(crand φr) for a particular value of the plastic parameter (η*). For larger values of this parameter, c and φ will keep a constant residual value. Accordingly, the functions representing these values take the following form:in standard numerical codes, for instance, FLAC (Itasca Consulting Group, 2009).

To fully define material behavior according to the theory of incremental plasticity, it is necessary to define the plastic parameter. Following Alonso et al. (2003) in an approach also adopted by other authors (Lee and Pietruszczak, 2008; Park et al., 2008), we used the plastic shear strain as the plastic parameter:

Fig. 4. (a) Peak and residual failure criteria of a SS material and stress-strain response of a test on a sample and (b) evolution models for cohesion and friction.

This plastic parameter is easy to be computed and can be related to other plastic parameters, like those suggested by Vermeer and de Borst (1984) or used in FLAC (Itasca Consulting Group, 2009).Using the plastic potential defined by Eq. (9), the flow rule imposes

for some˙λ to be computed, related to the type of loading (Lubliner,2006).

2.2.2. Estimating rock mass parameters for strain-softening behavior

Once it is accepted that rock mass behavior is SS, all the parameters that control this behavior need to be provided. For Mohr-Coulomb SS as defined above, the parameters needed to fully characterize a rock mass include peak strength criteria (cpand φp),elastic parameters (E and ν), residual strength criteria (crand φr)and post-failure parameters (η* and ψ).

Traditional rock mechanics approaches, based on rock mass classification systems and laboratory tests on rock samples, are used to obtain the needed parameters including peak strength criteria and elastic parameters. Post-failure parameters can be obtained, following recently suggested guidelines, as described below.

Classic rock mass characterization starts from compressive tests on rock samples in the laboratory, which provide the values for unconfined compressive strength (σci) and the Hoek-Brown parameter m.

Field work enables the quality of the rock mass to be assessed according to one or various classification systems. For the sake of simplicity, we chose the GSI (Marinos and Hoek, 2000). For characterization purposes also needed was an estimate of the density of the materials and, to obtain an estimate of the field stress, the depth of the excavation to be investigated.

With this information, typically available from any excavation study, and using the approach proposed by Hoek et al. (2002) and implemented in the freeware code RocLab (Rocscience Inc., 2009),one can obtain a reasonable estimate for the peak strength parameters (cpand φp) and the elastic parameters (E and ν).

The residual strength parameters can be obtained based on an estimate of the residual value for the GSIres. Based on Cai et al.(2007), Alejano et al. (2012) suggested roughly estimating this value according to

Post-failure parameters are estimated following the approach by Alejano et al. (2010). The value of η*, i.e. the plastic parameter marking the transition to the residual state, is estimated for every level of confinement stress σ3as

where

where speakis the Hoek-Brown parameter of the intact rock mass,E is the Young’s modulus for the rock mass, and Kψdepends on ψ. This dilation, considered constant for this application, can be estimated as

Note, in regard to Fig. 4, that in the axial branch of the stress-strain curve, the value of η* can be related to the negative slope of the post-failure component of this curve M (related to ω and E according to M = - ωE).

An aureole of material surrounding the tunnel (formed of blast damaged material) is considered here for a particular application.In order to quantify the parameters controlling the behavior of this damaged material, a simple (although somewhat rough) approach based on the characterization by Hoek et al. (2002) can be applied.

A parameter D, reflecting the level of disturbance, is included as an indicator of the damage suffered by the rock mass, roughly estimated in accordance with the blasting procedure used. This parameter takes a value between 0 (undamaged material) and 1(very damaged material), usually 0, 0.5 or 0.8. Thus, D = 0 when

excellent quality-controlled blasting or mechanical excavation produces minimal disturbance to the rock mass surrounding the excavation, D = 0.5 when squeezing problems result in significant floor heave and D = 0.8 when a high level of damage is associated with very poor quality blasting that produces severe local damage in the surrounding rock mass (Hoek et al., 2002).

Once D is defined, the same characterization procedure as defined previously for an undamaged rock mass can be followed,with the new D value included in the expressions provided by Hoek et al. (2002).

3. Numerical solution

Using a weak formulation, the (k + 1)-th step in the unloading problem stated in the previous section indicates an incremental displacement field Δk+1ursuch that the corresponding incremental stress Δk+1σ solves for any virtual displacement field υ. This formulation (corresponding to the principle of virtual work) can then be particularized to axisymmetric solutions.

Alternatively, a weak formulation of the axisymmetric problem can be directly obtained. The weak formulation for the axisymmetric problem reads (in the second case) as follows: find a radial displacement (incremental) field Δk+1ursolution of

for any virtual displacement field υ.

This problem is discretized (in space) using a standard Lagrange P1 finite element method (FEM). Then, after dividing the interval (a, Rf) using a (non-uniform) grid {rj}j=0,N(with

a = r0＜ r1＜ . . . ＜ rN-1＜ rN= Rf), we look for the (approximated) val

ues of the radial displacement (incremental) field Δk+1urat the nodes ({rj}j=0,N, Δk+1ur,j) as the solution for the system of nonlinear equations, which can be obtained from the weak formulation when the virtual displacement field is taken as each one of the functions in the finite element basis {φj(r)}j=0,N(de fined by{φj(ri)} = δij, where δijrepresents the Kronecker delta).

In the FEM discretization, the integrals over each element will be approximated using Gaussian quadrature. If a one-point scheme is used (as should be done when using a P1 FEM discretization of the displacement field), quadrature leads to the rectangle formula using the central point. As a result, radial and circumferential(incremental) stresses, Δk+1σrand Δk+1σθ, are only computed at the centre of each element.

To determine radial and circumferential (incremental) stresses,Δk+1σrand Δk+1σθ, and the components of the tangent stiffness matrix, Kk+1t, at the points used in the numerical quadrature (i.e. the central points of the elements), the constitutive law must be integrated for some prescribed deformation path. The (incremental)displacement field is estimated for each iteration in the solution of the non-linear problem and the associated deformation path (along the unloading step) is computed at the quadrature points. The constitutive law is then integrated along this path to determine the(incremental) stresses at these points. A residual of the weak formulation is obtained to assess the quality of the approximation, and a new estimate is computed if the residual is not small enough.

3.1. Integrating constitutive equations

As explained in the previous sections, in order to solve the kth step in the unloading problem, an incremental displacement field Δk+1urmust be found such that the corresponding incremental stresses, Δk+1σrand Δk+1σθ, solve the equilibrium equation.When using a finite element discretization, incremental stresses(associated with a given incremental discrete displacement field and characterized by incremental displacement at the mesh nodes)must be computed at each numerical integration (quadrature)node.

This prediction is trivial if material behavior is elastic. Thus, in order to integrate the constitutive equations for the material, an elastic trial is first computed (assuming that material behavior is elastic along the whole increment in the displacement field associated with the unloading step Δk+1ur). If the final stress given by the elastic trial remains in the elasticity domain, the trial is accepted and the incremental stresses Δk+1σrand Δk+1σθare computed.Otherwise, the total incremental displacement field is decomposed into an elastic part (associated with the fraction of the incremental displacements needed by the elastic trial to reach the yield surface)and a plastic part (the rest).

More precisely, given the incremental strain field Δk+1ε associated with the incremental displacement field Δk+1u, the elastic trial is computed as Δk+1σe= CΔk+1ε (where C stands for the tensor of elastic constants). If the stress field updated using Δk+1σelies outside the elastic domain, Δk+1ε is decomposed as

where Δk+1εe= αΔk+1ε and the scalar α is determined so that the updated stress field (and internal variables) associated with the elastic prediction for Δk+1εelies on the yield surface. In this case,the rate-independent plasticity constitutive law for the plastic part Δk+1εp(starting from the yield surface) must be integrated.

As previously mentioned, to find the value of the incremental stresses (at each quadrature node), the integral in Eq. (4) can be used, where the derivative of the corresponding stress will be given by the constitutive law (in this case, the SS behavior described in Section 2.2). Thus, the relationship between stress and (plastic)

As previously mentioned, to find the value of the incremental stresses (at each quadrature node), the integral in Eq. (4) can be used, where the derivative of the corresponding stress will be given by the constitutive law (in this case, the SS behavior described in Section 2.2). Thus, the relationship between stress and (plastic)strain given by the flow rule (see Eq. (13)) can be inverted (Lubliner,2006):where Cepstands for the elasto-plastic stress-strain tensor(dependent on the stress tensor and the plastic parameter η). In our computation, since the incremental strain field is known (and,hence, initial and final values of ε are prescribed, and a simple path with constant strain rates can be constructed), the previous equation can be seen as a system of ordinary differential equations in the stress field σ and in η (where the definition of η, given by Eq. (12)and written in differential form, must be incorporated in this system). As a consequence, standard integration schemes for ordinary differential equations can be used to obtain Δk+1σ (Sloan, 1987).

Alternatively, the incremental stresses can be directly computed by seeking a final plastic stage associated with the imposed incremental displacement and compatible with the constitutive law.Some details and considerations of benefits and drawbacks of both approaches will be shown in the following sub-sections.

3.1.1. Explicit algorithm: sub-stepping scheme

Algorithms to integrate the constitutive law (at each quadrature node) using the numerical integration of Eq. (2) are usually known as sub-stepping techniques (Sloan, 1987; Sloan et al., 2001). The plasticity evolution parameter (i.e. the independent variable in the system of ordinary differential equations) can be arbitrarily chosen to vary in (0, 1). Note that, in the present computation, the integration interval corresponds to the transition associated with the k-th unloading step.

The integration interval must be discretized in a number of substeps that can be determined in advance or (in more elaborated algorithms) that can be computed along the integration using some kind of step control strategy (based on an estimate of the local truncation error). In order to keep the low computational cost, explicit schemes are usually considered.

A rather simple alternative for solving Eq. (2), for the plastic part of the incremental strain tensor, with the additional equation for the plastic parameter η, is

Here we use the modified Euler scheme with constant time steps, with Δεp(τ) = Δk+1εp(τ).

This approach for the incremental stress computation can be safely used only when rock mass behavior (as described in Section 2.2) is true SS behavior. It is important to note that strong SS (i.e. when df/dη is negative and has a large absolute value) can result in a singular tensor Cepin Eq. (22). This makes it impossible to use Eq. (2) to calculate the incremental stress field; it also causes the jump in the stress field that characterizes EB behavior.It should be observed that, in the context of the material behavior described in Section 2.2, the strength of the SS rock mass, for a given set of peak and residual parameters, is related to the residual plastic parameter η*. Thus, a critical value ηcrit, depending on the peak and residual parameters, can be defined so that SS behavior only results if η* ＞ ηcrit(see Alonso et al. (2003) for further details).In the particular case of the Mohr-Coulomb failure criterion, ηcritis given by

where

where Kppis as in Eq. (10) referred to peak friction angle φp.

3.1.2. Implicit algorithm

As explained above, when material behavior is EB (detectable by comparing the plastic parameter associated with the transition to the residual state η with the corresponding critical value ηcrit,as shown in Eq. (23)), the incremental stress field can no longer be computed by integrating the system of ordinary differential equations.

In this case, the final stress associated with EB behavior is computed by seeking final stress corresponding to the residual state and compatible with the (imposed) incremental plastic strain. This approach will define an implicit algorithm (in contrast with the explicit algorithm resulting from the use of explicit schemes to integrate the system of ordinary differential equations) and will require the use of an iterative method to find the solution for the resulting non-linear problem (in the values of the incremental stress field at the quadrature node under consideration).

Concerning the solution for the resulting non-linear system, a very effective strategy in most cases is using the Newton method with line search. Nevertheless, if the drop from peak values to residual values is large, convergence may be difficult with this method(despite additional robustness provided by line search). In these cases, a more robust technique, based on Powell’s hybrid method(Powell, 1970), is used.

3.2. Solution for the non-linear problem

At the k-th step of the unloading problem, the finite element discretization of the equilibrium equations leads to the solution of a non-linear problem in the displacement field at the nodes of the mesh. This equation is the discrete counterpart of Eq. (18), which represents the principle of virtual work. More precisely, a vector of nodal incremental displacements Δk+1u must be found as the solution for a system of non-linear equations that can be written as

where Δk+1F represents the balance of internal forces determined from the integration of the constitutive law at each quadrature point, starting from an initial state described by the fields (σk, ηk)corresponding to the previous unloading step and Δk+1b stands for the external forces associated with regional stress and with internal pressure at the excavation wall.

Since the continuation method, associated with the computation of the unloading curve, could provide (if the GRC is regular)a good initial guess for the incremental displacement field at each step (several previously computed steps can be used to elaborate this guess), using Newton’s method can result in a fast global algorithm. At the (j + 1)-th iteration of Newton’s method, once the j-th approximation of the incremental displacement field Δk,ju has been computed, a new approximation is computed to solve the linear system:

The tangent matrix Ktk+1,j, which is related to the derivative of the incremental stresses at the quadrature nodes with respect to the incremental displacements at the mesh nodes, must be computed from the information obtained in the integration of the constitutive laws. For instance, it can be extracted from the last substep in the case of the explicit sub-stepping technique described in the previous sections.

The previous linear system in the correction of the incremental displacement field Δk+1,j+1u - Δk+1,ju is usually rewritten as

Table 1Rock mass strain softening parameters for the intact and damaged rock masses.

where rk+1,j= Δk+1b - Δk+1F(σk, ηk; Δk+1,ju) represents the residual of the non-linear system of equations (i.e. the unbalance between internal and external forces in the j-th iteration) which is used to check convergence.

The convergence of the Newton algorithm relies on the regularity of the solution for the unloading problem. This algorithm works well for most but not all SS materials.

In brittle or guasi-brittle materials, abrupt changes in stresses,may appear and a irregular behavior of the unloading problem can be anticipated. In these cases, as done in the integration of the constitutive equation using an implicit algorithm, a more robust algorithm based on Powell’s hybrid method is used.

4. Application to a tunnel surrounded by blast-damaged rocks

Blasting is commonly used for drift and tunnel advancing in civil and mining engineering fields, particularly when the rock is hard and in a good or very good quality rock mass. However, blasting not only breaks the rock the engineers want to remove but also damages or disturbs the surrounding rocks. This undesirable effect has been studied widely from different perspectives (Holmberg and Persson, 1979; Raina et al., 2000; Ouchterlony et al., 2002). However, although a good number of studies analyze the nature of this damage, few assess the impact on tunnel design. Noteworthy is the works by Saiang and Nordlund (2009), who considered a reduction in the elastic parameters and analyzed their influence on tunnel response. A rather simple way to account for blast-induced damage is the approach proposed by Hoek et al. (2002), based on quantifying what is called the disturbance parameter D (Section 2.2). This method, together with the rock mass characterization approach described here, permits the estimation of all the parameters needed to model the behavior of the damaged material found around the excavated tunnel.

The extent of the blast-damaged zone also needs to be known for modeling purposes. Even if it is clear that the degree of damage decreases from the excavation surface to the interior of the rock mass, approximate estimates of the extent of damage can becomputed using observation and classic approaches (Anl?ggnings AMA-98, 1999).

Table 2Parameters of the explicit algorithm for GSI = 25 and 40.

Table 3Parameters of the implicit Powell algorithm for GSI = 60.

Another possible approach, developed by García-Bastante et al.(2012), shows how an estimate can be based on Langefor’s theory of blast calculation, starting from the explosive energy, the coupling factor, the rock mass constant and the gas isentropic expansion factor.

Accounting for the typical values for these parameters in actual blasting tunnel advance, the extent of the damaged zone may vary between 0.5 m and 2.5 m but usually takes values in the range of 1-2 m.

For illustrative purposes, we modeled 3 m radius tunnels excavated in different quality rock masses in which different degrees of blast damage were induced. All the rock masses were composed of standard rock. The laboratory value for uniaxial compressive strength was 75 MPa, for m = 10 and an estimated depth of 600 m(σ0= 15 MPa). Rock masses with low (GSI = 25), low to average(GSI = 40) and average to high (GSI = 60) geotechnical quality were modeled. For each type of rock mass, we computed the value of the plastic radius for the unloaded case and the GRC for four different levels of blast-induced damage associated with different values of D (0, 0.5, 0.8 and 1). These damage levels can be related to different situations as in Hoek et al. (2002).

For this general approach, the extent of the damaged zone was considered to vary in accordance with the damage level, from 1 m for the case D = 0.5 to 2 m for the case D = 1 (very damaged rock).All the parameters needed to model these excavations were calculated following the procedure described in Section 2.2. Recall that dilatancy is considered to be constant, the parameters obtained for the different rock masses and blasting-induced damage levels are presented in Table 1. The algorithm parameters for GSI = 25 and 40 are shown in Table 2 and those for GSI = 60 in Table 3. Results for the different GRCs and plastic radii Rpare shown in Fig. 5.

Note that, as would be expected, in all cases the plastic radii and final displacements increase with higher levels of blasting-induced damage. Note also that the displacement increments are larger for poorer quality rock masses, where the influence of the damage may be more relevant in terms of support needs.

These results indicate noticeable differences in expected displacements and plastic radii, depending on whether blasting-induced damage around a drift or tunnel is considered and, if considered, its magnitude. The approach described above enables the influence of this damage to be assessed within the CCM approach.

Alejano et al. (2009b), using real data for a Peruvian mine,showed that poor blasting techniques increased the need for support and reinforcement in mining drifts. This increased need for support and reinforcement can be quantified by means of the approach described here and within the framework of the application of the CCM to estimate support needs (e.g. Oreste, 2003).Note that the fact of computing the GRC of a drift or tunnel while accounting for blasting-induced damage can lead to better and less costly support design and can help decide whether blasting controlled techniques are a suitable option.

Fig. 5. GRCs in terms of piand displacement urand final plastic radius, for the tunnels excavated in different quality rock masses (GSI = 60) surrounded by different types (size and damage level) of blasting-induced damaged material.

5. Application to a tunnel surrounded by a ring-like zone of reinforced rock

Indraratna and Kaiser (1990a, 1990b) proposed an analytical model to represent the behavior of a reinforced rock mass near a circular underground opening, considering a proper interaction mechanism between the ground and the grouted (or friction) bolts.To analyze the influence of the bolt pattern on tunnel response,they introduced the bolt density parameter, β, as a dimensionless parameter that reflects the relative density of bolts on the tunnel perimeter, and considered the shear stresses opposing the rock mass displacements near the tunnel wall. β is defined as where d is the bolt diameter; SLand STare the spacings of the bolt pattern in the tunnel advance direction and in the tunnel section,respectively; andˉλ relates the mobilized shear stress acting on the grouted bolt to the stress acting normal to the bolt. This value is analogous to the coefficient of friction or the bond angle used to analyze reinforced earth and split-set bolts.

The dimensionless ratio β/ˉλ is the inverse function of the bolt spacing for a given bolt and tunnel geometry and it generally ranges between 0.1 and 0.4. Once this parameter is defined, the above named authors proposed obtaining equivalent Hoek-Brown parameters for the reinforced rock mass by means of the expres

sions:

Indraratna and Kaiser (1990a) described how to derive the analytical model and illustrated the effect of bolts on the stress and displacement field near an opening. Verification of the theory by laboratory simulations and field measurements was subsequently described in detail. This type of geometrical distribution of rock and reinforced rock proposed to include reinforcement in GRC behavior can be used within the frame of our methodology.

As an example application, a laboratory simulation performed by Indraratna and Kaiser (1990b) was studied by means of the proposed methodology and compared to results obtained by these same authors. A small 130 mm radius tunnel was excavated in artificial rock with the following properties: E = 1500 MPa, ν = 0.25,φ = 32°, σc= 3.5 MPa, s = 0.9 and Kψ= 2.

The tunnel was reinforced with 100 mm brass bolts, with an estimated friction factor ofˉλ = 0.36, subjected to an isotropic field stress of 14 MPa. Results were obtained for the unsupported tunnel(β = 0) and for open and closed bolting meshes (β = 0.145 and 0.291,respectively).

Table 4Rock mass strain softening parameters for the intact and damaged rock masses.

Table 5Parameters of the implicit Powell algorithm for GSI = 60.

For these cases, following the approach by Indraratna and Kaiser(1990a), the parameters presented in Table 4 were obtained. Note that in this case peak and residual frictions were considered equal and brittle behavior was considered for cohesion, in such a way that the residual value of the compressive strength was computed from the expression:

Remind that cohesion and compressive strength could be related according to

The algorithm parameters are presented in Table 5. The GRCs for the different cases together with the values of the plastic radius for the fully unloaded case are shown in Fig. 6. The results compare quite well those obtained by Indraratna and Kaiser (1990a) in their models and are checked against actual laboratory results.

Even if the results have been compared to those for a small laboratory model, they still seem to indicate that the approach to include the effect of reinforcement in the CCM is suitable for real engineering applications as suggested by Oreste (2003).

Note that other ground improvement techniques such as grouting or the use of forepole umbrellas can also be dealt with in the framework of this approach.

Fig. 6. GRCs in terms of piand displacement and plastic radius, for the case of a 0.13 m tunnel excavated in a rock mass surrounded by different types (bolt density parameter β) of rings of reinforced material.

6. Conclusions

We have described a methodology for solving the axisymmetric problem associated with the computation of the GRC and plastic radius around tunnels excavated in non-homogeneous axisymmetric rock masses. The method is a robust and fast approach that can be used when self-similar solutions are not possible due to the non-homogeneity of the rock mass.

The method solves the sequence of equilibrium problems associated with tunnel unloading. A non-homogeneous stress distribution associated with the regional stress field needs to be computed and used as the initial condition in the problem. Assuming axial symmetry of this solution, the equilibrium equation in radial coordinates is solved by means of 1D FEM over an adapted

mesh.

Coupling with the elastic (far- field) solution is used as a boundary condition at an external boundary located at a sufficient distance from the tunnel wall.

A robust strategy to integrate the constitutive law is adopted in combination of an explicit technique (based on a sub-stepping algorithm) for the general cases (EPP and SS) and an implicit (and more time consuming) technique for rock masses behaving in a brittle way, i.e. when the stress drop occurs brusquely.

The Appendix of the manuscript contains both a verification of the developed code and a validation of the numerical aspects of the algorithm. In this sense, the proposed numerical method to deal with circular excavations in non-homogeneous, axisymmetric strain-softening rock masses can be considered as validated.Our method has been moreover verified for homogeneous and non-homogeneous cases and two applications of practical interest are presented. The first application deals with the computation of GRCs around tunnels surrounded by an aureole of blast-damaged materials. Results for different quality rock masses submitted to different damage levels point to an increased need for support and reinforcement in the case of poor blasting techniques. The second application refers to a tunnel surrounded by a ring of rock bolt-reinforced material. The resulting GRC shows diminishing deformation associated with the stabilizing effect of the rock bolts.Numerical results are in good agreement with laboratory tests. The two examples demonstrate the usefulness of the proposed numerical approach, capable of making predictions in an efficient way.

Future work includes comparing numerical results with in situ monitoring data.

Acknowledgements

The authors thank the Spanish Ministry of Science and Technology for financial support awarded under Contract Reference Numbers BIA2009-09673 and MTM2010-21235-C02-02. Ailish M.J.Maher is gratefully acknowledged for language edition of a version of the manuscript.

Appendix A. Verification

In this Appendix, computation of the GRC using the proposed methodology will be considered for verification purposes in following three test cases:

(1) A circular excavation in a homogenous EPP rock mass, where the solution is compared with the self-similar solution computed using the method proposed by Alonso et al. (2003).

(2) A circular excavation in a rock mass where an aureole of damaged zone is presented (both materials exhibit EPP behavior), where two different plasticity laws (Mohr-Coulomb and Drucker-Prager) are considered and where the solution is compared with analytical and 2D numerical solutions.

(3) A circular excavation in an EB homogeneous rock mass, where the solution is compared with the self-similar solution computed using the method proposed in Alonso et al. (2003).

Table A.1Rock mass parameters for homogeneous case.

A.1. Test case 1

A circular tunnel of 2 m radius is excavated in a homogeneous strain softening rock mass submitted to an isotropic field stress of 4.07 MPa, equivalent to 160 m in depth. This stress field is considered as the external boundary condition. Parameters used in this test case are shown in Table A.1 and the resulting GRC is shown in Fig. A.1, where the GRC computed from the self-similar solution is also represented. As can be observed in Fig. A.1, both solutions agree (note that no self-similarity is imposed here but naturally arises from the problem symmetry).

Concerning the computation time, although the proposed methodology must solve a sequence of (1D) equilibrium problems instead of integrating an initial value problem associated with a reduced system of ordinary differential equations (as is the case of the self-similar solution), the computational cost is not very much increased. For instance, in this example using an implementation with GNU/OCTAVE the computation time using the approach in Alonso et al. (2003) is 1.3 s, while the corresponding time using the proposed methodology rises to 6 s.

Fig. A.1. GRCs in terms of piand displacement for the homogeneous rock mass verification case.

A.2. Test case 2

In this test case, first considered is a Mohr-Coulomb material with EPP behavior. In this case, σ0= 4.07 MPa, Rc= 2.2 m and a = 2 m (the excavation radius). So a damage aureole of 0.2 m is considered. Material properties (both for the undisturbed and disturbed materials) are summarized in Table A.2. Computed radial and circumferential stress fields corresponding to the complete unloading (pi= 0) are plotted in Fig. A.2, also illustrating the analytical solution. As shown, the computed and analytical solutions agree.

Table A.2Rock mass parameters for EPP non-homogeneous case.

Table A.3Rock mass parameters for Drucker-Prager model.

The second case computed refers to a Drucker-Prager material with EPP behavior and associated flow rule, defined as in Regueiro and Borja (1999). In this case, σ0= 15 MPa, Rc= 3.113 m and a = 2 m.So a damage aureole of 1.113 m is considered.

Material properties are summarized in Table A.3. The solution computed with the proposed methodology is compared with that obtained with the FEM code TAHOE using a 2D model. TAHOE(Sandia National Laboratories, 2003) is an open source FEM code developed for the simulation of complex material behavior. The resulting GRC is represented in Fig. A.3.

Fig. A.2. Stress field for a tunnel section for the non-homogeneous EPP rock mass verification case.

For this case, where a rather fine mesh has been used to exhibit the large difference in computation time when comparing the proposed method with a general solver for a 2D model, the computation time for the proposed method rises from a few seconds in the previous test cases to one minute; meanwhile, the computation time used by TAHOE to solve a quarter of the 2D geometry,for a mesh size similar to the size used in the proposed method, is around five minutes.

Table A.4Rock mass parameters for homogeneous-brittle case.

Fig. A.3. GRCs in terms of piand displacement for the Drucker-Prager model.

Fig. A.4. GRCs in terms of piand displacement for the homogeneous-brittle rock mass verification case.

A.3. Test case 3

A circular tunnel of 3 m radius is excavated in a homogeneous EB rock mass submitted to an isotropic field stress of 14.07 MPa,equivalent to 560 m in depth. This stress field is considered as the external boundary condition. Again, the computed solution is compared with that obtained using a self-similarity approach (as done in Alonso et al., 2003).

The parameters used in this test case are shown in Table A.4 and the resulting GRC is represented in Fig. A.4. The method described in this paper is shown to be able to deal with this singular behavior by switching to an implicit integration of the constitutive law;meanwhile, the self-similar solution is again recovered as can be observed in Fig. A.4.

Plotted also is the self-similar solution computed using the approach in Alonso et al. (2003), where elasto-brittle transition is much easier to detect and deal with.

Alejano LR, Alonso E, Rodríguez-Dono A, Fdez-Manín G. Application of the convergence-confinement method to tunnels excavated in rock masses exhibiting Hoek-Brown strain-softening behavior. International Journal of Rock Mechanics and Mining Sciences 2010;47(1):150-60.

Alejano LR, Alonso E, Rodríguez-Dono A, Ordo?ez C, Cordova D. Estimate of support and reinforcement cost increase associated to poor blasting practices in drifting. In: Proceedings of the 9th Int. Symp. on Rock Fragmentation by Blasting -Fragblast 9, Sept. 2009. London: Taylor and Francis Group; 2009b. p. 731-42.

Alejano LR, Rodríguez-Dono A, Alonso E, Fdez-Manín G. Ground reaction curves for tunnels excavated in different quality rock masses showing several types of post-failure behavior. Tunneling and Underground Space Technology 2009a;24(6):689-705.

Alejano LR, Rodríguez-Dono A, Veiga M. Plastic radii and longitudinal deformation profiles of tunnels excavated in strain-softening rock masses. Tunneling and Underground Space Technology 2012;30:169-82.

Alonso E, Alejano LR, Varas F, Fdez-Manín G, Carranza-Torres C. Ground reaction curves for rock masses exhibiting strain-softening behaviour. International Journal for Numerical and Analytical Methods in Geomechanics 2003;27(13):1153-85.

Anl?ggnings AMA-98. General materials and works description for construction work section CBC: Bergschakt. Stockholm: Svensk Byggtj?anst; 1999 (in Swedish).

Cai M, Kaiser PK, Tasaka Y, Minami M. Determination of residual strength parameters of jointed rock masses using the GSI system. International Journal of Rock Mechanics and Mining Sciences 2007;44:247-65.

Carranza-Torres C, Fairhurst C. The elasto-plastic response of underground excavations in rock masses that satisfy the Hoek-Brown failure criterion. International Journal of Rock Mechanics and Mining Sciences 1999;36(6):777-809.

Carranza-Torres C. Self-similarity analysis of the elasto-plastic response of underground openings in rock and effects of practical variables. Minneapolis, USA:University of Minesota; 1998, Ph.D. Thesis.

Carranza-Torres C. Elasto-plastic solution of tunnel problems using the generalized form of the Hoek-Brown failure criterion. International Journal of Rock Mechanics and Mining Science 2004;41(Supp. 1):1-11.

Duncan-Fama ME. Numerical modelling of yield zones in weak rocks. In: Hudson JA,editor. Comprehensive rock engineering. Pergamon: Oxford; 1993. p. 49-75.

García-Bastante F, Alejano L, González-Cao J. Predicting the extent of blast-induced damage in rock masses. International Journal of Rock Mechanics and Mining Sciences 2012;56:44-53.

Guan Z, Jiang Y, Tanabasi Y. Ground reaction analyses in conventional tunnelling excavation. Tunnelling and Underground Space Technology 2007;22(2):230-7.

Hill R. The mathematical theory of plasticity. New York: Oxford University Press;1950.

Hoek E, Carranza-Torres C, Corkum B. Hoek-Brown failure criterion-2002 edition.In: Proceedings of the 5th North American Rock Mechanics Symposium; 2002.p. 267-73.

Holmberg R, Persson PA. Design of tunnel perimeter blast hole patterns to prevent rock damage. In: Jones MJ, editor. Proc. Tunneling. London, UK: IMM; 1979. p.280-3.

Indraratna B, Kaiser P. Analytical model for the design of grouted rock bolts.International Journal for Numerical and Analytical Methods in Geomechanics 1990a;14:227-51.

Indraratna B, Kaiser P. Design of grouted rock bolts based on the convergence control method. International Journal of Rock Mechanics and Mining Sciences and Geomechanics Abstracts 1990b;27:269-81.

Itasca Consulting Group. FLAC2D-fast Lagrangian analysis of continua. Minneapolis,USA: Itasca Consulting Group; 2009.

Lee YK, Pietruszczak S. A new numerical procedure for elasto-plastic analysis of a circular opening excavated in a strain-softening rock mass. Tunnelling and Underground Space Technology 2008;23(5):588-99.

Lubliner J. Plasticity theory. Upper Saddle River, NJ: Pearson Education; 2006.

Marinos P, Hoek ET. GSI: A geologically friendly tool for rock mass strength estimation. In: Proc. GeoEng2000 Conference. Melbourne: University of Melbourne;2000. p. 1422-42.

Oreste P. Analysis of structural interaction in tunnels using the convergence confinement approach. Tunnelling and Underground Space Technology 2003;18:347-63.

Ouchterlony F, Olsson M, Bergqvist I. Towards new Swedish recommendations for cautious perimeter blasting. Fragblast - International Journal for Blasting and Fragmentation 2002;6(2):235-61.

Panet M. Understanding deformations in tunnels. Comprehensive Rock Engineering 1993;1:663-90.

Park KH, Tontavanich B, Lee JG. A simple procedure for ground response curve of circular tunnel in elastic-strain softening rock masses. Tunnelling and Underground Space Technology 2008;23(2):151-9.

Powell MJD. A hybrid method for nonlinear equations. In: Rabinowitz P, editor.Numerical methods for nonlinear algebraic equations. New York: Gordon and Breach; 1970.

Raina AK, Chakraborty AK, Ramulu M, Jethwa JL. Rock mass damage from underground blasting, a literature review, and lab- and full scale tests to estimate crack depth by ultrasonic method. Fragblast - International Journal for Blasting and Fragmentation 2000;4:103-25.

Regueiro RA, Borja RI. A finite element model of localized deformation in frictional materials taking a strong discontinuity approach. Finite Elements in Analysis and Design 1999;33(4):283-315.

Rocscience Inc. RocLab 1.0. Toronto: Rocscience Inc; 2009.

Saiang D, Nordlund E. Numerical analyses of the influence of blast-induced rock around shallow tunnels in brittle rock. Rock Mechanics and Rock Engineering 2009;42:421-48.

Sandia National Laboratories. TAHOE - users guide 3.4.1. Albuquerque, New Mexico:Sandia National Laboratories; 2003.

Sloan SW, Abbo AJ, Sheng D. Re fined explicit integration of elasto-plastic models with automatic error control. Engineering Computations 2001;18(1/2):121-54.

Sloan SW. Substepping schemes for the numerical integration of elastoplastic stress-strain relations. International Journal for Numerical Methods in Engineering 1987;24:893-911.

Varas F, Alonso E, Alejano LR, Fdez-Manín G. Study of bifurcation in the problem of unloading a circular excavation in a strain-softening material. Tunnelling and Underground Space Technology 2005;20(4):311-22.

Vermeer PA, de Borst R. Non-associated plasticity for soils, concrete and rock. Heron 1984;29(3):1-62.

Wang S, Yin X, Tang H, Ge X. A new approach for analyzing circular tunnel in strain softening rock masses. International Journal of Rock Mechanics and Mining Sciences 2010;47(1):170-8.

Zhang Q, Jiang BS, Wang S, Ge XR, Zhang HQ. Elasto-plastic analysis of a circular opening in strain-softening rock mass. International Journal of Rock Mechanics and Mining Sciences 2012;50:38-46.