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

Generally speaking, the subaqueous tunnels can be regarded as the shallow-buried ones. Consequently, the classical problem of an elastic half planewith 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.


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.

Description of the Problem
The problem deals with an elastic half plane with a round tunnel; see Figure 1.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 ().The elastic region  represents rock or soil medium in which the tunnel is located.The radius of the hole is denoted by  0 , the depth of its centre below the surface by ℎ, and the covering thickness by ; see Figure 1.The ratio /ℎ will be considered as the essential geometrical parameter.
The complex variable method for the solution of twodimensional linear elastic problems [9] involves two analytic functions of a complex variable,  =  + , where  is the real part of  and  is the imaginary part of .The two complex potentials marked as  1 (),  1 () are bound to be analytic all over the elastic region .The stresses are related to these complex functions by the following equations [9]: where   and   are the normal stresses in the  and  directions, respectively, and   represents the shear stress.
The horizontal and vertical displacements,   and   , in this elastic region  can be expressed as follows (for the plane strain problem): where  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: Thus the boundary conditions of this specific case are formulated as follows: where  1 () is a given function of complex variable  along the upper boundary and  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  1 () depends on the actual stress distribution along the surface boundary.The integrand   +   represents the component-wise manner of the prescribed surface force distribution function (), and the lower limit of integration  1 is the initial point of the prescribed surface force distribution.

Conformal Mapping to an Annulus
We conformally map the region  in the -plane onto an annular region  in the transformed -plane (see Figure 2).The ring  is bounded by the circles || = 1 and || = , where  < 1.The appropriate conformal transformation is given by [10] where The depth ℎ and radius  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 ( = 0) and the cavity ( 2 + ( + ℎ) 2 =  2 0 ) boundaries in Figure 1, respectively.The complex potentials  1 () and  1 (), which must be analytic throughout the region  in the -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: The coefficients   ,   ,   , and   must be determined by the boundary conditions.Both series expansions are supposed to converge throughout the annular region  up to its boundaries.The point  = − corresponding to infinity in the -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  = −: (10)  4) and ( 5).When transforming the two conditions in terms of the variable , the term   1 () needs particular concerns.Because the term   1 () can be defined as  1 /, the corresponding term   1 () can be determined as follows via the chain rule:

Boundary Conditions
It now follows that 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().Then,  =  −1 .This gives 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.

The Surface Boundary Condition.
There is a given distribution of surface force, denoted as (), 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: = 1 :  () +  () ()   () +  () =  () + .
The radius  = 1 on this boundary; then the expression (13) can be simplified substantially: The function (), in general, can be transformed into a Fourier series, provided that the distribution function of surface load, () =   +  , has been prescribed.If the range of the surface load is intended to be [ 1 ,  2 ] (see Figure 1), the right side of ( 4) can be integrated as follows: Along the outer boundary ||, the radius  = 1, so that  =  =  = exp().Then the expression (16) can be converted as follows: Apparently, as a periodic function, the fluctuation cycle of () is 2.In addition, this function satisfies the Dirichlet's condition of convergence.Consequently, this function can be written as a Fourier series, where 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 : The unknown coefficients   and   can be expressed with the other ones   and   via (20).Once   and   have been figured out, the determination of   and   is explicit and straightforward.Fortunately, the coefficients   and   can be obtained from the tunnel boundary condition.

The Cavity Boundary Condition
. By means of ( 5), the free boundary condition for the cavity can be transformed as follows in the -plane: =  :  () +  () On this boundary,  = , so the expression (13) can be transformed as follows: A considerably complicated system of equations for the coefficients   ,   ,   , and   can be derived after substituting ( 8), (9), and ( 22) into (21).Furthermore, the coefficients   and   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 positive powers : From these two equations the coefficients can be calculated using a recursion method.If the values  1 and  1 are determined, the coefficients  2 and  2 can be figured out and so on.The initial solution,  1 and  1 , can be obtained from the conditions that the coefficients of  0 and  −1 must be zero.This gives The constant  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: The coefficients  1 and  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   (1, 0) in the -plane (as described in Section 3) requires that all coefficients tend to zero for  → ∞, and this is not automatically ensured.
Because the iterative equations ( 24) are linear and because the corresponding homogeneous equations require   = −  , 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  1 can be determined by (26), assuming  1 = 0 first.It can be expected that, for a very large value of  ( = 1000 or 10, 000), a constant limiting value, except zero, will be obtained for the coefficient   .Subtracting this limiting value from   and −  , the precise values for these coefficients can be found easily.
The residual coefficients   and   can be figured out from (20) eventually.Meanwhile, the complex potentials () and () are completely determined.

Approximately Analytical Solution in Physical Plane
In order to simplify expressions, the parameter  is introduced as follows: Then the conformal mapping relationship (see (6)) can be transformed into the following formula: The corresponding point  in the image plane can be calculated for an arbitrary point  =  +  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: In order to avoid the laborious calculations, the whole aforementioned calculation process is strongly suggested to be coded by Fortran language.

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 Table 1.
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.

Displacement Analysis.
The settlement curves at different depths in the overburden are shown in Figure 4, 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 crosssection 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.

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.

Figure 3 :
Figure 3: Contour plots for the principal stresses.

Figure 4 :Figure 5 :
Figure 4: Settlement curves at different depths in the overburden.

Table 1 :
Calculation preferences for a specific underwater tunnel.