Coding the “New Implicit Method” in MatLab for preliminary tunnel design Implementação no MatLab do “Novo Método Implícito” para pré-dimensionamento de túneis

Resumo The tunnel design is a subject that demands a three-dimensional analysis dealing with two different players: lining and rockmass. The traditional methods of getting the necessary parameters for design do not consider a correct interaction between the two players, but the New Implicit Method (NIM) takes this characteristic in its core and develops some formulations for elastic, plastic and viscous rockmasses. Understanding that the stiffness of the lining and the distance between the lining and the tunnel face change the convergence of the tunnel, a code in MatLab for NIM is validated through Finite Element (FEM) with its results being presented. The validation of this method was compared with FEM analysis and the results obtained an average accuracy of 12% what represents a good approximation regarding geotechnical issues.


Introduction
The design of tunnel support is one of the most complex tasks in the tunnel design, and most of the projects remain empirically analyzed and designed. However, it is getting unsuitable, since it is being proved that the interaction of the soil and lining matters in a way to change the stresses taken by the concrete as well as the deformation caused by the rockmass. This deformation of the rockmass over the lining is the main parameter and it is going to be called convergence. The Association Française des Tunnels et de L'Espace Souterrain (AFTES) [1] classifies the tunnel support design methods into four types: (i) purely empirical methods indicating the most appropriate type of support for a situation defined from various geotechnical classification systems; (ii) methods for determining the loads acting on the support, regardless of support type and deformation; (iii) support design methods which consider loads exerted by the ground as input data with some integration of support stiffness, deformation and the reactions of the surrounding ground; (iv) and more recent methods taking full account of the ground/ support interaction. The first method consists on the reproduction of succeeded tunnel supports and does not provide an assessment of the safety factor. The second and third method derives from the strength of the materials and consider the tunnel a semicircular beam reacting to external loads. The support reactions are given as Winkler springs to reproduce the strength of the soil. The last classification is derived from the continuous media theory and this is the reason why it was chosen to be applied in this work. This classification integrates the traditional Convergence-Confinement Method (CV-CF) for analytical solution and some numerical solutions with Finite Element which can describe better the 3D effect around the tunnel face. The New Implicit Method (NIM) was developed by BERN-AUD and ROUSSET [3] as evolution of the Convergence-Confinement Method. The NIM fully integrates the influence of the lining on the final convergence at equilibrium, since it shows that stronger the lining is, the smaller is the convergence (dimensionless displacement). Regarding that, the paper's main objectives are: (1) to present the New Implicit Method and (2) to introduce the NIM in MatLab for preliminary tunnel design. Every relevant and recurrent parameter shown in the text are summarized on Appendix I at the end.

Interaction problem
The determination of the load transferred to the support requires an analysis of the interaction of the load-deformation characteristics of the elements comprising the system ground/support. The preliminary design of tunnels foresees the necessary strength of the lining since it is able of presenting the convergences of the rockmass. Let us consider a circular deep tunnel (radius R i ), driven in a homogeneous and isotropic rockmass [ Figure 1], initially subjected to a geostatic stress field (σ 0 ): (1) The Equation (1) gives the initial stresses (before the ex-caThe Equation (1) gives the initial stresses (before the excavation). It is dependent on the P ∞ which is the pressure exerted by the soil mass above the tunnel. According to PARK [9] since it is a circular deep tunnel, it is accepted to consider the same pattern for vertical and horizontal axes. Parameters ρ and g are the unit mass and gravity acceleration respectively. With the assumption that the lining is a ring of constant thickness e, made of a homogeneous and isotropic material; and, moreover, that the lining is set at a distance d 0 from the plane and vertical tunnel face, the problem becomes an axisymmetric problem dependent only on the radial distance r to the tunnel axis and the time t. In the other hand, the deep section of the tunnel can be treated as a Plane Strain problem since the lining speed is constant (time dependent), so the parameters will depend only on the radial distance r as well. Firstly, the most used method to predict the stresses upon the lining is the Convergence-Confinement method [1] which is a procedure that allows the load imposed on a support installed behind the face of a tunnel to be estimated. Per CARRANZA and FAIRHURST [5], when a section of support is installed in the immediate vicinity of the tunnel face, it does not carry the full load to which it will be  Coding the "New Implicit Method" in MatLab for preliminary tunnel design subjected eventually; a part of the load is redistributed around the excavation face itself. However, as the tunnel and face advance, this 'face effect' decreases and the support must carry a greater proportion of the load that the face had carried earlier.
The method explained in the previous paragraph is defined by two curves [ Figure 2]: the convergence curve is the dashed line, also known as the ground reaction curve that gives the internal pressure exerted at the wall versus the convergence of the tunnel's wall and it depends on the characteristics of the rockmass and the adopted plastic criteria; the second curve is the confinement one (or support characteristic curve), which illustrates the relationship between the pressure P i versus the closure of the lining U i s defined as the straight line in [ Figure 2]. The Confinement curve starts at the convergence of the installation of the lining U o and its slope is defined as the lining strength K s . The parameters of Figure 2 A are: U o and P 0 : the convergence (radial deformation) and the pressure at the distance d 0 from the tunnel face; U eq and P eq : the convergence and the pressure at a section far from the tunnel face at equilibrium, where both pressure and convergence are stabilized; R p : the plastic radius is the limit between the elastic and plastic zone around the tunnel. Figure 2 illustrates the relationship between the pressure P i versus the closure of the lining U i s (convergence). Since it is shown in Figure 2 (A) that the rockmass presents an elastoplastic behavior, there is a plastic zone around the tunnel (B). The convergences U o and U eq are shown in (B) and the pattern of initial increase and final stabilization of the rockmass displacement is shown as well. Eventually, the Convergence-Confinement curve is drawn and the safety factor can be obtained by comparing the concrete strength with the stress at a distance far from the tunnel face that is imposed by the rockmass upon the lining structure.

The lining stiffness
PANET [8] brings equations for the lining stiffness (K s ). There are two main equations which describe the lining stiffness. They are derived from a thin shell and a thick pipeline. This work uses the normal solution only, considering that in deep tunnel, there is a uniform convergence only. For other cases, Panet [8] shows the lining stiffness for bending resistance. Considering a thin shell and R i / e >10, the lining stiffness is: In the other hand, when a thick pipeline is used (R i / e <10), the lining stiffness is:

The New Implicit Method (NIM)
The Convergence at equilibrium is the final displacement imposed on the lining ring and it has been calculated by the NIM, the Pressure at equilibrium is defined by the NIM too. The Method derived from the CV-CF method contains four main steps and for every rockmass behavior (elastic, elastoplastic and viscoplastic) the steps follow the same sequence (see Figure 3). The New Implicit Method (NIM) was developed by BER-NAUD and ROUSSET [3] and consists on deducing the curve of convergence with a geometrical transformation from unsupported to supported tunnel. The transformation

Figure 2
Tunnel equilibrium in a diagram Pi x Ui (a) and equilibrium of a tunnel (b) (BERNAUD and ROUSSET, 1996, [3]) is a function a 0 (x) of the unsupported tunnel and it was approximated by: is the convergence at the tunnel face of the unsupported tunnel and U i f (∞)=U ∞ is the convergence of the unsupported tunnel far from the tunnel face. The NIM showed that a s (x) depends on the K s (lining stiffness) that had not been predicted by the CV-CF method; so, a function for the supported tunnel is proposed: Then, for a given constitutive law of rockmass, the application of the NIM can be summarized by the following steps: a) Choose the reduced lining function α^* (K s ' ).Then, for the supported tunnel, Equation (5) becomes: b) Solve the convergence of the tunnel face U f according to the rockmass behavior. c) The solution of the problem is then obtained by solving the following equations:

Figure 3
Sequence of solution of a tunnel convergence U eq and Pressure P eq

Figure 4
Calculation of a tunnel with the NIM Coding the "New Implicit Method" in MatLab for preliminary tunnel design (8) d) In the System (8), the convergence U 0 is a function of U eq (that explains the name of 'implicit method'). By (7) and (8) we can write: (9) Figure 4 sums up how the NIM method works integrating the two charts which correlate the convergence, the pressure and the location of the studied section by setting the lining supported function and getting the value for U 0 and its respective pressure. Another use of the simplified method is to determine the ratio λ 0 (or the convergence U 0 ) at the moment of lining installation and then use this value as a good approach to calculate the tunnel in a 2D plane strain condition with a more complex geometry and loading (anisotropic loading for instance). Figure 5 shows how the pressure P i is obtained by using the ratio λ 0 .

The Numerical Code GEOMEC91
The tunnel face advance process is studied here by the 'deactivation/reactivation' method developed in the numerical code 'GEOMEC91' [4] It is a 2D axisymmetric analysis of tunneling that considers the sequences of excavation and lining placement as it is illustrated in Figure 6. The excavation process is simulated by the 'deactivation' of the stiffness of the elements representing the ground removed during the excavation step (face advance). This is accomplished by a significance reduction of their Young's Modulus and Poisson's ratio. Tunnel support is placed using the reactivation process at a specified distance, d 0 , from the advancing face. This is calculated by substituting the lining characteristics in the corresponding lining elements. At this moment, the lining elements are stress and strain free. BERNAUD [4] validated this numerical code with the experimental values of the Boom clay in Belgium, and for a 230-m deep tunnel the in-situ convergence was about 2.20% and the numerical solution returned a convergence of 2.04% presenting an average difference of 0.16% (between GEO-MEC91 and experimental measures) showing the accuracy of the CEOMEC91 used here as the validation tool of the New Implicit Method (see Figure 6). QUEVEDO [10] and FIORE et al [6] compared the solution of the GEOMEC91 with the numerical solution given by software ANSYS and showed an accuracy of GEOMEC91 used here. The model used by Fiore presented an accuracy

Figure 5
Application of the ratio λ 0 in a plane strain model

Figure 6
Typical mesh of GEOMEC91 (BERNAUD and ROUSSET, 1996 [4]) of 5% when compared to the GEOMEC91 solution. The values for the convergence and stress at equilibrium are two of the driving characteristics to build a reliable and safe tunnel, and as the text follows, the solution of the incognitos cited before are going to be given for different behavior of the rockmasses depending upon the geology, type and creep characteristic of the rock around the lining.

Elastic rockmass
Some dimensionless parameters are going to be introduced to reduce the number of parameters shown in text. For an elastic solution, the dimensionless parameters are: Where: P ∞ : Initial geostatic stress as defined in Equation (1); E: Modulus of Young of the rockmass; K s : stiffness of the lining. The steps to follow by the NIM to find the convergence and pressure at equilibrium are as follow: a) Get the entry parameters: P ∞ , E, K s , d 0 , ν (Poisson coefficient of the rockmass) and R i ; b) For an unsupported tunnel, we have: (11) c) For elastic, and plastic solutions, the function a s (x) is given from Equation (7) in function of a dimensionless lining function (α*), and it is illustrated in Figure 7 proving that this function is dependent on the lining stiffness: The approximate method by BERNAUD and ROUSSET [3] proposed an average function for the reduced lining function, so α* = α* average . See Figure 7: (13) d) So, the function a s (x') from Equation (7) is going to be changed by replacing the average function of α* from Equation (13): (15) e) Calculate the convergence and stress at equilibrium: (17) f) The convergence at d 0 is given by Equation (9).

Plastic rockmass
The rockmasses subjected to plastic strain presenting cohesion and/or friction angle will be treated here with Tresca (ϕ = 0) and Coulomb (ϕ > 0) criteria. However, some other criteria like Hoek-Brown is valid and presented in the original text of NIM [3]. Coding the "New Implicit Method" in MatLab for preliminary tunnel design

Tresca criterion
This criterion is valid for rockmasses whose main geotechnical characteristic is given by the cohesion. The steps to follow by the NIM to find the convergence and pressure at equilibrium are as follow: a) Get the entry parameters: P ∞ , E, K s , d 0 , C, ν = 0.5 and R i . b) For Tresca, the Friction Angle (ϕ) is zero and a new parameter is introduced: c) The reduced dimensionless lining function is given in Equation (13). d) Calculate the convergence U ∞ and at the tunnel face U f , as follows: (20) e) The a s (d 0 ') is the same as the shown one in Equation (14); The solution of this criterion valid for P eq < P ∞ -C; if P eq > P ∞ -C the rockmass is elastic. So, solving the System (8), the convergence and the stress at equilibrium are as follow: g) The plastic radius is: For Mises, the solution is given by the System (8) with C replaced by .

Coulomb criterion
This criterion is valid for rockmasses whose main geotechnical characteristic is given by cohesion and friction angle when . The steps to follow by the NIM to find the convergence and pressure at equilibrium are as follow: a) Get the entry parameters: P ∞ , E, K s , d 0 , C, ϕ, ν = 0.5 and R i ; b) For Coulomb, two new parameters depending on the friction angle (ϕ) are introduced: When the Friction Angle is zero, the solution is for Tresca criterion shown in Equation (22). c) For Coulomb, the reduced lining function (13) is changed by the insertion of the influence of the friction angle on it: (28) d) The Equation (28) changes the Equation (7)  (32) g) The solution of the System (8) gives the convergence and stress at equilibrium:

Viscoplastic rockmass
The NIM for viscoplastic rockmasses was developed as the same way as the elastic and plastic rockmasses, so the solution is given by getting a transformation equation from unsupported to supported tunnel. The reduced velocity of excavation V* is the new parameter for viscoplasticity and it is written as: Where: η: viscosity constant.
The transformation function obtained from the unsupported tunnel is changed and it is given as follows: (37) The parameters A and B in Equation (37) Since it is a viscoplastic solution, a criterion for the plasticity surface is required, the solution shown in the next topics is given for Tresca criterion presented before.

Tresca solution
a) The convergence U ∞ is: The solution of the System (8) is: (52) c) The plastic radius is given by Equation (25). Figure 8 illustrates the influence of the reduced velocity on the convergence profile of the tunnel, showing that the faster it gets, the smaller is the convergence. According to BERNAUD and ROUSSET [3], when the velocity trends to infinity, the solution trends to the elastic convergence, when there is no enough time to generate residual strain.

Results and validation
The code is developed in MatLab and gives the .gui executable. The code has been developed with the equations of the New Implicit Method shown before, and afterwards some response charts and solutions are given to compare the solution of NIM with GEOMEC91, the finite element tool described in Item 4.1.

Figure 8
The influence of the reduced velocity on the convergence. (BERNAUD and ROUSSET, 1994 [3]) Coding the "New Implicit Method" in MatLab for preliminary tunnel design

Elastic solution
The first results are shown for elasticity used for instantaneous strain pattern, we adopted: The solution obtained for the convergence through the Finite Element Analysis (U eq FE (%)) was taken from the GEOMEC91 solution introduced by [4]. The errors between the FEM and NIM are shown and presented good approximation as shown in Table 1. U eq FE,NIM : solution for the convergence at a far distance from d 0 through Finite Element Method and NIM method analysis respectively. Figure 9 shows the results obtained through the .gui executable for the second d 0 ' shown in the Table above, showing the charts developed automatically by the executable for stiffness varying from 0 to 20 000 MPa.
The executable generates three charts which are the solutions for many supporting stiffness, it also generates a .xlsx file with the numerical solution for each supporting stiffness, the dashed line represents the solution for P ∞ ' = 0.008, d 0 ' = 0.67 and K s ' = 72.

Plastic solution
The results shown in Table 2 concerns the plasticity with a Tresca and a Coulomb criterion. The plastic solution is used in over consolidated soils which present residual strain. For this analysis, we adopted:

Figure 9
Solution for an elastic rockmass using .gui MatLab The results for Coulomb and Tresca criteria present some errors smaller than 12% which is acceptable if compared to the results of the Finite Element analysis, indicating that the NIM is an adequate method for predicting the convergence and stress at equilibrium.

Viscoplastic solution
The results for viscoplastic cases are using the Tresca criterion as the plasticity surface. Some soils, mainly clays, present a creep behavior what increases the strain after some time and this behavior is important to assess. The results shown here have been compared to the ones given by GEOMEC91 as well. The results for many viscoplastic examples are given in Table 3.
As it is shown above, the method works well for every type of condition. Figure 10 shows the results of the executable for a viscoplatic rockmass, the three charts given by the .gui are in function of the reduced velocity instead of the lining stiffness as it was for the elastic and plastic solutions shown before, the dashed line shows the solution for N s = 5, d 0 ' = 1, V* = 500 and K s ' = 0.72 matching with the solution shown is Table 3.

Conclusion
A supported tunnel behaves as a three-dimensional structure for which strain and stress fields of the surrounding ground are strongly influenced. This paper develops an easy tool executable with MatLab to solve the problem of convergence and stress at equilibrium for elastic, elastoplastic and elastoviscoplastic tunnels.
It is important to know beforehand some geotechnical characteristics of the rockmass and the loading pattern to select the stress-strain criteria correctly. It is known that if the rockmass does not present permanent strain after being stressed it is going to behave as an elastic model, however the existence of creep and permanent deformations implies a plastic or a viscoplastic trend which give different solutions. The problem was proved to be slightly different from the one solved by the Convergence-Confinement method, since the New Implicit Method highlights that the stiffness of the support is a great player on the convergence U 0 (at a distance d 0 from the tunnel face). A point to be added is to take the permanent deformation of the supporting, for instance the concrete one, since it may change the curvature of the confinement curve.
The accuracy of the new method is satisfactory for a  Coding the "New Implicit Method" in MatLab for preliminary tunnel design geotechnical preliminary study, regarding the numerical limitations. The average error stays at a maximum of 12%, with most of the solutions under 8% when compared to the numerical solution. Generally, most of the differences are being for the safety, since it is increasing the value of the stress at equilibrium which is imposed over the lining. The results are good for values of K s ' greater than one. In addition to that, the results for smaller d 0 ' present better convergence than the ones of greater distance of lining placement. These smaller values represent most of the excavated tunnels, such as the value of d 0 ' ≅ 0,16 used in Tunnel Paraíso in Sao Paulo (Mafra, 2011, [7]). So, for conditions of K s ' > 1 and d 0 ' < 1, most of the solutions converge better to the numerical solution. The numerical code GEOMEC91 presents an approximation of actual behavior. However, considering the validation of GEOMEC91 with clays shown before and that it is a realistic tool to preview the behavior of deep tunnels, it was used to evaluate the solutions given by the New Implicit Method encoded here in MatLab.

Figure 10
Solution for a viscoplastic rockmass through .gui MatLab

APPENDIX I
The most relevant parameters used in this paper are summarized in the following list. The parameters subscripted with (L used in E L , ν L and so on) are defined regarding the lining characteristics. The superscript (' used in Ks',d 0 ' and so on) indicates that such parameter is dimensionless.
The other remaining parameters are listed on work presented before. They are referred as soon as they show up, and the list above are with the most recurrent and important parameters to take note.