Stability Analysis of Circular Tunnels in Cohesive-Frictional Soil Using the Node-Based Smoothed Finite Element Method (NS-FEM)

Abstract

In this chapter, the stability of a circular tunnel and dual circular tunnels in cohesive-frictional soils subjected to surcharge loading is investigated by using the node-based smoothed finite element method (NS-FEM). In the NS-FEM, the smoothing strain is calculated over smoothing domains associated with the elements’ nodes. The soil is assumed as a uniform Mohr-Coulomb material, and it obeys an associated flow rule. By using the second-order cone programming (SOCP) for solving the optimization problems, the ultimate load and failure mechanisms of the circular tunnel are considered. This chapter discusses the influence of the soil weight γD/c, the tunnel diameter ratio to its depth H/D, the vertical and horizontal spacing ratio (L/D, S/D) of two tunnels and soil internal friction angle ϕ on the stability numbers σs/c are calculated. The stability numbers obtained from the present approach are compared with the available literature for tunnels.

1. Introduction

In recent years, underground systems have become essential for the rapid development of many big cities. Underground infrastructures as an underground railway and gas pipeline have become increasingly popular in many metropolises to meet public demand. During tunnels’ construction, the overburden depth needs to be investigated carefully and plays an important role in constructing process and may reduce construction costs. Therefore, engineers need a practical approach to determine more precise the collapse load and failure mechanism in the circular tunnels’ preliminary design stage.

The first studied on the stability of a circular tunnel was performed at Cambridge University in the 1970s. Atkinson and Pott [1], Atkinson and Cairncross [2] investigated a series of centrifuge model tests of tunnels in dry sand and Mohr-Coulomb material subjected to surcharge loading. Cairncross [3] and Seneviratne [4] conducted a series of centrifuge model tests to determine the deformation around a circular tunnel in stiff clay and soft clay. Mair [5], Chambon and Corte [6] also conducted some centrifugal model tests to estimate circular tunnels’ stability in soft clay and sandy soil. Recently, Kirsch [7] and Idinger et al. [8] performed a small-scale tunnel model in a geotechnical centrifuge to investigate shallow tunnel face stability dry sand.

Some decades ago, several researchers have studied the stability of a tunnel in cohesive material using the upper and lower bound theorems, for example, the works of Davis et al. [9], Mühlhaus [10], Leca and Dormieux [11]. Recently, Zhang et al. [12] proposed a new 3D failure mechanism using the upper bound limit analysis theory to determine the tunnel face’s limit support pressure. In engineering practice, based on the 3D finite element method, Tosun [13] investigated the performance of concrete lining to compare with those obtained from observation and measurements during the excavation of rock masses and installing the temporary support system.

In recent decades, the finite element method using the triangular element (FEM-T3) has been rapidly developed to solve important geotechnical problems. Sloan and Assadi [14] first applied a finite element procedure for linear analysis to evaluate a square tunnel’s undrained stability in a soil whose shear strength increases linearly with depth. Then, Lyamin and Sloan [15], Lyamin et al. [16] and Yamamoto et al. [17, 18] used finite element limit analysis (FELA) to calculate the failure mechanisms of circular and square tunnels in cohesive-frictional soils. Recently, Yamamoto et al. [19], Xiao et al. [20] proposed an efficient method to calculate the stability numbers and failure mechanisms of dual circular tunnels in cohesive material; however, the nonlinear optimization procedure required large computational efforts.

However, one of the drawbacks of FEM-T3 elements is the volumetric locking phenomenon, which is often occurred in the nearly incompressible materials. To overcome this, Chen et al. [21] proposed a stabilized conforming nodal integration using the strain smoothing technique. Then, Liu and his co-workers [22, 23, 24, 25] applied this technique to standard FEM and proposed a class of smoothed finite element method (S-FEM). Typical S-FEM models include the cell-based S-FEM model (CS-FEM) [26], node-based S-FEM model (NS-FEM) [27, 28], and edge-based S-FEM model (ES-FEM) [29, 30, 31, 32]. Several papers demonstrate that the NS-FEM performs well in heat transfer analysis [33, 34], fracture analysis [35], acoustic problems [36, 37], axisymmetric shell structures [38], static and dynamic analysis [39, 40, 41]. Recently, T. Vo-Minh and his co-workers [42, 43, 44, 45] applied an upper bound limit analysis using NS-FEM and second-order cone programming (SOCP) to determine the twin circular and dual square tunnels’ stability cohesive-frictional soils.

This chapter aims to summarize our research papers [42, 43, 44, 45] using the node-based smoothed finite element method (NS-FEM) to estimate the collapse load and failure mechanism of single and two circular tunnels in cohesive-frictional soil subjected to surcharge loading. The numerical results are available for cases with ϕ ≤ 30°, and geotechnical engineers can use them in the preliminary design stage.

2. Upper bound limit analysis for a plane strain with Mohr-coulomb yield criterion using NS-FEM

2.1 A brief on the node-based smoothed finite element method

In the NS-FEM, the problem domain Ω is discretized by N_e triangular elements with totally N_n nodes and N_n smoothing domains Ω^(k) associated with the node k such that $\mathrm{\Omega }={\sum }_{k=1}^{N_n}{\mathrm{\Omega }}^{\left(k\right)}$ and Ωⁱ ∩ Ω^j = ∅, i ≠ j. The smoothing domain of the node k in NS-FEM is constructed based on the elements connected to the nodes k, as illustrated in Figure 1. The requirement of the smoothing domain is non-overlap and not required to be convex. Therefore, the smoothing domain is created by connecting the mid-edge-points sequentially to the surrounding triangles’ centroids.

Figure 1.

Triangular elements and smoothing cells associated with nodes.

The smoothed strain associated with the node k in the matrix form can be calculated by

$${\tilde{\mathbf{\epsilon }}}_k=\sum _{I\in N^{\left(k\right)}}{\tilde{\mathbf{B}}}_I\left({\mathbf{x}}_k\right){\mathbf{d}}_I$$E1

where N^(k) is a group of nodes associated with smoothing domain Ω_k^S, d_I is the nodal displacement vector and ${\tilde{\mathbf{B}}}_I\left({\mathbf{x}}_k\right)$ is the smoothed strain displacement matrix on the smoothing domain Ω _k^S that can be determined as

$${\tilde{\mathbf{B}}}_I\left({\mathbf{x}}_k\right)=\left[\begin{array}{l}{\tilde{b}}_{\mathit{Ix}}\left({\mathbf{x}}_k\right)0\\ 0{\tilde{b}}_{\mathit{Iy}}\left({\mathbf{x}}_k\right)\\ {\tilde{b}}_{\mathit{Iy}}\left({\mathbf{x}}_k\right){\tilde{b}}_{\mathit{Ix}}\left({\mathbf{x}}_k\right)\end{array}\right]$$E2

$${\tilde{b}}_{\mathit{Im}}\left({\mathbf{x}}_k\right)=\frac{1}{A^{\left(k\right)}}\underset{{\mathrm{\Gamma }}^{\left(k\right)}}{\int }{\mathbf{n}}_m^{\left(k\right)}\left(\mathbf{x}\right){\mathbf{N}}_I\left(\mathbf{x}\right)\mathrm{d\Gamma },\left(m=x,y\right)$$E3

where $A^{\left(k\right)}=\underset{{\mathrm{\Omega }}^{\left(k\right)}}{\int }d\mathrm{\Omega }$ is the area of the cell Ω_k^S, ${\mathbf{n}}_m^{\left(k\right)}\left(\mathbf{x}\right)$ is a matrix with components of the normal outward vector on the boundary Γ_k, N_I(x) is the FEM shape function for node I.

By using Gauss integration over each sub-boundary Γ_k of Ω _k^S, Eq. (3) can be rewritten as

$${\tilde{b}}_{\mathit{Im}}\left({\mathbf{x}}_k\right)=\frac{1}{A^{\left(k\right)}}\sum _{j=1}^{n_{\mathit{eg}}}{\mathbf{N}}_I\left({\mathbf{x}}_j^{\mathit{GP}}\right)n_{\mathit{jm}}^{\left(k\right)}{l}_j^{\left(k\right)},\left(m=x,y\right)$$E4

where n_eg is the total number of the sub-boundary segment of Γ_k, x_j^GP is the Gauss point of the sub-boundary segment of Γ_k which has length l_j^(k) and outward unit normal n_jm^(k).

2.2 A brief on the upper bound theorem

Consider a two-dimensional structure made of rigid-perfectly plastic materials with the domain Ω bounded by a continuous boundary ${\mathrm{\Gamma }}_{\dot{\mathbf{u}}}\cup {\mathrm{\Gamma }}_{\mathbf{t}}=\mathrm{\Gamma },{\mathrm{\Gamma }}_{\dot{\mathbf{u}}}\cap {\mathrm{\Gamma }}_{\mathbf{t}}=\varnothing .$ The structure subjected to body forces f and external tractions g on ${\mathrm{\Gamma }}_{\mathbf{t}}$ and the boundary ${\mathrm{\Gamma }}_{\dot{\mathbf{u}}}$ prescribed by the displacement velocity vector $\dot{\mathbf{u}}.$ The strain rates can express as:

$$\dot{\mathbf{\epsilon }}={\left[{\dot{\epsilon }}_{\mathit{xx}}{\dot{\epsilon }}_{\mathit{yy}}{\dot{\gamma }}_{\mathit{xy}}\right]}^T=\nabla \dot{\mathbf{u}}$$E5

The linear form of the external work rate can be calculated by

$$W_{\mathit{ext}}\left(\dot{\mathbf{u}}\right)=\underset{\mathrm{\Omega }}{\int }\mathbf{f}.\dot{\mathbf{u}}d\mathrm{\Omega }+\underset{{\mathrm{\Gamma }}_t}{\int }\mathbf{g}.\dot{\mathbf{u}}d\mathrm{\Gamma }$$E6

The internal plastic dissipation of the two-dimensional domain Ω can be written as

$$W_{\mathrm{int}}\left(\dot{\mathbf{\epsilon }}\right)=\underset{\mathrm{\Omega }}{\int }D\left(\dot{\mathbf{\epsilon }}\right)d\mathrm{\Omega }=\underset{\mathrm{\Omega }}{\int }\mathbf{\sigma }.\dot{\mathbf{\epsilon }}.d\mathrm{\Omega }$$E7

A space of kinematically admissible velocity field is denoted by

$$U=\left\{\dot{\mathbf{u}}\in {\left(H^1\left(\mathrm{\Omega }\right)\right)}^2,\dot{\mathbf{u}}=\overline{\dot{\mathbf{u}}}\begin{array}{c}\mathrm{on}{\mathrm{\Gamma }}_{\dot{\mathbf{u}}}\end{array}\right\}$$E8

We define a convex set that contains admissible stress fields

$$\mathbb{S}=\left\{\mathbf{\sigma }\in \sum \left|\psi \left(\mathbf{\sigma }\right)\le 0\right.\right\}$$E9

where Σ is symmetric stress tensors, ψ (σ) is the yield function.

The upper bound theorem states that there exists a kinematically admissible displacement field $\dot{\mathbf{u}}\in U,$ such that

$$W_{\mathrm{int}}\left(\dot{\mathbf{\epsilon }}\right)<{\alpha }^{+}{W}_{\mathrm{ext}}\left(\dot{\mathbf{u}}\right)+W_{\mathit{ext}}^0\left(\dot{\mathbf{u}}\right)$$E10

where ${\alpha }^{+}$ is the limit load multiplier of the load g,f and $W_{\mathit{ext}}^0\left(\dot{\mathbf{u}}\right)$ is the work of additional load g_o,t_o not subject to the multiplier.

Defining $C=\left\{\dot{\mathbf{u}}\in U|W_{\mathit{ext}}\left(\dot{\mathbf{u}}\right)=1\right\},$ the upper bound limit analysis becomes the optimization problem to determine the collapse multiplier ${\alpha }^{+}$

$${\alpha }^{+}=min\underset{\mathrm{\Omega }}{\int }D\left(\dot{\mathbf{\epsilon }}\right)d\mathrm{\Omega }-W_{\mathit{ext}}^0\left(\dot{\mathbf{u}}\right)$$E11

$$\mathit{st}\left\{\begin{array}{l}\dot{\mathbf{u}}=0\mathrm{on}{\mathrm{\Gamma }}_u\\ W_{\mathit{ext}}\left(\dot{\mathbf{u}}\right)=1\end{array}\right.$$E12

For plane strain in geotechnical problems, the Mohr-Coulomb yield function can be expressed as

$$\psi \left(\mathbf{\sigma }\right)=\sqrt{{\left({\sigma }_{\mathit{xx}}-{\sigma }_{\mathit{yy}}\right)}^2+4{\tau }_{\mathit{xy}}^2}+\left({\sigma }_{\mathit{xx}}+{\sigma }_{\mathit{yy}}\right)sin\varphi -2ccos\varphi$$E13

For an associated flow rule, the plastic strain rates vector is given by

$$\dot{\mathbf{\epsilon }}=\dot{\mu }\frac{\partial \psi \left(\mathbf{\sigma }\right)}{\partial \mathbf{\sigma }}$$E14

where $\dot{\mu }$ is the plastic multiplier.

Makrodimopoulos and Martin [46] used the Mohr-Coulomb failure criterion and associated flow rule to determine the power of plastic dissipation in geotechnical problems as follows

$$D\left(\dot{\mathbf{\epsilon }}\right)={\mathit{cA}}_i{t}_{i}cos\varphi$$E15

where A_i is the area of the element of node i, c is the cohesion, ϕ is the internal friction angle of soil, t_i is a vector of additional variables defined by

$$t_i\ge \sqrt{{\left({{\dot{\epsilon }}^i}_{\mathit{xx}}-{{\dot{\epsilon }}^i}_{\mathit{yy}}\right)}^2+{\left({\dot{\gamma }}_{\mathit{xy}}^i\right)}^2}$$E16

By using the smoothed strains rates ${\dot{\tilde{\mathbf{\epsilon }}}}_i$ in Eq. (1), the upper bound limit analysis problem for the plane strain using NS-FEM can be discretized in the simple form as follows

$${\alpha }^{+}=min\left(\sum _{i=1}^{N_n}{\mathit{cA}}_i{t}_{i}cos\varphi -W_{\mathit{ext}}^0\left(\dot{\mathbf{u}}\right)\right)$$E17

$$\mathit{st}\left\{\begin{array}{l}\dot{\mathbf{u}}=0\mathrm{on}{\mathrm{\Gamma }}_u\\ W_{\mathit{ext}}\left(\dot{\mathbf{u}}\right)=1\\ {{\tilde{\epsilon }}^i}_{\mathit{xx}}+{{\tilde{\epsilon }}^i}_{\mathit{yy}}=t_{i}sin\varphi ,i=1,2,..\dots \dots ,N_n\\ t_i\ge \sqrt{{\left({{\dot{\epsilon }}^i}_{\mathit{xx}}-{{\dot{\epsilon }}^i}_{\mathit{yy}}\right)}^2+{\left({\dot{\gamma }}_{\mathit{xy}}^i\right)}^2},i=1,2,..\dots \dots ,N_n\end{array}\right.$$E18

where α⁺ is a stability number, and A_i is the area of the smoothing domain of node i, N_n is the total number of nodes in the domain. The fourth constraint in Eq. (18) is a form of quadratic cones. As a result, the conic interior-point optimizer of the academic MOSEK package [47] is used for solving this problem.

3. Numerical examples

In this chapter, GiD [48] software was used to generate triangular elements with reduced element size close to the tunnel’s periphery. The domain’s size is assumed sufficiently large to eliminate the boundary effects and the plastic zones to be contained entirely within the domain. The computations were performed on a Dell Optiplex 990 (Intel Core™ i5, 1.6GHz CPU, 8GB RAM) in a Window XP environment. The NS-FEM approach has been coded in the Matlab language.

Example 1: Stability of a circular tunnel in cohesive-frictional soil.

Figure 2 shows the analysis model for a plane strain circular tunnel. The tunnel has diameter D and depth H. The soil behavior is modeled as a Mohr-Coulomb material with cohesion c, friction angle ϕ and unit weight γ. The surcharge loading σ_s is applied over the ground surface with smooth and rough interface conditions. Figure 3 illustrates the half of typical meshes of the upper bound limit analysis problem. The horizontal displacement of nodes is freeing or fixing along the ground surface, respectively, to describe the smooth or rough interface conditions between the loading and the soil. The following equation can calculate the upper bound limit analysis using NS-FEM:

Figure 2.

Model of a circular tunnel subjected to surcharge loading.

Figure 3.

Typical meshes of the circular tunnel using NS-FEM (H/D = 1).

$${\alpha }^{+}=\frac{{\sigma }_s}{c}=f\left(\frac{H}{D},\frac{\mathit{\gamma D}}{c},\varphi \right)$$E19

Figure 4a shows the plastic dissipation distribution of circular tunnel for shallow tunnel H/D = 1, ϕ = 15^o. A failure mechanism originates from the middle part of the tunnel and extends up to the ground surface. The power dissipation for medium and deep tunnels are illustrated in Figure 4b and c. The collapse surface develops from the bottom of the tunnel and extends up to the ground surface. It is noticeable that the failure mechanisms obtained by this proposed procedure are identical to those derived from rigid blocks and the results of Yamamoto et al. [18]. However, the values of stability number obtained from assuming rigid-block mechanism are greater than those of this proposed numerical procedure.

Figure 4.

Plastic dissipation distributions of circular tunnels. (a) H/D = 1, ϕ = 15^oγD/c = 1. (b) H/D = 2, ϕ = 15^oγD/c = 1. (c) H/D = 5, ϕ = 15^oγD/c = 1.

The stability numbers using NS-FEM for different values of ϕ, γD/c and H/D are listed in Tables 1 and 2 for smooth and rough interface conditions. The positive result means that the tunnel collapsed when subjected to compressive stress on the ground surface as this value. The negative value implies that theoretically, only normal tensile stress can be applied to the ground surface to sustain the tunnels’ stability. The stability results approximate zero are indicated by $"-"$; it means that the tunnel failure due to gravity occurs. The stability number results derived from this proposed method agree with the average values of the upper bound and lower bound reported by Yamamoto et al. [18], and illustrated in Figures 5 and 6.

ϕ	H/D	γD/c
		0	0.5	1	1.5	2	2.5	3
0	1	2.44	1.85	1.26	0.64	0.02	−0.62	−1.29
	2	3.46	2.33	1.18	0.02	−1.15	−2.35	−3.57
	3	4.13	2.48	0.80	−0.89	−2.60	−4.33	−6.08
	4	4.64	2.47	0.28	−1.92	−4.15	−6.39	−8.65
	5	5.04	2.35	−0.33	−3.07	−5.79	−8.56	−11.31
5	1	2.94	2.30	1.67	1.03	0.38	−0.28	−0.95
	2	4.42	3.18	1.94	0.68	−0.57	−1.84	−3.12
	3	5.47	3.64	1.81	−0.03	−1.88	−3.75	−5.62
	4	6.29	3.88	1.46	−0.96	−3.41	−5.86	−8.34
	5	6.99	3.98	1.00	−2.02	−5.06	−8.12	−11.23
10	1	3.65	2.95	2.26	1.56	0.87	0.15	−0.55
	2	5.88	4.49	3.10	1.70	0.30	−1.12	−2.55
	3	7.64	5.57	3.48	1.37	−0.76	−2.92	−5.13
	4	9.10	6.34	3.55	0.72	−2.16	−5.13	—
	5	10.40	6.95	3.46	−0.12	−3.79	−0.20	—
15	1	4.69	3.92	3.15	2.36	1.59	0.80	0.02
	2	8.31	6.71	5.09	3.45	1.80	0.13	−1.58
	3	11.50	9.06	6.58	4.04	1.44	−1.24	—
	4	14.38	11.08	7.70	4.21	0.59	−3.29	—
	5	17.09	12.92	8.65	4.15	−0.58	−0.21	—
20	1	6.36	5.46	4.57	3.67	2.77	1.85	0.95
	2	12.77	10.81	8.83	6.81	4.76	2.66	0.48
	3	19.25	16.20	13.05	9.79	6.42	2.87	−0.95
	4	25.65	21.47	17.09	12.49	7.63	2.39	—
	5	32.21	26.81	21.20	15.14	8.70	1.48	—
25	1	9.26	8.16	7.09	5.98	4.89	3.76	2.66
	2	21.97	19.42	16.82	14.14	11.39	8.53	5.55
	3	37.34	33.19	28.81	24.22	19.40	14.29	8.81
	4	54.33	48.44	42.14	35.43	28.23	20.43	11.82
	5	73.50	65.64	57.41	48.19	38.44	27.47	15.29
30	1	14.91	13.47	12.09	10.62	9.21	7.72	6.25
	2	44.59	40.88	37.03	33.03	28.64	24.49	19.93
	3	89.47	82.93	75.92	68.42	60.42	51.86	42.65
	4	146.15	136.32	125.60	113.94	101.20	87.35	72.18
	5	218.76	204.65	190.30	172.99	155.12	134.21	112.11
35	1	27.87	25.74	23.73	21.57	19.48	17.27	15.20
	2	114.91	108.36	101.47	94.18	86.50	78.47	69.98
	3	290.04	276.58	261.98	246.33	229.38	211.06	191.18
	4	551.96	530.01	505.28	478.04	447.58	414.36	377.71
	5	946.26	910.14	875.14	829.11	783.95	727.41	670.19

Table 1.

Stability numbers σ_s/c of a circular tunnel (smooth interface).

ϕ	H/D	γD/c
		0	0.5	1	1.5	2	2.5	3
0	1	2.51	1.92	1.33	0.72	0.10	−0.55	−1.22
	2	3.53	2.40	1.25	0.08	−1.09	−2.31	−3.53
	3	4.20	2.54	0.87	−0.82	−2.54	−4.27	−6.03
	4	4.70	2.53	0.33	−1.88	−4.11	−6.36	−8.63
	5	5.10	2.41	−0.29	−3.01	−5.75	−8.51	−11.29
5	1	3.03	2.40	1.76	1.11	0.46	−0.20	−0.89
	2	4.51	3.27	2.02	0.76	−0.51	−1.79	−3.07
	3	5.57	3.73	1.89	0.03	−1.83	−3.70	−5.59
	4	6.40	3.98	1.54	−0.89	−3.35	−5.82	−8.31
	5	7.09	4.08	1.07	−1.96	−5.02	−8.10	−11.23
10	1	3.78	3.08	2.38	1.67	0.96	0.24	−0.48
	2	6.04	4.63	3.22	1.80	0.38	−1.06	−2.51
	3	7.80	5.70	3.59	1.46	−0.69	−2.87	−5.10
	4	9.28	6.50	3.68	0.82	−2.09	−5.10	—
	5	10.59	7.11	3.58	−0.03	−3.75	−0.29	—
15	1	4.89	4.10	3.31	2.51	1.71	0.91	0.11
	2	8.58	6.94	5.28	3.61	1.93	0.22	−1.52
	3	11.79	9.30	6.77	4.19	1.55	−1.16	—
	4	14.72	11.37	7.94	4.39	0.70	−3.25	—
	5	17.49	13.26	8.89	4.33	−0.49	−0.25	—
20	1	6.66	5.74	4.82	3.89	2.96	2.03	1.08
	2	13.24	11.23	9.18	7.10	5.00	2.84	0.61
	3	19.80	16.67	13.44	10.11	6.65	3.03	−0.83
	4	26.39	22.10	17.61	12.91	7.94	2.55	—
	5	33.15	27.66	21.89	15.63	9.00	1.64	—
25	1	9.76	8.64	7.52	6.38	5.23	4.08	2.91
	2	22.95	20.30	17.58	14.79	11.93	8.98	5.91
	3	38.59	34.31	29.81	25.08	20.13	14.88	9.26
	4	56.16	50.07	43.56	36.59	29.17	21.12	12.18
	5	76.22	68.28	59.59	50.09	39.84	28.58	15.91
30	1	15.87	14.41	12.92	11.43	9.92	8.36	6.81
	2	46.90	43.02	38.99	34.79	30.43	25.89	21.13
	3	93.45	86.72	79.52	71.83	63.67	54.96	45.53
	4	151.97	141.69	130.44	118.30	105.10	90.62	74.69
	5	230.82	216.96	201.32	183.92	164.66	143.29	119.64
35	1	30.01	27.85	25.66	23.42	21.16	18.85	16.49
	2	122.01	115.10	107.88	100.34	92.42	84.10	75.30
	3	309.29	295.79	281.07	265.32	248.11	229.58	209.79
	4	581.11	557.87	532.05	503.36	471.31	436.51	398.22
	5	1038.50	1000.41	964.58	920.20	871.60	817.40	757.44

Table 2.

Stability numbers σ_s/c of a circular tunnel (rough interface).

ϕ	L/D	γD/c
		0	0.5	1	1.5	2	2.5	3
0	1.5	2.43	1.77	0.78	−0.34	−1.60	−2.92	−4.31
	2	2.21	1.63	0.48	−1.01	−2.72	−4.59	—
	3	2.43	1.84	0.33	−1.71	−3.93	—	—
	4	2.43	1.84	0.06	−2.28	—	—	—
	5	2.43	1.84	−0.91	—	—	—	—
5	1.5	2.94	2.25	1.28	0.17	−1.06	−2.37	−3.75
	2	2.66	2.07	0.99	−0.42	−2.11	—	−0.21
	3	2.94	2.30	1.22	−0.76	−3.07	—	—
	4	2.94	2.30	1.09	−1.24	—	—	—
	5	2.94	2.30	0.45	−2.98	—	—	—
10	1.5	3.63	2.92	1.93	0.84	−0.41	−1.76	−3.18
	2	3.29	2.64	1.65	0.32	−1.33	—	−0.35
	3	3.63	2.94	2.25	0.54	−1.82	—	—
	4	3.63	2.94	2.25	0.28	−2.57	—	—
	5	3.63	2.94	2.25	−0.91	—	—	—
15	1.5	4.65	3.89	2.85	1.73	0.47	−0.95	−2.47
	2	4.22	3.49	2.58	1.29	−0.26	—	−0.20
	3	4.65	3.89	3.13	2.31	0.20	—	—
	4	4.65	3.89	3.13	2.31	−0.12	—	—
	5	4.65	3.89	3.13	2.31	−1.78	—	—
20	1.5	6.32	5.37	4.26	3.07	1.78	0.30	−1.34
	2	5.70	4.85	3.97	2.71	1.23	−0.59	—
	3	6.32	5.43	4.53	3.64	2.74	0.60	—
	4	6.32	5.43	4.53	3.64	2.74	0.60	—
	5	6.32	5.43	4.53	3.64	2.74	0.60	—
25	1.5	9.12	7.92	6.66	5.32	3.91	2.37	0.58
	2	8.27	7.23	6.19	5.02	3.51	1.82	−0.27
	3	9.18	8.10	7.02	5.94	4.84	3.73	2.61
	4	9.18	8.10	7.02	5.94	4.84	3.73	2.61
	5	9.18	8.10	7.02	5.94	4.84	3.73	2.61
30	1.5	14.29	12.78	11.22	9.60	7.93	6.15	4.22
	2	13.29	11.92	10.53	9.13	7.66	5.85	3.85
	3	14.71	13.32	11.94	10.54	9.08	7.61	6.16
	4	14.71	13.32	11.94	10.54	9.08	7.61	6.16
	5	14.71	13.32	11.94	10.54	9.08	7.61	6.16
35	1.5	25.47	24.14	21.32	19.14	16.89	14.56	12.11
	2	24.77	22.72	20.64	18.55	16.43	14.27	12.05
	3	27.52	25.49	23.43	21.35	19.24	17.08	14.85
	4	27.52	25.49	23.43	21.35	19.24	17.08	14.85
	5	27.52	25.49	23.43	21.35	19.24	17.08	14.85

Table 3.

Stability numbers σ_s/c of two vertical circular tunnels (H/D = 1).

The stability numbers at the no-interaction points for dual vertical circular tunnels are highlighted in bold.

Figure 5.

The variation of stability numbers σ_s/c for different values of H/D (smooth interface).

Figure 6.

The variation of stability numbers σ_s/c for different values of H/D (rough interface).

Example 2: Stability of dual circular tunnels in cohesive-frictional soil.

Figure 7 illustrates the main geometrical parameters of two circular tunnels under plane strain conditions. The dual circular tunnels with the same diameter D, the cover depth of tunnel H, the vertical L and horizontal S distances between two tunnel centres. Continuous loading σ_s is applied to the ground surface. The soil is assumed to be rigid perfectly plastic and modeled by a Mohr-Coulomb yield criterion with cohesion c, friction angle ϕ, and unit weight γ. The typical mesh of dual circular tunnels is shown in Figure 8.

Figure 7.

Model of two circular tunnels subjected to continuous loading.

Figure 8.

The typical mesh of two circular tunnels using NS-FEM (H/D = 1, S/D = 3.5, L/D = 0.5).

The stability number σ_s/c is defined as a function of ϕ, γD/c, S/D, L/D and H/D to investigate two circular tunnels’ stability. The following equation can calculate the upper bound limit analysis using NS-FEM:

$${\alpha }^{+}=\frac{{\sigma }_s}{c}=f\left(\frac{H}{D},\frac{L}{D},\frac{S}{D},\frac{\mathit{\gamma D}}{c},\varphi \right)$$E20

3.1 Stability of two horizontal circular tunnels (L/D = 0, S/D ≠ 0)

Figure 9a–c show the distribution of power dissipation of shallow tunnel H/D = 1, ϕ = 15^o and γD/c = 1 at different values of S/D, i.e., S/D = 1.5, 2.0 and 3.5. In Figure 9a and b, 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. When the distance between two tunnels increases continuously and exceeds a certain value (S_c), i.e., S ≥ S_c = 3.5D as shown in Figure 9c, the failure mechanism becomes that of two single individual tunnels.

Figure 9.

Power dissipation of dual circular tunnels in the case H/D = 1. (a) γD/c = 1, S/D = 1.5, ϕ = 15^o. (b) γD/c = 1, S/D = 2, ϕ = 15^o. (c) γD/c = 1, S/D = 3.5, ϕ = 15^o.

Figure 10 shows the power dissipation of moderate tunnels H/D = 3, ϕ = 15^o and γD/c = 1 at different values of S/D, i.e., S/D = 2.0, 3.5 and 7.0. In Figure 10a and b, a slip failure between two circular tunnels enlarges to the top and bottom of tunnels, and a large surface originates the bottom of the tunnel extends up to the ground surface. When the distance between two tunnels increases continuously and exceeds a critical spacing (S_c), i.e., S ≥ S_c = 7D as shown in Figure 10c, the failure mechanism becomes that of two single individual tunnels and no influence on the failure mechanism of each tunnel.

Figure 10.

Power dissipation of dual circular tunnels in the case H/D = 3. (a) γD/c = 1, S/D = 2, ϕ = 15^o. (b) γD/c = 1, S/D = 3.5, ϕ = 15^o. (c) H/D = 3, γD/c = 1, S/D = 7, ϕ = 15^o.

The stability number results using NS-FEM for different values of ϕ, γD/c, S/D and H/D are listed in Table A1 of our research paper [43] and shown in Figures 11 and 12. The results derived from this proposed method agree well with the average values of the upper bound and lower bound reported by Yamamoto et al. [19]. The stability numbers at the no-interaction points for dual circular tunnels are highlighted in bold. When the spacing between the tunnels exceeds these points, the results obtained from NS-FEM tend to become constant. Therefore, the horizontal distance between two circular tunnels S/D plays an important role in the behavior of the failure mechanism. The increase in stability number is due to the effects of interaction between two circular tunnels.

Figure 11.

The variation of stability numbers σ_s/c for H/D = 1 (smooth interface).

Figure 12.

The variation of stability numbers σ_s/c for H/D = 3 (smooth interface).

3.2 Stability of two circular tunnels at different depth (L/D ≠ 0, S/D ≠ 0)

Figure 13a shows the typical power dissipation of two circular tunnels for shallow tunnel H/D = 1, ϕ = 15^o, L/D = 0.5, S/D = 1.5. It is noticed that a little slip failure occurs between two tunnels, and a large failure from the middle part of the tunnels extends up to the ground surface. When the distance between two tunnels increases continuously and exceeds a critical spacing (S_c), i.e., S ≥ S_c = 4D as shown in Figure 13b, only the tunnel near the ground surface failure and no interaction between dual circular tunnels.

Figure 13.

Power dissipation of dual circular tunnels in the case H/D = 1. (a) H/D = 1, L/D = 0.5, S/D = 1.5, γD/c = 1.5, ϕ = 15°. (b) H/D = 1, L/D = 0.5, S/D = 4, γD/c = 1.5, ϕ = 15°.

Figure 14a shows failure mechanism for moderate depth tunnel H/D = 3,ϕ = 15^o, L/D = 1, S/D = 3. 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 bottom of the tunnel extends up to the ground surface. When the distance between two tunnels is large enough and exceeds a critical spacing (S_c), i.e., S ≥ S_c = 8D as illustrated in Figure 14b, only the deeper tunnel failure and no influence on the failure mechanism of each tunnel.

Figure 14.

Power dissipation of dual circular tunnels in the case H/D = 3. (a) H/D = 3, L/D = 1, S/D = 3, γD/c = 1.5, ϕ = 15°. (b) H/D = 3, L/D = 1, S/D = 8, γD/c = 1.5, ϕ = 15°.

The stability number results using NS-FEM for different values of ϕ, γD/c, S/D, L/D and H/D are listed in Tables 2–4 of our research paper [45]. The results derived from this proposed method agree well with the average values of the lower bound and upper bound solution reported by Xiao et al. [20]. The stability numbers at the no-interaction points for dual circular tunnels are highlighted in bold. When the spacing between the tunnels exceeds these points, the results obtained from NS-FEM tend to become constant. Therefore, the horizontal distance between two circular tunnels S/D plays an important role in the failure mechanism’s behavior. The comparison of stability numbers between NS-FEM and Xiao et al. [20] solution is shown in Figures 15 and 16.

Figure 15.

The variation of stability numbers σ_s/c for H/D = 1, L/D = 1 (smooth interface).

Figure 16.

The variation of stability numbers σ_s/c for H/D = 3, L/D = 1 (smooth interface).

3.3 Stability of two vertical circular tunnels (L/D ≠ 0, S/D = 0)

Figure 17 shows the variation of the failure mechanisms in the case H/D = 1, L/D = 1.5, γD/c = 1 with different internal friction angle ϕ. In Figure 17a, a large slip surface develops from the middle of the above tunnel, and a slip failure originates from the bottom of the below tunnel extends up to the ground surface. In Figure 17b–d, the failure mechanism’s width decreases with increasing ϕ and the slip surface only originates from the above tunnel extends up to the ground surface.

Figure 17.

Power dissipation of dual vertical circular tunnels in the case H/D = 1, L/D = 1.5, γD/c = 1. (a) ϕ = 0°. (b) ϕ = 10°, (c) ϕ = 20°. (d) ϕ = 30°.

The variation of the failure mechanisms in the case H/D = 1, L/D = 3, γD/c = 1 with the different internal friction angle is illustrated in Figure 18. When a small friction angle ϕ ≤ 5^o, a large slip surface develops between two circular tunnels, and a small slip failure originates from the bottom of the below tunnel extends up to the ground surface, shown in Figure 18a and b. With increasing friction angle ϕ ≥ 10^o shown in Figure 18c, the failure mechanism only originates from the above tunnel extends to the ground surface and no failure mechanism on the below tunnel. Therefore, when the vertical distance between two tunnels exceeds a critical spacing (L_c), i.e., L ≥ L_c = 3D, the failure mechanism occurs with the above shallow tunnel and no effect on the deep tunnel. It means that the below tunnel is more stable when the soil internal friction angle ϕ increase and the slip surface only occurs from the top tunnel.

Figure 18.

Power dissipation of dual vertical circular tunnels in the case H/D = 1, L/D = 3, γD/c = 1. (a) ϕ = 0°. (b) ϕ = 5°. (c) ϕ = 10°, ϕ = 20°, ϕ = 30°.

The stability number results of two vertical using NS-FEM for different values of ϕ, γD/c, L/D in the case H/D = 1 are summarized in Table 3. The stability numbers at the no-interaction points for dual vertical circular tunnels are highlighted in bold. When the spacing between the tunnels exceeds these points, the results obtained from NS-FEM tend to become constant. Therefore, the vertical distance between two circular tunnels L/D plays an important role in the failure mechanism’s behavior.

4. Conclusions

Based on the upper bound limit analysis using NS-FEM, some concluding remarks can be shown as follow:

1. The stability numbers of a circular tunnel decrease continuously with increasing γD/c and increase with rising parameters H/D and ϕ.

2. A typical failure mechanism of two circular tunnels in cohesive-frictional soils consist of two parts: a small slip surface between the tunnels enlarges to the top and bottom of tunnels, and a large surface from the outside edge of tunnels extends up to the ground surface.

3. The stability results increase with increasing horizontal distance S/D for shallow dual tunnels (H/D = 1). In this case, the stability results increase with increasing horizontal distance S/D until it reaches a critical value S = 3.5D – 4D. The stability number tends to become constant, and this value is exactly equal to that of a single isolated tunnel. For the cases medium and deep tunnels H/D = 3, H/D = 5, the stability number increases until it reaches the approximate values of S = 8D and S = 12D, the stability number becomes constant and arrives at the maximum value.

4. In the case of two circular tunnels at a different depth, the horizontal distance S/D ratio plays an essential role in the behavior of dual circular tunnels’ failure mechanisms in cohesive-frictional soils. When the S/D ratio between two tunnels exceeds a certain value, the stability number tends to become constant, while the vertical distance L/D had no significant effect on the stability number. The failure mechanism becomes only a single tunnel and depends on the soil parameters ϕ and γ D/c.

Acknowledgments

This research was partially supported by the Foundation for Science and Technology at Ho Chi Minh City University of Technology (HUTECH). This support is gratefully acknowledged.

Conflict of interest

The authors declare no conflict of interest.