Vol. 3, No.

12 Modern Applied Science

22
Finite Element Analysis of the Stability of Tunnel Surrounding Rock
with Weak Rock Layer
Yangsong Zhang
Nanjing University of Science and Technology, Nanjing 210094, China
Tel: 86-25-8431-5773 E-mail: zysnj318@126.com
Xiaozhao Li
Nanjing University, Nanjing 210093, China
Tel: 86-25-8359-7196 E-mail: lixz@nju.edu.cn
The research is financed by the Natural Science Foundation of China (No. 40872172). (Sponsoring information)
Abstract
The evaluation of local stability of the tunnel surrounding rock mass with weak rock layer is very important for the
economy and security of the tunnel construction. Based on the rock tunnel geological conditions in the construction, the
finite element model of the tunnel surrounding rock mass stability analysis with weak rock layer is established in the
article by selecting the ABAQUS program linear D-P materials model as the gneiss yield criterion. The computation
result indicates that the explicit dynamic analysis program could simulate the whole process of the instabilization and
failure of the weak rock layer, and avoid the convergence existing in the general finite element computation, and the
analysis model and program could offer effective measure for the analysis of local stability of the tunnel surrounding
rock mass with weak rock layer.
Keywords: Rock tunnel, Weak rock layer, Explicit finite element analysis, Stability
1. Introduction
The empirical method and the numerical analysis method are two main methods to analyze the stability of the tunnel
surrounding rock mass at present. The empirical method is to mainly qualitatively analyze the stability of tunnel
surrounding rock mass according to the quality classification result of rocks, and the quality classification methods
extensively used in the world include the Geomechanics Classification System (RMR, Bieniawski, 1989), the Q-system
(Barton et al., 1974 & Barton, 2002), the Geological Strength Index (GSI, Hoek and Brown, 1997), and the Chinese
national standard BQ classification method which is extensively used in China (The Ministry of Water Resources
People Republic of China, 1995). These classification methods are the results of researches in many tunnel cases, and
have been applied in many evaluations about the stability of practical tunnel surrounding rock mass. But these methods
could not quantitatively compute the stress redistribution aroused in tunnelling works (Basarir et al., 2005), and could
not reasonably consider the local stability of unequal tunnel surrounding rocks with weak rock layer. On the other hand,
the finite element method could offer the displacement of the surrounding rock mass after tunnel is excavated, the stress
redistribution, the plastic zone distribution, and other useful information, which can help to quantitatively evaluate the
stability of the surrounding rock mass (Zhu Weishen et al., 2003 & Basarir et al., 2005). However, similarly, the
reasonability of the results computed by many numerical methods such as finite element method depends on whether
the space distribution of rock medium in the established numerical analysis models is close to the actual situation to
large extent, and the reasonability of input parameters in various rock mass computations.
For the surrounding rock mass that the rock mass have good quality and high strength, when analyzing the stability of
tunnel surrounding rock mass by the finite element method, the security coefficient is often used to denote the factor of
safety of stability, but to solve the factor of safety of the stability of tunnel surrounding rock mass, the strength
reduction method is mainly adopted at present (Zheng Yingren and Zhao Shangyi, 2005). But when computing the local
stability of surrounding rock mass with weak rock layer, and the strength parameter of rock mass is not reduced, the
weak rock layer may be failured locally, and it is very important to mainly study when the weak rock layer begins to
become into instabilization and failure, and the influences of the whole stability of tunnel surrounding rock mass after
Modern Applied Science December, 2009
23
the weak rock layer is failured, and the local supporting opportunity, and the supporting effect. There are not systematic
research results about these problems at present. In the article, taking geological conditions of Yuntaishan Road Tunnel
in Lianyungang of Jiangsu, China, as the example, the finite element computation analysis method about the local
stability of the tunnel surrounding rock with weak rock layer is discussed as follows.
2. Explicit finite element analysis of the whole process of instabilization and failure of tunnel surrounding rock
mass
2.1 Instabilization criteria of tunnel surrounding rock mass
At present, the criterion that the finite element method is adopted to compute extreme status of instabilization and
collapsed of rock-soil mass includes following kinds.
(1) The distribution of plastic zone. The size of the plastic zone of surrounding rock mass in excavaiing is one standard
to evaluate the stability of surrounding rock mass (Ferrero Anna Maria, Migliazza Maria, Giani Gian Paolo, 2004). But
for tunnels, the concept of the size of the plastic zone of surrounding rock mass or whether it is connected or not is not
specific like the slope and the specific borderlines are difficult to be confirmed.
(2) The computating that the finite element is not convergent indicates that the rock-soil mass has been destroyed
(Dawson E.M., Roth W. H., 1999 & Griffiths D.V., Lane, P.A., 1999). According to this standard, different finite
element programs and different convergence standards will produce different computation results. And in the finite
element simulation analysis of tunnelling, the result of the general finite element program is not convergent often, so the
further change of the surrounding rock mass status after local destroying can not be known.
(3) The sign that the rock-soil mass is collapsed should be that the slide moves infinitely, and the slide surface of
rock-soil mass and the displacement change and development infinitely. If the numerical method can compute the
deformation of surrounding rock mass, and follow the whole process of the deformation and instability of surrounding
rock mass, so this rule is very good. And it can monitor various points of the surrounding rock mass of different classes
in different parts of the tunnel, which can not only know the local instability, but know the total instability induced by
the local instability. It specially suits for the analysis and judgment of the local stability of tunnel surrounding rock mass
with weak rock layer.
2.2 Constitutive model of rock mass
Mohr-Coulomb strength rule and the Drucker-Prager strength rule are two kinds of usual rules in the rock-soil mass
computation. The strength parameter of the Mohr-Coulomb is easily obtained usually, but in the extreme analysis of
rock-soil mass, its computation result often is not convergent at the sharp corner. In the article, the linear
Drucker-Prager model is adopted. The linear Drucker-Prager yielding function and the plastic potential energy
respectively are (ABAQUS Inc, 2006)
0 tan = = d p t F | (1)
tan p t G = (2)
Here,
(
(

|
|
.
|

\
|
|
.
|

\
|
+ =
3
1
1
1
1
2
1
q K K
q t

, S r
2
9
( = )
3
S S , I p S + = o
) , (
i
f u | is the slope of the linear yielding plate on the pt stress plane, and it is usually called as the friction angle of
materials.
d is the cohesive strength of materials, and
0
tan
3
1
1
c
d o | |
.
|

\
|
= ,
0
c
o is the uniaxial compressive yielding stress, and
) , (
i
f K u is the ratio between the uniaxial extension yielding stress and the uniaxial compressive yielding stress, so it
control the dependence of the middle main stress of the yielding plane.
) , (
i
f u is the dilation angle in the pt plane.
The Mohr-Coulomb rule can be denoted as (ABAQUS Inc, 2006)
0 cos 2 sin ) (
3 1 3 1
= + + o o o o c (3)
For the triaxial compression, the formula (1) can be rewrote as
0
tan
6
1
1
tan
3
1
1
) (
tan
3
1
2
tan
0
c 3 1 3 1
=
+

+
+
+ o
|
|
o o
|
|
o o (4)
Vol. 3, No. 12 Modern Applied Science
24
By Comparing the formula (3) with the formula (4), | and
0
c
o corresponding with the Drucker-Prager compression
cone and the Mohr-Coulomb rules are

|
sin 3
sin 6
tan

=
(5)

o
sin 1
cos
2
0
c

= c
(6)
In the same way, | and
0
c
o corresponding with the Drucker-Prager tensile cone and the Mohr-Coulomb rules
respectively are
) sin 3 (
sin 6
tan

|
+
=
(7)
|
|
o
tan 3
) tan 6 ( cos
0
c

= c
(8)
| tan
3
1
1
1
+
= K
(9)
Fit the average value of the compression cone and the tensile cone of the Drucker-Prager yielding rule with the
Mohr-Coulomb rule, and take the average | tan value of the formula (5) and the formula (7) to compute ) , (
i
f K u
according to the formula (9), and 0.778 ) , (
i
f K u 1, and in the same way, ) , (
i
f K u controls the dependence of the
middle main stress of the yielding plane. Large numerous of comparisons and computations indicated that the
parameters computed by the linear Drucker-Prager yielding rule and above methods and Mohr-Coulomb had good
fitting precision.
2.3 Analysis approaches
In general finite element analysis, there are three methods to simulate the tunnelling, i.e. (1) gradually reducing the
elasticity modulus of the excavating rock mass, (2) adopting the killing element or removing element method, (3)
substituting the rock mass in excavating tunnels by the fixed boundary conditions, and first computing the initial stress
to balance the reverse stress of above borderlines, and in subsequent computation, putting the borderline reverse
function on the responding nodes, and removing corresponding fixed borderline conditions on the tunnelling plane, and
computing the ground stress balance, and in subsequent computation, gradually reduce the reverse stress on the
tunnelling face to simulate the stress release in tunnelling according to certain proportion. In ABAQUS/Explicit analysis,
above two methods can not obtain satisfactory excavating modeling effect. Though the third method is complex, but it
fits for the ABAQUS/Explicit analysis, and it is convenient to accurately control the stress release.
According to the surrounding rock mass condition, each step to excavate tunnel is 3.0m, and the detailed simulation
process includes following steps.
(1) Compute the reverse forces of various nodes on the walls of the excavating tunnel. It only needs to implement the
fixed borderline conditions on the tunnelling face in the step of the initial geostatic stress balance, so the reverse forces
of various nodes can be obtained.
(2) Implement the reverse forces of various nodes in the computation of first step on the corresponding nodes to
compute the initial geostatic stress balance again.
(3) Gradually reduce the reverse forces of various nodes on the tunnel surfaces to zero step by step, to simulate the
stress release of rock mass in each time. The final simulation computation includes the release of the reverse forces of
the nodes on the surface of surrounding rock mass in front of the tunnel.
The first step is implemented in ABAQUS/Standard, and the second step and the third step are implemented continually
in ABAQUS/Explicit.
3. Example analysis
3.1 Geological conditions and physical mechanical parameters of surrounding rock masses
The computation antitype is the surrounding rock mass section of the tunnel, mainly includes the solid metamorphic
rocks, and the saturated uniaxial compressive strength of rocks is R
C
=50~100MPa. For two main joints strikes are NE
and NW in the rock mass, the joint plane is planar and coarse, and presents close-parting opening status without filling.
The schistose structural plane occurrence average dip is 155
0
, and the average obliquity angle is 28
0
. Some parts with
Modern Applied Science December, 2009
25
the green schist weak rock layer unveil, and the saturated uniaxial comprehensive strength of rocks is R
C
=5~17MPa,
and the schistose structural plane of rock mass grows. The weak rock distribute as the intercalated status along the
schistose structural plane direction. The field survey shows that when the weak rock layer is distributed on the top of the
tunnel, if the local support measure is not adopted, the local instabilization and failure of surrounding rock mass easily
occur. According to the field survey, the geometric relationship model between the weak rock distribution and the tunnel
is seen in Figure 1.
The physical mechanical parameters confirmed by the modified basic quality index [BQ] of surrounding rock masses
are seen in Table 1, and the Table 1 also shows the internal friction angle | and the compressive initial yielding stress
0
c
o of the linear Drucker-Prager rule.
For the linear Drucker-Prager rule in the ABAQUS program, to ensure that the yielding plane is out-extrusive, the value
of K should equal or exceed 0.778, so when value of k confirmed by the formula (6) is less than 0.778, K=0.778. Take
nonassociated flow rules, suppose the expansion angle is ) , (
i
f u =/4 (Hock and Brown, 1997).
3.2 Finite element analysis model
In the computation, for the model, the length is 200m (x axis), the height is 120m (z axis), the depth is 45m (y axis), the
span of the tunnel is 16m, and the height of the tunnel is 9.05m, and the simulated excavating depth of the tunnel is 15m,
and whole sections are dug. The length from the top of the tunnel to the top of the model is 41m, and the rock mass of
200m covering on it adopts the average press on the surface of the model to be simulated. As seen in Figure 2, the
model is divided into 121714 4-node linear tetrahedron elements and 21815 nodes.
3.3 Results of computation and analysis
According to the numerical simulation results, through observing the vertical displacement of the tunnel top and the
distribution of the plastic zone of surrounding rock masses, following rules can be summarized.
3.3.1 Vertical displacement of tunnel top
The distribution of the final vertical displacement (U3) of tunnel surrounding rock mass obtained in the finite element
computation is seen in Figure 3. From Figure 3, the vertical displacement value from the top of the tunnel in the weak
rock mass ([BQ]=240) distribution zone obviously exceeds the corresponding values surrounding the solid rock zone
([BQ]=500). The maximum of the vertical displacement of the tunnel top is at the node A, and it is 17.23mm in the
weak rock layer, and the minus denotes the displacement direction downwards, which is observe of the direction of Z
axis. But the displacements in the solid rock zone of the tunnel top are small, for example, the vertical displacement of
the node C far from the weak rock distribution is 6.38mm.
In the finite simulation computation of tunnelling, the load release proportion factor (LPD) and the computation time
have same meanings. With the increase of the excavating step computation time, the node force on the excavating
borderline gradually reduces to zero. So from the result of the finite element computation, the relationship between the
node displacement on the excavating plane and the load release proportion factor can be easily obtained.
Figure 4 shows the relationship curve between the vertical displacement and the load release proportion factor of three
places on the top of tunnel. Obvious inflexion occurs in the curve of the node A (weak rock layer) with the increase of
the load release proportion factor when LPD is close to 0.9, i.e. the displacement increases quickly, so according to the
result of the finite element computation, comparing with the extreme status rule 2, if any support measure is adopted,
the weak rock layer on the top of tunnel will be collapsed. That accords with actual situation, in the tunnelling
construction of Yuntaishan in Lianyungang, in the solid II-class surrounding rock mass, because the quality of rock
mass was good, the systematic primary support had not been adopted according to relative standards, and the stability
was good when the weak rock layer didnt exist, and the deformation of surrounding rock mass was small, but for the
same solid II-class surrounding rock mass, when the top of tunnel had the weak rock layer with big depth (1~3m) and
effective local support had not been adopted, several collapses occurred in the weak rock surrounding layer on the top
of tunnel in certain time after the tunnel was excavated.
From Figure 4, the vertical displacements of the node B and the node C in the solid rock zone are obviously less than
the vertical displacement of the node A, and because the node B is close to the edge of the weak rock layer, the vertical
displacement is larger than the node C far from the weak rock mass because of the influence of the weak rock plastic
failure. The vertical displacements of the node C and the node B increase slowly, which indicates that they have not
achieved the extreme status.
3.3.2 Plastic zone of surrounding rock mass
Figure 5 shows the distribution of the plastic zone of surrounding rock mass. From Figure 5, the plastic zone occurs in
the top of tunnel and the upper weak rock layer on the excavating plane. The solid rock distribution zones all have not
achieved the plastic failure. From the development process of the plastic zone, when the load proportion factor is 0.9,
obvious plastic variance begins in the weak rock layer, which accords with the time when the inflexion of the vertical
Vol. 3, No. 12 Modern Applied Science
26
displacement-load proportion factor curve of the node A, and that indicates that the quick increase of displacement is
the induced by the obvious plastic variance in surrounding rock mass.
4. Conclusions
Combining with the linear Drucker-Prager yielding rule, the ABAQUA/Explicit program is adopted to analyze the
stability of the tunnel surrounding rock mass with the weak rock layer. By the research, following conclusions can be
obtained.
(1) The numerical analysis program in the article can be used to compute the whole process of the deformation and
failure of the tunnel surrounding rock mass, and it can avoid the problem that the results of the finite element
computation is not convergent, so it specially suits for the computation of the stability of the tunnel surrounding rock
mass with weak rock layer.
(2) It is reasonable and feasible to apply the ABAQUS/Explicit program to compute and adopt the rule that the
deformation increases quickly as the extreme status of the instabilization and failure of tunnel surrounding rock mass.
On the one hand, the computation program can ensure to obtain the inflexion that the displacement quickly increases,
and on the other hand, the results of the computation and judgment are coherent with the actual situation.
(3) In the article, the quality classification standard of rock mass is adopted to evaluate the strength parameters of
Mohr-Coulomb rule, and these parameters are converted into the strength parameters of the linear Drucker-Prager
yielding rule by the method in the article. Then they are used to compute the stability and the displacement of tunnel
surrounding rock mass, and the computation results is coherent with the actual survey result.
References
ABAQUS Inc. (2006). ABAQUS Users Manual, Version 6.6.
Barton, N., Lien, R., Lunde, J. (1974). Engineering classification of rock masses for the design of tunnel support. Rock
Mech. No.6(4). P.189-236.
Basarir, H., Ozsan, A., Karakus, M. (2005). Analysis of support requirements for a shallow diversion tunnel at Guledar
dam site, Turkey. Engineering Geology. No.81. P.131-145.
Bieniawski, Z.T. (1989). Engineering Rock Mass Classification: A Complete Manual for Engineers and Geologists in
Mining, Civil, and Petroleum Engineering. Wiley, New York. 251 pp.
Dawson E.M. and Roth W. H. (1999). Drescher A. Slope stability analysis by strength reduction. Geotechnique. No.49
(6). P.835-840.
Ferrero Anna Maria, Migliazza Maria, Giani Gian Paolo. (2004). Analysis of tunnel stability: comparison between
continuous and discontinuous approaches. Int J Rock Mech. Min. Sci. No.41(3). CD-ROM. P.1-6.
Hoek E and Brown ET. (1997). Practical estimates of rock mass strength. Int J Rock Mech Min Sci Geomech Absrr.
No.34(8). P.1165-1186.
The Ministry of Water Resources People Republic of China. (1995). The Chinese Standards-Standard for engineering
classification of rock masses (GB50128-94). Beijing: China Planning Press (in Chinese).
Zheng, Yingren and Zhao, Shangyi. (2005). Limit state finite element method for geotechnical engineering analysis and
its applications. China civil engineering journal. No.38(1). P.91-104 (in Chinese).
Zhu, Weishen, Li, Shucai, Li, Shuchen, et al. (2003). Systematic numerical simulation of rock tunnel stability
considering different rock conditions and construction effects. Tunnelling and Underground Space Technology. No.18.
P.531-536.
Table 1. Physical mechanic parameters based on the modification of the basic quality index [BQ] of surrounding rock
mass
Type of rock [BQ]
Class
of
rock
mass
Deformation
modulus
E
mass
(GPa)
Passion
ratio
Density of
gravity
(kN/m
3
)
Mohr-Coulomb
criterion
Drucker-Prager
criterion
(
-
) c (MPa) (
-
)
0
c
o (MPa)
Solid rock 500 II 26.5 0.23 26.0 55 1.8 59.1 7.13
Weak rock 240 V 1.25 0.36 22.0 26 0.19 41.59 0.504
Modern Applied Science December, 2009
27
Figure 1. Weak Rock layer on the Top of the Tunnel (the lines in the model present the range of weak rock layer)
Figure 2. Finite Mesh of the Tunnel Computation Model
Figure 3. Vertical Displacement U3 of Tunnel Surrounding Rock mass (Unit: m) (Nodel A: weak rock layer,
Nodel B: solid rocks near to weak rock layer, Nodel C: solid rocks far from weak rock layer)
0.0?
0.0!
0.0!
0.00
0
0 0.? 0.+ 0.b 0.8 !
lD
|
3

m
\od' ^ \od' b \od' C
Figure 4. Relationship Curve between Vertical Displacement of Surrounding Rock Mass
on the Top of Tunnel and the Load Proportion factor (LPD)