Open access

Finite Element Analysis for the Problem of Tunnel Excavation Successive Phases and Lining Mounting

Written By

Simona Roatesi

Submitted: 16 November 2011 Published: 10 October 2012

DOI: 10.5772/46168

From the Edited Volume

Finite Element Analysis - New Trends and Developments

Edited by Farzad Ebrahimi

Chapter metrics overview

3,561 Chapter Downloads

View Full Metrics

1. Introduction

Numerical application to geomechanics is an area of research that has grown rapidly since its origins in the late 1960s. This growth led to new methodologies and analysis approaches that are nowadays commonly employed in geotechnical, mining or petroleum engineering practice.

The circular tunnel boundary problem is frequently encountered in rock and soil mechanics, geotechnics and generally in mining and petroleum engineering. Consequently, there is a great interest in solving boundary problems involving the excavation of an underground opening. Moreover, time dependent behavior of geomaterials benefits by a special attention in constitutive and numerical approach. Regarding the constitutive modeling and analytical approach [1] concerns both with vertical and horizontal underground openings excavated in an elasto-viscoplastic rock mass, governed by the constitutive laws which go by his name and for which analytical solution for the displacement and viscoplastic strain are derived; [2] provides analytical results for circular and non-circular openings with thermal effects as well; [3-5] present results in approaching the problem both analytically (convergence-confinement method) and numerically by FEM.

The viscoplastic approach in FEM has been considered by many authors, in terms either of numerical stability of schemes used in viscoplasticity, see [6-8], or practical applications using FEM solutions in different areas. For instance, tunneling using FEM for viscoplastic materials has been investigated widely, as well, e.g. [9] approaches the problem with outstanding results in computational geomechanics, with special reference to earthquake engineering, in numerical modeling of dynamic soil and pore fluid interaction and earthquake-induced liquefaction and multiphase pollutant transport in partially saturated porous media; Zienkiewicz provides outstanding numerically approaches in a wide range of problems, e.g., [10] is devoted to computational geomechanics; [11] focuses a great interest in the multi-disciplinarily aspect of the problem, taking into account also adjacent phenomenon occurring at the rock-support interaction; [12] deals with ductile damage and fracture FE modeling of viscoplastic voided materials for high strain rate problems,[13] provides a finite-strain viscoplastic law coupled with anisotropic damage both theoretical and numerical approach, [14] develops a FE procedure to model the tunnel installation and the liner to predict the likely extent of damage to surface structures caused by nearby shallow tunneling, [15] deals with FE modelling of excavation and advancement process of a shield tunnelling machine, [16] develops a FE micromechanical-based model for hydro-mechanical coupling for tunnelling application, [17] applies a model based on plastic damage evolution and permeability to excavation-disturbed zone simulation of the mudstone shield tunnel; [18] analyses tunnel depth effect on the stress and strain state around the tunnel; [19], [20] studies the face tunnel influence in the analysis of a circular tunnel with a time-dependent behaviour, etc.

This chapter deals with a FE code implementation of an elasto-viscoplastic constitutive law. The numerical calculations are performed with a finite element code called CESAR made in LCPC-Paris [21]. The viscoplastic module is coded and implemented in the FE code CESAR by the author. The validation of the numerical code is performed in different steps: with the boundary problem solution of the triaxial laboratory tests, with analytical solution of creep step, and of supported and unsupported underground openings in viscoplastic rock mass. This chapter is focused on further complex applications, such as the tunnel excavation successive phases and lining mounting which are approached using the viscoplastic constitutive module.


2. The elasto-viscoplastic constitutive equation

The constitutive law implemented in finite element code is proposed by Cristescu. The hypotheses for the constitutive equation formulation are (see [1]):

The material is considered homogeneous and isotropic. Thus, the constitutive functions depend only on the invariants of the stress and strain tensors. The stress tensor and the strain tensor are denotedσandε, respectively and the stress tensor principal components are denoted as:σ1,σ2,σ3,ε1,ε2,ε3. Among the stress invariants, those with an important physical meaning are:

the mean stress:σ=(σ1+σ2+σ3)3;

the equivalent stress σ¯2=σ12+σ22+σ32σ1σ2σ2σ3σ3σ1(or the octahedral shear stress τ=23σ¯=23IIσwith IIσbeing the second invariant of the stress deviator).

The displacements and the rotations are assumed small, so thatε˙=ε˙E+ε˙I, where ε˙Eand ε˙Ι being the elastic strain rate and the irreversible strain rate respectively.

The component of the elastic strain rate satisfies the Hooke law


with K, G being the bulk and shear modulus, respectively andIbeing the identity tensor.

The component of the irreversible strain rate satisfies:


where: k represents the viscosity coefficient, that can depend on the stress state and the strain state, and probably on a damage parameter d describing the history of the micro cracking the rock was subjected to, and the bracket <> represents the positive part of respective function:


The irreversible stress work is used as a hardening parameter or internal state variable, split into volumetric and deviatoric parts and it is given by the quantity:


We introduce the damage parameterd, see [1], defined by:


describing the energy released by micro-cracking during the entire dilatancy period. In (2)tmax represents the time for which WvIis maximum. The failure threshold is considered to be the total energy released by micro cracking during the entire dilatancy process and it is characterized by the following parameter (constant):

H(σ) represents the loading function, generally a function of stress tensorσ, withH(σ,σ¯)=WI(t)the creep stabilization boundary equation, the function H depending on the two stress invariants mentioned above.F(σ) represents a viscoplastic potential, that establishes the orientation ofε˙Ι.

The initial yield stress of the material is zero or very close to zero.

The applicable domain for the constitutive equation is considered for compressive stresses (positive) and bounded by the failure surface which may be incorporated in the constitutive equation.

We consider therefore, the following constitutive equation:


For instance, for the model describing the Borod coal behavior (see [1]), the constitutive functions and material constants used are:

H(σ,σ¯):=a0(σ¯/σ*)2a2σ/σ*+a1+b0σ¯/σ*+{c0sin(ωσ/σ*+φ)+c1if0σσ0c0+c1ifσ0σ (the yield surface)

where a0 = 7.65 × 10−4 MPa; a1 = 0.55; a2 = 8.159 × 10−3; b0 = 0.001 MPa; c0 = 4.957 × 10−4 MPa; c1 = 4.8955 × 10−4 MPa; ω =171.9270;σ0=1,0996MPa; k =6 × 10−6; s-1 (the viscosity coefficient); df = 4.48 × 103 Jm−3; σ*=1 MPa and the elastic constants are: E = 798.385 MPa, υ = 0,327.

The relation (5) shows that the viscoplastic potential equals the loading function, fact that determines an associated flow rule.

For the model describing the saturated sand behaviour, the constitutive functions and material constants are:

H(σ,σ¯):=a(σ¯)7(σσ¯33)5+bσ¯+cσ,(the yield surface)E9
F(σ,σ¯):=2h1σ322f+α+2fh1(23σ322f+α+2(1+α3)32(2f+α)3)+(h1σ¯[(1+α3)σ¯]12(2f+α)322fh1[(1+α3)σ¯]32(2f+α)52)ln[(2f+α)σ]12+[(1+α3)σ¯]12[(2f+α)]12[(1+α3)σ¯]12+g0σ¯+g1σ¯3(the plastic potential)E10

wherea=4.834×10-7 (kPa)1,b=1.33×103 ,c=1.058×103 ,h1=2.34×10 (kPa)12 ,g0=0.005 , g1=0.62×10-6(kPa)-2, f=0.562, α=1.34,k=3.6×10-6s1 (the viscosity coefficient) and the elastic constants are: E = 205300 kPa, υ = 0.33.

For this model, as the viscoplastic potential differs to the loading function, determines a non-associated flow rule.


3. The numerical integration of the elasto-viscoplastic equation

In this paragraph the way to implement the elasto-viscoplastic law presented in previous paragraph in finite element code is described. The interface with the program is performed by the calling of the subroutine that carries out the numerical integration of the proposed constitutive law, at a Gauss integration point level. The implicit form for the constitutive law is used and two solving methods, namely: Euler semi implicit method (θscheme) and Runge - Kutta method of fourth order.

The numerical integration of equation (4) is performed alternatively by two methods:

Semi implicit Euler method (θscheme) in which the evaluation of the viscoplastic strain increment is done with the rule: Δεvp=(ε˙1vp(1θ)+ ε˙2vpθ)Δt,θ[0,1]whereε˙1vp=f(σ,W), ε˙2vp=f(σ+Δσ1,W+ΔW) and Δσ1,ΔW1are the stress increment and associated irreversible work increment calculated with the standard scheme, corresponding to the viscoplastic strain incrementΔεvp1; σ,Ware the known values of stress and irreversible work at the beginning of the step and functionf(σ,W)=k(σ,d)1WI(t)H(σ,σ¯)Fσ.

Runge-Kutta method performs more evaluations of the function inside each time step dt, propagating thus on an interval a solution that is a combination of a few Euler type steps (each of it implies an evaluation of function f) and further using the obtained information to match an expansion in Taylor series of higher orders.

In particular, Runge-Kutta method of fourth order employed to integrate the proposed constitutive equation uses four evaluations of function f for the considered time step:

withΔσ1,Δσ2,Δσ3,ΔW1,ΔW2,ΔW3the intermediate evaluations of the stress increment, and the irreversible work increment respectively.

Runge-Kutta method of fourth order is used for the numerically integration of equation (4) to evaluate the viscoplastic strain increment and it is also used to evaluate the irreversible work (1) from the evolution differential equation of the hardening parameterΔWi.


4. The numerical solution of nonlinear problems

The classical solution algorithms used in finite element method are displacement type, incremental and iterative ones that often present some convergence difficulties related to the applied loading and the searching of the limit loading.

Generally, the solution algorithms of nonlinear problems will depend on the specific problem under consideration. At the same time, the algorithms are made up of two levels: one global level, of the general scheme of solving, that enables the calculation of the displacement field at the nodal points of the discretized structure; and a local level, of the integration scheme of a nonlinear constitutive equation, which enables at a point, as starting with the stress and strain history at that point, a new state of stress to be calculated.

The two schemes have a great reciprocal influence upon the solving process. Concerning the integration of the constitutive law, it may be done through explicit or implicit schemes.

The treating of a nonlinear problem with finite element method leads to the solving of a system of equations that may be put in the form:


where: u represents the displacement vector in the nodal points of the structure under consideration, R(u) is the nodal forces vector corresponding to the stresses at moment t, P designates the total loading applied to the structure and λ(t)represents the loading factor applied at moment t.

The solution of the system of equations (6) is in fact the pair (u,λ(t)) associated to the displacement response of the structure at the loading which it is subjected to. Generally it is impossible to obtain a solution with a direct solution technique, so an iterative process has to be adopted. Often, the used iterative process is based by the linearization of the nonlinear equations to be solved.

The iterative process assumes that the constitutive equation can enable the calculation of the exact value of the inequilibrium (residual) with respect to the only unknown entity, which is the present displacement u1. For this reason it is necessary to perform an incremental loading of the structure.

The algorithm for solving a static nonlinear problem with the initial stress method used in the framework of the finite element program is:


j = 1loop on the increments

The incrementation of the loading P

i = 1loop on the iterations (We assumed i state known and i+1 state has to be calculated).

Setting up the termΔFi=PijBtσidv that has to be calculated as accurate as possible.

The calculation of the global stiffness matrixK=BtDBdv, where the matrix D needs not be quite accurate and can be:

elastic matrix (the initial stress method)

the constitutive matrix Devp(the tangent stiffness method)

The solving of the algebraic system KΔui=ΔFi(nonlinear if K=K(u))

Δui=K1ΔFisolving the system of equationsui+1=ui+Δuiupdating displacementΔεi=BΔu1displacement-strain relationsσi+1=σi+Δσicalculation of the incremental stresses with the constitutive law; Δσcan be calculated as:


viscoplasticity: *Δσi=De(ΔεiΔεivp), withΔεivp=Λf/σ,proportional with the derivative of the viscoplastic potential.

*Δσi=DevpΔεi (explicit form).

From the law Δεi=f(Δσi)(implicit form), the numerical integration may be done with θ scheme, with Runge-Kutta method, the consistent matrix, etc.

R(ui+1)=Btσi+1dvcalculation of the residue φi+1=PjR(ui+1)equilibrium check

Convergence test:

no i = i+1 (next iteration)

yes j = j+1 (next step)

where u represents the vector of nodal displacements, v is the vector of out of balance nodal forces, B designates the matrix of the shape functions for strain,Pj represents the vector of the known nodal applied forces, R(ui)is the vector of the equivalent nodal forces, due to the stressesσi. These nodal forces are consistent with the current value of unknown u.

Returning again on the fact that it can be also deduced from the algorithm above, namely that two calculation levels can be distinguished:

local level: the application of the constitutive law for a material point (in fact, Gauss integration point in finite element analysis) for the calculation of the constitutive matrix D,(Δσ=DΔε) and the stressesσi+1.

global level: the application of the iterative process upon the vector of nodal displacements.

In order to ensure the convergence of the iterative process, the non-equilibrium of the structure (the residual), the displacement variation and the work variation done during the iteration have to be tested at every iteration. So, the testing is done upon the following quantities:

testing the residue |ψ(ui)||ψ(u0)|εR;

testing the displacement|Δui||ui|εD;

testing the work |ψ(ui)T.Δui||ψ(u0)T.Δu0|εW.

It has to be mentioned the fact that an iteration has no physical meaning: in fact, at the beginning of the iteration, the equilibrium is satisfied, but, at the end of one iteration, the constitutive law is satisfied in preference to the equilibrium. Therefore, only at the end of an incremental stage, the solution can converge to a physical sense for the studied structure.


5. Comparison of the numerical and analytical solution for the creep step

In order to test the subroutine which performs the numerical integration of the elasto-viscoplastic constitutive equation, a comparison of the numerical solution with the analytical formula for the creep step is performed. For this purpose, an analytical formula is used for the calculation of the strain rate in the case of the application of a number of successive steps at constant stress, namely: it is assumed that at momentt0the stress state is increased suddenly to the value of σ(t0)and it is kept constant.

5.1. Determination of the analytical solution for creep step

To establish the formulae that describe the creep deformation, we will write firstly the formula which supplies the variation of WI(t)when all stress components are constant; it is easily obtained by integration of the constitutive law (4) (see [1]):


whereσ(t)=σ(t0)=constant and WI(t0)is the initial value of WIfort=t0, denoted furtherWIP.

This is the relation that describes the variation in time of WI(t)under a constant stress. The variation in time of WI(t)is longer or shorter depending especially on the value of the viscosity parameter k.

Introducing the left hand side of the relation (7) into the constitutive equation (4), we get the total strain under constant stress:


using the initial conditions:


where σ(t)=σ(t0) is the initial stress state, reached instantaneously and taken with respect to the state at t0 when the loading is applied (it represents the elastic response). It is observed that the value of irreversible strain rate during the creep is governed by the expression 1WI(t)H(σ,σ¯).

If we wish to have a stress path made up of small stress incrementsΔσ, followed by time intervals of constant stress, we get for the strain increments some formulae similar with (8), namely:


whereσ1andσ2being the principal stresses and where during the time intervalt(t0,t] the state of stress is constant and WI(t0) is calculated at momentt0, just before of a new stress increment Δσthat occurs att0.

5.2. Comparison of the numerical results with the analytical solution

The numerical solution is obtained for the same loading path as in the calculation of the analytical solution. Since the calculation of the analytical solution uses the hypothesis of an instantaneous loading, for the numerical solution the loading is applied in a very small time intervalΔt, namely 1010s. Similar to the calculation of the analytical solution in the hypothesis of a constant stress, for the numerical solution, at every time step, null loading increment is applied.

To perform a comparison between the analytical solution and the numerical results obtained using the Euler method (θscheme) and Runge-Kutta method of fourth order, we present as example, one step loading path, with a null initial state of strain and an initial stress state of σ10=0.7kPa,σ20=0.1kPa,σ30=0.1kPa.The numerical results obtained using Runge-Kutta method of fourth order shows the superiority of this method.

For instance, in the case of Δt = 0.01 s, the analytical solution for the strain components gives:

ε1analytic=  3.08839 × 10-6,ε2analytic =-7.99371 × 10-7,ε3analytic=-7.99371 × 10-7;E18

while the same components using θ scheme are:

ε1Euler=  3.08841 × 10-6,ε2Euler =-7.99383 × 10-7,ε3Euler =-7.99383 × 10-7;E19

and using Runge-Kutta of fourth order methods, are:

ε1RK=  3.08840 × 10-6,ε2RK =-7.99378 × 10-7,ε3RK =-7.99378 × 10-7.

6. Comparison of the numerical solution of the triaxial test boundary problem and the experimental data

The numerical solution is performed using the elasto-viscoplastic constitutive law presented above and represents the triaxial test for a cylindrical sample of saturated sand(σ1increasing,σ2=σ3) under axi-symmetric hypothesis, see [20].

A quarter of a sample was considered because of the symmetries, requiring only one quadrilateral finite element with 8 nodal points. So, the mesh is a very simple one, like in Figure 1, the width being denoted with a and the length l. In this case, we consider a = l = 1.

Figure 1.

The mesh for the simulation of the triaxial test for a cylindrical sample

The stress state in element being uniform, the 8 nodal points of the element present the same state of stress and strain. For validation we deal only with 8th nodal point element, which has the advantage that since its height being equal with unity, its vertical displacement is equal with the longitudinal strain, and its horizontal displacement is equal with the radial strain.

The initial boundary problem is then:

(inner face)z=0andr[0,a]ur(t)=0,σrz(t)=0,t0.(left lateral surface)z[0,l]andr=0ur(t)=0,σrz(t)=0,t0.(right lateral surface)z[0,l]andr=aσrr(t)=f(t),σrz(t)=0,t0.(outer face)z=landr[0,a]σzz(t)=g(t),σrz(t)=0,t0.E21

We consider for the functions f and g some particular forms that can simulate the triaxial test in two stages, namely:

the hydrostatic stage when both lateral pressure and the vertical one are increased (in steps for this test) until a certain value is reached.

the deviatoric stage, when the lateral pressure remains constant at while the vertical one is increased (in steps as well) until failure.

The following functions are used for the two stages:

f(t)={σhtTfor0tT(the hydrostatic stage)σhforTt(the deviatoric stage)g(t)={σhtTfor0tT(the hydrostatic stage)σh+Σh(tTτ)forTt(the deviatoric stage)E22

with a linear variation during the hydrostatic stage.Σh,τ,T are constants that characterize the loading.

We are going to consider as well a quadratic variation with respect to the time of the loading through the functions:

f(t)={σht2T2for0tT(the hydrostatic stage)σhforTt(the deviatoric stage)g(t)={σht2T2for0tT(the hydrostatic stage)σh+Σh(tTτ)2forTt(the deviatoric stage)E23

Let us consider now a hydrostatic stage of a classical triaxial test (we consider for the functions f, g, f', g' the upper branch of the function form defined above). We will try to find again some of the important constitutive features of the material with the numerical solution.

The comparison of the results for a linear incremental in time loading (considering the functions f, respective g) and a quadratic incremental in time loading (considering the functions f', respective g') are made for two different duration of creep time step, 10s and 30s respectively. In both cases, the volumetric strain corresponding to the quadratic loading is bigger than in the case of the linear one in time, the difference being more marked at smaller loading rates.

A comparison with respect to the loading rates for 7 loading steps, and 14 loading steps hydrostatic test respectively is made too. In both cases one can observe grater strains in the case of smaller loading rates, the strains being as great as the test performed with more loading steps.

For the deviatoric stage of the classical triaxial test (the lower branches of the functions f, g, f', g') in which the stress state that was reached in the hydrostatic stage is maintained constant and only the vertical component of the stress is increased, a comparison between the numerical results and the experimental data is achieved for saturated sand (see [22]).

There are two sets of tests with the confining pressure of 14.7 kPa and 29.4 kPa respectively. Figures 2a,b represents the experimental data, while Figures 3 a, b presents the numerical results for the two tests mentioned above, respectively.

Figure 2.

a), b). Experimental data for 2 tests with confining pressure of 14.7 kPa, 29.4 kPa respectively

Figure 3.

a),b). Numerical solutions for 2 tests with confining pressure of 14.7 kPa,29.4 kPa, respectively


7. Comparison of numerical and analytical solution for underground openings

In this paragraph, a comparison between the numerical results and an available analytical solution for the problem of underground openings in viscoplastic rock mass is performed and it represents the next step of validation of the numerical code. The approach presented here concerns the applications of supported and unsupported underground openings in viscoplastic rock mass with the assumption of constant primary stress in the whole domain.

7.1. Problem formulation

The proposed boundary problem is as follows: it is assumed that the rock mass is an infinite body in which circular opening is made. Therefore only plane strain condition is considered.

Assuming further that the underground opening is at a certain depth where the horizontal and vertical components of the primary (initial) stressσhandσvare known and, generally, distinct. It is also assumed that in the neighbourhood of the opening these components are constant and equal with their corresponding value for the opening axis depth. The influence of the ground surface is neglected, so that the opening is imagined as a cylindrical cavity of infinite length, excavated in an infinite space.

Cylindrical coordinate system is chosen for convenience with axis Oz being the symmetry axis of the opening (Figure 4).

Figure 4.

The domain of the boundary value problem

A pressure loading p(t) is considered on the surface r = a of the opening, that could be due to different causes. Therefore, the boundary conditions are:


with Ox and Oy the horizontal and the vertical axis respectively.

The conditions (9) can be written in the cylindrical coordinates as follows (see [1]):


The superscripts S, P, R mean the secondary, primary and relative components of the stress, displacement and strain respectively.

Further, the fundamental equations of the problem are presented.

The equilibrium equations written in the relative components are:


The components of the small strains:


and the compatibility equations in cylindrical coordinates has to be added, as well:


Due to the plane strain hypothesis, we have z=0,uzR=0,and then, according to the components of the small strains, zero values are founded for the relative strains εzzR=εrzR=εθzR=0.

A viscoplastic constitutive equation as presented in the previous paragraph is considered.

7.2. Analytical solution

It is assumed that the opening (or only a part of it) is excavated in a very short time interval t(0,t0)and then it is exploited in a much longer time intervalt(t0,t1). It is also assumed that during the first time interval the response of the rock is "instantaneous" and during the second one different time effects are possible, such as: creep and a slow variation in time of stress.

Consequently, if the tunnel is excavated att0, then immediately after excavation, the stress, the strain and the displacement of the rock are given by the elastic solution. Concerning the second interval (t0,t1) a possible solution will be present, obtained under a number of assumptions that would simplify a lot the analytical solution thus making it amenable to analysis.

The main hypothesis is that the same distribution of stresses for the elasto-viscoplastic model is assumed as in the elastic model: by considering that in the interval (t0,t1) the stress componentsσrr,σθθ,σrθremain constant and equal with those given by the elastic solution (see [1]).

Due to the fact of the plane strain hypothesis εzzR=0 it has to be accepted thatσzzvaries during this interval and therefore it satisfies the differential equation:


This differential equation may be integrated numerically using the initial conditionsσzz=σh. It was found that a very small variation in time (negligible) ofσzzwas exhibited and a very fast stabilization of this variation. Therefore, it can be assumed that all stress components could be constant during the rock creep around the opening.

In order to arrive at the expression that describes the creep strain, the equation (7) that supplies the variation of WI(t)when all stress components are constant is used.

The time taken for WI(t)to reach the asymptotic value depends mainly on the value of the constitutive parameter k, the viscosity. A reasonable correct value for k can be obtained only observing the convergence of a tunnel wall in time.

Then one can obtain, similarly as in the previous paragraph, the formulae (10) for the strain variation by integrating the constitutive equation (4) taking into account (7) and using the initial conditions (9).

In the case of σh=σv one can deduce a formula for the wall opening convergence as:

 uR  0.E30

Figure 5.

Domain used in numerical formulation

7.3. The numerical solution

In the case of numerical approach, the domain to be discretized and the boundary conditions, see Figure 5 and the FEM mesh, see Figure 6, are:

Figure 6.

FEM mesh

In this case a quarter of the domain needs to be considered as the principal components of the primary stressσhandσvare assumed constant over the entire domain.

For this problem a mesh of 6 nodded triangles was used. The presence of the lining is simulated either by an internal pressure acting on the tunnel walls or by introducing some elements (8 nodded quadrilaterals) in contact with the tunnel walls, with the mechanical characteristics of the lining material (Figure 6).

In the following examples the lining was considered made of concrete with Young modulus E = 200000 kPa and Poisson coefficient υ = 0.3.

There were considered the cases whenσh=σv= 2000 kPa and the caseσhσv, namelyσh= 1500 kPa andσv= 2000 kPa. The tunnel radius is 1 m.

It is well known that the choice of the time step is quite important in viscoplasticity problems. This choice must be correlated with the viscosity parameter k as well, so that the value of k dt to ensure a good convergence of the numerical scheme.

In the model under consideration of saturated sand it was observed that values forkΔtwhich exceeds the magnitude order of103produce a divergence of the numeric calculation with the present method (the initial stress method and the initial stress method combined with different methods of acceleration). So, some stronger methods have to be implemented, in order to offer a faster convergence of the numerical calculus, e.g. the method of consistent tangential operator, backward Euler method, etc.).

The next Figures present other aspects of the numerical solution both in the case of an internal pressure acting on the tunnel wall, and for the case of concrete lining, that cannot be represented by the analytical solution.

7.3.1. Numerical solution for the case of internal pressure acting on the tunnel wall

The following Figures represents the numerical solution for the case of a tunnel subjected to a hydrostatic primary stress ofσh=σv= 2000 kPa and an internal pressure on the tunnel walls p = 1000 kPa which simulates the lining.

Figures 7a and 7b represent the evolution of the Euclidian norm of the viscoplastic strain |εvpij| for the first time step and after 130 time steps with one time step 10000 s, respectively.

Figure 7.

a),b). Contours of |εvpij|evolution of tunnel with an internal pressure in the case of a hydrostatic primary stress after the first time step and after the 130th time step, respectively

A tunnel driven in a rock mass with a non-hydrostatic primary stress of σh= 1500 kPa andσv= 2000 kPa, exhibits for instance, a development of a viscoplastic zone in the tunnel wall that increases from 1.13E-4 for |εvpij|after the first time step Figure 8 to 5.13E-3 for |εvpij|after 130 time steps.

The radial stress component σrrpresents a very slow variation in time and it is presented in Figure 9.

Figure 8.

Contours of |εvpij|after the first time step

Figure 9.

Contours of σrrafter the first time step

For the same boundary problem the evolution of the hoop stressσθθis designed in Figures 10a for the first time step and 10b for the 130th time step.

Figure 10.

a) Contours σθθafter the first time step,b) Contours σθθafter the 130th time step

7.3.2. Numerical solution for the case concrete lining

In this subparagraph, the numerical solution for a concrete lining as a distinct material and mesh elements group is presented.

For a case of a lined tunnel by a concrete lining with Young modulus E = 200000 kPa and Poisson coefficient υ = 0.3 for the same type of rock mass subjected to a primary stress ofσh= 1500 kPa and σv= 2000 kPa the diagrams of variation are presented as follows: the total displacement utotal=ux2+uy2in Figures 11a and 11b, and the norm of the viscoplastic strain|εvpij| in Figures 12a and 12b respectively.

The evolution of the equivalent stress σ¯and hoop stressσθθ is also studied. Both in the case of the equivalent stress and of the hoop stressσθθ, an increasing of the smaller value zone from the roof, and a decreasing of the big value zone from the tunnel wall are observed.

Figure 11.

a),b) Contours of utotalevolution in a supported tunnel excavated in a rock mass with non-hydrostatic primary stress, after the first time step and after 130th time step, respectively

Figure 12.

a),b) Contours of|εvpij|evolution in a supported tunnel excavated in a rock mass with non-hydrostatic primary stress, after the first time step and after the130th time step, respectively

7.4. Comparisons between the numerical and analytical solutions

Although the analytical solution was obtained under the assumptions of constant stress during the creep, it represents a good benchmark for the numerical calculation, supplying important information both quantitatively and qualitatively.

On the other hand, the numerical solution introduces certain truncation due to its discretisation errors. Nevertheless, the comparison between the two solutions could lead to some constructive conclusions for both solutions.

In Figure 13a the variation with respect to the distance in radial direction of the radial displacement in a tunnel excavated in condition of primary stress σh= σv= 2000 kPa, with an internal pressure on the tunnel wall p = 1000 kPa, at a chosen time, namely after 5 time steps ( one time step is considered 10000 ) is presented. A good agreement in the neighbourhood of the opening is observed, then the two solutions begin to differ, taking into account on one hand the assumptions of constant stress for the analytical solution and on the other hand the errors accumulated during the numerical process, due to the course mash at great distance. Figure 4b represents a comparison between the two displacements at two locations (the polar radius r = a and the polar angle θ = 0) at tunnel boundary, under the same condition as above, but with respect to time. One can observe that for small numbers of time steps, the numerical solution is superior, but for larger number of time steps (i.e. more than 10 time steps), the analytical solution becomes larger.

Figure 13.

a),b). Tunnel outline total displacement with respect to distance and time, respectively, for the numerical and analytical solutions

In the case of the viscoplastic strain, the analytical solution predicts greater values than the numerical one, even from the beginning of time analysis. That can be observed in Figure 14, for the case of a tunnel excavated in a rock mass with a non-hydrostatic primary stress ofσh= 1500 kPa,σv = 2000 kPa, and the internal pressure p = 1000 kPa.

Figure 14.

The time evolution of viscoplastic strain of tunnel walls for numerical and analytical solutions

7.5. Final remarks

Concerning the analytical solution for the problem under consideration the main remark is that it is difficult to obtain the analytical solution and it can be obtained using some (apparently) severe assumptions such as the same stress distribution is used for linear elastic and elasto-viscoplastic condition and has to be invariant in time. On contrary, the numerical solution via Finite Element does not impose such restriction and the results are exhibiting a much slower variation with time during creep. Thus this could attenuate the severity of the assumption of the analytical solution, assumption that could not be considered totally unrealistic.

Since the analytical solution is the exact solution of the governing equations under some restrictions, it can be used to benchmark the FE analysis before putting it in general applications.

Some remarks can also be drawn taking into account the features exhibited by the numerical solution.

Comparison between the numerical solutions obtained in the case of hydrostatic and nonhydrostatic primary stress respectively, shows that a much greater increase in time for the viscoplastic strain in the nonhydrostatic case, and a larger zone in the rock mass is exhibiting viscoplastic behaviour too.

A very slower decrease in time of the radial stress is observed, while the hoop stress has a tendency to increase in time.

In a more realistic case of the lining being simulated as a distinct zone of material with specific mechanical characteristics (concrete in the cases under consideration), instead of being applied as an internal pressure on the tunnel walls, some remarks could be made too. With the concrete lining modelled as elements, the displacement is greatly reduced and the viscoplastic strain is much less than the case of the internal pressure being applied on the tunnel walls. The amount accumulated in time is 5 times smaller than the other ones.

All these features are in good accordance with the practical observation, so, more involved boundary problems could be envisaged to be solved with this code.


8. Numerical solution of the three phases tunnel excavation and lining mounting problem

A finite element solution for the problem of a circular tunnel excavated in a homogeneous isotropic elasto-viscoplastic rock mass is presented. The numerical model consists of the successive phases of the excavation and support mounting, emphasizing the role of two important factors of the analysis, namely the time and the tunnel face influence and taking into account the 3D aspect of the problem.

The behavior of rock mass is considered viscoplastic, while the concrete lining is elastic. A possible behavior of sliding at the rock-support interface, which requests some additional contact elements of the mesh, is neglected. The rock mass is homogeneous, for the simplicity of data input, though introducing another rock/soil layer, with different mechanical characteristics represents no difficulty.

It is convenient to lead calculation to a 2D, plane strain or axisymmetrical, if it is possible, since it is less costly in data input and running time than a 3D analysis.

8.1. Formulation of the problem

Let us consider the following boundary problem: the rock mass is an infinite body in which a circular opening is made, assuming then that the underground opening is at a certain depth characterized by a hydrostatic primary (initial) stress, σPPIwhere P = γh, h is the depth at which the tunnel is dug, γ is the specific gravity of the rock and Ιis the unity tensor.

As the tunnel possesses a circular geometry, the rock-mass and the lining mechanical properties do not depend on the angular coordinate θ and the far stress field in situ is hydrostatic (primary stress componentsσv,σh are assumed equal), the problem is an axisymmetrical one in Orz plane (Figure 15).

Consequently, the primary stress componentsσv,σh are assumed equal. The boundary conditions are:

OnBCi.e.z[zB,zC]andr=a:σrr=0and  σrz=0.E32

Figure 15.

Domain and boundary conditions for the problem study in Orz plane along the tunnel axis.

Cristescu's elasto-viscoplastic constitutive law is used for the rock-mass and elastic behavior for the concrete lining.

An important factor of the analysis is the time effect and it is involved by two different aspects: the rheological behavior of the rock mass on and the excavation history. Moreover, the tunnel support mounting determines a problem of interaction, getting thus a more involved calculus. Another important factor of the analysis of ground-support interaction during the tunnel excavation is the face tunnel. For instance, since the behavior of the rock-mass is viscplastic, rock pressure on the lining increases in time. On the other hand, closer the lining is installed to the tunnel face, more the pressure at the rock-lining interface increases with the advancing of the tunnel face.

The state of stress and strain around a lined tunnel depends explicitly on:

The mechanical and geometrical characteristic of the rock-mass and the support;

The excavation conditions, such as excavation rate, generally the excavation phases;

The support mounting conditions, namely the support mounting time after the excavation and the distance between the lining and the tunnel face.

Concerning the geometry and the loading, the successive phases of the tunnel excavation and support mounting is a three-dimensional problem. However, there are certain cases when the problem may be simplified assuming the hypothesis that close to the tunnel face, on the tunnel walls, r=a, the decompression of the primary stress component is occurring gradually.

The calculation of tunnel excavation and lining mounting is a complex problem. On one hand, the excavation is a three-dimensional problem that imposes taking into consideration the tunnel face influence that means a gradual decompression of the primary stress σh of the rock mass on the opening walls. On the other hand, the support mounting determines the problem to be a massive-lining interaction one, mainly based on the behavior of the rock-lining interaction.

The lining is often installed quickly enough after the excavation and at a relatively small distance from the tunnel face that induces a complex combination of the effects mentioned above.

Consequently, it is very important to take into consideration both time and tunnel face effect in the excavation and lining mounting problem. Essentially in this study is the calculation of lining pressure and the convergence of the tunnel walls.

8.2. The numerical solution of the successive phases tunnel excavation

The domain discretization of the boundary value problem is presented in Figure 16. As usual, in the tunnel surface region the mesh must be quite refine, while elsewhere a minimum possible number of elements is considered. It is used an 8 nodded-quadrilateral mesh, at least 2 layers of quadrilateral elements in the concrete lining.

We consider the tunnel radius a = 1.2 m and the lining thickness 0.2 m. The depth at which the tunnel is excavated is 273.5 m. One phase duration is 12 hours and respectively 1 day.

For the rock mass the Borod coal is used, whose material constants were presented previously. For the concrete the following material constants are used: Young modulus E = 20000 MPa, Poisson coefficientν = 0.3 and the volumetric weight γ = 0.02 MN/m2.

Figure 16.

The domain discretization of the boundary value problem.

It is used a numbering of elements group as they are activated/deactivated in the excavation and lining installing processes, as follows:

1st group is the rock-mass considered infinite

2nd group corresponds to the already mounted lining

3rd group is in the first phase the rock-mass that is going to be excavated and in the second phase is replaced by the concrete

4th group is in the first phase the rock-mass that is going to be excavated

5th group is in the first two phases the rock-mass that is going to be excavated in the third phase and eventually replaced by concrete in a possible fourth phase

6th group is in the first two phases the rock-mass that is going to be excavated in the third phase

Numerical model concerns three successive phases’ tunnel excavation and lining mounting. Let us detail the progression of phases of the example:

Phase 1: tunnel excavation and calculation of strain, stress, damage parameter, without lining mounting on new excavated zone (Excavation of element group 3 and 4). Structure elements corresponding to the support are inactive (null mechanical characteristics). Output storage for the following phase (phase 2) is performed as well.

Phase 2: lining installing at a certain given time T0, on the unexcavated zone in phase 1. Displacement and stress initialization starting from the previous phase output storage and output storage for the following phase (phase 3) are performed as well.

Phase 3: tunnel face advancing on a distance of unit radius, namely 1.2 m. Realization of a new excavation is simulated by inactivation of element group 5 and 6, considering null mechanical characteristic. Displacement and stress initialization with the 2nd phase state is performed as well.

In the following, we present some important results of the calculation concerning, for instance, the normal stress, the equivalent stress or the damage parameter df.

In Figures 17a, b, c isovalues zones for the damage parameter corresponding to the three phases are presented, respectively. It is observed that in tunnel face zone the damage is maximum (white area).

Figure 17.

a),b),c) Isovalues zones for damage parameter for the first, second, third phase, respectively

Isovalues zones for normal stress corresponding to the three phases are presented in Figures 18a, 18b, 18c. Great concentrations are observed in tunnel face zone (white area).

Figures 19a, b, c present isovalues zones for equivalent stress corresponding to the three phases, observing small tractions in the second, respectively the third phase. That signifies the possibility of fracture by exceeding the traction resistance, since it is known that it is very low for the rocks.

Figure 18.

a),b),c) Isovalues zones for normal stress for the first, second, respectively third phase.

Figure 19.

a),b),c). Isovalues zones for equivalent stress the first, second, respectively third phase.

8.3.Parametric. analysis: Influence of tunnel depth, lining stiffness and lining mounting time

8.3.1. Support rigidity influence

The numerical solution is used to analyze the influence of different parameters that intervene in the excavation process. In this chapter we present the influence of lining rigidity and depth at which the tunnel is excavated.

The conclusion is that the displacement and the damage through dilatancy, as it was incorporated in the constitutive law, are decreasing with the increasing of the lining rigidity and increasing with depth increasing. In this paragraph we exhibit these features of the solution by several Figures and observations.

The support rigidity can be calculated by the following formula [1], [20]:


The previous calculation was performed for a value of Young modulus E=20000 Mpa, Poisson coefficient ν = 0.3, external radius b = 1.2 m and internal radius c = 1 m. We perform the calculation for E=2500 Mpa, too. Figures 20a, b present the isovalues zones of damage parameter d defined in relation (2) and show that the dilatancy, as it was incorporated in the constitutive law, is decreasing with the increasing of the lining rigidity. The same results are obtained by author by analytical means in [20], [18].

Figure 20.

a),b). Isovalues zones for damage parameter for E=2500 Mpa, E=20000 Mpa, respectively

Figure 21.

a),b),c). Isovalues zones of dilatancy for tunnel depth of 150 m, 273.5 m, 850 m, respectively

Figure 22.

a), b). Isovalues zones for the normal stress for tunnel depth of 150 m, 850 m, respectively

8.3.2. Tunnel depth influence

To analyze the depth influence on the processes we perform the calculation for different values of the depth at which the tunnel is excavated. At the previous calculation performed initially for a depth of 273.5 m, we add another two calculations for 150 m and 850 m depth respectively. The conclusion is that both displacement and damage increase with depth such that: the maximum radial displacement is 6.65 cm at a depth of 150 m, 6.83 cm at a depth of 273.5 m and 7.31 cm at a depth of 850 m. Concerning the damage, for instance, let us observe the Figures 21a, b, c presenting the isovalues zones of the damage parameter for a depth of 150 m, 273.5 and 850 m, respectively. The damage is more pronounced at a depth of 850 m (the white zone is more extended).

Figures 22a, b present isovalues zones of the normal stressσrrfor instance, for the same depths of 150 m and 850 m, respectively. It is observed that the normal stress decreases with the depth increase in compression and consequently, it appears the risk of traction at smaller depth, which may induces rock or lining fracture.


9. Conclusions

The complex problem of a lined tunnel excavation in a viscoplastic rock mass is approached in this chapter both numerically by a FE code proposed by the author and by an analytical approach too. A good agreement between the numerical and the analytical solutions is obtained.

Both solutions are studied then for further specific features. All these features are in good agrrement with the practical observation, so, more involved boundary problems could be developed with this code and further improvements for the analytical solution.

A special study is devoted to the finite element solution for the simulation of a tunnel excavation with successive tunnel face advancing and the lining mounting. Due to the symmetry of the geometry and loadings, the problem is treated as an axisymmetrical one with an additional emphasis of the three-dimensional aspect of the problem, namely the tunnel face advancing and its proximity influence. So, the approach of a tunnel calculation in two-dimensional analysis along the tunnel axes, simulating thus the three-dimensional aspect of the problem, is more realistic than the classical cross section analysis and obviously less costly than an actual three-dimensional analysis. The parametric analysis performed in this study by the numerical solution is in good accordance with the results obtained by the author by analytical means [19], [20].


  1. 1. CristescuN.1989Rock rheologyKluwer Academic Publishers, Dordrecht, Holland.
  2. 2. SimionescuO.1999Mathematical methods in underground structure termomechanics, Romanian Academy Publ., Bucharest.
  3. 3. BernaudD.1991Tunnels profonds dans les milieux viscoplastique: approches expérimentale et numérique. Thése. École National des Ponts et Chaussées, France.
  4. 4. PanetM.1974Stabilité et soutenement des tunnels. La méc.des roches appl. aux ouvrages de génie civil,Ch.IX,ENPC.
  5. 5. SulemJ.PanetM.GuenotA.1987An analytical solution for time-dependent displacements in a circular tunnel, Int.J. Rock Mech. Min. Sci.& Geomech. Abstr., 243155
  6. 6. Zienkiewicz O C, Cormeau I C1974Visco-plasticity- plasticity and creep in elastic solids- A unified numerical solution approach, Int.J.Num.Meth.Engng.,8,821845
  7. 7. Cormeau, I.C., Numerical stability in quasi-static elasto/visco-plasticity1975Int. J. Num. Meth.Engng., 9109
  8. 8. IonescuI. R.SofoneaM.1993Functional and numerical methods in viscoplasticityOxford Univ. Press Oxford.
  9. 9. ChanA. H.OuJ.2008Three-Dimensional Numerical Analysis of a Dynamic Structure, Saturated Soil and Pore Fluid Interaction Problem, Trends in Eng. Comput. Techn., 335353
  10. 10. ZienkiewiczO. C.ChanA. H. C.PastorM.SchreflerB. A.ShiomiT.1999Computational Geomechanics with Special Refernces to Earthquake Engineering, J. Wiley & Sons, Chicester, 1999.
  11. 11. GiodaG.1993Finite element analysis of time dependent effects in tunnels, Advanced school on visco-plastic behaviour of geomaterials, CISM, Udine, Oct.
  12. 12. LauroF.BennaniB.DrazeticP.OudinJ.NiX.1997Ductile damage and fracture finite element modelling of elasto-viscoplastic voided materialsComputational Materials Science73295307
  13. 13. LinR. C.BrocksW.2006On a finite-strain viscoplastic law coupled with anisotropic damage: theoretical formulations and numerical applicationsArchive of Appl.Mech.,75315 EOF
  14. 14. Augarde C E, Burd H J2001D FE analysis of lined tunnels, Int. J. Num. Analyt. Meth. in Geomechanics, 25243262
  15. 15. KomiyaK.SogaK.AkagiH.HagiwaraT.BoltonM. D.1999FE, modeling of excavation and advancement process of a shield tunnelling machine, Soil &Foundations, 393
  16. 16. Leem J, Kemeny J M (1999) Finite element micromechanical-based model for hydro-mechanical coupling, The 37th U.S. Symp. on Rock Mech. (USRMS), June 7-9, Balkema, Vail, CO. Wang J, Wang W (2009) Multi-field Coupled Model&Num.Sim. During Excav.TunnelModeling, WMSO08,pp.334-337.
  17. 17. RoateşiS.AriciucM.2005Study on depth dependency in tunnel calculation, Scientifical Bulletin, UTCB, 1, 1630
  18. 18. RoateşiS.1996The tunnel face influence in the analysis of a circular tunnel with a time-dependent behaviour, Rev.Roum.Sci.Tech. Mec.Appl., 3-4247263
  19. 19. RoateşiS.2005Mathematical modeling and FEM problems in tunnel calculation,Romanian Academy Publisher, Bucharest.
  20. 20. CESAR − LCPC Documentation, Version 3.0.x, Publication LCPC, 1992.
  21. 21. Lade P V, Nelson R B, Ito Y M1987Nonassociated flow and stability of granular materialsJ.Engng.Mech., 113, 13021318

Written By

Simona Roatesi

Submitted: 16 November 2011 Published: 10 October 2012