Abstract
To study the tunnel stability at various static water pressures and determine the mechanical properties and deformation behavior of surrounding rock, a modified effective stress formula was introduced into a numerical integration algorithm of elastic-plastic constitutive equation, that is, closest point projection method (CPPM). Taking the effects of water pressure and seepage into account, a CPPM-based formula was derived and a CPPM algorithm based on Drucker-Prager yield criterion considering the effect of pore water pressure was provided. On this basis, a CPPM-based elastic-plastic numerical analysis program considering pore water pressure was developed, which can be applied in the engineering of tunnels and other underground structures. The algorithm can accurately take the effects of groundwater on stability of surrounding rock mass into account and it can show the more pronounced effect of pore water pressure on stress, deformation, and the plastic zone in a tunnel. The stability of water flooding in Fusong tunnel was systematically analyzed using the developed program. The analysis results showed that the existence of groundwater seepage under tunnel construction will give rise to stress redistribution in the surrounding rock mass. Pore water pressure has a significant effect on the surrounding rock mass.
1. Introduction
As a natural material, rock has numerous fractures, voids, and pores (i.e., fissures are rich in soluble rock), within which groundwater accumulates and is transported which affects deformation and failure of the rock. Groundwater is recognized as a dominant factor that influences the stability of rock engineering. In the construction of water-rich rock tunnels with a large buried depth, the surrounding rock mass is disturbed and damaged by construction, resulting in frequent accidents including flooding, mud slides, and leakage as a result of the coupling effects of high water pressure and high geostress. In the past, these misfortunes have resulted in serious economic damage and casualties [1–3]. Statistics concerning water-laden railway tunnels in China revealed that more than 80% of these tunnels either under construction or in operation had suffered various degrees of adversities including flooding and mudslides during construction and a full 30% of these tunnels are subject to groundwater pressure [4, 5]. As an additional consequence, these water-related mishaps induced by construction disturbance will also lead to a decline of groundwater levels with concomitant dry-up of springs and wells, leading to environmental deterioration including flow cutoff of rivers and brooks and withering of vegetation. This will directly affect local industrial and agricultural production as well as people’s lives. Therefore, the study on stability of tunnels and other underground engineering under high water pressure is of great significance.
The deformation and stress response of the rock mass surrounding a tunnel under high water pressure has been studied by many researchers in China and abroad. Using field testing, Bruno and Nakagawa [6] experimentally determined that water pore pressure had dual effects on crack growth and coalescence. Xing et al. [7] analyzed the influence of high confining pressure and high water pressure on deformability, strength, and the ductile-to-brittle transition of brittle rock using complete stress-strain triaxial compression tests. Xu et al. [8] conducted triaxial compression tests to study the effects of water pore pressure on the mechanical properties of sandstone under water-filling conditions. Bao et al. [9] studied the effect of high water pressure on the rock deformation and rock acoustic emission using experiment and numerical analysis. She and Cui [10] conducted rock creep tests under high pore pressure. Using analytical and numerical calculations, Li et al. [11] deduced an elastic-plastic analytical solution for a deep tunnel in a water seepage field. Based on effective stress principles and pore pressure in each element, Wang [12] and Wang et al. [13] simulated the deformation caused by constant pore pressure acting on rock and the localization process under various confining pressures and pore pressures using FLAC software. Fang et al. [14] reviewed analytical solutions for water pressure distribution under steady seepage in a tunnel and comprised a fluid-solid coupling simulation. These research results are valuable reference for the study of the instability mechanism of water flooding in deep tunnels of water-rich zone and the development of prevention measures. However, there are still no theoretical system and comprehensive conclusions due to the complexity of this problem. Therefore, further study and discussion are necessary.
In this paper, a modified effective stress formula was introduced into the CPPM and the integration algorithm for elastic-plastic constitutive equation. Then elastic-plastic numerical analysis program was developed considering the effect of water pore pressure in tunnels and other underground engineering structures. The stability of the tunnel under various water pressure conditions was analyzed based on the developed program. The variation laws of stress, strain, and plastic zone were also studied. This developed program is also used in the analysis of true engineering named Fusong tunnel of China.
2. Elastic-Plastic Numerical Program Considering Surrounding Rock Mass Subject to Water Pore Pressure
2.1. Closest Point Projection Method with Water Pore Pressure
The closest point projection method is a plastic finite element method algorithm, which consists of two parts, the initial elastic trial step and the plastic corrector step. They are driven by the total strain increment and plastic parameter increment, respectively. The plastic strain and the internal variables are kept fixed in the elastic trial step. The total strain remains unchanged in the plastic corrector step [15, 16].
Pore pressure distribution designated from the possible mechanical deformation is sufficient for solving the problem of static groundwater. The constant pore pressure is established by considering the action of pore pressure for the stress calculation by solving the effective stress of all nodes. Water pore pressure is constant throughout the whole calculation process and the pore pressure affects the rock primarily through effective stress. The pore pressure distribution is consistent with the initial condition and not the influence of mechanical deformation. Some scholars proposed improved effective stress formula according to rock mass as follows [17]:where is effective stress tensor, is the total stress tensor, and is equivalent pore pressure coefficient; it is decided by rock pore and crack growth degree. is pore water pressure and is equivalent pore pressure.
The effect of water pressure on the strength of the material is represented by Mohr’s circle. The stress condition is shown in Figure 1 before the introduction of water pore pressure . All the normal components of the stress tensor are effectively reduced through the pore pressure after the introduction of pore pressure as shown in Figure 1. Mohr’s circle moves towards the left under the action of the water’s pore pressure and to the destruction of the material.
The closest point projection method for the Drucker-Prager criteria, the updated stress at the step, is as follows:where is the return vector, is trial stress, is the plastic multiplier, and is the flow vector. There are two possible forms of the return mapping algorithm according to the flow vector at the smooth cone or at the apex.
Returning to the smooth cone, the accumulated strain iswhere is the accumulate strain at the step, is the increment of accumulate strain, and is the constant relevant to the friction angle.
The consistency condition iswhere is the deviator stress tensor, is hydrostatic pressure, is a constant related to the friction angle, and is hardened modulus.
Solving formula (4), the calculated , updated stress at the step is
Returning to the apex singularity as shown in Figure 2, stress is only relevant with a body strain increment and an accumulated plastic strain :
The consistency condition is
Solving (7), we obtain the updated stress :
The above calculation assumes the return to the apex singularity but actually also solves to judge whether to return to the smooth cone or the apex singularity, as shown in Figure 3:
Stress returning to the cone of the yield surface is effective. Otherwise, it returns to the apex of the yield surface.
The total stress of (1) at the closest point projection method of the Drucker-Prager criteria can be obtained from the calculation formula of the closest point projection method considering the action of pore water pressure; then we can calculate the effective stress .
Assuming the water is incompressible, nonsteady seepage continuity equation is as follows according to Darcy’s law:where is permeability head. , , and are the spatial coordinates, is the time coordinate, , , and are the permeability coefficients whose directions are parallel to the -, -, and -axis, respectively, and is the water storage rate.
The initial condition of unconfined nonsteady seepage definite condition with free surface is as follows:
Head boundary and flow boundary are as follows:where is the known head boundary, is the known head boundary value, is the known flow boundary, and is the known flow boundary value.
2.2. Self-Development of the Program of Elastic-Plastic FEM
Self-development of the elastic-plastic program of Drucker-Prager criteria under the action of pore water pressure is based on the original elastic finite element program. The program was developed with an object-oriented C++ language through the debug in the Visual 6.0 environment. Dealing with grid division by ANSYS software and stress, displacement and other values are graphically displayed by the Tecplot software. Completion of the whole process of the solution is accomplished using the finite element program, which can analyze special engineering problems, considering the effect of water pressure based on the CPPM under the Drucker-Prager constitutive model as shown in Figure 4.
Making use of an incremental-iterative method in the nonlinear finite element, a piecewise linear solution is obtained that approximates the nonlinear solution. The total load is divided into a number of load increments . The resulting finite element incremental equation is
Solving (13) using the Newton-Raphson iteration scheme, the process is as follows:(1)Initialize the displacement , stress , and strain , respectively.(2)At load step (), have the load increment , to solve .(3)At the iteration of load step , where maxiter, have , to solve internal forces .(4)Use the formula for convergence test; if there is convergence, then return to step , to the next load step calculation; otherwise proceed to step .(5)Calculate the tangent stiffness matrix by solving the elastic stiffness matrix in the initial elastic stage, using the consistent tangent stiffness matrix in the plastic stage, solving the elastoplastic constitutive matrix .(6)Solve the constitutive equation (13), and obtain the displacement increment .(7)For each step, , , , and , the strain increment is .(8)The CPPM of the Drucker-Prager criteria is presented to solve the effective stress , solving the elastic-plastic factor . Using the Newton-Raphson iterative scheme can obtain approximate square convergence speed once again as follows:①At the time , calculate the strain increment and the internal variable , and obtain the predicted hydrostatic pressure and partial stress tensor in the elastic predicted state.②Judge, in the plastic state; if , then obtain the stress at the time . Exit, into the next load step calculation, or the plastic corrector step.③Assuming a return to the cone of yield surface, set the initial value .④Make use of the Newton-Raphson iteration scheme to solve , and judge whether there is convergence; if , obtain the updated stress ; otherwise, continue to the Newton-Raphson iteration.⑤Judge whether to return to the location of the trial stress; the above assumes a return to the cone of yield surface; if , predict if the stress needed to return to the cone of the yield surface is effective. Otherwise predict the stress needed to return to the apex singularity.⑥Set the initial value , solve using the Newton-Raphson iteration scheme, and judge whether there is convergence. If , get the updated stress ; otherwise, continue to Newton-Raphson iteration.(9)Calculate the stress after at the iteration step, turn to step , and continue to iteration.
2.3. Validation of Program and Numerical Simulation
In order to verify the reliability of the program, a plane-strain FEM model of a circular tunnel was established [18]. As shown in Figure 5, the model contains 613 elements and 671 nodes with a size of 10 m × 10 m. The simulated tunnel with a radius of 1 m is subject to a uniformly distributed load at the upper surface of the model. The elastic-plastic Drucker-Prager model is applied in the simulation. It is assumed that rock mass is homogeneous and fractures are not developed. The density of the rock mass is kg/m3 and the distributed load KN/m2. In addition, the simulation requires the specification of Young’s modulus , Poisson’s ratio , cohesion , a friction angle , and an angle of dilatation . All calculations were based on values of GPa, , MPa, = 30°, and = 30°. The equivalent coefficient of the pore water pressure, , was set to be 1.0.
When there is no pore water pressure, the uniformly distributed load is applied incrementally and the load increment will be applied over 5 steps. The load increment is 500 KN/m2 for each step. Calculated values of displacement using the program are compared with the values obtained from ANSYS software under the same condition. As shown in Figure 6, the displacement values at some nodes around the tunnel area are given.
(a) Displacement comparison in -direction
(b) Displacement comparison in -direction
According to displacement comparison at all nodes, the maximum relative error is 8.9% in -direction and 2.03% in -direction between the values obtained from the program and ANSYS. The conformity shown in these plotted data reflects the reliability of the program.
3. Mechanical Characteristics and Plastic Zone of the Circular Tunnel with Pore Water Pressure
3.1. Project Profile
Fusong tunnel is located in Jingyu county, Baishan city, in the east of Jilin province of China. The tunnel consists of two bores with centers separated by 13~35 m. The left bore has a total length of 1625 m with its starting and ending chainage of ZK275+170 and ZK276+795, respectively. The right bore has a total length of 1600 m with its starting and ending chainage of RK275+180 and RK276+780. The maximum excavation width of the tunnel is about 12 m. The tunnel height is 7.6 m. The area where tunnel is located belongs to the Changbai mountains foothills area, which is an overlying Quaternary loose formation. Outcropped strata in the area consist primarily of Jurassic lake-river facies clasolite, andesite and Tertiary-Quaternary basalt. In addition, Archaeozoic and Upper Proterozoic strata are sporadically distributed. Groundwater in the tunnel area mainly consists of pore phreatic water in a loose accumulation layer and fissure water in bedrock which is recharged from precipitation and discharged into a gully. The area has a continental climate, with a dry and severely cold winter and a hot summer. There is a relatively abundant rainfall from June to August with an average annual precipitation of 939.4 mm. There is a fault near RK275+390 at the entrance of Fusong tunnel where a seasonal river runs over this section of the tunnel. Obviously, the influence of the river water needs to be considered with respect to the tunnel stability during the rainy season. Therefore, a section in right tunnel from RK275+173 to RK275+390, with a total length of 204 m, was selected for study.
3.2. Seepage Analysis of the Tunnel
3.2.1. Seepage Analysis Model and Analysis Process
For the study of the tunnel construction stability at the entrance taking seepage into consideration, a practical rule of thumb is that the variation of stress and displacement induced by unloading in the circular tunnel which is excavated in a homogeneous elastic infinite domain is less than 1% over a radius being 5 times tunnel radius and 5% over a radius being 3 times tunnel radius. Therefore, in accordance with the specific requirements of the project, discretization error, and computational error in the FEM, the calculation scale was at least a radius being 3-4 times tunnel radius in each direction. The numerical model had dimensions of 204 m × 70 m × 90 m (length × width × height) and full-face excavation was applied in the numerical simulation (Figure 7).
In order to consider the effect of seepage on strength of the tunnel’s surrounding rocks, physical and mechanical parameters of the geomaterials in FEM model, including permeability and the saturated unit weight of soils, were specified, as shown in Table 1. The boundary conditions of seepage in the FEM model are shown in Figure 8.
The hydraulic head of the nodes in the vertical direction varied from 52 to 69 m and the surrounding rock mass below the hydraulic head was considered to be saturated before excavation. It was also assumed that the geologic model was impermeable at anterior, rear, and bottom surfaces. After excavation was completed, the boundary condition of the linings considered to be drained was applied. The process of analysis considering the tunnel construction procedures subject to seepage was composed of the following steps:(1)Initial seepage field analysis(2)Initial stress field analysis based on step (3)Tunnel excavation which was assumed to be made in one step under seepage and stress fields(4)Seepage field changes of surrounding rock mass after tunnel excavation.
3.2.2. Analysis of Seepage Characteristics in the Tunnel
The seepage characteristics in the tunnel at lower water levels were analyzed. Analysis of seepage characteristics in the tunnel shows that before excavation the pore pressure at the bottom agrees with total hydraulic head value and pore pressure increases from top to bottom. Figure 9 is the seepage velocity of tunnel section ZK275+290.
(a) -direction
(b) -direction
(c) -direction
(d) -direction
It is shown that groundwater maintains a balance, existing in the form of a hydrostatic pressure. After excavation the seepage field changed (linings were set to be drained), the pore pressure decreased dramatically, and the seepage velocity increased rapidly around the tunnel. Meanwhile, new permeable boundary was generated around the tunnel and pore pressure was zero. The seepage velocity vectors in the -direction were also analyzed, which means groundwater seeps into the tunnel along radial direction and the seepage flow along tunnel axis is weak due to a small hydraulic gradient. The seepage velocity of water in surrounding rock mass at ZK275+290 section is not generally large, with maximum value of up to 1.1 × 10−4 m/s.
3.3. Comparison Analysis of Tunnel Stability Subject to Seepage or No Seepage
The most commonly used tunnel surrounding rock safety indexes are as follows: limit displacement of surrounding rock; the scope of damage zone of surrounding rock; safety factor of surrounding rock. We analyzed the displacements, stress, and damage zone (plastic zone) of surrounding rock considering and not considering seepage based on the numerical simulation.
3.3.1. Influence of Seepage on Displacement of Surrounding Rock Mass
As shown in Figure 10, displacements of the surrounding rock mass at the ZK275+290 section are presented, which consider coupling or no coupling of seepage and stress. Coupling analysis of seepage and stress was considered during the tunnel construction as shown in Figures 10(a) and 10(c), while only rock self-weight was considered in Figures 10(b) and 10(d). Nephograms in Figures 10(b) and 10(d) indicate that displacements at both the tunnel roof and bottom are greater after excavation, horizontal and vertical displacements are distributed symmetrically about the center of the tunnel and the vertical displacement at the vault and bottom increased with a decrease in the buried depth of the calculated nodes. The horizontal displacement exhibited the same characteristics at either side of the arch waist, with a displacement value of approximately 1.7 mm. The maximum value of the horizontal displacement appeared at the arch waist. Therefore, a one-step tunnel excavation has an influence on the surrounding rock mass mainly in the form of a vertical deformation. Comparing between the calculation displacements and warning value, the stability degree of surrounding rock could be obtained. The limit displacement of surrounding rock of the section is decided by experience analogies and standard as 25 mm.
(a) -direction displacement considering seepage
(b) -direction displacement without considering seepage
(c) -direction displacement considering seepage
(d) -direction displacement without considering seepage
Nephograms shown in Figures 10(a) and 10(c) indicate that the maximum displacement of the surrounding rock mass appeared on tunnel floor and heaving displacement at the inverted arch was as much as 22.7 mm, while the settlement displacement was up to 12.2 mm. The surrounding rock displacements meet safety requirement. The heaving displacement was about 1.83 times the value that occurred when no seepage was considered. Comparison between Figures 10(a) and 10(b) shows that the following seepage increased the horizontal displacement of both sides of the arch waist 7 times when compared with no seepage. These findings show that the seepage has an adverse impact on stability of the tunnel. Therefore, strengthening of the support measures to prevent collapse of tunnel needs to be considered.
3.3.2. Influence of Seepage on the Stress of the Surrounding Rock Mass
As shown in Figure 11, the maximum and minimum principal stresses are developed depending on whether there is seepage or no seepage. Comparison of the nephograms shown in Figure 11 indicates that the seepage significantly contributed to an increase of the maximum principal stress and it also was found that the maximum principal stress below the tunnel floor is apparently greater than that at the tunnel roof when seepage is considered. However, the maximum principal stress in the surrounding rock mass was uniformly distributed. All these findings show that seepage stress redistribution in surrounding rock mass was induced by seepage and the influence was more significant on tunnel floor than on tunnel roof. Furthermore, it also was found that, in the case of no seepage, the local tensile stress appeared below tunnel floor due to seepage while there was only compressive stress in the surrounding rock mass.
(a) Maximum principal stress considering seepage
(b) Maximum principal stress considering no seepage
(c) Minimum principal stress considering seepage
(d) Minimum principal stress considering no seepage
(e) Shear stress considering seepage
(f) Shear stress considering no seepage
Developed shear stresses are shown in Figures 11(d) and 11(e). When there was no seepage, it is found that shear stress appeared along the arch waist of the tunnel and at both ends of the inverted arch of the tunnel. The maximum shear stress with X-shaped distribution appeared at the left arch foot. Due to influence of seepage, the distribution range of shear stress was wider and the value of the shear stress was larger. Therefore, seepage will cause redistribution and an increase of shear stress in the surrounding rock mass.
3.3.3. Influence of Seepage on the Equivalent Plastic Zone in the Surrounding Tunnel Rock Mass
The plastic zone produced during the tunnel excavation can be represented by an equivalent plastic strain, which can be used to illustrate the extent of surrounding rock damage. Figure 12 depicts the equivalent plastic strain of the surrounding rock mass after completion of tunnel excavation.
(a) Equivalent plastic zone under seepage
(b) Equivalent plastic zone under no seepage
As can be seen, the plastic zone with seepage is larger and more widely distributed than the plastic zone without seepage and even extends to the region below the tunnel. When seepage was not considered, the plastic zone occurred only at the arch foot. Therefore, it is important to detect water in advance during construction and antiseepage measures need to be promptly applied in order to prevent a wide range of the plastic zone. Particular attention should be paid to the support at the arch foot.
3.4. Analysis of Tunnel Stability Considering Seepage at Various Hydraulic Heads
There was a discordant contact zone between calcareous mudstone and calcareous siltstone at the ZK275+290 chainage section in the Fusong tunnel. Over this contact zone, there was a seasonal river that has a hydraulic connection with the water in the contact zone. Thus, the stability of tunnel’s surrounding rock mass and potential water flooding through the heading face can occur during seasonal flash floods. The scope of surrounding rock damage zone can reflect the safety of surrounding rock. The larger the damage zone, the more dangerous the surrounding rock. Consequently, the distribution characteristics of the plastic zone under seepage pressures of 0, 0.1 MPa, 0.25 MPa, 0.3 MPa, 0.5 MPa, and 0.6 MPa were analyzed based on the assumption that the rock mass in the plastic zone was homogenous with perfect elastic-plastic behavior. Simulations based on these assumptions showed that there will be a rapid increase in the displacement of geomaterials and penetration of the plastic zone in heading face when the geomaterials undergo plastic yield. This scenario will cause major mishaps such as major flooding, mudslides, and a potential failure of the tunnel structure. Figure 13 shows the plastic zones in surrounding rock mass and heading face of the tunnel at various different seepage pressures.
(a) Pore water pressure is 0 MPa
(b) Pore water pressure is 0.1 MPa
(c) Pore water pressure is 0.25 MPa
(d) Pore water pressure is 0.3 MPa
(e) Pore water pressure is 0.5 MPa
(f) Pore water pressure is 0.6 MPa
As shown in Figure 13(a), the distribution of the plastic zone when the seepage pressure is zero indicates that the area of plastic zone is small in the excavated area under just the weight of the earth and is mainly concentrated at arch waist and arch springing. The tensile stress appears at the bottom of heading face, but it is distributed in a small range in the middle of the heading face. As has been noted, the area of the plastic zone at the heading face accounts for about 25% of the total area of the heading face. With an increase in seepage pressure, the area of plastic zone gradually increases and the plastic zone at the side wall of the excavated section extends towards the heading face. Under these conditions, the depth and range of the plastic zone in the heading face also increased and the range of tensile stress at the bottom of the heading face continued to expand. When the seepage pressure increased to 0.25 MPa, the range of the tensile stress in the heading face increased by 50% and rock mass in the heading face deformed plastically over a small depth in the plastic zone. In addition, there was no penetration between the plastic zones in excavated section and the heading face. But when seepage pressure increased to 0.3 MPa, the plastic zone in excavated section penetrated into the heading face. When the seepage pressure increased to 0.5 MPa, the plastic zone in excavated section completely penetrated into the heading face. Consequently, the region of the tensile stress appeared below the heading face which accounted for more than 30% of the total area of the heading face. At this time, water flooding may occur through the heading face which will result in an integral-sliding failure if the seepage pressure is any greater.
4. Conclusions
In light of the frequent disasters including flooding, mudslides, and general leakage that can result from high water pressure during construction of tunnels and other underground structures in water-rich regions, a mathematical algorithm was adopted to model the problem. This algorithm consisted of a modified effective stress formula that was introduced into the CPPM together with an adopted integration algorithm of an elastic-plastic constitutive equation and elastic-plastic finite element program which can adopt water pore pressure. The resulting study showed that pore water pressure has a significant influence on stress, deformation, and the plastic zone of a rock mass that surrounds a tunnel. With an increase of pore pressure, changes in the plastic zone of the surrounding rock mass are more pronounced.
The developed program was used to study the effect of seepage pressure on stress, deformation, and the plastic zone of surrounding rock mass of the tunnel, which goes through a discordant contact zone at the tunnel entrance. Some of conclusions of this study were as follows:(1)Pore water pressure has a significant influence on stress, deformation, and the plastic zone of the surrounding rock mass of Fusong tunnel. With an increase of pore pressure, the areas of the plastic zone on tunnel heading face and completed section of the tunnel increased. Therefore, the influence of groundwater on tunnel stability needs to be considered when the tunnel is located at high water levels and passes through a permeable discordant contact zone.(2)Existence of groundwater flow during tunnel construction results in stress redistribution in the surrounding rock mass. When considering the influence of the seepage, the settlement displacement at the tunnel top and heaving displacement at the bottom caused by tunnel excavation will increase. It was also found that the horizontal and vertical displacements of the surrounding rock mass are distributed symmetrically about the center of the tunnel. Vertical displacement at the tunnel vault and bottom increased with a decrease in buried depth of the calculated nodes.(3)Comparison between displacements under circumstances of seepage and no seepage indicates that tensile stress appeared in the surrounding rock mass at the bottom of the tunnel due to seepage, which is detrimental to the stability of the tunnel. Therefore, tunnel sections in areas with high hydraulic head seepage require extra measures to enhance support and prevent collapse of the tunnel in the design phase.(4)Stability analysis of the tunnel during construction considering different groundwater levels reveals that, during the construction of the tunnel crossing a fault in rainy season, potential water flooding at heading face must be considered if there is a connection between water in discordant contact zone and surface water. So it is suggested to adopt full-face curtain grouting to ensure the safety of the heading face.
Competing Interests
The authors declare that there is no conflict of interests regarding the publication of this paper.
Acknowledgments
The authors deeply appreciate support from the National Natural Science Foundation (51678101 and 51578447).