Theoretical Solution for Tunneling‐Induced Stress Field of Subdeep Buried Tunnel

Abstract

In the traditional Kirsch solution of stress field induced by tunneling in rock mass, the body force was not taken into consideration, and therefore the Kirsch solution is limited to demonstrate stress redistribution of deep‐buried tunnel. In order to consider the effect of body force on the stress redistribution induced by tunneling, a new secondary stress field solution for tunnel between shallow and deep tunnel (called subdeep tunnel) is developed with elastic mechanics and complex function employed. The stress field from theoretical solution is verified by numerical models, and the results showed good agreements with each other. This solution can be the basic theory in the analysis of the stress and field of subdeep tunnel, which have not been valuated through theoretical study yet.

1. Introduction

Tunnels are always classified as shallow and deep in tradition, and the classifying standard is the buried depth of tunnel, also known as h_q, which is the limit height that a pressure arch would form in the surrounding rock, and h_q is dependent on the quality of rock mass and tunnel height (Xu et al., 2000). This classification standard is just an empirical approach and is not based on mechanical behavior of rock mass. When tunneling in intact rock mass near the ground surface, the rock mass have partial self‐loading capacity which on the other hand is not enough to reach stability for rock mass. And in the presented research, a new type of tunnel classified as subdeep tunnel is described as: primary stress field of tunnel varies along the depth while the redistribution of secondary stress field of tunnel does not reach the ground surface or just have limited influence.

In general, the tunnel support design follows two kinds of strategies based on the tunnel classification of depth, and essentially, the support design strategy is a response to the stress redistribution and rock failure modes. For shallow tunnel, the surround rock has a very limited capacity to form a pressure arch and is very easy to collapse (Yang and Yang, 2008), and therefore the support system are designed to bear the whole weight of loosen zone above the tunnel (Peng and Liu, 2009; Terzaghi, 1943), and the loosen zone is determined by rock properties and tunnel section shapes (Wang et al., 2014). However, in this design approach self‐loading capacity of surround rock is not considered in subdeep tunnel, especially in intact rock mass, which would cause a significant waste of materials and always brings uneconomic results. On the other hand, according to the definition of subdeep tunnel, pressure arch which requires a significant burying depth cannot form and thus self‐stability cannot be expected, which makes researches on the pressure arch (Li, 2006; Poulsen, 2010; Sansone and Silva, 1998) and support design not suitable for subdeep tunnel.

For a better understanding of mechanical behavior of surround rock for subdeep tunnel, which is the basis for a reasonable support design, a theoretical solution of secondary stress field is put forward with elastic mechanics and complex function adopted. And the derivation process is described in detail as well as the results. The stress field including radial, tangential, and shear stresses are described and discussed. For an intuitive comparison, a numerical model is also built and the results from both theoretical and numerical models are presented, which is a powerful verification for the theoretical solution. Besides, the reasonable depth of subdeep tunnel is suggested by the theoretical results. And this solution for subdeep tunnel can be a very useful theoretical basis for safe and economic tunneling.

2. Theoretical model and primary stresses for subdeep tunnel

2.1. Subdeep tunnel model

In the stress analysis of subdeep tunnel, the following assumptions are made:

• The rock mass is elastic, homogeneous, and intact;

• Only self‐weight induced stress field exists in the tunnel site;

• The whole section of the tunnel is excavated in one step; and

• The tunnel is long enough and the model can be treated as planer strain.

As shown in Figure 1, a circular tunnel with a radius of a is excavated at depth of h_t. The vertical stress ${\sigma }_z$ at the bottom boundary is caused by the self‐weight of rock mass and ${\sigma }_z=0$ at the ground surface; horizontal stress ${\sigma }_x$ varies along z. And this stress field model can be decomposed into two parts, one is the primary stress field as shown in Figure 1(b), and the other is induced by stress‐released at tunnel boundary, as shown in Figure 1(c).

Figure 1.

Mechanical model and decomposition of the secondary stress field of subdeep tunnel. (a) Secondary stress field. (b) Primary stress field. (c) Stress released.

The primary stress field in Figure 1(b) can be described as follows according to Heim's hypothesis (Fegert, 2013) and Gold Nick's hypothesis (Dessler, 1982):

$\{\begin{array}{l}{\sigma }_z=\gamma (h_t-z)\hfill \\ {\sigma }_x=k_0{\sigma }_z=k_0\gamma (h_t-z)\hfill \\ {\tau }_{xz}=0\hfill \end{array}$E1

where z is the coordinate from tunnel center; $k_0$ is the lateral stress coefficient; $h_t$ is the depth from the tunnel center to ground; and $k_0=\frac{\nu }{1-\nu }$.

2.2. Analysis of the stress to release at tunnel boundary

The stress in polar coordinates in Figure 1 can be converted from Eq. (1) in rectangular coordinates as follows:

$\{\begin{array}{l}{\sigma }_r={\sigma }_x{\cos }^2\theta +{\sigma }_z{\sin }^2\theta +{\tau }_{xz}\sin 2\theta =\frac{{\sigma }_z+{\sigma }_x}{2}-\frac{{\sigma }_z-{\sigma }_x}{2}\cos 2\theta +{\tau }_{xz}\sin 2\theta \hfill \\ {\sigma }_{\theta }={\sigma }_x{\sin }^2\theta +{\sigma }_z{\cos }^2\theta -{\tau }_{xz}\sin 2\theta =\frac{{\sigma }_z+{\sigma }_x}{2}+\frac{{\sigma }_z-{\sigma }_x}{2}\cos 2\theta -{\tau }_{xz}\sin 2\theta \hfill \\ {\tau }_{r\theta }=\frac{{\sigma }_z-{\sigma }_x}{2}\sin 2\theta +{\tau }_{xz}\cos 2\theta \hfill \end{array}$E2

Where ${\sigma }_r$ = radial stress, ${\sigma }_{\theta }$= tangential stress, and ${\tau }_{r\theta }$= shear stress.

By submitting Eq. (1) into Eq. (2), the primary stress field can be expressed as:

$\{\begin{array}{l}{\sigma }_r=\frac{1+k_0}{2}\gamma (h_t-z)-\frac{1-k_0}{2}\gamma (h_t-z)\cos 2\theta \hfill \\ {\sigma }_{\theta }=\frac{1+k_0}{2}\gamma (h_t-z)+\frac{1-k_0}{2}\gamma (h_t-z)\cos 2\theta \hfill \\ {\tau }_{r\theta }=\frac{1-k_0}{2}\gamma (h_t-z)\sin 2\theta \hfill \end{array}$E3

where $z=r\sin \theta$.

When $r=a$, which means at the tunnel boundary:

$\{\begin{array}{l}{\sigma }_{ra}^0=\frac{1+k_0}{2}\gamma (h_t-a\sin \theta )-\frac{1-k_0}{2}\gamma (h_t-a\sin \theta )\cos 2\theta \hfill \\ {\sigma }_{\theta }^0{}_a=\frac{1+k_0}{2}\gamma (h_t-a\sin \theta )+\frac{1-k_0}{2}\gamma (h_t-a\sin \theta )\cos 2\theta \hfill \\ {\tau }_{r\theta a}^0=\frac{1-k_0}{2}\gamma (h_t-a\sin \theta )\sin 2\theta \hfill \end{array}$E4

Then, Eq. (4) is the stress to release in the following steps shown in Figure 1.

2.3. Complex function for elastic mechanics

With complex function, stress components in Figure 1 in polar coordinate can be expressed as

$z:\{\begin{array}{l}{\sigma }_{\theta }+{\sigma }_r=4\mathrm{Re}{{\chi }^{\prime }}_1(z)\\ {\sigma }_{\theta }-{\sigma }_r+2\text{i}{\tau }_{r\theta }=2\left[\overline{z}{{\chi }^{″}}_1(z)+{{\psi }^{\prime }}_1(z)\right]{\text{e}}^{2\text{i}\lambda }\end{array}$E5

where Re is the real part of complex function, $\overline{z}$ is the conjugate expression of z, ${\chi }_1\left(z\right)$ and ${\psi }_1\left(z\right)$ are complex potential function, ${\chi }_1^{\text{'}}\left(z\right)$ and ${\chi }_1^{\text{'}\text{'}}\left(z\right)$ are the first and second derivatives of ${\chi }_1\left(z\right)$, respectively, and ${\psi }_1^{\text{'}}\left(z\right)$ is the first derivative of ${\psi }_1\left(z\right)$.

2.4. Stress release

Similar to the solving process of an infinite plate with a hole (Wang, 2008), as the stress field approximates to zero, the solving process are described as follows:

1. Submitting Eq. (4) into the Fourier coefficient equations;

2. Through the obtained equations, solving the constant coefficients; and

3. As the problem has been assumed to be a planar strain work, which means $k_0=\frac{\nu }{1-\nu }={\nu }^{\text{'}}$, the following equations can be obtained based on steps (1) and (2):

   $\{\begin{array}{l}{{\chi }^{\prime }}_1(z)=\frac{a_1}{z}+\frac{a_2}{z^2}+\frac{a_3}{z^3}=\frac{\left(1+{\nu }^{\prime }\right)\gamma a^2}{8}\frac{\text{i}}{z}+\frac{\left(1-k_0\right)\gamma h_t{a}^2}{2}\frac{1}{z^2}-\frac{\left(1-k_0\right)\gamma a^4}{4}\frac{\text{i}}{z^3}\hfill \\ \overline{{{\chi }^{\prime }}_1(z)}=\frac{{\overline{a}}_1}{\overline{z}}+\frac{{\overline{a}}_2}{{\overline{z}}^2}+\frac{{\overline{a}}_3}{{\overline{z}}^3}=-\frac{\left(1+{\nu }^{\prime }\right)\gamma a^2}{8}\frac{\text{i}}{\overline{z}}+\frac{\left(1-k_0\right)\gamma h_t{a}^2}{2}\frac{1}{{\overline{z}}^2}+\frac{\left(1-k_0\right)\gamma a^4}{4}\frac{\text{i}}{{\overline{z}}^3}\hfill \\ {{\chi }^{″}}_1(z)=\frac{-a_1}{z^2}+\frac{-2a_2}{z^3}+\frac{-3a_3}{z^4}=-\frac{\left(1+{\nu }^{\prime }\right)\gamma a^2}{8}\frac{\text{i}}{z^2}-\frac{2\left(1-k_0\right)\gamma h_t{a}^2}{2}\frac{1}{z^3}+\frac{3\left(1-k_0\right)\gamma a^4}{4}\frac{\text{i}}{z^4}\hfill \\ {{\psi }^{\prime }}_1(z)=\frac{b_1}{z}+\frac{b_2}{z^2}+\frac{b_3}{z^3}+\frac{b_4}{z^4}+\frac{b_5}{z^5}=\frac{\left(3-{\nu }^{\prime }\right)\gamma a^2}{8}\frac{\text{i}}{z}+\frac{\left(1+k_0\right)\gamma h_t{a}^2}{2}\frac{1}{z^2}\text{ +}\frac{\left({\nu }^{\prime }-k_0\right)\gamma a^4}{4}\frac{\text{i}}{z^3}\hfill \\ +\frac{3\left(1-k_0\right)\gamma h_t{a}^4}{2}\frac{1}{z^4}-\left(1-k_0\right)\gamma a^6\frac{\text{i}}{z^5}\hfill \\ {\chi }_1(z)=\frac{\left(1+{\nu }^{\prime }\right)\gamma a^2}{8}\text{i}\ln z-\frac{\left(1-k_0\right)\gamma h_t{a}^2}{2}\frac{1}{z}+\frac{\left(1-k_0\right)\gamma a^4}{8}\frac{\text{i}}{z^2}+c_1\hfill \\ {\psi }_1(z)=\frac{\left(3-{\nu }^{\prime }\right)\gamma a^2}{8}\text{i}\ln z-\frac{\left(1+k_0\right)\gamma h_t{a}^2}{2}\frac{1}{z}-\frac{\left({\nu }^{\prime }-k_0\right)\gamma a^4}{8}\frac{\text{i}}{z^2}-\frac{\left(1-k_0\right)\gamma h_t{a}^4}{2}\frac{1}{z^3}+\frac{\left(1-k_0\right)\gamma a^6}{4}\frac{\text{i}}{z^4}+c_2\hfill \\ \overline{{\psi }_1(z)}=\frac{\left(3-{\nu }^{\prime }\right)\gamma a^2}{8}\text{i}\ln \overline{z}-\frac{\left(1+k_0\right)\gamma h_t{a}^2}{2}\frac{1}{\overline{z}}+\frac{\left({\nu }^{\prime }-k_0\right)\gamma a^4}{8}\frac{\text{i}}{{\overline{z}}^2}-\frac{\left(1-k_0\right)\gamma h_t{a}^4}{2}\frac{1}{{\overline{z}}^3}-\frac{\left(1-k_0\right)\gamma a^6}{4}\frac{\text{i}}{{\overline{z}}^4}+{\overline{c}}_2\hfill \end{array}$E6

By submitting Eq. (6) into Eq. (5), the relation between stress components can be expressed by:

$\{\begin{array}{c}{\sigma }_{\theta }+{\sigma }_r=2\mathrm{Re}\left(\frac{\left(1+{\nu }^{\prime }\right)\gamma a^2}{4}\frac{\text{i}}{z}+\left(1-k_0\right)\gamma h_t{a}^2\frac{1}{z^2}-\frac{\left(1-k_0\right)\gamma a^4}{2}\frac{\text{i}}{z^3}\right)\hfill \\ {\sigma }_{\theta }-{\sigma }_r+2\text{i}{\tau }_{r\theta }=2\left[\overline{z}{{\chi }^{″}}_1(z)+{{\psi }^{\prime }}_1(z)\right]{\text{e}}^{2\text{i}\theta }\hfill \\ =2[-\frac{\left(1+{\nu }^{\prime }\right)\gamma a^2}{8}\frac{\overline{z}\text{i}}{z^2}-\frac{2\left(1-k_0\right)\gamma h_t{a}^2}{2}\frac{\overline{z}}{z^3}+\frac{3\left(1-k_0\right)\gamma a^4}{4}\frac{\overline{z}\text{i}}{z^4}+\frac{\left(3-{\nu }^{\prime }\right)\gamma a^2}{8}\frac{\text{i}}{z}\\ +\frac{\left(1+k_0\right)\gamma h_t{a}^2}{2}\frac{1}{z^2}\text{+}\frac{\left({\nu }^{\prime }-k_0\right)\gamma a^4}{4}\frac{\text{i}}{z^3}+\frac{3\left(1-k_0\right)\gamma h_t{a}^4}{2}\frac{1}{z^4}-\left(1-k_0\right)\gamma a^6\frac{\text{i}}{z^5}]{\text{e}}^{2\text{i}\theta }\hfill \end{array}$E7

where z, which is a complex valuable, can be expressed as:

$\{\begin{array}{l}\frac{\text{i}}{z}=\frac{\text{i}}{r{\text{e}}^{\text{i}\theta }}=\frac{{\text{ie}}^{-\text{i}\theta }}{r}=\frac{\sin \theta +\text{i}\cos \theta }{r}\hfill \\ \frac{\text{1}}{z^2}=\frac{\text{1}}{r^2{\text{e}}^{\text{i2}\theta }}=\frac{{\text{e}}^{-\text{i2}\theta }}{r^2}=\frac{\cos 2\theta -\text{i}\sin 2\theta }{r^2}\hfill \\ \frac{\text{i}}{z^3}=\frac{\text{i}}{r^3{\text{e}}^{\text{i3}\theta }}=\frac{{\text{ie}}^{-\text{i3}\theta }}{r^3}=\frac{\sin 3\theta \text{+i}\cos 3\theta }{r^3}\hfill \end{array}$E8

By submitting Eq. (8) into the Eq. (7), a more simple relation between stress components are obtained, and by decomposing the imaginary part in the obtained equations, the subsidiary stress induced by the stress release process in Figure 1 is obtained. The obtained subsidiary stress is based on the infinite plate assumption, and it is an approximate value. For a more accurate solution, the subsidiary stress at ground surface should be released again, which is a very complicated process. By solving the obtained subsidiary stress components with Eq. (3), the stress field solution for subdeep tunnel can be expressed as:

$\{\begin{array}{l}{\sigma }_r=\frac{1+k_0}{2}\gamma (h_t-r\sin \theta )-\frac{1-k_0}{2}\gamma (h_t-r\sin \theta )\cos 2\theta +\gamma \{-\frac{\left(1+k_0\right)h_t{a}^2}{2r^2}\hfill \\ \begin{array}{cc}& +\left[\frac{\left(3+{\nu }^{\prime }\right)a^2}{4r}+\frac{\left({\nu }^{\prime }-k_0\right)a^4}{4r^3}\right]\sin \theta +\left(1-k_0\right)h_t\left(\frac{2a^2}{r^2}-\frac{3a^4}{2r^4}\right)\cos 2\theta -\left(1-k_0\right)a\left(\frac{5a^3}{4r^3}-\frac{a^5}{r^5}\right)\sin 3\theta \}\end{array}\hfill \\ {\sigma }_{\theta }=\frac{1+k_0}{2}\gamma (h_t-r\sin \theta )+\frac{1-k_0}{2}\gamma (h_t-r\sin \theta )\cos 2\theta +\gamma \{\frac{\left(1+k_0\right)h_t{a}^2}{2r^2}\hfill \\ \begin{array}{cc}& -\left[\frac{\left(1-{\nu }^{\prime }\right)a^2}{4r}-\frac{\left({\nu }^{\prime }-k_0\right)a^4}{4r^3}\right]\sin \theta +\left(1-k_0\right)h_t\frac{3a^4}{2r^4}\cos 2\theta +\left(1-k_0\right)a\left(\frac{a^3}{4r^3}-\frac{a^5}{r^5}\right)\sin 3\theta \}\end{array}\hfill \\ {\tau }_{r\theta }=\frac{1-k_0}{2}\gamma (h_t-r\sin \theta )\sin 2\theta \hfill \\ \begin{array}{cc}& +\gamma \left\{\left[\frac{\left(1-{\nu }^{\prime }\right)a^2}{4r}\text{+}\frac{\left({\nu }^{\prime }-k_0\right)a^4}{4r^3}\right]\cos \theta +\left(1-k_0\right)h_t\left(\frac{a^2}{r^2}-\frac{3a^4}{2r^4}\right)\sin 2\theta +\left(1-k_0\right)a\left(\frac{3a^3}{4r^3}-\frac{a^5}{r^5}\right)\cos 3\theta \right\}\end{array}\hfill \end{array}$E9

3. Distribution of secondary stress field solution for subdeep tunnel

The secondary stress field of subdeep tunnel calculated from Eq. (9), and the $v$, and $h_t$ are set at 0.3 and 6a, respectively. The results are shown in Figure 2.

Figure 2.

Second stress distribution in subdeep tunnels. (a) Distribution of σ_r. (b) Distribution of σ_θ. (c) Distribution of τ_rθ.

Stress distribution along r different directions in polar coordinate are also illustrated in Figure 3 in a normalized form.

Figure 3.

Stress distribution along r different directions. (a) Distribution of σ_r along r. (b) Distribution of σ_θ along r. (c) Distribution of τ_rθ along r.

From Figures 2 and 3, the distributions of secondary stress field induced by tunneling for subdeep tunnel are different from the deep‐buried tunnels, which has the following characteristics of its own:

1. The distribution of secondary subdeep tunnel just shows axial symmetry, not like the deep‐buried tunnel, whose stress field shows bidirectional symmetry. This difference is caused by the variation of ${\sigma }_z$ along z in the primary field.

2. The tangential stress near the tunnel boundary shows significant concentration at the side wall of the tunnel, to which close attention should be paid.

3. Figure 3 shows the variation of stress induced by tunneling along polar radius, and it is quite clear that, the stresses are approximate to 0, which proves the simplification of zero stress boundary to be reasonable and correct. Besides, subdeep tunnel can be determined as tunnels buried at depth of more than 2.5 D (D is the diameter of tunnel), and less than the depth of deep tunnel.

From the stress analysis, secondary stress field for subdeep tunnel is demonstrated and the limit depth to distinguish shallow tunnel and subdeep tunnel is obtained, which can be the guideline for tunnel support in a more reasonable, economic, and safer way.

4. Numerical solution and its comparison with theoretical results

To verify the theoretical solution, a numerical model is developed. The parameters for theoretical and numerical simulation are shown in Table 1. And a numerical model is developed in Flac3D as shown in Figure 4.

Figure 4.

Numerical model of subdeep tunnels.

Parameters of theoretical solution				Parameters of numerical simulation
Item	Symbol	Value	Unit	Item	Symbol	Value	Unit
Tunnel’s radius	a	3	m	Tunnel’s radius	a	3	m
Depth of tunnel center	ht	6 a	m	Depth of tunnel	d	15	m
Volume weight	γ	2.10E+04	N/m^3	Density	ρ	2.10E+03	kg/m^3
Elasticity modulus	E	9.00E+07	Pa	Bulk modulus	K	7.50E+07	Pa
Poisson’s ratio	ν	0.30	‐	Shear modulus	G	3.46E+07	Pa

Table 1.

Parameters for theoretical and numerical simulation.

Theoretical solution of horizontal stresses (σ_xx) of surrounding rock in subdeep tunnel is shown in Figure 5(a) in a rectangular coordinate system, and its numerical result is shown in Figure 5(b).

Figure 5.

Theoretical and numerical solution of horizontal stress. (a)Theoretical solution. (b) Numerical solution.

Theoretical solution of vertical stresses (σ_zz) of surrounding rock in subdeep tunnel is shown in Figure 6(a), and the numerical result is shown in Figure 6(b).

Figure 6.

Theoretical solutions and numerical value of vertical stress. (a) Theoretical solution. (b) Numerical solution.

Theoretical solution of shear stresses (τ_xz) of surrounding rock in subdeep tunnel is shown in Figure 7(a), and its numerical result is shown in Figure 7(b).

Figure 7.

Theoretical and numerical solution of shear stress. (a) Theoretical solution. (b) Numerical solution.

Contours of horizontal and vertical stress field from theoretical and numerical solution show good agreement with each other as shown in Figures 5 and 6, which proves the theoretical analysis reasonable and correct. On the other hand, shear stress from theoretical and numerical solution shows difference to some extent. In Figure 7, shear stress, from numerical solution, has a much smaller distribution area than that from numerical model, a probable explanation is that, the media in theoretical model is elastic while that in numerical model is elastic‐plastic, and stress concentration near the tunnel boundary may redistribute after rock failure, leading to a smoother distribution of shear stress in numerical model.

5. Conclusions

In the presented research, the stress field for tunnel buried at depth between deep tunnel and shallow tunnel is analyzed with elastic mechanics and complex function, and conclusions can be drawn as:

1. The subdeep tunnel is proposed and through theoretical analysis, and stress fields show that essential difference exists between the mechanical behavior of deep tunnel and subdeep tunnel, and the depth to distinguish shallow tunnel and subdeep tunnel is suggested as 2.5 times of the tunnel diameter.

2. By numerical analysis, the theoretical solution is proved to be reasonable and with high accuracy. And this theoretical solution can be a very good guideline for the support design for subdeep tunnel with economy and safety.

3. This theoretical solution also has its limitations, and if used for deep tunnels, the vertical stress in rock would be incorrect totally. Besides, this solution is not suitable for tunnels buried in depth less than 2.5 D, which does not satisfy the boundary stress condition.

4. Through the comparison of theoretical analysis and numerical model, the theoretical results are proved to be effective in the determination of subdeep tunnel, which can be very helpful in the design of subtunnel support for economy and safety purpose.