Stability of twin circular tunnels in cohesivefrictional soil using the nodebased smoothed finite element method (NSFEM)
Thien M. Vo^{1} , Tam M. Nguyen^{2} , An N. Chau^{3} , Hoang C. Nguyen^{4}
^{1, 2, 3}Department of Civil Engineering, Ho Chi Minh City University of Technology, Ho Chi Minh City, Vietnam
^{1, 4}Department of Civil Engineering, HCMC University of Technology, Ho Chi Minh City, Vietnam
^{1}Corresponding author
Journal of Vibroengineering, Vol. 19, Issue 1, 2017, p. 520538.
https://doi.org/10.21595/jve.2016.17896
Received 10 September 2016; received in revised form 5 November 2016; accepted 23 November 2016; published 15 February 2017
This paper presents an upper bound limit analysis procedure using the nodebased smoothed finite element method (NSFEM) and second order cone programming (SOCP) to evaluate the stability of twin circular tunnels in cohesivefrictional soils subjected to surcharge loading. At first stage, kinematically admissible displacement fields of the tunnel problems are approximated by NSFEM using triangular elements (NSFEMT3). 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.
Keywords: limit analysis, dual circular tunnels, stability, NSFEM, SOCP.
1. 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 FEMbased nonlinear analysis methods to calculate the failure mechanisms of circular and square tunnels in cohesivefrictional 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 their centers and then considered its influence on tunnels stability. The upper bound limit analysis was also applied in Osman’s study [16] to evaluate the stability of twin circular tunnels for largescale engineering problems. Recently, Yamamoto et al. [17] proposed an efficient method to assess the stability of dual circular tunnels in cohesive material, however, the nonlinear optimization procedure required large computational efforts.
Due to the simplicity, the standard finite element method using triangular element (FEMT3) are popularly. One of the marked drawbacks of FEMT3 elements is volumetric locking phenomenon, which is often occurred in the nearly incompressible materials. To overcome this, many methods were suggested as reduce integration methods [18], Bbar methods [19], enhanced assumed strain [2022], average nodal technique [23] and so on. In this study, NSFEM that is original from integration methods is used to overcome the challenge of volumetric locking. More precisely, the idea of the nodal integration in meshfree method using the strain smoothing technique was proposed by Chen and his collaborators [24, 25]. Then, Liu GR. and NguyenThoi T. [26] applied this technique to standard FEM to provide a softening effect to improve the solution of FEM, which is called smoothed finite element method (SFEM), including cellbased smoothed FEM (CSFEM) [27], nodebased smoothed FEM (NSFEM) [28], facebased smoothed FEM (FSFEM) [29], and edgebased smoothed FEM (ESFEM) [30]. In these SFEM models, the finite element mesh is used similarly as in FEM models. However, these SFEM models evaluate the weak form based on smoothing domains created from the entities of the element mesh such as cells/elements, or nodes, or faces, or edges. These smoothing domains can be located inside the elements (CSFEM), or cover parts of adjacent elements (NSFEM, FSFEM and ESFEM). These smoothing domains are linear independent and ensure stability and convergence of the SFEM models. The theoretical aspect of the SFEM is clearly presented in [31, 32]. Several further developments of SFEMs for limit and shakedown analysis have been investigated in [3337].
Recently, the nodebased smoothed finite element method (NSFEM) has been employed for upper and lower bound limit problems due to following advantages: (i) total degreesoffreedom significantly decreased, leading to a fast convergence for solutions, (ii) volumetric locking phenomenon is prevented by using NSFEM method in solving undrained (imcompressible) problems because elements do not have enough necessary degreesof freedom to find solutions with the condition of constant volume, (iii) by using of smoothed strains in NSFEM, the integration is conducted in the edges of smoothed cells, as a results, there is no need to compute first derivation of shape functions [28].
In upper bound limit analysis, the internal plastic dissipation is minimized to determine the ultimate load bearing capacity of structures. The yield criterion can be formed in a secondorder cone programming SOCP. To solve the resulting conic problems, the MATLAB (version 7.8.0) and the Mosek (version 6.0) [38] are used to give all solutions in this paper. The Mosek optimization toolbox can solve only convex optimization problems such as linear, quadratic and conic programming. Largescale SOCP problems can be solved effectively using primaldual algorithms based on the interiorpoint method. This algorithm was proved to be effectively optimization technique for limit analysis of structures.
Yamamoto et al. [17] investigated the upper and lower bound limit analysis for clay using FEMT3 and nonlinear optimization in 2012. The upper bound limit analysis using FEMT3 and linear optimization was presented by Jagdish Prasad Sahoo and Jyant Kumar [39] in 2013. Recently, Wilson et al. [40] studied twintunnels in the undrained condition and Tresca yield criterion by using FEMT3 and secondorder cone optimization programming (SOCP). The aim of this research is to present a numerical procedure using the nodebased smoothed finite element method (NSFEM) and secondorder cone programming (SOCP) to find the collapse load as well as failure mechanism of twin circular tunnels in cohesivefrictional soil subjected to surcharge loading. To evaluate the accuracy of this suggested procedure, the obtained results are compared with those of Yamamoto et al. [17].
This paper is outlined as follows. In Section 2, the problem definition is described. The upper bound limit analysis formulation is presented in Section 3. In Section 4, the brief on the nodebased smoothed finite element method NSFEM is introduced. NSFEM formulation for plane strain with MohrCoulomb yield criterion is presented in Section 5. In Section 6, some numerical examples are performed and discussed to demonstrate the effectiveness of presented method. Some concluding remarks are made in Section 7.
2. Problem definition
The twin circular tunnels which have diameter $D$, depth $H$ and distance $S$ 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 MohrCoulomb yield criterion with cohesion $c\text{'}$, friction angle $\varphi \text{'}$ and unit weight $\gamma $. 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 ${\sigma}_{s}/{c}^{\text{'}}$ is defined by using a function of $\varphi \text{'}$, $\gamma D/c\text{'}$, $S/D$ and $H/D$, as shown in the following equation:
In order to investigate the stability load factor ${\sigma}_{s}/{c}^{\text{'}}$, the variation in considered parameters are $H/D=$ 15, ${\varphi}^{\text{'}}=$ 0°20°, ${\gamma D/c}^{\text{'}}=$03 and tunnel spacing $S/D=$ 1.2512.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.
Fig. 1. Twin circular tunnels subjected to surcharge loading
3. Upper bound limit analysis formulation
Let us consider a rigidperfectly plastic body of area $\mathrm{\Omega}\in {R}^{2}$ with boundary $\mathrm{\Gamma}$, which is subjected to body forces $f$ and to surface tractions g on the free portion ${\mathrm{\Gamma}}_{t}$ of $\mathrm{\Gamma}$. The constrained boundary ${\mathrm{\Gamma}}_{u}$ is fixed and ${\mathrm{\Gamma}}_{u}\cup {\mathrm{\Gamma}}_{t}=\mathrm{\Gamma}$, ${\mathrm{\Gamma}}_{u}\cap {\mathrm{\Gamma}}_{t}=\mathrm{\varnothing}$. Let $\dot{\mathbf{u}}=[\begin{array}{ll}\dot{u}& \dot{v}{]}^{T}\end{array}$ belongs to a space $U$ of kinematically admissible velocity fields, where $\dot{u}$, $\dot{v}$are the velocity components in $x$ and $y$ directions, respectively. The strain rates$\dot{\epsilon}$can be expressed by relations:
The external work rate associated with a virtual plastic flow $\dot{u}$ is given by:
The internal plastic dissipation of the twodimensional domain $\mathrm{\Omega}$ can be written as:
The upper bound theorem states that there exists a kinematically admissible displacement field $\dot{\mathbf{u}}\in U$, such that:
where ${\lambda}^{+}$is the axact limit load multiplier of the load $g$, $f$ and ${W}_{ext}^{0}\left(\dot{\mathbf{u}}\right)$ is the work of additional ${g}_{0}$, ${t}_{0}$ not subjected to the multiplier.
Letting $C=\{\dot{\mathbf{u}}\in U{W}_{ext}\left(\dot{\mathbf{u}}\right)=1\},$ the collapse load multiplier ${\lambda}^{+}$^{}can be determined by the following formulae:
4. Brief on the nodebased smoothed finite element method (NSFEM)
The idea of the nodal integration in meshfree method using the strain smoothing technique was proposed by Chen and coworkers [24, 25]. Later, Liu and collaborators [2631] applied this technique to standard FEM to provide a softening effect to improve the solution of FEM, which is called smoothed finite element method (SFEM), such as ESFEM, CSFEM, NSFEM, and so on. The main difference between SFEM and the standard FEM is the strain field. In standard FEM, the displacement field is assumed, and the strain is calculated from the straindisplacement relation. In SFEM, a strain smoothing is calculated from the strain in FEM by a smoothed function. In this paper, the formulation of smoothing function using NSFEM is presented.
Fig. 2. Triangular elements and smoothing cells associated with the nodes in the NSFEM
In NSFEM, the integration based on nodes and strain smoothing technique [24, 25] is used. The problem domain $\mathrm{\Omega}$ is divided into ${N}_{n}$ smoothing cells ${\mathrm{\Omega}}^{\left(k\right)}$ associated with the node $k$ such that $\mathrm{\Omega}=\sum _{k=1}^{{N}_{n}}{\mathrm{\Omega}}^{\left(k\right)}$ and ${\mathrm{\Omega}}^{i}\bigcap {\mathrm{\Omega}}^{j}=\varnothing $, $i\ne j$ and ${N}_{n}$ is the total number of field nodes located in the entire problem domain. Each triangular element will be divided into three quadrilateral subdomain and each quadrilateral subdomain is attached with the nearest field node.
A strain smoothing formulation on the cell ${\mathrm{\Omega}}^{\left(k\right)}$ is now defined by the following operation:
where ${\mathrm{\Phi}}_{k}\left(x\right)$ is a smoothing function that satisfies positive and normalized to unity:
The smoothing function ${\mathrm{\Phi}}_{k}\left(x\right)$ is assumed constant:
where ${A}^{\left(k\right)}=\underset{{\mathrm{\Omega}}^{\left(k\right)}}{\int}d\mathrm{\Omega}$ is the area of the cell ${\mathrm{\Omega}}^{\left(k\right)}$ and the smoothed strain on the domain ${\mathrm{\Omega}}^{\left(k\right)}$ can be expressed as:
where ${\mathrm{\Gamma}}^{\left(k\right)}$ is the boundary of the domain ${\mathrm{\Omega}}^{\left(k\right)}$ as shown in Fig. 2 and ${\mathbf{n}}^{\left(k\right)}$ is a matrix with components of the outward normal vector on the boundary ${\mathrm{\Gamma}}^{\left(k\right)}$ given by:
The smoothed strain on the cell ${\mathrm{\Omega}}^{\left(k\right)}$ associated with node $k$ can be calculated by:
where ${N}^{\left(k\right)}$ is the set containing nodes that are directly connected to node $k$, ${\mathbf{d}}_{I}$ is the nodal displacement vector and the smoothed strain gradient matrix ${\stackrel{~}{\mathbf{B}}}_{I}\left({\mathbf{x}}_{k}\right)$ on the domain ${\mathrm{\Omega}}^{\left(k\right)}$ can be determined from:
where:
When a linear compatible displacement field along the boundary ${\mathrm{\Gamma}}^{\left(k\right)}$ is used, one Gauss point is sufficient for line integration along each segment of boundary ${\mathrm{\Gamma}}^{\left(k\right)}$ of ${\mathrm{\Omega}}^{\left(k\right)}$, the above equation can be determined from:
where $M$ is the total number of the boundary segment of ${\mathrm{\Gamma}}_{i}^{\left(k\right)}$, ${x}_{i}^{GP}$ is the Gauss point of the boundary segment of ${\mathrm{\Gamma}}_{i}^{\left(k\right)}$ which has length ${l}_{i}^{\left(k\right)}$ and outward unit normal ${n}_{ih}^{\left(k\right)}$.
5. NSFEM formulation for plane strain with MohrCoulomb yield criterion
The soil is assumed to be perfectly plastic, and it obeys the MohrCoulomb failure criterion and associated flow rule. The MohrCoulomb 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 $\dot{\mu}$:
Therefore, the power of dissipation can be formulated as a function of strain rates for each domain is presented in [41]:
where ${A}_{i}$ is the area of the element of node $i$, ${t}_{i}$ is a vector of additional variables defined by:
The change volume after deformation in cohesivefrictional soil can be calculated from:
Introducing an approximation of the displacement and the smoothed strains rates ${\dot{\stackrel{~}{\epsilon}}}_{i}$ can be calculated from Eq. (13), the upperbound limit analysis problem for plane strain using NSFEM can be determined by minimizing the objective function:
where ${N}_{n}$ is the total number of nodes in domain.
The fourth constraint in Eq. (24) is a form of quadratic cones.
6. 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 $H/D=$ 1 and $S/D=$ 2, the typical finite element meshes of 1672 triangular elements are employed in numerical analysis as shown in Fig. 3.
Fig. 3. Typical NSFEM meshes for twin circular tunnels ($H/D=$ 1, $S/D=$2)
Fig. 4. Deformed meshes ($H/D=$ 1, $S/D=$ 2)
Fig. 5. Velocity plot ($H/D=$ 1, $S/D=$ 2)
Fig. 6. The comparison between rigidblock mechanism and failure mechanism obtained by this method NSFEM. For the case: $H/D=$ 1, $\gamma D/{c}^{\text{'}}=$ 1, $S/D=$ 1.5, ${\varphi}^{\text{'}}=$ 10, smooth interface
a) Power dissipation ${\sigma}_{s}/c\text{'}=$0.74
b) Velocity plot
c) Rigidblock mechanism ${\sigma}_{s}/c\text{'}=$0.83
Figs. 69 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 $\varphi \text{'}$ and spacing ratio $S/D$ 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 $\varphi \text{'}$ and small distance between dual tunnels $S/D$. 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. [17].
Fig. 7. The comparison between rigidblock mechanism and failure mechanism obtained by this method NSFEM. For the case: $H/D=$ 1, $\gamma D/{c}^{\text{'}}=$ 1, $S/D=$2, ${\varphi}^{\text{'}}=$ 20°, smooth interface
a) Power dissipation ${\sigma}_{s}/c\text{'}=$2.45
b) Velocity plot
c) Rigidblock mechanism ${\sigma}_{s}/c\text{'}=$ 3.89
Fig. 8. The comparison between rigidblock mechanism and failure mechanism obtained by this method NSFEM. For the case: $H/D=$ 3, $\gamma D/{c}^{\text{'}}=$ 1, $S/D=$ 2, ${\varphi}^{\text{'}}=$10°, smooth interface
a) Power dissipation ${\sigma}_{s}/{c}^{\text{'}}=$ 0.85
b) Velocity plot
c) Rigidblock mechanism ${\sigma}_{s}/c\text{'}=$1.39
Figs. 8, 9 show the failure mechanism of twin circular tunnels for moderate depth, small friction angles $\varphi \text{'}$ medium spacing between the tunnels $S/D$. The slip surface between the tunnels enlarges to the top and bottom of tunnels, a large surface originates the bottom of tunnel extends up to the ground surface. The values of stability number obtained from rigidblock mechanism are slightly greater than those from NSFEM upper bound solution. The errors stability numbers calculated from NSFEM limit analysis and Yamamoto et al. [17] in the cases are shown in Figs. 69 are 2.6 %, 0.4 %, 4.4 % and 3.7 %, respectively.
In Figs. 69, the stability numbers of twin circular tunnels increase when increasing the value of the ratio $S/D$. When $S/D$ is small enough for two tunnels to interact, the failure mechanism enlarges as $H/D$ and $\gamma D/{c}^{\text{'}}$ increase. When the distance between two tunnels $S/D$ is large enough, there will be no influence on the failure mechanism of each tunnel. Therefore, the ratio $S/D$ 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.
Fig. 9. The comparison between rigidblock mechanism and failure mechanism obtained by this method NSFEM. For the case: $H/D=$ 3, $\gamma D/{c}^{\text{'}}=$1, $S/D=$ 3.5, ${\varphi}^{\text{'}}=$10°, smooth interface
a) Power dissipation ${\sigma}_{s}/c\text{'}=$1.56
b) Velocity plot
c) Rigidblock mechanism ${\sigma}_{s}/c\text{'}=$3.02
Fig. 10. Numerical results from NSFEM limit analysis (smooth interface)
a) Power dissipation: ($H/D=$ 1, $\gamma D/{c}^{\text{'}}=$ 1, $S/D=$ 3.5, ${\varphi}^{\text{'}}=$10°)
b) Power dissipation: ($H/D=$ 3, $\gamma D/{c}^{\text{'}}=$ 1, $S/D=$7, ${\varphi}^{\text{'}}=$10°
c) Power dissipation: ($H/D=$ 5, $\gamma D/{c}^{\text{'}}=$ 1, $S/D=$10, ${\varphi}^{\text{'}}=$10°
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 $H/D=$ 1, 3,5 and $\gamma D/{c}^{\text{'}}=$ 1, ${\varphi}^{\text{'}}=$10°, there is no interaction of twin circular tunnels when $H/D=$ 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. 1113. 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 $\gamma D/{c}^{\text{'}}$, 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 $H/D=$ 5. In this figure, when the ratio $S/D$ 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 $S/D$ 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.
Fig. 11. Comparisons of the stability numbers between present method and Yamamoto et al. [17] For the case $H/D=$ 1, smooth interface
a)$\varphi \text{'}=$ 5°
b)$\varphi \text{'}=$ 10°
c)$\varphi \text{'}=$ 15°
d)$\varphi \text{'}=$ 20°
The values of stability numbers obtained by using NSFEM 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 NSFEM tend to become constant. In cases of $H/D=$ 5 and $\gamma D/{c}^{\text{'}}=$ 3, the stability numbers that approximate zero are indicated by “–“, it means that the tunnels collapse under the weight of the soil.
Fig. 12. Comparisons of the stability numbers between present method and Yamamoto et al. [17]. For the case $H/D=$3, smooth interface
a)$\varphi \text{'}=$ 5°
b)$\varphi \text{'}=$ 10°
c)$\varphi \text{'}=$ 15°
d)$\varphi \text{'}=$ 20°
It is clearly that this numerical method gives a very good solution because most of obtained results are between those derived from lower bound (LB) and upper bound (UB) solutions given by Yamamoto et al. [17]. More interestingly, the total number of elements in the meshes used in the upper bound limit analysis ranges from 1672 to 3075 triangular elements (NSFEM), while there is a significantly larger number of elements employed in Yamamoto’s models (7680 triangular elements and 11412 stress/velocity discontinuities). The convergence of the method is examined in the stability analysis of dual circular tunnels for the case of $H/D=$1, $\gamma D/{c}^{\text{'}}=$ 1, $S/D=$ 1.25, ${\varphi}^{\text{'}}=$5° with a various number of elements in simulated model, and the results are summarized in Table 1. In comparison with those of Yamamoto et al. [17], it is recognized that the numerical procedure using NSFEM and SOCP not only reduces an appreciable amount of variables in optimization problem, but also gives a very better solution.
Fig. 13. Comparisons of the stability numbers between present method and Yamamoto et al. [17]. For the case $H/D=$5, smooth interface
a)$\varphi \text{'}=$ 5°
b)$\varphi \text{'}=$ 10°
c)$\varphi \text{'}=$ 15°
d)$\varphi \text{'}=$ 20°
Fig. 14. The convergence of stability factors of twin circular tunnels. For the case: $H/D=$1, $S/D=$ 1.25 and ${\varphi}^{\text{'}}=$5°
Table 1. The computational efficiency of present method using NSFEM and SOCP. For the case: $H/D=$1, $\gamma D/{c}^{\text{'}}=$1, $S/D=$1.25, $\varphi \text{'}=$ 5°
${\sigma}_{s}/c\text{'}$

0.4483

0.4081

0.3983

0.3906

0.3882

0.3871

0.3864

0.3856

Ne

523

914

1460

1935

3075

3975

5340

6408

Nvar

1555

2615

4060

5130

8270

10595

14110

16850

Iteration

17

18

20

20

20

21

20

20

Mosek time (s)

0.41

0.69

1.19

1.58

2.89

4.34

6.57

8.11

Ne = no. of elements, Nvar = no. of variables

7. Conclusions
In this paper, a numerical procedure based on the nodebased smoothed finite element method (NSFEM) is proposed to evaluate the stability of a plane strain twin circular tunnels in cohesivefrictional 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 $S/D$ 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.
Acknowledgements
This research was partially supported by the Foundation for Science and Technology at HCMUT University of Technology (Grant No. TNCS2015KTXD09). We sincerely appreciate Dr. Canh V. Le for sharing the upper bound analysis code. This support is gratefully acknowledged.
References
 Atkinson J. H., Cairncross A. M. Collapse of a Shallow Tunnel in a MohrCoulomb Material. Role of Plasticity in Soil Mechanics, Cambridge, 1973, p. 202206. [Search CrossRef]
 Cairncross A. M. Deformation Around Model Tunnels in Stiff Clay. Ph.D. Thesis, University of Cambridge, 1973. [Search CrossRef]
 Seneviratne H. N. Deformations and PorePressures Around Model Tunnels in Soft Clay. Ph.D. Thesis, University of Cambridge, 1979. [Search CrossRef]
 Mair R. J. Centrifugal Modelling of Tunnel Construction in Soft Clay. Ph.D. Thesis, University of Cambridge, 1979. [Search CrossRef]
 Davis E. H., Gunn M. J., Mair R. J., Seneviratne H. N. The stability of shallow tunnels and underground openings in cohesive material. Geotechnique, Vol. 30, Issue 4, 1980, p. 397416. [Publisher]
 Mühlhaus H. B. Lower bound solutions for circular tunnels in two and three dimensions. Rock Mechanics and Rock Engineering, Vol. 18, 1985, p. 3752. [Publisher]
 Leca E., Dormieux L. Upper and lower bound solutions for the face stability of shallow circular tunnels in frictional material. Geotechnique, Vol. 40, Issue 4, 1990, p. 581606. [Publisher]
 Sloan S. W., Assadi A. Undrained stability of a square tunnel in a soil whose strength increases linearly with depth. Computers and Geotechnics, Vol. 12, Issue 4, 1991, p. 321346. [Publisher]
 Lyamin A. V., Sloan S. W. Stability of a plane strain circular tunnel in a cohesive frictional soil. Proceedings of the J. R. Booker Memorial Symposium, Sydney, Australia, 2000, p. 139153. [Search CrossRef]
 Lyamin A. V., Jack D. L., Sloan S. W. Collapse analysis of square tunnels in cohesivefrictional soils. Proceedings of the First AsianPaciﬁc Congress on Computational Mechanics, Sydney, Australia, 2001, p. 405414. [Publisher]
 Yamamoto K., Lyamin A. V., Wilson D. W., Sloan S. W., Abbo A. J. Stability of a single tunnel in cohesivefrictional soil subjected to surcharge loading. Canadian Geotechnical Journal, Vol. 48, Issue 12, 2011, p. 18411854. [Publisher]
 Yamamoto K., Lyamin A. V., Wilson D. W., Sloan S. W., Abbo A. J. Stability of a circular tunnel in cohesivefrictional soil subjected to surcharge loading. Computers and Geotechnics, Vol. 38, Issue 4, 2011, p. 504514. [Publisher]
 Soliman E., Duddeck H., Ahrens H. Two and threedimensional analysis of closely spaced doubletube tunnels. Tunnelling and Underground Space Technology, Vol. 8, 1993, p. 1318. [Publisher]
 Wu B. R., Lee C. J. Ground movements and collapse mechanisms induced by tunneling in clayey soil. International Journal of Physical Modelling in Geotechnics, Vol. 3, 2003, p. 1529. [Publisher]
 Chehade F. H., Shahrour I. Numerical analysis of the interaction between twin tunnels: inﬂuence of the relative position and construction procedure. Tunnelling and Underground Space Technology, Vol. 23, 2008, p. 210214. [Publisher]
 Osman A. S. Stability of unlined twin tunnels in undrained clay. Tunnelling and Underground Space Technology, Vol. 25, 2010, p. 290296. [Publisher]
 Yamamoto K., Lyamin A. V., Wilson D. W., Sloan S. W., Abbo A. J. Stability of dual circular tunnels in cohesivefrictional soil subjected to surcharge loading. Computers and Geotechnics, Vol. 50, 2013, p. 4154. [Publisher]
 Hughes T. J. R. Reduced and selective integration techniques in the ﬁnite element analysis of plates. Nuclear Engineering and Design, Vol. 46, 1978, p. 203222. [Publisher]
 Hughes T. J. R. The ﬁnite Element Method. Dover Publications, PrenticeHall, 2000. [Search CrossRef]
 Piltner R., Taylor R. L. Triangular ﬁnite elements with rotational degrees of freedom and enhanced strain modes. Computers and Structures, Vol. 75, 2000, p. 361368. [Publisher]
 Simo J. C., Rifai M. S. A class of mixed assumed strain methods and the method of incompressible modes. International Journal for Numerical Methods in Engineering, Vol. 29, 1990, p. 15951638. [Publisher]
 Cardoso R. P. R., Yoon J. W., Mahardika M., Choudhry S., Alves de Sousa R. J., Fontes Valente R. A. Enhanced assumed strain (EAS) and assumed natural strain (ANS) methods for onepoint quadrature solidshell elements. International Journal for Numerical Methods in Engineering, Vol. 75, 2008, p. 156187. [Publisher]
 Bonet J., Burton A. J. A simple average nodal pressure tetrahedral element for incompressible and nearly incompressible dynamic explicit applications. Communications in Numerical Methods in Engineering, Vol. 14, 1998, p. 437449. [Publisher]
 Chen J. S., Wu C. T., Yoon S. A stabilized conforming nodal integration for Galerkin meshfree method. International Journal for Numerical Methods in Engineering, Vol. 50, 2001, p. 435466. [Publisher]
 Yoo J. W., Moran B., Chen J. S. Stabilized conforming nodal integration in the naturalelement method. International Journal for Numerical Methods in Engineering, Vol. 60, 2004, p. 861890. [Publisher]
 Liu G. R., NguyenThoi T. Smoothed ﬁnite Element Methods. CRC Press, New York, 2010. [Publisher]
 Liu G. R., Dai K. Y., NguyenThoi T. A smoothed ﬁnite element for mechanic problems. Computer and Mechanics, Vol. 39, 2007, p. 859877. [Publisher]
 Liu G. R., NguyenThoi T., NguyenXuan H., Lam K. Y. A node based smoothed ﬁnite element method (NSFEM) for upper bound solution to solid mechanics problems. Computer and Structures, Vol. 87, 2009, p. 1426. [Publisher]
 NguyenThoi T., Liu G. R., Lam K. Y., Zhang G. Y. A facebased smoothed ﬁnite element method (FSFEM) for 3D linear and nonlinear solid mechanics problems using 4node tetrahedral elements. International Journal for Numerical Methods in Engineering, Vol. 78, 2009, p. 324353. [Publisher]
 Liu G. R., NguyenThoi T., Lam K. Y. An edgebased smoothed ﬁnite element method (ESFEM) for static, free and forced vibration analyses of solids. Journal of Sound and Vibration, Vol. 320, 2009, p. 11001130. [Publisher]
 Liu G. R., NguyenThoi T., Dai K. Y., Lam K. Y. Theoretical aspects of the smoothed ﬁnite element method (SFEM). International Journal for Numerical Methods in Engineering, Vol. 71, 2007, p. 902930. [Publisher]
 Liu G. R., NguyenXuan H., NguyenThoi T. A theoretical study of SFEM models: properties, accuracy and convergence rates. International Journal for Numerical Methods in Engineering, Vol. 84, 2010, p. 12221256. [Publisher]
 NguyenXuan H., Rabczuk T., NguyenThoi T., Tran T. N., NguyenThanh N. Computation of limit and shakedown loads using a nodebased smoothed finite element method. International Journal for Numerical Methods in Engineering, Vol. 90, 2012, p. 287310. [Publisher]
 Le C. V., NguyenXuan H., Askes H., Bordas S., Rabczuk T., NguyenVinh H. A cellbased smoothed finite element method for kinematic limit analysis. International Journal for Numerical Methods in Engineering, Vol. 83, 2010, p. 16511674. [Publisher]
 NguyenXuan H., Liu G. R. An edgebased finite element method (ESFEM) with adaptive scaledbubble functions for plane strain limit analysis. Computer Methods in Applied Mechanics and Engineering, Vol. 285, 2015, p. 877905. [Publisher]
 NguyenXuan H., Rabczuk T. Adaptive selective ESFEM limit analysis of cracked planestrain structures. Frontiers of Structural and Civil Engineering, Vol. 9, 2015, p. 478490. [Publisher]
 NguyenXuan H., Wu C. T., Liu G. R. An adaptive selective ESFEM for plastic collapse analysis. European Journal of Mechanics – A/Solid, Vol. 58, 2016, p. 278290. [Search CrossRef]
 The MOSEK Optimization Toolbox for MATLAB Manual. http://www.mosek.com. [Search CrossRef]
 Jagdish Prasad Sahoo, Jyant Kumar Stability of long unsupported twin circular tunnels in soils. Tunnelling and Underground Space Technology, Vol. 38, 2013, p. 326335. [Publisher]
 Wilson D. W., Abbo A. J., Sloan S. W., Lyamin A. V. Undrained stability of dual circular tunnels. International Journal of Geomechanics, Vol. 14, 2014, p. 6979. [Publisher]
 Makrodimopoulos A., Martin C. M. Upper bound limit analysis using simplex strain elements and secondorder cone programming. International Journal for Numerical and Analytical Methods in Geomechanics, Vol. 31, 2006, p. 835865. [Publisher]
 GiD 11.0.4, Reference Manual. International Center for Numerical Methods in Engineering (CIMNE), http://www.cimne.com. [Search CrossRef]
 Chen W. F. Limit Analysis and Soil Plasticity. Elsevier, Amsterdam, 1975. [Search CrossRef]
Cited By
Computers and Geotechnics
T. VoMinh, L. NguyenSon

2021

Tunnelling and Underground Space Technology
T. VoMinh, L. NguyenSon, G. NguyenVan, T. ThaiPhuong

2021

International Journal of Computational Methods
Liming Zhou, Shuhui Ren, Yan Cai, Feng Li

2020

International Journal for Numerical Methods in Engineering
Shao‐Wei Wu, Chen Jiang, Chao Jiang, Gui‐Rong Liu

2020

Lecture Notes in Mechanical Engineering
T. VoMinh, T. NguyenMinh, A. ChauNgoc

2018
