Topics of Analytical and Computational Methods in Tunnel Engineering

Abstract

In this chapter, a selection of tunneling topics is presented, following the evolution of methods and tools from analytical to computational era. After an introductory discussion of the importance of elasticity and plasticity in tunneling, some practical topics are presented as paradigms to show the successful application of them in achieving a solution. The circular and horseshoe tunnel sections served as the basis of the elastic analysis of deep tunnels. Practical aspects such as influence zone and elastic convergences in both cases are examined. In the case of circular tunnels, the estimation of plastic zone formation is discussed for a selection of strength criteria. After a detailed discussion of the influence of surface proximity, the elastic and plastic analysis of shallow tunnels is examined in some detail. The presentation is completed by a short presentation of computational methods. An overview of recent developments and a classification of the methods are presented, and then some problems for the case of anisotropic rocks have been presented using finite element method (FEM). The last topic is the application of artificial intelligence (AI) tools in interpreting data and in estimating the relative importance of parameters involved in the problem of tunneling-induced surface settlements. In the conclusions a short discussion of the main topics presented follows.

1. Introduction

Underground work history goes back to prehistoric times [1]. The oldest known tunnel is the one underpassing the River Euphrates in Babylon constructed 4000 years ago. Hezekiah, King of Judea, built a tunnel 2700 years ago, whereas Eupalinos built the Eupalinion tunnel constructed in Samos Island, Greece, 2600 years ago, both for water supply purposes. For further information see [2, 3, 4, 5]. Shelters, underground works for warfare purposes, mineral resource exploitation, traffic tunnels and conveyance tunnels like water supply tunnels, etc. are the main categories of underground structures [3].

Tunnel engineering could be characterized as an art, not fine art of course, with the meaning of know-how based on past experience, lessons learned from tunnel disasters [6], knowledge of geological conditions and behaviour of the ground, innovations to overcome encountered difficulties, scientific knowledge from applied mechanics and advances in mechanized excavation tools and methods. In recent years, the experience from South Africa and North Europe rock engineering practice resulted in the introduction of comprehensive classification systems, namely, RMR [7, 8] and Q system [9]. In theory, Lamé in 1852 and Kirsch in 1898 [10] paved the road of scientific development of solutions to the problem of stresses and strains around underground holes [10]. Terzaghi [11] examined shallow tunneling through sands and cohesive soils. For a comprehensive presentation of the fundamentals of theoretical rock mechanics, see Jaeger et al. [12]. For a detailed coverage of underground excavations in rock, see Hoek and Brown [13]. Experimental methods such as photoelasticity [14, 15] and Moiré [16] and computational methods such as finite elements [17], finite difference, direct and indirect boundary elements [13, 18], distinct elements [19] and meshfree methods [20, 21] have been developed and extensively applied successfully. On the other hand, the progress regarding the material behaviour resulted in advanced constitutive laws [22] and in the more realistic modeling of continua and discontinua media [23]. Finally, the information systems and in particular the artificial intelligence offered new tools like artificial neural networks (ANN) and other machine learning methods.

2. Analytical methods: elastic problems

2.1 The case of deep tunnels

Deep tunnels can be considered as the limiting case of a cylinder with thick walls when the radius of the exterior wall tends to infinity. In that case the ground surface has no any influence on the stress field around the opening.

Jeffery [24] solved the problem of stress distribution inside the walls of a cylinder bounded by two non-concentric circles, with the distance of their centres denoted by d. In this way, Jeffery’s approach solves the problems of an infinite plate containing a circular hole as the limit of a thick cylinder with concentric boundaries (d = 0) and the problem of a semi-infinite plate containing a circular hole as a particular case of two eccentric boundaries. So, he gave the general solution of deep and shallow tunnels, respectively, as particular cases of a general problem.

For a historical account of the mathematical solutions of stress distribution around underground openings, see Gerçek [25].

2.1.1 Cyclic section

Assuming a homogeneous isotropic linear elastic medium (HILE medium) and considering the geometry of the problem and the equations of equilibrium, we obtain the following set of equations (Eqs. (1)–(4)) ( Figure 1 ):

Figure 1.

The cylinder with thick walls.

$${\sigma }_r=-\frac{2G}{1-2\nu }\left\{\frac{A}{2}-\left(1-2\nu \right)\frac{B}{r^2}\right\},\Rightarrow {\sigma }_r=C-\frac{D}{r^2},$$ E1

$${\sigma }_{\theta }=-\frac{2G}{1-2\nu }\left\{\frac{A}{2}+\left(1-2\nu \right)\frac{B}{r^2}\right\},\Rightarrow {\sigma }_{\theta }=C+\frac{D}{r^2},$$ E2

$${\tau }_{\mathit{r\theta }}=0$$ E3

$$u_r=-\frac{1}{2G}\left\{\left(1-2\nu \right)\mathit{Cr}+\frac{D}{r}\right\}$$ E4

The constants C and D are defined as:

$$\begin{array}{l}C=-\frac{\mathit{GA}}{1-2\nu },\\ D=-2\mathit{GB}\end{array}$$

Considering the boundary conditions:

1. for $r=a$ , ${\sigma }_r=p_a$ .

2. $r=b$ , ${\sigma }_r=p_b$ .

therefore:

$$\begin{array}{l}{p}_a=C-\frac{D}{a^2},\\ p_b=C-\frac{D}{b^2}\\ D=\frac{\left(p_b-p_a\right)a^2{b}^2}{\left(b^2-a^2\right)},\\ C=\frac{\left(b^2{p}_b-a^2{p}_a\right)}{\left(b^2-a^2\right)}\end{array}$$

When b tends to infinity, we obtain the special case of a cyclic opening excavated in an infinite medium under hydrostatic stress field p, by superposing the initial stresses induced by the excavation:

$$\begin{array}{l}{\sigma }_r=p\left(1-\frac{a^2}{r^2}\right),\\ {\sigma }_{\theta }=p\left(1+\frac{a^2}{r^2}\right),\\ {\tau }_{\mathit{r\theta }}=0,\\ u_r=-\frac{1}{2G}\frac{a^2}{r}p\end{array}$$ E5

For the general case of stress field, denoting the ratio of horizontal to vertical in situ principal stress with k, Kirsch [10] obtained the following equations:

E6

It is interesting to note that in the case of hydrostatic or isotropic stress field (k = 1), we can estimate an influence zone accepting a ± 0.04% deviation from the initial stress field. Then by introducing this condition in the above equations, we obtain the radius of influence zone as 5α, where α is the radius of the excavation.

2.1.2 The case of horseshoe section

Let us examine now the case of horseshoe section ( Figures 2 and 3 ). This problem is very useful in the case of excavations without using a tunnel boring machine (TBM). Figure 2 shows a typical metro section, whereas Figure 3 shows the simplified section used in the analysis.

Figure 2.

Typical horseshoe tunnel section.

Figure 3.

A simplified horseshoe tunnel section.

First, the maximum principal stress σ₁ contours obtained by a numerical analysis are shown Figure 4 . A slight influence of the asymmetric shape of the tunnel’s section can be observed. Schürch and Anagnostou [26] examined the influence of rotational symmetry violation on the applicability of the ground response curve.

Figure 4.

Maximum principal stress contours around a deep horseshoe tunnel section.

To obtain a solution to this problem, we adopted Muskhelishvili approach, using functions of complex variables and conformal mapping techniques [27].

Figure 5.

Conformal mapping of an infinite region surrounding a hole onto the unit circle [28].

Adopting Gerçek solution [28], we introduce the following mapping function, which transforms the infinite region surrounding the opening onto the interior of the unit circle ( Figure 5 ):

$$z=\omega \left(\zeta \right)=R\left(\frac{1}{\zeta }+\sum _{k=1}^3{a}_k\cdot {\zeta }^k\right)$$ E7

In Eq. (7) R is a real constant, and the complex coefficients α_k are defined as:

E8

In order to find the optimum values of the above coefficients to fit best the given horseshoe section, a (123x3) system of equations has been solved [29, 30]. The values of the coefficients are:

$$\begin{array}{c}{\mathrm{\alpha }}_1={\mathrm{b}}_3=0\\ {\mathrm{\alpha }}_2={\mathrm{b}}_2=0.047,\\ \begin{array}{c}{\mathrm{\alpha }}_3=0.029,\\ {\mathrm{b}}_1=0.104,\\ \mathrm{R}=3.756\end{array}\end{array}$$ E9

In [29, 30] an extended investigation of this problem is presented resulting in the calculation of the set of coefficients for 49 sections. In Figure 6 , the final section, which is the best for our problem, is presented superimposed to the original horseshoe shape.

Figure 6.

The final section that best fits the initial horseshoe section [29, 30].

Now, we can calculate the stresses around the opening as well as the strains for initial stress fields with k = 0, 0.333, 1 and 3 [29, 30]. In Figures 7 and 8 , the variation of σ_θ at the boundary of the excavation is shown for k = o.333 ( Figure 7 ) and k = 1 ( Figure 8 ).

Figure 7.

Hoop stress σ_θ variation along the boundary of the excavation for k = 0.333 [29].

Figure 8.

Hoop stress σ_θ variation along the boundary of the excavation for k = 1 [30].

A further result of practical significance obtained by this method is the calculation of the convergence of the excavation boundaries and the extent of the zone of influence. In Tables 1 and 2 , these values are presented in comparison with the corresponding values for the circular section. We can observe the greater deviation at the bottom of the horseshoe section because of the greater radius of this segment compared to the circle.

Initial stress field	Convergences of circular section [mm]	Horseshoe-induced roof convergences [mm]	Horseshoe-induced bottom convergences [mm]
k = 1	7.8	7.6	8.2
k = 0.333	7.1	7.0	7.9
k = 1	5.6	5.9	7.2
k = 3	1.1	2.5	5.0

Table 1.

Convergences of circular and horseshoe sections for different stress fields [29, 30].

Initial stress field	Influence zone (horseshoe section)	Influence zone (circular section)
k = 0	8.6R	7.5R
k = 0.333	7.8R	7R
k = 1	5.6R	5R
k = 3	7.4R	7R

Table 2.

Influence zone extends around circular and horseshoe sections for different stress fields.

R is the radius of equivalent circular section [29, 30].

In all cases presented in this paragraph, we assumed the initial stress field p = 1 MPa, modulus of elasticity E = 1GPa and Poisson ratio ν = 0.3. Using these values we are able to calculate the corresponding values for any p value and every kind of rock. The only constraint is that the value of Poisson ratio is equal to 0.3. In the following paragraph, we shall discuss the influence of ν upon the stresses.

2.2 The case of shallow tunnels: the circular section

Assuming that the tunnel is close to the ground surface, we can observe an influence of the boundary as the vertical stresses depend on the depth. So, the initial stresses at the roof and the bottom levels of the tunnel section are not equal. Bray [31] presented an extensive study on this problem. Some practical rules are summarized below.

By assuming an acceptable deviation, as in the case of the influence zone estimation, we may distinct three cases.

2.2.1 Case I

The tunnel depth, measured from the section’s centre, is Z > 25R, where R is the tunnel radius. In this case Kirsch equations apply. In Figure 9 , the distribution of stresses around a circular deep tunnel is shown [32]. There is no influence of the boundaries. Note that only a window of the stress field is shown.

Figure 9.

Stress distribution around a circular deep tunnel (Z = 25R, k = 1) [32].

2.2.2 Case II

The tunnel depth is 7R < Z < 25R. In this case we must add a correction term to the Kirsch equations depending on the Poisson ratio. That term is given by Bray as [31]:

$${\sigma }_{\theta }=\mathit{\gamma h}\left[1+k+2\left(1-k\right)cos2\theta \right]-\mathit{\gamma \alpha }sin\theta \left[\frac{3-4\nu }{2\left(1-\nu \right)}+2\left(1-k\right)cos2\theta \right]$$ E10

where α is the tunnel’s radius. This expression agrees with an equation given by Savin [33].

In Figure 10 , the distribution of stresses around a circular shallow tunnel is shown. The influence of the surface boundary starts to be noticeable [32].

Figure 10.

Stress distribution around a circular shallow tunnel (Z = 7R, k = 1) [32].

2.2.3 Case III

The tunnel depth is Z < 7R. This case is more difficult because the influence of the boundary becomes greater. Then Mindlin’s closed-form solutions apply [34, 35]. Mindlin obtained solutions to this problem for three cases of in situ stress fields at depth z (remote from the tunnel):

2.2.3.1 Case I

p_z = wz.

p_h = wz.

i.e. isotropic, or hydrostatic, gravitational pressure and w being the unit weight of mass.

2.2.3.2 Case II

p_z = wz.

p_h = {ν/(1-ν)}wz.

This is the case where the lateral deformation is a constraint remote from the tunnel.

2.2.3.3 Case III

p_z = wz.

p_h = 0.

This is the case of nonlateral constraint of the mass remote from the tunnel.

The solutions of Mindlin were in terms of bipolar coordinates α and β. The expression giving the stress on the circular boundary for the first case is [35]:

$$\begin{array}{cc}{\left[{\sigma }_{\beta }\right]}_{\alpha ={\alpha }_1}& =\frac{2w\text{A}\left(\cosh {\alpha }_1-\cos \beta \right)}{\sinh {\alpha }_1}\{\frac{1-\cosh {\alpha }_1\cos \beta }{{\left(\cosh {\alpha }_1-\cos \beta \right)}^2}\hfill \\ & -\coth {\alpha }_1-\frac{\left(7-8\nu \right)}{4\left(1-\nu \right)}\frac{\cos \beta }{\sinh {\alpha }_1}+2e^{-{\alpha }_1}\cos \beta {\displaystyle \sum }_{n=2}^{\infty }{R}_n\cos n\beta \}\end{array}$$ E11

where α₁ is the value of α corresponding to the boundary of the tunnel:

$$\begin{array}{l}\frac{Z}{R}=cosh{\alpha }_1\\ A=Ztanh{\alpha }_1=Rsinh{\alpha }_1\\ R_n=N_n-{\mathit{ne}}^{-n{\alpha }_1}\\ N_n=\frac{{\mathit{ne}}^{-n{\alpha }_1}\left(sinhn{\alpha }_{1}coshn{\alpha }_1-nsinh{\alpha }_{1}cosh{\alpha }_1\right)}{{sinh}^{2}n{\alpha }_1-n^2{sinh}^2{\alpha }_1}\end{array}$$ E12

Poulos and Davis [35] gave in figures and tables values of σ_β for 1 < (Z/R) < 4.

For Cases II and III above, the solution for σ_β is obtained by adding in the solution given by Eq. (11) a further expression [33]. From the above analysis, the importance of taking into account the influence of the ground surface on the stress distribution on tunnel boundary is obvious. A further important notice is the influence of Poisson ratio on the expressions of stresses. From Eq. (10), it is clear that this influence is important when Z < 25R.

Malvern [36] in a rigorous analysis concluded that when the resultant force on each boundary is zero, then the stress distribution is independent of the elastic constants. Otherwise the stress distribution will depend on ν. When the boundary conditions include displacement constraints, then the stress distribution will, also, depend on ν. These conclusions are important to be known for the stress analysis of tunnels by using computational methods like finite elements or boundary elements. Frocht [14] examined the influence of stresses from ν by extensive photoelastic experimental study. In this short presentation of the elastic solutions for the shallow circular tunnel, we restricted to the stress field influence. For an extensive presentation of elastic solutions for the important problem of surface settlements induced by tunneling, which is of great practical significance, see [37, 38, 39]. Another important topic in the case of shallow tunnels is the problem of stress field due to seismic loading. Recently, Pelli and Sofianos [40] published a paper addressing the very important topic of the stress field around shallow tunnels under seismic loading of SW waves.

3. Analytical methods: plastic problems

3.1 The case of deep tunnels

We now turn our interest in plastic problems. This is a case of great theoretical and practical significance, as every underground excavation must be analysed against failure. For a detailed coverage, see, almost, every standard textbook of soil and rock mechanics. For a more detailed coverage, see [41, 42, 43, 44]. Here, we examine special topics related to tunnel engineering. In order to analyse the competence of an excavation in soil or rock, we must extend our analysis beyond elasticity to the plastic behaviour of the surrounding medium. To this end, we have to introduce to our analysis a failure criterion. The most common criterion is the Mohr-Coulomb [11]. Extensive research resulted in a remarkable advancement in understanding the material behaviour and the implementation of advanced models in computational methods. A milestone was the introduction of critical state theory in the 1960s [45]. Further research in rock mechanics resulted in the introduction of the Hoek-Brown criterion, which is widely used in rock engineering practice [13]. Yu [46], in an extensive review, presented a wide range of failure criteria applicable on, almost, every kind of materials under every possible type of loading. In this review Yu discussed as well the soil and rock materials.

Let us examine the problem of practical significance of the formation of plastic zone around a deep tunnel of circular section. Closed-form solutions can be found for every criterion of soil and rock literature in the case of isotropic or hydrostatic stress field (k = 1), because of the axisymmetric type of this case. Recently, Vrakas and Anagnostou presented an extension of the small strain analyzes to obtain finite strain solutions [47]. For the general case of stress fields (k ≠ 1), closed-form solutions have been obtained so far for the Tresca criterion [48] and for the Mohr-Coulomb criterion [49]. Here, the elliptic paraboloid criterion developed by Theocaris [50, 51] is introduced to solve the problem of plastic zone around a circular deep tunnel in rock. First, the mathematical expression of the criterion is given, and a comparison of it with Griffith and Hoek-Brown criteria is presented.

This criterion is expressed through the three principal stresses. It is an energy criterion having as parameter the absolute value of the ratio R of uniaxial compressive strength over the uniaxial tensile strength. The elliptic paraboloid criterion is developed for applications in applied mechanics in general, so the tensile stresses are positive, and the compressive stresses are negative. Then:

$${\left({\sigma }_1-{\sigma }_2\right)}^2+{\left({\sigma }_2-{\sigma }_3\right)}^2+{\left({\sigma }_3-{\sigma }_1\right)}^2+2\left(R-1\right)\left({\sigma }_1+{\sigma }_2+{\sigma }_3\right){\sigma }_{\mathrm{t}}=2R{{\sigma }_{\mathrm{t}}}^2$$ E13

$$\text{where}\text{R}=\frac{\left|{\sigma }_c\right|}{{\sigma }_{t.}}.$$

Eq. (13) is the addition of two components, which express two parts of the total elastic energy. Indeed, the first part is, up to a constant multiplicative factor, the distortion energy expressed through the deviatoric stresses. This part represents the energy of the elastic change of shape [52] given by Eq. (14):

$${\left({\sigma }_1-{\sigma }_2\right)}^2+{\left({\sigma }_2-{\sigma }_3\right)}^2+{\left({\sigma }_3-{\sigma }_1\right)}^2$$ E14

The second part depends on the hydrostatic or spherical part of the stresses and represents the energy of elastic change in volume. In this way, the elliptic paraboloid criterion combines both the change of shape and volume. The latter is important for soil and rock materials as their strength depends on the confining pressure.

For the isotropic case of field stress, conditions of axisymmetry are valid, i.e. σ ₁ = σ ₂. Then from Eq. (13), we obtain:

$${\left({\sigma }_3-{\sigma }_1\right)}^2+\frac{R-1}{R}\left(2{\sigma }_1+{\sigma }_3\right){\sigma }_c-\frac{{\sigma }_c^2}{R}=0$$ E15

In the case of axisymmetric conditions, the criterion becomes paraboloid of revolution. In Figure 11 , a comparison of Griffith criterion with Hoek-Brown criterion being shown. It is obvious that the deviation between their representation for m value being equal to 7.88, close to the restriction of Griffith’s corresponding value being 8 [52].

Figure 11.

Comparison of Griffith criterion with Hoek-Brown criterion for different m values [52].

Now, we can proceed to a comparison of elliptic paraboloid criterion with Griffith. We assume an R value of 8 ( Figure 12 ).

Figure 12.

Comparison of elliptic paraboloid and Griffith criteria (R = 8) [52].

From the above figures, we may conclude that Griffith and elliptic paraboloid criteria agree, although their assumptions are completely different. Griffith assumed that failure initiates from pre-existed cracks because of the stress concentration at the tips of the cracks. This assumption is fundamental in fracture mechanics. On the other hand, the assumption of elliptic paraboloid criterion is an extension of von Mises criterion. Both criteria are identical in the case of equal strengths under compression and tension (R = 1). This assumption is close to the experimental results for metals. Theocaris applied this criterion in igneous and metamorphic rocks [45].

We, now, come to examine the problem of plastic zone formation around a deep tunnel of circular section under isotropic stress field (k = 1) [53] ( Figure 13 ).

Figure 13.

The geometry of the problem. The plastic zone around the circular tunnel is shown [53].

To solve the problem, we introduce the equation of equilibrium in polar coordinates (Eq. (16)) and the flow rule (Eq. (17)). Then, by introducing the expression of elliptic paraboloid criterion in Eq. (17), we obtain Eq. (18) [53]:

$$\frac{d{\sigma }_r}{\mathit{dr}}+\frac{{\sigma }_r-{\sigma }_{\theta }}{r}=0$$ E16

$$\frac{\partial f}{\partial {\sigma }_2}=0$$ E17

$$\frac{\partial \left[{\left({\sigma }_1-{\sigma }_2\right)}^2+{\left({\sigma }_2-{\sigma }_3\right)}^2+{\left({\sigma }_3-{\sigma }_1\right)}^2+2\cdot \left(R_b-1\right)\cdot \left({\sigma }_1+{\sigma }_2+{\sigma }_3\right)\cdot {\sigma }_{\mathit{tb}}-2\cdot R_b\cdot {{\sigma }_{\mathit{tb}}}^2\right]}{\partial {\sigma }_2}=0$$ E18

After complicated algebraic manipulations, we obtain the condition for plastic zone formation around the tunnel (Eq. (19)) [53]:

$$p_o<\frac{\left(R_i-1\right)\cdot {\sigma }_{\mathit{ti}}}{2}-p_i-\sqrt{-\left(R_i-1\right)\cdot {\sigma }_{\mathit{ti}}\cdot p_i+\frac{1}{3}\cdot R_i\cdot {{\sigma }_{\mathit{ti}}}^2+\frac{{\left(R_i-1\right)}^2\cdot {{\sigma }_{\mathit{ti}}}^2}{3}}$$ E19

Finally, the radius of the plastic zone r_c is obtained (Eq. (20)) [53]:

$$r_e=r_i\cdot e^{\frac{L-K}{2\cdot \left(R_b-1\right)\cdot {\sigma }_{\mathit{tb}}}-\frac{1}{2}\cdot ln\frac{K}{L}}$$ E20

L and K above are constants given by Eq. (21), (22):

$$K=-\left(R_b-1\right)\cdot {\sigma }_{\mathit{tb}}-\sqrt{\frac{4}{3}\cdot R_b\cdot {{\sigma }_{\mathit{tb}}}^2-4\cdot \left(R_b-1\right)\cdot {\sigma }_{\mathit{tb}}\cdot {\sigma }_{\mathit{re}}+\frac{4\cdot {\left(R_b-1\right)}^2\cdot {{\sigma }_{\mathit{tb}}}^2}{3}}$$ E21

$$L=-\left(R_b-1\right)\cdot {\sigma }_{\mathit{tb}}-\sqrt{\frac{4}{3}\cdot R_b\cdot {{\sigma }_{\mathit{tb}}}^2-4\cdot \left(R_b-1\right)\cdot {\sigma }_{\mathit{tb}}\cdot p_i+\frac{4\cdot {\left(R_b-1\right)}^2\cdot {{\sigma }_{\mathit{tb}}}^2}{3}}$$ E22

Hoek and Brown [9] obtained the radius of the plastic zone for Hoek-Brown criterion as follows [13]:

$${\mathrm{r}}_{\mathrm{e}}={\mathrm{r}}_{\mathrm{i}}{\mathrm{e}}^{\left\{\mathrm{N}-\frac{2}{{\mathrm{m}}_{\mathrm{r}}{\mathrm{\sigma }}_{\mathrm{c}}}{\left({\mathrm{m}}_{\mathrm{r}}{\mathrm{\sigma }}_{\mathrm{c}}{\mathrm{p}}_{\mathrm{i}}+{\mathrm{s}}_{\mathrm{r}}{\mathrm{\sigma }}_{\mathrm{c}}^2\right)}^{\frac{1}{2}}\right\}}$$ E23

$$\mathrm{N}=\frac{2}{{\mathrm{m}}_{\mathrm{r}}{\mathrm{\sigma }}_{\mathrm{c}}}\left({\mathrm{m}}_{\mathrm{r}}{\mathrm{\sigma }}_{\mathrm{c}}{\mathrm{p}}_0+{\mathrm{s}}_{\mathrm{r}}{\mathrm{\sigma }}_{\mathrm{c}}^2-{\mathrm{m}}_{\mathrm{r}}{\mathrm{\sigma }}_{\mathrm{c}}^2\mathrm{M}\right)$$ E24

$$\mathrm{M}=\frac{1}{2}{\left({\left(\frac{\mathrm{m}}{4}\right)}^2+\mathrm{m}\frac{{\mathrm{p}}_0}{{\mathrm{\sigma }}_{\mathrm{c}}}+\mathrm{s}\right)}^{\frac{1}{2}}-\frac{\mathrm{m}}{8}$$ E25

In Table 3 a comparison of the above criteria is shown. We can conclude that, except for low values of in situ stresses, both criteria are close in their predictions of radial stresses at the elastic–plastic interface. On the contrary, their predictions regarding the extent of the plastic zone differ with Hoek-Brown criterion being somehow more conservative.

Elliptic paraboloid			Hoek-Brown
P0 (MPa)	σre (MPa)	re (m)	σre (MPa)	re (m)
4	2.10	10.25	1.96	13.45
3	1.35	8.29	1.31	11.22
2	0.64	6.55	0.71	9.07
1.5	0.31	5.76	0.44	7.99
1	0.01	5.02	0.20	6.84

Table 3.

Comparison of elliptic paraboloid and Hoek-Brown criteria for σ_c = 1 MPa, σ_t = 0.125 MPa, R = 8.00, m = 8 and s = 1 [53].

3.2 The shallow tunnel case

The shallow tunnel case is complicated because of the influence of the proximity of ground surface on the stress field and the influence of gravity. As we notice for the case of deep tunnels, there are closed-form solutions for the plastic zone formation around circular tunnels for isotropic stress field and for the general case of Tresca [48] and Mohr-Coulomb [49] criteria. For the case of shallow tunnels, there were no closed-form solutions until 2009 when a solution was published for the Mohr-Coulomb solution under the assumption of isotropic stress field (k = 1) and no gravitational stresses, based on bipolar coordinates ( Figure 14 ).

Figure 14.

Bipolar coordinates used for the solution of shallow tunnel problem [54].

In the following the solution published in [54] and further applied [55] is presented omitting the detailed mathematical analysis:

$$\begin{array}{cc}{P}_{\mathit{cr}}& =\frac{2{\kappa }^2}{2\left(d_i^2-r_i^2{cos}^2\beta \right)+{\kappa }^2\left(\lambda -1\right)}\left[P_0\left(\frac{d_i^2-r_i^2{cos}^2\beta }{{\kappa }^2}\right)-\frac{Y}{2}\right]\\ & =\frac{2{\kappa }^2}{2\left({\kappa }^2+r_i^2{sin}^2\beta \right)+{\kappa }^2\left(\lambda -1\right)}\left[P_0\left(\frac{{\kappa }^2+r_i^2{sin}^2\beta }{{\kappa }^2}\right)-\frac{Y}{2}\right]\end{array}$$ E26

where r_i is the tunnel radius and d_i (= κcothα_i) is the depth of the tunnel axis from the surface. Considering the conditions of the problem, following [54] the final form of the equation giving the shape of plastic zone is given below:

$${\left(\frac{r_c}{r_i}\frac{d_i-r_{i}cos\beta }{d_c-r_{c}cos\beta }\right)}^{1-\lambda }=\frac{\left[2M_0+{\kappa }^2\left(\lambda -1\right)\right]\left[Y+P_i\left(\lambda -1\right)\right]}{2M_0\left[Y+P_0\left(\lambda -1\right)\right]}$$ E27

where $M_0={\kappa }^2+r_c^2{sin}^2\beta$ _.

In Figure 15 , a parametric study of the plastic zone shape is shown [54].

Figure 15.

Plastic zone formation for different P_i/P₀ values for P₀ = 500 kPa, c = 100 kPa, φ = 25⁰, tunnel’s depth = 20 m and r = 5 m [54].

Eq. (27) can be solved, also, by using MATLAB [56]. In Figure 16 , an example is shown.

Figure 16.

Plastic zone around shallow tunnel with MATLAB. Example SMG1 [54]: D = 10 m, r = 5 m, γ = 25 kN/m³, P₀ = 250 kPa, P_i = 50 kPa, c = 60 kPa, φ = 25⁰.

The analytical solution of this problem apart from its theoretical interest could be used in conjunction with computational analysis in practical problems [57]. In [55] a very important case is presented where the task was for a new tunnel to underpass the monumental Chandpole Gate in Jaipur, India. From the closed-form solution, the critical internal pressure was calculated to start with a further computational search for the optimum EPBM pressure in order to minimize the surface settlements induced by the excavation. For a more deep analysis of this problem, it seems that we have to proceed using semi-analytical methods. For the current status of this research, see [58, 59]. Important developments have been presented by Schofield [60] and Mair [61] in their Rankine Lectures of 1980 and 2008, respectively.

4. Applications of computational methods in tunnel engineering

In the last section of the chapter, a short, not complete, survey of computational methods in tunnel engineering is presented. It is a rather chronological account of methods and tools based on author’s personal experience.

To start with, in Figure 17 a very useful and clear classification of methods of analysis is presented [62, 63, 64]. For an extensive review of numerical analysis, see Potts [65].

Figure 17.

The four basic methods, in two levels, comprising eight different approaches to rock mechanic modeling [62, 63, 64].

4.1 Level 1: basic numerical methods: 1:1 mapping

The available numerical methods (FEM, BEM, DEM) belong to category “C” making “one-to-one mapping”. The meaning of this term is that they make a direct modeling of geometry and physical mechanisms [63]. In the 1970s finite element method (FEM) codes, based on differential equation formulation, were developed running in mainframe computers. Towards the end of the 1970s, the boundary element method (BEM) was developed, based on integral equation formulation. Bray and his co-workers introduced a 2D indirect formulation of BEM, and they developed a code included in [13].

Based on the example presented in Hoek and Brown ([13], p. 499), a twin cavern problem was analysed using the indirect BEM formulation [66]. In order to check the numerical analysis, a photoelastic model was, also, analysed. The BEM code was modified, to plot isochromes and isoclinics [14] obtained from photoelasticity, for comparison. In Figure 18 , plotting of isochromes and those obtained from the experiment is shown. The agreement between the two is remarkable [66].

Figure 18.

Plotting of isochromes of the twin cavern case. (a) Curves obtained from the numerical analysis and (b) curves obtained from the experiment [66].

In a later stage (2002), with the advances in information systems and methods, the numerical methods to solve engineering problems advanced exponentially. A 3D FEM code was developed, and an extensive study was undertaken for the modeling of rock mass and underground excavations [52]. In that code several failure criteria were implemented. It may be the only 3D FEM code running the elliptic paraboloid criterion described in Section 3.1. In order to benefit from the 3D capabilities of the code, a more complicated problem was analysed. The cavern section included in [13] was modified to include, also, a tunnel of circular section cutting the cavern at a right angle ( Figures 19 and 20 ). The rock mass is assumed to be anisotropic with two families of joints [51]. In formulating the problem, the following assumptions have been made: modulus of elasticity E = 7 GPa, ν = 0.25, γ = 25 kN/m³, tunnel’s depth Z = 400 m, vertical stress at infinity p = 10 MPa and in situ stress ratio k = 0.8. The model had 1441 nodes with 4323 degrees of freedom. About 236 isoparametric hexahedral elements with 20 nodes were used.

In the following figures, some characteristic results are presented. In Figure 19 , the formation of plastic zones around the tunnels in the conjunction area is presented. Because of the anisotropy of the problem, the plastic zones are not symmetric. In Figure 20 , the convergences of the excavation boundaries are shown. Finally, in Figure 21 , the displacement and stress contours are shown. The r and s coefficients are defined as r = E’/E and s = G’/G. All excavations are assumed to be unlined.

Figure 19.

Plastic zone formation around the excavations according to elliptic paraboloid failure criterion [52].

Figure 20.

Convergences at the excavation boundaries [52].

Figure 21.

Displacement contours (k = 0, β = 30°, r = 0.05, s = 0.013) and σ₃ stress contours (k = 0, β = 30°, r = 0.60, s = 0.91) [52].

4.2 Level 2: system approaches: non-1:1 mapping

In the last paragraph of this chapter, the application of artificial intelligence (AI) methods and tools in tunnel engineering is shortly discussed. This category of methods belongs to Level 2 methods achieving a non-1:1 mapping ( Figure 17 ). In this category of methods, the rock or soil mass is mapped indirectly by a network of nodes [67]. In the 1990s some applications of AI in geotechnical engineering were published [68, 69, 70, 71, 72, 73, 74]. From 2001 a great number of application of AI methods in geotechnical engineering have been published. There are two categories of approaches: the supervised learning and the unsupervised learning. Backpropagation, which was the first method applied in geotechnical engineering, belongs to the first category. Interconnected nodes corresponding to parameters involved in the problem represent the physical problem. The output of the training process is taken as a known target [67]. On the contrary, in unsupervised learning methods, the system has to extract knowledge from the data resulting in the underlying interconnections of the parameters of the problem [75]. Currently, a great number of publications are available using more advanced and sophisticated AI methods. A critical factor is the quality of data and the engineering judgment of the users. There are cases where the overtraining of the system resulted in a “blind” learning of the data guiding to wrong conclusions. The ability to “include creative ability, perception and judgment” [67] is, still, not achieved.

Here, an attempt to train a backpropagation network system using surface settlements induced by underground excavations is presented [76]. The total amount of the available data is 90 records, coming from different types of ground profiles and referring to tunnels constructed with different methods of excavation [77]. The provided information includes tunnel size and depth (Zo), maximum settlement (w), settlement trough width (i), volume (Vs), ground description, geological properties and method of working. An indicative part of the available data is given in Table 4 .

Tunnel	Tunnel data		Maximum recorded settlement	Ground geotechnical properties		Volumes	Settlement trough width	Tunneling method and soil conditions
	Depth zo (m)	2R (m)	w (mm)	cu (kN/m^2)	γzo/cu	Vs (%)	i (m)
London transport fleet line green park	29.30	4.15	6.17	270	2.1	1.40	12.6	Shield construction; Stiff, fissured overconsolidated London clay
NWA sewerage scheme, Hebburn	7.50	2.01	7.86	75	2.0	2.42	3.9	Shield construction; Laminated clay overlain by stony clay
NWA sewerage scheme, Willington Quay Syphon, Contract 32	13.37	4.25	81.50	33	9.4	13.10	9.1	Shield construction; Silty

Table 4.

Indicative examples of case history data [76, 77].

Case history data on some tunnels and tunneling conditions [76, 77].

The range of values of the data, the type of geological profiles and the excavated methods are given in Table 5 [77].

	zo (m)	2R (m)	w (mm)
Min	3.40	1.26	1.50
Max	37.49	30.50	280.00
Geology	Homogeneous formations
	Heterogeneous formations
Excavation method	Shield construction
	Hand excavation
	Slurry shield
	Shieldless excavation
	Full face tunneling machine
	Hydro-shield
	Full face blasting

Table 5.

Range of values of the available data and the categories of geological profiles and excavation methods [76, 77].

In training the ANN, two approaches followed. In the first approach, all the available data were included. The tunnel depth and diameter have been used as input data and the surface settlement as output. Qualitative “data” as the geological conditions and the excavation method were not initially taken into consideration as parts of the training process. The training was not successful because we mixed the nonhomogeneous data. Therefore we proceeded to the second approach using all available data and information. To this end, we assigned a number to distinct different geological conditions and excavation method. In this training four inputs were used and the training was successful. In Figure 22 , the relative importance of the parameters involved to the problem is shown. From this figure, we may conclude that the geometry of the tunnel, i.e. radius and depth, is the main factor with geological profile and excavation method having an important contribution too [76].

Figure 22.

Relative importance of the parameters contributing to induced settlements due to tunneling [76].

Further tunneling applications using AI tools can be found in the literature. In Zaré and Lavasan [78] an objective system approach is adopted to quantify the interaction of parameters involved in the problem of tunnel face stability. The method is based on a backpropagation ANN approach. A different approach of objective system approach, which adopted the unsupervised type of learning based on self-organizing maps, is, also, published with applications in rock engineering problems [75, 79, 80, 81, 82, 83, 84, 85].

5. Conclusions

In this chapter, a rather subjective view of tunnel engineering is presented. Nevertheless, the progress made in analytical and computational methods was followed with an effort to be well documented. Starting from a well-known problem of the elastic stress field around a circular deep tunnel, we investigated the influence of tunnel’s shape in stress and displacement fields around a tunnel of the horseshoe section, as well as the extent of the influence zone for several stress fields. To this end Muskhelishvili’s complex variable formulations of stress functions were used. The next problem was the case of the shallow circular tunnel for which established elastic solutions are presented. Proceeding to the more difficult problem of the plastic analysis, we examined the case of deep tunnels, and then we presented a closed-form solution of the plastic zone formation for the shallow circular tunnel. This is a topic still needing further investigation because of its mathematical difficulty. Computational methods were the last part of the chapter. A classification of methods was presented followed by a problem of deep tunnels analyzed using 3D finite element analysis. The increasing exploitation of artificial intelligence tools in analyzing geotechnical problems was the last topic. The presentation was based on the tunneling-induced surface settlements as a paradigm. The cited references listed below are a basic and indicative selection of literature for further reading.

Acknowledgments

The author acknowledges the contribution of his former PhD, MSc and MEng students at the National Technical University. Their theses are listed in the references, and they are referred in the text.