Numerical Modelling of Tunnels: G. Swoboda University of Innsbruck, Innsbruck, Austria

NUMERICAL MODELLING OF TUNNELS
G. Swoboda
University of Innsbruck, Innsbruck, Austria
ABSTRACT
Until a few years ago tunnel construction was based exclusively on experience. Nu-
merical methods, however, constitute a very valuable supplement. Shown here are the
approximations necessary for 2D analysis. The most important load cases, such as
dead load or water pressure, are also illustrated. The damage tensor theory needed for
realistic simulation of jointed rock is also presented. In future 3D analysis will take on
increasing significance, for which reason the pertinent models are also dealt with.
C. S. Desai et al. (eds.), Numerical Methods and Constitutive Modelling in Geomechanics

© Springer-Verlag Wien 1990
278 G. Swoboda
1 Introduction
This chapter is intended as a guide to those engineers responsible for th(~ dP.sign of
tunnels and aims to give them an overview of Lhe state of the art of practical 1111 -
merical models for tunneling. For many years the design of tunrwls wa:; hasPd soi('I.Y
on experience and sometimes on very simple analytical considerations. This d(H'S not
mean that a tunnel cannot be designed without the use of a numeri<:<d rnodd. Th(' fa.d.
is that such computer models are very helpful tools in undcrsta.ndi ng tlw nwcha.n ica.l
background of the excavation process.
Every tunnel design begins by selecting from engineering dra.wings the most critical
cross section for construction of the tunnel. This can be a function of the cxca.va.t.cd
area, the distance to other tunnels, the height of overburden, the influence of buildings
on the surface and the geological conditions.
E = 20000 kN/m 2
j<;;;:----"t---.._::,¥--;---:1f<;::-""'k;;:--'!E-;---~~----1!~- - v = 0 . 35
ko " 0 . 5
c " 10 kN/m 2
~ " zo•
E : 60000 kN/m 2
v " 0.7 ko" 0 . 7
c " 30 kN/m 2
zo•¢> "
~-----*":::...._~---\-'~f-;~ .---"--- : ----------

DC~
'"- E = 40000 kN/m 2

i.. v" 0 . 4 k 0 = 0 . 7
8 c = 25 kN/m 2
i.
9
il> = 20
Figure 1: Numerical model for a subway tunnel
Fig. 1 shows a typical numerical model for a subway tunnel. Modem mesb gener-
ation allows such a numerical model to be modified in a wide range due to geometric
and geological modifications between different cross sections.
The question of safety is so important for subway tunnels that hardly any design
or construction work can be performed without numerical model. Numerical study of
stress redistribution was of great importance for the progress of modern cross sections
and new excavation procedures. Fig. 2 shows different complex subway cross sections.
Numerical Modelling of Tunnels 279
I. 6 10m
Figure 2: Complex cross section for subway tunnels
2 Two-dimensional models
Numerical models in tunneling started with a full-load two-dimensional model (Fig. 3)
that was based on a rigid supporting shield driven in undisturbed primary state of
stress. Thus protected, the primary support (shotcrete lining) was erected without
stress. Only then was the imaginary supporting shield removed, with static pressure
being completely redistributed to the shotcrete. This type of model is very conservative
and stays on the safe side. It was not able to be used for tunnels with great overburden.
For this reason, the 3D support effect, as simulated by the strain relief in front of
the face, was taken into account.
Three models arc widely used to simulate the t hree-dimensional effect. These are
the
- stiffness reduction model
- load reduction model
- modified stiffness reduction model
'fhc stiffness reduction model uses a support core with a modified modulus of elasticity.
Fo r excavation, the full-stress field in front of the tunnel face is applied to this core.
Modification of the original modulus of elasticity E to the modulus of the support core
is done' with
(1)
Wlwn the o value is infinite, a = oo, we find the original initial stress field in front
of the face after stress relief. The stiffness factor a is calibrated by displacement
llH'<ts urcment based on the relation:
(2)
280 G. Swoboda
Stage 0 Stage 1 Construction Construction

Stage A Stage B
Initial state of stress rigi d support ing shiel d
removed shield
Eifio
0•0
E•O
Figure 3: Full load model
which means that the relation of the displacement on the front of the face uA to the
final displacement UB must be the same in the field (UA, UB) as in the numerical mode
(uA, UB) .
In the load reduction method, mainly used in Germany, the entire load from the
initial stresses is applied. At the same time, however, a support load (JN is applied
and the effective load on (1 - (J)N is reduced. In order to prevent uplift of the system,
the load factor has to be assumed variably over the circumference. The load factor can
be ~:tltin 0 < (3 < 1.
The modified stiffness reduction method attempts to create a synthesis of these two
methods. Here the load factor (3 is applied supplementary to the stiffness reduction
factor, with the result that the factor ab is within 0 < ab < 1.
3 Loading conditions
The final forces in the primary support system (shotcrete lining) are based on the
following influences:
Construction Stage 0 Construction Stage A Construction Stage 8

Initial state of stress Initial stress relieve Final construction stage
Stiffness reduction Support load

method method
E do E II do E
@
,,l(.(
/ \ I /TX
I ,_
fEs·aE\ (ds·~cfol
". . .t . . . . 'UA \ ...... j,..,/
\ I
UA Ue
d•f(t) ~. f(t)
0-sd <ro O-s~-s1
St i ffness factor Load factor
dy•yh
Figure 4: Two-dimensional simulation of the three-dimensional effect
3.1 Stress redistribution due to excavation

Tunnel excavation is simulated in a numerical model using the following steps:
- The elements are eliminated in the excavation area. This can be achieved by
reducing the elasticity modulus to zero or by deleting the elements from the list
of elements.
- At the same time, the primary stresses a 0 present in the element have to be
eliminated by integration over the clement.
This integration is performed element by element according to [7] with
(3)
The nodal forces {!} thus determined for an element are equilibrium loads. The
total vector { F} is the sum of all the elements E involved in the excavation.
E
{F} = L:Ue}
e=l
Since only the nodes are located on the excavation boundary, this total vector has
to be divided into two vectors
282 G. Swoboda
Displacement in front of face Additional displacement

u. 1.0 u.r"•:..___,.......,-,.----:t;;....,
Uc:I•O Ur
Cross section A Cross section B
(c) IEs· <!E I

Support
stiffness
Measured displacement
Excavation
(d)
Ua
Figure 5: Relation between the displacement in front of the face and final displacement
{F} = {~}
(5)
{F;} = {0}
whereby {Fa} represents the loads on the excavation boundary and { F;} those in the
interior, which, however, are not to be considered. Only now is the equilibrium load
vector { F} transformed into a load vector {Fa}, whereby the sum corresponds to the
dead load of the excavated material.
The dead load of the excavated material acts in the global y direction and gives rise
to the following problem. The size of the displacement from this load is dependent on
the size of the chosen FEM mesh. However, this is a problem that also occurs when
attempting an analytical solution.
In order to examine the question of loads in a continuum, a surface load is considered
in place of the internal dead load because it makes the problem easier to illustrate. In
the two-dimensional halfspace (Fig. 7) the displacement of a certain point P( x, y) can
be determined according to Girkmann [8]
Numerical Modelling of Tunnels
excavated area ---,

--.... .
\
\
X
Figure 6: Simulation of excavation
X X
a • P(x,y)
y P(O,y,O)
Q y
Figure 7: Line load in two- and three-dimensional halfspace
v = (1 + fl) F. [(1- fl) ln x2 + y2 - y2 ] (6)

E 7r h a2 x2 + y2
whereby h represents the thickness and a the distance to a reference point Q, where the
vertical displacement is v = 0. If the distance is infinite, the displacement V is infinite.
This means that we have a problem identical to the FEM model, namely displacement
dependent on a reference point.
Since this problem only occurs in the loading "dead load of the excavated material"
and, due to the type of driving employed, the loading only affects a. short dista.nn~,
it should be examined whether similar problems occur in three-dimensional ha.lfspa.cc.
Since finite displacements occur as a result of concentrated loads after Boussincsq,
this question was already investigated by Girkmann [9]. The displa.rcnwnt. of a. line
284 G. Swoboda
load p, which corresponds to the concentrated load F in the change from the two- to
three-dimensional model reads
(l+J-L)P [ l b+/x2+y2+b2 2 bxz ]

v= 2E7r . 4(1-Jt) n JxZ+yz + (x2+y2)/x2+y2+b2 (7)
whereby b is the length of the line load.
b [m] log
25 50 75 100 100 1000 10000 b [m]
10
10
20
50
30
40
100
v [mm] v [mm]
Figure 8: Displacement of a. line load on a. three-dimensional halfspace
Fig. 8 depicts the displacements for line loads of varied lengths and for points at
various depths. H is especially from the logarithmic diagram of the line load that we
sec that the displacement is a function of the load's length. Therefore, it can be said
that the displacements of a11 infinitely long line load arc infinitely large, whereby only
an infinite line load can be compared with a two-dimensional example.
Figure !l: Nonlinearity in the origin of the stress-strain curve
'!'he reason why it is so difficult to obtain an unmistakable solution is that the

nonlinearity iu the origin (Fig. 9) is not taken into consideration in the constitutive
law in general use today. For points whose initial strain of the theoretical stress-strain
diagram is not exceeded, integration of the analysis has to be interrupted if the strains
are less than the initial strain of the material co. In order for the numerical model to
take account of the dead load, since this is a very short line load and thus the strains
very quickly drop below the initial strains, the mesh size has to be reduced accordingly.
The program is easily revised to meet this situation, in that the elasticity modulus is set
at infinite for all mesh elements near the boundary where the strain e0 is not reached.
As indicated in Fig. 10, this brings about a modification in the support conditions.
''
\, 0
\
\
\
\
\
' ......... _____ ,..

''
new support
boundary
supports
Figure 10: l'vlodification of the boundary condition
3.2 Water pressure

The calculation of water load presents a special problem in the construction of shallow
tunnels in loose material. With earth pressure it can be assumed that the earth pressure
is reduced by the spatial carrying effect. A similar assumption is in no way possible
for water load, however since it will always exert its full load on the concrete arch.
When calculating this load case it must be decided whether the soil is damp, zone
a in Fig. 11, or wa.ter-saturatcd, zones b a.nd c. The soil stresses in the first case are
O"y = 1h
(8)
O"x = A O'y
whereby 1 represents the specific weight of the damp soil and .X the coefficient of earth
pressure at rest.. The stresses O"x and O'y are identical to the effective stresses 0"9 in the
soil's grain. Thus, t.hc occurring expansion can be determined by removing the soil
stress from the following clastic relationship:
{s} = [D]{a} = [D] { (9)
For zone b, which is under water pressure, the soil stress is:
286 G. Swoboda
Soilpressure or
BEAM- Element
Pore water pressure
JOINT- Element
Gap if llaw,r <or
Figure 11: Computational model for the loading case water
iiy = b' + 1.0) h =

jj g = -y'h ow (10)
iix = A.Uy
whereby -y' is the specific weight under uplift. The soil stresses are divided here into
the effective soil stresses i79 and the pore-water pressure ow.
If the water level drops due to construction, the geomechanical changes in stresses have
two causes:
(i) Change in effective soil stresses
Removal of the pore-water pressure in zone b causes an increase in the effective
soil stress from G-9 to ii9 • For the related expansion it must be assumed that
a
the expansions related to 9 have already taken place. This means the newly
occurred expansions in zone b are
{c} = [D]{~a9 } = [D] { (11)

(ii) Pore-water pressure on a joint face

If the drop in the water level causes water to accumulate along a joint face (Fig. 12,
Stage II) for example sheet piling, the full pore-water pressure ow will become
effective here.
The rcdisposition of forces described in (i) and (ii) produce a new primary stress
state in the soil. This is the decisive stress for the earth pressure during ecxavation,
meaning that pressure on the concrete arch is O"T. If the water pressure becomes
effective after completion of the tunnel, the reversible processes described in (i) and
(ii) will take place. This means the effective soil stress will be reduced and the pore-
water pressure removed from the joint face.
(iii) Pore-water pressure on the concrete arch
As shown in Fig. 12, the pore-water pressure aw will act on the joint face of
the tunnel arch identical to load (ii). In this way the tangential earth pressure
changes by the amount of ~O"IV,T· If the tensile stress ~aw,T becomes larger than
the original earth pressure, tensile stresses would occur between the concrete arch
and the soil, which however, cannot be the case. This is solved by the structure
of the calculation model shown in Fig. 11. A nonlinear JOINT element is inserted
between the concrete arch (DEAM elements) and the soil (LST elements).
When the tensile failure condition
O"T - ~0"1\",T >0 -t ET = 0 (12)

is reached the tangential modulus of elasticity equals 0. In this way the gap
between the soil and the concrete arch can be simulated.
4 Jointed rock
Uuti! now numerical models have not been applied to jointed rock because it is much
more difficult to predict the orimtation of a joint system than any soil parameter.
llut iu tunnels with pilot tunnel excavation systems the joint system is well known
and can be considered in the numerical model.
Goodman's [15] works first presented an element that describes the rock jointing.
This dement can only be used to describe individual faults and is also numerically
very i ust.able. Improvements for this clement have been proposed by many authors. In
n•n·nt. years, models have been developed to break down the discontinuum into blocks.
In the rigid body model (Rll~f) by Cundall [19], the rock mass is broken down into rigid
blocks. This is a model that docs not take consideration of the elastic deformations
and is certainly not suitable for <·a.lcula.ting shallow tunnel structures. In the further
developments by Asai, the clastic properties are simulated by springs in the nodes. In
288 G. Swoboda
STAGE I STAGE II Lowered ground-water table
a,
"·~-~- a,Ow
a,
ag•Y' h
STAGE Ill Filling with water
Figure 12: Loading case water
places where the block's displacement share is very large, a usable formulation of a
model around the discontinuum definitely has to be described.
The approach taken here, namely that of the "decoupled finite element method",
uses finite elements as elastic blocks. According to the faults, these are coupled by
constraint elements [21]. These elements can give consideration to the opening and
closing of the joint.
4.1 Decoupled finite element theory

4.1.1 Theory
When a continuum consists of two separate systems A and B (Fig. 13), the equation
system (Fig. 14) breaks down into two independent blocks.
[K]{a}- {F} = 0 (13)

Numerical Modelling of Tunnels
Figure 13: Coupled system
with
(14)
Figure 14: Decoupled constraint element
Coupling of the region A with the region B is done by means of two-node constraint
elements (Fig. 14). For this purpose, pairs of nodes with separate local displacements
are defined along the joint
( J[))
or in a global system
{ai} = [T]{ af} (W)

with
290 G. Swoboda
(T] = [ co~ ¢ sin ¢ ]

- sm cP cos cP
Further, coupling forces {).} are defined in the pairs of nodes in the local coordinate
system
(17)
+"
Figure l.j: States of the constraint element
The virtual internal work of the different node displacements is:
L - {ad
A;,a = 8 ({a2} L )T {>.} (18)
or, after transforming the displacements into the global system
8 ( (- [T]I (T] ){a} f {>.}

8( [C]{ a} f {>.} (19)
8{a}T(C]T{>.}
The corresponding external work is
(20)
whereby { F} is the nodal forces from the external load. Beyond this, the additional
internal forces also perform virtual work
A;,>.= 8{>.}T {(C]{a}} (21)

The corresponding external work reads
(22)
whereby {ar} is the initial displacerncnts of the nodes.
By equating internal and external work the following results:
11; = A,.
(23)
A;,, -1- l1;,.\ = A,.," -1- A,,.,
and by inserting (19), (20), (21) and (22)
(24)
the stiffness matrix of the constraint joint clement is received
o
[ (C] [CjT ] { a } = { F } (25)
0 A ar
From this it is seen that this element has no clastic stiffness; only additional constraints
are introduced, by means of which coupling can be forced.
In order to be able to orient the local system of coordinates, two additional nodes,
3 and 4, are introduced, as shown in Fig. 15. These are also used to formulate the
failure criterion for decoupling in stresses and not in forces.
The calculation is performed incrementally in load steps, whereby the fracture
criteria are examined. This means that the criteria will be examined for each of these
conditions (I<) starting with the condition of the last iteration (/{- 1). This gives the
following forces and displacements:
joint displacement
coupling forces
(26)
loading
{F}K = {F}K-I + {F}

whereby the added values of the pertinent load steps are without indices.
The friction law developed by Coulomb is used as the failure criterion, which gives
the maximum coupling forces:
shear
(27)
292 G. Swoboda
tension
(~8)
lu these equations, c is the cohesion, <p the angl<• of friction, a z t.lw tension failun•
stresses, L is the iufluence of the clement and d the thickness of l.h<' element. The
decision matrix in Table 1 gives t.he stiffn<·ss matrix d<'scrilwd in [8] for th<' conditions
fix state, free stale and slip stair. The conditions ar<' illust.rat.cd in Fig. 15.
Titble I: D<'cision Tahk
Iteration n Fix Slip Free

11- 1
'\"
N< - Fl'
z >/'N<- FKz
Fix )..K
N
>FK
z
< pK
)..1\
s- t s > FJ(
>..K
t
>..K
N-
< FJ(
z
>..K
N-
< FJ(
z
Slip )..K
N >FK
z
a.FtK < o a.F{ ~ o
Free fl.~< 0 >0

fl.KN-
4.1.2 Mesh generation of jointed rock

Jointed rock is described in the finite element mesh as rock along given possible fault
lines, where constraint elements make it possible for the failure mechanisms described
above to take place. In the framework of interactive mesh generation by means of a
digitizer, the NETDIG program [22] includes the possibility of introducing individual
fault lines after completion of the mesh topology. For this purpose, the node numbers on
both sides of the fault line are reassigned, whereby one node retains its original number
and the second, newly introduced node receives the highest node number available. The
reference nodes in Fig. 15 have to be included in the generation.
Special problems arise in this algorithm when two fault lines intersect. In this case,
appropriate constraint element combinations have to be generated, as shown in Fig. 16,
in order to permit all elements to move along the given lines. A total of up to eight
intersecting lines can be generated in the NETDIG program system.
Constraint Element
'\
I
I
I
/':
Element Joints
:/'::
/':
Figure 16: Crossing fault lines of jointed rock
4.1.3 Calculation of a tunnel in discontinuum

According to the New Austrian Tunneling Method (NATM) [1], driving is performed
in partial excavation, with rapid securement of the cavity using shotcrete. This effects
an activation of the surrounding rock and economic calculation of the cavity's securing
measures. After excavation of the top heading, the primary stresses in the undisturbed
soil or rock are redistributed transverse to the direction of driving and also via the face
by arch action. This redistribution of loads causes displacements in the roof, which in
shallow tunnels can extend to the surface in the form of settlement. The area behind
the face also undergoes displacement that is caused on the one hand by the arch action
following excavation, and on the other hand by the free surface of the face, which does
away with the support action in the tunnel's longitudinal direction. Instead of plane
displacement, plane stress is approximated.
Due to the extensive calculations entailed, the work is generally calculated on a
plane model, whereby the driving is approximated by individual sections transverse to
the tunnel axis. Literature makes reference to three simulation procedures that take
consideration of the influences in the tunnel's longitudinal direction. For example,
there is the possibility of variation of the modulus of elasticity of the shotcrete arch,
or calculation according to the partial-load or support-load method [23], [24].
In the work at hand, simulation is performed using the stiffness reduction method,
that has bceu seen to be the simplest possibility. This method takes account of the.
disphcemcnts caused by strain relief by considering the future excavation area as a
support core. The stresses here arc applied on the excavation's periphery as nodal
forces. The support core, which has a certain defined stiffness, acts as an elastic bed
for the surrounding elements. If the support core has a theoretically infinitely large
stiffness, preliminary relief of strain will preclude any displacements. Thus, the state
of stress from the previous loads remains unchanged. Realistic values for the stiffness
of the support core, that were confirmed in comparisons with calculations using in-situ
nwasurements, arc in a. range of 0.2 to 3 times the rock's module of elasticity.
Fig. 17 shows the excavation procedure as well as the analyzed cross section of the
individual construction steps.
294 G. Swoboda
cr
I
Figure 17: Excavation sequence
Tunnel driving is subsequently examined in a discontinuous rock mass with unin-

terrupted fault system, as shown in Fig. 18. In Fig. 19, construction step 3, a chimney-
like cave-in can be seen, a failure mechanism often observed in such jointed rock. The
load applied by the failure body presses particularly on the foot of the top heading's
shotcrete lining, causing the roof to heave slightly in construction step 4.
4.2 Damage tensor theory

Many rock types can be described as a homogeneous body with a number of cracks
located in one plane. The planes have a regular distance and can be overloaded by
other planes with different orientations. The damage tensor theory of rock is part of
the system of a continuum theory. The basis is the geological investigation of the fault
system of the rock.
4.2.1 The damage tensor- a linear representation
Let us assume that we have a three-dimensional element with one plane of cracks
(Fig.20). In the plane of cracks are the axes x 1 and x 2 of the cartesian coordinate
system xll x 2 and x 3. The x 3 axis coincides with the direction of {n }, the orientation
of the damaged rock. The share of damage or the remaining rock bridges are described
by the damage factor w [33].
If we assume a general plane in the undamaged rock ABC with the area S, the
normal vector of this plane is
(29)
-
--
~
"----
Figure 18: Tunnel in rock with horizontal and vertical jointing, arrangement of con-
straint elements.
In this relation, {v} is the normal vector and {e;} (i = 1,···,3) are the canonical
base vectors of the chosen coordinate system.
The effect of damage can be considered by reducing the effective area of OAB, in
the plane of cracks, by the factor (1- w), while the area of the planes OBC and OAC
remains unchanged.
Thus the reduced effective area will be described by OA* B* . Geometrically, we
change from the undamaged plane ABC to the damaged plane A* B*C* by reduc-
ing sides OA and OB by the factor ~while side OC is enlarged by the factor
1/ ~- The area of the thus obtained triangle A* B*C* is equal to s•, and its
normal vector is { v•}.
The normal vector of the damaged plane is
{S*} = S*{v*} = SJ{ed + S2{ez} + (1 - w) · S3{e3} (30)

and with the help of equation (29) we find
{S*} = S{v}- wS · [{e3 } 0 {e3 }]· {v}

whereby 0 is the tensor product:
[A]= [{n} 0 {n}] => A;i = n; · n i

If a tensor [!l], the so-called damage tensor, is defined as
[!l] = w · [{n} 0 {n}]
then {s·} is
{S*} = S ·[I- !l] · {v} =[I- !l] · {S} (31)

296 G. Swoboda
E
E
0
"'
'
'"' '
100 mm
' ......................_____-..,;
Figure 19: Fault system and displacements in construction steps 3 and 4
4.2.2 The damage tensor in rock mechanics

Let us consider a cube with a side length L and a volume V (Fig. 12) that has been
cut out of the rock. Through this cube we cut several parallel planes at intervals of l,
namely the mean distance of the individual joint planes. Next, we imagine the cracks
in the rock as being in these planes. The rock between them is homogeneous and
t
isotropic. Since represents the number of planes, the largest possible total area of
all cracks in the considered joint pattern is calculated as follows:
2 L V
S=L · - = - (32)
1 1
The k-te crack of this joint pattern ("i ") has the area ak,i and the normal vector { n ;}.
Its related area is thus
k = 1, ... ,N (33)
Now the damage tensor of this special crack can be calculated:
[nk,;] = wk,i [{n;} ® {ni}] (34)

The following is thus true for the entire joint pattern:
plane of
damage
Figure 20: A continuum theory of creep and creep damage
1·
L
N N
[0;] = l:)Ok,i] = ~ [{n;} ® {n;}] ak,i (35)
k=1 k=1
If the cube is cut through by several joint patterns, [Ototatl is obtained as the sum
of the individual damage tensors:
M 1 M N
[Ototatl = L[O;] = V L l; [{n;} ® {n;}] L ak,i (36)
•=1 •=1 k=l
Based on the damage factor of each set of cracks the damage tensor can also be
found by the following relation
N
[Oeotatl =L wk[{nk} ® {nk}]
k=1
where N is the number of sets of cracks.
4.2.3 Introduction of net stresses

Let us consider a 3D element V containing several cracks parallel to the y axis. The
separation factor is w (0 :::; w < 1). A normal stress of (J"x is at work on the element.
For reasons of equilibrium, a larger stress must be present inside the element (Point 1)
since, due to the cracks, the effective resisting area is smaller than the area on which
the load is imposed. We term these strains net strain (]";. They are expressed as
(1";(1-w) = (J"x (37)

or
298 G. Swoboda
z z
L--J~"
L
r
f L y
l X X
Figure 21: Joint system and joint planes
{38)
When these net stresses are applied to an isotropic, homogeneous element, the
clement is deformed under the given load, as the below element with cracks shows.
y
I I I I I I I
•1 'v'
I I I I I I I
I I I I I I I
)Risse
I I I I I I I
~--------------------------------7 X
Figme 22: Damaged rock
4.2.4 Mathematical formulation

011cc the st.r<•s:><'s 011 a homogeneous and isotl'Opic system have been calc.ula.ted, the net
stresses [a*] are determined. If [a] is the Cauchy stress tensor and {P} the resulting
vl'dor of force 011 !.he tetrahedron surface ABC in Fig. 18, the following holds true:
{p} =[a]· {JJ} (39)
{P} = S · {p} =.'-.'·[a]· {11} =[a]· {S}
{P} =[a*]· {S*} => [a*]· {S*} =[a]· {S}
=> [a*] = [a]· [I- nr 1 (40)

Since the effective area of the element is reduced by the cracks in the rock, these
net stresses must be larger than the stresses in a homogeneous system. The linear
representation [I- nr 1 describes these mechanical effects of the joints in the rock.
Now, however, it must be determined whether the stresses acting normally on the
cracks are positive or negative. Tensile stresses must be absorbed exclusively by the
effective area. Thus the net stresses are precisely calculated according to the above-
given formula. The compressive stresses, in contrast, can also be transferred through
the joint since some of the cracks are closed. The effective area is thus larger than
that described by the joint tensor. The same is also true for the shear stresses acting
parallel to the joint. The prevailing state of stress thus effects the size of the effective
area of any element.
If a matrix [T] is defined, whose columns are the eigenvector {</>;} of the damage
tensor [n] and if a diagonal matrix [0'] contains the eigenvalues w;, namely
[T] = [ {</>I} { </>2} {</>a}]; 1</>d = 1 i = 1, ... ,3
the following holds true:
Using [T] the Cauchy stress tensor [a] describing the stresses in the global system of
coordinates can be transformed to the direction of the damage tensor's principal axes:
[a'] = [Tf[a][T] (·H)

300 G. Swoboda
In order to be able to express the above cases in mathematical terms, [a'] must be
broken down into a tensor [a~] that contains the normal stresses and a tensor [a;] with
only shear stresses:
(42)
In its explicit form this relation reads as follows:
If the crack area were completely plane, no shear force could be transferred across
it. In this case the effective area would have to be calculated with [I- 0']. Since,
however, its surface is very uneven, a rough effect occurs that permits part of the shear
stresses to be absorbed. For this reason the effective area is modified as follows:
[I - 0'] f-t [I - CtO'J ;
As already stated, this is correspondingly true for compressive stresses:
[I- O'J f-t [I- CnO'J in case a < 0; 0 ::; Cn ::; 1
Contrarily, no modification is necessary for tensile stresses. The net stresses trans-
formed to the principal axis of the damage tensor are thus calculated as follows:
whereby [H](-) is an operator that permits us to distinguish between positive and

negative normal stresses:
zn case a'·
IJ >0
zn case a~.
lJ-
<0
Now the stresses [a*'] must be transformed to the original global system of coordi-
nates:
If [a] is the originally calculated Cauchy stresses, the differential stress tensor [7f] is
defined as follows:
4.2.5 Analysis in a finite element program

The equation to be solved in the framework of the principle of virtual work generally
reads as follows:
(44)
or
[K]{u} = {F}
After inserting the above relation, the equation reads:
or
[K]{u} = {F} + {F*}

It is noteworthy that the stiffness matrix [K] only takes the characteristic values
for isotropic and homogeneous rock. All mechanical effects of the cracks are covered
by the term {F*}.
(46)
When one single element is isolated from the entire system, {F*}e is the joint load
vector that causes the element to deform according to the differential stresses [t/1], or
that causes these differential stresses to occur in the element.
4.2.6 Practical procedure

(I) The equation [K]{ u'} = {F} is solved and the stresses of an isotropic sysstem are
calculated from the joint displacements {u'}.
(2) The· net. stresses [a*] arc calcula.tl'd from the Cauchy stress tensor [a] (under con-
sideration of the mentioned influences of the particular state of stress) with the
help of !.he damage tensor [OJ.
(:!} The differential stress tensor [1/•] = [a*]- [a] is calculated followed by the joint load
vcd.or {F*}.
('I) Once the syst.cm in question is statically determined, the final joint displacements
{u} of the jointed systems are obtained by solving the equation
[K]{u} = {F} + {F*}

302 G. Swoboda
(5) In the other case (see Point 4), solving the equation will only provide an approxima-
tion for the desired joint displacements. Since the introduction of jointed elements
to the system (namely in the form of an additional joint load vector {F*}) is equiv-
alent to changing the stiffness matrix, the stresses in the system are differently
distributed than in the isotropic calculation (see Point 1) that was the prereq-
uisite for calculation of the joint load vector. For this reason, the desired joint
displacement must be calculated by iterative means. The stress tensor [&], which
corresponds to a load of {F} + {F*}, is calculated according to Point 4. Since the
joint load vector {F*}, however, is not present in the actual system, its part of the
total stress [&] must not be considered in the new calculation of the net stresses.
For this reason, the differential stress tensor [~old] (from which {F*} was originally
calculated) must be subtracted from [&]. This permits the new net stresses, the
new differential stresses and thus the new joint load vector to be calculated as
follows:
[a] = [&]- [~old]
[a*] = f([a], [11])
[4-new] ·- [a*] - [a]
{F*}new = J({~new})
(6) The strain condition present can now be given on the element for every point with
the help of the [B] matrix:
{c:}(x,y,z) = [BJ(x,y,z){u}
This strain condition is independent of whether or not cracks exist in the element
in question because the strain condition depends solely on the joint displacement and
the chosen shape function.
Once the matrix of elasticity [D]e is known, the stresses at every point are also
known:
{a}e = {a}(.,,y,z) = [Dr{c:}e = [nnar{u}
As far as the normal stresses are concerned, these must occur in this form in the
rock bridges of the jointed rock. The shear stresses, however, are still those of the
isotropic system, but the desired net shear stresses can be calculated from these very
simply.
If [a] is the stress tensor equivalent to {a V, then
[a'] = [T]1'[a)[T]
[a'] = [a~] + [a~]
[a*'] = [a;]· [I- Ct11T 1 + [a~]
[a*] = [T][a*'][TjT
Unfortunately, this procedure does not always bring the desired results. Some
systems have been seen to have a critical separation factor of the joints, above which
the solution no longer converges.
4.2. 7 Calculation of the damage tensor

The sphere diagram is generally accepted for geomechanical representation of the fault
system, namely the sum of all the 3D data of a network of joint surfaces. For this
reason, calculation of the damage tensor proceeds from the assumption that the fault
system is given in this form.
On first inspection the network of joints appears to be very complicated. When it
is statistically analyzed by a suitable procedure, however, it is seen to contain main
directions. Individual joints can be ascribed to joint patterns with various orientations.
A joint pattern is understood as a number of joint planes that mainly have the same
(or a similar) position in space.
Minen any several joint patterns there are (usually) two or three that are more
pronounced in the frequency of their individual components and thus statistically sig-
nificant as compared to other joint patterns or individual joints. These are known as
main joint patterns and are mathematically defined by the damage tensor.
Figure 23: Sphere diagram [1]
In the sphere diagram (or in its representation) every joint is clearly defined by its
pencf.ration point diagram (Fig. 21 ). If we imagine the plane as being pushed to the
304 G. Swoboda
center of the sphere and if we take the normal vector to this plane, the desired point
is defined as the one where the surface of the sphere is penetrated. This point can be
more simply obtained by having the plane approach the sphere only until it touches
the sphere. This point is the desired penetration point.
In order to calculate the damage tensor, we need the penetration point diagram of
a maximum of three main joint patterns, whereby the particular damage tensor is used
to determine these. The tensor that describes the entire network of joints is the sum
of these. N
Figure 24: Orientation of the joint system
First, we need the normal vector of the main joint pattern in question, which is
easily obtained as shown in Fig. 22:
~{
- sin f3 sin a }
{n} - cosf3 (47)
sinf3 cos a
It is automatically assumed that the tunnel axis runs north-south. If this is not the
case, the sphere diagram must be turned so that a gives the "strike" of the main joint
pattern with reference to the tunnel axis.
{n} and the separation factor w give the 3-D damage tensor [0]:
[n] = w. [{n} ® {n}] (48)

If the rock is under a plane strain condition, it should also be described as a plane
problem since this is simpler than a 3D calculation. [OJ must therefore be converted to
a 2D damage tensor that describes rock behavior as similar as possible to the original
tensor (in a 3D calculation under plane strain condition).
For a unit stress condition
[a]= [ 1 1 0 l
the pertinent net stresses [a*] are ca.kulat.cd analogous t.o Cha.pl.t'r -I.~A .:
[a*] = [a]· [I- nr• (49)
From the 3D stress tensor [a*], those components arc isolated that exclusively cause
plane displacements in the area. in question:
The desired plane damage tensor is then given as follows:
A check of the formula's accuracy is simple: in the 2D tensor the unit strain con-
dition is given as [1]2x 2 • The net stresses [u*] are thus given as:
[i7*] = [u] · [I- nr 1 = [I] · [I- nr 1 = [I- nr 1 (50)

By inserting [0]2x 2 , the desired net stresses are obtained for the chosen plane strain
condition:
[i7*] = [a*hx2
4.2.8 Comparison of the results of this work with conventional methods

It is customary for the rock with its joints to be simulated by a material with an
isotropic material behavior. For this reason, the results of such calculations are em-
ployed here in order to get a feel for the mechanical effect of the damage tensor.
The compressive test described in Fig. 23 is an example of this. The joint planes
are inclined 45° away from the x axis and have a mean separation factor of 0.3.
The damage tensor thus reads:
cos45° }
{n} = { sin45°
[OJ= w · [{n} ® {n}] = [ ~:~~ ~:~~]

The finite element network has the following appareance:
306 G. Swoboda
p a 10000 kNtm• p •10000 kN/m 2
IOIOI [l]Iil]
, , , E ·10000 kN/m~f -
'' ' y •0.25
'' ' ' E nE

~ '-.2'____
llb
'''' E
0 ' ' ' ' 003 dy ,/\(45
l ~ ~'~': ~d·
--·· /~
m
·---
r- 0.2~
1 X "'(45°
r- 0.2~
I X
Figure 25: Comparison between damage theory and anisotropy
Fig. 25 shows the results of a damage tensor analysis for various combinations of
C1 and Cn. Since the net shear stresses do not cause any additional displacements,
the factor C1 has no influence on the size of the displacements. In this example these
are solely dependent on Cn and [0]. The maximum values of the displacements are
summarized in a table at the end of this work.
The transverse isotropic model roughly describes the displacement behavior of the
jointed rock when using the factor (1 - Cnw) for n (the relation of the moduli of
elasticity). The stresses calculated are, however, incorrect because the calculations
were based entirely on different prerequisites.
0.1 0.2
0.4 0.4
0.3 0.3
0.2 0.2
0.1 0.1
0.1 0.2 Ct • 1.0,Cn • 1.0 n • 0.7
Figure 26: Numerical model Figure 27: Displacement for different

damage constants [25]
5 Three-dimensional numerical models

Particularly in t.uml<'ls with large displacements, however, numerous aspects of the
construction works cannot b(• resolved without examining the static influence of the
driving fa.ce. The inOucnce of time, iu particular, cannot be realistically accounted
for in a 2D model. As shown in Fig. 28, complex interaction between driving velocity
v, rheological displacement of the ground s and the time-dependent change in the
elasticity modulus of the shotcrete combined with its viscoelasticity must be polued in
order to define the stresses on the tunnel face. Since these processes take place directly
behind the face, the use of 3D models is an absolute necessity.
5.1 Numerical modelling of the shotcrete lining

A sandwich shell element (Fig. 29) is used to simulate the shotcrete lining. Its multi-
layer system permits very good simulation of shotcrete application.
The isoparametric sandwich shell with quadratic shape function (ISSQ) [30] is based
on Ahmad's concept for degenerated isoparametric elements [29], whereby the element
geometry is described by the mid-surface. According to Fig. 30, every node of the shell
can be determined by means of the mid-surface's coordinates, the shell thickness t and
the v3 vector, a vector normal to mid-surface:
(51)
with
308 G. Swoboda
Interpolation of the coordinates from the element coordinates is performed with the
help of the shape function
x(e,7J,O}
{ y(e,7J,() =
z(e, 7], ()
(52)
whereby
{a}® {b} = {a1b1 a2b2 a3b:3 · • ·V
Interpolation of the displacements is performed in the same way as for the coordi-
nates
(53)
whereby a and f3 represent the rotations around the local axes XL and YL and thus
the displacement from the element node displacements '
u(e, 77, () }
{ v(e, 7], () =
w(e,7],()
(54)
The strains in the element are as follows
f.xL UIO,xL aXL

f.yL vlo,yL -{3YL
{€£} = "'fxLYL = I
U O,yL
+ V IO,xL +zL aYL- f3xL (55)
"'fxLZL Wlo,xL +a 0
"'fYLZL w 1o,YL- f3 0
1 6uo
U O,yL = 6yL
s = f(t)
v = f(t) E = f(t)
Figure 28: Influences of time: driving velocity, viscoplasticity of rock, viscoelasticity

of shotcrete
(
Figure 29: Sandwich shell element for simulation of the shotcrete lining
The displacement function used belongs to the family of Lagrange Polynomes and
is or(\ contiunity
N• ~((;(1 + ((;) 1]11i (1 + 1J17d i = 1,2,3,4

N; ~((;(1 + ((;)(1 - 77 2 ) i = 6.8 (56)
N• h77;(1 + 1J11;)(1- ( 2 ) i = 5, 7
N• (1- (2)(1- 112) i = 9
Using (56) and (55) , the strains can be found by differentiation:

310 G. Swoboda
(x,y ,z)bottom
Figure 30: Mapping of three-dimensional geometry to the mid-surface
{uo}
{vo}
{tL}=[B] {wo} (57)
{a}
{~}
[B] =
bT1 0 0
0 bT2 0
bT2 bT1 0
0 0 bT1
.eT
(bf 0
0 (bf
(bi (bf
bT3 0
0
0
0
0
I
· 0T · 2 ·
[ "'• 0 t
V1y @t
V1z@ t
-v, 0
-V2y@ t
-V2z ® t
t l (58)
0 0 b{ 0 bT
3 0
3n column• 2n column•
n is the number of integration points, and t the element thickness in the integration
points. Furthermore,
{bl}T = (N]~<fn + (N]~,J12

{b2}T = (NJ~.d21 + [N]~,J22 (59)
{b3}T = [N]~ /33
l
The coefficient fii results from
[F] = [ ~:~ ~:~ ~:: = (Q]T[Jt1 (60)

ht !32 !33
whereby [0] is the direction cosine matrix and [J]- 1 the inverted Jacobian matrix.
The stiffness matrix of the multi-layered shell elements is thus found as follows
Shotcrete:
Vfscoplastic Model
E,,,~ ·~
Time Time
Figure 31: Multilayered shotcrete- soil system
[k] = {;
I 1 1[B]T[Q]<kl[B] st
111111
1
detiJI de d1J d(k {61)
The matrix [QJ(k) is the stress-strain equation of the individual layers. Integration
is now only performed with the shell thickness sk and is subsequently added for all
layers l.
In order to account for the viscoelasticity of the shotcrete [:ll], the element load
vector { Fvd} must be determined from the strain increment {Alvei}
{AfveL} = (1- e-
Eb,ta
~ ·
At
)(2Jt-- {t:ve~(to)})
3 {uft}
b,ta
{62)
This is performed with the following known equation
1~
{M"") ~ E. ~ {t~\ }~ f.,inJT[D]{t>.wel) dVol (63)

My
whereby TJ represents the viscosity and Eb,ta the time-dependent modulus of elasticity
of the shotcrete.
The following equations for the viscoeleastic additional forces can be set up for the
existing element consisting of /layers:
312 G. Swoboda
Numerical Model
81 E E 1.50 ·10' kNfm'
I'= 0.25
1 = 20.0 kN/m'
B2 E = 0.50 · 10' kNfm'

I' =0.33
1 = 21.0 kN/m'
c = 1.0 kNfm'
., = 27.5°
83 E = 0.75 ·10' kN/m'

I'= 0.28
7 = 21.0 kNfm'
c = 1.0 kNfm'
., = 27.5°
E = 0.50 · 10' kN/m'

I'= 0.40
1 = 21.0 kNfm'
~ ~ 10000 Tag< · k.'V/m'
c = 1.0 kNfm'
n= l
., = 27.5°
Figure 32: Finite element model
{M}
(64b)
The load {tl.Fv.t} that simulates the viscoelastic displacements of the shot crete is
redefined in every time step and added to the static loads.
Fig. 31 shows the application of the multi-layered model. The shotcrete lining
was applied in three layers in this case, whereby each layer has different viscoelastic
properties Eb,ta and 1]b,ta according to the time of application to, t 1 or t 2 •
The time-dependent displacements of the soil were treated with the flow rule after
Desai [32]. The two regions were coupled by shifting the solid element and the shell
element u, v, w . The rotations a: and fJ were not coupled with the soil region.
5.2 Simulation of tunnel excavation

The method used to simulate excavation was not the simplified method, according
to which the excavation is simulated in a numerical model by displacing the model
without changing the system. This method is not suitable for simulating very jointed
rock or for analyzing special fault zones.
Working from a plain strain condition, the excavation is simulated step-by-step by
eliminating the clements. The assumption is that the elements concerned in one round
are those from n to m. This means that the modulus of elasticity is assumed to be
zero for these elements (65) and the stresses in the eliminated elements integrated after
(GG) and applied as extreme loads { Fex}
E; ~ 10-2o 1. = n, ni (65)
{Fe.r} = f])Bf{cro}d11 (66)

1.=n
The application was performed on a one-track tunnel profile (Fig. 32) for which
an appropriate time-displacement curve was available. The soil was loose, mainly
quartenary deposits. The dri\·ing velocity was 0.666 m/day. The viscosity of the soil
was modified as a parameter. The values "7 = 10000,25000 and 50000 days· KNfm 2
were calculated. The exponc11t for the viscosity rule was n = 1.0.
Figur<' :l:l: C01npa.rison of measur-ements and viscoplast.ic calculations. Roof displace-

nwnt.
314 G. Swoboda
The following characteri'itic values were chosen for the shotcrete:
Time Eb Eb,ta .,.,

[KN/m2] [KN/m2] [days· KN/m2]
0- 2days 15. 106 6. 106 6 ·108
> 2days 25 ·106 10. 106 10. 108
Fig. 33 shows a comparison of the displacements calculated and those measured in

the roof. The reason why the displacements in the numerical model have a pronounced
step function is that in the excavation, for example of the top heading, this process of
several hours takes place suddenly in the numerical model instead of slowly, step by
step under actual site conditions. Furthermore, the measurement curves are smoothed
since the measurements are not performed more often than twice a day. From the
longitudinal displacement line for the roof and invert, as shown in Fig. 34, it can be
seen that the displacements resulting from excavation start just in front of the face.
Thus, the redistribution of stresses takes place in a very small region before and after
the face.
10 time-steps = 1 day
.........•.·············. ·············...··.-.
' r
30 time-steps = 3 days
. ·.·.· •'•'•' =·=·=··=·=···..·:-;.::·············:·:···:·:·.·=········· ···.···.·=-=····.-.•,•,•.-•.·.··:
50 time-steps = 5 days
·.···=······:=:::-:-:;:···:;::·-:-:.:-.···:·:···...·····-.·····.·.·.············ ;..............
Figure 34: Excavation sequence and longitudinal roof and invert displacement
For various displacement velocities of the soil, and thus various viscosities, a com--
parison was made of the displacements in the roof (Fig. 35) and the longitudinal normal
forces Ny (Fig. 3(i). Particularly in th<' cas<' of nom1a.l forn•s, th<' displan•tw•nt. vt'locities
can sometimes cause very difl"erent forc<'s to O<TIIr. 'l'ht' longitudinal support effect is
especially pronounced in the base of t.h<' crown. Th<'s<' for-ct•s rt'::mlt from the cantilever
dfecl. of t.he lining in t.he t.op lwa.ding.
w [mm]
20
11 = 50000
15
11 = 25000
10 11 = 10000
5
0 10 20 30
Time-Steps
Figure 35: Parameter study of roof displacement. Influence of the soil viscosity
m
N
N N
I
I
_fj?J"'
_-:,22
_ze'3
-251
-286
-274
I I I I
~ ~ al $
17 = 50000 Days kN/m2 17 = 25000 Days kN/m2 17 = 10000 Days kNfm2
Figure 36: Longitudinal normal forces in the measuring cross section

316 G. Swoboda
References
[1] L. MULLER, Salzburg: Der Felsbau, III. Bd., Tunnelbau, F.Enke 1987
[2] M. BAUDENDISTEL: Zur Bemessung von Tunnelauskleidungen in wenig festem
Gebirge. Rock Mechanics, Suppl. 2, 279-312 (1973)
[3] G. SWOBODA, F. LAABMAYR: Beitrag zur Weiterentwicklung der Berechnung
flachliegender Tunnclbauten im Lockergestein. In Lessmann: "Moderner Tun-
nelbau bei der Miinchner U-Bahn". Wien: Springer 1978. Ubersetzung: Peking
1980, Tokio 1982.
[4] F. LAABMAYR, G. SWOBODA: Zusammenhang zwischen elektronischer Be-
rechnung und Messung; Stand der Entwicklung fiir seichtliegende Tunnel. Rock
Mechanics, Supplement 8 (1979), S 29 - 42.
[5] G. SWOBODA: Finite Element Analysis of the New Austrian Tunneling Method
(NATM). Proceedings of the 3rd International Conf. on Numerical Methods in
Geomechanics, Aachen 1979, p 581-586.
[6] G. SWOBODA: Special Problems During the Geomechanical Analysis of Tunnels.
Proceedings of the 4th Conference on Numerical Methods in Geomechanics 1982,
Edmonton, p 605-609.
[7] O.C. ZIENKIEWICZ: The Finite Element Method. London: McGraw Hill
1982.
[8] K. GIRKMANN: Flachentragwerke, Wien: Springer 1983.
[9} K. GIRKMANN: Zum Halbraumproblem von Michell. lng. Archiv 14, (1943)
106-112.
[10] K. KOVARI, C. AMSTAD, J. KOPPEL: Neue Entwicklungen in der Instru-
mentierung von Untertagbauten und anderen geotechnischen Konstruktionen,
Schweiz, Ingenieur und Architekt, H 41. (1979).
(11] G. SWOBODA, W. MERTZ, G. BEER: Application of coupled FEM-BEM anal-
ysis for three-dimensional tunnel analysis. in: Boundary Elements (ed. Q.H.
Du), p. 537-550, Peking: Pergamon Press, 1986.
[12] G. SWOBODA, W. MERTZ, G. BEER: Rheological analysis of tunnel excavation
by means of coupled finite Element (FEM) -boundary element (BEM) analysis.
Numerical and Analytical Methods in Geomechanics. (in press)
[13] G. SWOBODA: Boundary Element, ein Bindeglied zwischen analy- tischen und
numerischen Methoden. Baupraxis- Baustatik, Stuttgart 1987.
[14] G. SWOBODA, F. LAABMAYR, I. MADER: Grundlagen und Entwicklung bei
Entwurf und Berechnung im seichtliegenden Tunnelbau. Teil2, Felsbau 4, (1986)
S.184-187.
[15] R.E.GOODMAN, R.L.TAYLOR and T.BREKKE: A model for the mechanics of

jointed rock. Journal of Soil Mechanics and Foundation Division, ASCE 94, p.
637-658. (1968).
[16] C.S. DESAI, M.M. ZAMAN, J.G. LIGHTNER, H.J. SIRIWARDANE: Thin-layer
element for interfaces and joints. Internal Journal for Numerical and Analytical
Methods in Geomechanics S, p. 19-43 (1984).
[17] W. \VITTKE. Static analysis for underground openings in jointed rock. in Ch.S.
Desai, J.T. Christian: Numerical Methods in Geotechnical Engineering. :McGraw
Hill, New York 1977.
[1 8] O.C. ZIENKIEWICZ, C. DULLAGE: Ana.lysis of non-linear problems in rock
mechanics with particular reference to jointed rock systems. Proceedings of the
2nd International Congress on Rock Mechanics Sec., 8-14 (1970).
[19] T. MAINI, P. CUNDALL, J. MARTI, P. BERESFORD, N. LAST, M. ASGIAN:
Computer modelling of jointed rock mass. Technical Reprint N-78-4. U.S. Army
Engineers, Vicksburg, 1978.
[20] T. ASAI, M. NISIIIMURA, T. SAITO, l\L TERADA: Effects of rock bolting in
discontinuous rock mass. 5th International Conference on Numerical Methods in
Geomcchanics, Nagoya, p. 1273-1280 (1985).
[21] l\LG. KATONA: A simple contact-friction interface clement with application to
buried culverts. Internal .Journal for Numerical and Analytical Methods in Ge-
omecha.nics, 7, p. 371-38·1 (1983).
[22] NETDIG: Digitalisierung von Finite-Element-Netzen V 2.1 Manual University of
Innsbruck 1985.
[23] M. BAUDENDISTEL: Zum Eutwurf von Tuuneln mit grofiem Ausbruchsquer-
sclmitt. Rock 1\1cchanics, Supplement 8, p. 75-100 (1979).
[24] K. SCHIKORA, T. FIN1\: Berechnungsmethoden moderner bergma.nnischer Bau-
weisen beim U-I3almbau. Bauingenienr, 57, p. 193-198 (1982).
[25] T. KYOYA, Y. ICIIII\AWA, T. KAWAMOTO: A damage mechanics theory for
discontinous rock mass. 5th Int. Conf. on Num. 1\lcth. in Geom., Nagoya, 1985,
pA69-480
[26] C. SWOBODA, G. BEER: SUidtischer Tunnelbau- Rechenmodelle und Resultat-
interpretation als Grundlagc fiir Planung und Bauausfiihrung. Finite Elemente
i 11 dcr Bau praxis, l\1 iill(·hen, 1984.
[27] G. BEER, G. SWOBODA: On the efficient Analysis of Shallow Tunnels. Comp.
and Geomech. 1985, 1, pp. 15-31.
[28] G. SWOBODA: FINAL, Finite Element Analyse lincarer und nichtlinearer Struk-
turen, Programmvcrsion 6.1. Universitat lnnsbruck, 1988.
318 G. Swoboda
[29] S. AHMAD, B. IRONS, O.C. ZIENKIEWICZ: Analysis of Thick and Thin Shell
Structures by Curved Finite Elements. Inter. J. Num. Meth. Eng., 1970, 2, pp.
419 - 451.
[30] A. SCHMID: Beitrag zur Berechnung geschichteter Schalcntragwerke mittels dcr
Methode der Finiten Elemente. Dissertation, Innsbruck, 1988.
[31] W. MERTZ: Numerische Modelle drei-dimensionaler plastischer, visko-plastischer
und visko-elastischer Spannungen und Verformungen in der Geomcchanik unter
besonderer Beriicksichtigung seichtliegender Tunnel. Dissertation, lnnsbruck,
1988.
[32] C.S. DESAI: Somasundaram S., Frantziskonis, G., A hierarchical approach for
constitutive modelling of geological materials. Inter. Jrnl. Num. and Anal.
Meth in Geom., 1986, 11.
[33] F. PACHER: Kennziffern des Flachengefiiges. Geol. und Bauw. 24, 223ff, 1959.

Numerical Modelling of Tunnels: G. Swoboda University of Innsbruck, Innsbruck, Austria

Uploaded by

Copyright:

Available Formats

Numerical Modelling of Tunnels: G. Swoboda University of Innsbruck, Innsbruck, Austria

Uploaded by

Document Information

Original Title

Copyright

Available Formats

Share this document

Share or Embed Document

Sharing Options

Did you find this document useful?

Is this content inappropriate?

Copyright:

Available Formats

Numerical Modelling of Tunnels: G. Swoboda University of Innsbruck, Innsbruck, Austria

Uploaded by

Copyright:

Available Formats

NUMERICAL MODELLING OF TUNNELS

C. S. Desai et al. (eds.), Numerical Methods and Constitutive Modelling in Geomechanics

~-----*":::...._~---\-'~f-;~ .---"--- : ----------

'"- E = 40000 kN/m 2

Figure 1: Numerical model for a subway tunnel

Figure 2: Complex cross section for subway tunnels

Stage 0 Stage 1 Construction Construction

Figure 3: Full load model

Construction Stage 0 Construction Stage A Construction Stage 8

Stiffness reduction Support load

Figure 4: Two-dimensional simulation of the three-dimensional effect

3.1 Stress redistribution due to excavation

Displacement in front of face Additional displacement

Cross section A Cross section B

(c) IEs· <!E I

excavated area ---,

Figure 7: Line load in two- and three-dimensional halfspace

v = (1 + fl) F. [(1- fl) ln x2 + y2 - y2 ] (6)

(l+J-L)P [ l b+/x2+y2+b2 2 bxz ]

whereby b is the length of the line load.

Figure 8: Displacement of a. line load on a. three-dimensional halfspace

Figure !l: Nonlinearity in the origin of the stress-strain curve

'!'he reason why it is so difficult to obtain an unmistakable solution is that the

' ......... _____ ,..

Figure 10: l'vlodification of the boundary condition

3.2 Water pressure

{s} = [D]{a} = [D] { (9)

Pore water pressure

Figure 11: Computational model for the loading case water

iiy = b' + 1.0) h =

{c} = [D]{~a9 } = [D] { (11)

(ii) Pore-water pressure on a joint face

O"T - ~0"1\",T >0 -t ET = 0 (12)

STAGE I STAGE II Lowered ground-water table

STAGE Ill Filling with water

Figure 12: Loading case water

4.1 Decoupled finite element theory

[K]{a}- {F} = 0 (13)

Figure 13: Coupled system

Figure 14: Decoupled constraint element

{ai} = [T]{ af} (W)

(T] = [ co~ ¢ sin ¢ ]

Figure l.j: States of the constraint element

The virtual internal work of the different node displacements is:

8 ( (- [T]I (T] ){a} f {>.}

The corresponding external work is

A;,>.= 8{>.}T {(C]{a}} (21)

the stiffness matrix of the constraint joint clement is received

{F}K = {F}K-I + {F}

Titble I: D<'cision Tahk

Iteration n Fix Slip Free

Free fl.~< 0 >0

4.1.2 Mesh generation of jointed rock

Figure 16: Crossing fault lines of jointed rock

4.1.3 Calculation of a tunnel in discontinuum

Figure 17: Excavation sequence

Tunnel driving is subsequently examined in a discontinuous rock mass with unin-

4.2 Damage tensor theory

4.2.1 The damage tensor- a linear representation

{S} = S{v*} = SJ{ed + S2{ez} + (1 - w) · S3{e3} (30)

{P} =[a]· {S} => [a]· {S} =[a]· {S}