Available online at www.CivileJournal.

org

Civil Engineering Journal

Vol. 2, No. 6, June, 2016

Application of Hyperstatic Reaction Method for Designing of Tunnel

Permanent Lining, Part I: 2D Numerical Modelling
Rahim Hassani a, Rouhollah Basirat b*
a
MSc in Rock Mechanics Engineering, Amirkabir University of Technology, Tehran, Iran
b
PhD Student, Rock Mechanics Engineering, Tarbiat Modares University, Tehran, Iran

Received 14 May 2016; Accepted 21 June 2016

Abstract
The increase of bored tunnels in the entire world has raised the question how to design the tunnel structure in an
efficient way. This paper proposes a numerical approach to the Hyperstatic Reaction Method (HRM) for analysing
permanent tunnel linings. The permanent tunnel lining is known as main structure of tunnel maintenance during the time.
The HRM is one of the analysis methods for tunnel lining in long term. In this paper, two dimensional numerical
modelling is performed by considering hyperstatic reaction concepts. Loading is done after the calculation of long term
loads, and ground reaction is simulated by springs. Designing is done for Manjil-Rudabar freeway project, Tunnel No. 2.
The numerical analyses were performed for Operational Design Earthquake (ODE) and Maximum Design Earthquake
(MDE) loading conditions. A new simplified approach is used for considering the effect of earthquake loading on the
tunnel lining. Then, an interaction diagram between axial force and bending moment used for investigating the capacity
of tunnel lining. The thickness of tunnel lining and armature are calculated for three sections based on induced forces in
tunnel lining. These forces were different in every section according to the load combinations, rock mechanics properties,
lining properties, and overburden. The numerical results showed that the forces in tunnel lining for MDE condition is
approximately 50% more than ODE condition in earthquake loading. This numerical processing presented that the HRM
is a proper, fast, and practical method for designing and analysing the tunnel lining.
Keywords: Hyperstatic Method; Tunnel Lining; Numerical Modeling; Static Analysis; Dynamic Analysis.

1. Introduction
Numerical tunnel analysis is generally conducted for preliminary design with the two-dimensional beam-spring
model, which consists of linear beam and spring elements to simulate the behaviour of lining and ground, respectively.
Although this model is not only old but there are also even several advanced numerical models, lots of design
experiences in the old model have been accumulated. Furthermore, model complexities are not the most important key
factor for tunnel analysis. When factors such as input parameters, boundary conditions, and ground loads are well-
estimated, a simple model can give a good prediction which is comparable to or more accurate than results from
complex models.
These facts enable the beam-spring model to usefully and approximately examine the state of the tunnel, and
consequently the model has been continuously used. These limitations of the model could be improved by introducing
the Winkler-based beam element derived from the beam on foundation problem. This element considers ground
resistance to be distributed all over the element length, unlike the spring element. It was reported that the use of the
Winkler-beam foundation element enhances the convergence rate of the ground-structure interaction problem [1].
There has also been research related to expansions and refinements of the model. These studies considered the
nonlinear behaviour of the lining or ground in the beam-spring model [2-4].
The support ground interaction influences the stress state in the structure and this interaction depends on the

* Corresponding author: r.basirat@modares.ac.ir

244
Civil Engineering Journal Vol. 2, No. 6, June, 2016

mechanical characteristics of the ground. As these are only generally known with a certain approximation, it is often
necessary to carry out parametric or probabilistic type analyses in order to be able to completely describe the
uncertainty on the stress state of the support structure. These types of analyses need many calculations and the
hyperstatic reaction method (HRM) results to be particularly suitable for this purpose; due to the short time it requires
[4].
Figure 1 shows the beam-spring model with the following parameters: stress σ, strain ε, Young’s modulus E, spring
force F, spring deformation δ, and spring coefficient k. Note that the compression indicates positive values in the
material constitutive relation.
This type of model is preferred due to its simplicity and the ability to reasonably simulate structure–ground
interaction.
When the ground is in tension, it loses load resistance capacity. This phenomenon is reflected by using compression
only spring elements; a truss element may substitute for the spring element. The use of these elements leads to an
iterative process. The process terminates when there are no springs in tension.

Figure 1. Typical beam-spring model

Orest (2005) presented a probabilistic numerical approach for the design of primary tunnel supports, according to
the HRM. This in turn allowed an estimation of the costs of the different support systems to be made in relation to the
reliability level of the economic evaluations [5].
Do et al (2014) developed a specific implementation using a FEM framework for segmental tunnel lining. The
numerical results presented in the paper showed that the proposed HRM can be used to effectively estimate the
behaviour of a segmental tunnel lining [6].
In this paper, the applications of HRM are presented for designing of permanent tunnel lining including two
dimensional (2D) analyses. A real case study (Manjil-Rudabar freeway project, Tunnel No. 2) is designed by this
method.

2. Winkler-Based Beam Element

The Winkler-based beam element was derived from the governing differential equation for lateral deformation of a
beam. The original differential equation should be expanded to consider axial deformations due to tunnel geometry. In
the weak form, which can be interpreted as the well-known principle of virtual displacements, the differential equation
is:

∫ ( ) ( ) (1)

With the parameters as follows: x: local coordinate in element domain, v: lateral deformation, u: axial deformation, E l:
elastic modulus of the lining, Il: area moment of inertia of the lining, Al: cross-section area of the lining, η: subgrade
reaction modulus, q: arbitrary lateral load intensity, and p: arbitrary axial load intensity. With integrating by parts and
discretization, Eq. (1) was converted to the equilibrium equation of the Winkler-based beam element written in the
following form:
(2)
where Ks is the lining stiffness matrix, ; Kg is the ground stiffness matrix,
; d is the element nodal displacement vector; Ke is the total element stiffness matrix, K s+Kg; Fext is
the total element force vector, ; Bf is the flexural deformation displacement relation matrix,
; Ba is the axial deformation-displacement relation, ; Nf is the flexural shape function matrix using a
cubic Hermitian function; Na is the axial shape function matrix using a linear Hermitian function; and P is the nodal
force vector[7].

245
Civil Engineering Journal Vol. 2, No. 6, June, 2016

3. The Case Study (Rudbar-Manjil Freeway, Tunnel No. 2)

The Manjil-Rudabar freeway is one of the civil projects under construction in Iran. To complement the Qazvin –
Manjil freeway, the construction of two twin tunnels has been predicted. In this paper, the tunnel permanent lining is
designed for Tunnel No. 2. The length of the right and left tunnels are more than 1000 meters. The width of the right
and left tunnels are 14 and 12 m, respectively.
Based on geological longitudinal profile of Tunnel No. 2, there are some lithology units at the tunnel elevation
including pyroclastic andesitic rocks with tuff and tuff-breccia faces. In this paper, three sections are selected for
analysing. These sections are the weakest rocks with different overburdens. Figure 2 shows the geology longitudinal
profile map for the case study. The rock mechanics properties of three sections for numerical analysis are illustrated in
Table 1.

Figure 2. The geology longitude profile map for Manjil-Rudbar Freeway, Tunnel No.2

Table 1. Rock mechanics properties of tunnel sections

Parameter Unit Eat Edtb FZ
Elastic Modulus GPa 1.2 0.5 0.27
Density kN/m3 25 26 25
Cohesion kPa 450 180 200
Friction Angle Degree 35 30 20
Passion’s ratio - 028 0.3 0.35
Overburden m 40 60 140

4. Tunnel Lining Designing and 2D Numerical Simulation

4.1. Estimation of Spring Stiffness
The structural elements (i.e., beams) are usually modelled as linear elastic; their stiffness is a function of the
thickness and the elastic modulus of the constituting materials. Since the tunnel first lining is made of shotcrete and
steel ribs it is necessary to define an equivalent tunnel cross section and a modulus of deformability which take into
due account the different properties of shotcrete (continuous) and steel ribs (discontinuous).
The stiffness of the springs Kr and Kt (Figure 3) are usually evaluated from the rock mass data (Table 1) using very
simple relationships as those derived from Winkler theory [8-10]. The interface between lining and rock cannot
withstand tension; therefore, interface elements may be used or the springs deactivated when tensile stresses occur.
The radial and tangential spring stiffnesses, expressed in units of (subgrade reaction
coefficient), are estimated from [8]:

(3)

(4)

Where: and are radial and tangential spring stiffness, respectively.  is arc subtended by the beam element
(radian) and b is length of tunnel element considered.

Figure 3. Details of the rock–support interaction through Winkler’s spring’s criterion in the HRM

246
Civil Engineering Journal Vol. 2, No. 6, June, 2016

Since there are two θ angle in the tunnel sections (double-arch tunnel section), two spring stiffness are calculated. The
spring stiffness for different sections is illustrated in Table 2.

Table 2. Spring stiffness of tunnel sections

Unit θ1 θ2
Kr Kt Kr Kt
Eat 5026 2011 3351 1340
Edtb 2094 838 1396 559
FZ 1131 452.5 754 301.6

4.2. Loads and Load Combinations

Loads considered in the design of the lining are categorized according to their frequency of magnitude, continuity
and variation. Vertical and horizontal earth pressure, water pressure, dead weight of the lining, effects of surcharge
and other factors are fundamental ones, which continuously act on the lining without large variation and should be
always considered in the design of lining. In the design of the lining, it is necessary to select the substantial loads out
of the loads described above and to decide the appropriate magnitude of design load for each selected one [11].

4.2.1. Loads and Boundary Conditions

a) Vertical and horizontal soil/rock pressures
Methods giving loads exerted by the ground on the support determine the extent of the failure zone. Purely static
considerations then determine the reaction that needs to be exerted by the support to keep the failure zone stable.
These methods implicitly assume that severe convergence has occurred for failure mechanisms to occur; the
corresponding displacements are not necessarily acceptable for the support structure. Methods differ in the way they
define the failure zone.
However, what is known as the HRM deserves special mention. The support is modelled as bars and the ground
reaction. It is an attempt to address the interaction between the ground and support. The load on the support comes
from the actions needed to maintain the failure zone in equilibrium and the reaction of the ground to yielding of the
support [12].
The method requires the definition of the active loads that apply directly to the support structure. These loads can be
estimated using different methods that are known in the scientific literature [13].
Other passive loads, which can be developed and act on the tunnel lining sections at which the tunnel lining moves
towards the ground surrounding the tunnel, are due to the reaction of the ground to the displacement of the tunnel
lining.
Figure 4 shows the vertical and horizontal pressure (rock load) based on failure zone with boundary conditions in
HRM. These loads calculated by different approaches and are presented in Table 3.

Table 3. The vertical and horizontal pressure (rock load)

Unit Vertical Load (ton) Horizontal Load (ton)
Eat 15 8
Edtb 24 12
FZ 44 22

Figure 4. The vertical and horizontal pressure (rock load) based on failure zone with boundary conditions in HRM

b) Earthquake loading
Varies researches show that earthquake load is affected on the shallow tunnels [14]. Seismic safety with
consideration of the intended purposes of the tunnel shall be appropriately examined considering the importance of the
structure, the magnitude of the earthquake ground motion, geographical and ground conditions, structure and shapes,
and so on.
A new simple method is used for considering earthquake load in this paper. Figure 5 shows the loading condition for

247
Civil Engineering Journal Vol. 2, No. 6, June, 2016

earthquake loading. Δdlining is induced displacement by earthquake loading which can be calculated by analytical
methods that are presented in reference [15]. Considering properties in Table 1, Δdlining for the maximum and
operational design earthquakes (ODE and MDE) conditions is illustrated in Table 4.

Figure 5. Loading condition for earthquake loading

Table 4. Induced displacement by earthquake load

Unit ODE (mm) MDE (mm)
Eat 9.3 17.4
Edtb 14.4 26.8
FZ 12.0 22.4

c) Dead Load
Concrete load is considered 2.5 m3 for every one cubic meter of concrete.

4.2.2. Load Combinations

The ultimate limit state is generally verified based on every limit state for the various combinations of loads. On
the other hand, a combination of loads is selected for the limit state for cracking, deformation and other factors in the
check of serviceability limit state. As the characteristic values are decided for every combination of loads, the load
combination factor is unnecessary [11].
Design loading criteria for underground structures has to incorporate the additional loading imposed by ground
shaking and deformation. Once the ground motion parameters for MDE and ODE have been determined, load criteria
are developed for the underground structure using the load factor design method. This section presents the seismic
design loading criteria for MDE and ODE. Given the performance goals of the MDE the recommended seismic
loading combinations using the load factor design method for cut-and-cover tunnel structures is as follows [14]:
(5)
where U, D, EX and EQ are required structural strength capacity, effects due to dead loads, effects due to excavation
loads, and effects due to design earthquake motion, respectively.
For the ODE the seismic design loading combination depends on the performance requirements of the structural
members. The following loading criteria are recommended:
(6)

5. Model Verification
In the first step the numerical method is verified by an analytical method suggested by Japan Society of Civil
Engineering (JSCE). A common load distribution model for this method is shown in Figure 6, where vertical soil
reaction is uniform and horizontal soil reaction is distributed in a triangular between 45 to 135 degree from the crown
on both side. Horizontal deformation of a ring at the spring line, which will determine the magnitude of horizontal soil
reaction, is different, depending on whether the soil reaction derived from the dead weight of the lining is considered
or not. The axial force and bending moment is calculated according to Table 5.

248
Civil Engineering Journal Vol. 2, No. 6, June, 2016

Figure 6. The distribution of loads used in JSCE [11]

Table 5. Equations for member forces in JSCE [11]

Load Bending moment Axial force

Vertical load

Horizontal load

Horizontal
Triangular load

Soil reaction
( )

( ) ( )

Dead weight
( )
( (

) )

Without considering soil reaction derived from dead weight of lining:

( )
Horizontal
deformation of a
Considering soil reaction derived from dead weight of lining:
ring at spring
line ( )

Where EI: Flexural rigidity in unit width

A circle tunnel is considered as a verification model. The radius and thickness of tunnel lining are 3 and 0.5 m,
respectively. The horizontal and vertical loads are considered 15 and 8 tons.
The axial force and bending moment in the tunnel lining are shown in Figure 7 for analytical and numerical
methods . Figure 7 shows a meaningful agreement in numerical simulation with analytical method.
The maximum of difference in forces is less than 10%.

249
Civil Engineering Journal Vol. 2, No. 6, June, 2016

120 -15
100 -10
80 -5
60 0
40 5
20 10
θ 90
0 15

Bending Moment
Axial Force
(ton)

(ton.m)
Analytical
Analytical
Numerical
Numerical

Figure 7. The axial force and bending moment in the tunnel lining

6. Numerical Results
Figure 8 shows the numerical model in section Eat. Sixty seven elements are used in this section. The displacement
in three directions is fixed and the rotation is free.

Figure 8. The elements and joints in the numerical model

6.1. State of Axial Force and Bending Moment

The 2D numerical simulation of tunnels lining is divided into two groups including without and with invert. In this
project, there are two groups based on rock mechanics properties. The tunnel lining is without invert in E at and Edtb
sections and it has invert for FZ section. The numerical modelling was performed based on section 3.2.2 under static
and dynamic conditions (ODE & MDE) for three tunnels sections by considering hyperstatic reaction concepts. Two
loading are considered for analysing: a) dead load plus excavation load, and b) earthquake load. Then, forces are
accumulated in nodes based on superposition principle. Figures 9 to 12 show the state of axial force and bending
moment under ODE loading in tunnels lining for E at and FZ sections. The calculated thickness of lining is presented
for three sections in Table 5, according to analyses and considering the proper safety factor.

Table 6. The thickness of lining and armatur for three sections

Unit Thickness (cm) Armature
Eat 40 2×φ10@20
Edtb 40 2×φ18@20
FZ 60 2×φ25@10

250
Civil Engineering Journal Vol. 2, No. 6, June, 2016

Figure 9. The state of a) axial force and b) bending moment in tunnel lining under static loading (E at section)

Figure 10. The state of a) axial force and b) bending moment in lining under earthquake loading (Eat section)

Figure 11. The state of a) axial force and b) bending moment in tunnel lining under static loading (FZ section)

Figure 12. The state of a) axial force and b) bending moment in tunnel lining under earthquake loading (FZ section)

251
Civil Engineering Journal Vol. 2, No. 6, June, 2016

The maximum axial force and bending moment for ODE and MDE conditions are illustrated in Table 2. It is clear that
the maximum of axial force and bending moment occur in FZ section (i.e. with close invert), and more forces are
created in MDE condition for the earthquake loading.

Table 7. Maximum axial force and bending moment for all load combinations
ODE MDE
Static Load Dynamic Static Load Dynamic Static Load Dynamic Static Load Dynamic
Unit
(D+EX) Load (D+EX) Load (D+EX) Load (D+EX) Load
N M N M N M N M
Eat 319.6 3 1.8 5.3 282.3 19.9 2.2 6.4
Edtb 214.5 9.7 1.8 5.3 195.9 15.1 3.4 Fault Zone - MDE - Left Side
9.9
1.0Dead+1.0Rockload+1.0EQ
FZ 380.5 125.1 10.7 24.9- MDE - Edtb 205.4
Right Tunnel 72.7 19.9 48.8
Right Tunnel - MDE - EatBending moment in
* N= Axial 1.0Dead+1.0Rockload+1.0EQ
force in ton, M= 1.0Dead+1.0Rockload+1.0EQ

6.2. Interaction Diagram between Axial Force-Bending Moments

The bearing capacity of the column cross section can be determined from the interaction diagram moment-axial
force (P-M).The P-M interaction diagram is a suitable tool for designing and calculating the ultimate capacity of
tunnel lining sections in load combination conditions of axial force with bending moment. Figure 13 shows P-M
diagrams for different sections (Eat, Edtb, and FZ). According to Figure 13, tunnels lining is stable under ODE and
MDE conditions.
P-M Diagram
P-M Diagram P-M Diagram
500 600 1000

2× Φ 10@ 20 , D=40cm 2× Φ 25@ 10 D=60cm

500 2× Φ 18@ 20 D=40cm
400 ODE 800 ODE
ODE
MDE
MDE
400 MDE
300 600
300
Pn (ton)

Pn (ton)
400
Pn (ton)

200
200

100
200
100

0 0
0
-40 -20 0 20 40 -50 -40 -30 -20 -10 0 10 20 30 40 50 -200 -100 0 100 200
-100 -200
-100

Mn (ton-m) -200
-400
Mn (ton-m)
Mn (ton-m)

Figure 13. P-M diagrams for different sections (Eat, Edtb, and FZ)

7. Conclusion
In this paper, the application of HRM is presented for designing permanent tunnel lining in which two dimensional
analyses can be considered. Designing is done for Manjil-Rudabar freeway project, Tunnel No. 2. The numerical
analyses were performed for ODE and MDE loading conditions. Then interaction diagram between axial force and
bending moment was used for investigating the capacity of tunnel lining. The thickness of tunnel lining and armature
was calculated for three sections based on induced forces in tunnel lining. These forces were different in every section
according to the load combinations, rock mechanics properties, lining properties and overburden. The numerical
results showed that the forces in tunnel lining for MDE condition is approximately 50% more than ODE condition in
earthquake loading (Figure 14). This numerical processing presented that HRM is a proper, fast, and practical method
for tunnel engineers.
25 60
Axial Force (ODE)
Axial Force (MDE)
50
20 Bending Moment (ODE)
Bending Moment (MDE)
Bending Moment (ton.m)

40
Axial Force (ton)

15

30

10
20

5
10

0 0
Eat Edtb FZ
Figure 14. The axial force and bending moment in earthquake loading (ODE and MDE condition)

252
Civil Engineering Journal Vol. 2, No. 6, June, 2016

8. References
[1] Limkatanyu, S., Kuntiyawichai, K., Spacone, E., Kwon, M.H. “Nonlinear winkler-based beam element with
improved displacement shape functions”. KSCE J. Civil Eng. 17(1), 192–201, 2013.
[2] Do, N.A., Dias, D., Oreste, P., Djeran-Maigre, I. “The behavior of the segmental tunnel lining studied by the
hyperstatic reaction method”. Eur. J. Environmental Civil Eng. 18(4), 489–510, 2014.
[3] Barpi, F., Barbero, M., Peila, D. Numerical modelling of ground tunnel support interaction using bedded-beam-
spring model with Fuzzy parameters. Gospodarka Surowcami Mineralnymi 27, 71–87, 2011.
[4] Orestre, P. “A numerical approach to the hyperstatic reaction method for the dimensioning of tunnel support”.
Tunnelling and Underground Space Technology. 22(2), 185–205, 2007.
[5] Oreste, P. “A probabilistic design approach for tunnel supports”, Computers and Geotechnics 32, 520–534, 2005.
[6] Do, N. A., Dias, D., Oreste, P., and Maigre, D. I., “A new numerical approach to the hyperstatic reaction method
for segmental tunnel lining”, International Journal of Numerical and analytical method in Geomechanics, DOI:
10.1002/nag.2277, 2014.
[7] Kim, J. S, Kim, M. K, · Sam Dong Jung, S.D, “Two-dimensional numerical tunnel model using a Winkler-based
beam element and its application into tunnel monitoring systems”. Cluster Computer, DOI: 10.1007/s10586-014-
0418-4, 2015.
[8] DAUB (German Committee for Underground Construction). Recommendation for static analysis of shield
tunnelling machines. Tunnel, vol. 7, pp. 44–58, 2005.
[9] Association Française des Tunnels et de l’Espace Souterrain, Consideration on the usual method of tunnel lining
design, Tunnels & Ouvrages Souterrains, vol. 14, mars/avril, 1976.
[10] US Army Corps of Engineering (USACE) Tunnels and shafts in rock. Engineering Manual 1110-2-2901,
Washington (USA), 1997.
[11] Japan Society of Civil Engineering, 2006, Standard Specifications for Tunnelling-2006: Shield Tunnel.
[12] Panet, M. AFTES, Recommendation on the Convergence-Confinement Method, 2001.
[13] Mashimo H, Ishimura T. “Evaluation of the load on a shield tunnel lining in gravel”. Tunnelling and
Underground Space Technology; 18:233–241, 2003.
[14] Hashash YMA, Hook JJ, Schmidt B, Yao JI. “Seismic design and analysis of underground structures”. Tunnelling
and Underground Space Technology;16:247–93, 2001.
[15] Park, K.-H., Tantayopin, K., Tontavanich, B., and Owatsiriwong, A. “Analytical solution for seismic-induced
ovaling of circular tunnel lining under no-slip interface conditions: A revisit”, Tunnelling and Underground Space
Technology, 24(2): p. 231-235, 2009.

253