Finite Element Analysis of the Stability of Tunnel Surrounding Rock with Weak Rock Layer

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.


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 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.

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.

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) (2) Here, f is the slope of the linear yielding plate on the p-t stress plane, and it is usually called as the friction angle of materials.
d is the cohesive strength of materials, and , 0 c is the uniaxial compressive yielding stress, and 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.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.

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.

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 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 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 = /4 (Hock and Brown, 1997).

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.

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.

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 didn't 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.

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 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.

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.
angle in the p-t plane.The Mohr-Coulomb rule can be denoted as(ABAQUS Inc, 2006)    For the triaxial compression, the formula (1) can be rewrote as 0 formula (3) with the formula (4), and 0 c corresponding with the Drucker-Prager compression cone and the Mohr-Coulomb rules are sin the Drucker-Prager tensile cone and the Mohr-Coulomb rules respectively 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 )

Figure 1 .
Figure 1.Weak Rock layer on the Top of the Tunnel (the lines in the model present the range of weak rock layer)

Table 1 .
Physical mechanic parameters based on the modification of the basic quality index [BQ] of surrounding rock mass