Stress solution for an arbitrarily shaped tunnel with a neighboring circular cavity

Abstract

This paper gives an analytical stress solution for the plane strain problem that the deep rock mass contains an arbitrarily shaped tunnel and a circular cavity subjected to in situ stresses. Based on the discovery that the mapping function for an arbitrarily shaped hole can also manage a tunnel and a circular hole mapping into two circles, the plate bounded by the two holes is conformally mapped onto an annulus. The complex variable method is used to establish the stress boundary condition, and then the two analytical functions are solved by the power series method. Comparison with ANSYS solution is made to validate the correctness of deductions and solving process. To consider practical underground engineering, the cases of different lateral stresses, distances of the two holes, sizes, and positions of the cavity are discussed.

1. Introduction

Tunnel is a kind of underground structure built to utilize the underground space, which is widely seen in hydropower engineering, transportation engineering, and some other types of engineering. It is necessary to conduct a fundamental stress analysis of the tunnel for its functioning safety. A majority of tunnels are long and straight, and the boundary condition changes little along the tunnel axis, because of which it is treated as a plane strain problem. Furthermore, the influence of ground boundary can be ignored while the tunnels are deeply buried so that it is an infinite field problem.

However, the early studies on this kind of plane structure encountered problems that the conventional methods are not applicable because of the irregular tunnel shape. It was not until Muskhelishvili [1] proposed the complex variable method that the problem of an arbitrarily shaped tunnel was well solved. Since twin tunnels are typically closely arranged in parallel, their stress fields interact with each other which makes the mathematical form of the stress function extremely complicated. At present, only the problems of two holes with special shapes and arrangement have been solved analytically.

Howland and Knight [2] derived the solution for the cases of several equal circular holes. Green [3] presented the solving method for arbitrary circular holes, and Hoang and Abousleiman [4] acquired the solution for two unequal circular holes, which they both took advantage of biharmonic functions. Ling [5], Iwaki and Miyao [6] and Radi [7] used the bipolar coordinate method, which is proposed by Jeffery [8], to solve the cases of two circular holes. Ling solved for two holes of equal size under remote tension, while Iwaki and Miyao solved for two unequal holes with additional stress on the hole boundaries. And Radii used the obtained solution to evaluate the translational material forces and the energy release rate. The two methods mentioned here are only suitable for circular holes.

Davies and Hoddinott [9] first used the complex variable method to obtain the approximate solution of two unequal circular holes, in which they added a solution to the single hole solution. Haddon [10] also took advantage of the complex variable method to solve the same problem, but he acquired an analytical stress solution by simultaneously establishing the boundary conditions of two holes. The superiority of the complex variable method over others lies in the conformal mapping technique, with which Lu et al. [11] derived the solution for the case of an elliptical hole and a circular hole. In 2018, Zeng et al. [12] obtained the solution for two elliptical holes, the following year Zeng and Lu [13] found a general form of mapping function for two arbitrarily shaped holes, each of which should be symmetrical about the line joining their centers and successfully derived the stress solution. Although the complex variable method shows its capability for solving two non-circular holes, there is no solution for two arbitrarily shaped tunnels until now because Zeng and Lu’s solution is confined in two holes having a common symmetrical axis.

The only method existing for two arbitrarily shaped tunnels is the Schwartz alternating method. First, some scholars used it to solve for two circular holes and two elliptical holes [14–16], but it is too complicated to be used for problems of two arbitrarily shaped tunnels until Zhang and Lu [17] improved the method. They expanded the redundant stress in the Fourier series to solve for two U-shaped tunnels, and then Zhang et al. [18] and Zhang et al. [19] used the improved method to solve for multiple elliptical holes. Nevertheless, this method is still extraordinarily complicated to obtain a high-precision solution.

In this paper, the analytical stress solution for an arbitrarily shaped tunnel and an adjacent circular cavity, which are deeply buried and subjected to in-situ stresses, is presented using conformal mapping technique and complex variable method. The tunnel and cavity are mapped into an annulus simultaneously, which is accomplished by a mapping function for a single arbitrarily shaped hole. Then, the stress boundary conditions are established by the complex variable method and stress is obtained by calculating out the two analytical functions.

2. Statement of the problem

An arbitrarily shaped tunnel is placed underground deep and horizontally, whose axis is parallel to the z-axis, and a circular cavity is parallel to the tunnel (see Figure 1). The coordinate origin is located in the center of the tunnel and the x-axis coincides with the line connecting the centers of the two holes. Since the two cavities are deeply buried, it can be assumed that the stresses are uniform, which the horizontal stress is P and the vertical stress is Q.

Figure 1.

Cross-section of a tunnel and a circular hole in an infinite and elastic rock mass.

In Muskhelishvili’s complex variable method, stress function for plane elastic problem can be expressed by two analytical functions, which are ${\phi }_1$ (z) and ψ₁(z) in this paper. Then, stresses can be obtained by the following equations:

$${\sigma }_x+{\sigma }_y=4\mathrm{Re}\left[{{\phi }^{\prime }}_1\left(z\right)\right],$$

(1)

$${\sigma }_y-{\sigma }_x+2i{\tau }_{xy}=2\left[\overline{z}{{\phi }^{″}}_1\left(z\right)+{{\psi }_1}^{\prime }\left(z\right)\right].$$

(2)

The general stress boundary condition is mostly expressed in terms of integration:

$${\phi }_1\left(z\right)+z\overline{{\phi \prime }_1\left(z\right)}+\overline{{\psi \prime }_1\left(z\right)}=i{\int }_L({\overline{X}}_n+i{\overline{Y}}_n)\mathrm{ds}+c_A.$$

(3)

The rock mass should lie left while integrating. ${\overline{X}}_n$ and ${\overline{Y}}_n$ are the given stress components on the hole boundary in x-directions and y-directions, respectively. Usually, ${\overline{X}}_n$ and ${\overline{Y}}_n$ are the components of water pressure; they are equal to zero when there is no water in the two cavities. If the water pressure is P₁ in the tunnel and P₂ in the circular cavity, the stress boundary conditions of the two holes can be denoted as:

$${\phi }_1\left(z\right)+z\overline{{\phi \prime }_1\left(z\right)}+\overline{{\psi \prime }_1\left(z\right)}=P_{1}z,z\in L_1,$$

(4)

$${\phi }_1\left(z\right)+z\overline{{\phi \prime }_1\left(z\right)}+\overline{{\psi \prime }_1\left(z\right)}=P_{2}z+c_A,z\in L_2.$$

(5)

For constant does not contribute to the stress, $c_A$ is set as zero for tunnel boundary. And $c_A$ is set as an unknown quantity for the circular cavity boundary to discuss the general case.

As to the problem considered in this paper, the two analytical functions have the following forms:

$${\phi }_1\left(z\right)=\mathrm{\Gamma }z+{\phi }_0\left(z\right),$$

(6)

$${\psi }_1\left(z\right)=\mathrm{\Gamma }\prime z+{\psi }_0\left(z\right),$$

(7)

where $\mathrm{\Gamma }=(P+Q)/4$ , $\mathrm{\Gamma }\prime =(Q-P)/2$ , and ${\phi }_0(z)$ and ${\psi }_0(z)$ are single-valued analytical functions in the discussing region.

3. Conformal mapping

It is impossible to calculate out the two analytical functions from the stress boundary conditions in their current form. Conformal mapping should be conducted to simplify the problem. Referring to our earlier research [13], two steps are needed to map the region (shown in Figure 1) into an annulus.

In the first step, the infinite plate with a tunnel and a circular hole (see Figure 2(a)) should be mapped into an infinite plate with two circular holes (see Figure 2(b)). However, the form of the mapping function we proposed is no longer suitable for this problem because there does not exist a common symmetrical axis. Through our research, we found that the mapping function which maps an arbitrarily shaped hole into a circular hole can manage the discussed problem here, whose form is as follows:

$$z={\omega }_1(\zeta \prime )={\lambda }_1\zeta \prime +{\lambda }_0+\sum _{k=1}^{\infty }{\lambda }_{-k}\zeta {\prime }^{-k}.$$

(8)

Figure 2.

First step of the conformal mapping. (a) z plane with a tunnel and a circular hole and (b) $\zeta \prime$ plane with two circular holes.

For almost all the cases of single tunnel problems, the coefficients of the above formula, namely, λ₁, λ₀, and λ _−k (k = 1, 2, 3, …), are the real numbers or pure imaginary numbers. Without a common symmetrical axis for the two holes, the coefficients are complex numbers in this paper. Although the mapping function is much more complicated, it also can be figured out by optimization methods while λ _−k is truncated at a countable number k.

Set the radius of the image circle mapped from the tunnel be 1. Let the radius of the image circle mapped from the cavity be r₂ and the distance between the edges of the two circles be S₁ (see Figure 2(b)). r₂, S₁, λ₁, λ₀, and λ _-k are undetermined which can be found out through the correspondence between the coordinates before and after the mapping, and the solving process is as follows.

Select m points on each hole boundary in Figure 2(a) (the inflection point should be included), assume that they are, respectively, mapped into points $\zeta {\prime }_1,\zeta {\prime }_2,\dots ,$ and $\zeta {\prime }_{2m}$ in the image plane, which can be denoted by angles ${\beta }_1,{\beta }_2,\dots ,$ and ${\beta }_{2m}$ in local coordinate (see Figure 2(b)). ${\beta }_1,{\beta }_2,\dots ,$ and ${\beta }_{2m}$ are unknowns too, then there are 2 m + k + 4 unknowns including unknowns r₂, S₁, λ₁, λ₀, and λ _-k . Since they are hardly solved through equality constraints directly, we use the technique of least square method embedding mixed penalty function method to find them out. For the points mapped before and after having the correspondence relationship of equation (8), one obtains:

$${\lambda }_0+{\lambda }_1\zeta {\prime }_i+\sum _{k=1}^{\infty }{\lambda }_{-k}\zeta {\prime }_i^{-k}=z_i(i=1,2,\dots ,2m).$$

(9)

Let:

$$A=\left(\begin{array}{c}1{\zeta \prime }_1{\zeta }_1^{\prime -1}{\zeta }_1^{\prime -2}\dots {\zeta }_1^{\prime -k}\hfill \\ \multicolumn{1}{c}{1{\zeta \prime }_2{\zeta }_2^{\prime -1}{\zeta }_2^{\prime -2}\dots {\zeta }_2^{\prime -k}}\\ \multicolumn{1}{c}{⋮⋮⋮⋮\ddots ⋮}\\ \multicolumn{1}{c}{1{\zeta \prime }_{2m}{\zeta }_{2m}^{\prime -1}{\zeta }_{2m}^{\prime -2}\dots {\zeta }_{2m}^{\prime -k}}\end{array}\right),B={({\lambda }_0{\lambda }_1{\lambda }_{-1}{\lambda }_{-2}\dots {\lambda }_{-k})}^T,\mathrm{and}Z={(z_1{z}_2\dots z_{2m})}^T.$$

(10)

$\zeta {\prime }_1,\zeta {\prime }_2,\dots ,$ and $\zeta {\prime }_{2m}$ can be expressed by r₂, S₁ ${\beta }_1,{\beta }_2,\dots ,$ and ${\beta }_{2m}$ ,

$$\{\begin{array}{c}{\zeta \prime }_i=e^{i{\beta }_i}(i=1,2,\dots ,m)\hfill \\ \multicolumn{1}{c}{{\zeta \prime }_i=1+S_1+r_2{e}^{i{\beta }_i}(i=m,m+1,\dots ,2m)}\end{array}.$$

(11)

Noticing that equation (9) is a set of linear equations with the form $A_{2m\times (k+2)}{B}_{k+2}=Z_{2m}$ , so if $\zeta {\prime }_1,\zeta {\prime }_2,\dots ,$ and $\zeta {\prime }_{2m}$ are given, one can find out the coefficients of mapping function, namely, $B$ , by the least square method in the complex domain. If the objective function ( $F_1$ ) of the least square method is the square of relative error of observations $z_i$ and estimates ${\hat{z}}_i$ ,

$$F_1=\sum _{i=1}^{2m}{\left[\frac{|z_i-{\hat{z}}_i|}{\left|z_i\right|}\right]}^2.$$

(12)

B can be figured out by:

$$B={\left({\overline{A}}^T\mathrm{WA}\right)}^{-1}{\overline{A}}^T\mathrm{WZ},$$

(13)

where

$$W=\left[\begin{array}{cccc}\frac{1}{z_1{\overline{z}}_1}\hfill & 0\hfill & \cdots \hfill & 0\hfill \\ \multicolumn{1}{c}{0}& \frac{1}{z_2{\overline{z}}_2}\hfill & \cdots \hfill & 0\hfill \\ \multicolumn{1}{c}{⋮}& ⋮\hfill & \ddots \hfill & ⋮\hfill \\ \multicolumn{1}{c}{0}& 0\hfill & \cdots \hfill & \frac{1}{z_{2m}{\overline{z}}_{2m}}\hfill \end{array}\right].$$

(14)

For r₂, S₁ ${\beta }_1,{\beta }_2,\dots ,$ and ${\beta }_{2m}$ are the design variables, $A$ is known in each optimization process, then the solution $\hat{B}$ of the linear equations can be derived directly by equation (13).

However, linear equations may have no solution, the solution derived by the least square method is approximate, namely, that the equality constraints may do not hold. Hence, we use the penalty function method to optimize design variables X to make all of the equality constraints hold, where $X=(r_2,S_1,{\beta }_1,{\beta }_2,\dots ,{\beta }_{2m})$ . And the objective function in the penalty function method is:

$$F\left(X\right)=\sum _{j=1}^{2m}{(\sum _{i=1}^{k+2}{A}_{\mathrm{ji}}{\hat{B}}_i-z_j)}^2.$$

(15)

The angles in local coordinates should satisfy the following inequality constraints which need to be included in the optimization:

$$\begin{array}{c}\mathrm{s}.\mathrm{t}.G_i={\beta }_{i+1}-{\beta }_i>0,(i=1,2,...,m-1),\hfill \\ G_{i+m-1}={\beta }_{i+m+1}-{\beta }_{m+i}>0,(i=1,2,...,m-1),\hfill \\ G_{2m-1}=2\pi -{\beta }_m+{\beta }_1>0,\hfill \\ G_{2m}=2\pi -{\beta }_{2m}+{\beta }_{m+1}>0.\hfill \end{array}$$

(16)

Add the inequality constraints penalty term to the objective function to obtain the penalty function:

$$P(X,r_f)=\sum _{j=1}^{2m}{(\sum _{i=1}^{k+2}{A}_{\mathrm{ji}}{\hat{B}}_i-z_j)}^2+r_f\sum _{j\in I_1}\ln \frac{1}{G_j\left(X^{\left(0\right)}\right)}+\sum _{j\in I_2}{\frac{1}{r_f}\{min[0,G_j\left(X^{\left(0\right)}\right)]\}}^2,$$

(17)

where $I_1=\{j|G_j(X^{\left(0\right)})>0;1\le j\le 2m\}$ , $I_2=\{j|G_j(X^{\left(0\right)})\le 0;1\le j\le 2m\}$ , $X^{\left(0\right)}$ is the initial value, penalty factor $r_f$ is a set of determined decreasing positive values, and we construct a penalty function $P(X,r_f)$ for every penalty factor $r_f$ . The steps of optimization calculation are as follows:

1. Give initial values of $X^{\left(0\right)}$ , $\epsilon$ , and ${r_f}_{\left(1\right)}$ , set i = 1.

2. Use Davidon-Fletcher-Powell algorithm to minimize $P(X,{r_f}_{\left(i\right)})$ and obtain $X({r_f}_{\left(i\right)})$ .

3. If $|X_i-X_{i-1}|\le \epsilon$ , output $X_i$ and $F\left(X_i\right)$ , otherwise go to Step 4.

4. Let $r_f={r_f}_{(i+1)}$ , set i = i + 1, go to Step 2.

The equality constraints (equation (9)) will hold through optimizing $X$ , and the coefficients of the mapping function will be found synchronously. Figure 3 is an example of the mapping function obtained by the technique. A total of 30 points are selected on each boundary and the mapping function retains 12 negative power terms, and the obtained coefficients are shown in Table 1.

Figure 3.

Contrast between actual hole shape and mapped hole shape.

Table 1.

Coefficients of the mapping function for holes in Figure 3.

r 2	0.68970	λ −3	−0.24210 + 0.03580i	λ −9	−0.00344 + 0.00171i
S 1	0.58611	λ −4	0.01903 + 0.09437i	λ −10	−0.00352 + −0.00631i
λ 0	−0.00237 −0.25385i	λ −5	−0.00895 + 0.00178i	λ −11	−0.00729 + 0.00388i
λ 1	2.48055 + 0.12137i	λ −6	0.00635 + 0.02231i	λ −12	0.00262 + 0.00406i
λ −1	−0.23681 + 0.01295i	λ −7	0.01576 −0.00555i
λ −2	−0.01742 −0.18554i	λ −8	−0.00603 −0.01452i

The mapping function maps Figure 2(b) into Figure 2(a). The giving actual hole shape is in agreement with the mapped shape, which is shown in Figure 3.

Then, the first step of mapping is completed. The second step is to map the two circles into a ring (see Figure 4), and the function is known to be:

$$\zeta \prime ={\omega }_2(\zeta )=\frac{\zeta -d}{d\zeta -1},$$

(18)

where $d=(1+{\gamma }_1{\gamma }_2+\sqrt{({\gamma }_1^2-1)({\gamma }_2^2-1)})/({\gamma }_1+{\gamma }_1)$ and ${\gamma }_1=1+S+2r_2$ , ${\gamma }_2=1+S$ . The radius of the inner circle is $R_0=({\gamma }_1{\gamma }_2-1-\sqrt{({\gamma }_1^2-1)({\gamma }_2^2-1)})/({\gamma }_1-{\gamma }_2)$ and the radius of the outer circle is 1. The angles (θ₁ and θ₂) of the two circles are anticlockwise measured from the positive direction of the ξ-axis to keep the same measuring direction with the angles (α₁ and α₂) of holes in Figure 2(a), respectively.

Figure 4.

An annulus in image plane.

Substituting equation (18) into equation (8) obtains the complete mapping function:

$$z=\omega \left(\zeta \right)={\lambda }_1\frac{\zeta -d}{d\zeta -1}+{\lambda }_0+\sum _{k=1}^{\infty }{\lambda }_{-k}{\left(\frac{\zeta -d}{d\zeta -1}\right)}^{-k}.$$

(19)

4. Solving the two analytical functions

Once z plane is transformed into ζ plane, the two single-valued analytical functions ${\phi }_0(z)$ and ${\psi }_0(z)$ can be written as:

$${\phi }_0\left(z\right)={\phi }_0\left(\zeta \right)=\sum _{j=0}^{\infty }{a}_j{\zeta }^j+\sum _{j=1}^{\infty }{b}_j{\zeta }^{-j},$$

(20)

$${\psi }_0\left(z\right)={\psi }_0\left(\zeta \right)=\sum _{j=0}^{\infty }{c}_j{\zeta }^j+\sum _{j=1}^{\infty }{d}_j{\zeta }^{-j},$$

(21)

where a_j, b_j, c_j, and d_j are unknowns.

The two analytical functions will be much more easily solved in ζ plane. Taking advantage of equations (19)–(21), the two stress boundary conditions are transformed as:

$$(2\mathrm{\Gamma }-P_1)\omega \left(\sigma \right)+\mathrm{\Gamma }\prime \overline{\omega \left(\sigma \right)}+{\phi }_0\left(\sigma \right)+\frac{\omega \left(\sigma \right)}{\overline{\omega \prime \left(\sigma \right)}}\overline{{\phi \prime }_0\left(\sigma \right)}+\overline{{\psi }_0\left(\sigma \right)}=0,$$

(22)

$$(2\mathrm{\Gamma }-P_2)\omega \left(R_0\sigma \right)+\mathrm{\Gamma }\prime \overline{\omega \left(R_0\sigma \right)}+{\phi }_0\left(R_0\sigma \right)+\frac{\omega \left(R_0\sigma \right)}{\overline{\omega \prime \left(R_0\sigma \right)}}\overline{{\phi \prime }_0\left(R_0\sigma \right)}+\overline{{\psi }_0\left(R_0\sigma \right)}=c_A,$$

(23)

where σ represents the coordinate on the outer boundary, and R₀σ represents the coordinate on the inner boundary in Figure 4.

If the stress boundary conditions $(\mathrm{\Gamma },\mathrm{\Gamma }\prime ,P_1\mathrm{and}{P}_2)$ are given and the mapping function is found, then equations (22) and (23) are better written as:

$${\phi }_0\left(\sigma \right)+\frac{\omega \left(\sigma \right)}{\overline{\omega \prime \left(\sigma \right)}}\overline{{\phi \prime }_0\left(\sigma \right)}+\overline{{\psi }_0\left(\sigma \right)}=-(2\mathrm{\Gamma }-P_1)\omega \left(\sigma \right)-\mathrm{\Gamma }\prime \overline{\omega \left(\sigma \right)},$$

(24)

$${\phi }_0\left(R_0\sigma \right)+\frac{\omega \left(R_0\sigma \right)}{\overline{\omega \prime \left(R_0\sigma \right)}}\overline{{\phi \prime }_0\left(R_0\sigma \right)}+\overline{{\psi }_0\left(R_0\sigma \right)}-c_A=-(2\mathrm{\Gamma }-P_2)\omega \left(R_0\sigma \right)-\mathrm{\Gamma }\prime \overline{\omega \left(R_0\sigma \right)}.$$

(25)

The fractions in equations (24) and (25) should be eliminated to obtain the equation set which could directly solve the coefficients (a_j, b_j, c_j, and d_j). It is obvious that $\omega (\sigma )$ , $\omega (R_0\sigma )$ , $\omega (\sigma )/\overline{\omega \prime (\sigma )}$ , and $\omega (R_0\sigma )/(\overline{\omega \prime (R_0\sigma )}$ contain fractions, and the fractions are too complicated to be eliminated by simple multiplication, but they can be transformed into Fourier series for $\omega (\sigma )$ is a period function with a period of 2π. The coefficients of the Fourier series can be numerically figured out through solving the linear equation set [13]. Then, $\omega (\sigma )$ , $\omega (R_0\sigma )$ , $\omega (\sigma )/\overline{\omega \prime \left(\sigma \right)}$ , and $\omega (R_0\sigma )/\overline{\omega \prime (R_0\sigma )}$ can be written as:

$$\omega (\sigma )=\sum _{j=1}^n{e}_j{\sigma }^{-j}+\sum _{j=0}^n{f}_j{\sigma }^j,$$

(26)

$$\omega (R_0\sigma )=\sum _{j=1}^n{g}_j{\sigma }^{-j}+\sum _{j=0}^n{h}_j{\sigma }^j,$$

(27)

$$\frac{\omega (\sigma )}{\overline{\omega \prime (\sigma )}}=\sum _{j=1}^n{l}_j{\sigma }^{-j}+\sum _{j=0}^n{m}_j{\sigma }^j,$$

(28)

$$\frac{\omega (R_0\sigma )}{\overline{\omega \prime (R_0\sigma )}}=\sum _{j=1}^n{s}_j{\sigma }^{-j}+\sum _{j=0}^n{t}_j{\sigma }^j.$$

(29)

What should be mentioned is that the coefficients ( $e_j$ , $f_j$ , $g_j$ , $h_j$ , $l_j$ , $m_j$ , $s_j$ , and $t_j$ ) in our earlier research [13] are real numbers while they are complex numbers in this paper, because the coefficients of the mapping function $\omega (\zeta )$ are complex numbers. Generally, the more the number(n) of the power items is, the smaller the error between the original function and the series is. The number(n) of the negative and positive power items is set as 500 in this paper.

As to equations (24) and (25), the following equations can be easily deduced:

$$\overline{{\psi }_0\left(\sigma \right)}=\sum _{j=0}^{\infty }{\overline{c}}_j{\sigma }^{-j}+\sum _{j=1}^{\infty }{\overline{d}}_j{\sigma }^j,$$

(30)

$$\overline{{\phi \prime }_0\left(\sigma \right)}=\sum _{j=1}^{\infty }j{\overline{a}}_j{\sigma }^{-j+1}-\sum _{j=1}^{\infty }j{\overline{b}}_j{\sigma }^{j+1},$$

(31)

$$\overline{{\psi }_0\left(R_0\sigma \right)}=\sum _{k=0}^{\infty }{\overline{c}}_k{R}_0^k{\sigma }^{-k}+\sum _{k=1}^{\infty }{\overline{d}}_k{R}_0^{-k}{\sigma }^k,$$

(32)

$$\overline{{\phi \prime }_0\left(R_0\sigma \right)}=\sum _{k=1}^{\infty }k{\overline{a}}_k{R}_0^{k-1}{\sigma }^{-k+1}-\sum _{k=1}^{\infty }k{\overline{b}}_k{R}_0^{-k-1}{\sigma }^{k+1}.$$

(33)

Substituting equations (26)–(33) into equations (24) and (25), equations (24) and (25) can be written as:

$$\begin{array}{c}\sum _{j=0}^{\infty }{a}_j{\sigma }^j+\sum _{j=1}^{\infty }{b}_j{\sigma }^{-j}+(\sum _{j=1}^n{l}_j{\sigma }^{-j}+\sum _{j=0}^n{m}_j{\sigma }^j)(\sum _{j=1}^{\infty }j{\overline{a}}_j{\sigma }^{-j+1}-\sum _{j=1}^{\infty }j{\overline{b}}_j{\sigma }^{j+1})+\sum _{j=0}^{\infty }{\overline{c}}_j{\sigma }^{-j}+\sum _{j=1}^{\infty }{\overline{d}}_j{\sigma }^j\hfill \\ =-\sum _{j=1}^n[(2\mathrm{\Gamma }-P_1)e_j+\mathrm{\Gamma }\prime f_j]{\sigma }^{-j}-(2\mathrm{\Gamma }-P_1+\mathrm{\Gamma }\prime )f_0-\sum _{j=1}^n[(2\mathrm{\Gamma }-P_1)f_j+\mathrm{\Gamma }\prime e_j]{\sigma }^j,\hfill \end{array}$$

(34)

$$\begin{array}{c}\sum _{k=0}^{\infty }{a}_k{R}_0^k{\sigma }^k+\sum _{k=1}^{\infty }{b}_k{R}_0^{-k}{\sigma }^{-k}+(\sum _1^n{s}_k{\sigma }^{-k}+\sum _0^n{t}_k{\sigma }^k)\times (\sum _{k=1}^{\infty }k{\overline{a}}_k{R}_0^{k-1}{\sigma }^{-k+1}-\sum _{k=1}^{\infty }k{\overline{b}}_k{R}_0^{-k-1}{\sigma }^{k+1})\hfill \\ +\sum _{k=0}^{\infty }{\overline{c}}_k{R}_0^k{\sigma }^{-k}+\sum _{k=1}^{\infty }{\overline{d}}_k{R}_0^{-k}{\sigma }^k-c_A=-\sum _{j=1}^n[(2\mathrm{\Gamma }-P_2)g_j+\mathrm{\Gamma }\prime h_j]{\sigma }^{-j}-(2\mathrm{\Gamma }-P_2+\mathrm{\Gamma }\prime )h_0-\hfill \\ \sum _{j=1}^n[(2\mathrm{\Gamma }-P_2)h_j+\mathrm{\Gamma }\prime g_j]{\sigma }^j.\hfill \end{array}$$

(35)

The coefficients (a_j, b_j, c_j, and d_j) can be obtained from equations (34) and (35) for the equalities of the two equations should hold for every σ. Therefore, the coefficient of every power of σ in the left sides of the two equations should equal the coefficient of the same power on the right sides. Since the right sides of equations (34) and (35) are known, an equation set related to the coefficients (a_j, b_j, c_j, and d_j) can be established.

The left sides of equations (34) and (35) contain infinite power items for ${\phi }_0(\zeta )$ and ${\psi }_0(\zeta )$ are the infinite series. They can be truncated at finite items to solve the coefficients because ${\phi }_0(\zeta )$ and ${\psi }_0(\zeta )$ are the convergent series. Then, we can obtain the solution of arbitrary degree of accuracy by truncating at different items. In this paper, a_j, b_j, c_j, and d_j are truncated at N, which means that a_j, b_j, c_j, and d_j are equal to zero while j > N.

Based on the above assumption, the unknowns that should be calculated out are a₀ + c_0 −c_A, a₀ + c₀, a₁∼a_N, b₁∼b_N, c₁∼c_N, and d₁∼d_N. Since a₀ and c₀ do not contribute to the stress, constant terms a₀ + c_0 −c_A and a₀ + c₀ do not have to be calculated out. Then, the equation set related to a₁∼a_N, b₁∼b_N, c₁∼c_N, and d₁∼d_N can be established by equating the coefficients of ${\sigma }^1$ ∼ ${\sigma }^N$ and ${\sigma }^{-1}$ ∼ ${\sigma }^{-N}$ in the two sides of equations (34) and (35), obtaining 4N equations [13]. Equations (34) and (35) contain the unknowns and their conjugates, so the formed 4N equations become extremely complicated when they multiplied by $e_j$ , $f_j$ , $g_j$ , $h_j$ , $l_j$ , $m_j$ , $s_j$ , and $t_j$ , resulting in that the equations cannot be solved directly. Thus, each equation should be divided into real part and image part, obtaining 8N equations, and merge them to solve the unknowns simultaneously:

$$\mathrm{UX}=V,$$

(36)

where

$$\begin{array}{l}X=[\mathrm{Re}\left(a_1\right),\mathrm{Re}\left(a_2\right),\dots ,\mathrm{Re}\left(a_i\right),\mathrm{Re}\left(b_1\right),\mathrm{Re}\left(b_2\right),\dots ,\mathrm{Re}\left(b_i\right),\mathrm{Re}\left(c_1\right),\mathrm{Re}\left(c_2\right),\dots ,\mathrm{Re}\left(c_i\right),\\ \mathrm{Re}\left(d_1\right),\mathrm{Re}\left(d_2\right),\dots ,\mathrm{Re}\left(d_i\right)\\ \mathrm{Im}\left(a_1\right),\mathrm{Im}\left(a_2\right),\dots ,\mathrm{Im}\left(a_i\right),\mathrm{Im}\left(b_1\right),\mathrm{Im}\left(b_2\right),\dots ,\mathrm{Im}\left(b_i\right),\mathrm{Im}\left(c_1\right),\mathrm{Im}\left(c_2\right),\dots ,\\ \mathrm{Im}\left(c_i\right),\mathrm{Im}\left(d_1\right),\mathrm{Im}\left(d_2\right),\dots ,\mathrm{Im}\left(d_i\right)].\end{array}$$

(37)

Once the two analytical functions are solved, the stress components in the image plane, described by orthogonal curvilinear coordinates, can be obtained from equations (1) and (2):

$${\sigma }_{\rho }+{\sigma }_{\theta }=4\mathrm{Re}\left[\frac{\phi \prime (\zeta )}{\omega \prime (\zeta )}\right].$$

(38)

$${\sigma }_{\theta }-{\sigma }_{\rho }+2i{\tau }_{\rho \theta }=\frac{2{\zeta }^2}{{\rho }^2\overline{\omega \prime \left(\zeta \right)}}[\overline{\omega \left(\zeta \right)}\frac{\phi ″\left(\zeta \right)\omega \prime \left(\zeta \right)-\omega ″\left(\zeta \right)\phi \prime v}{{[\omega \prime \left(\zeta \right)]}^2}+\psi \prime \left(\zeta \right)].$$

(39)

5. Comparison with ANSYS solution and stress result around the holes

Highly complicated equation deductions and solving large sparse linear equations are necessary to derive the analytical stress solution. Comparison with ANSYS solution is made to validate the correctness of deductions and solving process. Despite stress solution obtained from ANSYS being approximate, the stress distributing regularity is roughly right. So, the comparison is meaningful.

The example used for comparison is shown in Figure 5(a), in which the height of the straight side wall and the radius of the semicircle are both 2 m, the radius of the cavity is 1.5 m, and the horizontal stress P and vertical stress Q are set as 0.5 and 1.0 MPa, respectively. The ANSYS plane model is a big square with its side measuring 80 m to simulate the infinite domain. Using Plane2 element to mesh the area obtains 503,603 nodes and 251,122 elements in the numerical simulation process. The tangential stresses of the points of every degree on the two boundaries calculated by two methods are shown in Figure 5(b). The numerical solution and the analytical solution are in good agreement except the two corner points. The two biggest differences appear at 231° and 309°, and the relative error is about 10%. It is mainly because the tangential stress varies too fast at the corner while the ANSYS constant-stress element has an area.

Figure 5.

Comparison between numerical solution and analytical solution. (a) ANSYS model and (b) contrast of tangential stress on the boundaries.

To demonstrate the stress result around the holes of this case, stress cloud maps of ${\sigma }_{\theta }$ , ${\sigma }_{\rho }$ , and ${\tau }_{\rho \theta }$ in the orthogonal curvilinear coordinates are presented in Figure 6.

Figure 6.

Stress cloud maps of (a) ${\sigma }_{\theta }$ , (b) ${\sigma }_{\rho }$ , and (c) ${\tau }_{\rho \theta }$ .

The largest stress concentrations appear around the tunnel corners where the tangential stress is much larger than the two other stress components. And it can be found from Figure 6 that the area of great stress concentration of left corner is bigger than that of right corner.

6. Examples

Using the method proposed above, the stresses in the plane are calculated out for some situations. But only the tangential stress is showed in figures because of its importance.

6.1. Effect of lateral stress coefficient

The geometry is chosen the same as which is in section 5. The imposed vertical stress Q is set as 1.0 and three different horizontal stresses P: 0.5, 1.0, and 1.5 are considered to investigate the influence of axial stress on tangential stress. The obtained stress distribution on the boundaries of the two holes is demonstrated in Figure 7 by every 1°.

Figure 7.

Tangential stress for different lateral stress coefficients: (a) left hole and (b) right hole.

Figure 7(a) shows the tangential stress on the boundary of the tunnel. The stress values at the two corner points, namely, 231° and 309°, are much larger than other locations. And they increase as the horizontal stress increases, for instance, the stresses at right corner are 8.86, 10.69, and 12.54 when horizontal stresses are 0.5, 1.0, and 1.5, respectively. The tangential stress at right corner is bigger than the left one for the existence of the circular hole. The tangential stresses at the vault (about ${\alpha }_1\in [45°,135°]$ ) and floor (about ${\alpha }_1\in [231°,309°]$ ) of the tunnel increase with the increase in horizontal stress. Tensile stress will appear at the tunnel floor when horizontal stress is small, which is unfavorable for rock material. The tangential stresses at the side wall of the tunnel decrease as the horizontal stress increases.

Figure 7(b) shows the tangential stress on the boundary of the circular cavity. The stress distribution of the circular hole is nearly symmetric with respect to ${\alpha }_2=180°$ . The maximum tangential stress is located near ${\alpha }_2=180°$ where the circle boundary is close to the tunnel and it decreases little when horizontal stress increases. The low stress zone is near the upper and lower position of the circle where tangential stress is close to zero. And the low stress zone moves toward ${\alpha }_2=180°$ when horizontal stress increases.

There are four points on the either boundary at which tangential stress is zero under only horizontal stress. This phenomenon exists in all of the two hole problems [12,13].

6.2. Effect of distance of the two holes

Three distances of the two boundaries are investigated here: S = 1, S = 3, and S = 9 when the other shape parameters are kept the same as section 5. Generally, the effect of distance differs with respect to applied load, so two uniaxial loads—horizontal stress and vertical stress—are discussed here.

Figure 8(a) illustrates the tangential stress distribution on the tunnel boundary for different distances while the plate is under horizontal stress. The stress remains largely unchanged when the cavity moves away from tunnel. Tangential stresses are tensile at range ${\alpha }_1\in [0,39°)\cup (136°,231°)\cup (309°,360°)$ , namely, on the side wall of the tunnel. Figure 8(b) displays the tangential stress distribution on cavity boundary for different distances while the plate is under horizontal stress. The absolute value of tension stress and compression stress both decreases because of the existence of the tunnel. The maximum stress occurs at the top and bottom of the circle. With the increase in horizontal stress, stress at ${\alpha }_2\in (55°,155°)\cup (215°,315°)$ increases and stress at about ${\alpha }_2\in (155°,215°)$ decreases.

Figure 8.

Tangential stress for different distances of the two holes. (a) Left hole, P = 1.0, (b) right hole P = 1.0, (c) left hole, Q = 1.0, and (d) right hole Q = 1.0.

The change of tangential stress on the tunnel due to different distances while under vertical stress is presented in Figure 8(c). Similarly, it varies little except the stress fluctuation at ${\alpha }_1\in (125°,231°)$ becoming larger when the two holes are very close. In this load condition, tangential stresses at ${\alpha }_1\in (67°,112°)\cup (231°,309°)$ are tensile, namely, on the vault and floor. Figure 8(d) demonstrates the tangential stress on cavity boundary of different distances while under vertical stress. It is obvious that the stress at the boundary near the tunnel decreases as the distance increases, and stresses change small at other place of the hole.

6.3. Effect of the size of circular cavity

Sometimes in karst area, tunnel and circular karst cave are arranged in parallel and the karst cave will enlarge over time. It is interest to investigate the changement of stress due to the size of circle. In this part, the tunnel geometry and distance of the two holes are kept the same as section 5 while three circle radii—0.5, 1.5, and 3.0—are considered.

Figure 9(a) demonstrates the tangential stress around tunnel boundary for different circle radius while under horizontal stress. With the enlargement of the circular hole, the stresses at the right side wall increase, the stresses at right side of arch and stresses at the floor decrease. In particular, the tangential stresses at left corner point reduce greatly, which are 6.27, 5.52, and 3.31 here. The stress distribution on the boundary of round cave is shown in Figure 9(a). With the increasing of the circle size, tangential stresses at about the range ${\alpha }_1\in (35°,135°)\cup (225°,330°)$ , namely, the upper part and lower part of the circle, increase and tangential stresses at the right side boundary decrease.

Figure 9.

Tangential stress for different sizes of the cavity. (a) Left hole, P = 1.0, (b) right hole P = 1.0, (c) left hole, Q = 1.0, and(d) right hole Q = 1.0.

Figure 9(c) and (d) presents the tangential stress distribution on the tunnel boundary and circular hole boundary while under vertical stress, respectively. In this load condition, tangential stress on both boundaries changes little while the circle size increases except tangential stress at tunnel corner increasing remarkably, for instance, the stresses at left corner are 6.45, 7.61, and 9.63 for radius are 0.5, 1.0, and 1.5, respectively.

6.4. Effect of the position of circular cavity

As to practical tunnel engineering, the cavity may be located in any position in the plane. To study the effect caused by the position of circle, four cases are calculated (see Figure 10), which the angles of the connecting line (O₁O₂) of the two hole centers with respect to x-axis are −90°, −30°, 30°, and 90°, and the distance of two hole centers and the sizes of tunnel and cavity are the same as section 5.

Figure 10.

The tunnel close to a circular cavity of different positions.

Figure 11(a) and (c) demonstrates the tangential stress on the tunnel boundary when the plane is applied horizontal stress and vertical stress, respectively. The stress almost remains unchanged except little change appearing at the right side wall and floor of the tunnel.

Figure 11.

Tangential stress for different positions of circular cavity. (a) Left hole, P = 1.0, (b) right hole P = 1.0, (c) left hole, Q = 1.0, and (d) right hole Q = 1.0.

Figure 11(b) and (d) displays the tangential stress on the cavity boundary when the plane is applied horizontal stress and vertical stress, respectively. To compare the stress variation of the same point, the angles are measured from the same position in the circle shown in Figure 10. The stress distributions of cavity differ largely from its position. While under horizontal stress, the maximum tangential stress occurs at the top or bottom of the circle and the maximum tensile stress is located at ${\alpha }_2=180°$ . And with the increase in orientation angle, the maximum tangential stress first decreases then increases. While under vertical stress, the stress distributions differ a lot for orientation angle. The maximum tangential stresses are 1.82, 3.89, 3.87, and 2.06 for α = −90°, α = −30°, α = 30°, and α = 90°, respectively. The absolute values of tensile stress are also larger when α = −30° and α = 30°, which are 1.23 and 1.32, respectively.

7. Conclusion

This paper is based on a new discovery that the mapping function of an arbitrarily shaped hole can also manage a tunnel and a circular hole mapping into two circles. With the technique of least square method embedded in mixed penalty function method, the coefficients of the mapping function can be found through optimization. Combined with Korn’s function, the two holes can be finally mapped into a ring. Then, stress boundary conditions are established by two analytical functions which are solved by the power series method. And ANSYS simulation is conducted to compare with the analytical solution, which they are in good agreement. Some numerical examples are calculated out to study the effect by parameters, and tangential stress distributions on the boundaries are plotted in the figures. The method is applicable to problem of an arbitrarily shaped hole near to a neighboring circular cavity. However, since the mapping functions should be analytic, the hole boundary must be smooth curves and the solution near the sharp corner is approximate which is acceptable for tunnel engineering.

Undoubtedly, the corners of the tunnel are the places where the largest stress concentration is located in. And the stress at the two corner increases as the horizontal stress increases. The side wall of tunnel will produce tensile stress when the plane is subjected to horizontal stress and the vault and floor of the tunnel appear tensile stress when the plane is under vertical stress. The circle size has great effect on the boundary stress. With the enlargement of the circular hole, the maximum stress at the corners decreases while under horizontal stress and it is just the opposite while under vertical stress. The distance of the two holes slightly affects the stress on the tunnel while the stress on the left part of the circular hole is larger for small distance under vertical stress. The position of circle also influences the stress on the tunnel little, but the stress distribution on the circle varies considerably for its different position.