Open AccessArticle

Random Responses of Shield Tunnel to New Tunnel Undercrossing Considering Spatial Variability of Soil Elastic Modulus

School of Civil Engineering and Architecture, Zhejiang Sci-Tech University, Hangzhou 310018, China

Research Centre of Coastal and Urban Geotechnical Engineering, Zhejiang University, Hangzhou 310058, China

Laboratory of Geotechnics, Department of Civil Engineering, Ghent University, 9052 Ghent, Belgium

⁴

Deltares, Geo-Engineering, P.O. Box 177, 2600 MH Delft, The Netherlands

Author to whom correspondence should be addressed.

Appl. Sci. 2024, 14(9), 3949; https://doi.org/10.3390/app14093949

Submission received: 31 March 2024 / Revised: 1 May 2024 / Accepted: 3 May 2024 / Published: 6 May 2024

(This article belongs to the Section Civil Engineering)

Download

Browse Figures

Versions Notes

Abstract

This paper investigates the effect of spatial variability of soil elastic modulus on the longitudinal responses of the existing shield tunnel to the new tunnel undercrossing using a random two-stage analysis method (RTSAM). The Timoshenko–Winkler-based deterministic method considering longitudinal variation in the subgrade reaction coefficient and the random field of the soil elastic modulus discretized by the Karhunen–Loeve expansion method are combined to establish the RTSAM. Then, the proposed RTSAM is applied to carry out a random analysis based on an actual engineering case. Results show that the increases in the scale of fluctuation and the coefficient of variation of the soil elastic modulus lead to higher variabilities of tunnel responses. A decreasing pillar depth and mean value of the soil elastic modulus and an increasing skew angle strengthen the effect of the spatial variability of the soil elastic modulus on tunnel responses. The variabilities of tunnel responses under the random field of the soil elastic modulus are overestimated by the Euler–Bernoulli beam model. The results of this study provide references for the uncertainty analysis of the new tunneling-induced responses of the existing tunnel under the random field of soil properties.

Keywords:

new tunnel undercrossing; random tunnel response; random two-stage analysis method; spatial soil variability; random field; soil elastic modulus

1. Introduction

The construction of the underground transportation system has been gradually increasing in major cities. As the most advanced tunnel construction method at present, the shield method has become the first choice for urban tunnel construction in various regions of the world due to its advantages in low environmental impact, high construction efficiency, and high degree of mechanization. Therefore, a large number of existing shield tunnels have been constructed in urban underground space. Shield tunnels are generally assembled from reinforced concrete segments connected with steel bolts. Owing to their structural characteristics, existing shield tunnels are vulnerable to new tunnel construction underneath [1,2,3]. Because of the new tunnel undercrossing, longitudinal tunnel deformations may develop and threaten the safety and serviceability of the tunnel during long-term operation [4,5]. Therefore, the evaluation of tunneling-induced longitudinal deformations of existing shield tunnels has become a growing concern.

The longitudinal performance of a shield tunnel due to new tunneling can be studied by numerical models implemented in such commercial software as ABAQUS, FLAC, or PLAXIS, which are effective tools used to consider complex interaction between the tunnels and the ground [6,7,8]. However, the selection of constitutive parameters in the numerical simulation is difficult, and the modeling and computation generally cost substantial time. The complexity and inconvenience of the numerical model limit its application in engineering practice. By contrast, the one-dimensional two-stage analysis method is efficient and convenient for rapidly evaluating the effects of new tunnel excavation on longitudinal behaviors of shield tunnels in the primary design stage [9,10,11,12]. The existing tunnel in these analysis methods is generally assumed as a one-dimensional continuous beam lying on various foundation models. The greenfield soil settlement induced by new tunneling at the existing tunnel location is first determined through empirical equations or analytical solutions [13,14,15,16]. Then, the soil settlement can be used as an input and imposed to the tunnel for obtaining tunnel responses. Zhang and Huang [17] introduced a simplified two-stage analysis method based on the Euler–Bernoulli beam theory and the Winkler foundation model to analyze the responses of existing shield tunnels to undercrossing construction. The shearing deformation of shield tunnels is considered in further research by adopting the Timoshenko beam for better simulation of longitudinal tunnel behavior. Li et al. [18] derived a Timoshenko beam solution to calculate the deformations of existing tunnels owing to new tunneling underneath. Zhang et al. [19] proposed an analytical method based on Timoshenko beam theory and a Kerr foundation model for new tunnel construction underneath and perpendicular to the existing shield tunnels. Franza et al. [20] used a Timoshenko beam-based elastic-continuum model to assess the longitudinal responses of existing tunnels to single and twin tunneling.

The long-term geological deposition and weathering process result in apparent spatial variability of soil properties [8,21,22,23,24,25], which may lead to significant influences on the longitudinal tunnel responses. For the one-dimensional analysis method, the spatial soil variability is reflected by the variations in the soil parameters along the longitudinal direction of the tunnel [26,27,28]. However, there are limited studies available at present using the one-dimensional two-stage analysis method to investigate the new tunneling–existing tunnel interaction in spatially variable soil. Compared with other mechanical parameters of the soil, the soil elastic modulus has a relatively large spatial variability [29,30]. Therefore, the main objective of this paper is to perform a random analysis based on a proposed random two-stage analysis method (RTSAM) to explore the effect of the spatial variability of soil elastic modulus on the tunnel responses to new tunneling underneath.

This paper is structured in the following way. First, the establishment and verification of a deterministic two-stage analysis method for the new tunneling–existing tunnel interaction considering longitudinal variation in the subgrade reaction coefficient are presented. Next, the random field of the soil elastic modulus is introduced and discretized to develop the RTSAM. After that, the coupled effect of the variability of soil elastic modulus and the new tunnel undercrossing on the existing tunnel is studied through random analysis. Finally, the main conclusions are given.

2. Deterministic Analysis Method

The deterministic one-dimensional two-stage analysis method is shown in Figure 1 for assessing the longitudinal deformations of a shield tunnel due to new tunneling underneath. The existing shield tunnel is assumed as a continuous Timoshenko beam, and the soil around the tunnel is modeled as the Winkler foundation. Because of the spatial variability of soil elastic modulus, the subgrade reaction coefficient of the Winkler foundation varies along the longitudinal direction and cannot be assumed to be a constant. The relationship between the soil elastic modulus and the subgrade reaction coefficient can be described through the equation introduced by Yu et al. [31]:

k (x) = \frac{3.08}{ρ} \frac{E_{s} (x)}{(1 - ν^{2}) D_{e}} \sqrt[8]{\frac{E_{s} (x) D_{e}^{4}}{{(E I)}_{eq}}}

(1)

ρ = \{\begin{cases} 2.18 & (Z / D_{e} \leq 0.5) \\ (1 + \frac{1}{1.7 Z / D_{e}}) & (Z / D_{e} > 0.5) \end{cases}

(2)

where k(x) and E_s(x) are the subgrade reaction coefficient and the soil elastic modulus varying along the longitudinal tunnel axis, respectively; D_e and Z are the diameter and buried depth of the existing tunnel, respectively; ρ is the depth parameter; ν is the Poisson’s ratio of the soil; (EI)_eq is the equivalent longitudinal flexural stiffness of the existing tunnel. It is noted that the original equation in the study of Yu et al. [31] is modified here to make the unit of k(x) be kN/m³.

2.1. Greenfield Soil Settlement Caused by New Tunneling

The widely used semi-analytical solution developed by Loganathan and Poulos [16] is used to describe the greenfield soil settlement due to the new tunneling. The soil settlement U(x) at the axis depth of the existing tunnel induced by a volume loss of V_L is given in Equation (3).

\begin{array}{l} U (x) = R^{2} \{- \frac{Z - H}{x^{2} + {(Z - H)}^{2}} + \frac{(3 - 4 ν) (Z + H)}{x^{2} + {(Z + H)}^{2}} - \frac{2 Z [x^{2} - {(Z + H)}^{2}]}{{[x^{2} + {(Z + H)}^{2}]}^{2}}\} \\ V_{L} \exp \{- [\frac{1.38 x^{2}}{{(H + R)}^{2}} + \frac{0.69 Z^{2}}{H^{2}}]\} \end{array}

(3)

where R is the radius of the new tunnel and H is the buried depth of the new tunnel.

Equation (3) can be further modified to consider the skew angle between the new tunnel and the existing tunnel [32]:

\begin{array}{l} U (x) = R^{2} \{- \frac{Z - H}{{(x \sin a)}^{2} + {(Z - H)}^{2}} + \frac{(3 - 4 ν) (Z + H)}{{(x \sin a)}^{2} + {(Z + H)}^{2}} - \frac{2 Z [{(x \sin a)}^{2} - {(Z + H)}^{2}]}{{[{(x \sin a)}^{2} + {(Z + H)}^{2}]}^{2}}\} \\ V_{L} \exp \{- [\frac{1.38 {(x \sin a)}^{2}}{{(H + R)}^{2}} + \frac{0.69 Z^{2}}{H^{2}}]\} \end{array}

(4)

where a is the skew angle between the new and the existing tunnel.

2.2. Soil–Existing Tunnel Interaction

The greenfield tunneling-induced soil settlement is imposed on the existing tunnel, and the force analysis of an element of the tunnel is shown in Figure 2. The following equations can be derived according to the static equilibrium of the element:

Q - k (x) U D_{e} d x + k (x) w D_{e} d x = Q + d Q

(5)

\begin{array}{l} M + (Q + d Q) d x + k (x) U D_{e} \frac{{(d x)}^{2}}{2} \\ - k (x) w D_{e} \frac{{(d x)}^{2}}{2} = M + d M \end{array}

(6)

where M and Q are the bending moment and shear force of the existing tunnel, respectively; w is the settlement of the tunnel.

The higher-order terms of Equations (5) and (6) can be neglected, and the following equations are derived:

\frac{d Q}{d x} = - k (x) U D_{e} + k (x) w D_{e}

(7)

Q = \frac{d M}{d x}

(8)

The relationship between the displacements and the internal forces of the Timoshenko beam can be given as [33]

M = - {(E I)}_{eq} \frac{d θ}{d x}

(9)

Q = {(κ G A)}_{eq} (\frac{d w}{d x} - θ)

(10)

where θ denotes the shear angle and (κGA)_eq is the equivalent longitudinal shearing stiffness of the existing tunnel.

The relationship between the settlement and the shear angle is obtained by combining Equations (7)–(10):

\frac{{(E I)}_{eq}}{{(κ G A)}_{eq}} \frac{d^{2} θ}{d x^{2}} - θ = - \frac{d w}{d x}

(11)

\frac{d θ}{d x} = \frac{d^{2} w (x)}{d x^{2}} - \frac{k (x) w D_{e}}{{(κ G A)}_{eq}} + \frac{k (x) U D_{e}}{{(κ G A)}_{eq}}

(12)

Taking the first derivative of Equation (10) and the second derivative of Equation (12), the following equations can be derived:

\frac{{(E I)}_{eq}}{{(κ G A)}_{eq}} \frac{d^{3} θ}{d x^{3}} - \frac{d θ}{d x} = - \frac{d^{2} w}{d x^{2}}

(13)

\begin{array}{l} \frac{d^{3} θ}{d x^{3}} = \frac{d^{4} w (x)}{d x^{4}} - \frac{D_{e}}{{(κ G A)}_{eq}} (\frac{d^{2} k (x)}{d x^{2}} w + 2 \frac{d k (x)}{d x} \frac{d w}{d x} + k (x) \frac{d^{2} w}{d x^{2}}) \\ + \frac{D_{e}}{{(κ G A)}_{eq}} (\frac{d^{2} k (x)}{d x^{2}} U + 2 \frac{d k (x)}{d x} \frac{d U}{d x} + k (x) \frac{d^{2} U}{d x^{2}}) \end{array}

(14)

The governing differential equation for the tunnel settlement is given as Equation (15) by combining Equations (12)–(14). It is noted that when k(x) is assumed as a constant value, Equation (15) degenerates into the general differential equation of the Timoshenko–Winkler model used in existing studies [18,32] which ignore the longitudinal variation in the subgrade reaction coefficient.

\begin{array}{l} {(E I)}_{eq} \frac{d^{4} w}{d x^{4}} - \frac{{(E I)}_{eq} D_{e}}{{(κ G A)}_{eq}} (\frac{d^{2} k (x)}{d x^{2}} w + 2 \frac{d k (x)}{d x} \frac{d w}{d x} + k (x) \frac{d^{2} w}{d x^{2}}) + k (x) w D_{e} \\ = k (x) U D_{e} - \frac{{(E I)}_{eq} D_{e}}{{(κ G A)}_{eq}} (\frac{d^{2} k (x)}{d x^{2}} U + 2 \frac{d k (x)}{d x} \frac{d U}{d x} + k (x) \frac{d^{2} U}{d x^{2}}) \end{array}

(15)

Then, the bending moment and shear force can be expressed as follows:

M = - \frac{{(E I)}_{eq} D_{e}}{{(κ G A)}_{eq}} (\frac{{(κ G A)}_{eq}}{D_{e}} \frac{d^{2} w}{d x^{2}} + k (x) U - k (x) w)

(16)

Q = - \frac{{(E I)}_{eq} D_{e}}{{(κ G A)}_{eq}} (\frac{{(κ G A)}_{eq}}{D_{e}} \frac{d^{3} w}{d x^{3}} + \frac{d k (x)}{d x} U + k (x) \frac{d U}{d x} - \frac{d k (x)}{d x} w - k (x) \frac{d w}{d x})

(17)

Equation (15) is solved by adopting the finite difference method (FDM), and the finite discretization of the tunnel is illustrated in Figure 1. The tunnel is separated into n real segments and four virtual segments along the longitudinal direction. Therefore, the total number of discretization nodes is n + 5. Based on the finite difference theory, the central differential form for the governing differential equation can be given as follows:

\begin{array}{l} {(E I)}_{eq} \frac{6 w_{i} - 4 (w_{i + 1} + w_{i - 1}) + (w_{i + 2} + w_{i - 2})}{l^{4}} \\ - \frac{{(E I)}_{eq} D_{e}}{{(κ G A)}_{eq}} (\frac{(k_{i + 1} - 2 k_{i} + k_{i - 1}) w_{i}}{l^{2}} + \frac{(k_{i + 1} - k_{i - 1}) (w_{i + 1} - w_{i - 1})}{2 l^{2}} + \frac{k_{i} (w_{i + 1} - 2 w_{i} + w_{i - 1})}{l^{2}}) + k_{i} w_{i} D_{e} \\ = k_{i} U_{i} D_{e} - \frac{{(E I)}_{eq} D_{e}}{{(κ G A)}_{eq}} (\frac{(k_{i + 1} - 2 k_{i} + k_{i - 1}) U_{i}}{l^{2}} + \frac{(k_{i + 1} - k_{i - 1}) (U_{i + 1} - U_{i - 1})}{2 l^{2}} + \frac{k_{i} (U_{i + 1} - 2 U_{i} + U_{i - 1})}{l^{2}}) \end{array}

(18)

where w_i is the settlement of the tunnel at node i; k_i is the subgrade reaction coefficient at node i; U_i is the greenfield soil settlement at node i; l = S/n with S as the tunnel length, i = 0, 1, 2, 3, …, n.

The two ends of the tunnel are assumed to be unconstrained, and the boundary conditions are given as

M_{0} = - \frac{{(E I)}_{eq}}{{(κ G A)}_{eq}} ({(κ G A)}_{eq} \frac{w_{1} - 2 w_{0} + w_{- 1}}{l^{2}} + k_{0} U_{0} D_{e} - k_{0} D_{e} w_{0}) = 0

(19)

M_{n} = - \frac{{(E I)}_{eq}}{{(κ G A)}_{eq}} ({(κ G A)}_{eq} \frac{w_{n + 1} - 2 w_{n} + w_{n - 1}}{l^{2}} + k_{n} U_{n} D_{e} - k_{n} D_{e} w_{n}) = 0

(20)

Q_{0} = - \frac{{(E I)}_{eq} D_{e}}{{(κ G A)}_{eq}} (\begin{array}{l} \frac{{(κ G A)}_{eq}}{D_{e}} \frac{w_{2} - 2 w_{1} + 2 w_{- 1} - w_{- 2}}{2 l^{3}} + \frac{k_{1} - k_{- 1}}{2 l} U_{0} \\ + k_{0} \frac{U_{1} - U_{- 1}}{2 l} - \frac{k_{1} - k_{- 1}}{2 l} w_{0} - k_{0} \frac{w_{1} - w_{- 1}}{2 l} \end{array}) = 0

(21)

Q_{n} = - \frac{{(E I)}_{eq} D_{e}}{{(κ G A)}_{eq}} (\begin{array}{l} \frac{{(κ G A)}_{eq}}{D_{e}} \frac{w_{n + 2} - 2 w_{n + 1} + 2 w_{n - 1} - w_{n - 2}}{2 l^{3}} + \frac{k_{n + 1} - k_{n - 1}}{2 l} U_{n} \\ + k_{n} \frac{U_{n + 1} - U_{n - 1}}{2 l} - \frac{k_{n + 1} - k_{n - 1}}{2 l} w_{n} - k_{n} \frac{w_{n + 1} - w_{n - 1}}{2 l} \end{array}) = 0

(22)

The settlement of the virtual nodes located at both ends of the tunnel can be determined through the boundary conditions and expressed in the following way:

w_{- 1} = (\frac{k_{0} l^{2}}{{(κ G A)}_{eq}} + 2) w_{0} - w_{1} - \frac{k_{0} U_{0} D_{e} l^{2}}{{(κ G A)}_{eq}}

(23)

w_{n + 1} = (\frac{k_{n} l^{2}}{{(κ G A)}_{eq}} + 2) w_{n} - w_{n - 1} - \frac{k_{n} U_{n} D_{e} l^{2}}{{(κ G A)}_{eq}}

(24)

\begin{array}{l} w_{- 2} = (\frac{k_{0}^{2} D_{e}^{2} l^{4}}{{(κ G A)}_{eq}^{2}} + \frac{4 k_{0} D_{e} l^{2}}{{(κ G A)}_{eq}} + 4 - \frac{(k_{1} - k_{- 1}) D_{e} l^{2}}{{(κ G A)}_{eq}}) w_{0} \\ - (4 + \frac{2 k_{0} D_{e} l^{2}}{{(κ G A)}_{eq}}) w_{1} + w_{2} + \frac{k_{0} D_{e} l^{2} (U_{1} - U_{- 1})}{{(κ G A)}_{eq}} + \frac{U_{0} D_{e} l^{2} (k_{1} - k_{- 1})}{{(κ G A)}_{eq}} \\ - (\frac{2 D_{e} l^{2}}{{(κ G A)}_{eq}} + \frac{k_{0} D_{e}^{2} l^{4}}{{(κ G A)}_{eq}^{2}}) k_{0} U_{0} \end{array}

(25)

\begin{array}{l} w_{n + 2} = (\frac{k_{n}^{2} D_{e}^{2} l^{4}}{{(κ G A)}_{eq}^{2}} + \frac{4 k_{n} D_{e} l^{2}}{{(κ G A)}_{eq}} + 4 - \frac{(k_{n - 1} - k_{n + 1}) D_{e} l^{2}}{{(κ G A)}_{eq}}) w_{n} \\ - (4 + \frac{2 k_{n} D_{e} l^{2}}{{(κ G A)}_{eq}}) w_{n - 1} + w_{n - 2} + \frac{k_{0} D_{e} l^{2} (U_{n - 1} - U_{n + 1})}{{(κ G A)}_{eq}} + \frac{U_{0} D_{e} l^{2} (k_{n - 1} - k_{n + 1})}{{(κ G A)}_{eq}} \\ - (\frac{2 D_{e} l^{2}}{{(κ G A)}_{eq}} + \frac{k_{n} D_{e}^{2} l^{4}}{{(κ G A)}_{eq}^{2}}) k_{n} U_{n} \end{array}

(26)

The matrix form for the tunnel settlement can be obtained through combining the discretized formulations of each node:

w_{(n + 1) \times 1} = {(K_{1} - K_{2} + K_{3} - K_{4} - K_{5})}^{- 1} (f_{1} - f_{2} - f_{3} - f_{4} - f_{5})

(27)

where w_(n+1)×1 is the vector of the tunnel settlement (w_(n+1)×1 = {w₀, w₁, w₂, …, w_n₋₁, w_n}^T); K₁, K₂, K₃, K₄ and K₅ are the stiffness matrices; f₁, f₂, f₃, f₄, f₅ are the loading vectors for the greenfield soil settlement. The explicit expressions for the above vectors and matrices are presented in Appendix A. Moreover, it is noted that the calculation process in the above analysis method only involves algebraic operations so high efficiency can be guaranteed in the following random analysis.

2.3. Equivalent Longitudinal Flexural Stiffness and Shearing Stiffness

The equivalent longitudinal flexural stiffness and the equivalent shearing stiffness of a shield tunnel is given as follows [34,35,36]:

{(E I)}_{eq} = \frac{l_{s}}{l_{s} - λ l_{b} + λ l_{b} \frac{\cos ψ + (π / 2 + ψ) \sin ψ}{\cos^{3} ψ}} E_{c} I

(28)

ψ + \cot ψ = \frac{π}{2} + \frac{π n_{b} E_{b} A_{b}}{E_{c} A_{c}}

(29)

{(κ G A)}_{eq} = ξ l_{s} {(\frac{λ l_{b}}{n_{b} κ_{b} G_{b} A_{b}} + \frac{l_{s} - λ l_{b}}{κ_{c} G_{c} A_{c}})}^{- 1}

(30)

where l_s and l_b are the lengths of the shield segment and the longitudinal steel bolt, respectively; λ is the influencing factor for circumferential joints; ψ is the angle of neutral axis; I is the moment of inertia of the tunnel cross-section; n_b is the number of the longitudinal steel bolts; E_c and E_b are the elastic modulus of the concrete shield segment and the steel bolt, respectively; A_b and A_c are the cross-section area of the bolt and the tunnel, respectively; ξ is the modified factor of the equivalent shearing stiffness, κ_b and κ_c are the shear coefficient of the bolt and shield segment, and they are set to 0.9 and 0.5, respectively; G_b and G_c are the shear modulus of the bolt and segment, respectively.

2.4. Verification

The developed FDM-based deterministic analysis method is verified with the Timoshenko–Winkler-based closed-form analytical solution which is capable of evaluating the longitudinal responses of a tunnel lying on heterogeneous strata with two discontinuities [36]. A scenario of an existing tunnel under-crossed by a new tunnel is hypothesized for verification, as shown in Figure 3. The subgrade reaction coefficient k(x) in the scenario is given in Equation (31). The diameters of the new tunnel and the existing tunnel are both set to 6 m. The buried depths of the existing and new tunnels are assumed as 12 m and 21 m, respectively, and the skew angle between the existing and new tunnel is 90°. The volume loss caused by the new tunneling is set to 1%, and the Poisson’s ratio of the soil is 0.3. The equivalent longitudinal flexural and shearing stiffness of the existing tunnel are assumed as 1.2 × 10⁵ MN·m² and 2 × 10³ MN, respectively.

k (x) = \{\begin{cases} 10 \begin{matrix}  \end{matrix} (x \leq - 5 m) \\ 5 \begin{matrix}  \end{matrix} (- 5 m < x \leq 5 m) \\ 10 \begin{matrix}  \end{matrix} (x > 5 m) \end{cases} \begin{matrix}  \end{matrix} {MN / m}^{3}

(31)

Figure 4 compares the new tunneling-induced longitudinal settlement, bending moment, and shear force of the existing tunnel calculated by the analytical solution and the proposed analysis method. The results obtained from the proposed method agree well with those by the analytical solution, thereby verifying the accuracy of the method. It is noted that the tunnel responses affected by spatial variability of soil properties, as in the field, are more complex and cannot be evaluated with the analytical solution that is used here, considering only few discontinuities in the strata [26]. Therefore, the FDM-based analysis method proposed in Section 2.2 should be used in the following analysis.

3. Establishment of RTSAM

The establishment of the RTSAM in this study is achieved in two steps. First, the random field theory is utilized to describe the spatial variability of soil elastic modulus along the longitudinal direction of existing tunnel. Then, the random field of soil elastic modulus discretized by the Karhunen–Loeve (K-L) expansion method and the deterministic analysis method developed in Section 2 are combined to form the RTSAM and calculate the new tunneling-induced random responses of the existing tunnel.

3.1. Random Field of Soil Elastic Modulus

The spatial variability of soil properties can be reasonably described by the random field method in practical applications [21,23,29,37,38]. In this study, it is assumed that the soil elastic modulus in the longitudinal direction E_s(x) is characterized by a one-dimensional stationary lognormal random field with mean value μ_E, standard deviation σ_E, and scale of fluctuation (SOF) θ_E. The coefficient of variation (COV) of soil elastic modulus COV_E is expressed as σ_E/μ_E. The use of lognormal distribution aims to avoid the generation of unreal negative values of soil elastic modulus. The following equation can be given for the random field of the soil elastic modulus with lognormal distribution:

G_{E} (x) = \exp (μ_{\ln E} + σ_{\ln E} G (x))

(32)

where G_E(x) is the lognormal random field of the soil elastic modulus; G(x) is the correlated standard normal random field; μ_lnE and σ_lnE are the mean and standard deviation of the correlated random field with normal distribution, respectively, that can be given as

σ_{\ln E} = \sqrt{\ln (1 + \frac{σ_{E}^{2}}{μ_{E}^{2}})}

(33)

μ_{\ln E} = \ln (μ_{E}) - \frac{σ_{l n E}^{2}}{2}

(34)

A common exponential autocorrelation function is adopted to represent the correlation between the soil elastic modulus of two arbitrary points:

ρ_{E} = \exp (- \frac{2 |x_{1} - x_{2}|}{θ_{E}})

(35)

where ρ_E is the autocorrelation function; x₁ and x₂ are the longitudinal locations of the random field.

The lognormally distributed random field G_E(x) can be discretized through the K-L expansion method, that is, approximately expressed by a series including finite random variables [37,39,40] as shown in Equation (36). The K-L expansion method is implemented in this study because of its high efficiency and accuracy.

G_{E} (x) \approx \exp (μ_{\ln E} + σ_{\ln E} \sum_{i = 1}^{M_{k}} \sqrt{λ_{i}} ξ_{i} φ_{i} (x))

(36)

where ξ_i is a set of independent random variables following the standard normal distribution; φ_i(x) and λ_i are the eigenfunctions and eigenvalues of the autocorrelation function, respectively; M_k is the number of truncation terms of the K-L expansion that is usually obtained based on truncation error ε_r [37,41]:

ε_{r} = 1 - \frac{\sum_{i = 1}^{M_{k}} λ_{i}}{\sum_{i = 1}^{M_{a}} λ_{i}} < ε_{a}

(37)

where M_a is the total number of discretized nodes of the random field and ε_a is the allowable value of the truncation error.

3.2. Calculation Procedure of RTSAM

Figure 5 shows in detail the calculation procedure of the proposed RTSAM, which can be divided into eight steps:

(1) Determine the greenfield soil settlement caused by the new tunneling and the structural parameters of the existing tunnel.

(2) Use the mean value, coefficient of variation, scale of fluctuation, and autocorrelation function to characterize the random field of soil elastic modulus.

(3) Discrete the existing tunnel into finite difference segments along its longitudinal direction and obtain the position coordinates of the corresponding finite difference nodes.

(4) Determine the number of truncation terms of K-L expansion M_k and calculate truncation error ε_r combining with Equation (37) and the coordinates of the finite difference nodes, which are coincided with the discretized nodes of the random field.

(5) Determine whether truncation error ε_r is less than allowable error ε_a (ε_a = 0.1 in this study) to ensure the accuracy of the random field. If ε_r ≥ ε_a, return to Step (4), and the value of M_k should be increased. If ε_r < ε_a, the truncation error meets the requirement. Proceed to Step (6).

(6) Obtain a set of standard normal distribution samples of size M_k through the Monte Carlo simulation (MCS).

(7) Generate a realization of the random field of soil elastic modulus by utilizing the K-L series expansion method, that is, a set of randomly generated soil elastic modulus at the coordinates of the finite difference nodes.

(8) Use the discretized random field of soil elastic modulus and the finite difference matrix equation (Equation (27)) to obtain the random longitudinal response of the existing tunnel to new tunneling underneath.

4. Random Analysis

Random analysis is performed based on an actual engineering case to analyze the effect of the spatial variability of soil elastic modulus on the responses of the existing tunnel to the new tunnel undercrossing. Zhang et al. [19] reported the responses of an overlying metro tunnel to the excavation of a new tunnel in Shanghai, China. Both the new tunnel and the existing tunnel were constructed by the shield method and have an outer diameter of 6.2 m. The cross-sectional view of the new and existing tunnel is shown in Figure 6. The buried depths of the existing and new tunnels are 9.1 m and 20.1 m, respectively, and the skew angle between the existing and new tunnel is 90°. The vertical pillar depth P is 4.8 m. The parameters of segments and bolts for typical metro shield tunnels in Shanghai are listed in Table 1 referring to Wu et al. [35]. Joint parameters λ and ξ are assumed to be 0.4725 and 10, respectively [36,42]. Therefore, the equivalent longitudinal flexural and shearing stiffness of the existing tunnel are 1.36 × 10⁵ MN·m² and 4.2 × 10⁴ MN, respectively, which are close to the parameter values used in the study of Zhang et al. [19]. The volume loss caused by the new tunneling is set to 0.6%, and the Poisson’s ratio of the soil is 0.3. The mean value μ_E, the coefficient of variation COV_E, and the SOF θ_E of the soil elastic modulus are assumed to be 16.49 MPa and 0.5 and 18.6 m (3D), respectively [19].

Figure 7 illustrates the longitudinal settlement, bending moment, and shear force of the existing tunnel calculated from 500 random field realizations. The results from the deterministic analysis with a constant soil elastic modulus of 16.49 MPa are also included. As shown in Figure 7, the deterministic results are contained in the distribution ranges of the longitudinal tunnel responses on the basis of the random field. It can be seen that the random results at the extreme points of the response curves exhibit larger variability due to the spatial variability of the soil elastic modulus. The higher uncertainties in these extreme values may result in higher possibilities of maximum tunnel responses exceeding allowable limits. Furthermore, though the greenfield soil settlement due to new tunneling is symmetric with respect to the centerline of the new tunnel, the curves of the longitudinal tunnel responses may show asymmetric profiles because of the asymmetric longitudinal distribution of the soil elastic modulus. The asymmetry is more apparent for the curves of the internal force and cannot be reflected when the longitudinal variation in the soil elastic modulus is ignored. The settlements of the existing tunnel obtained from the measurement data are also shown in Figure 7a. It can be found that most measurement data fall within the random results, which proves the applicability of the proposed RTSAM in practical engineering to a certain extent. The results indicate that the proposed RTSAM can offer the random responses of the existing tunnel under the random field of the soil elastic modulus and provide a simplified methodology for the uncertainty analysis of the longitudinal tunnel responses to new tunneling underneath.

The effect of the random field of the soil elastic modulus on the maximum tunnel responses is explored in the following sections. The dimensionless parameters for the maximum tunnel responses are defined as follows:

(1): normalized maximum tunnel settlement w_m = w_max/U_max, where w_max and U_max are the maximum tunnel settlement and greenfield soil settlement, respectively;
(2): normalized maximum bending moment M_m = M_maxD²/((EI)_eqU_max), where M_max is the maximum bending moment of the existing tunnel;
(3): normalized maximum shear force Q_m = Q_maxD/((κGA)_eqU_max), where Q_max is the maximum shear force of the existing tunnel.

Figure 8 illustrates the evolution of calculated mean values and COVs of these dimensionless parameters with the increase in the number of MCS. As shown in Figure 8, the number of MCS can be set as 5000 in the following random analysis to obtain stable statistical results.

4.1. Effect of SOF θ_E

Figure 9 shows the changes in the mean values of the normalized maximum tunnel responses (w_m, M_m and Q_m) under different SOF θ_E values. It can be seen that the variation in θ_E has a slight influence on the mean values of the normalized responses. With the increase in θ_E from 0.5D to 5D, the mean values decrease slightly at a decreasing rate. The mean values of the normalized responses almost remain stable when θ_E is larger than 5D. The variations in the COVs of the normalized maximum tunnel responses with increasing θ_E are presented in Figure 10. As shown in Figure 10, the COVs of the normalized internal forces (M_m and Q_m) are always larger than those of the normalized tunnel settlement (w_m), indicating that the variability of the normalized internal forces is more apparent under the random field of the soil elastic modulus. The larger θ_E results in larger COVs of the normalized maximum tunnel responses. This phenomenon can be explained according to the logic of the study of Pan and Dias [37], which investigated the stability of the tunnel face in spatially random soils. When θ_E becomes smaller, the correlation between the soil elastic modulus of two arbitrary points becomes weaker, which makes the value of soil elastic modulus at different points under a single realization of the random field vary more greatly. The average soil elastic modulus within the influence range of the new tunneling changes less apparently among different random field realizations because of the decrease in θ_E, and the normalized maximum tunnel responses obtained through different realizations become closer to each other. When θ_E is relatively large, the average soil elastic modulus under different realizations varies significantly, which leads to an increase in the variability of the normalized maximum tunnel responses. The COVs of the normalized maximum tunnel responses increase remarkably when θ_E varies in the range from 0.5D to 10D. However, the variation of θ_E has a limited effect on the COVs of the normalized responses when θ_E is larger than 10D.

The empirical cumulative distribution functions (CDFs) of w_m, M_m, and Q_m under different θ_E (θ_E = 1D, 3D, 10D) are shown in Figure 11. The distribution ranges of w_m, M_m, and Q_m under θ_E = 1D are the narrowest, while the distribution ranges under θ_E = 10D are the widest. Therefore, the increasing θ_E leads to higher probabilities of the occurrence of large w_m, M_m, and Q_m. For example, with the increase in θ_E from 1D to 10D, the probability of M_m > 0.17 increases from 0.03 to 0.1. The findings presented here align with the results shown previously, i.e., the COVs of w_m, M_m, and Q_m under θ_E = 10D are the largest. Moreover, it is also seen that the CDF curves rotate clockwise about their intersection point as θ_E increases.

4.2. Effect of Coefficient of Variation COV_E

Figure 12 illustrates the changes in the mean values of the normalized maximum tunnel responses with increasing COV_E under different θ_E (θ_E = 1D, 3D, 10D). The effect of the COV of the soil elastic modulus on the mean values of the normalized responses is not very obvious. The increasing COV_E induces slight decreases in the mean values of w_m and M_m. Furthermore, it is observed that the increase in COV_E amplifies the effect of the variation in θ_E on the mean values of the normalized maximum tunnel responses. The variations in the COVs of the normalized maximum tunnel responses with increasing COV_E under different θ_E (θ_E = 1D, 3D, 10D) are shown in Figure 13. The increase in COV_E results in an approximately linear increase in the COVs of the normalized responses. It is worth noting that the maximum slope is observed at the maximum θ_E (10D), whereas the minimum slope can be seen at the minimum θ_E (1D). For example, the slope of variation in M_m with increasing COV_E under θ_E = 10D is about 1.7 times larger than that under θ_E = 1D. Therefore, the results imply that the increasing θ_E strengthens the effect of the variation in COV_E on the variability of the tunnel responses.

4.3. Effect of Mean Value μ_E

Figure 14 presents the variation in the mean values of normalized maximum tunnel responses with increasing μ_E under different θ_E (θ_E = 1D, 3D, 10D). The mean values of w_m, M_m, and Q_m all gradually increase with the increase in μ_E from 6.49 MPa to 46.49 MPa. According to the two-stage analysis method, the response of the existing tunnel is controlled by two components, i.e., the greenfield soil settlement caused by the new tunnel and the soil–existing tunnel interaction, which can be understood as the ability of the existing tunnel to resist the greenfield soil movement. Vorster et al. introduced the concept of the “tradeoff” between the greenfield measurements and the interaction analysis [43]. The relative importance of the greenfield soil settlement and the soil–structure interaction may change significantly under different construction or soil conditions. The increase in θ_E weakens the effect of the soil–tunnel interaction on the tunnel responses in the random field, and the existing tunnel tends to follow the greenfield settlement profile more and more under each single realization, which leads to the increases in the mean values of w_m, M_m, and Q_m.

Figure 15 illustrates the variations in the COVs of the normalized maximum tunnel responses with increasing θ_E under different θ_E (θ_E = 1D, 3D, 10D). The increases in μ_E result in reduction in the COVs of the normalized tunnel responses, indicating that the normalized maximum tunnel responses under different random field realizations become closer to each other. Therefore, it can be concluded that the effect of spatial variability of the soil elastic modulus becomes less obvious for the soil with higher values of μ_E. Moreover, the decrease in θ_E slightly reduces the slope of variation in the COVs of normalized tunnel responses, which means that the decreasing θ_E weakens the effect of the changes in μ_E on the variability of the normalized tunnel responses.

4.4. Effect of Pillar Depth P

The changes in the mean values of the normalized maximum tunnel responses with increasing P under different θ_E (θ_E = 1D, 3D, 10D) are shown in Figure 16. Previous case studies indicate that the pillar depth typically varies between 0D and 0.8D [5]. It is seen from Figure 16 that the mean values of w_m gradually increase, and the mean values of M_m and Q_m decrease when P increases from 0.2D to 1.8D. The profiles of the greenfield soil settlement and the tunnel settlement become flattened and wide due to the increase in P, which results in smaller internal forces of the existing tunnel. Figure 17 presents the variations in the COVs of the normalized maximum tunnel responses with increasing P under different θ_E (θ_E = 1D, 3D, 10D). As shown in Figure 17, the increase in P leads to an approximately linear decrease in the COV of the normalized responses. For example, the COV of w_m, M_m, and Q_m under θ_E =3D decreases by 48%, 37%, and 37%, respectively, due to the increase in P from 0.2D to 1.8D. Therefore, the increasing P weakens the effect of the spatial variability of the soil elastic modulus. It is also observed that the variation in θ_E has nearly no effect on the slope of variation in the COVs of the normalized tunnel responses.

4.5. Effect of Skew Angle a

Figure 18 presents the changes in the mean values of the normalized maximum tunnel responses with increasing a under different θ_E (θ_E = 1D, 3D, 10D). It is seen that as a increases from 30° to 90°, the mean values of w_m decrease, and the mean values of M_m and Q_m increase at a decreasing rate. The effect of a on the COVs of normalized maximum tunnel responses under different θ_E (θ_E = 1D, 3D, 10D) is presented in Figure 19. The decrease in a leads to an apparent negative influence on the COVs of w_m, M_m, and Q_m. For example, the COV of M_m decreases by 47% with the reduction in a from 90° to 30°. The decreasing a induces wider greenfield settlement profiles as seen in Equation (3). Therefore, it can be deduced that the effect of spatial variability of the soil elastic modulus becomes less remarkable when the profile of the greenfield soil settlement widens. Furthermore, the effect of variation in a on the variability of the normalized maximum tunnel responses becomes more prominent with increasing θ_E because the larger slope can be observed under larger θ_E.

4.6. Effect of Modified Factor ξ

The effect of the modified factor of equivalent shearing stiffness ξ on the mean values of the normalized maximum tunnel responses under θ_E = 3D is shown in Figure 20. For comparison, Figure 20 also includes the mean values calculated by the Euler–Bernoulli beam model (ξ = ∞). As ξ increases from 0.5 to 10, the mean value of w_m decreases by 8%, and the mean value of Q_m decreases by 93%. On the contrary, the mean value of M_m increases by 42% because of the increase in ξ. The results imply that the mean values of the normalized maximum tunnel responses vary at a decreasing rate and gradually come closer to the calculation results according to the Euler–Bernoulli beam model for larger values of ξ.

The COVs of the normalized maximum tunnel responses for different ξ values under θ_E = 3D are presented in Figure 21. It can be found that the COV of w_m, M_m, and Q_m increases by 13%, 31%, and 35%, respectively with the increase in ξ from 0.5 to 10, implying that the increasing equivalent shearing stiffness of the existing tunnel results in higher variabilities of the normalized tunnel responses. It can also be deduced that the existing tunnel tends to displace following the greenfield soil settlement and is barely affected by the spatial variability of the soil elastic modulus when the tunnel is infinitely flexible. The comparison of the calculation results from the Timoshenko beam model and the Euler–Bernoulli beam model indicates that the variabilities in w_m and M_m may be significantly overestimated by the Euler–Bernoulli beam model. Furthermore, the calculation results show that the COV of the normalized maximum tunnel responses can be fitted with the following hyperbolic equation:

{COV}_{n} = \frac{a_{h}}{1 + b_{h} / ξ}

(38)

where COV_n is the COV of the normalized maximum tunnel responses and a_h and b_h are two coefficients of the hyperbolic equation. The hyperbolic fitting curve shows good agreement with the calculated results and can be used to establish the mathematical relationship between COV_n and ξ. For w_m and M_m, the fitting coefficient ah corresponds to the value of COV_n under ξ = ∞ (0.0697 for w_m and 0.1489 for M_m) and is basically identical to the COV_n calculated by the Euler–Bernoulli beam model (0.0703 for w_m and 0.1481 for M_m). These findings further verify the effectiveness of the fitting curve and provide new understanding of the differences between the Timoshenko beam and the Euler–Bernoulli beam under the random field of the soil elastic modulus.

5. Conclusions

The effect of spatial variability of soil elastic modulus on shield tunnel responses to new tunneling underneath is investigated in this paper using the random two-stage analysis method. A Timoshenko–Winkler-based deterministic analysis method considering the variation in the subgrade reaction coefficient along the longitudinal direction is first established. The RTSAM is then developed by combining the deterministic analysis method and the discretized random field of the soil elastic modulus. Random analysis is performed to investigate the influence of calculation parameters on the responses of the existing tunnel under the random field. The following conclusions can be summarized:

(1): The spatial variability of the soil elastic modulus induces apparent variabilities of the longitudinal tunnel responses caused by new tunneling and may cause asymmetric tunnel responses. With the increase in θ_E, the mean values of the normalized maximum tunnel responses (w_m, M_m and Q_m) slightly decrease, but the COVs of w_m, M_m, and Q_m increase significantly. The variations in θ_E within 10D lead to a strong effect in the variation in the COVs of w_m, M_m, and Q_m. However, the effect of variations in θ_E becomes limited when θ_E is larger than 10D. The increasing θ_E causes higher probabilities of the occurrence of large w_m, M_m, and Q_m.
(2): The increase in COV_E barely affects the mean values of w_m, M_m, and Q_m, but it significantly increases the COVs of w_m, M_m, and Q_m. The larger the θ_E, the more obvious the effect of variations in COV_E on the COVs of w_m, M_m, and Q_m. Furthermore, the reductions in the pillar depth and the mean value of the soil elastic modulus and the increase in the skew angle lead to larger COVs of w_m, M_m, and Q_m, i.e., they amplify the effect of spatial variability of the soil elastic modulus on tunnel responses.
(3): There is a reduction in the mean values of w_m and Q_m and an increase in the mean values of M_m due to the increase in the equivalent shear stiffness of the existing tunnel. Higher variabilities in w_m, M_m, and Q_m are observed when the equivalent shearing stiffness increases, and the variation in COVs of w_m, M_m, and Q_m can be fitted well by the hyperbolic equation. The Euler–Bernoulli beam model overestimates the variabilities of the tunnel responses caused by spatial variability of the soil elastic modulus.

This study provides theoretical references for the uncertainty analysis of the longitudinal responses of the existing tunnel affected by the new tunneling underneath and the spatial variability of the soil elastic modulus. However, it should be noted that the nonlinear deformations of the soil and the tunnel structure are not considered in this study and need to be investigated in further research.

Author Contributions

Conceptualization, X.G. (Xiaolu Gan) and N.L.; methodology, X.G. (Xiaolu Gan) and N.L.; software, X.G. (Xiaolu Gan); validation, X.G. (Xiaolu Gan); formal analysis, X.G. (Xiaolu Gan) and N.L.; investigation, X.G. (Xiaolu Gan); resources, X.G. (Xiaonan Gong); data curation, X.G. (Xiaolu Gan); writing—original draft preparation, X.G. (Xiaolu Gan); writing—review and editing, A.B. and N.L.; visualization, X.G. (Xiaolu Gan); supervision, X.G. (Xiaonan Gong); project administration, N.L.; funding acquisition, X.G. (Xiaolu Gan). All authors have read and agreed to the published version of the manuscript.

Funding

This research was funded by Zhejiang Provincial Natural Science Foundation of China, grant number LQ24E080008, and China Scholarship Council, grant number 202006320252.

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

The data that support the findings of this study are available from the corresponding author upon reasonable request.

Conflicts of Interest

The authors declare no conflicts of interest.

Abbreviations

D	Diameter of new tunnel	n_b	Number of longitudinal steel bolt
R	Radius of new tunnel	E_c	Elastic modulus of shield segment
D_e	Diameter of existing tunnel	E_b	Elastic modulus of steel bolt
H	Buried depth of new tunnel	A_c	Cross-section area of shield segment
Z	Buried depth of existing tunnel	A_b	Cross-section area of steel bolt
a	Skew angle between new and existing tunnel	ξ	Modified factor of equivalent shear stiffness
V_L	Volume loss induced by new tunneling	κ_c	Shear coefficient of shield segment
U	Greenfield soil settlement	κ_b	Shear coefficient of bolt
U_max	Maximum greenfield soil settlement	G_c	Shear modulus of shield segment
k	Subgrade reaction coefficient	G_b	Shear modulus of steel bolt
ν	Poisson’s ratio of soil	μ_E	Mean value of soil elastic modulus
E_s	Soil elastic modulus	σ_E	Standard deviation of soil elastic modulus
ρ	Depth parameter	θ_E	Scale of fluctuation of soil elastic modulus
(EI)_eq	Equivalent longitudinal flexural stiffness of existing tunnel	COV_E	Coefficient of variation of soil elastic modulus
(κGA)_eq	Equivalent longitudinal shearing stiffness of existing tunnel	G_E	Lognormal random field of soil elastic modulus
M	Bending moment of existing tunnel	G	Correlated standard normal random field of soil elastic modulus
Q	Shear force of existing tunnel	μ_lnE	Mean value of correlated standard normal random field
w	Settlement of existing tunnel	σ_lnE	Standard deviation of correlated standard normal random field
θ	Shear angle of existing tunnel	ρ_k	Autocorrelation function
w_max	Maximum tunnel settlement	ξ_i	Independent standard normal distribution random variable
M_max	Maximum tunnel bending moment	λ_i	Eigenvalues of autocorrelation function
Q_max	Maximum tunnel shear force	φ_i	Eigenfunctions of autocorrelation function
w_m	Normalized maximum tunnel settlement	M_k	Number of truncation term of K-L expansion
M_m	Normalized maximum tunnel bending moment	ε_r	Truncation error
Q_m	Normalized maximum tunnel shear force	M_a	Total number of discretized nodes in random field
w	Vector of tunnel settlement	ε_a	Allowable truncation error
K₁, K₂, K₃, K₄, K₅	Stiffness matrices	COV_n	Coefficient of variation of normalized maximum tunnel response
f₁, f₂, f₃, f₄, f₅	Loading vector for greenfield soil settlement	a_h, b_h	Coefficients of hyperbolic equation
l_s	Length of shield segment
l_b	Length of longitudinal steel bolt
λ	Influencing factor for circumferential joint
ψ	Neutral axis angle
I	Moment of inertia of tunnel cross-section

Appendix A. Vectors and Matrices of Equation (27)

K_{1} = \frac{1}{l^{4}} {[\begin{matrix} A_{1} & A_{2} & 2 & 0 \\ A_{3} & 5 & - 4 & 1 \\ 1 & - 4 & 6 & - 4 & 1 \\ 1 & - 4 & 6 & - 4 & 1 \\ ⋮ & ⋮ & ⋮ & ⋮ & ⋮ \\ 1 & - 4 & 6 & - 4 & 1 \\ 1 & - 4 & 6 & - 4 & 1 \\ 1 & - 4 & 5 & B_{3} \\ 0 & 2 & B_{2} & B_{1} \end{matrix}]}_{(n + 1) \times (n + 1)}

(A1)

where,

A_{1} = \frac{k_{0}^{2} D_{e}^{2} l^{4}}{{(κ G A)}_{eq}^{2}} - \frac{(k_{1} - k_{- 1}) D_{e} l^{2}}{{(κ G A)}_{eq}} + 2

A_{2} = - \frac{2 k_{0} D_{e} l^{2}}{{(κ G A)}_{eq}} - 4

A_{3} = \frac{k_{0} D_{e} l^{2}}{{(κ G A)}_{eq}} - 2

B_{1} = \frac{k_{n}^{2} D_{e}^{2} l^{4}}{{(κ G A)}_{eq}^{2}} - \frac{(k_{n - 1} - k_{n + 1}) D_{e} l^{2}}{{(κ G A)}_{eq}} + 2

B_{2} = - \frac{2 k_{n} D_{e} l^{2}}{{(κ G A)}_{eq}} - 4

B_{3} = \frac{k_{n} D_{e} l^{2}}{{(κ G A)}_{eq}} - 2

K_{2} = \frac{k_{i} D_{e}}{{(κ G A)}_{eq} l^{2}} {[\begin{matrix} C_{1} & 0 \\ 1 & - 2 & 1 \\ 1 & - 2 & 1 \\ 1 & - 2 & 1 \\ ⋮ & ⋮ & ⋮ & ⋮ & ⋮ \\ 1 & - 2 & 1 \\ 1 & - 2 & 1 \\ 1 & - 2 & 1 \\ 0 & C_{2} \end{matrix}]}_{(n + 1) \times (n + 1)}

(A2)

where,

C_{1} = \frac{k_{0} D_{e} l^{2}}{{(κ G A)}_{eq}}

C_{2} = \frac{k_{n} D_{e} l^{2}}{{(κ G A)}_{eq}}

K_{3} = \frac{D_{e}}{{(E I)}_{eq}} {[\begin{matrix} k_{0} & 0 \\ k_{1} \\ k_{2} \\ ⋱ \\ k_{i} \\ ⋱ \\ k_{n - 2} \\ k_{n - 1} \\ 0 & k_{n} \end{matrix}]}_{(n + 1) \times (n + 1)}

(A3)

K_{4} = \frac{(k_{i + 1} - k_{i - 1}) D_{e}}{2 l^{2} {(κ G A)}_{eq}} {[\begin{matrix} C_{3} & 2 & 0 \\ - 1 & 0 & 1 \\ - 1 & 0 & 1 \\ ⋱ \\ - 1 & 0 & 1 \\ ⋱ \\ - 1 & 0 & 1 \\ - 1 & 0 & 1 \\ 0 & - 2 & C_{4} \end{matrix}]}_{(n + 1) \times (n + 1)}

(A4)

where,

C_{3} = - \frac{k_{0} D_{e} l^{2}}{{(κ G A)}_{eq}} - 2

C_{4} = \frac{k_{n} D_{e} l^{2}}{{(κ G A)}_{eq}} + 2

K_{5} = \frac{k_{i + 1} - 2 k_{i} + k_{i - 1}}{l^{2}} \frac{D_{e}}{{(κ G A)}_{eq}} {[\begin{matrix} 1 & 0 \\ 1 \\ 1 \\ ⋱ \\ 1 \\ ⋱ \\ 1 \\ 1 \\ 0 & 1 \end{matrix}]}_{(n + 1) \times (n + 1)}

(A5)

f_{1} = \frac{D_{e}}{{(E I)}_{eq}} {\{\begin{matrix} k_{0} U_{0} & k_{1} U_{1} & \dots & k_{n - 1} U_{n - 1} & k_{n} U_{n} \end{matrix}\}}_{_{(n + 1)}}^{T}

(A6)

f_{2} = \frac{D_{e}}{{(κ G A)}_{eq} l^{2}} {\{\begin{matrix} (k_{1} - 2 k_{0} + k_{- 1}) U_{0} \\ (k_{2} - 2 k_{1} + k_{0}) U_{1} \\ ⋮ \\ (k_{n} - 2 k_{n - 1} + k_{n - 2}) U_{n - 1} \\ (k_{n + 1} - 2 k_{n} + k_{n - 1}) U_{n} \end{matrix}\}}_{1 \times (n + 1)}

(A7)

f_{3} = \frac{D_{e}}{{(κ G A)}_{eq} l^{2}} {\{\begin{matrix} (U_{1} - 2 U_{0} + U_{- 1}) k_{0} \\ (U_{2} - 2 U_{1} + U_{0}) k_{1} \\ ⋮ \\ (U_{n} - 2 U_{n - 1} + U_{n - 2}) k_{n - 1} \\ (U_{n + 1} - 2 U_{n} + U_{n - 1}) k_{n} \end{matrix}\}}_{1 \times (n + 1)}

(A8)

f_{4} = \frac{D_{e}}{2 l^{2} {(κ G A)}_{eq}} {\{\begin{matrix} (k_{1} - k_{- 1}) (U_{1} - U_{- 1}) \\ (k_{2} - k_{0}) (U_{2} - U_{0}) \\ ⋮ \\ (k_{n} - k_{n - 2}) (U_{n} - U_{n - 2}) \\ (k_{n + 1} - k_{n - 1}) (U_{n + 1} - U_{n - 1}) \end{matrix}\}}_{1 \times (n + 1)}

(A9)

f_{5} = {\{\begin{matrix} D_{1} & D_{2} & 0 & \dots & 0 & D_{3} & D_{4} \end{matrix}\}}_{_{(n + 1)}}^{T}

(A10)

where,

D_{1} = \frac{2 k_{0} U_{0} D_{e}}{{(κ G A)}_{eq} l^{2}} + \frac{k_{0} D_{e} (U_{1} - U_{- 1})}{{(κ G A)}_{eq} l^{2}} + \frac{U_{0} D_{e} (k_{1} - k_{- 1})}{{(κ G A)}_{eq} l^{2}} - \frac{(k_{1} - k_{- 1}) k_{0} U_{0} D_{e}^{2}}{2 {(κ G A)}_{eq}^{2}}

D_{4} = \frac{2 k_{n} U_{n} D_{e}}{{(κ G A)}_{eq} l^{2}} + \frac{k_{n} D_{e} (U_{n - 1} - U_{n + 1})}{{(κ G A)}_{eq} l^{2}} + \frac{U_{n} D_{e} (k_{n - 1} - k_{n + 1})}{{(κ G A)}_{eq} l^{2}} + \frac{(k_{n + 1} - k_{n - 1}) k_{n} U_{n} D_{e}^{2}}{2 {(κ G A)}_{eq}^{2}}

D_{2} = - \frac{k_{0} U_{0} D_{e}}{{(κ G A)}_{eq} l^{2}}

D_{3} = - \frac{k_{n} U_{n} D_{e}}{{(κ G A)}_{eq} l^{2}}

References

Cooper, M.L.; Chapman, D.N.; Rogers, C.D.F.; Chan, A.H.C. Movements in the piccadilly Line tunnels due to the heathrow express construction. Geotechnique 2002, 52, 243–257. [Google Scholar] [CrossRef]
Chen, R.P.; Lin, X.T.; Kang, X.; Zhong, Z.Q.; Liu, Y.; Zhang, P.; Wu, H.N. Deformation and stress characteristics of existing twin tunnels induced by close-distance EPBS under-crossing. Tunn. Undergr. Space Technol. 2018, 82, 468–481. [Google Scholar] [CrossRef]
Gan, X.L.; Yu, J.L.; Gong, X.N.; Zhu, M. Characteristics and Countermeasures of Tunnel Heave due to Large-Diameter Shield Tunneling Underneath. J. Perform. Constr. Facil. 2020, 34, 04019081. [Google Scholar] [CrossRef]
Jin, D.L.; Yuan, D.J.; Li, X.; Zheng, H. An in-tunnel grouting protection method for excavating twin tunnels beneath an existing tunnel. Tunn. Undergr. Space Technol. 2018, 71, 27–35. [Google Scholar] [CrossRef]
Gan, X.L.; Yu, J.L.; Gong, X.N.; Liu, N.W.; Zheng, D.Z. Behaviours of existing shield tunnels due to tunnelling underneath considering asymmetric ground settlements. Undergr. Space 2022, 7, 882–897. [Google Scholar] [CrossRef]
Chakeri, H.; Hasanpour, R.; Hindistan, M.A. Bahtiyar ünver, Analysis of interaction between tunnels in soft ground by 3D numerical modeling. Bull. Eng. Geol. Environ. 2011, 70, 439–448. [Google Scholar] [CrossRef]
Avgerinos, V.; Potts, D.M.; Standing, J.R. Numerical investigation of the effects of tunnelling on existing tunnels. Géotechnique 2017, 67, 808–822. [Google Scholar] [CrossRef]
Chen, F.Y.; Wang, L.; Zhang, W.G. Reliability assessment on stability of tunnelling perpendicularly beneath an existing tunnel considering spatial variabilities of rock mass properties. Tunn. Undergr. Space Technol. 2019, 88, 276–289. [Google Scholar] [CrossRef]
Huang, M.S.; Chen, Z.R.; Li, Z. A simplified analysis method for the influence of tunneling on grouped piles. Tunn. Undergr. Space Technol. 2009, 24, 410–422. [Google Scholar] [CrossRef]
Liang, R.Z.; Xia, T.D.; Huang, M.S.; Lin, C.G. Simplified analytical method for evaluating the effects of adjacent excavation on shield tunnel considering the shearing effect. Comput. Geotech. 2017, 81, 167–187. [Google Scholar] [CrossRef]
Liang, R.Z.; Kang, C.; Xiang, L.M.; Li, Z.C.; Lin, C.G.; Gao, K.; Guo, Y. Responses of in-service shield tunnel to overcrossing tunnelling in soft ground. Environ. Earth. Sci. 2021, 80, 183. [Google Scholar] [CrossRef]
Ding, Z.; Zhang, M.B.; Zhang, X.; Wei, X.J. Theoretical analysis on the deformation of existing tunnel caused by under-crossing of large-diameter slurry shield considering construction factors. Tunn. Undergr. Space Technol. 2023, 133, 104913. [Google Scholar] [CrossRef]
Peck, R.B. Deep excavations and tunneling in soft ground. In Proceedings of the 7th International Conference on Soil Mechanics and Foundation Engineering, State of the Art Report, Mexico City, Mexico, 25–29 August 1969; Edinburgh University Library Catalogue: Edinburgh, Scotland, 1969; pp. 225–290. [Google Scholar]
Sagaseta, C. Analysis of undrained soil deformation due to ground loss. Géotechnique 1987, 37, 301–320. [Google Scholar] [CrossRef]
Verruijt, A.; Booker, J.R. Surface settlements due to deformation of a tunnel in an elastic half plane. Géotechnique 1996, 46, 753–756. [Google Scholar] [CrossRef]
Loganathan, N.; Poulos, H.G. Analytical Prediction for Tunneling-Induced Ground Movements in Clays. J. Geotech. Geoenviron. Eng. 1998, 124, 846–856. [Google Scholar] [CrossRef]
Zhang, Z.G.; Huang, M.S. Geotechnical influence on existing subway tunnels induced by multiline tunneling in Shanghai soft soil. Comput. Geotech. 2014, 56, 121–132. [Google Scholar] [CrossRef]
Li, P.; Du, S.J.; Shen, S.L.; Wang, Y.H.; Zhao, H.H. Timoshenko beam solution for the response of existing tunnels because of tunneling underneath. Int. J. Numer. Anal. Methods Geomech. 2016, 40, 766–784. [Google Scholar] [CrossRef]
Zhang, D.M.; Huang, Z.K.; Li, Z.L.; Zong, X.; Zhang, D.M. Analytical solution for the response of an existing tunnel to a new tunnel excavation underneath. Comput. Geotech. 2019, 108, 197–211. [Google Scholar] [CrossRef]
Franza, A.; Viggiani, G.M.B. Role of shear deformability on the response of tunnels and pipelines to single and twin tunneling. J. Geotechn. Geoenviron. Eng. ASCE. 2021, 147, 04021145. [Google Scholar] [CrossRef]
Vanmarcke, E.H. Probabilistic modeling of soil profiles. J. Geotech. Eng. Div. 1977, 103, 1227–1246. [Google Scholar] [CrossRef]
Phoon, K.K.; Kulhawy, F.H. Characterization of geotechnical variability. Can. Geotech. J. 1999, 36, 612–624. [Google Scholar] [CrossRef]
Mollon, G.; Phoon, K.K.; Dias, D. Validation of a new 2D failure mechanism for the stability analysis of a pressurized tunnel face in a spatially varying sand. J. Eng. Mech. ASCE 2011, 137, 8–21. [Google Scholar] [CrossRef]
Phoon, K.K.; Ching, J. Risk and Reliability in Geotechnical Engineering; CRC Press: Boca Raton, FL, USA, 2015. [Google Scholar]
Guo, X.F.; Du, D.C.; Dias, D. Reliability analysis of tunnel lining considering soil spatial variability. Eng. Struct. 2019, 196, 109332. [Google Scholar] [CrossRef]
Huang, H.W.; Gong, W.P.; Khoshnevisan, S.; Juang, C.H.; Zhang, D.M.; Wang, L. Simplified procedure for finite element analysis of the longitudinal performance of shield tunnels considering spatial soil variability in longitudinal direction. Comput. Geotech. 2015, 64, 132–145. [Google Scholar] [CrossRef]
Gong, W.; Juang, C.; Martin, J.; Tang, H.; Wang, Q.; Huang, H.W. Probabilistic analysis of tunnel longitudinal performance based upon conditional random field simulation of soil properties. Tunn. Undergr. Space Technol. 2018, 73, 1–14. [Google Scholar] [CrossRef]
Huang, W.M.; Wang, J.C.; Yang, Z.X.; Xu, R.Q. Analytical analysis of the longitudinal response of shield tunnel lining considering ring-to-ring interaction. Comput. Geotech. 2022, 146, 104705. [Google Scholar] [CrossRef]
Huang, H.W.; Xiao, L.; Zhang, D.M.; Zhang, J. Influence of spatial variability of soil Young’s modulus on tunnel convergence in soft soils. Eng. Geol. 2017, 228, 357–370. [Google Scholar] [CrossRef]
Zhang, J.Z.; Huang, H.W.; Zhang, D.M.; Zhou, M.L.; Tang, C.; Liu, D.J. Effect of ground surface surcharge on deformational performance of tunnel in spatially variable soil. Comput. Geotech. 2021, 136, 104229. [Google Scholar] [CrossRef]
Yu, J.; Zhang, C.R.; Huang, M.S. Soil-pipe interaction due to tunnelling: Assessment of Winkler modulus for underground pipelines. Comput. Geotech. 2013, 50, 17–28. [Google Scholar] [CrossRef]
Gan, X.L.; Yu, J.L.; Gong, X.N.; Liu, N.W.; Zhu, M. Probabilistic analysis for twin tunneling-induced longitudinal responses of existing shield tunnel, Tunn. Undergr. Space Technol. 2022, 120, 104317. [Google Scholar] [CrossRef]
Timoshenko, S.P. On the correction for shear of the differential equation for transverse vibration of prismatic bars. Philo. Mag. 1921, 41, 744. [Google Scholar] [CrossRef]
Wu, H.N.; Shen, S.L.; Liao, S.M.; Yin, Z.Y. Longitudinal structural modelling of shield tunnels considering shearing dislocation between segmental rings. Tunn. Undergr. Space Technol. 2015, 50, 317–323. [Google Scholar] [CrossRef]
Wu, H.N.; Shen, S.L.; Yang, J.; Zhou, A. Soil-tunnel interaction modelling for shield tunnels considering shearing dislocation in longitudinal joints. Tunn. Undergr. Space Technol. 2018, 78, 168–177. [Google Scholar] [CrossRef]
Cheng, H.Z.; Chen, R.P.; Wu, H.N.; Meng, F.Y.; Yi, Y.L. General solutions for the longitudinal deformation of shield tunnels with multiple discontinuities in strata. Tunn. Undergr. Space Technol. 2021, 107, 103652. [Google Scholar] [CrossRef]
Pan, Q.; Dias, D. Probabilistic evaluation of tunnel face stability in spatially random soils using sparse polynomial chaos expansion with global sensitivity analysis. Acta Geotech. 2017, 12, 1415–1429. [Google Scholar] [CrossRef]
Griffiths, D.V.; Fenton, G.A. Bearing capacity of spatially random soil: The undrained clay Prandtl problem revisited. Geotechnique 2001, 51, 351–359. [Google Scholar] [CrossRef]
Ghanem, R.G.; Spanos, P.D. Stochastic Finite Element—A Spectral Approach; Dover Publications: New York, NY, USA, 2003. [Google Scholar]
Cho, S.E. Probabilistic stability analysis of rainfall-induced landslides considering spatial variability of permeability. Eng. Geol. 2014, 171, 11–20. [Google Scholar] [CrossRef]
Jiang, S.H.; Li, D.Q.; Cao, Z.J.; Zhou, C.B.; Phoon, K.K. Efficient System Reliability Analysis of Slope Stability in Spatially Variable Soils Using Monte Carlo Simulation. J. Geotech. Geoenviron. Eng. 2016, 141, 04014096. [Google Scholar] [CrossRef]
Xu, L. Study on the Longitudinal Settlement of Shield Tunnel in Soft Soil; Tongji University: Shanghai, China, 2005. [Google Scholar]
Vorster TE, B.; Klar, A.; Soga, K.; Mair, R.J. Estimating the effects of tunneling on existing pipelines. J. Geotech. Geoenviron. Eng. 2005, 131, 1399–1410. [Google Scholar] [CrossRef]

Figure 1. Deterministic two-stage analysis method considering longitudinal variation in subgrade reaction coefficient: (a) cross-section view; (b) plan view.

Figure 2. Force analysis for an element of the existing tunnel.

Figure 3. Hypothetical scenario of an existing tunnel under-crossed by a new tunnel.

Figure 4. Comparison between tunnel responses obtained from proposed method and analytical solution: (a) Settlement, (b) Bending moment; (c) Shear force.

Figure 5. Calculation procedure of proposed RTSAM.

Figure 6. Relative position for the existing tunnel and the new tunnel in random analysis.

Figure 7. Longitudinal tunnel responses obtained from 500 realizations of the random field: (a) Settlement, (b) Bending moment; (c) Shear force.

Figure 8. The convergence of statistical results of normalized maximum tunnel responses: (a) Mean value, (b) Coefficient of variation.

Figure 9. Mean values of normalized maximum tunnel responses under different θ_E values: (a) Normalized maximum tunnel settlement, (b) Normalized maximum bending moment; (c) Normalized maximum shear force.

Figure 10. COVs of normalized maximum tunnel responses under different θ_E values.

Figure 11. CDFs of normalized maximum tunnel responses under different θ_E values: (a) Normalized maximum tunnel settlement, (b) Normalized maximum bending moment; (c) Normalized maximum shear force.

Figure 12. Effect of COV_E on mean values of normalized maximum tunnel responses under different θ_E values: (a) Normalized maximum tunnel settlement, (b) Normalized maximum bending moment; (c) Normalized maximum shear force.

Figure 13. Effect of COV_E on COVs of normalized maximum tunnel responses under different θ_E values: (a) Normalized maximum tunnel settlement, (b) Normalized maximum bending moment; (c) Normalized maximum shear force.

Figure 14. Effect of μ_E on mean values of normalized maximum tunnel responses under different θ_E: (a) Normalized maximum tunnel settlement, (b) Normalized maximum bending moment; (c) Normalized maximum shear force.

Figure 15. Effect of μ_E on COVs of normalized maximum tunnel responses under different θ_E values: (a) Normalized maximum tunnel settlement, (b) Normalized maximum bending moment; (c) Normalized maximum shear force.

Figure 16. Effect of P on mean values of normalized maximum tunnel responses under different θ_E values: (a) Normalized maximum tunnel settlement, (b) Normalized maximum bending moment; (c) Normalized maximum shear force.

Figure 17. Effect of P on COVs of normalized maximum tunnel responses under different θ_E values: (a) Normalized maximum tunnel settlement, (b) Normalized maximum bending moment; (c) Normalized maximum shear force.

Figure 18. Effect of a on mean values of normalized maximum tunnel responses under different θ_E values: (a) Normalized maximum tunnel settlement, (b) Normalized maximum bending moment; (c) Normalized maximum shear force.

Figure 19. Effect of a on COVs of normalized maximum tunnel responses under different θ_E values: (a) Normalized maximum tunnel settlement, (b) Normalized maximum bending moment; (c) Normalized maximum shear force.

Figure 20. Effect of ξ on mean values of normalized maximum tunnel responses: (a) Normalized maximum tunnel settlement, (b) Normalized maximum bending moment; (c) Normalized maximum shear force.

Figure 21. Effect of ξ on COVs of normalized maximum tunnel responses: (a) Normalized maximum tunnel settlement, (b) Normalized maximum bending moment; (c) Normalized maximum shear force.

Table 1. Typical structure parameters of shield tunnels used in the Shanghai metro line.

Parameter	Value
Segmental rings
Outer diameter (m)	6.2
Inner diameter (m)	5.5
Thickness (m)	0.35
Elastic modulus (MPa)	3.45 × 10⁴
Length (m)	1
Poisson’s ratio	0.2
Bolts
Longitudinal bolts number	17
Length (mm)	400
Diameter (mm)	30
Elastic modulus (MPa)	2.06 × 10⁵
Poisson’s ratio	0.3

Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

© 2024 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https://creativecommons.org/licenses/by/4.0/).

Share and Cite

MDPI and ACS Style

Gan, X.; Liu, N.; Bezuijen, A.; Gong, X. Random Responses of Shield Tunnel to New Tunnel Undercrossing Considering Spatial Variability of Soil Elastic Modulus. Appl. Sci. 2024, 14, 3949. https://doi.org/10.3390/app14093949

AMA Style

Gan X, Liu N, Bezuijen A, Gong X. Random Responses of Shield Tunnel to New Tunnel Undercrossing Considering Spatial Variability of Soil Elastic Modulus. Applied Sciences. 2024; 14(9):3949. https://doi.org/10.3390/app14093949

Chicago/Turabian Style

Gan, Xiaolu, Nianwu Liu, Adam Bezuijen, and Xiaonan Gong. 2024. "Random Responses of Shield Tunnel to New Tunnel Undercrossing Considering Spatial Variability of Soil Elastic Modulus" Applied Sciences 14, no. 9: 3949. https://doi.org/10.3390/app14093949

APA Style

Gan, X., Liu, N., Bezuijen, A., & Gong, X. (2024). Random Responses of Shield Tunnel to New Tunnel Undercrossing Considering Spatial Variability of Soil Elastic Modulus. Applied Sciences, 14(9), 3949. https://doi.org/10.3390/app14093949

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Menu

Random Responses of Shield Tunnel to New Tunnel Undercrossing Considering Spatial Variability of Soil Elastic Modulus

Abstract

1. Introduction

2. Deterministic Analysis Method

2.1. Greenfield Soil Settlement Caused by New Tunneling

2.2. Soil–Existing Tunnel Interaction

2.3. Equivalent Longitudinal Flexural Stiffness and Shearing Stiffness

2.4. Verification

3. Establishment of RTSAM

3.1. Random Field of Soil Elastic Modulus

3.2. Calculation Procedure of RTSAM

4. Random Analysis

4.1. Effect of SOF θ_E

4.2. Effect of Coefficient of Variation COV_E

4.3. Effect of Mean Value μ_E

4.4. Effect of Pillar Depth P

4.5. Effect of Skew Angle a

4.6. Effect of Modified Factor ξ

5. Conclusions

Author Contributions

Funding

Institutional Review Board Statement

Informed Consent Statement

Data Availability Statement

Conflicts of Interest

Abbreviations

Appendix A. Vectors and Matrices of Equation (27)

References

Share and Cite

Article Metrics

Article Access Statistics

Further Information

Guidelines

MDPI Initiatives

Follow MDPI

Article Menu

Random Responses of Shield Tunnel to New Tunnel Undercrossing Considering Spatial Variability of Soil Elastic Modulus

Abstract

1. Introduction

2. Deterministic Analysis Method

2.1. Greenfield Soil Settlement Caused by New Tunneling

2.2. Soil–Existing Tunnel Interaction

2.3. Equivalent Longitudinal Flexural Stiffness and Shearing Stiffness

2.4. Verification

3. Establishment of RTSAM

3.1. Random Field of Soil Elastic Modulus

3.2. Calculation Procedure of RTSAM

4. Random Analysis

4.1. Effect of SOF θE

4.2. Effect of Coefficient of Variation COVE

4.3. Effect of Mean Value μE

4.4. Effect of Pillar Depth P

4.5. Effect of Skew Angle a

4.6. Effect of Modified Factor ξ

5. Conclusions

Author Contributions

Funding

Institutional Review Board Statement

Informed Consent Statement

Data Availability Statement

Conflicts of Interest

Abbreviations

Appendix A. Vectors and Matrices of Equation (27)

References

Share and Cite

Article Metrics

Article Access Statistics

Further Information

Guidelines

MDPI Initiatives

Follow MDPI

4.1. Effect of SOF θ_E

4.2. Effect of Coefficient of Variation COV_E

4.3. Effect of Mean Value μ_E