Abstract
The undesirable effect on the stability for cross-river tunnel faces considering pore water pressure was observed with the consideration of the soil arch effect by using the discrete technology for the first time. In light of the upper bound of plastic theory, an improved failure mechanism of the deep-buried tunnel face was established. A new discrete technology approach taking account into the soil arching effect was proposed to estimate the stability for cross-river tunnel faces subjected to pore water pressure. The presented approach is validated by comparing with the existing solutions as well as showing great improvements. After verification, based on the failure mechanism, this paper discusses the impact of the changing water level and the soil parameters on the normalized supporting pressure and meanwhile analyzes the variation of the shape of collapsing domain of soils ahead of the tunnel face considering the soil arching effect. The results illustrate that soils with the bigger friction angle form the arch more easily during excavation, and with higher water height, the soil arching effect appears not as obvious as expected, particularly on those soils with the smaller friction angle.
1. Introduction
Due to the increasing development of the urbanization, the traffic congestion has become one of the reasons restricting the economic growth. Therefore, the tunnel is regarded as an effective way to alleviate this issue. Recently, as the massive utilization of the shallow-buried space, numerous tunnels are chosen to construct in the deep-buried space gradually. As the shield tunnel has advantages of short construction period and high security, it is employed popularly by many civil engineers [1–3]. Nevertheless, the soil arching phenomenon is beneficial to maintain the face stability during the excavation of deep-buried tunnels, and this phenomenon is investigated by many scholars [4, 5].
The soil arching effect could trace back to the trapdoor system experiment done by Terzaghi [6], who discovered that soils with vertical displacement would be resisted by the shearing action of surrounding soils, which could make earth pressure transmit around and lead to the stress redistribution phenomenon. Hence, this means that the total earth pressure on the face may be less than the gravity of the soils above the tunnel, which could contribute to the face stability of tunnels obviously. Afterwards, the studies about the soil arching effect are conducted by many researchers. Chen et al. [7] introduced a wedge limit equilibrium failure model of tunnel face and computed the lateral earth pressure coefficient to incorporate the presence of soil arching the stability assessment of tunnel faces in sand. By conducting the indoor experiment, Jacobsz [8] obtained different heights of soil arch for a series of conditions. In addition, He and Zhang [9] revealed the reasons for major principal stress rotation involved in the soil arch by means of a numerical approach called discontinuous deformation analysis (DDA). So far, a rotational kinematical failure mechanism proposed by Mollon et al. [10] is popular among scholars as it is validated to achieve great improvements compared to other traditional models. Based on this failure model, Zou et al. [11] calculated the lateral earth pressure coefficients of soils and compared the solutions derived by the upper-bound method and limit equilibrium theorem, respectively. Then, Chen et al. [12] introduced the influence of the anisotropic and nonhomogeneous soil properties into the framework of the upper-bound theory and proposed an advanced failure mechanism of tunnel faces.
However, the aforementioned works mainly focus on the face stability issue with the consideration of the soil arching effect in dry soils, while most urban shield tunnels often are excavated in the aquifer. It is necessary to consider the presence of pore water pressure in the stability analysis of tunnel faces. Currently, common methods to compute the pore pressure include the pore water coefficient method and the numerical method basically. The first approach is originally proposed by Bishop and Morgenstern [13] to evaluate the slope stability in saturated soils. They defined a coefficient to describe the seepage field, and this simplification has been validated by many scholars. However, Pan and Dias [14] stated that the comparison between results computed by this method with that derived by the numerical approach shows that outcomes obtained by this approach with reasonable coefficients are satisfactory, while it has the drawback of difficulty to know an appropriate value in advance, which could limit its application to a large extent. In fact, this assumption is not in line with the actual pore water pressure distribution. On the other hand, the numerical method denotes using numerical software to establish a model whose dimension is the same as the theoretical model. Then, by means of the hydraulic coupling calculation of this software, the pore pressure field ahead of the tunnel face could be obtained [14–16]. It should be noted that mesh density is required to satisfy enough accuracy. Hence, the disadvantage of this approach is low computation efficiency apparently.
In order to overcome these limitations, with the assistance of the upper-bound method, this paper establishes a rotational failure mechanism with the soil arching effect by using the discrete technology by Mollon et al. [10], to estimate the tunnel face stability under the impact of pore pressure. The finite difference software FLAC3D is adopted to obtain the numerical solutions of pore pressure, and the corresponding contour is presented as well. Then, we propose a linear interpolation approach to introduce the pore water pressure into the presented failure model to compute the limit supporting pressure. By comparing the presented solutions with the existing solutions, the proposed approach is verified and shows great improvements. At the end of this study, the effect of soil properties with changed water heights on the limit supporting pressure and collapsing domain of soils is also discussed.
2. An Improved Failure Mechanism of the Deep-Buried Tunnel
2.1. The Upper-Bound Limit Analysis Theory
The limit analysis approach consists of the upper-bound theorem and lower-bound theorem and is introduced firstly by Chen [17] into the geotechnics as well as used widely by many researchers to analyze the stability issue of geotechnical structures. Compared to the lower-bound theory, a kinematically admissible velocity field could be constructed more easily, and therefore, the upper-bound theory is more popular among scholars [18–20]. The upper-bound method states that theoretical collapsing load derived by equaling work rates done by external forces to the entire internal energy dissipation should not be less than the actual collapsing load. Therefore, the critical supporting pressure needed to maintain the face stability of tunnels could be computed after establishing a reasonable failure mechanism.
Before utilizing the upper-bound theory to conduct the stability issue, there are some involved assumptions to be emphasized: (1) the whole collapsing blocks are assumed to be rigid, i.e., no volumetric strain would occur; (2) the soil considered is thought as an ideal plastic material and should follow the associated flow rule, that is, the dilation angle of any point located on the slip surface equals to the internal friction angle; and (3) all materials analyzed comply with the Mohr-Coulomb (MC) failure criterion. Note that the soils involved in this analysis are sand and soil cohesion is therefore zero.
2.2. The Rotational Failure Mechanism in the Underlying Layer
In light of the kinematical approach, Mollon et al. [10] proposed a discrete technology to establish a rotational failure mechanism which shows notable advantages in the stability assessment of tunnel faces. Moreover, the principle of this failure mechanism is extended to the stability issue of other geotechnical structures by many researchers [21–23], and the results derived are satisfactory. To highlight the positive influence of the soil arch, this paper tends to combine it with the rotational failure mechanism to investigate the stability of deep-buried shield tunnel faces, and the proposed failure model is composed of the soil arching model in the upper layer and rotational failure mechanism constructed by the discrete method. The aim of this section is to elaborate the principle of generating the failure mechanism in the underlying layer.
As mentioned in the works of Mollon et al. [10], the whole potential collapsing domain ahead of the tunnel face rotates around the horizontal -axis through point O at the angular velocity , with representing the buried depth and representing the diameter of the tunnel, as shown in Figure 1. The dashed line refers to the original failure mechanism which should be replaced by the soil arching model in the analysis. Besides, the surface BM denotes the outer boundary of the rotational failure mechanism, and it should be the log spiral starting from point B at the tunnel invert to point M at the top of the tunnel. Owing to the postulation of the normal condition, the angle between the velocity of any point located in the surface BM and outer normal direction is . In the polar coordinate system, the log spiral could be expressed by the polar radius and rotation angle : in which is the distance from point O to point B, denotes the counterclockwise angle between OB and vertical direction, and represents the friction angle.
As seen in Figure 2(a), the rotational failure mechanism could be divided into two regions including part 1 and part 2 according to the generation rule. Firstly, the tunnel face is uniformly discretized into parts, and the whole model is also discretized by a set of planes through rotation point O. Specifically, any plane in part 1 could be determined by point O and all discrete points on the tunnel face, and in part 2, any plane could be generated by the former plane rotating around horizontal -axis through point O with an angle . The entire rotational failure mechanism could be generated by repeating the above process, which is called the “point by point” as well.
As sketched in Figure 3, a discrete point in the plane could be obtained by the given points and in the former plane based on the normality condition and then connect all adjacent points , , and to acquire the whole surface of the failure mechanism in the underlying layer. More details could be explored in Mollon et al. [10].
Obviously, the establishment of rotational failure mechanism needs to be terminated at the top plane AM of the tunnel, as shown in Figure 1. Hence, a linear interpolation technology is adopted herein to realize it, and the detailed process is that during the process of “point by point,” the location of point would be judged every moment; if the distance between point and the tunnel invert is higher than the tunnel height, the generation process would be closed and this point would be replaced by point obtained by using a linear interpolation technology, where the point is the intersection of the top surface of the tunnel and the line , as shown in Figure 2(b). The point is the midpoint of line . Based on this linear technology, the rotational failure mechanism derived could be seen in Figure 4.
2.3. The Soil Arching Effect in the Upper Layer
As mentioned before, the soil arching is great of interest for the stability evaluation of deep-buried tunnel faces during excavation. Therefore, this section introduced a theoretical model proposed by Wan et al. [4] to construct an improve failure mechanism considering the soil arching. As seen in Figure 5, due to the shearing resistance of adjacent motionless soils, the pressure of soils with relative vertical displacement would be transmitted around, which could lead to a stress redistribution phenomenon. Particularly, the domain with bigger displacement shapes like a rectangle while that with smaller displacement shapes like a curve. Wan et al. [4] stated that the soil arch consists of a rectangular arch and a parabolic arch, and the soils between both of them are in the potential failure state with the bottom of parabolic arch being the fracture surface. In the analysis, , , and are utilized herein to describe the heights of rectangular arch, potential failure domain, and parabolic arch, respectively, and and could written as follows: where is the width of rectangular arch and is the friction angle. Because of the stress redistribution induced by the soil arching effect, the forces on the top of the rectangular arch could be divided into the gravity of soils in the potential failure state and the residual vertical stress after stress transfer. Apparently, could be computed easily by the geometry relationship shown in Figure 5, i.e., , where represents the counterclockwise angle between the tangential line of point at the parabolic arch and horizontal direction and is taken here as [4]. Thereby, the gravity of the whole potential failure soils could be rewritten as . Notice that the dry unit weight of soils would be substituted by the saturated unit weight when soils computed are under groundwater.
On the other hand, to derive the vertical force applied on the top of the rectangular arch, a tiny unit parabolic arch element is taken in the analysis, with of thickness, as illustrated in Figure 6(a). As mentioned before, the vertical earth pressure transferred from the soils with vertical displacement to the surrounding soils could be computed by subtracting vertical stress by static earth pressure, namely, .
(a)
(b)
However, the trapdoor system tests of Terzaghi [6] discovered that the soil arching has less influence on the stress state of soils from 2 to 3 times the width of the trapdoor above the center line to the surface. Hence, the height of entire soil arch is taken as . In other words, the soils in this domain could balance themselves by forming the arch, and the height of the parabolic arch is therefore . Then, we could derive the equivalent uniform force acting on the tiny unit parabolic element: where denotes the uniform surcharge load on the ground surface. The horizontal stress on the sides of the element could be computed by dividing the horizontal supporting force on both sides of the structure by the element thickness , where could be derived based on the mechanical characteristic of the arch structure, i.e., . According to the Mohr stress circle of point sketched in Figure 6(b), we could obtain the relationship between and , namely, , and then could compute the vertical stress on the top of the rectangular arch by combining Eq. (3):
Consequently, the uniform vertical force could be derived:
After that, to obtain the equivalent earth pressure on the top of the rotational failure mechanism in the underlying layer, an arch element with of thickness is selected herein, as seen in Figure 7.
Since the soils analyzed comply with the Mohr-Coulomb failure criterion, one could obtain by the mechanical equilibrium conditions for the soil arch element in horizontal and vertical directions: in which and denote the cross-sectional area and perimeter of the tiny element, respectively; and represent the vertical and horizontal stress on the tiny element; and refers to the shear stress on the element. The conclusions of Terzaghi [6] illustrated that there may be a relationship between and , and the ratio of them could be expressed as a function associated with the soil shearing strength, i.e., , where is the lateral stress ratio. In this study, the lateral stress ratio could be calculated by with the inspiration of Wan et al. [4], where is and is the counterclockwise angle between major principal stress of point and horizontal direction, i.e., .
On the other hand, as the bottom section of the rectangular arch is the same as that of the rotational failure mechanism in the underlying layer although it is irregular, the circumference and area of this section could be calculated based on the “point by point” described in Section 2.2. As shown in Figure 8, the coordinates of all points in the plane are known after generating the proposed failure mechanism and could be written as , . Hence, the unit perimeter is , and the perimeter of the entire section is . Likewise, the unit area is , and the area of the whole section is . Finally, with the combination of Eq. (6) and boundary condition , the equivalent pressure in the underlying layer could be obtained as follows:
3. The Limit Supporting Pressure
After establishing the failure mechanism, the subsequent task is to calculate the external work rates and internal energy dissipation. In this study, the total external forces include pore water pressure, supporting pressure of tunnel faces, soil gravity, equivalent earth pressure on the rotational failure mechanism, and uniform surcharge load. Due to the postulate of rigid blocks, volumetric strain is not allowed to occur during plastic flow, and therefore, internal energy only dissipates along the failure surfaces.
3.1. Work Equation
As stated by the upper-bound approach, in a kinematically admissible failure mechanism, the relationship between the external work rates and the total internal energy dissipation is as follows: in which is the gravity rate of sand, is the supporting pressure rate, is the work rate of the equivalent earth pressure, and and are the pore water pressure rate and the total internal energy dissipation rate, respectively. The internal dissipation rate is indeed calculated as based on Mollon et al. [10], while it should equal to zero as the sand is considered in this paper; that is, the internal energy dissipation could not occur during collapsing. Therefore, Eq. (8) could be simplified as , and the expression of the limit supporting pressure could hence be derived as follows: where and denote the unit saturated weight of sand and unit water weight and and represent the nondimensional parameters, corresponding to the impact of the gravity of sand and the equivalent earth pressure, and both of them could be computed as where , , , and refer to the polar coordinates of center points located on the surfaces and ; and are the volume of the element corresponding to facets and , respectively; is the area of facet; and and are the polar coordinates of center point of the facet, as shown in Figure 9.
To derive the critical supporting pressure of the tunnel face, an algorithm called fminsearch embedded in the mathematical software MATLAB is employed in this analysis, with and being the input, where and are the polar coordinate of point E. By considering the fact that the supporting pressure is to resist the collapse, its value is regarded to be negative, and hence, the minimum of the supporting pressure searched is as the output. Note that the input parameters should follow Eq. (11). Consequently, the only unknown item involved in Eq. (9) is , which would be discussed in Section 3.2.
3.2. The Calculation of Pore Water Pressure Rate
Viratjandr and Michalowski [24] stated that the pore water pressure rate could be computed as the sum of that done by the seepage and buoyance: in which is the pore water pressure, denotes the volumetric strain increment of the soil skeleton, and represent the volume of the failure mechanism and area of failure surfaces, is the unit outer vector normal to the surface , and is the velocity in the kinematically admissible velocity field. As mentioned before, no volumetric strain requires Eq. (12) to be simplified as . Therefore, how to derive the distribution of pore water pressure is the core problem to be resolved in this section.
As it is a bit difficult to solve the three-dimensional Laplace equation basically, this paper employed a numerical software FLAC3D to obtain the pore water pressure in the steady seepage field by means of the hydromechanical coupling calculation. The corresponding numerical model is sketched in Figure 10, and only half the model is established because of the symmetry.
Furthermore, to eliminate the boundary effect, the length, width, and height are set as 80 m, 40 m, and 65 m, respectively. Besides, the grids around the tunnel face are densified to ensure the high accuracy. In terms of the boundary condition, the displacement of the bottom is prohibited, and the vertical displacement of other three sides is fixed. It should be emphasized that the pore water pressure at the tunnel face is fixed as zero due to the tunnel lining. The other used parameters are m, m, GPa, and . Based on the works of Pan and Dias [14], the results could be satisfactory if the soil permeability coefficient is not lower than . Figure 11 shows the pore water pressure contours with various water height in the vertical plane.
In a word, the pore pressure computed by FLAC3D could be exported into the proposed failure mechanism by using a linear interpolation method, and the processes are as follows: (1) the numerical model with given geometry parameters (e.g., tunnel diameter , buried depth , and water height ) is established by FLAC3D at first, and then, the steady seepage calculation is on to derive the needed pore pressure. (2) The pore pressure corresponding to the discrete point could be derived by the linear interpolation method.
Therefore, by substituting the derived into Eq. (12), the nondimensional parameter could be computed as Eq. (13) with combination of Eqs. (8) and (9). Finally, we could obtain the critical supporting pressure by taking Eq. (13) into Eq. (9).
4. Comparisons
To validate the presented failure model of the deep-buried tunnel face considering soil arching effect, this section firstly compares the shapes of the theoretical model with the experimental results of Chambon and Corte [25] with different ratios of the cover depth to the tunnel diameter. Note that in the experiment tests of Chambon and Corte [25], the soil friction angle, soil cohesion, and sand density vary from 38° to 42°, 0 kPa to 5 kPa, and 15.3 kN/m3 to 16.1 kN/m3 , respectively, due to the inherent uncertainty during the experiment operation. To facilitate the comparison, the means of these sand parameters are selected as the input to derive the solutions of the proposed method, that is, °, kPa, and kN/m3.
Figure 12 shows the comparison between the presented failure mechanisms and the collapse domain of soils ahead the tunnel face by the experiment of Chambon and Corte [25] under various conditions where , 1, and 2. It could be found that our theoretical failure model agrees well with the experimental result while the failure region seems not to extend around the bottom of the tunnel if the cover depth is relatively small. In other words, the proposed method could describe well the collapse domain of deep-buried tunnel face.
Apart from that, to verify the presented failure model and also illustrate the advantage further, the presented normalized limit supporting pressure with changed water height is compared with the existing solutions, as shown in Figure 13. The used analytical parameters are as follows: , , kN/m3, , m, and m.
As shown in Figure 13, the presented results are always higher than the other solutions of Zou and Qian [15], Pan and Dias [14], and Lee et al. [26], and differences between them are rising with the increasing water height. The maximum differences are found in the case of , namely, 2.38%, 20.61%, and 38.51%, respectively, which indicates that the pore water pressure has a potential influence on the soil arch. The explanation for this may be that the pore pressure could reduce the matric suction between soil particles. Hence, it is crucial to incorporate the soil arching effect into the stability assessment of deep-buried shield tunnels especially for those cross-river tunnels excavated in the aquifer. On the other hand, this section verified that the presented failure mechanism could be effective for the stability issue of tunnel faces.
5. Parameter Analyses
Based on the verification in Section 4, this section is aimed at extending the proposed theoretical model to the practical engineering by analyzing the effect of different friction angles of sand on the normalized limit supporting pressure with changed water height. Also, the discussion about the impact of soil properties on the shape soil arch is our focus as well. The relative parameters are as follows: , , m, and m.
As seen in Figure 14, with the rise of the water height, the normalized supporting pressure rises linearly approximately, which means that the limit value of supporting pressure in a specific case could be computed based on two given cases by the linear interpolation, to ensure the construction safety. Additionally, it could be observed that, for some soils with the smaller friction angle, the variation of is more sensitive to the change of water height, and it could be attributed to the fact that these soils form the arch with difficulty, being adverse to the tunnel face stability. By contrast, the advantage that soils with bigger friction angles could form the arch more easily could weaken the undesirable impact of the pore pressure on the tunnel face stability.
Figure 15 describes the different shapes of failure mechanisms with various water height. It could be seen that the height of “horn” in the top of the model would become lower increasingly with the rise of water height, and the end of the rotational block looks more and more gentle. In addition, the increasing water height also leads to the wider collapsing domain, which could be explained that the larger pore pressure leads to the fast collapse and destroys the soil arch to some extent.
6. Conclusions
With the assistance of the upper-bound method, by utilizing the discrete technology of Mollon et al. [10], we introduce an advanced failure mechanism with the soil arching effect. In order to assess the stability of tunnel faces under the impact of pore pressure, the work rate of pore pressure is derived by incorporating the results of numerical software into this failure mechanism. Afterwards, comparisons between the existing works verify the presented approach. In the end of this paper, the effect of the changing water heights on the normalized limit supporting pressure is discussed. The following conclusions are obtained: (1)For those soils with the bigger friction angle, soils ahead the tunnel face form the arch more easily during excavation. However, soils with smaller shear strengths are suffered from the pore water pressure heavily; that is, the variation of normalized supporting pressure is more obvious. Therefore, civil engineers should seriously concern about the surrounding distribution of pore water pressure to avoid the collapse(2)During excavation in some layers with higher water height, the phenomenon of soil arch is not as obvious as expected, particularly on those soils with smaller friction angle. In other words, the pore water pressure would lead to the wider and deeper failure domain of soils. The possible reason for this phenomenon could be that the rising pore water pressure could reduce the matric suction of soils
Data Availability
The data used to support the findings of this study are available from the corresponding author upon request.
Conflicts of Interest
The authors declare that they have no conflicts of interest.
Acknowledgments
This research was funded by the Key Engineering Science and Technology Projects of Jiangxi Provincial Transportation Department (grant number 2021C0002).