Acessibilidade / Reportar erro

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

Uma formulação tridimensional de elementos finitos para modelar um túnel revestido com concreto projetado e pré-moldado

ABSTRACT

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.

Keywords:
finite elements; tunnel; viscoelasticity; viscoplasticity; ANSYS

RESUMO

Este artigo apresenta um procedimento para a simulação numérica por elementos finitos tridimensionais do processo de escavação de um túnel com revestimento em concreto. O comportamento do maciço rochoso é representado pela lei viscoplástica de Perzyna. O concreto, projetado ou pré-moldado, é modelado como um material viscoelástico através das cadeias de Maxwell e Kelvin. As simulações por elementos finitos são realizadas pela incorporação de sub-rotinas dos modelos viscoelásticos para o concreto no programa ANSYS. O método para simular a escavação do túnel é pela ativação e desativação de elementos em etapas sucessivas. Na primeira parte do artigo duas validações do modelo são apresentadas. A solução analítica e a obtida pelo programa ANSYS, para a deformação alcançada após a estabilização, são comparadas para um túnel sem revestimento. Os resultados para um túnel revestido, considerando o maciço viscoplástico e o revestimento elástico, são comparados com os valores determinados pelo programa GEOMEC91. Na segunda parte, o mesmo túnel é analisado, considerando o revestimento de concreto viscoelástico, conforme os modelos das cadeias de Maxwell e Kelvin. Finalmente, é modelado o túnel experimental de Kielder, para o qual dados medidos in situ estão disponíveis.

Palavras-chave:
elementos finitos; túneis; viscoelasticidade; viscoplasticidade; ANSYS

1. 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][1] Ward WH, Tedd P, Berry NSM. "The Kielder experimental tunnel: final results." Géotechnique 1983; vol. 33, no. 3, pp. 275-291. , which in situ measurement data is available.

2. 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):

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):

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

If σ is outside or in the border of the elastic domain then F(σ)<0 and:

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):

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][2] Zienkiewicz OC, Cormeau I. "Visco-Plasticity and Creep in Elastic Solids - A Unified Numerical Solution Approach." International Journal for Numerical Methods in Engineering 1974; vol. VIII, pp. 821-845. :

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 F 0 is the reference stress.

3. The long term behavior of concrete structures

3.1 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.

3.2 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][3] Fédération Internationale du Béton, "Fib Model Code 2010. Final Draft" 2012; vol. 1, no. Bulletins 65, p. 350., 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):

where:

εci (t ') is the initial strain at loading; ε cc (t) is the creep strain at time t > t'; ε cs (t) is the shrinkage strain; and ε cT (t) is the thermal strain.

3.3 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][4] Creus GJ. Viscoelasticity - Basic Theory and Applications to Concrete Structures, Heidelberg: Springer-Verlag Berlin 1986; p. 169..

In order to simulate numerically the effects of the creep and aging of the concrete structures, Bazant [5][5] Bazant ZP. Mathematical Modeling of Creep and Shrinkage of Concrete, John Wiley & Sons, 1988. 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:

Where J(t,t') is the creep function, that represents strain in time t, caused by a constant unitary stress applied in the interval t-t'; Ec (t') is the elasticity modulus at time t'; φ (t,t') is the creep coefficient; Ec is the elasticity modulus of the concrete in the age of 28 days.

The term t' is important to include the effects of aging, allowing a correct response of the material.

3.4 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][3] Fédération Internationale du Béton, "Fib Model Code 2010. Final Draft" 2012; vol. 1, no. Bulletins 65, p. 350. is given by equation (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.

Figure 1
a) Maxwell model; b)Kelvin model

4. Viscoelastic formulation: basic rheological models

In accordance to Creus [4][4] Creus GJ. Viscoelasticity - Basic Theory and Applications to Concrete Structures, Heidelberg: Springer-Verlag Berlin 1986; p. 169. 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.

Figure 2
Generalized Maxwell chain model

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:

t5 →∞ , to simulate a coefficient of viscosity tending to infinity η5 →∞ 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][5] Bazant ZP. Mathematical Modeling of Creep and Shrinkage of Concrete, John Wiley & Sons, 1988. . 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.28 days. 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.

4.2 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.

Figure 3
Generalized Kelvin chain model

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][3] Fédération Internationale du Béton, "Fib Model Code 2010. Final Draft" 2012; vol. 1, no. Bulletins 65, p. 350.. These facts lead to a major numerical simplification and make the Kelvin model more advantageous than the Maxwell model [7][7] Bazant Z, Prasannan S. "Solidification theory for aging creep." Cement and Concrete Research 1988; vol. 18, pp. 923-932. . In accordance to Bazant [5][5] Bazant ZP. Mathematical Modeling of Creep and Shrinkage of Concrete, John Wiley & Sons, 1988. , it is advantageous to use progressively increasing time steps Δtr .

Eight loading ages and thirty discrete times are used for representing properly the creep functions at Kelvin model. From these creep curves calculated with discrete points, the values of E μ for each element of the chain are given from the least squares fitting method. 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).

of the time steps until the values of time greater than the shorter retardation time, keeping stability and a good accuracy.

4.3 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∙fcm (fcm 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.

4.4 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][3] Fédération Internationale du Béton, "Fib Model Code 2010. Final Draft" 2012; vol. 1, no. Bulletins 65, p. 350. 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.

Figure 4
Creep function J for t'=2.8 days by fib 2010 and computational models

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).

The use of an exponential algorithm allows the gradually growth

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

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.

5. 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.

6. Validation of FEM ANSYS code

6.1 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][8] Corbetta F. "Calculus Analytiques et Numériques de Tunnels Profonds," Thèse de Doctorat de l'École Nationale Supérieure des Mines de Paris, 1990. presented an analytical solution for convergence, in elasticity, at a section far from the tunnel face. It is given by the equation 18.

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][9] Bernaud D. "Tunnels profonds dans les milieux viscoplastiques: approaches experimentale et numérique," Thèse de Doctorat de l'École Nationale des Ponts et Chaussées, Paris, France, 1991. ) and the NIM (1D solution, [10][10] Bernaud D, Rousset G. "The New Implicit Method for Tunnel Analysis," International Journal for Numerical and Analythical Methods in Geomechanics 1996; no. 20, pp. 673-690. ). The results present a very good agreement.

6.2 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).

Figure 6
Usermat flowchart for Maxwell chain model

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][8] Corbetta F. "Calculus Analytiques et Numériques de Tunnels Profonds," Thèse de Doctorat de l'École Nationale Supérieure des Mines de Paris, 1990. at equation 19.

Figure 7
Usermat flowchart for Kelvin chain model

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.

7. 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; distance

Figure 8
Convergence curves at validation

Figure 9
Convergence curves at validation for unlined tunnel

Figure 10
Convergence curves at validation for lined tunnel

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][11] AFTES, "AFTES Recommendations," in Design of Sprayed Concrete for Underground Support, www.aftes.asso.fr, Ed., Association Française des Tunnels et de l'Espace Souterrain 2000; p. 36.
www.aftes.asso.fr...
.

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.

7.1 The model in the computational ANSYS code

The tunnel construction simulation is performed with 37 excavation steps (tunnel with circular section), with the radius R m=1 , in an isotropic rockmass, with a distance between the tunnel face and

Figure 11
3D model mesh

the lining and a step of excavation 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 P = 4 MPa . The concrete lining has a

constant thickness equal to

Due to the 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):

where σ is the stress, m is the parameter related with the strain rate of the material, g is the parameter related to the viscosity of the material and σ o is the yield stress. These parameters are shown in Table 1. is the rate of equivalent strain,

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][12] ANSYS, "Mechanical APDL Structural Analysis Guide," ANSYS, Inc., Version 14.5, 2012..

7.2 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

Table 1
Rockmass and concrete properties

Figure 12
Tunnel convergence after construction and after five years

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.

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

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.

7.3 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).

7.4 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

Table 2
Results relation between different models and time

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.

8. 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][1] Ward WH, Tedd P, Berry NSM. "The Kielder experimental tunnel: final results." Géotechnique 1983; vol. 33, no. 3, pp. 275-291. 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

Figure 14
Sections through the experimental Kielder tunnel. Freeman [13][13] Freeman BTJ. "The behaviour of fully-bonded rock bolts in the Kielder Experimental tunnel." Tunnels & Tunneling, June 1978, pp. 37-40.

Table 3
Rockmass and concrete properties of the experimental Kielder tunnel

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][1] Ward WH, Tedd P, Berry NSM. "The Kielder experimental tunnel: final results." Géotechnique 1983; vol. 33, no. 3, pp. 275-291. and Freeman [13][13] Freeman BTJ. "The behaviour of fully-bonded rock bolts in the Kielder Experimental tunnel." Tunnels & Tunneling, June 1978, pp. 37-40.. 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.

Figure 15
Results comparison over time

9. 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][7] Bazant Z, Prasannan S. "Solidification theory for aging creep." Cement and Concrete Research 1988; vol. 18, pp. 923-932. have been tested with the results of the fib Model Code 2010 [3][3] Fédération Internationale du Béton, "Fib Model Code 2010. Final Draft" 2012; vol. 1, no. Bulletins 65, p. 350.. 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 ANSYS 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.

10. References

  • [1]
    Ward WH, Tedd P, Berry NSM. "The Kielder experimental tunnel: final results." Géotechnique 1983; vol. 33, no. 3, pp. 275-291.
  • [2]
    Zienkiewicz OC, Cormeau I. "Visco-Plasticity and Creep in Elastic Solids - A Unified Numerical Solution Approach." International Journal for Numerical Methods in Engineering 1974; vol. VIII, pp. 821-845.
  • [3]
    Fédération Internationale du Béton, "Fib Model Code 2010. Final Draft" 2012; vol. 1, no. Bulletins 65, p. 350.
  • [4]
    Creus GJ. Viscoelasticity - Basic Theory and Applications to Concrete Structures, Heidelberg: Springer-Verlag Berlin 1986; p. 169.
  • [5]
    Bazant ZP. Mathematical Modeling of Creep and Shrinkage of Concrete, John Wiley & Sons, 1988.
  • [6]
    Pande GN, Owen DRJ, Zienkiewicz OC. "Overlay models in time-dependent non-linear material analysis." Computers and Structures 1977; vol. 7, no. 3, pp. 435-443.
  • [7]
    Bazant Z, Prasannan S. "Solidification theory for aging creep." Cement and Concrete Research 1988; vol. 18, pp. 923-932.
  • [8]
    Corbetta F. "Calculus Analytiques et Numériques de Tunnels Profonds," Thèse de Doctorat de l'École Nationale Supérieure des Mines de Paris, 1990.
  • [9]
    Bernaud D. "Tunnels profonds dans les milieux viscoplastiques: approaches experimentale et numérique," Thèse de Doctorat de l'École Nationale des Ponts et Chaussées, Paris, France, 1991.
  • [10]
    Bernaud D, Rousset G. "The New Implicit Method for Tunnel Analysis," International Journal for Numerical and Analythical Methods in Geomechanics 1996; no. 20, pp. 673-690.
  • [11]
    AFTES, "AFTES Recommendations," in Design of Sprayed Concrete for Underground Support, www.aftes.asso.fr, Ed., Association Française des Tunnels et de l'Espace Souterrain 2000; p. 36.
    » www.aftes.asso.fr
  • [12]
    ANSYS, "Mechanical APDL Structural Analysis Guide," ANSYS, Inc., Version 14.5, 2012.
  • [13]
    Freeman BTJ. "The behaviour of fully-bonded rock bolts in the Kielder Experimental tunnel." Tunnels & Tunneling, June 1978, pp. 37-40.

Publication Dates

  • Publication in this collection
    June 2016

History

  • Received
    15 Dec 2015
  • Accepted
    10 Mar 2016
IBRACON - Instituto Brasileiro do Concreto Instituto Brasileiro do Concreto (IBRACON), Av. Queiroz Filho, nº 1700 sala 407/408 Torre D, Villa Lobos Office Park, CEP 05319-000, São Paulo, SP - Brasil, Tel. (55 11) 3735-0202, Fax: (55 11) 3733-2190 - São Paulo - SP - Brazil
E-mail: arlene@ibracon.org.br