A tridimensional finite element approach to model a tunnel with shotcrete and precast concrete

Resumo This paper describes a numerical simulation with 3D finite elements of a tunnel. The viscoplastic law of Perzyna represents the rockmass behavior. The concrete, shotcrete or precast, is modeled as a viscoelastic material through the Maxwell and Kelvin chain models. Finite element simulations are performed by incorporating subroutines for viscoelastic concrete models in the ANSYS code. The method to simulate tunnel excavations is by activating and deactivating elements in sequential steps. In the first part of the paper two validations are performed. The analytical solution and the deformation achieved on the stabilization in the ANSYS code are compared with an unlined tunnel. A lined tunnel, with an elastic and viscoplastic rockmass combined with an elastic lining, is compared with the results of the GEOMEC91 code. In the second part, it is compared the same tunnel with two different concrete lining for two chain models. Finally, it is modeled the Kielder experimental tunnel, which in situ measured data is available.

A tridimensional finite element approach to model a tunnel with shotcrete and precast concrete

Introduction
Stabilizing underground openings such as tunnels excavated in rockmass remain a major concern of geotechnical engineers dealing with this kind of structures.In tunnels, the rockmass deformation and the ground pressure on lining depend on the stress and characteristics of the rockmass as well as of the tunnel geometry, the lining stiffness and the time of the lining installation.Pressure variation on lining and strain are caused by the advance of excavation and by the time-dependent properties of the rockmass and lining.The lining is mainly focused on limiting the consequences of a pressure relief of the surrounding ground due to excavation process, and more specifically to maintain the tunnel closure within an admissible value compatible with appropriate subsequent operating conditions of the structure.Therefore, adequate support systems are essential.In this paper it is studied a deep tunnel that presents a nonlinear behavior.The inelastic deformation time rate effects are always present at some degree.In the study of structural components under static loading conditions at normal temperatures it is accepted that the time rate effects are generally not important and the conventional theory of plasticity models the behavior appropriately.However, some structures as deep tunnels, present a time-dependent behavior.In order to take in account this phenomenon the rockmass will be modeled with a viscoplastic law and the lining with a viscoelastic one.Combined with the finite element ANSYS code, the proposed constitutive laws are carried out for simulating the advancement of a tunnel and calculating its convergence as the excavation proceeds.In the first part of the paper two validations are performed.For an unlined tunnel, the analytical solution and the deformation achieved on the stabilization in the ANSYS code are compared.For a lined tunnel, validation is performed by comparing the results of the GEOMEC91 code with the ANSYS code.The results of the ANSYS code agree perfectly with the analytical solution and with the GEOMEC91 code.Finally, it is simulated a real experimental tunnel, the Kielder tunnel [1], which in situ measurement data is available.

The constitutive law for the rockmass
The description of the constitutive law for the rockmass reveals two fundamentals conception: the elasticity of the material and the viscoplastic flow rule that describes the evolution of the viscoplastic stress.The elastic-viscoplastic formulation is based in the decomposition of the tensor of total strain as shown in the equation ( 1): (1) Where ε e , ε vp and ε 0 are the elastic strain tensor, the viscoplastic strain tensor and the initial one.The relation between stress and strain is given by the Hooke's law, shown in equation ( 2): (2) where σ 0 represents the vector of initial stress.
In the study of the state of stress in a given viscoplastic material, the following situations can appear: If σ is inside the elastic domain, then the criterion of viscoplasticity and (3) If σ is outside or in the border of the elastic domain then F(σ)<0 and: (4) The rate of total strain is the sum of an elastic part and a viscoplastic one, equation ( 4).The rate of viscoplastic strain is written as an evolution law given by equation ( 5): (5) The function g is zero when F(σ)<0 .
The law of viscoplastic evolution is described in the equations 6 and 7. and they are in conformity with the formulation given by [2]: (6) (7) Where F(σ,α) is the criterion of viscoplasticity; G(σ) is the viscoplastic potential; α is the parameter of hardening; . a is the rate of the parameter of hardening; , n η are the constants of viscosity; and 0 F is the reference stress.

Concrete as a tunnel lining
The aim of the systems of structural linings, in the perimeter of the tunnel, is to support the displacement of the walls and the face in order to ensure the stability of the tunnel and give the protection against the eventual local breaks.
The shotcrete are created in the USA in the beginning of 20th century and afterwards it have been used as a layer of protection against corrosion, fire and also for the lining of tunnels, slopes, etc.
The shotcrete is actually the dominant technique in the mining excavation.It is convenient in the tunnels where the section changes with the tunnel excavation.
The precast concrete lining is largely used in the soft soils.This system is convenient in the mechanized excavation where a fast progress for the long tunnel lengths is required.In this case the precast concrete is an economic technique and is efficient as a tunnel lining.

The long term behavior of the concrete
The viscoplastic materials show a time-dependent stress-strain relation.
In the concrete there are immediate strain that come from the elastic behavior and long-term strain that are due to the viscous and shrinkage behavior.The concrete has some features that transform it in a very much dependent on the loading, on the time and on the environment where it is placed (humidity and temperature).
In this way, according to the fib Model Code 2010 [3], the total strain at time t of a concrete member uniaxially loaded at time t' with a constant stress σ c (t') is given by equation ( 8):

The constitutive law for creep
In this analysis, creep is assumed to be linearly dependent to stress.Then linear superposition is applied to creep of concrete under variable stress [4].In order to simulate numerically the effects of the creep and aging of the concrete structures, Bazant [5] have presented many options that includes: an algorithm to obtain stress in concrete for a given strain history, the evaluation of creep rate for concrete with aging based in Maxwell chain and a methodology to evaluate the time dependency spectrum of Kelvin chain.The method was based in the expanded Dirichlet's series for creep rate.The retardation time gives a rule for creep rate that is required for the analysis of large structures.In generally, the numerical formulations for the evaluation of the time-dependent strain of the concrete, takes in account two diferents effects: creep and shrinkage.The fib Model Code 2010 presents some equations for both effects.Adding portions of a strain history, due to small stress increments, after a time t, the creep law for an uniaxial stress is given by equation (9).
Where t is the current age of the concrete (in days), t' is the age of the concrete at loading; σ is the total stress; ε is the total strain; ε 0 is the inelastic strain including the shrinkage and the thermal dilatation; and J(t,t') is the creep function that means:

Shrinkage of the concrete
Shrinkage is the phenomenon where the concrete reduces its volume by losing water in the concrete mass during the process of hardening.
The shrinkage is a given deformation independent of the current stress.It increases when the thickness of the structure in contact with the environment decreases.
In the numerical formulation, the shrinkage is considered as a given deformation, as observed in equation (8).Then the strain due to shrinkage is subtracted of the total one.The total shrinkage, given by the fib Model Code 2010 [3] is given by equation ( 11): (11) Where ε cas is the autogenous shrinkage and ε cds is the drying shrinkage.
The formulation presented in the fib Model Code 2010 is a function of the environment humidity, of the concrete age, of the concrete age in the beginning of the shrinkage and of the thickness of the concrete structure.
A tridimensional finite element approach to model a tunnel with shotcrete and precast concrete

Viscoelastic formulation: basic rheological models
In accordance to Creus [4] the viscoelastic models are formulated by the union of two basic rheological elements: the elastic elements (springs) and the viscous (dashpots) one.For the spring element it is admited that the response is instant and reversible.For the dashpot element it is considered a Newton viscosity, where the rate of strain is proportional to the stress.
The springs and the dashpots can be combined to obtain diferents configurations as shown in Figure 1.
The generalized models are given by the combination of the models shown before: the generalized Maxwell model is composed by the parallel union of several Maxwell elements; and the generalized Kelvin model is composed by a serie union of several Kelvin elements.

Generalized Maxwell chain model
The model of the generalized Maxwell chain, shown in Figure 2, represents both responses: the instantaneous and the long-term one.According to Pande [6], in this model it is assumed that the structure is composed of several superimposed layers.Each of these layers presents different mechanical properties.The layers have the same total strain, producing different stress field that contributes for the total stress field.The procedure of the application of the generalized Maxwell chain is provided for a spectrum of a long-term shrinkage.The strain-stress relation, based in a Maxwell chain, is given by the equation ( 12). ( Where σ μ is the stress in the spring μ; t μ is the time of relaxation related to the element µ of the chain; E μ is the elasticity modulus in the spring µ, that depends of the age t of the concrete; and n is the number of elements (adopted n=5).
In the case of the Maxwell chain, the element μ = 5 has an important relaxation time: , to simulate a coefficient of viscosity tending to infinity this assumption permits to say that the fifth element of the chain is an isolate spring that obeys the Hooke`s Law and does not present a viscous behavior.It allows that the deformation results in an asymptotic curve, such as the real concrete behavior.The others relaxation dates have been chosen according to equation ( 13).The integration of the equation (12), results in an exponential serie that represents the relaxation function of the Maxwell's chain.
Then, for a given time, the values of E μ can be obtained from the relaxation function with discrete values.The determination of the relaxation function of creep was based in the procedure proposed by Bazant [5].
In order to obtain an accurate representation of the relaxation spectrum, eight loading ages have been adopted, beginning with t' 1 = 2.8 days and ending with t' 8 = 8854.28days.For each loading time, thirty discrete dates are considered.
With the Bazant's procedure, the relaxation curve is computed for any loading time t', using the least squares fitting method.The aging of concrete is considered by the update of the elasticity modulus, of each element of the chain, as time increases.

Generalized Kelvin chain model
In this model, shown in Figure 3, the parameters of viscoelasticity, the elasticity modulus E μ and the coefficient of viscosity η μ doesn´t depend on the time.These parameters are calculated at the beginning of an experimental model calibration and they remain constant during the analysis, as it is performed in the classical viscoelasticity.The aging is considered only in the isolate spring.This effect is calculated as a function of the elasticity modulus ranging as time passes by the method of fib Model Code 2010 [3].These facts lead to a major numerical simplification and make the Kelvin model The differential equation for the Kelvin chain without aging is given by the equation ( 14). ( where E μ is the elasticity modulus of the element µ of the chain; η μ is the coefficient of viscosity of the element µ of the chain; and ε μ is the strain of the element µ of the chain. Integrating the equation ( 14) for a constant stress σ applied at time t', in the chain of n Kelvin elements, results in a viscous strain ε vp as shown in equation ( 15). ( where θ μ is the retardation time, given by the equation ( 16). ( For the calibration of the generalized Kelvin model, it was considered fixed time steps given by equation ( 17).

(17)
The use of an exponential algorithm allows the gradually growth of the time steps until the values of time greater than the shorter retardation time, keeping stability and a good accuracy.

fib Model Code 2010: creep formulation
In the case of the linear viscoelasticity, the material obeys the principle of superposition applied to the history of stress and strain.In this way, the fib Model Code 2010 presents a formulation for the creep function that is valid to concretes submitted to a maximum stress smaller than 0.4•f cm (f cm is the mean compressive strength of concrete at an age of 28 days in MPa) and submitted to a loading applied in the time t', as show in the equation (10).The formulation takes into account the temperature, the humidity, the size of the structural element, the type of cement and the concrete strength.

Maxwell and Kelvin chain models calibration for viscoelastic concrete
For eight different loading application times t', the creep curves proposed by the fib Model Code 2010 [3] are compared with those obtained by computation models of Maxwell and Kelvin chains described before.Some of the curves obtained in the calibration are shown in Figure 4 and Figure 5.A tridimensional finite element approach to model a tunnel with shotcrete and precast concrete At curves of Figure 4 and Figure 5, it is shown that creep effect is a function of the concrete age when loading is applied and also is a function of the load duration.Concrete is classified as a viscoelastic material with aging.The curves obtained from the fib Model Code 2010 and from the numerical models have presented a very good agreement.

Computational implementation in the ANSYS code
The ANSYS finite element code can adapt the program to the user needs by the incorporation of subroutines written in Fortran language.In these simulations, Maxwell and Kelvin chains for the shotcrete and precast are compared.The viscoelastic constitutive law (Maxwell and Kelvin chains) for concrete is incorporated by using the subroutine usermat for the 3D case.This subroutine can be modificated by the user in order to implement its own constitutive equations for a given material.This subroutine is used for each integration point, for each element and for each iteration as shown in Figure 6 and Figure 7.
In Figure 6, EMU is the elastic modulus matrix for the selected eight loading ages and the i subscript represents each direction: xx, yy, zz, xy, xz, yz.

Elastic rockmass with elastic lining case
The aim here is to compare a closed form solution of the 1D tunnel in elasticity with the numerical results (ANSYS) and with a simplified method (NIM).For an unsupported circular tunnel, Corbetta [8] presented an analytical solution for convergence, in elasticity, at a section far from the tunnel face.It is given by the equation 18. (18) n where P ∞ is the geostatic stress at the depth of the tunnel.
Figure 8 illustrates the results of the ANSYS code, the GEOMEC91 code (is a 2D axisymmetric code [9]) and the NIM (1D solution, [10]).The results present a very good agreement.

Viscoplastic rockmass and elastic lining case
The aim here is to compare the closed form solution in plasticity of a 1D tunnel with the numerical results (ANSYS) and with a simplified method (NIM).For an unlined tunnel at stabilization, the solution is given by the analytical plastic solution for a Tresca rockmass (with a cohesion C) presented by Corbetta [8] at equation 19. ( Figure 9 shows a very good agreement between the analytical plastic solution and the numerical solution of a viscoplastic tunnel calculated at stabilization with the ANSYS code.For a supported tunnel, an analytical solution is not available and the validation is performed between the numerical ANSYS and GEOMEC91 codes, as shown in Figure 10.Far from the tunnel face the difference between the convergence values at stabilization is about 5%.The two codes have different viscoplastic formulations that may explain the small differences between them.

Numerical analysis of tunnels using generalized Maxwell and Kelvin chain models
Underground structures design has as priority the local stability of the structure and around it.The induced displacement must be tolerable for a given structure.The numerical methods, in civil engineering, allow the analysis of complex structures.The stress and the strain of the tunnel depends mainly of the constructive method; excavation method; thickness of the lining; dis-  A tridimensional finite element approach to model a tunnel with shotcrete and precast concrete tance between the lining and tunnel face; constitutive law of the rockmass and of the lining.A 3D computational method is necessary to analyze these conditions and the requirements of an underground structure [11].
The aim in this section is to simulate the 3D tunnel construction procedure and the lining placement using the ANSYS program.
The viscoplastic rockmass behavior is combined with the viscoelastic lining behavior.

The model in the computational ANSYS code
The tunnel construction simulation is performed with 37 excavation steps (tunnel with circular section), with the radius (20 where σ is the stress, is the rate of equivalent strain, m is the parameter related with the strain rate of the material, g is the pa- rameter related to the viscosity of the material and o σ is the yield stress.These parameters are shown in Table 1.In the ANSYS code, the activation/deactivation of elements corresponds to the "Birth and Death" command.The element is deactivated by multiplying its stiffness by a reduction factor and removing the mass of the global matrix [12].

Numerical results
In order to simplify the nomenclature of each case it is used the following abbreviations: for Kelvin we have the codes KE with a letter J or M corresponding to the shotcrete and precast concrete.For   Maxwell, it is used the code MA with the letter J or M. Figure 12 illustrates the convergence of all cases.So it is possible to compare the Maxwell and the Kelvin models in two times: the end of the construction and at five years later.Due to the long term behavior of the rockmass (viscoplasticity) and of the lining (viscoelasticity), it is observed that the tunnel deformation continues as time passes until the final stabilization.
In the first ages the deformation of the concrete is important.It is observed, in the curves corresponding to the final construction, that the convergence far from the face presents an important increasing.
In Figure 12, it can be observed that the convergences of Maxwell case are greater than Kelvin ones; and shotcrete convergences are greater than precast concrete ones, as expected.
In order to evaluate the behavior as time passes in a point of the lining, it is chosen the point with coordinates x,y,z equal to (in centimeter) -70.71; 70.71; 2883.33,shown in Figure 13a.
All cases are illustrated in Figure 13b.Maxwell results present a greater deformation than Kelvin results.
The precast concrete curves are below the shotcrete ones due to the fact that precast concrete has a bigger stiffness at loading.

Comparison between the models of Kelvin and Maxwell
It is possible to describe the long term behavior of the Kelvin model by using properties of E μ and η μ independents of time, which appears to be a considerable simplification of the problem.Then, the total viscous deformations are calculated independently of the previous steps and easily converge even changing the size of the steps.In Maxwell model the total viscous deformations are determined with an explicit formulation, which is very instable and sensible to the size of steps.The calculations with the model of Maxwell have been performed with a time step of one day, while in the model of Kelvin the time step is increasing and so advancing much faster.The calculation with the model of Kelvin takes 35 minutes while the model of Maxwell takes 7 hours (Computer Intel i7-4770K 3.50GHz with a memory RAM of 8 Gb).

Comparison between the shotcrete and precast concrete
The difference between shotcrete and precast concrete is given by the age that they are incorporate in the model.It is considered  A tridimensional finite element approach to model a tunnel with shotcrete and precast concrete that the precast concrete is placed at 28 days old while the shotcrete is placed at 2.8 days old.These parameters, which specify if concrete will be precast or shotcrete and the method to use, are indicated in the ANSYS input file.This data enters to the usermat subroutine like data material parameters.
Figure 12 shows the differences between the many results: the convergence at five years is greater than 12%, 10%, 14%, 12% for the cases KE-J, KE-M, MA-J, MA-M, respectively, in relation with the convergence at the end of the construction.This figure also shows the difference between the shotcrete and the precast concrete.At the instance of the end of entire construction the relation between the shotcrete and the precast concrete is of 4% for the Kelvin model and it is equal to 3% for the Maxwell model.At five years old, the difference between the shotcrete and precast concrete is equal to 6% for the Kelvin model and is 5% for the Maxwell model.

The experimental Kielder Tunnel modelling
In this section, it is compared the convergence of the experimental tunnel Kielder with the numerical simulation performed with the code ANSYS.
Ward [1] describes the history and the process of construction of the Kielder tunnel, giving the experimental results.The experimental tunnel Kielder was constructed in England in 1974 and was closed in 1979.It has a length of 32 km and a circular section with a 3.3 m diameter.
In order to compare the behavior of many systems of lining, seven different linings were chosen for the layer of mudstone and two methods of excavation were used, as shown in Figure 14.A length of 11 m was chosen for each stretch.
The instrumentation was placed 30 cm above the lining in the middle of the tunnel.The instruments measured the rock displacements for a period of five years.
In this calculation the Kelvin model has been used to simulate the shotcrete.The physical and geometrical parameters, used in the calculations, are based in the works given by Ward [1] and Freeman [13].In the model, the rockmass has a viscoplastic behavior and the lining is viscoelastic: the characteristics are presented in Table 3.
The tunnel radius is equal to R=1.65 m and the thickness of the lining is equal to e=0.143 m.The mesh has 3904 hexahedral elements and 16599 nodes.The measuring instrument was placed 30 cm from the lining, in the middle of the tunnel.So the displacement of a node in the middle of the length (11 km) and 47 cm far from the lining is plotted in Figure 15.The curves of Figure 15 show that the long term numerical results are very similar to the experimental results.At short term, the   difference between both methods is more important and it may be due to the fact that some concrete parameters are not specified in the literature.

Conclusions
The numerical code ANSYS proved to be a useful tool in the simulation of the excavation of a tunnel and placement of a lining by using the activation/deactivation method.The models of Maxwell and Kelvin [7] have been tested with the results of the fib Model Code 2010 [3].So it was possible to show that the methods of Maxwell and Kelvin are able to represent the viscoelastic behavior of the concrete.
The Kelvin model is more convenient and more economic because the properties E μ and η μ are time-independent and the concrete aging is calculated independently.This model doesn't have the problem of convergence, while the Maxwell model is significantly unstable due to the explicit formulation and the parameters update.Therefore, the formulation and the total computation time is greater at the Maxwell case.
In the case of the supported tunnel with a viscoelastic lining two kinds of concrete have been studied: the shotcrete and the precast concrete.The difference between the two concretes is the age at loading, leading to different elasticity modulus.For a same model these two linings present a small difference, about 3 to 6% at the analyzed instants.
Comparing the final construction instance with those at five years old, the difference varies between 10% to 14%.
Finally everything of what was done previously is used to model a real tunnel, the Kielder tunnel.In this case the numerical calculations are made with the Kelvin model because of the convergence advantages and smaller computation time.The convergence of the tunnel at five years old given by the ANSYS code is compared with the in situ measurements: they have presented a very good agreement.The relevance of a tunnel simulation with a finite element code is to make the study easier, allowing to choose the most advantageous configuration concerning the project and the construction.The AN-SYS code allows obtaining a wide variety of numerical results, with several visual options of post-processing.This is very important in a project analysis.The behavior of shotcrete and precast lining depends on time.The incorporation of viscoelastic models by subroutines of the ANSYS code enables to represent the long term behavior of concrete with more accuracy.For these reasons, these techniques are increasingly attractive at representation of numerical models for artworks (tunnels in this case) taking advantage of pre and post-processing skills of a commercial code.The code ANSYS, allow incorporating additional models specially developed for determined situation analysis.

Figure 5 -
Figure 5 -Creep function J for t'=28 days by fib 2010 and computational models

Figure 8 -Figure 9 -
Figure 8 -Convergence curves at validation in an isotropic rockmass, with a distance between the tunnel face and The face is flat and transverse to the longitudinal axis.The excavation speed is 10 m/day.The tunnel is deep and it is considered a geostatic stress equal to symmetry, only one quarter of the total model has been simulated.The measures are: 30 m in the longitudinal axis Z of the tunnel, 20 m in the axis X and Y .Element SOLID186 has been used with radial mapping for meshing the model.The mesh consists of 7504 elements and 34344 nodes as shown in Figure 11.The convergence curves after the entire construction of the tunnel and at the stabilization, which simulates the situation at an infinity time, are obtained.The rockmass is a non-linear viscoplastic material with a Mises viscoplastic criterion.The ANSYS program uses the Perzyna formulation as given in the equation (20):

Figure 10 -Figure 11 -
Figure 10 -Convergence curves at validation for lined tunnel

Figure 12 -
Figure 12 -Tunnel convergence after construction and after five years

Figure 13 -
Figure 13 -a) Node location b) Node evolution over time

Figure 14 -
Figure14-Sections through the experimental Kielder tunnel.Freeman[13] + The term t' is important to include the effects of aging, allowing a correct response of the material.

Table 1 -
Rockmass and concrete properties

Table 3 -
Rockmass and concrete properties of the experimental Kielder tunnel P. V. FIORE | D. B. MAGHOUS | A. CAMPOS FILHO Figure 15 -Results comparison over time