# A Complex Variable Solution for Rectangle Pipe Jacking in Elastic Half-Plane

Authors: Xin-yuan Li; Guo-bin Liu
Journal: Mathematical Problems in Engineering
 (2017)

---

## Abstract

In mechanics, the solution of soil stresses and displacements field caused by shallow rectangular jacking pipe construction can be simplified as half-plane problem. Both the boundary conditions of the surface and the cavity boundary must be taken into account. It is the essential prerequisite for mechanical analysis of the pipe jacking with the complex variable theory that the mechanical boundary must be transformed from the half-plane with a rectangle cavity to the concentric ring. According to Riemann’s existence theorem and basic complex variable theory, a conformal mapping function is established. Both sides of boundary conditions equation are developed into Laurent series, and then the coefficients of complex stress function are solved by power series method. The derived solution is applied to an example and a comparison is made using FEM method to show the accuracy of the methods. The result shows the following:(1) the method presented in this paper is applicable to a shallow-buried rectangular tunnel; (2) the complex function method proposed in this paper is characterized by clear steps, fast convergence, and simple operation.


---

## Body

## 1. Introduction

Underground tunnel excavation work increases with economic development. A number of methods have been proposed to calculate the stresses and displacements caused by tunnel excavation [1, 2]. For instance, the complex variables method developed by Muskhelishvili [3] is particularly suitable for solving underground tunnel problems as it has the tool of conformal mapping, which can transform the complex boundary shape into a simple shape. Using this technique, numerous highly accurate analytical results have been obtained by scholars [4, 5]. The problems solved by the complex function method are mainly deep-buried tunnels. The complex variables method can be applied for solving stresses and displacement for single or multiple tunnels of arbitrary shapes at great depth [6].The complex variables method was also proposed for stress and displacement around a circular tunnel in an elastic half-plane [7–12]. However, solving shallow-buried tunnel of rectangular shape using the complex function method is still difficult. The main difficulty is that the conformal mapping function from the external domain of the rectangular excavation section in the half-plane to the concentric ring domain is highly complex [7] and cannot be applied directly as in Verruijt’s work.An attempt was made in this study to determine a conformal mapping function and boundary conditions handling method that provides an approach for the solution of the difficulty.


## 2. Conformal Mapping Function

The focus in this study is the stresses and displacements for a shallow-buried rectangular shape tunnel. To solve the problem using the complex variable method, the conformal transformation must be recurred. To transform the half-planex<h with the exclusion of the hole (as shown in Figure 1) into the concentric ring domain (as shown in Figure 2), the transformation is given as(1)z=0.6h0.3+∑j=1mCjζj,where Cj are coefficients to be determined, m is the number of the coefficients which depend on the shape of the cavity and mapping accuracy, and h is the distance between the horizon line and coordinate, as shown in Figure 1. If the tunnel is symmetrical about the x-axis, the coefficients Cj in (1) must be real numbers [13].Figure 1
Shallow cavern with rectangular excavation cross-section.Figure 2
The region after conformal mapping.Function (1) can be regarded as (2)z=0.6h0.3+w,(3)w=∑j=1mCjζj.The horizon linex=h is mapped to the circumference line Γ′2 whose center locates at coordinates with a radius of 0.3 in the w plane with function (2), whereas the tunnel is mapped to the outer boundary Γ′1 whose equation can be determined uniquely. If the specific mapping function (3) is adopted, the outer boundary Γ′1 in the w plane is mapped to the unit circle Γ′′1 in the ζ plane [13]. The circumference Γ′2 is mapped to internal circle Γ′′2 with the same function.According to the boundary correspondence theory, the function that can map the external boundary of thew plane into the unit circle of the plane in the ζ plane is the desired mapping function [14]. Cj can be solved with an iterative method between even and odd points [15]. The value of r can be determined according to the corresponding condition of Γ′2 and Γ′′2. An m order algebraic equation about r can be obtained through formula (3) by assuming point (0.3,0) in the w plane mapped to point A′′ in the ζ plane accurately (4)0.3=∑j=1mrjCj.AsCj have been obtained, the value of r can be obtained by solving the real (<1) root of this equation [13].It can be verified easily that the circleζ=r corresponds to the horizon line x=h and that the circle ζ=1 corresponds to the tunnel.


## 3. Basic Equations for Solving Complex Stress Function

The complex variables method for the rectangular shape tunnel excavation problem in the elastic half-plane can be resolved to solve the complex stress function according to surface boundary conditions and tunnel boundary conditions.Because the conformal transformation functionz(ζ) is analytic in the ring D′′ bounded by the circles ζ=1 and ζ=r, the functions φ(z) and ψ(ζ), which must be analytic throughout D′′, can be expanded in the form of Laurent series [16](5)φz=φζ=a0+∑k=1∞akζk+∑k=1∞bkζ-k,ψz=ψζ=c0+∑k=1∞ckζk+∑k=1∞dkζ-k.If the values ofak, bk, ck, dk are known, the expressions of φ(z) and ψ(ζ) can be obtained. Thus, the problem can be resolved to solve for the values of ak, bk, ck, dk according to the boundary conditions.The first boundary condition is that the upper boundaryx=h must be entirely free of stress. This yields(6)φz+zφ′z¯+ψz¯x=h=0.When this condition is transformed to theζ plane, its form is(7)φζ+zζz′ζ¯φ′ζ¯+ψζ¯=0.Becausez(ζ) is an analytic function, the expression z(ζ)/z′ζ¯ must be an analytic function in ring D′′. This means that they can be represented by their Laurent series expansions [8].Letζ=ρeiα=ρδ, where ρ and α are the radius and argument of the points on the ζ plane, respectively. We have ρ=r on the surface boundary; thus (8)zζz′ζ¯=∑j=0mCjrjδ-j2∑j=0mCjrjδj∑j=1mCjrj-1δ1-j=∑-∞∞Ekδk,where Ek are undetermined real coefficients. For specific problems, Ek can be obtained by comparing the coefficient of same power of δ on the left and right sides of the equation(9)∑j=1-mk+sDjEk-j=0k∈1-m-s,m-s,∑j=1-mmDjEk-j=0k∈1+m-s,-2m-1,∑j=1-mmDjEk-j=∑j=-k-mmCjC-k-jr-kk∈-2m,-m,∑j=1-mmDjEk-j=∑j=0-kCjC-k-jr-kk∈1-m,0,∑j=1-mmDjEk-j=0k∈1,1+s-m,∑j=k-smDjEk-j=0k∈2-m+s,m+s,where s is the number of the coefficients which depend on the calculation accuracy. We take s=50 in this paper(10)Dk=∑j=0k+m-1CjCj-k-1j-k-1rj-k-2k∈1-m,0∑j=kmCjCj-k+1j-k+1rj-kk∈1,m.We have2s+1 undetermined coefficients, but there are 2s+2m equations. They are overdetermined equations. Ek can be obtained using the least square method.By substituting (5) and (8) into (6), the boundary condition on x=h can be expressed as (11)a0+∑k=1∞akrkδk+∑k=1∞bkr-kδ-k+∑-∞∞Ekδk∑k=1∞-kbkr-k-1δk+1+∑-∞∞Ekδk∑k=1∞kakrk-1δ1-k+∑k=0∞ckrkδ-k+∑k=1∞dkr-kσk=0.A group of equations to solveak, bk, ck, and dk can be obtained by setting the coefficients of all powers of δ equal to zero. The result is(12)ckrk=-bkr-k+∑v=1∞vE-v-k-1bvr-v-1-∑v=1∞vEv-k-1avrv-1,dkr-k=-akrk+∑v=1∞vE-v+k-1bvr-v-1-∑v=1∞vEv+k-1avrv-1,a0-∑v=1∞vE-v-1bvr-v-1+∑v=1∞vEv-1avrv-1+c0=0.The coefficientsck and dk can be obtained if ak and bk have been determined. From the boundary condition at the tunnel boundary, another group of solving equations can be obtained.The second boundary value problem, in which the displacement is prescribed along the tunnel boundary, is considered. The boundary condition at the corresponding boundary in theζ plane is(13)2Gu¯+iv¯=κφζ-zζz′ζ¯φ′ζ¯-ψζ¯,where u¯ and v¯ are given displacements in x and y directions, respectively.It can be assumed thatu¯+iv¯ can be written as a Fourier series(14)u¯+iv¯=∑v=-∞∞Fveivθ,where(15)Fv=12π∫02πu¯+iv¯e-ivθdθ.θ is the argument of the points on the tunnel boundary in the z plane. It can be expected that for all problems of practical significance, such an expansion is possible.The tunnel boundary in theζ plane can be expressed as(16)ζ=δ.Based on the analytic functionz(ζ), the one-to-one corresponding relations between θ and α are determined as (17)θ=fαand we have(18)2G∑v=-∞∞Fveivθ=2G∑v=-∞∞Fveivfα.Each item of formula (18) can be expanded into a Fourier series with a period of 2π. Thus,(19)2GFveivfα=∑k=-∞∞Gvkeikα=∑k=-∞∞Gvkδk,where(20)Gvk=12π∫02π2GFveivfαe-ikαdαand thus(21)2Gu¯+iv¯=∑v=-∞∞∑k=-∞∞Gvkeikα.For the tunnel boundaryζ=δ, we have(22)zζz′ζ¯=∑j=0mCjδ-j2∑j=0mCjδj∑j=1mCjδ1-j=∑-∞∞Hkδk,where Hk are undetermined real coefficients. For specific problems, Hk can be obtained by comparing the coefficient of the same power of δ on the left and right sides of the equation.By substituting (21) and (22) into (13), we have(23)∑v=-∞∞∑k=-∞∞Gvkδk=κa0+∑k=1∞akδk+∑k=1∞bkδ-k+∑k=1∞∑v=1∞vH-v+k-1bvδk+∑k=0∞∑v=1∞vH-v-k-1bvδ-k-∑k=1∞∑v=1∞vHv+k-1avδk-∑k=0∞∑v=1∞vHv-k-1avδ-k-∑k=0∞ckδ-k-∑k=1∞dkδk.Another group of equations to solveak, bk, ck, and dk can be obtained by setting the coefficients of all powers of δ equal to zero. The result is(24)∑v=-∞∞Gvk=κak+∑v=1∞vH-v+k-1bv-∑v=1∞vHv+k-1av-dk,∑v=-∞∞Gv-k=κbk+∑v=1∞vH-v-k-1bv-∑v=1∞vHv-k-1av-ck,∑v=-∞∞Gv0=κa0+∑v=1∞vH-v-1bv-∑v=1∞vHv-1av-c0.a
k, bk, ck, and dk can be obtained by combining (12) and (24). It should be noted that, according to Verruijt's research, ak, bk, ck, and dk cannot be obtained only by (12) and (24); another equation is need. This conclusion is also verified in practical calculation. Based on Verruijt’s work [7], another equation can be obtained. The value of a0 can be determined by assuming that a0=0, calculating the limiting value of ak when k is sufficiently large (e.g., ak=a, when k=1000), and repeating this calculation for an initial value a0=1 (e.g., ak=b, when k=1000). Then, the correct value of a0 is (25)a0=aa-b.


## 4. Solving Formulas for Stresses and Displacements

The stress components caused by excavation,σx, σy, and σxy, can be derived by (26)σy+σx=2φ′z+φ′z¯,σy-σx+2iσxy=2z-φ′′z+ψ′z,while the displacement components,ux and uy, are given by (27)(27)2Gux+iuy=3-4νφz-zφ′z¯-ψz¯.Here we assume that there is no vertical displacement at infinity. Herein, the vertical displacement actually is the relative vertical displacement to infinity, given as(28)u′x=ux,u′y=uy-uy∞,where u′x and u′y are the actual horizontal and vertical displacement, respectively.


## 5. Example and Discussion

In this section, the derived solution is applied to an example and a comparison is made using ABAQUS finite element code to verify the solution. Take the rectangle shaped tunnel witha=2m, b=2m, and h=10m as an example (see Figure 1).When the pipe jacking machine is used to excavate the soil, a slight overexcavation is necessary. The space between the cavity and excavation machine is maintained to facilitate the steering of the pipe jacking machine, installing pipes, and providing space for pumping slurry [17]. Although slurry is injected, there is still soil loss. Referring to the circular pipe jacking engineering experience [18], the rectangular pipe jacking has the uniform and nonuniform convergence model as shown in Figure 3. The uniform convergence mode is that segment suspended in the slurry, and the top clearance between the segment and soil is equal to bottom clearance. In the mode of nonuniform convergence, the whole segment sinks and the bottom clearance is 0. The rectangular pipe jacking has a large cross-section, and the gap between the bottom of the segment and the soil is often very small due to the weight. Thus, the displacement mode is generally nonuniform deformation in practical engineering. We assumed u¯=2v¯ in this paper; the design requires that the ground loss rate must be controlled within 5%; we have(29)2a×u¯+4b×v¯-2u¯×v¯2a×2b=5‰.Figure 3
Uniform and nonuniform convergence model.This givesu¯=10 mm and v¯=5 mm.We use the following parameters for calculation:E=10 MPa and ν=0.3.Form=10, the coefficients obtained according to the method given in Section 2 are shown in Table 1.Table 1
Coefficients determined for conformal mapping function.
Number
C
j
r
(
1
)
2.5196
0.1174
(
2
)
-
0.3004
(
3
)
0.0953
(
4
)
0.2129
(
5
)
0.2725
(
6
)
-
0.1859
(
7
)
0.0790
(
8
)
-
0.0014
(
9
)
0.0123
(
10
)
-
0.0177The specific expression of functionθ=f(α) is involved during the solving procedure of the Fourier series. Its form should be obtained accurately through the conformal mapping function. Since the conformal mapping function is highly complex, the solving is very difficult. As a simplification, θ=f(α) can be replaced approximately with a polynomial function. The corresponding relationship between θ and α and the polynomial fitting results are shown in Figure 4. There are nine polynomial items in this study. According to the figure, the accuracy of the fitting result is relatively high for most points and only a few points contain errors. It can be assumed that the accuracy of the fitting result increases with an increase in polynomial item quantity.Figure 4
The corresponding relationship betweenθ and α.The coefficient values and fitting results are shown in Table2. According to the table, the accuracy can meet the requirements when adopting a nine times polynomial fitting.Table 2
Polynomial fitting results ofθ=f(α).
Power
Coefficient
R-square
0
-
0.0669
0.9979
1
4.772299
2
-
13.0071
3
19.10024
4
-
15.0771
5
6.870741
6
-
1.85365
7
0.291853
8
-
0.02478
9
0.000876The displacement boundary condition of the tunnel can be expressed in the Fourier series form ofσ with (21), after the polynomial is obtained.In the actual solving procedure, the item quantities ofφ(ζ) and ψ(ζ), Ek, Gk, and Hk are impossible to be taken as infinity. Thus, relatively large numbers are generally used. Values of 200, 100, 100, and 100, respectively, are used in the solving procedure of this study.Here, the solution is compared with that obtained from finite element analysis using ABAQUS software. The same parameters as analytical solution are adopted.A comparison of ground surface settlement caused by tunnel construction is presented in Figure5, demonstrating good agreement between the analytical and finite element solutions. There is a difference of about 1% between the two. The distribution of displacement caused by the excavation of rectangular pipe jacking is approximately normal distribution, and the reverse bending point occurs at approximately 1.5 times of the excavation depth. Surface displacement occurs mainly within 4 times of the excavation depth; the settlement decreases with the increase of the depth; minimum settlement occurs at surface, with a value of about of 70% of the settlement of the tunnel boundary.Figure 5
Deformation of the surface.The stressesσx and σy of the above-mentioned example are shown in Figure 6. The negative sign here indicates compressive stress. It can be seen from Figure 6 that the analytical results are in good agreement with ABAQUS results. There is a difference of about 2% between the two. As shown in the figure, for σx, the maximum (35 KPa) appears at the floor and side corners and the minimum (−40 KPa) appears at the roof and side corners, while, for the stress σy, the maximum (25 KPa) appears at the roof and side corners and the minimum (−30 KPa) appears at the floor and side corners.Figure 6
Stresses curves at different position.
(a)
Horizontal stresses along lines aty=0 m, y=2 m, and y=4 m
(b)
Vertical stresses along lines aty=0 m, y=2 m, and y=4 mIn general, the finite element result can be seen as the accurate solution. The analytical solution of this study and the ABAQUS finite element code have no significant differences for most points. Most discrepancies between them are found at tunnel corners, which may be due to the fact that the grid size in the ABAQUS model is not sufficiently small. The analytical solution in this study is sufficiently accurate for the engineering application.


## 6. Conclusions

A complex variable solution for rectangle pipe jacking in elastic half-plane is derived in this paper. The conformal mapping function that transformed the half-planex<h with the exclusion of the hole into the concentric ring domain can be solved by the two-step method proposed in this paper. Through the mapping function proposed in this paper, both sides of the boundary condition equation can be transformed into a Laurent series expression. Then, the coefficient of the complex stress function can be obtained using a power series method. The accuracy of the result is validated by FEM results; there is a difference of about 2% between them. The analytical solution proposed in this paper has clear steps, has rapid convergence, has high accuracy, is easy to program, and is sufficiently accurate for this engineering application.



---

*Source: 5713063-2017-12-07.xml*