A Complex Variable Solution for Lining Stress and Deformation in a Non-Circular Deep Tunnel II Practical Application and Verification

A new complex variable method is presented for stress and displacement problems in a non-circular deep tunnel with certain given boundary conditions at infinity. In order to overcome the complex problems caused by non-circular geometric configurations and the multiply-connected region, a complex variable method and continuity boundary conditions are used to determine stress and displacement within the tunnel lining and within the surrounding rock. The coefficients in the conformal mapping function and stress functions are determined by the optimal design and complex variable method, respectively. The new method is validated by FLAC3D finite difference software through an example. Both the new method and the numerical simulation obtained similar results for the stress concentration and the minimum radial displacement occurred at a similar place in the tunnel. It is demonstrated that the new complex variable method is reliable and reasonable. The new method also provides another way to solve non-circular tunnel excavation problems in a faster and more accurate way.


Introduction
With the rapid economic development in China, which has caused the expansion of road and railway networks from east to west and to areas in the northeast that are surrounded by mountains, the construction of tunnels is broadly used to improve existing transportation networks. Lining is the primary support adopted to ensure rock pressure. It has been of high interest in determining stress fields within lining using analytical methods. Analytical solutions for stress distribution within circular lining and around circular and elliptical holes have been proposed by many authors (Bobet [1,2]; Lee and Nam [3]; Timoshenko and Goodier [4]).
Many studies have been carried out to determine the stress and deformation of tunnels by applying numerical methods, which have been generally used to provide an understanding of how the stress and deformation of lining are influenced by different parameters. Numerical methods are largely employed to find the stress and deformation stages of lining in the preliminary stages of design.
Muskhelishvili's [5] complex variable method is one of the useful analytical approaches that fully expounded a basic theory of complex potential functions in order to address some issues of plane elasticity mechanics. Based on this method, Exadaktyol and Stavtopoulou [6,7] proposed a closed-form plane strain solution for stress and displacement around semicircular holes. Verruijt [8,9] calculated the stress and displacement components around a circular tunnel in an elastic half-plane. Zhao and Yang [10] obtained a general solution for deep square tunnels with different pressure coefficients.
Kargar et al. [11] made an attempt to study lining stress and deformation in a non-circular deep tunnel using the Cauchy integral formula of complex variable methods. Li and Chen [12] obtained analytical solutions for a non-circular tunnel lining in power series forms. However, the analytical solutions for stress and displacement found in the above-mentioned literature have seldom come to an available expression for non-circular deep tunnels, especially non-circular deep tunnels with a lining, except for Kargar et al. [11] and Li [12].
When considering the problem with circular support, the solution is much easier. When a non-circular lining is included, the problem considers a multiply-connected region and conformal mapping, which increases the complexity of the problem. Kargar's [11] method can calculate the stress of a non-circular deep tunnel with a lining, although the method needs to integrate. The integral is too complex and difficult to calculate. Li's [12] method only considered the non-circular tunnel lining without thinking about the surrounding rock around the lining. In order to simplify the calculation, Li [12] assumed that the surrounding rock has a certain given surface traction applied to the lining. In this study it will be shown that these difficulties can be surmounted at last for the case of a non-circular deep tunnel with certain boundary conditions at infinity, by a new complex variable method in power series forms.
Based on the research findings proposed, an attempt was made in this paper to find the stress distribution and radial displacement of the lining in a non-circular deep tunnel, considering the boundary conditions of the surrounding rocks by applying a complex variable power series method, which is more efficient, simple and accurate. Finally, the new method was validated by the FLAC finite difference software through an example.

Problem Statement and General Consideration
The problem refers to a non-circular tunnel with lining in an elastic geomaterial. The tunnel is located at great depth compared with the tunnel dimension; the problem is considered a single hole with support in an infinite plane. The infinite plate on the complex plane is divided into the two isotropic homogenous regions of S 1 and S 2 bounded by contours L 1 and L 2 . The regions S 1 and S 2 refer to the lining and the surrounding rock, respectively. The boundary of the tunnel lining inside (L 1 ) is free of stress, and the rock-lining interface (L 2 ) satisfies the continuity boundary conditions. It is assumed that the region S 1 in the z-plane can be mapped conformally onto a ring (O 1 ) in the ζ-plane. The surrounding rock, region S 2 in the z-plane, can be mapped conformally onto the region O 2 , which is the area outside the L 2 outline in the ζ-plane, see Figure 1. The general formula of the conformal mapping function is determined based on the Laurent series as follows: where R is positive real number reflecting the hole's size, and the J k are generally complex numbers, which are determined by the shape of the tunnel. In most situations it is accurate enough to only take the first few of J k of the series. θ and ρ are assumed to be two polar coordinates of point ζ in the ζ-plane.  [11] made an attempt to study lining stress and deformation in a noncircular deep tunnel using the Cauchy integral formula of complex variable methods. Li and Chen [12] obtained analytical solutions for a non-circular tunnel lining in power series forms. However, the analytical solutions for stress and displacement found in the above-mentioned literature have seldom come to an available expression for non-circular deep tunnels, especially non-circular deep tunnels with a lining, except for Kargar et al. [11] and Li [12].
When considering the problem with circular support, the solution is much easier. When a noncircular lining is included, the problem considers a multiply-connected region and conformal mapping, which increases the complexity of the problem. Kargar's [11] method can calculate the stress of a non-circular deep tunnel with a lining, although the method needs to integrate. The integral is too complex and difficult to calculate. Li's [12] method only considered the non-circular tunnel lining without thinking about the surrounding rock around the lining. In order to simplify the calculation, Li [12] assumed that the surrounding rock has a certain given surface traction applied to the lining. In this study it will be shown that these difficulties can be surmounted at last for the case of a non-circular deep tunnel with certain boundary conditions at infinity, by a new complex variable method in power series forms.
Based on the research findings proposed, an attempt was made in this paper to find the stress distribution and radial displacement of the lining in a non-circular deep tunnel, considering the boundary conditions of the surrounding rocks by applying a complex variable power series method, which is more efficient, simple and accurate. Finally, the new method was validated by the FLAC finite difference software through an example.

Problem Statement and General Consideration
The problem refers to a non-circular tunnel with lining in an elastic geomaterial. The tunnel is located at great depth compared with the tunnel dimension; the problem is considered a single hole with support in an infinite plane. The infinite plate on the complex plane is divided into the two isotropic homogenous regions of S1 and S2 bounded by contours L1 and L2. The regions S1 and S2 refer to the lining and the surrounding rock, respectively. The boundary of the tunnel lining inside (L1) is free of stress, and the rock-lining interface (L2) satisfies the continuity boundary conditions. It is assumed that the region S1 in the z-plane can be mapped conformally onto a ring (O1) in the ζ-plane. The surrounding rock, region S2 in the z-plane, can be mapped conformally onto the region O2, which is the area outside the L2 outline in the ζ-plane, see Figure 1. The general formula of the conformal mapping function is determined based on the Laurent series as follows: where R is positive real number reflecting the hole's size, and the J k are generally complex numbers, which are determined by the shape of the tunnel. In most situations it is accurate enough to only take the first few of J k of the series. θ and ρ are assumed to be two polar coordinates of point ζ in the ζ-plane.  In the complex variable method [5,13,14], the solution is expressed in terms of two functions ϕ and ψ, which must be analyzed in the region S 1 occupied by the elastic material. The stresses are determined based on the stress functions in the equations: where σ ρ , σ θ and τ ρθ are the radial, circumferential and tangential stress components, respectively. The displacements are given by where G is the shear modulus of the elastic material, and κ is determined by Poisson's ratio υ, κ = 3 − 4υ and κ = 3−υ 1+υ are plane strain and plane stress, respectively. In this paper, plane strain conditions are assumed. Based on the Kargar's [11] method and Chen [15], the stress functions ϕ 1 and ψ 1 in the region O 2 have the following expansions: where where h k , m k are generally complex numbers, which must be determined from boundary conditions. ϕ 0 (ζ) and ψ 0 (ζ) are holomorphic functions with ϕ 0 (∞) = 0 and ψ 0 (∞) = 0, Γ and Γ are real and complex constants with regard to the stress state at infinity, which are determined as follows: where σ 1 and σ 2 are the principal stress components at infinity; α is the angle made between the σ 1 direction and the x axis; K is the lateral pressure coefficient; γ and H are the unit weight of the surrounding rock and the depth of tunnel, respectively. The stress functions ϕ 2 (ζ) and ψ 2 (ζ), which are represent region O 1 , have the following Laurent series expansions: where a k , b k , c k , d k are generally complex numbers that must be determined from boundary conditions. Stress functions ϕ 1 (ζ), ψ 1 (ζ), ϕ 2 (ζ) and ψ 2 (ζ) should satisfy the continuity boundary conditions on R 2 circles and satisfy the boundary conditions on R 1 circles. The displacement at lining and rock should be equal on the boundary R 2 . It can be concluded as: where u R2 ρ1 and u R2 θ1 are the lining displacement components of boundary R 2 in the ρ and θ directions, and u R2 ρ2 and u R2 θ2 are the rock displacement components of boundary R 2 in the ρ and θ directions, respectively.
The stress at lining and rock should be equal on boundary R 2 . It can be concluded as: where f R2 1 and f R2 2 are displacement components of boundary R 2 in the lining and rock, respectively. The stress at the lining on the boundary R 1 should be 0. It can be concluded as: Based on Chen [15], the boundary conditions (13)-(15) can be rewritten. Equation (13) is expressed as: The stress boundary condition of Equation (14) is expressed as: The stress boundary condition of Equation (15) is expressed as: where t 1 and t 2 are the point of the boundary R 1 and R 2 in the ζ-plane, respectively. u ρ1 and u θ1 are the displacement components of the lining-rock interface in the ρ and θ directions. The expressions Γw(ζ) and Γ w(ζ) should not be incorporated into continuity Equation (16) since they define initial ground stress and deformation in the surrounding rock when tunnels are excavated. Equations (16) and (17) are concerned with the continuity of deformation and the stress field across the lining-rock interface due to the no-slip condition. Equation (18) is concerned that the tunnel lining inside (L 1 ) is entirely free of stress.

Solution
In order to eliminate the difficulties caused by the power series, Equations (16)-(18) are rewritten in the form: The form of an infinite polynomial times another infinite polynomial, such as ϕ 0 (t 2 )w (t 2 ) of Equation (19), is hard to calculated and merge so that the equation cannot be calculated. To acquire a solution for Equations (19)-(21), Li's method needs to be introduced, which converts the form of the multiplication of two infinite polynomials to an infinite matrix, which can be easily calculated and merged for a computer.
From the power series theory of the complex variable method written by Chen [15], the first factor of Equation (19) on the left can be rewritten as follows: where k is substituted by v to facilitate the calculation. σ is the angle of the point of the boundary; ρ is the radius of the boundary circle in the ζ-plane, which equals 1 at the lining-rock interface. σ and ρ can be related by ζ = ρσ, where σ = exp (iθ). L v can be calculated by R and J k , and Equation (1) is rewritten in the form: The first item on the right of the Equation (22) must be expanded as follows: where the positive power of σ in Equation (24) is obtained as: The zero and negative powers of σ in Equation (24) can be derived: The second item on the right of Equation (22) can be expanded as follows: There are only zero and negative powers of σ in the above Equation (23), which can be derived as follows: σ 0 The general formula of Equation (29) are determined based on the expanded functions of Equation (25) as follows: The general formula of Equation (30) is determined based on the expanded functions of Equations (26) and (28) as follows: Based on the power series of the complex variable method written by Chen [15], which have been presented in Equations (22)-(30), the other items of Equation (19) can be determined and separate the positive and negative exponents of Equation (19); the positive power system of Equation (19) can be expanded as follows: The negative power system of Equation (19) can be expanded as follows: The positive power system of Equation (20) is determined as follows: The negative power system of Equation (20) can be expanded as follows: where ρ is the radius of the circle, which equals 1 at the lining-rock interface. The positive power system of Equation (21) is determined as follows: The negative power system of Equation (20) can be expanded as follows: Math. Comput. Appl. 2018, 23, 43 where ρ is the radius of the circle, which equals 0.8633 at the lining-atmosphere interface. Equations (31)-(36) can be written as a set of linear equations:  (37) is linearly related when k → ∞ and ρ = 1. The number of conditions in the equation set is reduced from (6v + 3) to (6v + 2). Considering that ϕ 0 (∞) = 0 and ψ 0 (∞) = 0, the coefficients h k and m k can take any value if k → ∞ . It is indicated that h ∞ and m ∞ can take zero. Since the coefficients h ∞ and m ∞ are known, the number of unknown coefficients in the equation set is reduced from (6k + 4) to (6k + 2). The number of unknown coefficients equals the boundary conditions, in which all the coefficients have been determined uniquely. Equations (5), (6), (11) and (12) are obtained by the calculated coefficients, and the stress and deformation of the non-circular deep tunnel are obtained through Equations (2)-(4).

Implementation
In this section, the new complex variable method is applied to an example and a comparison is provided with FLAC finite difference software in order to verify the formula.

Fundamental Assumption
(a) The tunnel is assumed to have an infinite length; the surrounding rock mass is homogeneous, isotropic and linear elastic and without creep or viscous behaviors. (b) The tunnel's length and depth are assumed to far outweigh its diameter; the surrounding rock mass conforms to the plane strain condition (κ = 3 − 4υ).

Comparison of the New Analytical Solution with That of the Numerical Simulation Results
The tunnel distribution diagram is presented and a 3916 zones and 13887 grid-points finite mesh calculation model was used to simulate stress and displacement distribution in Figure 2. The horizontal displacement of the finite mesh calculation model is constrained by the left and right boundary, the vertical displacement is constrained by the bottom boundary, and the top boundary is free and unconstrained. The numerical model is concerned with continuity of deformation and the stress field across the lining-rock interface due to the no-slip condition. The tunnel lining inside should be entirely free of stress.
The calculation parameters are shown in Table 1.  According to Lv's [16,17] method and the geometry of the tunnel (in Figure 2), the conformal mapping function (Equation (1)) was determined and provided by a self-programming optimal design software as follows: where coefficient k = 4 is close enough to Equation (1). The radius ρ = 0.8633 is related to ζ by ζ = ρσ, where σ = exp (iθ). It is assumed that the tunnel lining inside (L1) can be mapped conformally onto a circle (R1). When the radius ρ = 1 the L2 can be mapped conformally onto R2.
As an example, the boundary condition across the lining-rock interface and tunnel lining inside can be determined by Equations (13)- (15). The boundary conditions at infinity can be expressed through Equations (9) and (10).
Based on Equation (37), a simple computer program written by MATLAB was applied to solve the problem.
As the coordinates of the analytical solution and numerical simulation are different, it is necessary to rewrite the results. The comparison of the rewritten results between the analytical solution and numerical simulation are shown below. Figure 3 shows the circumferential stress along the rock-lining interface predicted by the new analytical solution and the FLAC finite difference software. It can be observed that the maximum circumferential stress happens at a position of an 85-degree angle (i.e., the widest part of the tunnel) and, moreover, the circumferential stress of the new analytical solution has not declined rapidly, but creates a stress concentration at = 115°. The circumferential stress along the inner lining periphery is presented in Figure 4a, demonstrating good agreement between the new analytical solution and the numerical simulation apart from an 85-degree angle. Normal stress and shear stress along the inner lining periphery is illustrated in Figure 4b   According to Lv's [16,17] method and the geometry of the tunnel (in Figure 2), the conformal mapping function (Equation (1)) was determined and provided by a self-programming optimal design software as follows: where coefficient k = 4 is close enough to Equation (1). The radius ρ = 0.8633 is related to ζ by ζ = ρσ, where σ = exp (iθ). It is assumed that the tunnel lining inside (L 1 ) can be mapped conformally onto a circle (R 1 ). When the radius ρ = 1 the L 2 can be mapped conformally onto R 2 .
As an example, the boundary condition across the lining-rock interface and tunnel lining inside can be determined by Equations (13)- (15). The boundary conditions at infinity can be expressed through Equations (9) and (10).
Based on Equation (37), a simple computer program written by MATLAB was applied to solve the problem.
As the coordinates of the analytical solution and numerical simulation are different, it is necessary to rewrite the results. The comparison of the rewritten results between the analytical solution and numerical simulation are shown below. Figure 3 shows the circumferential stress along the rock-lining interface predicted by the new analytical solution and the FLAC finite difference software. It can be observed that the maximum circumferential stress happens at a position of an 85-degree angle (i.e., the widest part of the tunnel) and, moreover, the circumferential stress of the new analytical solution has not declined rapidly, but creates a stress concentration at α = 115 • . The circumferential stress along the inner lining periphery is presented in Figure 4a, demonstrating good agreement between the new analytical solution and the numerical simulation apart from an 85-degree angle. Normal stress and shear stress along the inner lining periphery is illustrated in Figure 4b Figure 5 shows the radial displacement of the tunnel along the rock-lining interface predicted by the analytical solution. It could be demonstrated that the radial displacement along the rock-lining interface was in good agreement with the displacement boundary condition, which proves its high accuracy.  Figure 5 shows the radial displacement of the tunnel along the rock-lining interface predicted by the analytical solution. It could be demonstrated that the radial displacement along the rock-lining interface was in good agreement with the displacement boundary condition, which proves its high accuracy.  Figure 5 shows the radial displacement of the tunnel along the rock-lining interface predicted by the analytical solution. It could be demonstrated that the radial displacement along the rock-lining interface was in good agreement with the displacement boundary condition, which proves its high accuracy. The radial displacement of the tunnel along the inner lining periphery is illustrated in Figure 6, which shows that the radial displacement of the analytical solution and numerical simulation are all zero at α = 75°. It was demonstrated that the numerical solution agrees well with the analytical The radial displacement of the tunnel along the inner lining periphery is illustrated in Figure 6, which shows that the radial displacement of the analytical solution and numerical simulation are all zero at α = 75 • . It was demonstrated that the numerical solution agrees well with the analytical solution. Considering the results of Figure 5, the analytical solution had good agreement with the displacement boundary condition, thus the analytical solution results were reasonable.  The discrepancies between the new analytical solution and the numerical simulation described above, especially Figure 4, may be due to the fact that the grid size in the numerical modeling is not small enough to produce accurate results. The stress information of the numerical simulation is stored in zones not in grid points, and the stress cannot be accurately expressed along a specified boundary.

Conclusions
In this paper, the stress and displacement of a non-circular deep tunnel and within their lining supports were studied using a new analytical solution, which is based on the basic theory of complex variables and plane elasticity [17], and the following conclusions can be made.
The analytic functions were exactly established to predict the stress and displacement distribution of the non-circular deep tunnel within their lining supports, but it is obviously not entirely true that the stress and displacement value is only determined by the in-situ stress boundary conditions and coefficient of the elasticity modulus, Poisson's ratio, lateral pressure and material density. The analytical solution for radial displacement was smaller than the numerical simulation results, however, and further study will be needed to develop these functions.
Due to the fact that the grid size in modeling was not small enough to produce accurate results and stress information was stored in zones, not in grid points in the numerical simulation, the stress cannot be accurately expressed along a specified boundary. But the analytical solution results were not affected by grid size and zones in the numerical modeling. The discrepancies between the new analytical solution and the numerical simulation described above, especially Figure 4, may be due to the fact that the grid size in the numerical modeling is not small enough to produce accurate results. The stress information of the numerical simulation is stored in zones not in grid points, and the stress cannot be accurately expressed along a specified boundary.

Conclusions
In this paper, the stress and displacement of a non-circular deep tunnel and within their lining supports were studied using a new analytical solution, which is based on the basic theory of complex variables and plane elasticity [17], and the following conclusions can be made.
The analytic functions were exactly established to predict the stress and displacement distribution of the non-circular deep tunnel within their lining supports, but it is obviously not entirely true that the stress and displacement value is only determined by the in-situ stress boundary conditions and coefficient of the elasticity modulus, Poisson's ratio, lateral pressure and material density.
The analytical solution for radial displacement was smaller than the numerical simulation results, however, and further study will be needed to develop these functions.
Due to the fact that the grid size in modeling was not small enough to produce accurate results and stress information was stored in zones, not in grid points in the numerical simulation, the stress cannot be accurately expressed along a specified boundary. But the analytical solution results were not affected by grid size and zones in the numerical modeling.
The curves of the stress value showed that the new analytical solution and numerical simulation were in reasonable agreement. Both solutions predicted a normal stress concentration at the lower and upper corners of the tunnel, and both maximum circumferential stress results occurred in the widest part of the tunnels. The normal and shear stress values of the tunnel along the inner lining periphery were almost zero, which proved its high accuracy.
Although numerical simulation is the main tool for solving tunnel excavation problems, especially non-circular tunnels, the complex variable method can provide another way to solve non-circular tunnel excavation problems in a faster and more accurate way.