On the Application of Numerical Analysis in Tunnelling
P. A. Vermeer, S. C. Möller & N. Ruse
Institute of Geotechnical Engineering, Pfaffenwaldring 35, 70569 Stuttgart,Germany
At present tunnels tend to be analyzed on the basis of 2D finite element computations , because 3D analyses are considered to be ex-
tremely time consuming. As a result, 3D analyses are presently the domain of researchers. Consulting engineers will only perform 3D-
FEM analyses when facing complex geometries, e.g. tunnel joints or connections to underground stations, but not for straight-ahead
tunneling. For judging possibilities, we distinguish between the 3 main focuses of tunnel analyses: tunnel heading stability, surface set-
tlements and structural forces in linings. No doubt, tunnel heading stability requires a 3D FE-analysis, as otherwise one can not possi-
bly capture the very significant arching in frictional ground. Here, it will be shown that safety factors can be computed relatively easily
by the so-called SSR-FEM. In contrast to stability analyses, settlement analyses tend to require very large 3D meshes, but a so-called
smart 3D analysis is proposed to overcome this difficulty. For the computation of structural forces in tunnel linings, a relatively simple
2D analysis is presented which matches the solution of a full 3D analysis reasonably well.
1 INTRODUCTION 2 TUNNEL HEADING STABILITY
The uses of underground space are many, varying from simple One of the main problems when constructing a tunnel is to ensure
storage caves and cellars to large complexes for storage of both the stability of the tunnel heading. In many European countries
liquids and solid materials and finally to many different kinds of this stability is analysed on the basis of the sliding wedge mecha-
underground traffic and transport routes. Tunnels built for high- nism in Fig. 2, being used for shield tunnelling as well as
ways and railroads may be shallow in urban areas or deep under- NATM-tunnelling. No doubt, this method gives an indication
neath major mountain ranges. Such deep tunnels may operate un- about the minimum support pressure needed in shield tunnelling
der considerable external groundwater pressure, whereas water and the factor of safety in NATM-tunnelling, but the model is
supply and sewage tunnels also operate under internal pressure. crude and thus bound to be replaced by finite element analyses.
Both the cut & cover method for very shallow excavations and Shield tunnelling in homogeneous soil can also be analysed by
deep underground excavations call upon the disciplines of soil means of empirical formulas for the minimum failure pressure pf.
and rock engineering. Such an equation is of the form
As for all construction activities empirical knowledge and engi-
neering judgement play a dominant role, but with the advance of p f = c ⋅ N c + γ ⋅ D ⋅ Nγ + q ⋅ N q (1)
computers and geomechanics numerical analysis is of growing
importance. Indeed, forces in support structures of complex where c can be a drained as well as an undrained cohesion, γ is
three-dimensional underground openings (Fig. 1) need to be as- the unit soil weight, D the tunnel diameter and q a possible sur-
sessed numerically. Moreover numerical analyses are used for face load. For undrained situations, the stability numbers Nc, Nγ
settlement predictions. On the other hand possibilities of finite and Nq depend on the relative cover C/D. For drained conditions,
element analyses for tunnel heading stability are rarely used. In these stability numbers depend on the soil friction angle.
this paper attention is focussed on stability as well as on settle-
ments and structural forces.
Fig. 1. Complex three dimensional tunnel geometries Fig. 2. Analytical wedge model for tunnel heading stability analysis
3 ON THE STABILITY NUMBERS
The face pressure Eq. (1) was intensively studied by Vermeer et
al. (2002). Despite its restriction to homogeneous soil conditions,
the formula is of great practical value as it offers good insight
into the face stability problem, at least in combination with
proper expressions for Nc, Nγ and Nq. For undrained shield
tunnelling, we proposed to use
N c = 5.86 ç ÷ and N γ = 0.5 + C D . (2)
For drained conditions, we derived
Nc = cot ϕ´ , Nγ = − 0.05 , Nq = 0 . (3)
9 tan ϕ´
At least for ϕ´ > 20° and a relative cover C/D > 1.5.
The face-pressure Eq. (1) not only of great importance for getting
insight into the problems of face stability, but also for validation
of numerical analyses. Indeed, the non-linear elastoplastic finite
element method is a powerful tool, but computer codes need to
be validated. Moreover, users will have to train themselves and
need benchmark problems with known solutions when using rela-
Fig. 3. Above: Typical strength-displacement curve for assessment of
tively new numerical methods such as the SSR-FEM.
safety factor by SSR-FEM. Below: Shadings of incremental displace-
ments at collapse, showing shell type sliding body.
4 SHEAR STRENGTH REDUCTION FEM (SSR-FEM)
The most versatile method for stability analysis is the SSR-FEM.
This method has been applied to slope stability analysis in two-
dimensional as well as three-dimensional situations. Numerical
comparisons have shown that SSR-FEM is a reliable and robust
method for assessing the safety factor and corresponding failure
mechanism. One of the main advantages is that both the safety
factor and the mechanism emerge naturally from the analysis. In
SSR-FEM, the factor of safety is defined as
tan ϕ real c
F = = real (4)
tan ϕ min cmin
where cmin and tan ϕmin are minimum values as needed for equi-
librium. These values are obtained by reducing the real shear Fig. 4. Soil profile of Rennsteig road tunnel in Thuringia, Germany
strength parameters stepwise down to failure in an elastoplastic
FE-analysis. The method originates from Zienkiewicz et al Relevant data for a stability analysis is also indicated in Fig 4.
(1973) and has been used most recently by Cai & Ugai (2003) for Using these data we assessed factors of safety for several soil pro-
anchored slopes. We have adopted this method for three- files; two resulting factors for the top-heading being shown in
dimensional stability analyses of tunnel headings (Vermeer et al, Fig. 4 and denoted by the symbol F. It should be noted that all
2002). Fig. 3 shows results from a particular tunnel stability SSR-FEM results were obtained by using a fully three-
problem. Considering homogeneous ground we have compared dimensional code (Plaxis 3D Tunnel). However, numerical effort
this method to the face-pressure Eq. (1) to find excellent agree- is limited as relatively small meshes can be used as indicated in
ment. Fig. 3. This is completely different from 3D settlement analyses,
where relatively large meshes are required in order to guarantee
5 CASE STUDY WITH SSR-FEM
The SSR-FEM is applicable both in shield and NATM tunnel- 6 STEP-BY-STEP EXCAVATION
ling. For NATM tunnelling, this is well demonstrated by the
Rennsteig tunnel in Thuringia. With a length of 8 km it is the To get a better understanding of mesh requirements for a 3D set-
longest motorway tunnel in Germany. The excavation of this tlement analysis we modelled a NATM tunnel. The sequential
double-tube tunnel was done by sequential construction of a top excavations of top-heading and invert were analysed with the
heading followed by bench and invert, as indicated in Fig. 4. elastoplastic Mohr-Coulomb (MC) model and results are shown
Fig. 5. Computed displacements for an NATM tunnel. The 1st figure
shows displacements after excavation of top-heading in graded shadings
from black to white. The 2nd figure shows the corresponding settlement
trough. The 3rd figure shows shadings of displacements after excavation
of the invert.
in Fig. 5. The way of modelling NATM tunnels is also known as Fig. 6. Settlement trough and ground-response curve from classical 2D
the Step-By-Step Method where settlements are obtained by a analysis with unloading factor β
stepwise removal of soil elements and installation of tunnel lin-
ing. From this analysis we obtained that meshes have to be larger Usually in practice unloading factors are either chosen according
than the steady-state settlement trough itself. Due to boundary to the experience of the judging engineer or by considering two
conditions of vertical rollers initial disturbances on the right side different conservative cases: e.g. β = 0.2 for settlements, as it will
of Fig. 5 affect settlements over a considerable excavation length. give relatively large settlements and β = 0.7 for structural forces
As a consequence one has to simulate a lot of excavation steps in in the lining, as it will give relatively large loads on the lining.
order to arrive at a reliable steady state solution. Therefore such a No doubt this method is not very satisfactory from a scientific
three dimensional settlement analysis is very time consuming point of view.
(Vermeer et al, 2002).
8 SMART 3D-ANALYSIS
To overcome shortcomings of 2D analyses as well as time con-
Because of the costs of a full 3D calculation, most tunnels are suming 3D analyses we developed a smart way of 3D analysis.
presently analysed by the use of 2D finite element computations. Based on this new approach we use a full three dimensional
In practice the most widespread 2D method is the so-called Load model but instead of excavating many slices we only compute
Reduction Method (Panet & Guenot, 1982). As indicated in Fig. two excavations. From the last excavation increment we can ex-
6, a prescribed percentage β of initial stress σ0 is left inside the trapolate the entire three dimensional settlement trough. This
tunnel as a support pressure to take into account the missing three method has been well described by Möller et al (2003). The final
dimensional arching effect. Most settlements (∆s1) will take place settlement value coming from this smart 3D calculation can now
in this phase. After installation of the lining this support pressure be used to find appropriate unloading factors for a 2D calcula-
is removed and some additional settlements (∆s2) occur due to tion. For a tunnel with sequential excavations we obtained
loading of the lining. Indeed, this way of modelling the trans- β = 0.3 for top heading and β = 0.7 for invert (Bonnier et al,
verse settlement trough is very efficient, as very low computer 2002) giving insight that unloading factors are also strongly de-
run time is required. Moreover, by comparing results we obtained pendent on the geometry of excavated tunnel sections.
that for appropriate unloading factors β the transverse settlement
trough exactly corresponds to the 3D solution. Although the
Load Reduction Method is straight forward to use and therefore 9 STEEPNESS OF SETTLEMENT TROUGH
very common in practice there is still one major shortcoming re-
maining about it. No secured method has yet been provided for Besides determining unloading factors for a two dimensional
determining proper β-values. Settlements and also loads on lin- analysis the fast way of tunnel analysis is a useful tool to analyse
ings will be computed directly dependent on this input parameter. the settlement trough using different soil models. To compare
results of different soil models with measurements from site we
modelled the Steinhaldenfeld subway tunnel in Stuttgart, Ger-
many shown in Fig. 7. The ground shown in Fig. 8 consisted
mainly of overconsolidated marl, which may be considered to be
a hard soil as well as a soft rock. As the 1000m long NATM tun-
nel was constructed in an urban area settlements were carefully
For modelling the settlement trough in practice, the Mohr-
Coulomb (MC) Model is a very common approach as models of
higher order are mostly the domain of researchers. Fig. 9 shows
results of MC and Hardening Soil (HS) Model (Brinkgreve et al,
2002) compared to measurements from site. It is observed that
the MC-model predicts a relatively wide and shallow trough,
which deviates significantly from the measured one. However HS
Fig. 7. Steinhaldenfeld tunnel in Stuttgart does a lot better as it takes into account the different behaviours
for primary loading and unloading-reloading of overconsolidated
10 BENDING MOMENTS AND NORMAL FORCES
An important topic of tunnelling concerns the assessment of
structural forces in linings. In order to predict realistic values of
bending moments and normal forces a full 3D analysis is needed.
To calculate realistic values for a circular tunnel, we divided a
block of 100x40x28m into 8840 volume elements with a total of
26809 nodes (Fig. 11). For the parameters of the MC-model, we
Fig. 8. Ground profile of Steinhaldenfeld tunnel
assumed 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 dila-
tancy angle of ψ=0 and K0=1-sinϕ for the coefficient of initial
lateral earth pressure. The NATM tunnel with a diameter of 8m
and a cover of 16m was modelled in a symmetric half with an un-
supported excavation length of 2m. The shotcrete lining has a
thickness of 30cm, a Young’s Modulus of 20MPa and a Pois-
son’s Ratio of ν=0.
Fig. 12 presents resulting normal forces after 80m of ‘step-by-
step installation’. Within a single lining segment of d=2m a sharp
drop of tangent normal force is observed from 1400kN/m at the
front to 200kN/m at the back. This is due to the fact that the un-
supported tunnel heading is mostly arching on the front and not
so much on the back of the last lining segment. The average tan-
gent normal force, i.e. the solid line, appears to have a steady-
state magnitude of 750kN per meter of tunnel length. The steady
Fig. 9. Cross section of observed and computed settlement trough of state is characterized by the horizontal part of the solid line The
Steinhaldenfeld tunnel average normal force following directly after the tunnel heading
is 600kN/m and still below the average value of 750kN/m. Fur-
ther excavation steps will disturb the 3D arching effect around
the tunnel heading and therefore lead to additional loading on the
lining. In the steady state the arching is predominantly 2D.
Fig. 11. Shadings of vertical displacements after 80m of stepwise excava-
Fig. 10. Observed and computed longitudinal settlement trough of Stein- tion. The steady-state settlement of s=4.5cm is reached after an excava-
haldenfeld tunnel tion of 35m.
Similar to normal forces bending moments follow a zigzagging
pattern as shown in Fig. 13. Considering average values, small
bending moments of -3kN/m are found next to the tunnel head-
ing. In the middle of the excavated tunnel bending moments
reach an average steady state value of -22kNm/m. Towards the
right boundary of the finite element mesh the lining is loaded
more heavily giving a bending moment of –43kNm/m. However,
this is a numerical effect that relates to smooth roller boundaries.
One has to note, that the two thinner lines above and below the
average lines in Fig. 12 and Fig. 13 mark the maximums and
minimums of results due to a coarser 3D mesh. Indeed, computed
Fig. 12: The step-by-step installation of the tunnel leads to zigzagging structural forces appear to depend significantly on the discretisa-
normal forces in the ring direction of the lining (compression is positive) tion of lining elements. Fig. 14 and Fig. 15 show values of the
two different mesh coarsenesses in more detail. The finer mesh
was generated with three elements per lining segment and the
coarser mesh was generated with only one element per segment.
Each element thereby contains two Gaussian stress points. It is
important to note, that the increase of the maximum values due to
the finer mesh is about 40% for normal forces and 10% for bend-
ing moments. This large increase relates possibly to the higher
flexibility of the finer mesh and obviously very much to the fact,
that the Gaussian stress points from the finer mesh are positioned
closer to the edges of the segment. This means, that the stress
points are more close to the actual maximums and minimums in
the front and the rear part of the segment.
Fig. 13: The step-by-step installation of the tunnel leads to realistic zig-
There is no doubt, that the zigzagging bending moments and
zagging bending moments in the ring direction of the lining
normal forces are realistic, but the cost of a ‘step-by-step’ simula-
tion is difficult to justify for many practical tunnel applications.
Therefore it is advisable to perform a 2D analysis in combination
with a smart 3D Analysis as discussed in previous sections. After
matching a 2D calculation to the 3D solution giving β=0.3 for the
present case, we computed bending moments and normal forces
as illustrated in Fig. 16. This figure can be used to compare the
2D data to the zigzagging 3D data of Fig. 12 and Fig. 13. It ap-
pears from Fig. 16a, that 2D normal forces are below the average
value of the ones from the 3D ‘step-by-step installation’. On the
other hand, Fig. 16a shows that 2D bending moments appear to
match the average ones from a 3D analysis quite well. This aver-
age value would seem to be a highly realistic value, as the zig-
zagging in Fig. 12 and Fig. 13 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.
Fig. 14: Computed normal forces for coarse and fine mesh
Fig. 15: Computed bending moments for coarse and fine mesh Fig. 16: 2D and 3D Bending moments and normal forces
Examples of numerical applications for tunnel heading stability
and settlements have been shown. For analysing tunnel heading
stability, the shear strength reduction method appeared to be a
powerful tool to obtain both factors of safety and mechanisms of
failure. However, as numerical results need to be validated, it is
recommended to compare factors of safety from shear strength
recuction method as much as possible to findings from the face
pressure Eq. (1) with stability numbers as given by Eps. (2) and
For computing settlements it has been shown, that transverse set-
tlement troughs can be modelled well with the load reduction
method, at least when using an appropriate constitutive soil
model. However appropriate values for unloading factors (β) re-
main uncertain. Therefore reliable settlements may only be mod-
elled three dimensional. A proposed smart 3D analysis can over-
come shortcomings of time consuming conventional 3D analyses.
Results of smart 3D analysis have been shown for different soil
models. Mohr-Coulumb Model gives too wide and too shallow
settlement troughs, but the Hardening Soil Model can improve a
lot, at least for overconsolidated soils.
Results for structural forces in tunnel linings from ‘step-by-step’
simulation have been shown. Both bending moments and normal
forces appear to have a zigzagging pattern within one lining seg-
ment. It has been shown, that relatively fine FE meshes have to
be used in order to compute exact maximums and minimums of
the observed zigzagging. Even though these results are believed
to be highly realistic the expense of such a full 3D analysis is dif-
ficult to justify at least for straight ahead tunnelling. Therefore
for the computation of structural forces in tunnel linings it is ad-
vocated to use a 2D analysis. Appropriate unloading factors may
come from the proposed “Smart 3D Analysis”.
Bonnier, P.G., Möller, S.C. & Vermeer, P.A. 2002. Bending
Moments and Normal Forces in Tunnel linings. Paper pre-
sented at the 5th European Conference of Numerical Methods
in Geotechnical Engineering, Paris, France.
Brinkgreve, R.B.J. & Vermeer, P.A. 2002. Plaxis 3D Manual.
Delft: A.A. Balkema.
Cai, F. & Ugai, K. 2003. Reinforcing mechanism of anchors in
slopes: a numerical comparison of results of LEM and FEM.
Numerical and Analytical Methods in Geomechanics (27):
Möller, S.C., Vermeer, P.A. & Bonnier, P.G. 2003. A fast 3D
tunnel analysis. Paper presented at the Second MIT Confer-
ence on Computational Fluid and Solid Mechanics, Boston,
Panet, M. & Guenot, A. 1982. Analysis of convergence behind
the face of tunnel. Institution of Mining and Metallurgy, Proc.
Vermeer, P.A., Bonnier, P.G. & Möller, S.C. 2002. On A Smart
Use of 3D-FEM in Tunneling. Paper presented at the 8th In-
ternational Symposium on Numerical Models in Geomechan-
ics, Rom, Italy.
Vermeer, P.A., Ruse, N. & Marcher, T. 2002. Tunnel heading
stability in drained ground. Felsbau: 8 – 18.
Zienkiewicz, O.C., Humpheson, C. & Lewis, R.W. 1975.
Assosciated and nonassociated visco-plasticity in soil mechan-
ics. Géotechnique (25) : 671 - 689.