BENDING MOMENTS AND NORMAL FORCES IN TUNNEL LININGS

TITLE OF THE PAPER IN FRENCH

P.G. Bonnier 1, S.C. Möller 2, P.A. Vermeer 2

1
Plaxis B.V., Delft Netherlands
2
Institute of Geotechnical Engineering, University of Stuttgart, Germany

ABSTRACT - Three-dimensional finite-element analyses are considered for the purpose of

predicting loads on linings and surface settlements. First of all computational results from proper
step-by-step tunnel installations are shown. Because such 3D analyses are extremely computer
time consuming, for loads on linings it is advocated to use a 2D FEM analysis. A new method is
presented where the so-called β-value can be determined from a simple 3D analysis.

RESUME – Abstract in French

1. Introduction

At present tunnels in praxis are mostly analysed by the use of 2D finite element computations,
because 3D analyses are considered to be very time consuming. Therefore consulting engineers
will only perform 3D analyses when facing complex geometries, e.g. tunnel joints or underground
stations, but not for straight ahead tunnelling. In fact 3D analyses can become very elaborate and
it is useful to retain simple 2D approaches. However it is smart to supplement 2D approaches by
some 3D calculations. For analysing tunnel construction it will be differentiated between the three
main focuses (Vermeer, 2001), i.e.

A) tunnel heading stability, B) loads on linings and C) surface settlements.

For tunnel heading stability the reader is referred to Vermeer and Ruse (2001) and this topic will
not be discussed. In this paper the focus will mainly be on loads on linings and surface
settlements. First of all loads on tunnel linings from a full 3D excavation of a circular NATM (New
Austrian Tunnel Method) tunnel will be shown, where the bending moments and the normal
forces are found to have a non constant, oscillating zigzag pattern. Hereafter results from a
similar 2D analysis will be compared to those of the 3D. In order to perform a 2D analysis which
matches the results from a full 3D analysis a new method will be presented, in which the so-called
β-value can be determined from a very simple 3D analysis, which does not need excessive
computer run time. This method can also be used to describe the complete shape of the 3D
settlement trough.

2. Bending moments and normal forces

An important topic of tunnelling are loads on linings. To investigate realistic values of bending
moments and normal forces a full 3D analysis is needed. In order to calculate realistic values
from a circular tunnel we divided a block of 100x40x28m into 8840 volume elements with a total
of 26809 nodes (Figure 1). For the parameters of the Mohr-Coulomb (MC) model, we took a
Young’s Modulus of E=42Mpa, a Poisson’s Ratio of ν=0.25, a cohesion of c=20kPa, a friction
angle of ϕ=20o, a dilatancy angle of ψ=0 and K0=1-sinϕ. The NATM tunnel with a diameter of 8m
and a cover of 16m was modelled in a symmetric half with an unsupported excavation length of
2m. Each computational phase consists thus of d=2m of excavation, in which one slice of soil
 Figure 1: Shadings of vertical displacements after 80m of stepwise excavation. The steady-state
settlement of s=4.5cm is reached after an excavation of 35m.

element is switched off. Within the same phase a ring of lining elements is switched on to support
the previous excavation. The shotcrete lining has a thickness of 30cm, a Young’s Modulus of
20MPa and a Poisson’s Ratio of ν=0. This way of modelling a NATM tunnel will be referred to as
‘step-by-step installation’ (Wittke, 1984).
Figure 2 presents normal forces after 80m of ‘step-by-step installation’. Within a single
shotcrete ring of width d=2m there is a sharp drop of the normal force from about 1400kN/m at
the front of the ring down to about 200kN/m at the back of the ring. This is logical as the
unsupported tunnel head is arching on the front and not on the back of this tunnel segment. The
average normal force, i.e. the solid line, appears to have a magnitude of about 750kN per metre
of tunnel length. At the tunnel heading the normal force has not yet reached the average value of
about 750kN/m. Instead of a lower value of about 600kN/m is obtained.
Just like the normal forces the bending moments show a zigzagging pattern that matches the
‘step-by-step installation’ with d=2m. For convenience we will focus on the average value, as
indicated by the solid line. Near the tunnel heading vanishing small bending moments of about -
3kN/m are found. However with the advance of the tunnel face the bending moment in Figure 3
reaches an average steady state value of about –22kNm/m. Beyond the steady state part on the
extreme right in Figure 3 the lining is more heavily loaded up to –43kNm/m. However, this is a
numerical effect that relates to the use of smooth roller boundaries to the sides of the mesh block.
One has to note, that the two thinner lines above and below the average lines in Figure 2 and
Figure 3 mark the maximums and minimums of the zigzag results due to a more coarse 3D mesh.

Figure 2: The step-by-step installation of the tunnel leads to zigzagging normal forces in the ring
direction of the lining (compression is positive)
Figure 3: The step-by-step installation of the tunnel leads to realistic zigzagging bending moments
in the ring direction of the lining

Figure 4 and Figure 5 show the compared values of the two different mesh coarsenesses. The
finer mesh was generated with three elements per slice, as against the coarser mesh was
generated with only one element per slice. Each element contains thereby two Gaussian stress
points. It is important to note, that in this case the increase of the maximum values due to the
finer mesh is about 40% for normal forces and 10% for bending moments. This large increase is

Figure 4: Computed normal forces for coarse and fine mesh

Figure 5: Computed bending moments for coarse and fine mesh

obviously connected to the higher flexibility of the finer mesh on the one hand, which is leading to
more deflection and on the other hand to the fact, that the Gaussian stress points from the finer
mesh are positioned closer to the edges of the slice. This means, that the stress points are more
close to the actual maximums and minimums in the front and the rear part of the slice. The
obtained values of bending moments and normal forces can be used to control the design of the
tunnel lining by means of interaction diagrams (Leca et. al., 2000).
There is no doubt, that the shown zigzagging steady-state bending moments and normal
forces after some 20m of tunnel excavation are realistic, but the cost of a ‘step-by-step’ simulation
is difficult to justify for many practical tunnel applications. Therefore it is advisable to perform a 2D
analysis which matches a full 3D analysis.

3. The 2D analysis

In order to perform a 2D analysis that matches a full 3D analysis the so-called β-method (cf.
Panet and Guenot, 1982) is used, where the lining is installed after a prescribed amount of
unloading. After installation of the lining initial stresses are reduced totally and the lining is loaded.
To set the exact β-value one has to compute the respective depth of the settlement trough from a
full 3D analysis. A new method for the fast determination of this settlement value will be described
in the section ‘The fast settlement analysis’. For the present case we find a settlement of 4.5cm.
This value can then be used to select the appropriate β-factor. One may simply perform a 2D
analysis to obtain the load-settlement curve of Figure 6, where a settlement of 4.5cm corresponds
to β=0.31, i.e. 31% of the initial supporting pressure should be retained before installing the ‘2D’
tunnel lining. After lining installation this 31% is also taken away and the lining will be loaded. For
very flexible linings an additional settlement will occur, but as a rule there will hardly be any
additional settlement due to the loading of the lining.

Figure 6: 2D-tunnel analysis in which soil elements are gradually removed

Figure 7: 2D and 3D Bending moments and normal forces

After matching a 2D calculation for the present case, we find the bending moments and normal
forces as illustrated in Figure 7. This figure can be used to compare the 2D data to the zigzagging
3D data of Figure 2 and Figure 3. It appears that 2D normal forces are below the average value of
the ones from the 3D ‘step-by-step installation’. If one considers the hardening of the shotcrete for
the 3D analysis, the peak values are probably too high. To find more realistic values one has to
model the lining with a stepwise increase of stiffness, where the final stiffness is reached only
after several phases. Thus the found value of the 2D calculation is still believed to be realistic. 2D
bending moments appear to match the average ones from a 3D analysis quite well. This average
value would seem to be a highly realistic value, as the zigzagging in Figure 2 and Figure 3 is most
probably excessively large. Shotcrete shows substantial creep and stress-relaxation so that one
may expect a considerable damping of all oscillations around the average values.

4. The fast settlement analysis

Observations from practice have shown that the distribution of the developing longitudinal
settlement trough due to tunnel excavation is a s-like-curve (Mair and Taylor, 1979; cf. Peck,
1969) as shown in Figure 8. If one wants to describe the surface settlement distribution in a
certain stage of tunnel excavation this curve is given as a function depending on the location
variable x: All settlements of this curve can be described with a function S(x). As the tunnel
construction is moving on, the settlement function is also becoming dependent on the time of the
progressing tunnel excavation. With the tunnel being excavated onward the profile of the
settlement function is translated, following the moving tunnel. Therefore the settlement function is
supplemented and written as S=S(x,t). If one supposes the tunnel construction velocity as a
constant movement, the following connection between tunnel excavation length l and tunnel
construction time t can be drawn:)v=l/t and t=l/v. From that conclusion it is clear, that tunnel
construction time is proportional to tunnel excavation length: t~l. Therefore the settlement function
can be written as

(1) S = S (x , l) . (1)

After Equation (1) the settlement function is dependent on the variables x and l. Thereby the
tunnel excavation length l=l0 can be kept at an constant value, while the occurred settlement
distribution is being observed as a function dependent on the variable x, as well as the
settlements can be described in one singular point P0, while the tunnel excavation length l is the
variable and the location x=x0 of the point P0 is kept at a constant value. Therefore for the
settlement function can be written additionally S=S(x,l0)=S(x0,l). In the second case x=x0 is a

Figure 8: Settlements due to tunnel excavation

defined place on the surface and P0 is a reading point, which records the settlements of the
developing through. Based on a local coordinate system x* (note Figure 8) above the moving
tunnel face, the distance to point P0 according to a constant value of x 0 and a variable l is found
as x*=x0-l, so that for the settlement function one can write S=S(x*). The settlement function S(x*)
is dependent on many complex influencing factors and therefore describable only with great
effort. Considerably simpler to describe is the gradient of this settlement function, which can be
obtained already from one small excavation step ?l. The gradient that follows from a settlement
?s(x*) in a surface point P0 due to an excavation of ?l is described as Sl=?s/?l. The integration of
this gradient over an excavation length from -∞ up to a certain stage of tunnel construction l0,
yields the settlement in the place x 0. The final settlement S∞ is the integration of all settlements
?s(x*), which will occur in the reading point P0, when a complete tunnel over an excavation length
from -∞ until ∞ is excavated:

∞
∂s
(2) S∞ = ∫ ∂l dl .
−∞
(2)

The relationship of x*=x-l and dx*/dl=-1 leads to dl=-dx*. From this relationship one can transform
Equation (2) into

−∞ ∞
∂s ( x*) ∂s ( x*)
(3) S∞ ( x*) = − ∫ dx * = ∫ dx * . (3)
∞ ∂l −∞ ∂l

The integral of Equation (3) can then be expressed by a numerical approximation, which is
representing the ‘step-by-step installation’ method from a FEM simulation, where the soil
elements inside the tunnel are also being switched off and the tunnel lining elements are being
switched on in a stepwise procedure:

∆s ( x*) ∆x *
(4) S∞ ( x*) ≈ ∑ ∆x * = ∑ ∆s ( x*) . (4)
∆l ∆l

As shown in Figure 9 ∆l represents the excavation length of an excavation step from a FEM
simulation and ∆s(x*) is the occurring settlement due to this excavation step. ∆x* is the distance
of the repeating interval where ∆s(x*) is being evaluated. For practicable application of Equation
(4) one chooses ∆x*=∆l and obtains

(5) S∞ ( x*) ≈ ∑ ∆S ( x*) . (5)

Figure 9: Development of the longitudinal settlement trough in an FEM analysis

The repeated summation of the settlements ∆s(x*) from one single excavation step ∆l responds to
a translation of the settlement curve. Figure 9 indicates clearly how the repeated summation of a
steady-state settlement crater in a 3D FEM calculation leads to the final distribution of the
longitudinal settlement curve. This summation can also be done with the whole shape of the
settlement crater in order to achieve the entire 3D settlement trough.
If one considers a vertical section through the array of curves and a surface point P0 in this
section, then one will find, that the distances between the settlement curves ∆s i are responding to
the settlements ∆s(x*), which this surface point experiences due to the proceeding tunnel
excavation. The sum of all these settlements results in the final settlement S8 . Assuming that the
settlement craters in consequence of equal excavation lengths in the steady-state region of the
longitudinal settlement curve remain the same, the conclusion can be made, that the settlements
∆s i from the vertical section respond to the settlement ∆s(x*) of the settlement crater from one
single excavation step: The same amount of settlements from the vertical section can be found
within the distances of ∆l on the settlement curve from one excavation step. Therefore the final
settlement S8 can already be obtained from one single settlement crater.

5. The ‘all-in-once installation’

The use of this theory in a FEM simulation leads to the conclusion, that one simply needs to
evaluate only one excavation step instead of many excavations in a complex ‘step-by-step
installation’ method. From the obtained settlement crater of one excavation step one can choose
a longitudinal section over the tunnel axis or in any other desired position of this crater. The
described summation will then lead to the final settlement in this section. For the present case this
method led to a final settlement over the tunnel axis of s=4,3cm, whereas from the full 3D
calculation a value of s=4,5cm was obtained, i.e. including an error of only 5%. For the application
of this method it has to be noted, that the computed settlement crater due to one single
excavation step needs to have a steady-state value (Vermeer et. al., 2002). This can only be
accomplished, if the so called boundary effects from the FEM model are eliminated. In order to
eliminate these boundary effects one needs to install the whole tunnel lining up to the steady-
state region of the 3D settlement trough in one calculation phase with an unsupported excavation
length of ∆l in front of the tunnel lining. Afterwards one performs a second phase, where one
more slice of soil element with the length of ∆l is switched off and the previous unsupported slice
gets supported by the tunnel lining. All displacements from the first phase are reset to zero and
the computed values will only respond to the settlement crater of the second phase, as shown in
Figure 10. This method is referred to as the ‘all-in-once installation’.
Finally it should be noted that a 3D analysis of volume loss does not require the relatively fine
mesh of Figure 1. On computing the settlement crater after two computational phases as
described above, the fine mesh is only needed at the tunnel face in the middle of the mesh block.
Towards the back and beyond the tunnel face the element size can be gradually increased.

Figure 10: Settlement crater from ‘all-in-once installation’

6. Conclusions

It has been shown that a full 3D analysis of tunnel excavation is needed for the prediction of loads
on linings, as the values of both bending moments and normal forces are found to have a zigzag
pattern, which can not be found in a 2D calculation. Therefore an extremely large number of
excavation steps is needed. Moreover it has been shown, that a relatively fine 3D mesh is
required, because values of both bending moments and normal forces are found to be much
higher than from the solution of a relative coarse mesh. Due to this requirement computer run
time gets excessive. Such full 3D analysis are not feasible in engineering practice.
On comparing 2D and 3D analyses one observes that a 2D analysis matches the values from
a full 3D analysis well. If one considers the hardening of the shotcrete and models the 3D
analysis with a stepwise increasing stiffness, then even more realistic values, as also obtained
from the 2D analysis, are believed to be found. This will be a topic of further research.
Such 2D analysis can respond to a full 3D analysis when the settlement values are taken the
same. It has been shown that the settlement value, which responds to a full 3D analysis, can be
obtained from a 3D analysis, that requires little computational effort. This method is based on
reaching a steady-state value after one excavation step. It is concluded that this method is
applicable for constitutive laws with both liner elastic and constant plastic soil behaviour. However
if there is non linear soil behaviour and also other boundary conditions, e.g. very shallow tunnels,
it is believed that a steady-state value will be reached only after several steps of excavation. The
determination of those influences will be another topic of further research. In order to arrive at the
same settlement in a 2D analysis one has to select appropriate β-values. It has been shown that
β-values can be obtained from a very simple 2D analysis.

7. References

Leca, E., Leblais, Y and Kuhnhenn, K. (2000). Underground works in soils and soft rock tunneling. Int.
Conf. On Geotech. and Geolog. Eng., Vol. 1 pp 220-268, Melbourne
Mair, R.J. and Taylor, R.M. (1997). Bored tunnelling in the urban environment. Proceedings of the 14th Int.
Conf. on Soil Mech. and Found. Eng., Vol. 4 pp 2353-2385, Hamburg
Panet, M. and Guenot, A. (1982). Analysis of convergence behind the face of tunnel. Institution of Mining
and Metallurgy, Proc. Tunneling 82 pp. 197-204, London
th
Peck, R.B. (1969). Deep excavation and tunnelling in soft ground. Proc. 7 International Conference on Soil
Mechanics and Foundation Engineering, State of the Art Volume. pp. 225-290, Mexico City.
Vermeer, P.A. (2001). On a smart use or 3D-FEM in tunnelling. Bulletin of the PLAXIS Users Assciation
(NL), N0 11, Delft
Vermeer, P.A. and Ruse, N. (2001). Tunnel Face Stability (in German). Geotechnik, 24 No 3 pp 186-193
Vermeer, P.A., Bonnier, P.G. and Möller, S.C. (2002). On a smart use of 3D-FEM in tunnelling. 8th Int.
Symp. Num. Mod. Geomech., Rome
Wittke, W. (1984). Rock Mechanics. Springer-Verlag, Berlin