Stability of twin circular tunnels in cohesive-frictional soil using the node-based smoothed finite element method ( NS-FEM )

This paper presents an upper bound limit analysis procedure using the node-based smoothed finite element method (NS-FEM) and second order cone programming (SOCP) to evaluate the stability of twin circular tunnels in cohesive-frictional soils subjected to surcharge loading. At first stage, kinematically admissible displacement fields of the tunnel problems are approximated by NS-FEM using triangular elements (NS-FEM-T3). Next, commercial software Mosek is employed to deal with the optimization problems, which are formulated as second order cone. Collapse loads as well as failure mechanisms of plane strain tunnels are obtained directly by solving the optimization problems. For twin circular tunnels, the distance between centers of two parallel tunnels is the major parameter used to determine the stability. In this study, the effects of mechanical soil properties and the ratio of tunnel diameter and the depth to the tunnel stability are investigated. Numerical results are verified with those available to demonstrate the accuracy of the proposed method.


Introduction
The rapid development of transportation in large cities requires the use of underground areas to construct infrastructures and facilities. In practice, two parallel circular tunnels are widely established to make separate way for each road or railway. In order to reduce construction costs, it is essential to give an accurate estimation of the stability of dual circular tunnels in the design of shallow tunnels. This issue has drawn a profound concern of many researchers, especially Atkinson and Cairncross [1], Cairncross [2], Seneviratne [3] and Mair [4]. In 1980, Davis et al. [5] proposed a theoretical solution for single circular tunnel in cohesive material under undrained condition using both upper bound and lower bound limit analysis. In the study of Mühlhaus [6], lower bound approach was applied to evaluate the stability of tunnels subject to uniform internal pressure. By using limit analysis, Leca and Dormieux [7] presented a more general solution for the face stability of shallow circular tunnels in frictional material.
In recent decades, the standard finite element method (FEM) has been rapidly developed to solve complicated engineering problems. A finite element procedure for linear analysis was first given by Sloan and Assadi [8] to evaluate the undrained stability of a square tunnel in a soil whose shear strength increases linearly with depth. Lyamin and Sloan [9], Lyamin et al. [10] and Yamamoto et al. [11,12] developed FEM-based nonlinear analysis methods to calculate the failure mechanisms of circular and square tunnels in cohesive-frictional soils. However, most of the theoretical studies were only focused on stability of single tunnel. There is a scarcity of evidences in literature which were suggested to determine the stability of twin tunnels. To investigate the stability of two parallel tunnels, a series of centrifuge model tests for clays was conducted by Wu and Lee [14], and their studies were aimed to formulate ground movements and failure mechanisms of soils surrounding two tunnels. Chahade and Shahrour [15] used Plaxis software to simulate the construction procedure of twin tunnels with various distance between

Problem definition
The twin circular tunnels which have diameter , depth and distance between the tunnels are illustrated as in Fig. 1. The ground deformation takes place under plane strain. The soil is assumed to be rigid perfectly plastic and modelled by a Mohr-Coulomb yield criterion with cohesion ′, friction angle ′ and unit weight . Drained loading conditions are also considered, and surcharge loading is applied to the ground surface. In order to assess the stability of the tunnel, a dimensionless load factor ⁄ is defined by using a function of ′, ⁄ ′, ⁄ and ⁄ , as shown in the following equation: In order to investigate the stability load factor ⁄ , the variation in considered parameters are ⁄ = 1-5, = 0°-20°, ⁄ = 0-3 and tunnel spacing ⁄ = 1.25-12. 5. To describe the smooth interface condition between the loading and the soil, the values of vertical velocities are equal to zero along the ground surface.

Upper bound limit analysis formulation
Let us consider a rigid-perfectly plastic body of area Ω ∈ with boundary Γ, which is subjected to body forces and to surface tractions g on the free portion Γ of Γ. The constrained boundary Γ is fixed and Γ ∪ Γ = Γ, Γ ∩ Γ = ∅. Let = [ ] belongs to a space of kinematically admissible velocity fields, where , are the velocity components in -anddirections, respectively. The strain rates can be expressed by relations: 523 The external work rate associated with a virtual plastic flow is given by: The internal plastic dissipation of the two-dimensional domain Ω can be written as: The upper bound theorem states that there exists a kinematically admissible displacement field ∈ , such that: where is the axact limit load multiplier of the load , and ( ) is the work of additional , not subjected to the multiplier. Letting = { ∈ | ( ) = 1}, the collapse load multiplier can be determined by the following formulae:

Brief on the node-based smoothed finite element method (NS-FEM)
The idea of the nodal integration in meshfree method using the strain smoothing technique was proposed by Chen and co-workers [24,25]. Later, Liu and collaborators [26-31] applied this technique to standard FEM to provide a softening effect to improve the solution of FEM, which is called smoothed finite element method (S-FEM), such as ES-FEM, CS-FEM, NS-FEM, and so on. The main difference between S-FEM and the standard FEM is the strain field. In standard FEM, the displacement field is assumed, and the strain is calculated from the strain-displacement relation. In S-FEM, a strain smoothing is calculated from the strain in FEM by a smoothed function. In this paper, the formulation of smoothing function using NS-FEM is presented. The problem domain Ω is divided into smoothing cells Ω ( ) associated with the node such that Ω = ∑ Ω ( ) and Ω ⋂Ω = ∅, ≠ and is the total number of field nodes located in the entire problem domain. Each triangular element will be divided into three quadrilateral sub-domain and each quadrilateral sub-domain is attached with the nearest field node.
A strain smoothing formulation on the cell Ω ( ) is now defined by the following operation: where Φ ( ) is a smoothing function that satisfies positive and normalized to unity: The smoothing function Φ ( ) is assumed constant: where ( ) = Ω ( ) is the area of the cell Ω ( ) and the smoothed strain on the domain Ω ( ) can be expressed as: where Γ ( ) is the boundary of the domain Ω ( ) as shown in Fig. 2 and ( ) is a matrix with components of the outward normal vector on the boundary Γ ( ) given by: The smoothed strain on the cell Ω ( ) associated with node can be calculated by: where ( ) is the set containing nodes that are directly connected to node , is the nodal displacement vector and the smoothed strain gradient matrix ( ) on the domain Ω ( ) can be determined from: where: When a linear compatible displacement field along the boundary Γ ( ) is used, one Gauss point is sufficient for line integration along each segment of boundary Γ ( ) of Ω ( ) , the above equation can be determined from: where is the total number of the boundary segment of Γ ( ) , is the Gauss point of the boundary segment of Γ ( ) which has length ( ) and outward unit normal ( ) .

NS-FEM formulation for plane strain with Mohr-Coulomb yield criterion
The soil is assumed to be perfectly plastic, and it obeys the Mohr-Coulomb failure criterion and associated flow rule. The Mohr-Coulomb yield function can be expressed in the form of stress components as: For an associated flow rule, the direction of the plastic strain rates vector is given by the gradient to the yield function, with its magnitude given by the plastic multiplier rate : Therefore, the power of dissipation can be formulated as a function of strain rates for each domain is presented in [41]: where is the area of the element of node , is a vector of additional variables defined by: The change volume after deformation in cohesive-frictional soil can be calculated from: Introducing an approximation of the displacement and the smoothed strains rates ̃ can be calculated from Eq. (13), the upper-bound limit analysis problem for plane strain using NS-FEM can be determined by minimizing the objective function: = 0, on Γ , ( ) = 1, ̃ + ̃ = sin , ‖ ‖ ≤ , = 1, 2, … , , where is the total number of nodes in domain. The fourth constraint in Eq. (24) is a form of quadratic cones.

Numerical results
Due to symmetry, only half of the problem is considered. In this paper, GiD [42] is employed for automatic mesh generation with three node triangular elements and adaptive refinement along the periphery of the tunnel. The size of domain is chosen sufficiently large enough to ensure that the failure mechanism only taking place inside the considered domain. For the case of ⁄ = 1 and ⁄ = 2, the typical finite element meshes of 1672 triangular elements are employed in numerical analysis as shown in Fig. 3. Figs. 6-9 illustrate the failure mechanisms and velocity plots for various cases of twin tunnels. Fig. 6(a) shows the power dissipation of twin circular tunnels for shallow tunnel in the case that friction angle ′ and spacing ratio ⁄ are relatively small. It is noticeable that a small slip failure occurs between twin tunnels and a large failure domain extends up to the ground surface. Fig. 7(a) shows the case for shallow depth, moderate friction angle ′ and small distance between dual tunnels ⁄ . In this figure, a small slip surface between two circular tunnels enlarges to the top and bottom of tunnels and a large surface from the middle part of the tunnel extends up to the ground surface. The power dissipations obtained in this study are quite identical to those which were derived from rigid block technique proposed by Chen [43] and solution of Yamamoto et al. In Figs. 6-9, the stability numbers of twin circular tunnels increase when increasing the value of the ratio ⁄ . When ⁄ is small enough for two tunnels to interact, the failure mechanism enlarges as ⁄ and ⁄ increase. When the distance between two tunnels ⁄ is large enough, there will be no influence on the failure mechanism of each tunnel. Therefore, the ratio ⁄ parameter plays an important role in the behaviour of the failure mechanism and the increase of stability number is due to the effects of interaction. The power dissipations in Fig. 10 show no interaction between dual circular tunnels, and the failure mechanism of two parallel tunnels looks like the case of individual single tunnel. For example, Figs. 10(a)-(c) show that the power dissipations in case the ratio ⁄ = 1, 3,5 and ⁄ = 1, = 10°, there is no interaction of twin circular tunnels when ⁄ = 3.5, 7, 10, respectively. To maximize stability number, the twin circular tunnels should be placed far enough apart to ensure no interaction between the tunnels and that stability number is equal to the single circular tunnel.
To evaluate the accuracy of this research, the obtained results and those of Yamamoto et al [17] were plotted in Figs. 11-13. These figures show that the results derived from this proposed method are very good because most of them are lower than upper bound solutions of Yamamoto et al. [17] and higher than Yamamoto's lower bound solution. As shown in Figs. 11 and 12, stability numbers decrease when increasing the value of the ratio ⁄ , its mean that the weight of the soil effects to the failure mechanism of twin circular tunnels. Fig. 13 presents the results of stability numbers for deep tunnels ⁄ = 5. In this figure, when the ratio ⁄ increases from 1.25 to 3, the stability numbers slightly decrease. This is because of the fact that when the two tunnels are very close together, the extra resistance gained by increasing the width of the pillar is not enough to be the counterbalance of the extra soil mass that it must support. Increasing ⁄ further, the stability numbers tend to become constant when the space between the tunnels exceeds a certain value. At this point, the failure mechanisms become two individual single tunnels. The values of stability numbers obtained by using NS-FEM and SOCP are summarized in Table A1 (in Appendix). As expected, when the space between twin tunnels exceeds a certain value, the stability load factor tends to become constant. The negative results imply that a tensile normal stress can be applied to the ground surface to ensure that there is no collapse occurred, but this cannot be seen in engineering practice. The positive one means that the tunnel will be collapsed when it is subjected a compressive stress on the ground surface as this value. The stability numbers at the nointeraction points for twin circular tunnels are highlighted in bold. When the spacing between the tunnels exceeds these points, the values obtained from NS-FEM tend to become constant. In cases of ⁄ = 5 and ⁄ = 3, the stability numbers that approximate zero are indicated by "-", it means that the tunnels collapse under the weight of the soil.  Table 1. In comparison with those of Yamamoto et al. [17], it is recognized that the numerical procedure using NS-FEM and SOCP not only reduces an appreciable amount of variables in optimization problem, but also gives a very better solution.

Conclusions
In this paper, a numerical procedure based on the node-based smoothed finite element method (NS-FEM) is proposed to evaluate the stability of a plane strain twin circular tunnels in cohesive-frictional soil subjected to surcharge loading. The obtained results are in well agreement with the average values of the lower bound and upper bound reported by Yamamoto et al. [17]. The ratio ⁄ plays an important role in the failure mechanism of twin circular tunnels. To maximize stability number, the distance between the twin circular tunnels should be large enough to ensure no interaction between them and the stability number approaches the value of the single circular tunnel. Various numerical examples for twin circular tunnel problems have been carried out showing that the presented method is able to provide accurate and stable solutions with minimal computational effort. It is promising to develop the proposed method for more complex and large scale problems. [34] Le C. V., Nguyen-Xuan H., Askes H., Bordas S., Rabczuk T., Nguyen-Vinh H. A cell-based smoothed finite element method for kinematic limit analysis.