# Elastic Analysis for Subaqueous Tunnel Surrounding Rock via the Complex Variable Method

Authors: Li-yuan Yu; Hong-wen Jing; Ying-chao Wang
Journal: Mathematical Problems in Engineering
 (2014)

---

## Abstract

Generally speaking, the subaqueous tunnels can be regarded as the shallow-buried ones. Consequently, the classical problem of an elastic half plane with a round cavity, loaded arbitrarily along the surface boundary, can be used to obtain the stress and displacement fields of the surrounding rock for this type of tunnels. The solution uses the complex variable method, with a conformal mapping onto a circular ring in the image plane. Because of the convergence of the complex potentials throughout the annular region, the coefficients in the Laurent series expansion form for complex functions can be determined by a system of liner recurrent equations, obtained from both the horizontal and the cavity boundary conditions.  The stresses and deformations of the surrounding rock can then be calculated via some relevant equations. The whole calculation program should be coded by Fortran language. As an example, the case of a specific underwater tunnel is considered in some detail eventually.


---

## Body

## 1. Introduction

Since the modern times, water bodies, such as rivers, lakes, and bays, have become the considerable restricting factor for the development of economy and society with the progress of land transportation. Traditionally, the common modes of transportation to overcome these obstacles were ferries and bridges. But, in recent decades, the subaqueous tunnel has progressed vigorously due to its special features.The rock covering thickness, a critical parameter to subaqueous tunnels, is no more than twice the tunnel span in general case [1]. Therefore, the subaqueous tunnel can be regarded as shallow-buried tunnel with ground load derived from the overlying water. Besides, effective reinforcement measures can substantially improve complex geological surrounding rock in mechanical properties such as continuity, homogeneity, and isotropy. For the abovementioned reasons, the solution of mechanical field for underwater tunnels can be simplified as a classical problem of an elastic half plane with a circular cavity, loaded arbitrarily on the horizontal boundary.Compared with deep-buried tunnels, the shallow-buried ones are much more difficult to address via mathematic methods owing to the surface boundary and ground load. The influence of the surface was taken into account by Jeffery [2] and Mindlin [3, 4] with the aid of the bipolar coordinates method. Nevertheless, this approach is incapable of providing the deformation field in the situation that the tunnel is quite adjacent to the ground surface. Bobet [5] proposed another elastic solution for ground deformations of a shallow tunnel in a saturated ground by using the general series form stress function in polar coordinate. Bobet’s solution was obtained with the boundary condition of uniform radial displacement at the tunnel opening, while Park [6] presented the elastic solution to predict the tunnelling-induced undrained ground deformation by imposing prescribed oval-shaped displacement at the tunnel periphery. The virtual image technique was firstly employed by Sagaseta [7] to consider the presence of top free surface for a shallow tunnel in an incompressible soil layer. Verruijt and Booker [8] extended Sagaseta’s method to acquire a simple closed-form analytical solution for a tunnel in a homogeneous elastic half space.The complex variable method, a once-popular means for problems of multiply connected regions [9], gained the capacity to solve the more general case of problems for regions bounded by two eccentric circles due to the creative attempt from Verruijt [10] and Strack and Verruijt [11]. The advantage of complex variable method in solving elastic problems is that both the stresses and displacements can be obtained simultaneously and that the approaches to handle stress or displacement boundary conditions are analogous. The method has recently yielded solutions, in terms of infinite series, for the problems of specified displacements [10] or stresses [12] along the tunnel boundary.Although many research results about shallow-buried tunnels have been presented as mentioned above, the problem derived from the underwater tunnel, that is, an elastic half plane with a circular cavity, loaded uniformly in vertical direction on the upper boundary, has never been analyzed. This math conundrum will be solved in this paper by the complex variable method, with a conformal mapping onto a circular ring.


## 2. Description of the Problem

The problem deals with an elastic half plane with a round tunnel; see Figure1. The boundary of the tunnel is free of stress and displacement, and loading takes place along the ground boundary of half plane, in the form of a given distribution of surface force, denoted by q
(
z
). The elastic region R represents rock or soil medium in which the tunnel is located. The radius of the hole is denoted by r
0, the depth of its centre below the surface by h, and the covering thickness by d; see Figure 1. The ratio r
/
h will be considered as the essential geometrical parameter.Figure 1
Elastic half plane with a round tunnel.The complex variable method for the solution of two-dimensional linear elastic problems [9] involves two analytic functions of a complex variable, z
=
x
+
i
y, where x is the real part of z and y is the imaginary part of z. The two complex potentials marked as φ
1
(
z
), ψ
1
(
z
) are bound to be analytic all over the elastic region R. The stresses are related to these complex functions by the following equations [9]:
(1)
σ
y
+
σ
x
=
4
Re
ϕ
1
′
(
z
)
,
σ
y
-
σ
x
+
2
i
τ
x
y
=
2
[
z
-
ϕ
1
′
′
(
z
)
+
ψ
1
′
(
z
)
]
,
where σ
x and σ
y are the normal stresses in the x and y directions, respectively, and τ
x
y represents the shear stress. The horizontal and vertical displacements, u
x and u
y, in this elastic region R can be expressed as follows (for the plane strain problem):
(2)
2
G
(
u
x
+
i
u
y
)
=
(
3
-
4
μ
)
ϕ
1
(
z
)
-
z
ϕ
1
′
(
z
)
¯
-
ψ
1
(
z
)
¯
,
where G denotes the shear modulus of the elastic material and μ represents Poisson’s ratio.Despite the definite surface force along both boundaries, it is usually more convenient to express the boundary condition in terms of the integral of the surface tractions, integrated along the boundary:(3)
[
ϕ
1
(
z
)
+
z
ϕ
1
′
(
z
)
¯
+
ψ
1
(
z
)
¯
]
s
=
i
∫
(
F
~
x
+
i
F
~
y
)
d
s
.
Thus the boundary conditions of this specific case are formulated as follows:
(4)
y
=
0
:
φ
1
(
z
)
+
z
ϕ
1
′
(
z
)
¯
+
ψ
1
(
z
)
¯
=
i
∫
l
1
z
(
q
x
+
i
q
y
)
d
z
=
Q
1
(
z
)
+
C
,
(5)
x
2
+
(
y
+
h
)
2
=
r
2
:
φ
1
(
z
)
+
z
φ
1
′
(
z
)
¯
+
ψ
1
(
z
)
¯
=
0
,
where Q
1
(
z
) is a given function of complex variable z along the upper boundary and C is an unknown integration constant. This constant can be assumed to be zero along the cavity boundary without loss of generality. The precise form of the continuous function Q
1
(
z
) depends on the actual stress distribution along the surface boundary. The integrand q
x
+
i
q
y represents the component-wise manner of the prescribed surface force distribution function q
(
z
), and the lower limit of integration l
1 is the initial point of the prescribed surface force distribution.


## 3. Conformal Mapping to an Annulus

We conformally map the regionR in the z-plane onto an annular region γ in the transformed ζ-plane (see Figure 2).Figure 2
The conformally mappedζ-plane.The ringγ is bounded by the circles |
ζ
|
=
1 and |
ζ
|
=
α, where α
<
1. The appropriate conformal transformation is given by [10]
(6)
z
=
ω
(
ζ
)
=
-
i
h
1
-
α
2
1
+
α
2
1
+
ζ
1
-
ζ
,
where
(7)
r
0
h
=
2
α
1
+
α
2
with
α
=
h
-
h
2
-
r
0
2
r
0
.
The depth h and radius r
0 for the tunnel are shown in Figure 1. It can be easily confirmed that the circles |
ζ
|
=
1 and |
ζ
|
=
α in Figure 2 correspond to the upper (
y
=
0
) and the cavity (
x
2
+
(
y
+
h
)
2
=
r
0
2
) boundaries in Figure 1, respectively.The complex potentialsφ
1
(
z
) and ψ
1
(
z
), which must be analytic throughout the region R in the z-plane, can be considered as functions of ζ thanks to the analytic transformation function ω
(
ζ
) in the annular region γ. Consequently, the two functions are both analytic in the region γ in the ζ-plane, and then they can be expanded into Laurent series in the same region including its boundaries:
(8)
ϕ
1
(
z
)
=
ϕ
1
(
ω
(
ζ
)
)
=
ϕ
(
ζ
)
=
∑
k
=
0
∞
a
k
ζ
k
+
∑
k
=
1
∞
b
k
ζ
-
k
,
(9)
ψ
1
(
z
)
=
ψ
1
(
ω
(
ζ
)
)
=
ψ
(
ζ
)
=
∑
k
=
0
∞
c
k
ζ
k
+
∑
k
=
1
∞
d
k
ζ
-
k
.
The coefficients a
k, b
k, c
k, and d
k must be determined by the boundary conditions. Both series expansions are supposed to converge throughout the annular region γ up to its boundaries. The point ζ
=
-
i corresponding to infinity in the z-plane merits some careful consideration. By virtue of the stress distributions of an elastic half plane with a normal point load acting on the upper boundary (see (10), derived from Flamant [9]), it can be inferred that the stresses tend towards zero at infinity and that the displacements will be finite at infinity; therefore the Laurent series will converge on the point ζ
=
-
i:
(10)
σ
x
=
-
2
P
π
x
3
(
x
2
+
y
2
)
2
,
σ
y
=
-
2
P
π
x
y
2
(
x
2
+
y
2
)
2
,
τ
x
y
=
-
2
P
π
x
2
y
(
x
2
+
y
2
)
2
.


## 4. Boundary Conditions

### 4.1. Transformation of the Boundary Equations

The boundary conditions are formulated in terms of the functionsφ
1
(
z
), ψ
1
(
z
), and z
φ
1
′
(
z
)
¯; see (4) and (5). When transforming the two conditions in terms of the variable ζ, the term z
φ
1
′
(
z
)
¯ needs particular concerns. Because the term φ
1
′
(
z
) can be defined as d
φ
1
/
d
z, the corresponding term φ
1
′
(
ζ
) can be determined as follows via the chain rule:
(11)
φ
1
′
(
ζ
)
=
d
ϕ
1
d
ζ
=
d
ϕ
1
d
z
d
z
d
ζ
=
ϕ
1
′
(
z
)
ω
′
(
ζ
)
.
It now follows that
(12)
z
φ
1
′
(
z
)
¯
=
ω
(
ζ
)
ϕ
1
′
(
ζ
)
¯
ω
′
(
ζ
)
¯
.
The mathematical difficulties involved in solving boundary value problems for a certain region depend on the factor ω
(
ζ
)
/
ω
′
(
ζ
)
¯. We have ζ
=
ρ
σ on a circle with radius ρ in the ζ-plane, where σ
=
exp
⁡
⁡
(
i
ϑ
). Then, ζ
¯
=
ρ
σ
-
1. This gives
(13)
ω
(
ζ
)
ω
′
(
ζ
)
¯
=
-
1
2
(
1
+
ρ
σ
)
(
σ
-
ρ
)
2
σ
2
(
1
-
ρ
σ
)
.
This factor appears to be relatively simple for the circular tunnel case. However, for problems with a complex shape hollow the factor may be so complicated that it actually prohibits analytic solution of the problem.


### 4.2. The Surface Boundary Condition

There is a given distribution of surface force, denoted asq
(
z
), along the ground boundary (see Figure 1). In virtue of (4), the transformed form of the boundary condition along the corresponding boundary in the ζ-plane can be obtained as follows:
(14)
|
ζ
|
=
1
:
φ
(
ζ
)
+
ω
(
ζ
)
ω
′
(
ζ
)
¯
φ
′
(
ζ
)
¯
+
ψ
(
ζ
)
¯
=
Q
(
ζ
)
+
C
.
The radius ρ
=
1 on this boundary; then the expression (13) can be simplified substantially:
(15)
ω
(
ζ
)
ω
′
(
ζ
)
¯
=
1
2
(
1
-
σ
-
2
)
.
The function Q
(
ζ
), in general, can be transformed into a Fourier series, provided that the distribution function of surface load, q
(
z
)
=
q
x
+
i
q
y, has been prescribed. If the range of the surface load is intended to be [
l
1
,
l
2
] (see Figure 1), the right side of (4) can be integrated as follows:
(16)
y
=
0
:
Q
1
(
x
)
=
i
∫
l
1
x
(
q
x
+
i
q
y
)
d
x
=
{
0
,
x
∈
(
-
∞
,
l
1
)
,
(
-
q
y
+
i
q
x
)
(
x
-
l
1
)
,
x
∈
[
l
1
,
l
2
]
,
(
-
q
y
+
i
q
x
)
(
l
2
-
l
1
)
,
x
∈
(
l
2
,
+
∞
)
.
Along the outer boundary |
ζ
|, the radius ρ
=
1, so that ζ
=
ρ
σ
=
σ
=
exp
⁡
⁡
(
i
ϑ
). Then the expression (16) can be converted as follows:
(17)
|
ζ
|
=
1
:
Q
(
ζ
)
=
{
(
-
q
y
+
i
q
x
)
(
σ
2
-
σ
1
)
,
θ
∈
[
0
,
θ
2
)
,
(
-
q
y
+
i
q
x
)
(
σ
-
σ
1
)
,
θ
∈
[
θ
2
,
θ
1
]
,
0
,
θ
∈
(
θ
1
,
2
π
)
.
Apparently, as a periodic function, the fluctuation cycle of Q
(
ζ
) is 2
π. In addition, this function satisfies the Dirichlet’s condition of convergence. Consequently, this function can be written as a Fourier series,
(18)
Q
(
ζ
)
=
∑
k
=
-
∞
+
∞
e
k
σ
k
,
where
(19)
k
=
0
,
e
0
=
-
q
y
+
i
q
x
2
π
[
σ
1
(
1
-
1
2
σ
1
)
-
σ
2
(
1
-
1
2
σ
2
)
]
k
>
0
,
e
k
=
-
q
y
+
i
q
x
2
k
π
0000000000
×
{
[
+
cos
⁡
(
k
σ
1
)
-
cos
⁡
(
k
σ
2
)
k
(
σ
1
-
σ
2
)
sin
k
00000000000000
+
cos
⁡
⁡
(
k
σ
1
)
-
cos
⁡
⁡
(
k
σ
2
)
k
]
000000000000
+
i
[
sin
(
k
σ
1
)
-
sin
(
k
σ
2
)
k
(
σ
1
-
σ
2
)
cos
⁡
k
0000000000000000
-
sin
⁡
(
k
σ
1
)
-
sin
⁡
(
k
σ
2
)
k
]
}
k
>
0
,
e
-
k
=
e
k
¯
.
On the basis of the Laurent series expansions (8) and (9) and the expression (15), the left-hand side of (14) can be elaborated out. The expression (18) completes the expansion of the right-hand side of (14). The resulting equation must be satisfied for all possible values of σ, so (20) can be obtained finally, leading to sums of positive, zero, and negative powers of σ:
(20)
positive
powers

:
c
k
=
e
k
-
b
k
¯
+
1
2
(
k
-
1
)
a
k
-
1
positive
powers
：
00
-
1
2
(
k
+
1
)
a
k
+
1
,
k
=
1,2
,
…
,
zero
power

:
c
o
=
e
o
¯
-
a
o
¯
-
1
2
a
1
-
1
2
b
1
+
C
¯
,
negative
powers

:
d
k
=
e
k
¯
-
a
k
¯
+
1
2
(
k
-
1
)
b
k
-
1
negative
powers
：
00
-
1
2
(
k
+
1
)
b
k
+
1
,
k
=
1,2
,
…
.
The unknown coefficients c
k and d
k can be expressed with the other ones a
k and b
k via (20). Once a
k and b
k have been figured out, the determination of c
k and d
k is explicit and straightforward. Fortunately, the coefficients a
k and b
k can be obtained from the tunnel boundary condition.


### 4.3. The Cavity Boundary Condition

By means of (5), the free boundary condition for the cavity can be transformed as follows in the ζ-plane:
(21)
|
ζ
|
=
α
:
φ
(
ζ
)
+
ω
(
ζ
)
ω
′
(
ζ
)
¯
φ
′
(
ζ
)
¯
+
ψ
(
ζ
)
¯
=
0
.
On this boundary, ζ
=
α
σ, so the expression (13) can be transformed as follows:
(22)
ω
(
ζ
)
ω
′
(
ζ
)
¯
=
-
α
σ
-
(
1
-
2
α
2
)
+
α
(
2
-
α
2
)
σ
-
1
-
α
2
σ
-
2
2
(
1
-
α
σ
)
.
A considerably complicated system of equations for the coefficients a
k, b
k, c
k, and d
k can be derived after substituting (8), (9), and (22) into (21). Furthermore, the coefficients c
k and d
k can be eliminated by using (20). The system of equations turns out to be less complicated, with only two levels of coefficients. The final result is given by
(23)
positive
powers

:
(
1
-
α
2
)
α
2
k
(
k
+
1
)
a
k
+
1
positive
powers
：
+
(
1
-
α
2
k
+
2
)
b
k
+
1
¯
positive
powers
：
=
(
1
-
α
2
)
α
2
k
k
a
k
+
(
1
-
α
2
k
)
b
k
¯
positive
powers
：
+
α
2
k
e
k
-
α
2
k
+
2
e
k
+
1
,
k
=
1,2
,
…
,
negative
powers

:
(
1
-
α
2
k
+
2
)
a
k
+
1
+
(
1
-
α
2
)
(
k
+
1
)
b
k
+
1
¯
negative
powers
：
=
α
2
(
1
-
α
2
k
)
a
k
+
(
1
-
α
2
)
k
b
k
¯
negative
powers
：
+
e
k
+
1
-
α
2
e
k
,
k
=
1,2
,
…
.
From these two equations the coefficients can be calculated using a recursion method. If the values a
1 and b
1 are determined, the coefficients a
2 and b
2 can be figured out and so on. The initial solution, a
1 and b
1, can be obtained from the conditions that the coefficients of σ
0 and σ
-
1 must be zero. This gives
(24)
σ
0
:
(
α
2
-
1
)
(
a
1
+
b
1
¯
)
+
C
¯
=
α
2
e
1
-
e
0
,
σ
-
1
:
(
α
2
-
1
)
(
a
1
+
b
1
¯
)
-
α
2
C
=
α
2
e
0
-
e
1
-
α
2
a
0
.
The constant a
0 does not produce stresses in the medium and corresponds to only an arbitrary rigid body displacement. Consequently, it can be regarded as zero. Then the following two solutions are readily available:
(25)
C
=
Re
(
e
1
-
e
0
)
+
i
α
2
+
1
α
2
-
1
Im
⁡
(
e
1
-
e
0
)
,
(26)
a
1
+
b
1
¯
=
α
2
e
0
-
e
1
+
α
2
C
α
2
-
1
.
The coefficients a
1 and b
1 cannot be determined only by (26). Fortunately, this difficulty can be overcome via the convergence of the complex potentials φ
(
ζ
) and ψ
(
ζ
) for all values of ζ in the ring γ. The convergence of the series expansions (8) and (9) at the point B
′
(
1,0
) in the ζ-plane (as described in Section 3) requires that all coefficients tend to zero for k
→
∞, and this is not automatically ensured. Because the iterative equations (24) are linear and because the corresponding homogeneous equations require a
k
=
-
b
k
¯, we can infer that an arbitrary constant can be added up to each of these coefficients without altering the solution of the equations given above. Accordingly, the coefficient b
1 can be determined by (26), assuming a
1
=
0 first. It can be expected that, for a very large value of k (k
=
1000 or 10,000), a constant limiting value, except zero, will be obtained for the coefficient a
k. Subtracting this limiting value from a
k and -
b
k
¯, the precise values for these coefficients can be found easily. The residual coefficients c
k and d
k can be figured out from (20) eventually. Meanwhile, the complex potentials φ
(
ζ
) and ψ
(
ζ
) are completely determined.


## 4.1. Transformation of the Boundary Equations

The boundary conditions are formulated in terms of the functionsφ
1
(
z
), ψ
1
(
z
), and z
φ
1
′
(
z
)
¯; see (4) and (5). When transforming the two conditions in terms of the variable ζ, the term z
φ
1
′
(
z
)
¯ needs particular concerns. Because the term φ
1
′
(
z
) can be defined as d
φ
1
/
d
z, the corresponding term φ
1
′
(
ζ
) can be determined as follows via the chain rule:
(11)
φ
1
′
(
ζ
)
=
d
ϕ
1
d
ζ
=
d
ϕ
1
d
z
d
z
d
ζ
=
ϕ
1
′
(
z
)
ω
′
(
ζ
)
.
It now follows that
(12)
z
φ
1
′
(
z
)
¯
=
ω
(
ζ
)
ϕ
1
′
(
ζ
)
¯
ω
′
(
ζ
)
¯
.
The mathematical difficulties involved in solving boundary value problems for a certain region depend on the factor ω
(
ζ
)
/
ω
′
(
ζ
)
¯. We have ζ
=
ρ
σ on a circle with radius ρ in the ζ-plane, where σ
=
exp
⁡
⁡
(
i
ϑ
). Then, ζ
¯
=
ρ
σ
-
1. This gives
(13)
ω
(
ζ
)
ω
′
(
ζ
)
¯
=
-
1
2
(
1
+
ρ
σ
)
(
σ
-
ρ
)
2
σ
2
(
1
-
ρ
σ
)
.
This factor appears to be relatively simple for the circular tunnel case. However, for problems with a complex shape hollow the factor may be so complicated that it actually prohibits analytic solution of the problem.


## 4.2. The Surface Boundary Condition

There is a given distribution of surface force, denoted asq
(
z
), along the ground boundary (see Figure 1). In virtue of (4), the transformed form of the boundary condition along the corresponding boundary in the ζ-plane can be obtained as follows:
(14)
|
ζ
|
=
1
:
φ
(
ζ
)
+
ω
(
ζ
)
ω
′
(
ζ
)
¯
φ
′
(
ζ
)
¯
+
ψ
(
ζ
)
¯
=
Q
(
ζ
)
+
C
.
The radius ρ
=
1 on this boundary; then the expression (13) can be simplified substantially:
(15)
ω
(
ζ
)
ω
′
(
ζ
)
¯
=
1
2
(
1
-
σ
-
2
)
.
The function Q
(
ζ
), in general, can be transformed into a Fourier series, provided that the distribution function of surface load, q
(
z
)
=
q
x
+
i
q
y, has been prescribed. If the range of the surface load is intended to be [
l
1
,
l
2
] (see Figure 1), the right side of (4) can be integrated as follows:
(16)
y
=
0
:
Q
1
(
x
)
=
i
∫
l
1
x
(
q
x
+
i
q
y
)
d
x
=
{
0
,
x
∈
(
-
∞
,
l
1
)
,
(
-
q
y
+
i
q
x
)
(
x
-
l
1
)
,
x
∈
[
l
1
,
l
2
]
,
(
-
q
y
+
i
q
x
)
(
l
2
-
l
1
)
,
x
∈
(
l
2
,
+
∞
)
.
Along the outer boundary |
ζ
|, the radius ρ
=
1, so that ζ
=
ρ
σ
=
σ
=
exp
⁡
⁡
(
i
ϑ
). Then the expression (16) can be converted as follows:
(17)
|
ζ
|
=
1
:
Q
(
ζ
)
=
{
(
-
q
y
+
i
q
x
)
(
σ
2
-
σ
1
)
,
θ
∈
[
0
,
θ
2
)
,
(
-
q
y
+
i
q
x
)
(
σ
-
σ
1
)
,
θ
∈
[
θ
2
,
θ
1
]
,
0
,
θ
∈
(
θ
1
,
2
π
)
.
Apparently, as a periodic function, the fluctuation cycle of Q
(
ζ
) is 2
π. In addition, this function satisfies the Dirichlet’s condition of convergence. Consequently, this function can be written as a Fourier series,
(18)
Q
(
ζ
)
=
∑
k
=
-
∞
+
∞
e
k
σ
k
,
where
(19)
k
=
0
,
e
0
=
-
q
y
+
i
q
x
2
π
[
σ
1
(
1
-
1
2
σ
1
)
-
σ
2
(
1
-
1
2
σ
2
)
]
k
>
0
,
e
k
=
-
q
y
+
i
q
x
2
k
π
0000000000
×
{
[
+
cos
⁡
(
k
σ
1
)
-
cos
⁡
(
k
σ
2
)
k
(
σ
1
-
σ
2
)
sin
k
00000000000000
+
cos
⁡
⁡
(
k
σ
1
)
-
cos
⁡
⁡
(
k
σ
2
)
k
]
000000000000
+
i
[
sin
(
k
σ
1
)
-
sin
(
k
σ
2
)
k
(
σ
1
-
σ
2
)
cos
⁡
k
0000000000000000
-
sin
⁡
(
k
σ
1
)
-
sin
⁡
(
k
σ
2
)
k
]
}
k
>
0
,
e
-
k
=
e
k
¯
.
On the basis of the Laurent series expansions (8) and (9) and the expression (15), the left-hand side of (14) can be elaborated out. The expression (18) completes the expansion of the right-hand side of (14). The resulting equation must be satisfied for all possible values of σ, so (20) can be obtained finally, leading to sums of positive, zero, and negative powers of σ:
(20)
positive
powers

:
c
k
=
e
k
-
b
k
¯
+
1
2
(
k
-
1
)
a
k
-
1
positive
powers
：
00
-
1
2
(
k
+
1
)
a
k
+
1
,
k
=
1,2
,
…
,
zero
power

:
c
o
=
e
o
¯
-
a
o
¯
-
1
2
a
1
-
1
2
b
1
+
C
¯
,
negative
powers

:
d
k
=
e
k
¯
-
a
k
¯
+
1
2
(
k
-
1
)
b
k
-
1
negative
powers
：
00
-
1
2
(
k
+
1
)
b
k
+
1
,
k
=
1,2
,
…
.
The unknown coefficients c
k and d
k can be expressed with the other ones a
k and b
k via (20). Once a
k and b
k have been figured out, the determination of c
k and d
k is explicit and straightforward. Fortunately, the coefficients a
k and b
k can be obtained from the tunnel boundary condition.


## 4.3. The Cavity Boundary Condition

By means of (5), the free boundary condition for the cavity can be transformed as follows in the ζ-plane:
(21)
|
ζ
|
=
α
:
φ
(
ζ
)
+
ω
(
ζ
)
ω
′
(
ζ
)
¯
φ
′
(
ζ
)
¯
+
ψ
(
ζ
)
¯
=
0
.
On this boundary, ζ
=
α
σ, so the expression (13) can be transformed as follows:
(22)
ω
(
ζ
)
ω
′
(
ζ
)
¯
=
-
α
σ
-
(
1
-
2
α
2
)
+
α
(
2
-
α
2
)
σ
-
1
-
α
2
σ
-
2
2
(
1
-
α
σ
)
.
A considerably complicated system of equations for the coefficients a
k, b
k, c
k, and d
k can be derived after substituting (8), (9), and (22) into (21). Furthermore, the coefficients c
k and d
k can be eliminated by using (20). The system of equations turns out to be less complicated, with only two levels of coefficients. The final result is given by
(23)
positive
powers

:
(
1
-
α
2
)
α
2
k
(
k
+
1
)
a
k
+
1
positive
powers
：
+
(
1
-
α
2
k
+
2
)
b
k
+
1
¯
positive
powers
：
=
(
1
-
α
2
)
α
2
k
k
a
k
+
(
1
-
α
2
k
)
b
k
¯
positive
powers
：
+
α
2
k
e
k
-
α
2
k
+
2
e
k
+
1
,
k
=
1,2
,
…
,
negative
powers

:
(
1
-
α
2
k
+
2
)
a
k
+
1
+
(
1
-
α
2
)
(
k
+
1
)
b
k
+
1
¯
negative
powers
：
=
α
2
(
1
-
α
2
k
)
a
k
+
(
1
-
α
2
)
k
b
k
¯
negative
powers
：
+
e
k
+
1
-
α
2
e
k
,
k
=
1,2
,
…
.
From these two equations the coefficients can be calculated using a recursion method. If the values a
1 and b
1 are determined, the coefficients a
2 and b
2 can be figured out and so on. The initial solution, a
1 and b
1, can be obtained from the conditions that the coefficients of σ
0 and σ
-
1 must be zero. This gives
(24)
σ
0
:
(
α
2
-
1
)
(
a
1
+
b
1
¯
)
+
C
¯
=
α
2
e
1
-
e
0
,
σ
-
1
:
(
α
2
-
1
)
(
a
1
+
b
1
¯
)
-
α
2
C
=
α
2
e
0
-
e
1
-
α
2
a
0
.
The constant a
0 does not produce stresses in the medium and corresponds to only an arbitrary rigid body displacement. Consequently, it can be regarded as zero. Then the following two solutions are readily available:
(25)
C
=
Re
(
e
1
-
e
0
)
+
i
α
2
+
1
α
2
-
1
Im
⁡
(
e
1
-
e
0
)
,
(26)
a
1
+
b
1
¯
=
α
2
e
0
-
e
1
+
α
2
C
α
2
-
1
.
The coefficients a
1 and b
1 cannot be determined only by (26). Fortunately, this difficulty can be overcome via the convergence of the complex potentials φ
(
ζ
) and ψ
(
ζ
) for all values of ζ in the ring γ. The convergence of the series expansions (8) and (9) at the point B
′
(
1,0
) in the ζ-plane (as described in Section 3) requires that all coefficients tend to zero for k
→
∞, and this is not automatically ensured. Because the iterative equations (24) are linear and because the corresponding homogeneous equations require a
k
=
-
b
k
¯, we can infer that an arbitrary constant can be added up to each of these coefficients without altering the solution of the equations given above. Accordingly, the coefficient b
1 can be determined by (26), assuming a
1
=
0 first. It can be expected that, for a very large value of k (k
=
1000 or 10,000), a constant limiting value, except zero, will be obtained for the coefficient a
k. Subtracting this limiting value from a
k and -
b
k
¯, the precise values for these coefficients can be found easily. The residual coefficients c
k and d
k can be figured out from (20) eventually. Meanwhile, the complex potentials φ
(
ζ
) and ψ
(
ζ
) are completely determined.


## 5. Approximately Analytical Solution in Physical Plane

In order to simplify expressions, the parameterm is introduced as follows:
(27)
m
=
1
-
α
2
1
+
α
2
h
.
Then the conformal mapping relationship (see (6)) can be transformed into the following formula:
(28)
ζ
=
(
z
+
i
m
)
(
z
-
i
m
)
.
The corresponding point ζ in the image plane can be calculated for an arbitrary point z
=
x
+
i
y in the physical plane, using (28). The complex potentials φ
(
ζ
) and ψ
(
ζ
) are then obtained by (8) and (9). The stress and displacement components can finally be determined via the following equations:
(29)
σ
y
+
σ
x
=
4
Re
[
φ
′
(
ζ
)
ω
′
(
ζ
)
]
,
σ
y
-
σ
x
+
2
i
τ
x
y
=
2
ω
′
(
ζ
)
{
ω
(
ζ
)
¯
[
φ
′
(
ζ
)
ω
′
(
ζ
)
]
′
+
ψ
′
(
ζ
)
}
,
2
G
(
u
x
+
i
u
y
)
=
(
3
-
4
γ
)
ϕ
(
ζ
)
-
ω
(
ζ
)
ω
′
(
ζ
)
¯
ϕ
′
(
ζ
)
¯
-
ψ
(
ζ
)
¯
.
In order to avoid the laborious calculations, the whole aforementioned calculation process is strongly suggested to be coded by Fortran language.


## 6. Example

As an example, a subaqueous tunnel with a radius of 5 m is considered. The exemplary tunnel, with a buried depth of 20 m, suffers from hydrostatic pressure due to the overlying water of 30 m in height. Although the scope of hydrostatic pressure along the ground boundary is quite extensive, the range of tenfold tunnel diameter can meet accuracy requirements for the aforementioned engineering problem. The precise parameters for this example are shown in Table1.Table 1
Calculation preferences for a specific underwater tunnel.
r
/
m
h
/
m
E/GPa
γ
l
1
/
m
l
2
/
m
q
y/MPa
5
25
10
0.3
−55
55
−0.3


### 6.1. Stress Analysis

The principal stresses can be obtained from the normal stresses and shear stress as follows:(30)
σ
1
σ
2
=
σ
x
+
σ
y
2
±
(
σ
x
-
σ
y
2
)
2
+
τ
x
y
2
.
The contour plots for the principal stresses σ
1 and σ
2 are shown in Figure 3 for a specific domain (70 m and 80 m in the horizontal and vertical orientation, resp.). As shown in the figure, the stress concentration is comparatively obvious within double the tunnel span. The major principal stress σ
1 reaches the maximum (−8.14 KPa) at the tunnel vault and the minimum (−200 KPa) at the arch foot part, while, for the minor principal stress σ
2, the maximum (−80 KPa) appears at the tunnel vault and the minimum (−800 KPa) appears at the middle of side walls. The negative sign here indicates compressive stress conventionally.Figure 3
Contour plots for the principal stresses.


### 6.2. Displacement Analysis

The settlement curves at different depths in the overburden are shown in Figure4, while Figure 5 represents the cross-section profiles before and after the tunnel excavation (notice that displacements are magnified a thousand times for a legible display). The tunnel cross-section area reduces by about 8.3% due to the excavation. The settlement abnormally decreases with increasing depth owing to without considering gravity. As the distance between the tunnel central line and the specific point increases, the settlement curves drop more evidently with increasing burial depth, and, for any particular curve, the settlement magnitude sees a gradual reduction despite being increasingly slight.Figure 4
Settlement curves at different depths in the overburden.Figure 5
Cross-section profiles of the tunnel before and after excavation.


## 6.1. Stress Analysis

The principal stresses can be obtained from the normal stresses and shear stress as follows:(30)
σ
1
σ
2
=
σ
x
+
σ
y
2
±
(
σ
x
-
σ
y
2
)
2
+
τ
x
y
2
.
The contour plots for the principal stresses σ
1 and σ
2 are shown in Figure 3 for a specific domain (70 m and 80 m in the horizontal and vertical orientation, resp.). As shown in the figure, the stress concentration is comparatively obvious within double the tunnel span. The major principal stress σ
1 reaches the maximum (−8.14 KPa) at the tunnel vault and the minimum (−200 KPa) at the arch foot part, while, for the minor principal stress σ
2, the maximum (−80 KPa) appears at the tunnel vault and the minimum (−800 KPa) appears at the middle of side walls. The negative sign here indicates compressive stress conventionally.Figure 3
Contour plots for the principal stresses.


## 6.2. Displacement Analysis

The settlement curves at different depths in the overburden are shown in Figure4, while Figure 5 represents the cross-section profiles before and after the tunnel excavation (notice that displacements are magnified a thousand times for a legible display). The tunnel cross-section area reduces by about 8.3% due to the excavation. The settlement abnormally decreases with increasing depth owing to without considering gravity. As the distance between the tunnel central line and the specific point increases, the settlement curves drop more evidently with increasing burial depth, and, for any particular curve, the settlement magnitude sees a gradual reduction despite being increasingly slight.Figure 4
Settlement curves at different depths in the overburden.Figure 5
Cross-section profiles of the tunnel before and after excavation.


## 7. Conclusion

The subaqueous tunnel can be idealized as a classical problem of an elastic half plane with a round cavity, loaded uniformly along the ground boundary in the normal direction. It has been shown that the complex variable method can be applied to this problem successfully. The method used consists of first finding out the appropriate conformal mapping function and writing the Laurent series expansion for both of the analytic functions. All the coefficients of the Laurent series expansions are determined by both the surface and the cavity boundary conditions, thanks to the convergence of the series expansions throughout the annular region in the image plane. The stress and displacement components for the elastic medium in the physical plane can finally be calculated via some relevant equations. As an example, the solution method has been illustrated by calculating the deformations and stresses for the case of a specific underwater tunnel.



---

*Source: 492745-2014-09-25.xml*