Study on the Horizontal Axis Deviation of a Small Radius TBM Tunnel Based on Winkler Foundation Model

Abstract

During the construction stage of the small radius TBM (tunnel boring machine) interval, the improper control of the boring parameters and the boring posture can cause the horizontal axis deviation of the shield tunnel. In order to address this issue, the TBM segments lining structure of the small radius interval is simplified as the continuous circular curved beam based on the longitudinal equivalent continuous model and Winkler elastic foundation beam theory. The theoretical model is solved through the transfer matrix method, and its applicability is verified by comparing it with the field monitoring data. It is found that the horizontal axis deviation of the completed tunnel increases with the total jack thrust, and the lateral displacement tends to be stable when the distance between the ring and the tail is far. The horizontal axis deviation has a negative relationship with the thrust difference or path difference when the jack thrust in the outside of the shield curve is larger than that of inside the shield curve. The horizontal axis deviation has a positive relationship with the thrust difference or path difference when the jack thrust in the outside of the shield curve is smaller than that of inside the shield curve.

1. Introduction

Tunnel boring machine (TBM) technology has been widely used for the construction of the urban subway tunnel. The longitudinal performance research of the TBM tunnel based on the equivalent continuous model and elastic foundation theory is a hotspot. The longitudinal affects the safe operation of the subway, which has also been widely studied through field monitoring, numerical analysis and analytical solution [1,2,3,4,5]. The research on the tunnel longitudinal performance focuses on the longitudinal settlement and horizontal axis deviation, which are caused by the coupling effects of multi-load [6,7,8,9].

Koizumi et al. [10] and Shiba et al. [11] were the pioneers developing the analytical solution for calculating the longitudinal settlement based on longitudinal bending stiffness. Elastic foundation models such as the Winkler model [12], the Pasternak model [13], and the Kerr model [14] are the most widely used theoretical calculation models in the longitudinal deformation research of tunnels. After that, the generalized longitudinal equivalent continuous model was proposed by Zhang et al. [15]. Li et al. [16] proposed the model considering the mechanical behavior of blots. Shiba et al. [17] and Talmon and Bezuijen [18] used the continuous elastic beam to simulate the longitudinal structure model of the shield tunnel, which can analyze the longitudinal property of the shield tunnel. Yu et al. [19] derived an analytical solution for the longitudinal bending stiffness of a segmental liner, which was mainly utilized in the shield tunnels. This method is verified by the simulation results from the finite element program.

The longitudinal settlement mechanism of the tunnel can also be analyzed by the two-dimensional and three-dimensional numerical simulation methods [20,21,22]. Huang et al. [23] developed a finite element method for analyzing the longitudinal performance of shield tunnels considering the longitudinal variation of geotechnical parameters. Yuan et al. [9] analyzed the deformation of the shield tunnel and the surrounding soil by building a three-dimensional model of shield tunnel. A numerical model was established to simulate the earth pressure balance (EPB) in excavation processes—the simulation results were compared with those obtained by the field measurement [24,25].

Also, field monitoring is an important method to analyze the longitudinal settlement of the shield tunnel. Ocak [26] studied the relationship between shield parameters and the shield surface settlement through field surface settlement measurements. Fargnoli et al. [27] collected the measurements of tunneling-induced settlements during the construction of the new Milan underground line 5, which were back analyzed with the method of the classical Gaussian empirical expression. Then, the detailed description of the EPB tunneling performance in the transverse and longitudinal directions could be achieved.

The research on the longitudinal performance of shield tunnels is mainly focused on the longitudinal settlement of tunnels. Limited research has been conducted on the problem of horizontal axis deviation of the small radius TBM interval construction.

In the design of urban subway lines, the sharp curve of the shield tunnel with a small curvature radius is usually adopted to avoid the adverse geological conditions and surrounding existing buildings [28,29]. When the curvature radius is less than 40 times of the TBM’s diameter, the construction difficulty and quality risk will be increased significantly. The damage caused by the horizontal axis deviation of the constructed tunnel is particularly prominent, which can easily affect the assembly quality of the segments. The horizontal axis deviation can result in tunnel mud leakage and local structural damage, which pose a severe threat to the safety of metro operation. Therefore, it is of great significance to study the effect of the construction load on the horizontal axis deviation of a small radius TBM interval tunnel.

The objective of this paper is to propose a solution of the tunnel horizontal axis deviation of small radius TBM interval with the transfer-matrix method, which is based on generalized longitudinal equivalent continuum model and Winkler elastic foundation beam theory. The applicability of the calculation model is verified through the monitored data in the field. The factors such as the jack thrust and shield posture affecting the tunnel horizontal axis deviation are analyzed.

2. Generalized Longitudinal Equivalent Continuous Model

2.1. Model Description

The prefabricated segment lining structure used in the TBM tunnel has a large number of longitudinal joints, which makes the deformation mechanism of the lining structure complicated. Various theoretical models for longitudinal deformation calculation have simplified the joint structures. Among these models, the longitudinal equivalent continuous model proposed by Shiba et al. [17] has the most extensive application. The tunnel cross-section is assumed as a homogeneous ring, and the tunnel longitudinal integral rigidity is also reduced when considering the joints. The tunnel is simplified as a uniform continuous beam with the equivalent stiffness on the elastic foundation.

Traditional models ignore the impacts of deformation of the tunnel cross-section and the limitation of the influential region of the longitudinal joints, which leads to the lower longitudinal effective rigidity ratio compared with the measured value. The generalized longitudinal equivalent continuous model [15] is adopted to calculate the longitudinal equivalent bending rigidity of the TBM tunnel segment lining structure in this paper. This model takes account of the influential region of the joints and lateral stiffness.

2.2. Basic Assumptions

Some assumptions are proposed in the generalized longitudinal equivalent continuous model:

(a)

The initial cross-section is uniform circular and develops into ellipse under the pressure of surrounding rock.

(b)

The cross-section before and after deformation is always plane under the longitudinal bending.

(c)

The stress distribution along the tunnel axial keeps unchanged on the neutral axis and inside the lining ring. One side of the neutral axis is under tension while the other side is under compression. The variation of the stress of lining structure along the thickness is ignored.

(d)

The bolts discretely distributed along the ring are transformed into continuously distributed bolt-ring, which can connect the annular gap. The stiffness of bolt-ring is uniformly distributed along the longitudinal length and radial thickness of the joints. The equivalent stiffness of bolt rings k_j is as follows [30].

$$k_j=\frac{n{k}_{ji}}{2\pi Rt}=\frac{n{E}_j{A}_j}{2\pi Rtl},$$

(1)

where k_ji is the elastic stiffness coefficient of a single bolt, n is the number of bolts inside the ring, R is the average of the internal and external diameter of the segment, t is the thickness of the segment, E_j is the elastic modulus of the bolt, A_j is the cross-section area of the bolt, and l is the interaction length of the bolt.

(e)

The influence coefficient of the annular gap λ is introduced, and the influence scope of the annular gap is λl. When λ < 1, the interaction length of ring direction joints within the influence scope of the annular gap is λl; when λ ≥ 1, the interaction length of ring direction joints within the influence scope of the annular gap is l.

(f)

Within the influence scope of the annular gap, the lining ring undertakes all the compressive stress caused by bending, and the tensile stress is jointly undertaken by equivalent bolt-ring and segment lining ring. Outside the influence scope of the annular gap, the tensile and compressive stresses caused by bending are both undertaken by the segment lining ring.

2.3. Longitudinal Equivalent Bending Stiffness

Under the longitudinal bending, when the tensions of the bolts are all smaller than the elastic limit P_y, the whole segment is in the elastic state. Figure 1 and Figure 2 show the longitudinal deformation and stress distribution of the segment inside and outside the influence scope of annular gaps. a and b are the long half axis and short half axis of the ellipse respectively, the ellipse is the cross-section of the segment under surrounding rock pressure, ϕ shows the position of neutral axis, c is the distance between neutral axis and y-axis, ds is the length of the micro-segmentation, dβ is the corresponding angle and x is the distance between neutral axis and the micro-segmentation, ε_t is the tension strain of the segment and ε_c is the compressive strain of the segment, M is the bending moment, E_c is the elastic modulus of tunnel reinforced concrete segment and δ is the displacement; ls is the segment length between the centerline of the two-segment rings and θ is the rotation angle.

Based on the basic assumptions, the equilibrium equation, the deformation compatibility equation, and the tunnel equivalent bending stiffness (EI)_eq can be obtained [31].

$${\left(EI\right)}_{eq}=\frac{M{l}_s}{\theta }=\frac{M{l}_s}{\frac{M}{\frac{E_c{Z}_1}{\lambda l}+\frac{E_c{k}_j{Z}_2}{\lambda E_c+\lambda k_{j}l}}+\frac{M\left(l_s-\lambda l\right)}{E_c{I}_c}}=E_c{I}_c\frac{l_s}{\frac{E_c{I}_c}{\frac{E_c{Z}_1}{\lambda l}+\frac{E_c{k}_j{Z}_2}{\lambda \left(E_c+k_{j}l\right)}}+\left(l_s-\lambda l\right)}.$$

(2)

The longitudinal stiffness effective rate η is [32],

$$\eta =\frac{{\left(EI\right)}_{eq}}{E_c{I}_c},$$

(3)

where [16],

$$\begin{gathered} Z_1=\frac{\left(a+b\right)b^{2}t\left(\frac{\pi }{2}-\phi \right)}{2}+\frac{{\left(a+b\right)}^{3}t{\sin }^2\phi \left(\frac{\pi }{2}-\phi \right)}{4}\hspace{0.33em}+\frac{\left(a+b\right)b^{2}t\sin 2\phi }{4}-\frac{{\left(a+b\right)}^{2}bt\sin 2\phi }{2} \\ {Z}_2=\frac{\left(a+b\right)b^{2}t\left(\frac{\pi }{2}+\phi \right)}{2}+\frac{{\left(a+b\right)}^{3}t{\sin }^2\phi \left(\frac{\pi }{2}+\phi \right)}{4}\hspace{0.33em}\hspace{0.33em}-\frac{\left(a+b\right)b^{2}t\sin 2\phi }{4}+\frac{{\left(a+b\right)}^{2}bt\sin 2\phi }{2} \end{gathered}$$

I_c is the inertia moment of tunnel segment cross-section; t is the thickness of lining ring.

3. Elastic Foundation Beam Model Used in the Small Radius TBM Interval

3.1. Model Description

The elastic foundation model is the most widely used theoretical model in the study of tunnel longitudinal deformation, mainly it has three categories: (a) Winkler foundation model [12], the foundation is regarded as a series of independent springs. The soil properties are manifested by the stiffness of the spring. (b) Pasternak model [13], an incompressible shear layer, only having horizontal shear deformation, is added in the Winkler foundation model. The shearing between springs is considered in the Pasternak model. (c) Kerr model [14], the spring layer, is added based on the Pasternak model. In contrast, the parameters in the Winkler foundation model are straightforward. When the compressible strata are thin and a hard layer exists, the calculation results based on the Winkler foundation model agree well with the real present situation. Therefore, the Winkler foundation model is selected to calculate the horizontal axis deviation of a small radius TBM tunnel subjected to construction loading. The lining structure of TBM segments in the small radius interval is simplified as the continuous circular curved beam on the Winkler foundation based on the equivalent elastic foundation beam theory (Figure 3). q_r is the radial linear load on the curved beam; q_t is the tangential linear load on the curved beam. Φ is the arbitrary angle on the circular curved beam; p is the concentrated force.

3.2. Basic Assumptions

Some basic assumptions are:

(a)

The small radius internal segment is regarded as the Euler-Bernoulli circular curved beam on the Winkler foundation.

(b)

The material of the curved beam is the linearly elastic body and conforms to the Hook’s law.

(c)

The curvature radius of the curved beam neutral axis keeps unchanged after deformation.

(d)

One end of the segment, away from shield tail and anchored in the surrounding rocks by grouting, is simplified as the fixed support.

(e)

The modified longitudinal equivalent continuous model is used to solve the longitudinal equivalent stiffness of the segment lining.

3.3. Model Solution

The uniform cross-section circular curved beam on the Winkler elastic foundation is depicted in Figure 4. The microelement body of the circular curved beam on the Winkler elastic foundation is shown in Figure 5. The equivalent equation of the microelement body, ignoring the high order, is shown below [32].

$$\frac{dQ}{dx}=p-\frac{N}{R}-q_r,$$

(4)

$$\frac{dN}{dx}=\frac{Q}{R}-q_t,$$

(5)

$$\frac{dM}{dx}=Q,$$

(6)

where Q is the shear force; N is the axial force.

For simplification, the effect of axial force on the deformation is neglected. Then, the deflection differential equation of the circular curved beam is obtained [34],

$$\left(\frac{d^{2}y}{d{x}^2}+\frac{y}{R^2}\right)+\frac{M}{EI}=0,$$

(7)

Combining Equations from (4) to (7), we can get deflection differential equation of the circular curved beam, which ignores the tangential deformation [34],

$$\frac{d^{5}y}{d{\mathsf{\Phi }}^5}+2\frac{d^{3}y}{d{\mathsf{\Phi }}^3}+{\mu }^2\frac{dy}{d\mathsf{\Phi }}=\frac{R^4}{EI}\left(\frac{d{q}_r}{d\mathsf{\Phi }}+q_t\right),$$

(8)

where EI is the bending stiffness of the tunnel, μ² is the coefficient, μ² = 1 + R⁴DK/EI, K is the coefficient of foundation bedding, D is the width of the curved beam cross-section, R is the curvature radius of the curved beam.

When there is no load on the beam span, the Equation (8) can be changed into,

$$\frac{d^{5}y}{d{\mathsf{\Phi }}^5}+2\frac{d^{3}y}{d{\mathsf{\Phi }}^3}+{\mu }^2\frac{dy}{d\mathsf{\Phi }}=0.$$

(9)

This is the homogenous fifth-order differential equation with constant coefficients, whose general solution is as below [34],

$$y\left(\mathsf{\Phi }\right)=C_1+\left(C_2\mathrm{ch}\alpha \mathsf{\Phi }+C_3\mathrm{sh}\alpha \mathsf{\Phi }\right)\cos \beta \mathsf{\Phi }+\left(C_4\mathrm{ch}\alpha \mathsf{\Phi }+C_5\mathrm{sh}\alpha \mathsf{\Phi }\right)sin\beta \mathsf{\Phi },$$

(10)

where,

$$\alpha =\sqrt{\frac{\left(\mu -1\right)}{2}},\beta =\sqrt{\frac{\left(\mu +1\right)}{2}};$$

C₁–C₅ are undetermined coefficients, obtained by loading and boundary conditions.

With the successive derivation of Equation (10), we can obtain analytical expressions of the rotation angle θ, the bending moment M, the shear force Q, the axial force N. With the initial condition of the beam, we can get the solutions of undetermined coefficients C₁–C₅, the general solution of the defection differential equation can be converted into the below equation [34],

$$y\left(\mathsf{\Phi }\right)=y_0{F}_{y1}\left(\mathsf{\Phi }\right)+{\theta }_0{F}_{y2}\left(\mathsf{\Phi }\right)+M_0{F}_{y3}\left(\mathsf{\Phi }\right)+Q_0{F}_{y4}\left(\mathsf{\Phi }\right)\hspace{0.33em}\hspace{0.33em}+N_0{F}_{y5}\left(\mathsf{\Phi }\right),$$

(11)

where

$$\begin{array}{l}\begin{array}{ll}\hfill & \left\{\begin{array}{c}{F}_{y1}\left(\mathsf{\Phi }\right)=\frac{2\alpha \beta }{1+2\alpha \beta }+\frac{1}{1+2\alpha \beta }{W}_1\left(\mathsf{\Phi }\right)+\frac{R^{4}kD}{EI\left(1+2\alpha \beta \right)}{W}_2\left(\mathsf{\Phi }\right)\\ F_{y2}\left(\mathsf{\Phi }\right)=\frac{\alpha }{{\alpha }^2+{\beta }^2}{W}_2\left(\mathsf{\Phi }\right)+\frac{\beta }{{\alpha }^2+{\beta }^2}{W}_3\left(\mathsf{\Phi }\right)\hfill \\ F_{y3}\left(\mathsf{\Phi }\right)=\frac{R^2}{\left(1+2\alpha \beta \right)EI}{W}_1\left(\mathsf{\Phi }\right)-\frac{R^{6}kD}{2\alpha \beta \left(1+2\alpha \beta \right){\left(EI\right)}^2}{W}_4\left(\mathsf{\Phi }\right)\\ F_{y4}\left(\mathsf{\Phi }\right)=\frac{R^3}{2EIa\left({\alpha }^2+{\beta }^2\right)}{W}_2\left(\mathsf{\Phi }\right)-\frac{R^3}{2EI\beta \left({\alpha }^2+{\beta }^2\right)}{W}_3\left(\mathsf{\Phi }\right)\\ F_{y5}\left(\mathsf{\Phi }\right)=-\frac{R^3}{2\alpha \beta EI}{W}_4\left(\mathsf{\Phi }\right)\hfill \end{array}\right.\hfill \\ \hfill & \left\{\begin{array}{c}\begin{array}{c}{W}_1\left(\mathsf{\Phi }\right)=ch\left(\alpha \mathsf{\Phi }\right)cos\left(\beta \mathsf{\Phi }\right)\\ W_2\left(\mathsf{\Phi }\right)=sh\left(\alpha \mathsf{\Phi }\right)cos\left(\beta \mathsf{\Phi }\right)\\ W_3\left(\mathsf{\Phi }\right)=ch\left(\alpha \mathsf{\Phi }\right)sin\left(\beta \mathsf{\Phi }\right)\\ W_4\left(\mathsf{\Phi }\right)=sh\left(\alpha \mathsf{\Phi }\right)sin\left(\beta \mathsf{\Phi }\right)\end{array}\end{array}\right.\hfill \end{array}\end{array}$$

When the arbitrary point on the curved span has the concentrated moment M_i, the concentrated force P_i, the radial distribution load q_ri, the tangential distribution load q_ti, the external load on the beam span can be regarded as the partial initial parameters. Therefore, additional deflection caused by the concentrated moment M_i is as below [34],

$$y{\left(\mathsf{\Phi }\right)}_{M_i}=M_i{F}_{y3}\left(\mathsf{\Phi }-{\mathsf{\Phi }}_{M_i}\right)\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\left(\mathsf{\Phi }\ge {\mathsf{\Phi }}_{M_i}\right),$$

(12)

The additional defection caused by the concentrated force P_i when Φ ≥ Φ_Pi is shown below [34],

$$y{\left(\mathsf{\Phi }\right)}_{P_i}=-P_i{F}_{y4}\left(\mathsf{\Phi }-{\mathsf{\Phi }}_{P_i}\right)\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\left(\mathsf{\Phi }\ge {\mathsf{\Phi }}_{P_i}\right),$$

(13)

The additional defection caused by the radial distribution load q_ri when Φ ≥ Φ_ai can be obtained [34],

$$y{\left(\mathsf{\Phi }\right)}_{q_{ti}}=-D{\displaystyle {\int }_{{\mathsf{\Phi }}_{ci}}^{\mathsf{\Phi }}{q}_t\left(\delta \right)F_{y5}\left(\mathsf{\Phi }-\delta \right)d\delta }\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\left(\mathsf{\Phi }\ge {\mathsf{\Phi }}_{ci}\right),$$

(14)

When Φ ≥ Φ_bi, the integral upper limit Φ = Φ_bi is selected.

The additional defection caused by the tangential distribution load q_ti when Φ ≥ Φ_ci can also be obtained [34],

$$y{\left(\mathsf{\Phi }\right)}_{q_t}=-D{\displaystyle {\int }_{{\mathsf{\Phi }}_c}^{\mathsf{\Phi }}{q}_t\left(\delta \right)F_{y5}\left(\mathsf{\Phi }-\delta \right)d\delta }\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\left(\mathsf{\Phi }\ge {\mathsf{\Phi }}_c\right),$$

(15)

When Φ > Φ_di, the integral upper limit Φ = Φ_di is selected.

With the initial conditions while Φ = 0 and the above external load, the deflection of the curved beam at any points can be achieved [34],

$$\begin{array}{ll}y\left(\mathsf{\Phi }\right)=& y_0{F}_{y1}\left(\mathsf{\Phi }\right)+{\theta }_0{F}_{y2}\left(\mathsf{\Phi }\right)+M_0{F}_{y3}\left(\mathsf{\Phi }\right)+Q_0{F}_{y4}\left(\mathsf{\Phi }\right)+N_0{F}_{y5}\left(\mathsf{\Phi }\right)+{\displaystyle \sum _{i=1}^{n_M}{M}_i{F}_{y3}\left(\mathsf{\Phi }-{\mathsf{\Phi }}_{M_i}\right)}\\ & -{\displaystyle \sum _{i=1}^{n_P}{P}_i{F}_{y4}\left(\mathsf{\Phi }-{\mathsf{\Phi }}_{P_i}\right)}-{\displaystyle \sum _{i=1}^{n_{q_r}}D{\displaystyle {\int }_{{\mathsf{\Phi }}_{ai}}^{\mathsf{\Phi }}{q}_{ri}\left(\delta \right)F_{y4}\left(\mathsf{\Phi }-\delta \right)d\delta }}-{\displaystyle \sum _{i=1}^{n_{q_t}}D{\displaystyle {\int }_{{\mathsf{\Phi }}_{ci}}^{\mathsf{\Phi }}{q}_{ti}\left(\delta \right)F_{y5}\left(\mathsf{\Phi }-\delta \right)d\delta }}\end{array}$$

(16)

With the same procedure, we can obtain the analytical expressions of rotation angle θ, bending moment M, shear force Q, and axial force N when some loads are acting on the beam span.

The TBM tunnel segment experiences leaving the shield tail, pea gravel backfilling, and grouting in turn. There are differences in the stiffness of the medium near the segment and the deformation of segment cross-section with loading. Therefore, the bending of the limited length variable cross-section circular curved beam on the inhomogeneous Winkler foundation needs further discussion.

Figure 6 shows the equal section circular curved beam subjected to the distributed load on the homogeneous Winkler foundation. The analytical expressions of the displacement and the internal force can be expressed in the following matrix form [34].

$$A{\left(\mathsf{\Phi }\right)}_{5\times 1}=B{\left(\mathsf{\Phi }\right)}_{5\times 5}A{\left(0\right)}_{5\times 1}+C{\left(\mathsf{\Phi }\right)}_{5\times 1}+D{\left(\mathsf{\Phi }\right)}_{5\times 1},$$

(17)

where,

$$\begin{gathered} A{\left(\mathsf{\Phi }\right)}_{5\times 1}={\left[y\left(\mathsf{\Phi }\right),\theta \left(\mathsf{\Phi }\right),M\left(\mathsf{\Phi }\right),Q\left(\mathsf{\Phi }\right),N\left(\mathsf{\Phi }\right)\right]}^T; \\ B{\left(\mathsf{\Phi }\right)}_{5\times 5}=\left[\begin{array}{ccccc}{F}_{y1}\left(\mathsf{\Phi }\right)& F_{y2}\left(\mathsf{\Phi }\right)& F_{y3}\left(\mathsf{\Phi }\right)& F_{y4}\left(\mathsf{\Phi }\right)& F_{y5}\left(\mathsf{\Phi }\right)\\ F_{\theta 1}\left(\mathsf{\Phi }\right)& F_{\theta 2}\left(\mathsf{\Phi }\right)& F_{\theta 3}\left(\mathsf{\Phi }\right)& F_{\theta 4}\left(\mathsf{\Phi }\right)& F_{\theta 5}\left(\mathsf{\Phi }\right)\\ F_{M1}\left(\mathsf{\Phi }\right)& F_{M2}\left(\mathsf{\Phi }\right)& F_{M3}\left(\mathsf{\Phi }\right)& F_{M4}\left(\mathsf{\Phi }\right)& F_{M5}\left(\mathsf{\Phi }\right)\\ F_{Q1}\left(\mathsf{\Phi }\right)& F_{Q2}\left(\mathsf{\Phi }\right)& F_{Q3}\left(\mathsf{\Phi }\right)& F_{Q4}\left(\mathsf{\Phi }\right)& F_{Q5}\left(\mathsf{\Phi }\right)\\ F_{N1}\left(\mathsf{\Phi }\right)& F_{N2}\left(\mathsf{\Phi }\right)& F_{N3}\left(\mathsf{\Phi }\right)& F_{N4}\left(\mathsf{\Phi }\right)& F_{N5}\left(\mathsf{\Phi }\right)\end{array}\right]; \\ A{\left(0\right)}_{5\times 1}={\left[y_0,{\theta }_0,M_0,Q_0,N_0\right]}^T; \\ C{\left(\mathsf{\Phi }\right)}_{5\times 1}={\left[C_y\left(\mathsf{\Phi }\right),C_{\theta }\left(\mathsf{\Phi }\right),C_M\left(\mathsf{\Phi }\right),C_Q\left(\mathsf{\Phi }\right),C_N\left(\mathsf{\Phi }\right)\right]}^T; \\ D{\left(\mathsf{\Phi }\right)}_{5\times 1}={\left[D_y\left(\mathsf{\Phi }\right),D_{\theta }\left(\mathsf{\Phi }\right),D_M\left(\mathsf{\Phi }\right),D_Q\left(\mathsf{\Phi }\right),D_N\left(\mathsf{\Phi }\right)\right]}^T; \end{gathered}$$

The elements in the matrix C(Φ)_5×1 are as follows [34],

$$\left\{\begin{array}{c}{C}_y\left(\mathsf{\Phi }\right)=-D{\displaystyle {\int }_{{\mathsf{\Phi }}_a}^{\mathsf{\Phi }}{q}_r\left(\delta \right)F_{y4}\left(\mathsf{\Phi }-\delta \right)d\delta }\\ C_{\theta }\left(\mathsf{\Phi }\right)=-D{\displaystyle {\int }_{{\mathsf{\Phi }}_a}^{\mathsf{\Phi }}{q}_r\left(\delta \right)F_{\theta 4}\left(\mathsf{\Phi }-\delta \right)d\delta }\\ C_M\left(\mathsf{\Phi }\right)=-D{\displaystyle {\int }_{{\mathsf{\Phi }}_a}^{\mathsf{\Phi }}{q}_r\left(\delta \right)F_{M4}\left(\mathsf{\Phi }-\delta \right)d\delta }\\ C_Q\left(\mathsf{\Phi }\right)=-D{\displaystyle {\int }_{{\mathsf{\Phi }}_a}^{\mathsf{\Phi }}{q}_r\left(\delta \right)F_{Q4}\left(\mathsf{\Phi }-\delta \right)d\delta }\\ C_N\left(\mathsf{\Phi }\right)=-D{\displaystyle {\int }_{{\mathsf{\Phi }}_a}^{\mathsf{\Phi }}{q}_r\left(\delta \right)F_{N4}\left(\mathsf{\Phi }-\delta \right)d\delta }\end{array}\right.;$$

The elements in the matrix D(Φ)_5×1 are as follows [34],

$$\left\{\begin{array}{c}{D}_y\left(\mathsf{\Phi }\right)=-D{\displaystyle {\int }_{{\mathsf{\Phi }}_a}^{\mathsf{\Phi }}{q}_t\left(\delta \right)F_{y5}\left(\mathsf{\Phi }-\delta \right)d\delta }\\ D_{\theta }\left(\mathsf{\Phi }\right)=-D{\displaystyle {\int }_{{\mathsf{\Phi }}_a}^{\mathsf{\Phi }}{q}_t\left(\delta \right)F_{\theta 5}\left(\mathsf{\Phi }-\delta \right)d\delta }\\ D_M\left(\mathsf{\Phi }\right)=-D{\displaystyle {\int }_{{\mathsf{\Phi }}_a}^{\mathsf{\Phi }}{q}_t\left(\delta \right)F_{M5}\left(\mathsf{\Phi }-\delta \right)d\delta }\\ D_Q\left(\mathsf{\Phi }\right)=-D{\displaystyle {\int }_{{\mathsf{\Phi }}_a}^{\mathsf{\Phi }}{q}_t\left(\delta \right)F_{Q4}\left(\mathsf{\Phi }-\delta \right)d\delta }\\ D_N\left(\mathsf{\Phi }\right)=-D{\displaystyle {\int }_{{\mathsf{\Phi }}_a}^{\mathsf{\Phi }}{q}_t\left(\delta \right)F_{N4}\left(\mathsf{\Phi }-\delta \right)d\delta }\end{array}\right.;$$

The curved beam is divided into n sections into the abruptly changed sections of bending stiffness and foundation bedding coefficient, in the action points of the concentrated force and moment, and the starting action points of the distributed loading. There are concentrated moments M_i and concentrated forces P_i at the junctures Φ = Φ_i (i = 0, 1, 2, …, n) of the adjacent curved beams. Each curved beam section has radial distributed load q_ri(Φ) and tangential distributed load q_ti(Φ).

For curved beam i, in the interval Φ_i−1⁺ ≤ Φ < Φ_i⁻,

$$\begin{array}{ll}{A}_i{\left({\mathsf{\Phi }}^{-}\right)}_{5\times 1}& {=B}_i{\left(\mathsf{\Phi }-{\mathsf{\Phi }}_{i-1}\right)}_{5\times 5}{A}_i{\left({\mathsf{\Phi }}_{i-1}^{+}\right)}_{5\times 1}+C_i{\left(\mathsf{\Phi }-{\mathsf{\Phi }}_{i-1}\right)}_{5\times 1}\\ & +D_i{\left(\mathsf{\Phi }-{\mathsf{\Phi }}_{i-1}\right)}_{5\times 1}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\left(i=1,2,\dots ,n-1,n\right)\end{array},$$

(18)

where the α, β, K, EI in the factor expression should be corresponding to the α_i, β_i, K_i, and (EI_eq)_i in section i.

The interface between section i−1 and section i, the cross-section Φ_i−1, the Equation based on the deformation consistency and force equilibrium is as follows,

$$A_i{\left({\mathsf{\Phi }}_{i-1}^{+}\right)}_{5\times 1}=A_{i-1}{\left({\mathsf{\Phi }}_{i-1}^{-}\right)}_{5\times 1}+E_{i-1}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\left(i=2,3,\dots ,n-1,n\right),$$

(19)

where,

$$E_{i-1}={\left[0,0,M_{i-1},P_{i-1},0\right]}^T\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\left(i=2,3,\dots ,n-1,n\right);$$

The calculation results at the end of curved beam i−1 can be regarded as the initial condition of the curved beam i (2 ≤ i ≤ n), then we can obtain,

$$A_i{\left({\mathsf{\Phi }}^{-}\right)}_{5\times 1}=\tilde{B}{\left({\mathsf{\Phi }}^{-}\right)}_{5\times 5}{A}_1{\left({\mathsf{\Phi }}_0^{+}\right)}_{5\times 1}+\tilde{C}\left({\mathsf{\Phi }}^{-}\right)+\tilde{D}\left({\mathsf{\Phi }}^{-}\right)\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\left(i=2,3,\dots ,n-1,n\right),$$

(20)

where,

$$\tilde{B}\left(\mathsf{\Phi }\right)=B_i\left(\mathsf{\Phi }-{\mathsf{\Phi }}_{i-1}\right)B_{i-1}\left({\mathsf{\Phi }}_{i-1}-{\mathsf{\Phi }}_{i-2}\right)\mathrm{\dots }{B}_1\left({\mathsf{\Phi }}_1-{\mathsf{\Phi }}_0\right);$$

$$\begin{array}{ll}\tilde{C}\left(\mathsf{\Phi }\right)=& C_i\left(\mathsf{\Phi }-{\mathsf{\Phi }}_{i-1}\right)+B_i\left(\mathsf{\Phi }-{\mathsf{\Phi }}_{i-1}\right)\left[C_{i-1}\left({\mathsf{\Phi }}_{i-1}-{\mathsf{\Phi }}_{i-2}\right)+E_{i-1}\right]\\ & +B_i\left(\mathsf{\Phi }-{\mathsf{\Phi }}_{i-1}\right)B_{i-1}\left({\mathsf{\Phi }}_{i-1}-{\mathsf{\Phi }}_{i-2}\right)\left[C_{i-2}\left({\mathsf{\Phi }}_{i-2}-{\mathsf{\Phi }}_{i-3}\right)+E_{i-2}\right]\\ & +\dots \\ & +B_i\left(\mathsf{\Phi }-{\mathsf{\Phi }}_{i-1}\right)B_{i-1}\left({\mathsf{\Phi }}_{i-1}-{\mathsf{\Phi }}_{i-2}\right)\dots B_3\left({\mathsf{\Phi }}_3-{\mathsf{\Phi }}_2\right)\times \left[C_2\left({\mathsf{\Phi }}_2-{\mathsf{\Phi }}_1\right)+E_2\right]\\ & +B_i\left(\mathsf{\Phi }-{\mathsf{\Phi }}_{i-1}\right)B_{i-1}\left({\mathsf{\Phi }}_{i-1}-{\mathsf{\Phi }}_{i-2}\right)\dots B_3\left({\mathsf{\Phi }}_3-{\mathsf{\Phi }}_2\right)B_2\left({\mathsf{\Phi }}_2-{\mathsf{\Phi }}_1\right)\times \left[C_1\left({\mathsf{\Phi }}_1-{\mathsf{\Phi }}_0\right)+E_1\right]\end{array};$$

$$\begin{array}{ll}\tilde{D}\left(\mathsf{\Phi }\right)=& D_i\left(\mathsf{\Phi }-{\mathsf{\Phi }}_{i-1}\right)\\ & +B_i\left(\mathsf{\Phi }-{\mathsf{\Phi }}_{i-1}\right)\left[D_{i-1}\left({\mathsf{\Phi }}_{i-1}-{\mathsf{\Phi }}_{i-2}\right)+E_{i-1}\right]\\ & +B_i\left(\mathsf{\Phi }-{\mathsf{\Phi }}_{i-1}\right)B_{i-1}\left({\mathsf{\Phi }}_{i-1}-{\mathsf{\Phi }}_{i-2}\right)\left[C_{i-2}\left({\mathsf{\Phi }}_{i-2}-{\mathsf{\Phi }}_{i-3}\right)+E_{i-2}\right]\\ & +\dots \\ & +B_i\left(\mathsf{\Phi }-{\mathsf{\Phi }}_{i-1}\right)B_{i-1}\left({\mathsf{\Phi }}_{i-1}-{\mathsf{\Phi }}_{i-2}\right)\dots B_3\left({\mathsf{\Phi }}_3-{\mathsf{\Phi }}_2\right)\times \left[D_2\left({\mathsf{\Phi }}_2-{\mathsf{\Phi }}_1\right)+E_2\right]\\ & +B_i\left(\mathsf{\Phi }-{\mathsf{\Phi }}_{i-1}\right)B_{i-1}\left({\mathsf{\Phi }}_{i-1}-{\mathsf{\Phi }}_{i-2}\right)\dots B_3\left({\mathsf{\Phi }}_3-{\mathsf{\Phi }}_2\right)B_2\left({\mathsf{\Phi }}_2-{\mathsf{\Phi }}_1\right)\times \left[D_1\left({\mathsf{\Phi }}_1-{\mathsf{\Phi }}_0\right)+E_1\right]\end{array};$$

When i = n and Φ₀⁺ = 0, Φn⁻ = Φ_max, the following equation can be achieved based on Equation (20),

$$A{\left({\mathsf{\Phi }}_{\max }\right)}_{5\times 1}=\tilde{B}{\left({\mathsf{\Phi }}_{\max }\right)}_{5\times 5}A{\left(0\right)}_{5\times 1}+\tilde{C}\left({\mathsf{\Phi }}_{\max }\right)+\tilde{D}\left({\mathsf{\Phi }}_{\max }\right).$$

(21)

The unknown initial parameters in matrix A(0)_5×1 can be solved with Equation (21) based on the boundary conditions of the two ends of the curved beam. Putting all the initial parameters at the initial section into the Equation (20), we can obtain the deflection, rotation angle, moment, shear force and axial force at the arbitrary angle Φ of the circular curved beam when the foundation bedding coefficient and bending stiffness stepped change.

4. Analysis of Engineering Example

4.1. Engineering Description

The entrance/exit tunnel MRDK0+457.8~688.2 of Min-Le parking lot belongs to the Shenzhen rail transit line 6 phase II, which is a small radius TBM internal with the curvature radius R = 300 m (as illustrated in Figure 7).

The reinforced concrete (C50) universal wedge-shaped segments with a ring width of 1.2 m are utilized for the tunnel lining structure. The outer and inner diameters are 6200 mm and 5400 mm, respectively, and the thickness is 400 mm. The rings are connected by ten longitudinal bending bolts (M24). There are a total of 192 rings in the interval. The No. 1 and 2 segments are wrapped by the shield tail brush and sealing materials, the elastic resistance coefficient of foundation at the shield tail k₁ = 100 kPa/m; The No. 3 to 5 segments have finished the backfilling of the bottom part by pea gravel, then the elastic resistance coefficient of foundation k₂ = 292 MPa/m; The No. 6 to 10 segments have already finished the backfilling of the sidewall by pea gravel, then the elastic resistance coefficient of foundation k₃ = 1701 MPa/m; The No. 11 to 96 segments are in the slightly weathered cataclasite formation, then the elastic resistance coefficient of foundation k₄ = 1993 MPa/m; The No. 97 to 120 segments are in the slighted weathered granite formation, the elastic resistance coefficient of foundation k₅ = 4770 MPa/m; The No. 121 to 192 segments are in the section where the grout is not solidified, the elastic resistance coefficient of foundation k₆ is linearly changed, and finally reaching to 10,309 MPa/m. The tangential resistance in the model is provided by the friction between the segments and the surrounding medium. The frictional coefficient between the shield tail brush and the segment is 0.3, and the frictional coefficient between the pea gravel and the segment is 0.5. The action sphere of the tangential resistance can be solved through condition N(Φ_qt) = 0.

The influence range coefficient of the annular gap λ is set as 0.2 in this paper, and the generalized longitudinal equivalent continuous model is utilized to calculate the longitudinal equivalent bending stiffness (EI)_eq in the tunnel. The grouting pea gravels haven’t formed the ring in the wall back of No. 1 to 10 segments, and the segments haven’t undertaken the surrounding rock pressure, then the bending stiffness (EI)_eq1 = 162.7 GPa·m⁴; The vertical uniform pressure of No. 11 to 96 segments is 107.7 kPa, the base uniform reaction is 110.8 kPa, the horizontal uniform pressure is 32.3 kPa, and the bending stiffness (EI)_eq2 = 163.1 GPa·m⁴; The vertical uniform pressure of No. 97 to 192 segments is 54.9 kPa, the base uniform reaction is 58.0 kPa, the horizontal uniform pressure is 8.2 kPa, and the bending stiffness (EI)_eq3 = 162.8 GPa·m⁴.

The signs of the force and moment are stipulated for the convenience of further study. When the horizontal jack thrust on the outside of the shield curve is larger than that on the inside of the shield curve, the horizontal force couple is positive; when the path of the shield on the outside of the shield curve is larger than that on the inside of the shield curve, the jack thrust deviation angle is positive When the segment’s lateral displacement points to the circle center of the circular curve, this direction is positive, and vice versa. Through the monitoring data, the average of the total auxiliary cylinder jack thrust is 5931.33 kN, the average of the jack thrust on the outside of the shield curve is larger than that on the inside of the shield curve. The average of the force couple caused by the difference of the jack thrust in the horizontal direction is 1708.20 kN·m. The average of the lateral component caused by the path difference of the auxiliary cylinder jack is 40.36 kN.

4.2. The Accumulative Value of the Segments’ Lateral Displacements

The calculated lateral displacements of the ring segments and the monitored lateral displacements of the ring segments from the field are shown in Figure 8, the parameters used in this calculation are shown in

Table 1. Initial parameters for calculation model.

Total Thrust F/(kN)	Lateral Force Couple M/(kN·m)	Deviation Angle α/(°)	Monitoring Data
			N0/(kN)	M0/(kN·m)	Q0/(kN)
5931.47	−1708.20	−0.390	5931.33	−1708.20	−40.36

. From Figure 8, we understand the following.

(a)

When the segment is in the shield tail, there is a big difference between (2.267 mm) the simulated displacements and the field monitoring results. This is probably due to the simplification of the wrapping influence of the shield tail brush and sealing grease on the shield tail segments during the theoretical simulation.

(b)

In the process of segments leaving the shield tail around 1 to 4 rings, the simulated accumulated lateral displacement agrees well with the lateral displacement from the field. Compared with the monitored data, the simulated results are a little underestimated. The curvature radius in this project has already reached the turning limitation of the TBM. The poorly controlled shield posture and the improper selection of assembly position can cause the squeezing action of the second ring in the shield tail from the shield shell or the tail brush. The squeezing action also causes the increase of lateral offset of some rings, which just leaving the tail and are backfilled. The first ring leaving the tail is most influenced by this effect. From the field data, we can know that the lateral displacement increases rapidly when the ring leaves the tail.

(c)

In the process of segments leaving the shield tail around 5 to 8 rings, the simulated displacements coincide quite well with the monitored data. When the backfilling finishes, the accumulated displacements from the field tend to be stable. However, the simulated displacements still increase with low speed.

(d)

When the segments leaving the tail around nine rings, there is a big difference between theoretical results and field data. When the pea gravels form the ring, the lateral accumulated displacement declines slowly under the pressure difference of surrounding rocks and tends to be stable from the 17th ring. For the theoretical simulation, the pressure difference of surrounding rocks from both sides of the shield line is ignored as the 3D tunnel is simplified as the 2D circular curved beam in the model. The simulated accumulated displacement increases slowly from the 11th ring and tends to be stable from the 15th ring.

Finally, the simulated displacements agree well with the field data; the simulated results are more conservative. Therefore, this model is applicable in the prediction of the horizontal accumulated displacement of the shield segments in the small radius interval.

5. Single-Factor Influence Analysis of the Segments’ Horizontal Displacements Based on the Model

5.1. The Effect of the Jack Total Thrust

The total jack thrust is set as N = 6000 kN, 7000 kN, 8000 kN, 9000 kN, and 10,000 kN, respectively, during the calculation in the model while other parameters are kept the same. The parameters in the theoretical model are summarized in

Table 2. The initial parameters under different total jack thrust.

Total Thrust F/(kN)	Lateral Force Couple M/(kN·m)	Deviation Angle α/(°)	N0/(kN)	M0/(kN·m)	Q0/(kN)
6000	−1000	−1	5999.09	−1000	−104.71
7000			6998.93		−122.17
8000			7998.78		−139.62
9000			8998.63		−157.07
10,000			9998.48		−174.52

. Then, the evolution law of the total thrust on the lateral accumulated displacements of a single ring is analyzed.

Figure 9 shows that the lateral accumulated displacement increases with the increasing distance between the ring and the shield tail. With the increase of the leaving distance (completed rings), the lateral accumulated displacements tend to be stable. The segment results in the deformation along the axis under the jack thrust. When the outside of the segment contacts with the surrounding medium, the friction is produced to resist the jack thrust. When the frictional resistance of the first ring inside the shield tail cannot resist the jack thrust, the remaining thrust can be delivered to the next rings. The transverse component of the remaining thrust can cause a horizontal axis deviation. Therefore, with the increase of the total jack thrust, the remaining thrust acting on the segments increase, and the horizontal axis deviation also increases.

Based on the stable criteria that the displacement difference of the adjacent ring is no more than 0.1 mm, the stable ring number that was leaving the tail increase from the 12th to the 17th when the total thrust increases from 6000 kN to 10,000 kN. In the stable state, the minimum lateral accumulated displacement is −15.347 mm, while the maximum lateral accumulated displacement is −28.321 mm, increasing by 84.54%.

5.2. The Effect of the Jack Thrust Misalignment

The lateral force couple M (2000 kN·m, 1000 kN·m, 0 kN·m, −1000 kN·m, and −2000 kN·m) caused by the thrust difference of the right and left jacks is selected as the variable parameters. The parameters for calculated are summarized in

Table 3. The initial parameters under different jack thrust misalignments.

Total Thrust F/(kN)	Lateral Force Couple M/(kN·m)	Deviation Angle α/(°)	N0/(kN)	M0/(kN·m)	Q0/(kN)
6000	−2000	−1	5999.09	−2000	−104.71
	−1000			−1000
	0			0
	1000			1000
	2000			2000

. The lateral accumulated displacement of a single ring is simulated based on the model.

Figure 10 illustrates that the lateral accumulated displacement increases with the distance of the shield tail, and tends to be stable when leaving the rail around 12 rings. The effect of the jack thrust on the shield ring will decay with the increasing distance, causing the decreasing lateral accumulated displacement. The direction of displacements is leaving the circular center; the displacement when the inside force > the outside force is larger than that when the outside force > the inside force. When the force couple is positive, the increasing force couple will prevent the lateral deviation of the shield efficiently.

5.3. The Effect of the Jack Thrust Deviation Angle

The deviation angle α caused by the path difference of lateral jacks is set as control variables. The lateral accumulated displacements of a single ring with different deviation angles (3°, 2°, 1°, 0°, −1°, −2°, and −3°) are simulated with the model, and the parameters used in the calculation are summarized in

Table 4. The initial parameters under different jack thrust deviation angles.

Total Thrust F/(kN)	Lateral Force Couple M/(kN·m)	Deviation Angle α/(°)	N0/(kN)	M0/(kN·m)	Q0/(kN)
6000	−1000	−3	5991.78	−1000	−314.02
		−2	5996.34		−209.40
		−1	5999.09		−104.71
		0	6000.00		0.00
		1	5999.09		104.71
		2	5996.34		209.40
		3	5991.78		314.02

.

As depicted in Figure 11, the lateral accumulated displacement increases with the distance between the ring and the tail and tends to be stable when the tail leaves the 12th ring. In the construction of the small radius TBM tunnel, the path of the external jack is always larger than that of the internal jack, which causes the deviation angle is negative. Therefore, the axial of the ring will have an external deviation unavoidably. The situation that the external path is larger than the internal path should be avoided in construction, especially in the process of a regripping cycle.

6. Conclusions

(a)

Based on the generalized longitudinal equivalent continuous model and the elastic foundation beam theory, the horizontal axis deviation of a small radius TBM tunnel is simplified as the deflection of the circular curved beam on the Winkler foundation.

(b)

The model has been solved by the transfer-matrix method, which can solve the deformation and internal force with high accuracy at any cross-section conveniently and fast. When the foundation bedding coefficient and bending stiffness stepped change, the exact solutions of the model can be achieved. When the foundation bedding coefficient and bending stiffness change continuously, the numerical solutions meeting accuracy requirements can be obtained by increasing the number of the order.

(c)

In a setting of Shenzhen rail transit line 6 phase II, the calculation model of segment axial horizontal deviation for the small radius TBM tunnel construction is established. The applicability of this model is verified by the field data.

(d)

The effect of the jack thrust on the lateral deviation of the tunnel segments decreases with the increase of the distance between the shield tail and the ring. The influence range depends on the total thrust and increases with the total thrust.

(e)

When the lateral thrust in the external side of the curved line is larger than that in the internal side of the curved line, the horizontal axis of the segments moves to the external side, and the offset has a positive relationship with the thrust difference (or path difference). When the lateral thrust in the external side of the curved line is smaller than that in the internal side of the curved line, the horizontal axis of the segments still moves to the external side. However, the offset decreases with the increase of thrust difference (or path difference).