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

Simona Roatesi

300 Finite Element Analysis – New Trends and Developments

*Systems and Signal Processing 25,* vol. 1393–1407, 2011.

[21] R. Gutkin*, et al.*, "On acoustic emission for failure investigation in CFRP," *Mechanical* 

Additional information is available at the end of the chapter

http://dx.doi.org/10.5772/46168

## **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 (convergenceconfinement 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

© 2012 Roatesi, licensee InTech. This is an open access chapter distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/3.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited. © 2012 Roatesi, licensee InTech. This is a paper distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/3.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

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

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

**σ**

iii. The component of the elastic strain rate satisfies the Hooke law

**ε σ**

into volumetric and deviatoric parts and it is given by the quantity:

0

We introduce the *damage parameter d*, see [1], defined by:

characterized by the following parameter (constant):

two stress invariants mentioned above.

constitutive equation.

*T*

iv. The component of the irreversible strain rate satisfies:

respective function:

 11 1 32 2 **I** *<sup>E</sup>*

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

( ) ( ,)1 (,)

*<sup>I</sup> Wt F k d*

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

*A AA A* / 2

The irreversible stress work is used as a hardening parameter or internal state variable, split

*I I*

describing the energy released by micro-cracking during the entire dilatancy period. In (2) *max <sup>t</sup>* represents the time for which *<sup>I</sup> Wv* is maximum. The failure threshold is considered to be the total energy released by micro cracking during the entire dilatancy process and it is

. *I I*

*H*( **σ** ) represents the loading function, generally a function of stress tensor **σ** , with , *<sup>I</sup> H Wt* the creep stabilization boundary equation, the function *H* depending on the

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

*F*( **σ** ) represents a viscoplastic potential, that establishes the orientation of .

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

We consider therefore, the following constitutive equation:

**σ ε** .

*I I II W T t t dt W T W T v d* (1)

*v max v dt W t W t* (2)

*<sup>f</sup> v max v failure d Wt Wt* (3)

*KG G*

*I*

*<sup>H</sup>* ,

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


iii. The component of the elastic strain rate satisfies the Hooke law

$$\dot{\mathbf{o}}^E = \left(\frac{1}{3K} - \frac{1}{2G}\right)\dot{\mathbf{o}}\mathbf{I} + \frac{1}{2G}\dot{\mathbf{o}}\mathbf{I}$$

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

iv. The component of the irreversible strain rate satisfies:

302 Finite Element Analysis – New Trends and Developments

tunnel with a time-dependent behaviour, etc.

**2. The elasto-viscoplastic constitutive equation** 

an important physical meaning are:


 2 2 3 3

hypotheses for the constitutive equation formulation are (see [1]):

123 <sup>3</sup> ;


constitutive module.

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 hydromechanical 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

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

The constitutive law implemented in finite element code is proposed by Cristescu. The

i. 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 3123 , , , , , . Among the stress invariants, those with

ii. The displacements and the rotations are assumed small, so that *E I* , where *<sup>E</sup>* and being the elastic strain rate and the irreversible strain rate respectively.

1 2 3 12 23 31 (or the octahedral shear stress

*II* with *II* being the second invariant of the stress deviator).

$$\dot{\mathbf{e}}^I = k(\mathbf{o}\sigma, d) \left\langle 1 - \frac{W^I(t)}{H(\overline{\mathbf{o}}, \sigma)} \right\rangle \frac{\partial F}{\partial \mathbf{o}} \,\,\,\, \dot{\mathbf{e}}$$

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:

$$
\left\langle A \right\rangle = \left( A + \left| A \right| \right) / \mathcal{Q} = A^{+} 
$$

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:

$$\mathcal{W}^I\left(T\right) = \bigcap\_{0}^{T} \mathbf{o}\left(t\right) \cdot \dot{\mathbf{e}}^I\left(t\right) dt = \mathcal{W}\_v^I\left(T\right) + \mathcal{W}\_d^I\left(T\right). \tag{1}$$

We introduce the *damage parameter d*, see [1], defined by:

$$d\begin{pmatrix} t \\ \end{pmatrix} = \mathcal{W}\_v^I \begin{pmatrix} t\_{\max} \\ \end{pmatrix} - \mathcal{W}\_v^I \begin{pmatrix} t \\ \end{pmatrix} \tag{2}$$

describing the energy released by micro-cracking during the entire dilatancy period. In (2) *max <sup>t</sup>* represents the time for which *<sup>I</sup> Wv* is 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):

$$d\_f = \mathcal{W}\_v^I \left( t\_{\max} \right) - \mathcal{W}\_v^I \left( t\_{failure} \right). \tag{3}$$

*H*( **σ** ) represents the loading function, generally a function of stress tensor **σ** , with , *<sup>I</sup> H Wt* 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 .


We consider therefore, the following constitutive equation:

$$\dot{\mathbf{e}} = \left(\frac{1}{3K} - \frac{1}{2G}\right) \mathbf{e} \mathbf{I} + \frac{1}{2G} \mathbf{e} + k(\mathbf{o}, d) \left\langle 1 - \frac{W^I(t)}{H(\overline{\mathbf{e}}, \mathbf{o})} \right\rangle \frac{\partial F}{\partial \mathbf{o}} \tag{4}$$

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

$$H(\boldsymbol{\sigma}, \overline{\boldsymbol{\sigma}}) \coloneqq a\_0 \frac{\left(\overline{\boldsymbol{\sigma}} / \boldsymbol{\sigma}\_\*\right)^2}{a\_2 \boldsymbol{\sigma} / \boldsymbol{\sigma}\_\* + a\_1} + b\_0 \overline{\boldsymbol{\sigma}} / \boldsymbol{\sigma}\_\* + \begin{cases} c\_0 \sin(\boldsymbol{\alpha} \boldsymbol{\sigma} / \boldsymbol{\sigma}\_\* + \boldsymbol{\phi}) + c\_1 & \text{if} \quad 0 \le \boldsymbol{\sigma} \le \boldsymbol{\sigma}\_0 \\ c\_0 + c\_1 & \text{if} \quad \boldsymbol{\sigma}\_0 \le \boldsymbol{\sigma} \end{cases} \tag{the yield surface}$$

$$F(\sigma, \overline{\sigma}) \coloneqq H(\sigma, \overline{\sigma}) \tag{5}$$

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

**3. The numerical integration of the elasto-viscoplastic equation** 

Runge - Kutta method of fourth order.

, ,1 , **σ σ**

the irreversible work increment respectively.

*W<sup>i</sup>* .

*W t <sup>I</sup> <sup>F</sup> fWkd <sup>H</sup>* .

<sup>1</sup>

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

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

*vp f W* ( ) **<sup>σ</sup>** , 2 1 <sup>1</sup>

corresponding to the viscoplastic strain increment <sup>1</sup>

**σ**

information to match an expansion in Taylor series of higher orders.

 

  1

*vp vp*

*vp*

*vp*

a. Semi implicit Euler method ( scheme) in which the evaluation of the viscoplastic

strain increment is done with the rule: 1 2 ( 1 ) ) , [0,1] *vp vp vp <sup>t</sup>* where

associated irreversible work increment calculated with the standard scheme,

stress and irreversible work at the beginning of the step and function

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

In particular, Runge-Kutta method of fourth order employed to integrate the proposed

,

*f WW*

,

*f WW*

,

*f WW*

constitutive equation uses four evaluations of function *f* for the considered time step:

,

**σ**

*f W*

**σ σ**

**σ σ**

**σ σ**

2 11

3 22

4 33

with 123 1 2 3 **σσσ** *WWW* , , the intermediate evaluations of the stress increment, and

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

*vp f WW W* ( , ) and , **σ σ σ** are the stress increment and

*vp* ; **σ**, *W* are the known values of

where *a*0 = 7.65 10−4 MPa; *a*1 = 0.55; *a*2 = 8.159 10−3; *b*0 = 0.001 MPa; *c*0 = 4.957 10−4 MPa; *c*1 = 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(\sigma,\overline{\sigma}) \coloneqq a \frac{\left(\overline{\sigma}\right)^{7}}{\left(\sigma - \frac{\overline{\sigma}^{3}}{3}\right)^{5}} + b\overline{\sigma} + c\sigma \tag{\text{the yield surface}}$$

$$F(\sigma,\overline{\sigma}) \coloneqq -\frac{2h\_{1}\sigma^{\frac{1}{2}}}{2f+\alpha} + 2fh\_{1}\left(\frac{2}{3}\frac{\sigma^{\frac{1}{2}}}{2f+\alpha} + \frac{2\left(1 + \frac{\alpha}{3}\right)^{\frac{1}{2}}}{\left(2f+\alpha\right)^{3}}\right) +$$

$$\begin{split} \left(h\_{1}\overline{\sigma}\frac{\left[\left(1 + \frac{\alpha}{3}\right)\overline{\sigma}\right]^{\frac{1}{2}}}{\left(2f+\alpha\right)^{\frac{1}{2}}} - 2fh\_{1}\frac{\left[\left(1 + \frac{\alpha}{3}\right)\overline{\sigma}\right]^{\frac{1}{2}}}{\left(2f+\alpha\right)^{\frac{2}{3}}}\right) \ln\frac{\left[\left(2f+\alpha\right)\sigma\right]^{\frac{1}{2}} + \left[\left(1 + \frac{\alpha}{3}\right)\overline{\sigma}\right]^{\frac{1}{2}}}{\left[\left(2f+\alpha\right)^{\frac{1}{2}} - \left[\left(1 + \frac{\alpha}{3}\right)\overline{\sigma}\right]^{\frac{1}{2}}} + \\\\ \overline{g\_{0}\overline{\sigma}} + g\_{1}\overline{\sigma}^{\frac{3}{2}} \qquad\qquad\qquad\qquad\qquad\qquad\qquad\left(\text{the plastic potential}\right) \end{split}$$

where -7 1 *<sup>a</sup>*=4.834 10 (kPa) , <sup>3</sup> *<sup>b</sup>*=1.33 10 , <sup>3</sup> *<sup>c</sup>*=1.058 10 , <sup>1</sup> 2 <sup>1</sup> *h* =2.34 10 (kPa) , <sup>0</sup> *g* =0.005 , -6 -2 <sup>1</sup> *<sup>g</sup>* =0.62 10 (kPa) , *f*=0.562, =1.34, -6 1 *<sup>k</sup>*=3.6 10 s (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 *nonassociated 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:

a. Semi implicit Euler method ( scheme) in which the evaluation of the viscoplastic strain increment is done with the rule: 1 2 ( 1 ) ) , [0,1] *vp vp vp <sup>t</sup>* where

 <sup>1</sup> *vp f W* ( ) **<sup>σ</sup>** , 2 1 <sup>1</sup> *vp f WW W* ( , ) and , **σ σ σ** are the stress increment and associated irreversible work increment calculated with the standard scheme, corresponding to the viscoplastic strain increment <sup>1</sup> *vp* ; **σ**, *W* are the known values of stress and irreversible work at the beginning of the step and function

$$f\left(\sigma, W\right) = k\left(\sigma, d\right) \left\langle 1 - \frac{W^I\left(t\right)}{H\left(\sigma, \overline{\sigma}\right)} \right\rangle \frac{\partial F}{\partial \sigma} \dots$$

304 Finite Element Analysis – New Trends and Developments

functions and material constants used are:

*Ha b*

*c*1 = 4.8955 10−4 MPa; =171.9270;

determines an *associated flow rule*.

material constants are:

2

0 0 \*

*a a c c*

*k d*

 11 1 ( ) ( ,)1 32 2 (,) **σI σ σ**

For instance, for the model describing the *Borod coal* behavior (see [1]), the constitutive

\* 0 \*1 0

where *a*0 = 7.65 10−4 MPa; *a*1 = 0.55; *a*2 = 8.159 10−3; *b*0 = 0.001 MPa; *c*0 = 4.957 10−4 MPa;

The relation (5) shows that the viscoplastic potential equals the loading function, fact that

For the model describing the *saturated sand* behaviour, the constitutive functions and

 

> 3 3 2 2

 

<sup>2</sup> <sup>2</sup> <sup>3</sup> ( , ): <sup>2</sup> 2 3 2 2

*<sup>h</sup> <sup>F</sup> fh f f <sup>f</sup>*

2 2 1

<sup>1</sup> *<sup>g</sup>* =0.62 10 (kPa) , *f*=0.562, =1.34, -6 1 *<sup>k</sup>*=3.6 10 s (the viscosity coefficient) and the elastic

For this model, as the viscoplastic potential differs to the loading function, determines a *non-*

 

1 1 21 3 3 3

2 ln 2 2 2 1

( ) ( , ):

7 <sup>5</sup> <sup>3</sup>

3

*c c*

**σ**

(the yield surface)

*Wt F <sup>I</sup>*

*F H* ( , ): ( , ) (5)

0=1,0996MPa; *k* =6 10−6; s-1 (the viscosity coefficient); *df* =

*H a bc* (the yield surface)

2 1

2

2

3

1 2

 

2 <sup>1</sup> *h* =2.34 10 (kPa) , <sup>0</sup> *g* =0.005 ,

3 2

*f*

1 3 1 2 2 1 2

1 3

= 0,327.

*KG G H* (4)

2 \*1 01 0 (/ ) sin( / ) if 0 ( , ): / / if

0 1 *g g* the plastic potential

= 0.33.

1 1

*h fh*

3


*associated flow rule*.

constants are: *E* = 205300 kPa,

3 5

*f f <sup>f</sup>*

where -7 1 *<sup>a</sup>*=4.834 10 (kPa) , <sup>3</sup> *<sup>b</sup>*=1.33 10 , <sup>3</sup> *<sup>c</sup>*=1.058 10 , <sup>1</sup>

1

4.48 103 Jm−3; \* =1 MPa and the elastic constants are: *E* = 798.385 MPa,

b. 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:

$$\begin{aligned} \dot{\mathbf{s}}\_{vp}^1 &= f\left(\boldsymbol{\sigma}\_\prime W\right) \\ \dot{\mathbf{s}}\_{vp}^2 &= f\left(\boldsymbol{\sigma} + \boldsymbol{\Delta}\boldsymbol{\sigma}^1\prime, W + \boldsymbol{\Delta}W^1\right) \\ \dot{\mathbf{s}}\_{vp}^3 &= f\left(\boldsymbol{\sigma} + \boldsymbol{\Delta}\boldsymbol{\sigma}^2\prime, W + \boldsymbol{\Delta}W^2\right) \\ \dot{\mathbf{s}}\_{vp}^4 &= f\left(\boldsymbol{\sigma} + \boldsymbol{\Delta}\boldsymbol{\sigma}^3\prime, W + \boldsymbol{\Delta}W^3\right) \end{aligned}$$

with 123 1 2 3 **σσσ** *WWW* , , the 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 *W<sup>i</sup>* .

## **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:

$$\Psi\left(\mathbf{u},\boldsymbol{\lambda}\right) = \mathbf{R}(\mathbf{u}) - \boldsymbol{\lambda}(t)\mathbf{P} \tag{6}$$

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

*ii i FP v* that has to be calculated as accurate as possible.

The calculation of the global stiffness matrix **K B DB**<sup>d</sup> *<sup>t</sup> <sup>v</sup>* , where the matrix **D** needs not be

*i ii* <sup>1</sup> calculation of the incremental stresses with the constitutive law; can be

*i ii* , with / *vp*

\* **D***evp*

From the law ( ) *i i f* (implicit form), the numerical integration may be done with θ

where **u** represents the vector of nodal displacements, **v** is the vector of out of balance nodal

the known nodal applied forces, ( ) **R** *<sup>i</sup> u* is the vector of the equivalent nodal forces, due to the stresses *<sup>i</sup>* . 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



*<sup>i</sup> f* ,proportional with the

represents the vector of

*i i* (explicit form).

Setting up the term **<sup>B</sup>** <sup>d</sup> *<sup>j</sup> <sup>t</sup>*


<sup>1</sup> **<sup>K</sup>** <sup>F</sup> *i i <sup>u</sup>* solving the system of equations

*i i* . - viscoplasticity: \* **D** ( ) *<sup>e</sup> vp*

derivative of the viscoplastic potential.

*i i u v* calculation of the residue

that two calculation levels can be distinguished:

**D,**( ) **D** and the stresses *<sup>i</sup>* <sup>1</sup> .

*i i u* equilibrium check

scheme, with Runge-Kutta method, the consistent matrix, etc.

forces, **B** designates the matrix of the shape functions for strain, **P***<sup>j</sup>*

 *i ii* <sup>1</sup> *u uu* updating displacement *<sup>i</sup>* **B** <sup>1</sup> *u* displacement-strain relations


The solving of the algebraic system **K** F *i i u* (nonlinear if **K=K(u)**)

quite accurate and can be:

calculated as:

1 1 **R B** () d *<sup>t</sup>*

1 1 P () **<sup>R</sup>** *<sup>j</sup>*

Convergence test:

no *i = i+*1 (next iteration) yes *j = j+*1 (next step)

displacements.


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 *u*1. 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:

$$\begin{cases} \boldsymbol{\mu\_1} = 0\\ \boldsymbol{\mu\_1} = \mathbf{P} \end{cases} \text{initialization}$$

*j* = 1 loop on the increments

The incrementation of the loading **P** 

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

Setting up the term **<sup>B</sup>** <sup>d</sup> *<sup>j</sup> <sup>t</sup> ii i FP v* that has to be calculated as accurate as possible.

The calculation of the global stiffness matrix **K B DB**<sup>d</sup> *<sup>t</sup> <sup>v</sup>* , where the matrix **D** needs not be quite accurate and can be:


The solving of the algebraic system **K** F *i i u* (nonlinear if **K=K(u)**)

<sup>1</sup> **<sup>K</sup>** <sup>F</sup> *i i <sup>u</sup>* solving the system of equations

*i ii* <sup>1</sup> *u uu* updating displacement

*<sup>i</sup>* **B** <sup>1</sup> *u* displacement-strain relations

 *i ii* <sup>1</sup> calculation of the incremental stresses with the constitutive law; can be calculated as:


306 Finite Element Analysis – New Trends and Developments

**4. The numerical solution of nonlinear problems** 

applied loading and the searching of the limit loading.

system of equations that may be put in the form:

the framework of the finite element program is:

applied at moment t.

equations to be solved.

loading of the structure.

*j* = 1 loop on the increments The incrementation of the loading **P** 

calculated).

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

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

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

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

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

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 *u*1. For this reason it is necessary to perform an incremental

The algorithm for solving a static nonlinear problem with the initial stress method used in

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

<sup>0</sup> initialization

 1 1

*u*

**P**

**u Ru P** , ( ) ()*t* (6)

with the stress and strain history at that point, a new state of stress to be calculated.


$$\mathbf{^\*}\ \Delta \sigma\_i = \mathbf{D}^{evp} \Delta \varepsilon\_i \tag{\text{explicit form}}.$$

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

 1 1 **R B** () d *<sup>t</sup> i i u v* calculation of the residue

 1 1 P () **<sup>R</sup>** *<sup>j</sup> i i u* 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, **P***<sup>j</sup>* represents the vector of the known nodal applied forces, ( ) **R** *<sup>i</sup> u* is the vector of the equivalent nodal forces, due to the stresses *<sup>i</sup>* . 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:


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:

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

*KG G* (9)

**σ σ**

> **σ σ**

(8)

(10)

0

**σ**

Introducing the left hand side of the relation (7) into the constitutive equation (4), we get the

*F H*

*<sup>H</sup> k F t t*

0 0 0 11 1 ; 0; 32 2 **ε ε I σ** *<sup>I</sup> t t*

where **σ σ** <sup>0</sup> *t t* is the initial stress state, reached instantaneously and taken with respect to the state at 0*t* 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

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

*<sup>t</sup> <sup>H</sup> k F t t <sup>t</sup> G H F*

1 exp , <sup>2</sup> <sup>1</sup>

1 exp , <sup>2</sup> <sup>1</sup>

where **σ**<sup>1</sup> and **σ**<sup>2</sup> being the principal stresses and where during the time interval <sup>0</sup> *t tt*, the state of stress is constant and <sup>0</sup> *W t <sup>I</sup>* is calculated at moment 0*<sup>t</sup>* , just before of a new stress

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

*<sup>t</sup> <sup>H</sup> k F t t <sup>t</sup> G H F*

 

*I*

1 1

<sup>1</sup> ,

2 2

*H*

**5.2. Comparison of the numerical results with the analytical solution** 

<sup>1</sup> ,

*H*

0

*W t F*

1 0

**σ σ**

 

*I*

0

*W t F*

2 0

**σ σ**

 

<sup>1</sup> ,

*H*

0

*W t <sup>I</sup> F*

**σ σ**

<sup>0</sup>

1 exp <sup>1</sup> **σ ε ε σ**

total strain under constant stress:

using the initial conditions:

expression

namely:

 1 . , *W t <sup>I</sup> H*

increment **σ** that occurs at <sup>0</sup>*t* .

0


$$\text{{-}s}\quad\text{{+}s}\quad\text{{+}s}\quad\text{{+}s}\quad\text{{-}\quad\text{{+}}\quad\text{{-}dC}\quad\text{{+}\${}^{\mu}\${}^{\nu}\${}^{\mu}\${}^{\nu}\${}^{\mu}\${}^{\nu}\${}^{\mu}\${}^{\nu}\${}^{\nu}\${}^{\mu}\${}^{\nu}\${}^{\nu}\${}^{\mu}\${}^{\nu}\${}^{\nu}\${}^{\mu}\${}^{\nu}\${}^{\nu}\${}^{\mu}\${}^{\nu}\${}^{\nu}\${}^{\mu}\${}^{\nu}\${}^{\nu}\${}^{\mu}\${}^{\nu}\${}^{\nu}\${}^{\mu}\${}^{\nu}\${}^{\nu}\${}^{\mu}\${}^{\nu}\${}^{\nu}\${}^{\mu}\${}^{\nu}\${}^{\nu}\${}^{\mu}\${}^{\nu}\${}^{\nu}\${}^{\mu}\${}^{\nu}\${}^{\nu}\${}^{\mu}\${}^{\nu}\${}^{\nu}\${}^{\mu}\${}^{\nu}\${}^{\nu}\${}^{\mu}\${}^{\nu}\${}^{\nu}\${}^{\mu}\${}^{\nu}\${}^{\nu}\${}^{\mu}\${}^{\nu}\${}^{\nu}\${}^{\nu}\${}^{\mu}\${}^{\nu}\${}^{\nu}\${}^{\nu}\${}^{\mu}\${}^{\nu}\${}^{\nu}\${}^{\mu}\${}^{\nu}\${}^{\nu}\${}^{\nu}\${}^{\mu}\${}^{\nu}\${}^{\nu}\${}^{\nu}\${}^{\mu}\${}^{\nu}\${}^{\nu}\$$$

$$\text{--} \qquad \text{testing the work} \quad \frac{\left| \Psi(\boldsymbol{\mu}\_{i})^{T} . \boldsymbol{\Delta\boldsymbol{\mu}\_{i}} \right|}{\left| \Psi(\boldsymbol{\mu}\_{0})^{T} . \boldsymbol{\Delta\boldsymbol{\mu}\_{0}} \right|} \leq \boldsymbol{\varepsilon}\_{W} . \text{ }$$

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 elastoviscoplastic 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 moment <sup>0</sup>*t* the stress state is increased suddenly to the value of **σ** <sup>0</sup>*t* 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 *W t <sup>I</sup>* when all stress components are constant; it is easily obtained by integration of the constitutive law (4) (see [1]):

$$1 - \frac{W^I\left(t\right)}{H\left(\sigma, \overline{\sigma}\right)} = \left(1 - \frac{W^I\left(t\_0\right)}{H\left(\sigma, \overline{\sigma}\right)}\right) \exp\left[\frac{k}{H} \frac{\partial F}{\partial \sigma} \cdot \sigma\left(t\_0 - t\right)\right] \tag{7}$$

where **σ σ** <sup>0</sup> *t t* constant and <sup>0</sup> *W t <sup>I</sup>* is the initial value of *W<sup>I</sup>* for <sup>0</sup> *t t* , denoted further *WIP* .

This is the relation that describes the variation in time of *W t <sup>I</sup>* under a constant stress. The variation in time of *W t <sup>I</sup>* 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:

$$\mathbf{e} = \mathbf{e}^{0} + \frac{\left(1 - \frac{W^{I}\left(t\_{0}\right)}{H\left(\sigma, \overline{\boldsymbol{\sigma}}\right)}\right) \frac{\partial F}{\partial \boldsymbol{\sigma}}}{\frac{1}{H} \frac{\partial F}{\partial \boldsymbol{\sigma}} \cdot \boldsymbol{\sigma}} \left\{1 - \exp\left[\frac{k}{H} \frac{\partial F}{\partial \boldsymbol{\sigma}} \cdot \boldsymbol{\sigma}\left(t\_{0} - t\right)\right]\right\} \tag{8}$$

using the initial conditions:

308 Finite Element Analysis – New Trends and Developments

 <sup>0</sup> ( ) ; ( ) *i*

( ).

*T i i T W*

*u u u u*

0 0

( ).

*u u*

*i u u*

increased suddenly to the value of **σ** <sup>0</sup>*t* and it is kept constant.

**5.1. Determination of the analytical solution for creep step** 

easily obtained by integration of the constitutive law (4) (see [1]):

  *R*

 ; *<sup>i</sup> D*

.




further *WIP* .

viscosity parameter *k*.

quantities:

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

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 elastoviscoplastic 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 moment <sup>0</sup>*t* the stress state is

To establish the formulae that describe the creep deformation, we will write firstly the formula which supplies the variation of *W t <sup>I</sup>* when all stress components are constant; it is

> 0 <sup>0</sup> 1 1exp , , **<sup>σ</sup>**

*Wt Wt I I k F t t*

where **σ σ** <sup>0</sup> *t t* constant and <sup>0</sup> *W t <sup>I</sup>* is the initial value of *W<sup>I</sup>* for <sup>0</sup> *t t* , denoted

This is the relation that describes the variation in time of *W t <sup>I</sup>* under a constant stress. The variation in time of *W t <sup>I</sup>* is longer or shorter depending especially on the value of the

**σ**

*H H <sup>H</sup>* (7)

$$t = t\_0; \quad \mathbf{e}^I = 0; \quad \mathbf{e} = \mathbf{e}^0 = \left(\frac{1}{3K} - \frac{1}{2G}\right)\mathbf{e}\_0\mathbf{I} + \frac{1}{2G}\mathbf{e}\_0\tag{9}$$

where **σ σ** <sup>0</sup> *t t* is the initial stress state, reached instantaneously and taken with respect to the state at 0*t* 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 1 . , *W t <sup>I</sup> 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:

$$\begin{split} e\_1(t) &= \frac{\sigma\_1(t)}{2G} + \frac{\left(1 - \frac{W^I(t\_0)}{H(\mathfrak{o}, \overline{\mathfrak{o}})}\right) \frac{\partial F}{\partial \sigma\_1}}{\frac{1}{H} \frac{\partial F}{\partial \mathfrak{o}} \cdot \mathfrak{o}} \left\{1 - \exp\left[\frac{k}{H} \frac{\partial F}{\partial \mathfrak{o}} \cdot \mathfrak{o}\left(t\_0 - t\right)\right]\right\}, \\ e\_2(t) &= \frac{\sigma\_2(t)}{2G} + \frac{\left(1 - \frac{W^I(t\_0)}{H(\mathfrak{o}, \overline{\mathfrak{o}})}\right) \frac{\partial F}{\partial \sigma\_2}}{\frac{1}{H} \frac{\partial F}{\partial \mathfrak{o}} \cdot \mathfrak{o}} \left\{1 - \exp\left[\frac{k}{H} \frac{\partial F}{\partial \mathfrak{o}} \cdot \mathfrak{o}\left(t\_0 - t\right)\right]\right\}, \end{split} \tag{10}$$

where **σ**<sup>1</sup> and **σ**<sup>2</sup> being the principal stresses and where during the time interval <sup>0</sup> *t tt*, the state of stress is constant and <sup>0</sup> *W t <sup>I</sup>* is calculated at moment 0*<sup>t</sup>* , just before of a new stress increment **σ** that occurs at <sup>0</sup>*t* .

#### **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 <sup>10</sup> 10 s. 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.

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

*z l r ut t t*

*r rz*

*zz rz*

*z l r a t ft t t*

*r rz*

*rr rz*

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.

left lateral surface 0, and 0 0, 0, 0. right lateral surface 0, and ( ), 0, 0.

*z r a ut t t*

*z l r a t gt t t*

We consider for the functions *f* and g some particular forms that can simulate the triaxial test



for 0 (the hydrostatic stage) for (the deviatoric stage)

with a linear variation during the hydrostatic stage. , *<sup>h</sup>* ,*T* are constants that characterize

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

for 0 (the hydrostatic stage)

for 0 (the hydrostatic stage)

for (the deviatoric stage)

for (the deviatoric stage)

for 0 (the hydrostatic stage)

for (the deviatoric stage)

steps for this test) until a certain value is reached.

 

*<sup>t</sup> t T f t <sup>T</sup>*

*<sup>T</sup> g t t T T t*

*T t*

*<sup>t</sup> t T*

is increased (in steps as well) until failure.

The following functions are used for the two stages:

*h*

*h*

*h*

*h h*

2 2

2 2

*<sup>t</sup> t T f t <sup>T</sup>*

*h*

*h h*

*T t*

2

*t T T t*

*<sup>t</sup> t T*

*h*

*h*

*<sup>T</sup> g t*

 

inner face 0 and 0, 0, 0, 0.

outer face and 0, ( ), 0, 0.

The initial boundary problem is then:

in two stages, namely:

the loading.

through the functions:

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 <sup>000</sup> 1 23 0.7kPa, 0.1kPa, 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:

$$\varepsilon\_1^{\text{analytic}} = 3.08839 \times 10^{-6}, \varepsilon\_2^{\text{analytic}} = 7.99371 \times 10^{-7}, \varepsilon\_3^{\text{analytic}} = 7.99371 \times 10^{7} \text{ yr}$$

while the same components using θ scheme are:

 -6 -7 -7 <sup>123</sup> = 3.08841 10 , =-7.99383 10 , =-7.99383 10 ; *Euler Euler Euler*

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

 -6 -7 -7 <sup>123</sup> = 3.08840 10 , =-7.99378 10 , =-7.99378 10 . *R K R K R K*

## **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* 1 2 <sup>3</sup> increasing, 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:

310 Finite Element Analysis – New Trends and Developments

method of fourth order shows the superiority of this method.

while the same components using θ scheme are:

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

1 2 <sup>3</sup> increasing, under axi-symmetric hypothesis, see [20].

**(0,1)**

**(0,.5)**

**(0,0)**

**Figure 1.** The mesh for the simulation of the triaxial test for a cylindrical sample

**problem and the experimental data** 

increment is applied.

gives:

<sup>000</sup>

interval *t* , namely <sup>10</sup> 10 s. 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

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

1 23 0.7kPa, 0.1kPa, 0.1kPa. The numerical results obtained using Runge-Kutta

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

 -6 -7 -7 <sup>123</sup> = 3.08839 10 , =-7.99371 10 , =-7.99371 10 ; *analytic analytic analytic*

 -6 -7 -7 <sup>123</sup> = 3.08841 10 , =-7.99383 10 , =-7.99383 10 ; *Euler Euler Euler*

 -6 -7 -7 <sup>123</sup> = 3.08840 10 , =-7.99378 10 , =-7.99378 10 . *R K R K R K*

**6. Comparison of the numerical solution of the triaxial test boundary** 

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*

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.

**1 2 3**

**(.5,0) (1,0)**

**(1,.5)**

**(1,1)**

**4 5**

**678**

**(.5,1)**

$$\begin{array}{l}\text{(inner face)}\quad z=0 \text{ and } r \in [0,a] \quad u\_r\left(t\right)=0, \ \sigma\_{rz}\left(t\right)=0, \ t \ge 0.\\ \text{(left lateral surface)}\quad z \in \left[0,l\right] \text{ and } r=0 \quad u\_r\left(t\right)=0, \ \sigma\_{rz}\left(t\right)=0, \ t \ge 0.\\ \text{(right lateral surface)}\quad z \in \left[0,l\right] \text{ and } r=a \quad \sigma\_{rr}\left(t\right)=f\left(t\right), \ \sigma\_{rz}\left(t\right)=0, \ t \ge 0.\\ \text{(outer face)}\quad z=l \text{ and } r \in \left[0,a\right] \text{ } \sigma\_{zz}\left(t\right)=\mathbf{g}\left(t\right), \ \sigma\_{rz}\left(t\right)=0, \ t \ge 0.\end{array}$$

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


The following functions are used for the two stages:

$$f\left(t\right) = \begin{cases} \sigma\_h \frac{t}{T} & \text{for } 0 \le t \le T \\\\ \sigma\_h \text{ for } T \le t & \text{(the deviatoric stage)} \end{cases}$$

$$g\left(t\right) = \begin{cases} \quad \sigma\_h \frac{t}{T} \text{ for } 0 \le t \le T \\\\ \sigma\_h + \Sigma\_h \left(\frac{t-T}{\tau}\right) \text{for } T \le t & \text{(the deviatoric stage)} \end{cases}$$

with a linear variation during the hydrostatic stage. , *<sup>h</sup>* ,*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) = \begin{cases} \sigma\_h \frac{t^2}{T'^2} \text{ for } 0 \le t \le T' & \text{(the hydrostatic stage)}\\ \sigma\_h \text{ for } T' \le t & \text{(the deviatoric stage)} \end{cases}$$

$$\text{g'}(t) = \begin{cases} \sigma\_h \frac{t^2}{T'^2} \text{ for } 0 \le t \le T' & \text{(the hydrostatic stage)}\\ \sigma\_h + \Sigma\_h \left(\frac{t - T'}{\tau}\right)^2 \text{ for } T' \le t & \text{(the deviatoric stage)} \end{cases}$$

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.

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

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

**7. Comparison of numerical and analytical solution for underground** 

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

(a) (b)

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

Assuming further that the underground opening is at a certain depth where the horizontal and vertical components of the primary (initial) stress *<sup>h</sup>* and *<sup>v</sup>* are 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

Cylindrical coordinate system is chosen for convenience with axis Oz being the symmetry

respectively

**openings** 

whole domain.

considered.

**7.1. Problem formulation** 

axis of the opening (Figure 4).

infinite length, excavated in an infinite space.

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**

312 Finite Element Analysis – New Trends and Developments

solution.

smaller loading rates.

loading steps.

(see [22]).

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

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

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

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

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

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

(a) (b)

results for the two tests mentioned above, respectively.

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 *<sup>h</sup>* and *<sup>v</sup>* are 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:

$$\begin{aligned} \left\{ \begin{array}{c} r = a \\ \text{any} \end{array} \right\} : \quad \sigma\_{rr}^{S} = p\_{\prime} \sigma\_{r0}^{S} = 0, \quad \begin{array}{c} r \to \infty \\ \text{any} \end{array} \right\} : \quad \sigma\_{\text{xx}/w}^{S} = \sigma\_{h\prime} \sigma\_{yy/w}^{S} = \sigma\_{v\prime} \tag{11} \end{aligned} \tag{11}$$

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

θ θ

*rz z*

 

*u u uu r r rr*

, ; <sup>2</sup>

1 1

*r r*

z

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

*r r*

 

*z z z r*

2 2 0,

*z zz rz zz*

 

1 1 <sup>1</sup> , , . <sup>2</sup>

*uu u u u rr r z z*

 

*z*

<sup>1</sup> 0.

*z*

components of the small strains, zero values are founded for the relative strains

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

It is assumed that the opening (or only a part of it) is excavated in a very short time interval <sup>0</sup> *t t* 0, and then it is exploited in a much longer time interval 0 1 *t tt*, . 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

Consequently, if the tunnel is excavated at <sup>0</sup>*t* , then immediately after excavation, the stress, the strain and the displacement of the rock are given by the elastic solution. Concerning the second interval 0 1 *t t*, 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

*z u*

*r*

*z*

0,

0,

and then, according to the

 

 

*rz r z rz r rr*

*r r rz r r rz z* 

*r rz rr z z*

2 2 2 2

*r r rr r z rz z r*

2 2 2 2

*z rz r z zz zz*

 

2

22 2 2

*rz z r z z r*

θ θ θθ θ zz

 

*u u z r*

; <sup>2</sup>

*r z rz*

θ

rr r

1

*r r*

2 2

*r θ r r*

2 0,

*rr <sup>r</sup> rr*

2 2

2

2 22 2

*rr r*

<sup>2</sup> <sup>2</sup> <sup>2</sup>

2 0,

*rz zz rr*

2 2

<sup>1</sup> <sup>2</sup>

Due to the plane strain hypothesis, we have 0, 0, *<sup>R</sup>*

2

*r rr*

 0 . *RRR zz rz z*

stress.

analysis.

**7.2. Analytical solution** 

2 22

*z r r z*

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

$$\begin{split} \boldsymbol{\sigma}\_{rr/\infty}^{S} &= \boldsymbol{\sigma}\_{rr}^{P} = \frac{1}{2} (\boldsymbol{\sigma}\_{h} + \boldsymbol{\sigma}\_{v}) + \frac{1}{2} (\boldsymbol{\sigma}\_{h} - \boldsymbol{\sigma}\_{v}) \cos 2\boldsymbol{\theta}\_{r} \\ \boldsymbol{\sigma}\_{00/\infty}^{S} &= \boldsymbol{\sigma}\_{00}^{P} = \frac{1}{2} (\boldsymbol{\sigma}\_{h} + \boldsymbol{\sigma}\_{v}) - \frac{1}{2} (\boldsymbol{\sigma}\_{h} - \boldsymbol{\sigma}\_{v}) \cos 2\boldsymbol{\theta}\_{r} \\ \boldsymbol{\sigma}\_{r0/\infty}^{S} &= \boldsymbol{\sigma}\_{r0}^{P} = -\frac{1}{2} (\boldsymbol{\sigma}\_{h} - \boldsymbol{\sigma}\_{v}) \sin 2\boldsymbol{\theta}\_{r} \end{split}$$

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:

$$\begin{split} & \frac{\hat{\mathcal{O}}\boldsymbol{\sigma}\_{rr}}{\hat{\mathcal{O}}r} + \frac{1}{r} \frac{\hat{\mathcal{O}}\boldsymbol{\sigma}\_{r0}}{\hat{\mathcal{O}}\boldsymbol{\Theta}} + \frac{\hat{\mathcal{O}}\boldsymbol{\sigma}\_{rz}}{\hat{\mathcal{O}}z} + \frac{\boldsymbol{\sigma}\_{rr} - \boldsymbol{\sigma}\_{\boldsymbol{\Theta}0}}{r} = 0; \\ & \frac{\hat{\mathcal{O}}\boldsymbol{\sigma}\_{r\boldsymbol{\Theta}}}{\hat{\mathcal{O}}r} + \frac{1}{r} \frac{\hat{\mathcal{O}}\boldsymbol{\sigma}\_{\boldsymbol{\Theta}0}}{\hat{\mathcal{O}}\boldsymbol{\Theta}} + \frac{\hat{\mathcal{O}}\boldsymbol{\sigma}\_{\boldsymbol{\Theta}z}}{\hat{\mathcal{O}}z} + 2 \frac{\boldsymbol{\sigma}\_{r\boldsymbol{\Theta}}}{r} = 0; \\ & \frac{\hat{\mathcal{O}}\boldsymbol{\sigma}\_{rz}}{\hat{\mathcal{O}}r} + \frac{1}{r} \frac{\hat{\mathcal{O}}\boldsymbol{\sigma}\_{z\boldsymbol{\Theta}}}{\hat{\mathcal{O}}\boldsymbol{\Theta}} + \frac{\hat{\mathcal{O}}\boldsymbol{\Theta}\_{zz}}{\hat{\mathcal{O}}z} + \frac{\boldsymbol{\sigma}\_{rz}}{r} = 0. \end{split}$$

The components of the small strains:

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

$$\begin{split} \varepsilon\_{\varepsilon \text{rr}} &= \frac{\widehat{\varepsilon}u\_{r}}{\widehat{\varepsilon}r}, \varepsilon\_{\text{r}0} = \frac{1}{2} \left[ \frac{1}{r} \frac{\widehat{\varepsilon}u\_{r}}{\widehat{\varepsilon}\Theta} + \frac{\widehat{\varepsilon}u\_{\text{o}}}{\widehat{\varepsilon}r} - \frac{u\_{\text{o}}}{r} \right]; \\ \varepsilon\_{\varepsilon \text{z}} &= \frac{1}{2} \left[ \frac{\widehat{\varepsilon}u\_{r}}{\widehat{\varepsilon}z} + \frac{\widehat{\varepsilon}u\_{\text{z}}}{\widehat{\varepsilon}r} \right]; \\ \varepsilon\_{\text{\varepsilon\text{0}}} &= \left[ \frac{1}{r} \frac{\widehat{\varepsilon}u\_{\text{o}}}{\widehat{\varepsilon}\Theta} + \frac{u\_{r}}{r} \right], \varepsilon\_{\text{x}0} = \frac{1}{2} \left[ \frac{1}{r} \frac{\widehat{\varepsilon}u\_{\text{z}}}{\widehat{\varepsilon}\Theta} + \frac{\widehat{\varepsilon}u\_{\text{o}}}{\widehat{\varepsilon}z} \right], \varepsilon\_{\text{xx}} = \frac{\widehat{\varepsilon}u\_{\text{z}}}{\widehat{\varepsilon}z}. \end{split}$$

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

 <sup>2</sup> <sup>2</sup> <sup>2</sup> 2 2 2 22 2 2 2 2 22 2 2 2 2 2 2 2 2 0, 2 2 0, 2 0, <sup>1</sup> <sup>2</sup> *rr <sup>r</sup> rr z zz rz zz rz zz rr rz r z rz r rr z r r r r r θ r r rr r z z z r z r r z r r rz r r rz z* 2 2 2 2 2 22 2 2 2 0, 0, <sup>1</sup> 0. *z r rz rr z z z rz r z zz zz r r r rr r z rz z r r rr rz z r z z r*

Due to the plane strain hypothesis, we have 0, 0, *<sup>R</sup> z u z* and then, according to the components of the small strains, zero values are founded for the relative strains 0 . *RRR zz rz z*

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

#### **7.2. Analytical solution**

314 Finite Element Analysis – New Trends and Developments

**Figure 4.** The domain of the boundary value problem

any

displacement and strain respectively.

The components of the small strains:

*r a*

different causes. Therefore, the boundary conditions are:

 : , 0,

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

/ 2

Further, the fundamental equations of the problem are presented.

The equilibrium equations written in the relative components are:

 

*S P*

*S P*

*S P*

*S S rr r*

*p*

A pressure loading *p(t)* is considered on the surface *r = a* of the opening, that could be due to

*z* 

*z*

 

*hv hv*

sin 2 .

any

The superscripts *S, P, R* mean the secondary, primary and relative components of the stress,

<sup>1</sup> 0;

<sup>1</sup> 2 0;

*z r z r z*

<sup>1</sup> 0.

*z*

 

*r r r*

*rr r rz rr*

*rr r*

*rz z zz rz*

*r r r*

 

 

1 1

 

1 1

*rr rr h v h v*

*r*

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

*r r hv*

/ 2 2 1

/ 2 2

/ / : ,,

cos2 , cos2 , (11)

*S S xx h yy v*

> It is assumed that the opening (or only a part of it) is excavated in a very short time interval <sup>0</sup> *t t* 0, and then it is exploited in a much longer time interval 0 1 *t tt*, . 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 at <sup>0</sup>*t* , then immediately after excavation, the stress, the strain and the displacement of the rock are given by the elastic solution. Concerning the second interval 0 1 *t t*, 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 0 1 *t t*, the stress components , , *rr r* remain constant and equal with those given by the elastic solution (see [1]).

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

In the case of numerical approach, the domain to be discretized and the boundary

In this case a quarter of the domain needs to be considered as the principal components of

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

In the following examples the lining was considered made of concrete with Young modulus

There were considered the cases when *<sup>h</sup>* = *<sup>v</sup>* = 2000 kPa and the case <sup>h</sup> *<sup>v</sup>* , namely *<sup>h</sup>* =

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

In the model under consideration of saturated sand it was observed that values for *k t* which exceeds the magnitude order of <sup>3</sup> 10 produce 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

the primary stress *<sup>h</sup>* and *<sup>v</sup>* are assumed constant over the entire domain.

characteristics of the lining material (Figure 6).

E = 200000 kPa and Poisson coefficient = 0.3.

1500 kPa and *<sup>v</sup>* = 2000 kPa. The tunnel radius is 1 m.

tangential operator, backward Euler method, etc.).

value of *k* d*t* to ensure a good convergence of the numerical scheme.

**7.3. The numerical solution** 

**Figure 6.** FEM mesh

conditions, see Figure 5 and the FEM mesh, see Figure 6, are:

Due to the fact of the plane strain hypothesis 0 *<sup>R</sup> zz* it has to be accepted that *zz* varies during this interval and therefore it satisfies the differential equation:

$$\frac{1}{3} \left( \frac{1}{3K} + \frac{1}{G} \right) \dot{\sigma}\_{zz} = -k \left( 1 - \frac{W^I}{H} \right) \frac{\partial F}{\partial \sigma\_{zz}}.$$

This differential equation may be integrated numerically using the initial conditions zz *<sup>h</sup>* . It was found that a very small variation in time (negligible) of zz was 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 *W t <sup>I</sup>* when all stress components are constant is used.

The time taken for *W t <sup>I</sup>* 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:

**Figure 5.** Domain used in numerical formulation

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

## **7.3. The numerical solution**

316 Finite Element Analysis – New Trends and Developments

Due to the fact of the plane strain hypothesis 0 *<sup>R</sup>*

observing the convergence of a tunnel wall in time.

during this interval and therefore it satisfies the differential equation:

components could be constant during the rock creep around the opening.

supplies the variation of *W t <sup>I</sup>* when all stress components are constant is used.

3 3

solution (see [1]).

initial conditions (9).

0

*R*

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 0 1 *t t*, the stress components , , *rr r* remain constant and equal with those given by the elastic

> 11 1 1 .

This differential equation may be integrated numerically using the initial conditions zz *<sup>h</sup>* . It was found that a very small variation in time (negligible) of zz was exhibited and a very fast stabilization of this variation. Therefore, it can be assumed that all stress

In order to arrive at the expression that describes the creep strain, the equation (7) that

The time taken for *W t <sup>I</sup>* 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

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

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

 

0

*H k F u u ttr F H*

<sup>1</sup> , 1 exp <sup>1</sup>

 

*IP*

*H*

**Figure 5.** Domain used in numerical formulation

*W F*

*K G H*

*zz*

*I*

*W F <sup>k</sup>*

*zz*

*zz* it has to be accepted that *zz* varies

with **u**R 0.

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 *<sup>h</sup>* and *<sup>v</sup>* are 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 *<sup>h</sup>* = *<sup>v</sup>* = 2000 kPa and the case <sup>h</sup> *<sup>v</sup>* , namely *<sup>h</sup>* = 1500 kPa and *<sup>v</sup>* = 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* d*t* to ensure a good convergence of the numerical scheme.

In the model under consideration of saturated sand it was observed that values for *k t* which exceeds the magnitude order of <sup>3</sup> 10 produce 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.

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

**Figure 8.** Contours of

*ij*

**Figure 9.** Contours of *rr* after the first time step

10a for the first time step and 10b for the 130th time step.

For the same boundary problem the evolution of the hoop stress is designed in Figures

*vp* after the first time step

## *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 *<sup>h</sup>* = *<sup>v</sup>* = 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 *ij vp* for the first time step and after 130 time steps with one time step 10000 s, respectively.

**Figure 7.** a), b). Contours of *ij vp* 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 *<sup>h</sup>* = 1500 kPa and *<sup>v</sup>* = 2000 kPa, exhibits for instance, a development of a viscoplastic zone in the tunnel wall that increases from 1.13E-4 for *ij vp* after the first time step Figure 8 to 5.13E-3 for *ij vp* after 130 time steps.

The radial stress component *rr* presents a very slow variation in time and it is presented in Figure 9.

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

**Figure 8.** Contours of *ij vp* after the first time step

318 Finite Element Analysis – New Trends and Developments

represented by the analytical solution.

**Figure 7.** a), b). Contours of *ij*

that increases from 1.13E-4 for *ij*

130 time steps.

Figure 9.

*ij*

walls *p* = 1000 kPa which simulates the lining.

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

The following Figures represents the numerical solution for the case of a tunnel subjected to a hydrostatic primary stress of *<sup>h</sup>* = *<sup>v</sup>* = 2000 kPa and an internal pressure on the tunnel

Figures 7a and 7b represent the evolution of the Euclidian norm of the viscoplastic strain

*vp* for the first time step and after 130 time steps with one time step 10000 s, respectively.

*vp* evolution of tunnel with an internal pressure in the case of a

*vp* after the first time step Figure 8 to 5.13E-3 for *ij*

*vp* after

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 *<sup>h</sup>* = 1500 kPa and *<sup>v</sup>* = 2000 kPa, exhibits for instance, a development of a viscoplastic zone in the tunnel wall

(a) (b)

The radial stress component *rr* presents a very slow variation in time and it is presented in

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

**Figure 9.** Contours of *rr* after 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.

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

**Figure 11.** a), b) Contours of total *u* evolution in a supported tunnel excavated in a rock mass with non-

(a) (b)

*vp* evolution 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

hydrostatic primary stress, after the first time step and after the130th time step, respectively

(a) (b)

**Figure 12.** a), b) Contours of *ij*

**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 *<sup>h</sup>* = 1500 kPa and *<sup>v</sup>* = 2000 kPa the diagrams of variation are presented as follows: the total displacement 2 2 total *x y u uu* in Figures 11a and 11b, and the norm of the viscoplastic strain *ij vp* 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 10.** a) Contours after the first time step, b) Contours after the 130th time step

In this subparagraph, the numerical solution for a concrete lining as a distinct material and

(a) (b)

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 *<sup>h</sup>* = 1500 kPa and *<sup>v</sup>* = 2000 kPa the diagrams of variation are presented as follows: the total

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.

total *x y u uu* in Figures 11a and 11b, and the norm of the viscoplastic strain

*7.3.2. Numerical solution for the case concrete lining* 

mesh elements group is presented.

displacement 2 2

*vp* in Figures 12a and 12b respectively.

*ij*

**Figure 11.** a), b) Contours of total *u* evolution in a supported tunnel excavated in a rock mass with nonhydrostatic primary stress, after the first time step and after 130th time step, respectively

**Figure 12.** a), b) Contours of *ij vp* evolution in a supported tunnel excavated in a rock mass with nonhydrostatic 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.

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

**Figure 14.** The time evolution of viscoplastic strain of tunnel walls for numerical and analytical

analytical solution, assumption that could not be considered totally unrealistic.

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 elastoviscoplastic 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

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

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

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

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

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.

solutions

applications.

exhibiting viscoplastic behaviour too.

tendency to increase in time.

solution.

**7.5. Final remarks** 

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 *<sup>h</sup>* = *<sup>v</sup>* = 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 *<sup>h</sup>* = 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**

322 Finite Element Analysis – New Trends and Developments

**7.4. Comparisons between the numerical and analytical solutions** 

steps (i.e. more than 10 time steps), the analytical solution becomes larger.

important information both quantitatively and qualitatively.

some constructive conclusions for both solutions.

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

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

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 *<sup>h</sup>* = *<sup>v</sup>* = 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

**Figure 13.** a), b). Tunnel outline total displacement with respect to distance and time, respectively, for

(a) (b)

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 *<sup>h</sup>*

= 1500 kPa, *v* = 2000 kPa, and the internal pressure *p* = 1000 kPa.

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

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

**Figure 15.** Domain and boundary conditions for the problem study in O*rz* plane along the tunnel axis.

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

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 rockmass 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

for the concrete lining.

gradually.

lining interaction.

above.

increases with the advancing of the tunnel face.

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

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

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

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-

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

h of the

the tunnel face influence that means a gradual decompression of the primary stress

## **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, **σ** *<sup>P</sup> PI* where *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 *<sup>v</sup>* , *<sup>h</sup>* are assumed equal), the problem is an *axisymmetrical* one in O*rz* plane (Figure 15).

Consequently, the primary stress components *<sup>v</sup>* , *<sup>h</sup>* are assumed equal. The boundary conditions are:

On . . [ , ] and : and 0. *AB i e z z z r a p <sup>A</sup> <sup>B</sup> rr rz* On . . [ , ]and : 0 and 0. *B C rr rz BC i e z z z r a* On . . and [0, ] : 0 and 0. *CD i e z z z r a C D rr rz* On . . [ , ] and 0 : 0 and 0. *DE i e z z z r u D E <sup>r</sup> rz* On . . and [0, ] : and 0. *E F F zz v rz EF i e z z z r r* On . . and : and 0. *F G F rr v rz FG i e z z z r r* On . . and [ , ] : and 0. *GA i e z z r r r A A G zzv rz*

into account the 3D aspect of the problem.

characteristics represents no difficulty.

**8.1. Formulation of the problem** 

*axisymmetrical* one in O*rz* plane (Figure 15).

On . . [ , ] and : and 0. *AB i e z z z r a p <sup>A</sup> <sup>B</sup> rr rz*

On . . [ , ]and : 0 and 0. *B C rr rz BC i e z z z r a* On . . and [0, ] : 0 and 0. *CD i e z z z r a C D rr rz* On . . [ , ] and 0 : 0 and 0. *DE i e z z z r u D E <sup>r</sup> rz* On . . and [0, ] : and 0. *E F F zz v rz EF i e z z z r r* On . . and : and 0. *F G F rr v rz FG i e z z z r r* On . . and [ , ] : and 0. *GA i e z z r r r A A G zzv rz*

tensor.

conditions are:

**mounting problem** 

All these features are in good accordance with the practical observation, so, more involved

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

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

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

It is convenient to lead calculation to a 2D, plane strain or axisymmetrical, if it is possible,

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, **σ** *<sup>P</sup> PI* where *P = h*, *h* is the

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 *<sup>v</sup>* , *<sup>h</sup>* are assumed equal), the problem is an

Consequently, the primary stress components *<sup>v</sup>* , *<sup>h</sup>* are assumed equal. The boundary

is the unity

since it is less costly in data input and running time than a 3D analysis.

depth at which the tunnel is dug, is the specific gravity of the rock and

boundary problems could be envisaged to be solved with this code.

**Figure 15.** Domain and boundary conditions for the problem study in O*rz* 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 rockmass 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:


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

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



In the following, we present some important results of the calculation concerning, for

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

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

(a) (b) (c)

Isovalues zones for normal stress corresponding to the three phases are presented in Figures

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

18a, 18b, 18c. Great concentrations are observed in tunnel face zone (white area).

and output storage for the following phase (phase 3) are performed as well.

instance, the normal stress, the equivalent stress or the damage parameter *df*.

state is performed as well.

maximum (white area).

very low for the rocks.

## **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:


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


326 Finite Element Analysis – New Trends and Developments

20000 MPa, Poisson coefficient

lining pressure and the convergence of the tunnel walls.

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

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

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* =

It is used a numbering of elements group as they are activated/deactivated in the excavation




Numerical model concerns three successive phases' tunnel excavation and lining mounting.


third phase and eventually replaced by concrete in a possible fourth phase

= 0.3 and the volumetric weight = 0.02 MN/m2.

**8.2. The numerical solution of the successive phases tunnel excavation** 

mesh, at least 2 layers of quadrilateral elements in the concrete lining.

**Figure 16.** The domain discretization of the boundary value problem.

and lining installing processes, as follows:

third phase


second phase is replaced by the concrete

Let us detail the progression of phases of the example:



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.

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

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

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

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

The previous calculation was performed for a value of Young modulus *E*=20000 Mpa,

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

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

(a) (b)

*Eb c*

2 2

= 0.3, external radius *b* = 1.2 m and internal radius *c* = 1 m. We perform

*b c*

2 2

( )

(1 ) 1 2

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

*q*

**mounting time** 

Poisson coefficient

*8.3.1. Support rigidity influence* 

rigidity and depth at which the tunnel is excavated.

solution by several Figures and observations.

obtained by author by analytical means in [20], [18].

**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*

328 Finite Element Analysis – New Trends and Developments

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

(a) (b) (c)

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

(a) (b) (c)

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

$$q = \frac{E(b^2 - c^2)}{(1 + \upsilon)\left[ (1 - 2\upsilon)b^2 + c^2 \right]}$$

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

330 Finite Element Analysis – New Trends and Developments

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

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

Figures 22a, b present isovalues zones of the normal stress *rr* for 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

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

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

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

[1] Cristescu N (1989) *Rock rheology*, Kluwer Academic Publishers, Dordrecht, Holland. [2] Simionescu O (1999) *Mathematical methods in underground structure termomechanics,*

*numérique.* Thése. École National des Ponts et Chaussées, France.

[3] Bernaud D (1991) *Tunnels profonds dans les milieux viscoplastique: approches expérimentale et* 

developed with this code and further improvements for the analytical solution.

(the white zone is more extended).

author by analytical means [19], [20].

*Military Technical Academy, Bucharest, Romania* 

Romanian Academy Publ., Bucharest.

**Author details** 

**10. References** 

Simona Roatesi

**9. Conclusions** 

obtained.

depth, which may induces rock or lining fracture.

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

**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 *rr* for 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**

330 Finite Element Analysis – New Trends and Developments

*8.3.2 Tunnel depth influence* 

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

(a) (b) (c)

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

(a) (b)

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

## **Author details**

Simona Roatesi *Military Technical Academy, Bucharest, Romania* 

## **10. References**

	- [4] Panet M (1974) *Stabilité et soutenement des tunnels. La méc.des roches appl. aux ouvrages de génie civil*,Ch.IX,ENPC.

**Chapter 15** 

© 2012Santos da Silva et al., licensee InTech. This is an open access chapter distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/3.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Attribution License (http://creativecommons.org/licenses/by/3.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

© 2012Santos da Silva et al., licensee InTech. This is a paper distributed under the terms of the Creative Commons

**Finite Element Modelling of** 

and Luciano Rodrigues Ornelas de Lima

http://dx.doi.org/10.5772/50376

technologic and design issues.

expected footbridge service life.

**1. Introduction** 

solutions.

Additional information is available at the end of the chapter

**the Dynamic Behaviour of Tubular Footbridges** 

Tubular hollow sections are increasingly used in off-shore structures, highway bridges, pedestrian foot-bridges, large-span roofs and multi-storey buildings due to their excellent properties and the associated advances in fabrication technology. The intensive use of tubular structural elements in Brazil, such as the example depicted in Figure 1, mainly due to its associated aesthetical and structural advantages, led designers to be focused on their

Nowadays in Brazil, there is still a lack of code that deals specifically with tubular design. This fact induces designers to use other international tubular design codes. Consequently, their design methods accuracy plays a fundamental role when economical and safety points of view are considered. Additionally, recent tubular joint studies indicate that further research is needed, especially for particular geometries. This is even more significant for some failure modes where the collapse load predictions lead to unsafe or uneconomical

Steel and composite tubular footbridges are currently subjected to dynamic actions with variable magnitudes due to the pedestrian crossing on the concrete deck [1-4]. These dynamic actions can generate the initiation of fractures or even their propagation in the structure. Depending on the magnitude and intensity, these adverse effects can compromise the structural system response and the reliability which may also lead to a reduction of the

Generally, fatigue assessment procedures are usually based on S-N curves which relate a nominal or geometric stress range S to the corresponding number N of load cycles to fatigue

José Guilherme Santos da Silva, Ana Cristina Castro Fontenla Sieira, Gilvan Lunz Debona, Pedro Colmar Gonçalves da Silva Vellasco


## **Finite Element Modelling of the Dynamic Behaviour of Tubular Footbridges**

José Guilherme Santos da Silva, Ana Cristina Castro Fontenla Sieira, Gilvan Lunz Debona, Pedro Colmar Gonçalves da Silva Vellasco and Luciano Rodrigues Ornelas de Lima

Additional information is available at the end of the chapter

http://dx.doi.org/10.5772/50376

## **1. Introduction**

332 Finite Element Analysis – New Trends and Developments

*Num. Meth.Engng*., Vol.9, pp.109−127.

*génie civil*,Ch.IX,ENPC.

Univ. Press Oxford.

Vol. 7, Issue 3, pp. 295-307.

in Geomechanics, 25: 243–262.

Modeling, WMSO08,pp.334–337.

Bulletin, UTCB, 1, pp.16-30.

Academy Publisher, Bucharest.

J.Engng.Mech., 113, pp. 1302-1318.

pp.155−164.

353.

1999.

325.

37-52.

[4] Panet M (1974) *Stabilité et soutenement des tunnels. La méc.des roches appl. aux ouvrages de* 

[5] Sulem J, Panet M, Guenot A (1987) *An analytical solution for time-dependent displacements in a circular tunnel*, Int.J. Rock Mech. Min. Sci.& Geomech. Abstr., Vol.24, No.3,

[6] Zienkiewicz O C, Cormeau I C (1974) *Visco-plasticity - plasticity and creep in elastic solids -* 

[7] Cormeau, I.C., Numerical stability in quasi-static elasto/visco-plasticity (1975) *Int. J.* 

[8] Ionescu I R, Sofonea M (1993) *Functional and numerical methods in viscoplasticity,* Oxford

[9] Chan A H, Ou J (2008) *Three-Dimensional Numerical Analysis of a Dynamic Structure, Saturated Soil and Pore Fluid Interaction Problem*, Trends in Eng. Comput. Techn., pp. 335-

[10] Zienkiewicz O C, Chan A H C, Pastor M, Schrefler B A, Shiomi T (1999) *Computational Geomechanics with Special Refernces to Earthquake Engineering*, J. Wiley & Sons, Chicester,

[11] Gioda G (1993) *Finite element analysis of time dependent effects in tunnels*, Advanced school

[12] Lauro F, Bennani B, Drazetic P, Oudin J, Ni X (1997) *Ductile damage and fracture finite element modelling of elasto-viscoplastic voided materials*, Computational Materials Science,

[13] Lin R C, Brocks W (2006) *On a finite-strain viscoplastic law coupled with anisotropic damage: theoretical formulations and numerical applications*, Archive of Appl.Mech.,Vol.75, 6-7, 315-

[14] Augarde C E, Burd H J (2001) *3D FE analysis of lined tunnels*, Int. J. Num. Analyt. Meth.

[15] Komiya K, Soga K, Akagi H, Hagiwara T, Bolton M D (1999) *FE, modeling of excavation and advancement process of a shield tunnelling machine*, Soil &Foundations, Vol. 39, No.3,

[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. [17] Wang J, Wang W (2009) *Multi-field Coupled Model&Num.Sim. During Excav.Tunnel* 

[18] Roateşi S, Ariciuc M (2005), *Study on depth dependency in tunnel calculation*, Scientifical

[19] Roateşi S (1996) *The tunnel face influence in the analysis of a circular tunnel with a time-*

[20] Roateşi S (2005) *Mathematical modeling and FEM problems in tunnel calculation*, Romanian

[22] Lade P V, Nelson R B, Ito Y M (1987) *Nonassociated flow and stability of granular materials*,

*dependent behaviour*, Rev.Roum.Sci.Tech. Mec.Appl., no 3-4, pp. 247-263.

[21] CESAR − LCPC *Documentation*, Version 3.0.x, Publication LCPC, 1992.

on visco-plastic behaviour of geomaterials, CISM, Udine, Oct.

*A unified numerical solution approach*, Int.J.Num.Meth.Engng.,8,pp. 821-845.

Tubular hollow sections are increasingly used in off-shore structures, highway bridges, pedestrian foot-bridges, large-span roofs and multi-storey buildings due to their excellent properties and the associated advances in fabrication technology. The intensive use of tubular structural elements in Brazil, such as the example depicted in Figure 1, mainly due to its associated aesthetical and structural advantages, led designers to be focused on their technologic and design issues.

Nowadays in Brazil, there is still a lack of code that deals specifically with tubular design. This fact induces designers to use other international tubular design codes. Consequently, their design methods accuracy plays a fundamental role when economical and safety points of view are considered. Additionally, recent tubular joint studies indicate that further research is needed, especially for particular geometries. This is even more significant for some failure modes where the collapse load predictions lead to unsafe or uneconomical solutions.

Steel and composite tubular footbridges are currently subjected to dynamic actions with variable magnitudes due to the pedestrian crossing on the concrete deck [1-4]. These dynamic actions can generate the initiation of fractures or even their propagation in the structure. Depending on the magnitude and intensity, these adverse effects can compromise the structural system response and the reliability which may also lead to a reduction of the expected footbridge service life.

Generally, fatigue assessment procedures are usually based on S-N curves which relate a nominal or geometric stress range S to the corresponding number N of load cycles to fatigue

© 2012Santos da Silva et al., licensee InTech. This is an open access chapter distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/3.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited. © 2012Santos da Silva et al., licensee InTech. This is a paper distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/3.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

failure. In this situation fatigue assessment refers to the nominal stress range Δσ in a tubular structural member.

Finite Element Modelling of the Dynamic Behaviour of Tubular Footbridges 335

On the other hand, the structural engineers experience and knowledge allied by the use newly developed materials and technologies have produced tubular steel and composite (steel-concrete) footbridges with daring structures. This fact have generated very slender tubular steel and composite pedestrian footbridges and consequently changed the serviceability and ultimate limit states associated to their design. A direct consequence of

Considering all aspects mentioned before, the main objective of this investigation is to present the finite element modelling of the dynamic behaviour of tubular composite (steelconcrete) footbridges submitted to human walking vibration. Based on the results obtained in this study, a fatigue assessment will be performed, in order to evaluate the tubular

The investigated structural model was based on a tubular composite (steel-concrete) footbridge, spanning 82.5 m. The structure is composed by three spans (32.5 m, 17.5 m and 20.0 m, respectively) and two overhangs (7.50 m and 5.0 m, respectively). The structural system consists of tubular steel sections and a concrete slab and is currently used for

The proposed computational model adopted the usual mesh refinement techniques present in finite element method simulations, based on the ANSYS program [13]. The finite element model has been developed and validated with the experimental results. This numerical model enabled a complete dynamic evaluation of the investigated tubular footbridge especially in terms of human comfort and its associated vibration serviceability limit states.

footbridges service life. Further research in this area is currently being carried out.

this design trend is a considerable increase of structural vibrations [1-4, 8-11].

**Figure 2.** Typical multiplanar K-joint with notations.

pedestrian crossing [1-2].

**Figure 1.** Example of a tubular steel pedestrian footbridge in Rio de Janeiro/RJ, Brazil.

The fatigue resistance is given according to a classification catalogue in the form of standardized S-N curves. Structural details classified in this catalogue, see e.g. Eurocode 3 Part 1.9 [5], correspond specifically to a situation of stress range, direction, crack position, detail dimension and weld quality which had been characteristic for the tests on which the classification is based [6-7].

The use of circular hollow section members as part of the structure of pedestrian footbridges is a relatively new constructional concept. During the last couple years several steel-concrete composite footbridges had been constructed in Brazil, as illustrated in Figure 1.

The typical cross-section of this type of pedestrian footbridge generally consists of a tubular spatial truss girder carrying the concrete deck slab, as presented in Figure 1. The deck slab is connected directly to the steel structure by either shear studs, concrete dowels or in some cases where no top chord exists. At the bottom chord of the tubular space truss four brace members have to be connected to the continuous bottom chord. This type of joint is usually named K-joint, as depicted in Figure 2 [6].

Steel and composite tubular footbridges can be subjected to the material imperfections of its structural elements, such as mechanical and metallurgic discontinuities. Such defects lead to cracking in the-se structural elements. When these elements are subjected to dynamic actions, the fatigue phenomenon occurs and produces stress concentrations and possible fractures. These fractures are directly responsible for reducing the local or global footbridge stabilities or even its life service [7].

On the other hand, the structural engineers experience and knowledge allied by the use newly developed materials and technologies have produced tubular steel and composite (steel-concrete) footbridges with daring structures. This fact have generated very slender tubular steel and composite pedestrian footbridges and consequently changed the serviceability and ultimate limit states associated to their design. A direct consequence of this design trend is a considerable increase of structural vibrations [1-4, 8-11].

**Figure 2.** Typical multiplanar K-joint with notations.

334 Finite Element Analysis – New Trends and Developments

structural member.

classification is based [6-7].

named K-joint, as depicted in Figure 2 [6].

stabilities or even its life service [7].

failure. In this situation fatigue assessment refers to the nominal stress range Δσ in a tubular

**Figure 1.** Example of a tubular steel pedestrian footbridge in Rio de Janeiro/RJ, Brazil.

composite footbridges had been constructed in Brazil, as illustrated in Figure 1.

The fatigue resistance is given according to a classification catalogue in the form of standardized S-N curves. Structural details classified in this catalogue, see e.g. Eurocode 3 Part 1.9 [5], correspond specifically to a situation of stress range, direction, crack position, detail dimension and weld quality which had been characteristic for the tests on which the

The use of circular hollow section members as part of the structure of pedestrian footbridges is a relatively new constructional concept. During the last couple years several steel-concrete

The typical cross-section of this type of pedestrian footbridge generally consists of a tubular spatial truss girder carrying the concrete deck slab, as presented in Figure 1. The deck slab is connected directly to the steel structure by either shear studs, concrete dowels or in some cases where no top chord exists. At the bottom chord of the tubular space truss four brace members have to be connected to the continuous bottom chord. This type of joint is usually

Steel and composite tubular footbridges can be subjected to the material imperfections of its structural elements, such as mechanical and metallurgic discontinuities. Such defects lead to cracking in the-se structural elements. When these elements are subjected to dynamic actions, the fatigue phenomenon occurs and produces stress concentrations and possible fractures. These fractures are directly responsible for reducing the local or global footbridge Considering all aspects mentioned before, the main objective of this investigation is to present the finite element modelling of the dynamic behaviour of tubular composite (steelconcrete) footbridges submitted to human walking vibration. Based on the results obtained in this study, a fatigue assessment will be performed, in order to evaluate the tubular footbridges service life. Further research in this area is currently being carried out.

The investigated structural model was based on a tubular composite (steel-concrete) footbridge, spanning 82.5 m. The structure is composed by three spans (32.5 m, 17.5 m and 20.0 m, respectively) and two overhangs (7.50 m and 5.0 m, respectively). The structural system consists of tubular steel sections and a concrete slab and is currently used for pedestrian crossing [1-2].

The proposed computational model adopted the usual mesh refinement techniques present in finite element method simulations, based on the ANSYS program [13]. The finite element model has been developed and validated with the experimental results. This numerical model enabled a complete dynamic evaluation of the investigated tubular footbridge especially in terms of human comfort and its associated vibration serviceability limit states.

This investigation is carried out based on correlations between the experimental results related to the footbridge dynamic response and those obtained with finite element models [1-2]. The structural system dynamic response, in terms of peak accelerations, was obtained and compared to the limiting values proposed by several authors and design standards [9,14].

Finite Element Modelling of the Dynamic Behaviour of Tubular Footbridges 337

This walking load model can be represented by the load static parcel, related to the individual weight, and a combination of harmonic forces whose frequencies are multiples or harmonics of the basic frequency of the force repetition, e.g. step frequency, fs, for human activities. This load model considers a space and temporal variation of the dynamic action over the structure and the time-dependent repeated force can be represented by the Fourier

> ( ) [1 cos 2 ]

In this load model, five harmonics were considered to represent the dynamic load associated to human walking [12]. Table 2 shows the dynamic coefficients and phase angles used in this load model. Figure 3 presents a dynamic loading function for one person walking with

Harmonic i Dynamic Coefficients α<sup>i</sup> Phase Angles φ<sup>i</sup> 1 0.37 0 2 0.10 π/2 3 0.12 π/2 4 0.04 π/2 5 0.08 π/2

 *i si Ft P if t* (1)

**2.1. Load model I (LM-I)** 

series, see Equation (1).

F(t) : dynamic load;

P : person's weight (800 N [1-4]);

fs : walking step frequency; : harmonic phase angle;

step frequency equal to 2 Hz.

i : harmonic multiple (i = 1,2,3…,n);

**Table 2.** Dynamic coefficients and phase angles [12].

**Figure 3.** LM-I: dynamic load function for one person walking (fs = 2.0 Hz).

αi : dynamic coefficient for the harmonic force;

Where:

t : time.

The peak acceleration values found in the present investigation indicated that the analysed tubular footbridge presented problems related with human comfort. Hence it was detected that this type of structure can reach high vibration levels that can compromise the footbridge user's comfort and especially its safety.

## **2. Human walking modelling**

Human loads comprise a large portion of the acting live loads in offices and residential building floors. In general, the human live loads are classified into two broad categories: in situ and moving. Periodic jumping due to music, sudden standing of a crowd, and random in-place movements are examples of in situ activities. Walking, marching, and running are examples of moving activities. As the main purpose of footbridges is the pedestrian's crossing, they must be safe and do not cause discomfort to users [1-4].

On the other hand, human activities produce dynamic forces and their associate vibration levels should not disturb or alarm their users. Therefore, this investigation describes four different load models developed to incorporate the dynamic effects induced by people walking on the footbridges dynamic response. It must be emphasized that the geometry of the human body walking is an organized leg motion that cause an ascent and descending movement of the effective body mass at each step [1-4].

The human body mass accelerations are associated to floor reactions, and are approximately periodic to the step frequency. The two feet produce this type of load, as function of the static parcel associated to the individual weight and three or four harmonic load components. These harmonics appear due to the interaction between the increasing loads, represented by one foot, and the simultaneous unload of the other foot [1-4].

However, it is also necessary to incorporate several other parameters in the human walking representation, like step distance and speed. These variables are related to the step frequency and are depicted in Table 1 [12]. Table 1 presents a detailed description of the excitation frequency values, dynamic coefficients, as well as the phase angles to be employed in the mathematical representation of the four dynamic loading models implemented and used in the present investigation.


**Table 1.** Characteristics of the human walking [12].

#### **2.1. Load model I (LM-I)**

This walking load model can be represented by the load static parcel, related to the individual weight, and a combination of harmonic forces whose frequencies are multiples or harmonics of the basic frequency of the force repetition, e.g. step frequency, fs, for human activities. This load model considers a space and temporal variation of the dynamic action over the structure and the time-dependent repeated force can be represented by the Fourier series, see Equation (1).

$$F(t) = P\{1 + \sum \alpha\_i \cos\left(2\,\pi\,\text{i}\,\,f\_s\,t + \phi\_i\right)\}\tag{1}$$

Where:

336 Finite Element Analysis – New Trends and Developments

footbridge user's comfort and especially its safety.

**2. Human walking modelling** 

This investigation is carried out based on correlations between the experimental results related to the footbridge dynamic response and those obtained with finite element models [1-2]. The structural system dynamic response, in terms of peak accelerations, was obtained and compared to the limiting values proposed by several authors and design standards [9,14].

The peak acceleration values found in the present investigation indicated that the analysed tubular footbridge presented problems related with human comfort. Hence it was detected that this type of structure can reach high vibration levels that can compromise the

Human loads comprise a large portion of the acting live loads in offices and residential building floors. In general, the human live loads are classified into two broad categories: in situ and moving. Periodic jumping due to music, sudden standing of a crowd, and random in-place movements are examples of in situ activities. Walking, marching, and running are examples of moving activities. As the main purpose of footbridges is the pedestrian's

On the other hand, human activities produce dynamic forces and their associate vibration levels should not disturb or alarm their users. Therefore, this investigation describes four different load models developed to incorporate the dynamic effects induced by people walking on the footbridges dynamic response. It must be emphasized that the geometry of the human body walking is an organized leg motion that cause an ascent and descending

The human body mass accelerations are associated to floor reactions, and are approximately periodic to the step frequency. The two feet produce this type of load, as function of the static parcel associated to the individual weight and three or four harmonic load components. These harmonics appear due to the interaction between the increasing loads,

However, it is also necessary to incorporate several other parameters in the human walking representation, like step distance and speed. These variables are related to the step frequency and are depicted in Table 1 [12]. Table 1 presents a detailed description of the excitation frequency values, dynamic coefficients, as well as the phase angles to be employed in the mathematical representation of the four dynamic loading models

Activity Velocity (m/s) Step Distance (m) Step Frequency (Hz)

Slow Walking 1.1 0.60 1.7 Normal Walking 1.5 0.75 2.0 Fast Walking 2.2 1.00 2.3

represented by one foot, and the simultaneous unload of the other foot [1-4].

crossing, they must be safe and do not cause discomfort to users [1-4].

movement of the effective body mass at each step [1-4].

implemented and used in the present investigation.

**Table 1.** Characteristics of the human walking [12].


In this load model, five harmonics were considered to represent the dynamic load associated to human walking [12]. Table 2 shows the dynamic coefficients and phase angles used in this load model. Figure 3 presents a dynamic loading function for one person walking with step frequency equal to 2 Hz.


**Table 2.** Dynamic coefficients and phase angles [12].

**Figure 3.** LM-I: dynamic load function for one person walking (fs = 2.0 Hz).

## **2.2. Load model II (LM-II)**

In this load model, the time-dependent repeated force also can be represented by the Fourier series, as expressed in Equation (1) and four harmonics were considered to represent the dynamic action associated to human walking [9]. This model also considers a space and temporal variation of the dynamic action over the structural system. Table 3 shows the dynamic coefficients and phase angles used in this model. Figure 4 presents a dynamic loading function for one person walking with step frequency equal to 2 Hz.

Finite Element Modelling of the Dynamic Behaviour of Tubular Footbridges 339

In Equation (2), the magnitudes ΔP1, ΔP2 and ΔP3 are associated with harmonic amplitudes. The first harmonic amplitude, ΔP1, is equal to 0.4P for fs equal to 2.0 Hz and 0.5P for fs equal to 2.4Hz. A simple interpolation between these two values can be used in intermediate cases. The second and third harmonic amplitudes, ΔP2 e ΔP3, were assumed to

The phase angles φ2 and φ3 depend on various other factors and should represent the most favourable used load combinations. In the present study the phase angles φ2 and φ3 were assumed to be equal to π/2 and phase angle φ1 was assumed to be equal to zero [15]. Figure 5 presents a dynamic loading function for one person walking with step frequency equal to

The fourth walking load model considered the same idea of the previous models. The main difference was the incorporation of the human heel effect in this particular load representation with the aid of Equations (3) to (6). The mathematical model behind this strategy was proposed by Varela [8] as well as a numerical approach to evaluate the floor

The proposed mathematical model, see Equations (3) to (6), used to represent the dynamic actions produced by people walking on floor slabs is not simply a Fourier series. This is due to the fact that the mentioned equations also incorporate the heel impact effect [8]. This loading model also considers a space and temporal variation of the dynamic action over the

Additionally, Load Model IV (LM-IV) also incorporates the transient effect due to the human heel impact [8]. The present investigation used a heel impact factor equal to 1.12 (fmi

**Figure 5.** LM-III: dynamic load function for one person walking (fs = 2.0 Hz).

i : harmonic phase angle;

**2.4. Load model IV (LM-IV)** 

structure reaction, as presented in Figure 6.

structure and is evaluated considering four harmonics.

be equal to 0.1P for fs equal to 2.0Hz [15].

t : time.

2 Hz.


**Table 3.** Dynamic coefficients and phase angles [9].

**Figure 4.** LM-II: dynamic load function for one person walking (fs = 2.0 Hz).

## **2.3. Load model III (LM-III)**

In this case a general expression is used to represent the excitation produced by an individual walking throughout time. These loads are produced with both feet, as function of a static part associated to the individual weight and three harmonics were considered to represent the dynamic action related to human walking [15], as illustrated in Equation (2). This dynamic loading model considers a space and temporal variation of the dynamic action over the footbridge.

$$\mathbf{F(t) = P + \Lambda P\_1 \sin \left(2\pi \mathbf{f\_s t - \phi p\right) + \Lambda P\_2 \sin \left(4\pi \mathbf{f\_s t - \phi p\right) + \Lambda P\_3 \sin \left(6\pi \mathbf{f\_s t - \phi p\right)}\right)}\tag{2}$$

Where:


i : harmonic phase angle; t : time.

338 Finite Element Analysis – New Trends and Developments

**Table 3.** Dynamic coefficients and phase angles [9].

**2.3. Load model III (LM-III)** 

over the footbridge.

F(t) : dynamic load;

P : person's weight (800 [1-4]); fs : walking step frequency;

Where:

In this load model, the time-dependent repeated force also can be represented by the Fourier series, as expressed in Equation (1) and four harmonics were considered to represent the dynamic action associated to human walking [9]. This model also considers a space and temporal variation of the dynamic action over the structural system. Table 3 shows the dynamic coefficients and phase angles used in this model. Figure 4 presents a dynamic

Harmonic i Dynamic Coefficients α<sup>i</sup> Phase Angles φ<sup>i</sup> 1 0.50 0 2 0.20 π/2 3 0.10 π 4 0.05 3π/2

loading function for one person walking with step frequency equal to 2 Hz.

**Figure 4.** LM-II: dynamic load function for one person walking (fs = 2.0 Hz).

In this case a general expression is used to represent the excitation produced by an individual walking throughout time. These loads are produced with both feet, as function of a static part associated to the individual weight and three harmonics were considered to represent the dynamic action related to human walking [15], as illustrated in Equation (2). This dynamic loading model considers a space and temporal variation of the dynamic action

F(t) = P + P1 sin (2 fs t - φ1) + P2 sin (4 fs t - φ2) + P3 sin (6 fs t - 3) (2)

**2.2. Load model II (LM-II)** 

In Equation (2), the magnitudes ΔP1, ΔP2 and ΔP3 are associated with harmonic amplitudes. The first harmonic amplitude, ΔP1, is equal to 0.4P for fs equal to 2.0 Hz and 0.5P for fs equal to 2.4Hz. A simple interpolation between these two values can be used in intermediate cases. The second and third harmonic amplitudes, ΔP2 e ΔP3, were assumed to be equal to 0.1P for fs equal to 2.0Hz [15].

The phase angles φ2 and φ3 depend on various other factors and should represent the most favourable used load combinations. In the present study the phase angles φ2 and φ3 were assumed to be equal to π/2 and phase angle φ1 was assumed to be equal to zero [15]. Figure 5 presents a dynamic loading function for one person walking with step frequency equal to 2 Hz.

**Figure 5.** LM-III: dynamic load function for one person walking (fs = 2.0 Hz).

## **2.4. Load model IV (LM-IV)**

The fourth walking load model considered the same idea of the previous models. The main difference was the incorporation of the human heel effect in this particular load representation with the aid of Equations (3) to (6). The mathematical model behind this strategy was proposed by Varela [8] as well as a numerical approach to evaluate the floor structure reaction, as presented in Figure 6.

The proposed mathematical model, see Equations (3) to (6), used to represent the dynamic actions produced by people walking on floor slabs is not simply a Fourier series. This is due to the fact that the mentioned equations also incorporate the heel impact effect [8]. This loading model also considers a space and temporal variation of the dynamic action over the structure and is evaluated considering four harmonics.

Additionally, Load Model IV (LM-IV) also incorporates the transient effect due to the human heel impact [8]. The present investigation used a heel impact factor equal to 1.12 (fmi

= 1.12). However, it must be emphasized that this value can vary substantially from personto-person. Figure 6 illustrates the dynamical load function for an individual walking at 2 Hz, based on Equations (3) to (6) and Tables 1 and 3 [9,12].

Finite Element Modelling of the Dynamic Behaviour of Tubular Footbridges 341

Fm : maximum Fourier series value, given by Equation (4);

The developed computational model adopted the usual mesh refinement techniques present in finite element method simulations, based on the ANSYS program [13]. The finite element model has been developed and validated with the experimental results [1-2]. This numerical model enabled a complete dynamic evaluation of the investigated tubular footbridge especially in terms of human comfort and its associated vibration serviceability limit states, see Figure 9. In this model, all steel tubular sections were represented by three-dimensional beam elements (PIPE16 and BEAM44) with tension, compression, torsion and bending capabilities. These elements have six degrees of freedom at each node: translations in the

nodal x, y, and z directions and rotations about x, y, and z axes, see Figure 9.

**Figure 6.** LM-IV: dynamic load function for one person walking (fp = 2.0 Hz).

On the other hand, the reinforced concrete slab was represented by shell finite elements (SHELL63), as presented in Figure 9. This finite element has both bending and membrane capabilities. Both in-plane and normal loads are permitted. The element has six degrees of freedom at each node: translations in the nodal x, y, and z directions and rotations about the

αi : dynamic coefficient for the harmonic force;

C1, C2 : coefficients given by Equations (5) and (6).

i : harmonic multiple (i = 1,2,3…,n);

**4. Finite element model** 

nodal x, y, and z axes.

fmi : human heel impact factor;

Tp : step period; fs : step frequency; i : harmonic phase angle; P : person's weight;

t : time;

### **3. Investigated structural model**

The structural model consists of tubular steel sections and a 100 mm concrete slab and is currently submitted to human walking loads [1-2]. The structure was based on a tubular composite (steel-concrete) footbridge, spanning 82.5 m. The structure is com-posed by three spans (32.5 m, 17.5 m and 20.0 m, respectively) and two overhangs (7.50 m and 5.0 m, respectively), as illustrated in Figures 7 and 8.

The steel sections used were welded wide flanges (WWF) made with a 300 MPa yield stress steel grade. A Young's modulus equal to 2.05 x 105 MPa was adopted for the tubular footbridge steel beams and columns. The concrete slab has a 20 MPa specified compression strength and a 2.13 x 104 MPa Young's Modulus.

$$F(t) = \begin{cases} \left(\frac{F\_{jm}F\_m - P}{0.04T\_p}\right)t + P & \text{if } 0 \le t < 0.04T\_p \\\\ \left.F\_{ml}\right|\_{m} \left[\frac{C\_1\left(t - 0.04T\_p\right)}{0.02T\_p} + 1\right] & \text{if } 0.04T\_p \le t < 0.06T\_p \\\\ F\_m & \text{if } 0.06T\_p \le t < 0.15T\_p \\\\ P + \sum\_{i=1}^{\text{ph}} P \exp\left[2\pi i f\_p \left(t + 0.1T\_p\right) + \phi\_i\right] & \text{if } 0.15T\_p \le t < 0.90T\_p \\\\ \left.10\left(P - C\_2\right)\left(\frac{t}{T\_p} - 1\right) + P & \text{if } 0.90T\_p \le t < T\_p \\\\ F\_m = P\left(1 + \sum\_{i=1}^{\text{ph}} \alpha\_i\right) \\\\ C\_1 = \left(\frac{1}{f\_{ml}} - 1\right) & \text{(5)} \end{cases} \tag{4}$$

$$C\_2 = \begin{cases} P\left(1 - \alpha\_2\right) & \text{if } nh = 3\\ P\left(1 - \alpha\_2 + \alpha\_4\right) & \text{if } nh = 4 \end{cases} \tag{6}$$

Where:


#### **4. Finite element model**

340 Finite Element Analysis – New Trends and Developments

**3. Investigated structural model** 

respectively), as illustrated in Figures 7 and 8.

strength and a 2.13 x 104 MPa Young's Modulus.

0.04

*mi m*

 

1

*p*

2

*i*

1

*nh*

m

*F t*

Where:

0.04

*p*

*p*

0.02

*Ct T*

10 1 0.90

*C*

 

 

2

*<sup>P</sup> if nh <sup>C</sup>*

*fF P t P if t T <sup>T</sup>*

*p*

2

*f F if T t T <sup>T</sup>*

*mi m p p*

( ) 0.06 0.15

*P sen if t T if T t T <sup>P</sup>*

1 *nh*

*m i i*

*<sup>t</sup> P C <sup>P</sup> if T t T <sup>T</sup>*

 

*F if T t T*

2 0.1 0.15 0.90

 1

 <sup>1</sup> <sup>1</sup> <sup>1</sup> *mi*

 

1 3 1 4

2 4

*p pi p p*

0 0.04

*p*

*p p*

*F P* (4)

*<sup>P</sup> if nh* (6)

*f* (5)

*p p*

(3)

1 0.04 0.06

based on Equations (3) to (6) and Tables 1 and 3 [9,12].

= 1.12). However, it must be emphasized that this value can vary substantially from personto-person. Figure 6 illustrates the dynamical load function for an individual walking at 2 Hz,

The structural model consists of tubular steel sections and a 100 mm concrete slab and is currently submitted to human walking loads [1-2]. The structure was based on a tubular composite (steel-concrete) footbridge, spanning 82.5 m. The structure is com-posed by three spans (32.5 m, 17.5 m and 20.0 m, respectively) and two overhangs (7.50 m and 5.0 m,

The steel sections used were welded wide flanges (WWF) made with a 300 MPa yield stress steel grade. A Young's modulus equal to 2.05 x 105 MPa was adopted for the tubular footbridge steel beams and columns. The concrete slab has a 20 MPa specified compression

The developed computational model adopted the usual mesh refinement techniques present in finite element method simulations, based on the ANSYS program [13]. The finite element model has been developed and validated with the experimental results [1-2]. This numerical model enabled a complete dynamic evaluation of the investigated tubular footbridge especially in terms of human comfort and its associated vibration serviceability limit states, see Figure 9. In this model, all steel tubular sections were represented by three-dimensional beam elements (PIPE16 and BEAM44) with tension, compression, torsion and bending capabilities. These elements have six degrees of freedom at each node: translations in the nodal x, y, and z directions and rotations about x, y, and z axes, see Figure 9.

**Figure 6.** LM-IV: dynamic load function for one person walking (fp = 2.0 Hz).

On the other hand, the reinforced concrete slab was represented by shell finite elements (SHELL63), as presented in Figure 9. This finite element has both bending and membrane capabilities. Both in-plane and normal loads are permitted. The element has six degrees of freedom at each node: translations in the nodal x, y, and z directions and rotations about the nodal x, y, and z axes.

Finite Element Modelling of the Dynamic Behaviour of Tubular Footbridges 343

**Figure 9.** Tubular footbridge finite element model

**5.1. Natural frequencies and vibration modes** 

**Table 4.** Tubular footbridge natural frequencies.

In a second phase, the steel-concrete composite tubular footbridge natural frequencies vibration modes and peak accelerations were determined with the aid of the numerical

It can be clearly noticed that there is a very good agreement between the structural model natural frequency values calculated using finite element simulations [1] and the experimental results [2], see Table 4. Such fact validates the finite element model here presented, as well as the results and conclusions obtained throughout this investigation. The

Tubular Footbridge Natural Frequencies (Hz) f01 f02 f03 Finite Element Model (see Figure 9) 1.61 2.12 5.39 Experimental Results 1.56 2.34 5.08 Error (%) 3.20 9.40 6.10

When the tubular footbridge freely vibrates in a particular mode, it moves up and down with a certain configuration or mode shape. Each footbridge natural frequency has an

simulations [1], based on the finite element method using the ANSYS program [13].

vibration modes of the tubular footbridge are depicted in Figures 10 to 12.

**Figure 7.** Investigated steel-concrete composite tubular footbridge.

**Figure 8.** Internal section of the investigated structural model.

The footbridge pier bearings were represented by a non-linear rotational spring element (COMBIN39). This element is a unidirectional element with non-linear generalized forcedeflection capability that can be used in any analysis.

The finite element model presented 71540 degrees of freedom, 11938 nodes and 15280 finite elements (BEAM44: 1056; PIPE16: 5642; SHELL63: 8580 and COMBIN39: 8), as presented in Figure 9. It was considered that both structural elements (steel tubular sections and concrete slab) have total interaction with an elastic behaviour.

## **5. Dynamic analysis**

Initially, the steel-concrete composite tubular footbridge natural frequencies, vibration modes and peak accelerations were determined based on experimental tests [2]. The peak acceleration values were obtained considering three types of human walking: slow walking, regular walking and fast walking.

**Figure 9.** Tubular footbridge finite element model

**Figure 7.** Investigated steel-concrete composite tubular footbridge.

**Figure 8.** Internal section of the investigated structural model.

deflection capability that can be used in any analysis.

slab) have total interaction with an elastic behaviour.

**5. Dynamic analysis** 

regular walking and fast walking.

The footbridge pier bearings were represented by a non-linear rotational spring element (COMBIN39). This element is a unidirectional element with non-linear generalized force-

The finite element model presented 71540 degrees of freedom, 11938 nodes and 15280 finite elements (BEAM44: 1056; PIPE16: 5642; SHELL63: 8580 and COMBIN39: 8), as presented in Figure 9. It was considered that both structural elements (steel tubular sections and concrete

Initially, the steel-concrete composite tubular footbridge natural frequencies, vibration modes and peak accelerations were determined based on experimental tests [2]. The peak acceleration values were obtained considering three types of human walking: slow walking, In a second phase, the steel-concrete composite tubular footbridge natural frequencies vibration modes and peak accelerations were determined with the aid of the numerical simulations [1], based on the finite element method using the ANSYS program [13].

## **5.1. Natural frequencies and vibration modes**

It can be clearly noticed that there is a very good agreement between the structural model natural frequency values calculated using finite element simulations [1] and the experimental results [2], see Table 4. Such fact validates the finite element model here presented, as well as the results and conclusions obtained throughout this investigation. The vibration modes of the tubular footbridge are depicted in Figures 10 to 12.


**Table 4.** Tubular footbridge natural frequencies.

When the tubular footbridge freely vibrates in a particular mode, it moves up and down with a certain configuration or mode shape. Each footbridge natural frequency has an

associated mode shape. It was verified that longitudinal amplitudes were predominant in the fundamental vibration mode (f01 = 1.61 Hz), related with Z axis direction, see Figure 10. In the second mode shape lateral displacements were predominant (f02 = 2.12 Hz), associated with X axis direction, as presented in Figure 11. On the other hand, in the third vibration mode (f03 = 5.39 Hz), the flexural effects were predominant, related to vertical amplitudes in the Y axis direction, as illustrated in Figure 12.

Finite Element Modelling of the Dynamic Behaviour of Tubular Footbridges 345

**Figure 12.** Vibration mode associated with the 3rd footbridge natural frequency (f03=5.39Hz).

The finite element modelling follows with the evaluation of the footbridge performance in terms of vibration serviceability due to dynamic forces induced by people walking. The first step of this investigation concerned in the determination of the tubular footbridge peak

The dynamic loading models (see Equations (1) to (6) and Figures 3 to 6), related to one, two and three people crossing the tubular footbridge on the concrete slab centre, in normal

The maximum accelerations (peak accelerations) were obtained utilizing an integration time step of 2x10-3 s (Δt = 2x10-3 s). In this investigation, seven sections of the structural model were analysed, see Figure 16. These maximum accelerations were compared to the limits recommended by design codes [9,14]. The structural damping coefficient adopted in this investigation was equal to 0.01 (ζ=1%), in accordance with the measured experimental

walking, see Figures 13 to 15, were applied on the investigated footbridge over 55.0 s.

**5.2. Determination of the tubular footbridge peak accelerations** 

accelerations, based on a linear time-domain dynamic analysis.

**Figure 13.** One person walking on the footbridge (regular walking).

**Figure 14.** Two people walking on the footbridge (regular walking).

damping [2].

**Figure 10.** Vibration mode associated with the 1st footbridge natural frequency (f01=1.61 Hz).

**Figure 11.** Vibration mode associated with the 2nd footbridge natural frequency (f02=2.12 Hz).

**Figure 12.** Vibration mode associated with the 3rd footbridge natural frequency (f03=5.39Hz).

## **5.2. Determination of the tubular footbridge peak accelerations**

344 Finite Element Analysis – New Trends and Developments

the Y axis direction, as illustrated in Figure 12.

associated mode shape. It was verified that longitudinal amplitudes were predominant in the fundamental vibration mode (f01 = 1.61 Hz), related with Z axis direction, see Figure 10. In the second mode shape lateral displacements were predominant (f02 = 2.12 Hz), associated with X axis direction, as presented in Figure 11. On the other hand, in the third vibration mode (f03 = 5.39 Hz), the flexural effects were predominant, related to vertical amplitudes in

**Figure 10.** Vibration mode associated with the 1st footbridge natural frequency (f01=1.61 Hz).

**Figure 11.** Vibration mode associated with the 2nd footbridge natural frequency (f02=2.12 Hz).

The finite element modelling follows with the evaluation of the footbridge performance in terms of vibration serviceability due to dynamic forces induced by people walking. The first step of this investigation concerned in the determination of the tubular footbridge peak accelerations, based on a linear time-domain dynamic analysis.

The dynamic loading models (see Equations (1) to (6) and Figures 3 to 6), related to one, two and three people crossing the tubular footbridge on the concrete slab centre, in normal walking, see Figures 13 to 15, were applied on the investigated footbridge over 55.0 s.

The maximum accelerations (peak accelerations) were obtained utilizing an integration time step of 2x10-3 s (Δt = 2x10-3 s). In this investigation, seven sections of the structural model were analysed, see Figure 16. These maximum accelerations were compared to the limits recommended by design codes [9,14]. The structural damping coefficient adopted in this investigation was equal to 0.01 (ζ=1%), in accordance with the measured experimental damping [2].

**Figure 13.** One person walking on the footbridge (regular walking).

**Figure 14.** Two people walking on the footbridge (regular walking).


Finite Element Modelling of the Dynamic Behaviour of Tubular Footbridges 347

**Figure 18.** LM-II: tubular footbridge acceleration response at section B. One pedestrian crossing the

**Figure 19.** LM-III: tubular footbridge acceleration response at section B. One pedestrian crossing the

concrete slab centre at resonance condition. Normal walking.

concrete slab centre at resonance condition. Normal walking.

**Figure 15.** Three people walking on the footbridge (regular walking).

**Figure 16.** Tubular composite (steel-concrete) footbridge investigated sections.

In sequence, Figures 17 to 20 illustrate the tubular footbridge dynamic response, along the time, related to the section B (see Figure 16), when one pedestrian crosses the footbridge in regular walking (resonance condition).

**Figure 17.** LM-I: tubular footbridge acceleration response at section B. One pedestrian crossing the concrete slab centre at resonance condition. Normal walking.

regular walking (resonance condition).

**Figure 15.** Three people walking on the footbridge (regular walking).

**Figure 16.** Tubular composite (steel-concrete) footbridge investigated sections.

In sequence, Figures 17 to 20 illustrate the tubular footbridge dynamic response, along the time, related to the section B (see Figure 16), when one pedestrian crosses the footbridge in

**Figure 17.** LM-I: tubular footbridge acceleration response at section B. One pedestrian crossing the

concrete slab centre at resonance condition. Normal walking.

**Figure 18.** LM-II: tubular footbridge acceleration response at section B. One pedestrian crossing the concrete slab centre at resonance condition. Normal walking.

**Figure 19.** LM-III: tubular footbridge acceleration response at section B. One pedestrian crossing the concrete slab centre at resonance condition. Normal walking.

Finite Element Modelling of the Dynamic Behaviour of Tubular Footbridges 349

(alim in m/s2 Investigated Sections )\*

(alim in m/s2 Investigated Sections )\*

(alim in m/s2 Investigated Sections )\*

Limit Accelerations

0.49

Limit Accelerations

0.49

Limit Accelerations

0.49

On the other hand, these maximum acceleration values increases when the normal walking associated to two and three people (see Figures 14 and 15) is considered in the analysis, as

> Tubular Footbridge Peak Accelerations (ap in m/s2)

A B1 B B2 C D E

Tubular Footbridge Peak Accelerations (ap in m/s2)

A B1 B B2 C D E

Tubular Footbridge Peak Accelerations (ap in m/s2)

A B1 B B2 C D E

Load Model I (LM-I) 1.38 0.17 0.53 0.17 0.53 0.39 0.72 Load Model I (LM-II) 1.50 0.18 0.58 0.18 0.58 0.43 0.79 Load Model I (LM-III) 1.38 0.17 0.55 0.17 0.55 0.40 0.74 Load Model I (LM-IV) 1.00 0.16 0.44 0.16 0.38 0.33 0.78

**Table 5.** Structural model peak accelerations corresponding to one individual walking.

Load Model I (LM-I) 2.14 0.19 0.74 0.19 0.72 0.61 1.23 Load Model I (LM-II) 2.32 0.21 0.81 0.21 0.80 0.67 1.35 Load Model I (LM-III) 2.14 0.20 0.76 0.19 0.75 0.63 1.27 Load Model I (LM-IV) 1.87 0.21 0.57 0.22 0.58 0.55 1.36

**Table 6.** Structural model peak accelerations corresponding to two people walking.

Load Model I (LM-I) 2.71 0.25 0.97 0.25 0.94 0.81 1.68 Load Model I (LM-II) 2.93 0.27 1.06 0.27 1.04 0.89 1.84 Load Model I (LM-III) 2.70 0.26 1.00 0.25 0.97 0.84 1.73 Load Model I (LM-IV) 2.61 0.31 0.76 0.32 0.75 0.73 1.86

**Table 7.** Structural model peak accelerations corresponding to three people walking.

presented in Tables 5 to 7.

Dynamic Loading Models

Dynamic Loading Models

Dynamic Loading Models

\*alim = 1.5%g = 0.15 m/s2: indoor footbridges [9,14] \*alim = 5.0%g = 0.49 m/s2: outdoor footbridges [9,14]

\*alim = 1.5%g = 0.15 m/s2: indoor footbridges [9,14] \*alim = 5.0%g = 0.49 m/s2: outdoor footbridges [9,14]

\*alim = 1.5%g = 0.15 m/s2: indoor footbridges [9,14] \*alim = 5.0%g = 0.49 m/s2: outdoor footbridges [9,14]

**Figure 20.** LM-IV: tubular footbridge acceleration response at section B. One pedestrian crossing the concrete slab centre at resonance condition. Normal walking.

Figures 17 to 20 present the vertical acceleration versus time graph for the tubular footbridge at section B (see Figure 16). These figures show that the vertical acceleration of the structure gradually increase along the time. In this particular case, the third harmonic with a 2.0 Hz step frequency (fs = 2.0 Hz), was the walking load resonant harmonic.

The maximum acceleration value found at section B (see Figure 16) was equal to 0.53 m/s2 (LM-I), 0.58 m/s2 (LM-II), 0.55 m/s2 (LM-III) and 0.44 m/s2 (LM-IV), as illustrated in Figures 17 to 20. These figures also indicate that from the moment that the pedestrian leaves the footbridge span (Section B, see Figure 16), when the time is approximately equal to 26 s, the structural damping minimises the dynamic structural model response, as presented in Figures 17 to 20. This assertive occurs only in dynamic loading models that consider the load spatial variation.

The peak acceleration analysis was focused on the steel-concrete composite tubular footbridge dynamic behaviour when the pedestrian normal walking was considered in this work. In sequence, Tables 5 to 7 present the maximum accelerations (peak accelerations: ap in m/s2), related to seven structural sections of the investigated footbridge (A, B, B1, B2, C, D and E), as illustrated in Figure 16.

The maximum acceleration values (peak accelerations) found in this investigation are respectively equal to 1.50 m/s2 (Section A), 0.18 m/s2 (Section B1), 0.58 m/s2 (Section B), 0.18 m/s2 (Section B2), 0.58 m/s2 (Section C), 0.43 m/s2 (Section D) and 0.79 m/s2 (Section C), corresponding to one individual crossing the composite footbridge in normal walking (resonance condition), as illustrated in Figure 13.

On the other hand, these maximum acceleration values increases when the normal walking associated to two and three people (see Figures 14 and 15) is considered in the analysis, as presented in Tables 5 to 7.


\*alim = 1.5%g = 0.15 m/s2: indoor footbridges [9,14]

348 Finite Element Analysis – New Trends and Developments

**Figure 20.** LM-IV: tubular footbridge acceleration response at section B. One pedestrian crossing the

with a 2.0 Hz step frequency (fs = 2.0 Hz), was the walking load resonant harmonic.

Figures 17 to 20 present the vertical acceleration versus time graph for the tubular footbridge at section B (see Figure 16). These figures show that the vertical acceleration of the structure gradually increase along the time. In this particular case, the third harmonic

The maximum acceleration value found at section B (see Figure 16) was equal to 0.53 m/s2 (LM-I), 0.58 m/s2 (LM-II), 0.55 m/s2 (LM-III) and 0.44 m/s2 (LM-IV), as illustrated in Figures 17 to 20. These figures also indicate that from the moment that the pedestrian leaves the footbridge span (Section B, see Figure 16), when the time is approximately equal to 26 s, the structural damping minimises the dynamic structural model response, as presented in Figures 17 to 20. This assertive occurs only in dynamic loading models that consider the

The peak acceleration analysis was focused on the steel-concrete composite tubular footbridge dynamic behaviour when the pedestrian normal walking was considered in this work. In sequence, Tables 5 to 7 present the maximum accelerations (peak accelerations: ap in m/s2), related to seven structural sections of the investigated footbridge (A, B, B1, B2, C, D

The maximum acceleration values (peak accelerations) found in this investigation are respectively equal to 1.50 m/s2 (Section A), 0.18 m/s2 (Section B1), 0.58 m/s2 (Section B), 0.18 m/s2 (Section B2), 0.58 m/s2 (Section C), 0.43 m/s2 (Section D) and 0.79 m/s2 (Section C), corresponding to one individual crossing the composite footbridge in normal walking

concrete slab centre at resonance condition. Normal walking.

load spatial variation.

and E), as illustrated in Figure 16.

(resonance condition), as illustrated in Figure 13.

\*alim = 5.0%g = 0.49 m/s2: outdoor footbridges [9,14]

**Table 5.** Structural model peak accelerations corresponding to one individual walking.


\*alim = 1.5%g = 0.15 m/s2: indoor footbridges [9,14]

\*alim = 5.0%g = 0.49 m/s2: outdoor footbridges [9,14]

**Table 6.** Structural model peak accelerations corresponding to two people walking.


\*alim = 1.5%g = 0.15 m/s2: indoor footbridges [9,14]

\*alim = 5.0%g = 0.49 m/s2: outdoor footbridges [9,14]

**Table 7.** Structural model peak accelerations corresponding to three people walking.

It must be emphasized that the footbridge overhang sections (Sections A and E, see Figure 16) have presented very high peak accelerations, in all investigated situations, as presented in Tables 5 to 7. This is explained due to the fact that the transient impact produced by the pedestrians on the entrance and exit of the investigated structure have generated high acceleration values.

Finite Element Modelling of the Dynamic Behaviour of Tubular Footbridges 351

The results found throughout this investigation have indicated that the dynamic actions produced by human walking could generate peak accelerations that surpass design criteria limits developed for ensuring human comfort. Hence it was detected that this type of structure can reach high vibration levels that can compromise the footbridge user's comfort

The analysis methodology presented in this paper is completely general and is the author's intention to use this solution strategy on other pedestrian foot-bridge types and to investigate the fatigue problem. The fatigue problem is a relevant issue and certainly much more complicated and is influenced by several design parameters and footbridge types.

José Guilherme Santos da Silva, Ana Cristina Castro Fontenla Sieira, Gilvan Lunz Debona, Pedro Colmar Gonçalves da Silva Vellasco and Luciano Rodrigues Ornelas de Lima

The authors gratefully acknowledge the support for this work provided by the Brazilian

[1] Debona GL. Modelagem do comportamento dinâmico de passarelas tubulares em aço e mistas (aço-concreto) (Modelling of the dynamic behaviour of steel-concrete composite tubular footbridges), MSc Dissertation (in Portuguese), Civil Engineering Post-Graduate Programme, PGECIV, State University of Rio de Janeiro, UERJ, Rio de

[2] Zúñiga JEV. Análise da resposta dinâmica experimental de uma passarela tubular mista, aço-concreto, submetida ao caminhar humano (Dynamic experimental analysis of a steel-concrete composite tubular footbridge submitted to human walking), MSc Dissertation (in Portuguese), Civil Engineering Post-Graduate Programme, PGECIV,

[4] Silva JGS da, Vellasco PCG da S, Andrade SAL de, Lima LRO de, Figueiredo FP. Vibration analysis of footbridges due to vertical human loads. Computers & Structures,

[5] Eurocode 3: Design of steel structures. Part 1.9: General rule - Fatigue. European

State University of Rio de Janeiro, UERJ, Rio de Janeiro, Brazil, pp. 1-135; 2011. [3] Figueiredo FP, Silva JGS da, Vellasco PCG da S, Andrade SAL de, Andrade, SAL de. A parametric study of composite footbridges under pedestrian walking loads.

Further research in this area is currently being carried out.

*State University of Rio de Janeiro (UERJ), Rio de Janeiro/RJ, Brazil* 

Science Foundation CAPES, CNPq and FAPERJ.

Janeiro, Brazil, pp. 1-154; 2011.

2007; 85:1693-1703.

Engineering Structures, 2008; 30:605-615.

Committee for Standardisation; 2005.

and especially its safety.

**Author details** 

**Acknowledgement** 

**7. References** 

Based on a quantitative analysis of the maximum accelerations, it was verified that the loading model I (LM-I) has produced the highest peck acceleration values practically in all investigated cases, as illustrated in Tables 5 to 7.

These peak accelerations presented in Tables 5 to 7 are related to a pedestrian normal walking situation. It must be emphasized that the limit acceleration value is equal to 0.49 m/s2 [9,14], when the outdoor footbridges are considered in the analysis.

Based on the finite element modelling of the steel-concrete composite tubular footbridge dynamic behaviour, the numerical results presented in Tables 5 to 7 indicated that the dynamic actions produced by human walking led to peak accelerations higher than the limiting values present in design code recommendations (Outdoor footbridges: alim = 5%g = 0.49 m/s² [9,14]), as depicted in Tables 5 to 7.

## **6. Final remarks**

This contribution covers the application of tubular structural elements in pedestrian footbridge design and tries to give an overview about the evaluation of tubular footbridges dynamic behaviour, objectifying to help practical structural engineers to deal with this kind of problem and to allow for a further application of tubular structural elements in pedestrian footbridge design.

The present investigation was carried out based on four dynamic loading models (LM-I to LM-IV) implemented objectifying to incorporate the dynamic effects induced by people walking on the footbridges dynamic response. In these models, the position of the human walking load was changed according to the individual position. However, a more realistic loading model (LM-IV) considered the ascent and descending movement of the human body effective mass at each step load (human walking load) and additionally also incorporates the transient effect due to the human heel impact.

The proposed analysis methodology considered the investigation of the dynamic behaviour, in terms of serviceability limit states, of a composite tubular footbridge, spanning 82.5 m. The structure is composed by three spans (32.5 m, 17.5 m and 20.0 m, respectively) and two overhangs (7.50 m and 5.0 m, respectively). The structural system is constituted by tubular steel sections and a concrete slab and is currently used for pedestrian crossing.

A computational model, based on the finite element method, was developed using the ANSYS program. This model enabled a complete dynamic evaluation of the investigated tubular footbridge especially in terms of human comfort and its associated vibration serviceability limit states.

The results found throughout this investigation have indicated that the dynamic actions produced by human walking could generate peak accelerations that surpass design criteria limits developed for ensuring human comfort. Hence it was detected that this type of structure can reach high vibration levels that can compromise the footbridge user's comfort and especially its safety.

The analysis methodology presented in this paper is completely general and is the author's intention to use this solution strategy on other pedestrian foot-bridge types and to investigate the fatigue problem. The fatigue problem is a relevant issue and certainly much more complicated and is influenced by several design parameters and footbridge types. Further research in this area is currently being carried out.

## **Author details**

350 Finite Element Analysis – New Trends and Developments

investigated cases, as illustrated in Tables 5 to 7.

0.49 m/s² [9,14]), as depicted in Tables 5 to 7.

acceleration values.

**6. Final remarks** 

pedestrian footbridge design.

serviceability limit states.

It must be emphasized that the footbridge overhang sections (Sections A and E, see Figure 16) have presented very high peak accelerations, in all investigated situations, as presented in Tables 5 to 7. This is explained due to the fact that the transient impact produced by the pedestrians on the entrance and exit of the investigated structure have generated high

Based on a quantitative analysis of the maximum accelerations, it was verified that the loading model I (LM-I) has produced the highest peck acceleration values practically in all

These peak accelerations presented in Tables 5 to 7 are related to a pedestrian normal walking situation. It must be emphasized that the limit acceleration value is equal to 0.49

Based on the finite element modelling of the steel-concrete composite tubular footbridge dynamic behaviour, the numerical results presented in Tables 5 to 7 indicated that the dynamic actions produced by human walking led to peak accelerations higher than the limiting values present in design code recommendations (Outdoor footbridges: alim = 5%g =

This contribution covers the application of tubular structural elements in pedestrian footbridge design and tries to give an overview about the evaluation of tubular footbridges dynamic behaviour, objectifying to help practical structural engineers to deal with this kind of problem and to allow for a further application of tubular structural elements in

The present investigation was carried out based on four dynamic loading models (LM-I to LM-IV) implemented objectifying to incorporate the dynamic effects induced by people walking on the footbridges dynamic response. In these models, the position of the human walking load was changed according to the individual position. However, a more realistic loading model (LM-IV) considered the ascent and descending movement of the human body effective mass at each step load (human walking load) and additionally also

The proposed analysis methodology considered the investigation of the dynamic behaviour, in terms of serviceability limit states, of a composite tubular footbridge, spanning 82.5 m. The structure is composed by three spans (32.5 m, 17.5 m and 20.0 m, respectively) and two overhangs (7.50 m and 5.0 m, respectively). The structural system is constituted by tubular

A computational model, based on the finite element method, was developed using the ANSYS program. This model enabled a complete dynamic evaluation of the investigated tubular footbridge especially in terms of human comfort and its associated vibration

steel sections and a concrete slab and is currently used for pedestrian crossing.

m/s2 [9,14], when the outdoor footbridges are considered in the analysis.

incorporates the transient effect due to the human heel impact.

José Guilherme Santos da Silva, Ana Cristina Castro Fontenla Sieira, Gilvan Lunz Debona, Pedro Colmar Gonçalves da Silva Vellasco and Luciano Rodrigues Ornelas de Lima *State University of Rio de Janeiro (UERJ), Rio de Janeiro/RJ, Brazil* 

## **Acknowledgement**

The authors gratefully acknowledge the support for this work provided by the Brazilian Science Foundation CAPES, CNPq and FAPERJ.

## **7. References**

	- [6] Kuhlmann U, H-P Günther, Saul R, Häderle, M.-U. 2003. Welded circular hollow section (CHS) joints in bridges. ISTS 2003: Proceedings of the 10th International Symposium on Tubular Structures, Madrid, Spain.

**Chapter 16** 

© 2012 Akour and Maaitah, licensee InTech. This is an open access chapter distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/3.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Attribution License (http://creativecommons.org/licenses/by/3.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

© 2012 Akour and Maaitah, licensee InTech. This is a paper distributed under the terms of the Creative Commons

**Finite Element Analysis of Loading** 

**Beyond the Yield Limit** 

Salih Akour and Hussein Maaitah

http://dx.doi.org/10.5772/50433

**1. Introduction** 

Additional information is available at the end of the chapter

to transfer the forces between sandwich panel components.

**1.1. Main principles of sandwich structures** 

**Area Effect on Sandwich Panel Behaviour** 

Research efforts continuously are looking for new, better and efficient construction materials. The main goal of these researches is to improve the structural efficiency, performance and durability. New materials typically bring new challenges to designer who utilizes these new materials. In the past decades various sandwich panels have been implemented in aerospace, marine, architectural and transportation industry. Light-weight, excellent corrosion characteristics and rapid installation capabilities created tremendous opportunities for these sandwich panels in industry. Sandwich panel normally consists of a low-density core material sandwiched between two high modulus face skins to produce a lightweight panel with exceptional stiffness as shown in Figure 1. Face skins act like flanges of an I-beam. These faces are typically bonded to a core to achieve the composite action and

Typical sandwich composite construction consists of three main components as illustrated in Figure 1. The sandwich consists of two thin, stiff and strong faces are separated by thick, light and weaker core. Faces and core materials are bonded together with an adhesive to facilitate the load transfer mechanism between the components, therefore effectively utilize all the materials used. The two faces are placed at a distance from each other to increase the moment

In sandwich structure, typically the core material is not rigid compared to face sheets; therefore, the shear deflection within the core is insignificant in most cases. The shear deflection in the faces can be also neglected. The effect of shear rigidity in the core is shown

of inertia, and consequently the flexural rigidity, about the neutral axis of the structure.

