A Variational Approach to Analyze the Settlement of Existing Tunnels Caused by Ground Surcharge

Abstract

This paper presents a variational approach to assess the settlement of the operational shield tunnels resulting from surface loading. The vertical additional force on the tunnel, induced by the surcharge, was computed using the Boussinesq solution. The structural behavior of the tunnel was modeled using the Timoshenko beam theory, which accounts for both bending and shear deformation mechanisms. Furthermore, the two-parameter Pasternak foundation model, which accounts for the continuity of foundation deformation, was used to model the interaction between the tunnel and the surrounding ground. A finite Fourier series was employed to approximate the vertical displacement and cross-sectional rotation angle of the tunnel. By conducting work and energy analyses, the energy balance equations for the tunnel and the soil were obtained. The governing equations were then formulated according to the minimum potential energy principle. The displacement and cross-sectional rotation angle of the tunnel were then expressed through the variational method. The accuracy of the proposed method was validated by comparison with in situ measurement data, confirming its effectiveness in predicting tunnel responses under a ground surcharge. Finally, a parametric study was conducted to evaluate the impact of various parameters on the settlement of the shield tunnel.

1. Introduction

As urban subway networks expand, the density of tunnels continues to increase, which makes them increasingly susceptible to ground surcharge effects [1,2,3]. In 2014, 16 reported incidents of ground surcharge occurred along the Shanghai subway system [4]. Ground surcharge exerts additional force on the existing tunnel, causing longitudinal deformation in shield tunnels. If this deformation exceeds safe limits, it can lead to serious issues such as segment cracking, joint leakage, and bolt failure, all of which significantly jeopardize the tunnel’s operational safety [5,6,7]. Therefore, it is essential to study the longitudinal deformation of shield tunnels caused by ground surcharge to ensure the safety and structural integrity of urban subway systems.

Various research methods have been developed to investigate the adverse impacts of ground surcharge on the settlement of shield tunnels, such as theoretical analysis [8,9,10,11], case studies [5,12], numerical simulations [13,14], and model experiments [15,16]. Among these, theoretical analysis stands out due to its conceptual clarity and ease of application, meaning it is frequently used for early-stage settlement prediction in shield tunnels affected by ground surcharge. Hou [17] simplified the shield tunnel in soft soil as an infinite-length elastic foundation beam and derived the analytical solution for the beam under a uniform distributed load. Jiang et al. [18] analyzed the mechanical behavior of shield tunnels under ground surcharge and obtained the stress distribution across the tunnel section. Wang et al. [11] developed a simplified analytical framework to address the transverse forces and structural deformation of shield tunnels when exposed to substantial surface loads. Dai et al. [19] applied the Euler–Bernoulli beam model on a Winkler foundation and derived the equilibrium equation for longitudinal deformation. Their findings, achieved through a finite differential methodology, revealed that this model tends to overestimate internal forces while underestimating deformation. Kang et al. [20] modeled the shield tunnel as a Timoshenko beam supported by a Winkler foundation and subsequently calculated its longitudinal deformation. Wang et al. [21] derived a difference solution for the tunnel under arbitrary distributed loads using a Timoshenko beam model, taking into account spatial soil variability and arbitrary elastic boundary conditions. Wu et al. [22] proposed an analytical solution based on the simplified Timoshenko beam model to analyze the dislocation between rings, which was used to assess shield tunnel behavior under external forces transferred from ground surcharge. Liang et al. [23] applied a two-stage analysis with a Timoshenko beam on a Pasternak foundation to derive the analytical solution for longitudinal displacement. Zheng et al. [24] presented a simplified solution to determine the longitudinal deformation of the tunnel under ground loading. The Pasternak foundation model accounts for interactions between Winkler foundation springs by incorporating a shear layer, overcoming the limitations of the Winkler model.

The variational approach, known for its clear concepts, simplicity in solutions, and ease of achieving high calculation accuracy, has been widely applied to geotechnical engineering problems in recent years [25,26,27,28,29]. Avramidis and Morfidis [28] analyzed the bending problem of a Timoshenko beam resting on a Kerr foundation. They derived the general governing equation using the variational method and solved it to obtain an analytical expression for beam deformation. Liu et al. [29] introduced a method to assess vertical pipeline settlement influenced by adjacent tunnel crossing construction utilizing the energy variational method, which offers a new perspective on pipeline settlement analysis. However, its application in the theoretical study of the settlement of the existing tunnel resulting from ground surcharge is still relatively rare.

This paper presents a variational approach to investigate the effect of ground surcharge on the settlement of an existing tunnel, addressing gaps in previous studies on this topic. The tunnel was modeled as a Timoshenko beam supported by a Pasternak foundation, which accommodates the continuous deformation of the foundation. Using the two-stage method, the Boussinesq solution was employed to compute the vertical additional force exerted on the existing tunnel due to the ground surcharge. Subsequently, the total potential energy equation was derived based on the relationship between work and energy. By applying the variational approach, an analytical method for forecasting the settlement of the tunnel was obtained. The calculated results were compared with in situ measurement data, confirming the validity of the proposed method. On this basis, the impact of factors such as the intersection angle between the load centerline and the axis of the existing tunnel, the burial depth, and the surface surcharge magnitude on the settlement were analyzed. The research results provide recommendations for determining a reasonable range of temporary ground surcharge.

2. Methods

2.1. Problem Statement

The relative positional relationship between the ground surcharge and the existing shield tunnel is illustrated in Figure 1. The dimensions of the ground surcharge are defined as follows: L for the length, B for the width, p for the ground surcharge load, h for the burial depth of the tunnel axis, D for the outer diameter of the tunnel, l for the tunnel length, and α for the intersection angle between the load centerline and the axis of the existing tunnel. To clearly describe the positional relationship between the existing tunnel and the ground surcharge, both a global coordinate system (XOY) and a local coordinate system (xoy) are established. The global coordinate system’s origin, denoted as O, is located at the center of the ground surcharge, where the X-axis is parallel to the longer side and the Y-axis is perpendicular to it. The local coordinate system’s origin, denoted as o, is positioned at a specific point within the tunnel, where the x-axis coinciding with the tunnel axis and the y-axis is perpendicular to it. The distance between O and o is represented as d.

2.2. Additional Force Induced by Surface Surcharge

As illustrated in Figure 2, the schematic model used for calculating the settlement of the tunnel subjected to ground surcharge is based on several key assumptions: (1) The foundation is considered an isotropic, semi-infinite elastic body [11,21,23]. (2) The effects of the tunnel itself are neglected when calculating the additional forces on the tunnel subjected to the ground surcharge. (3) The Timoshenko model is employed to analyze the structural response of the existing tunnel, accounting for both bending and shear deformations. (4) The Pasternak model is adopted to describe the interaction between the tunnel and the surrounding soil, considering both the shear deformation and compression of the foundation.

In analyzing the settlement of the existing tunnel subjected to ground surcharge, calculating the additional force exerted on the tunnel axis due to ground surcharge is a critical step. To achieve this, the Boussinesq formula was employed, which describes the vertical force distribution within a semi-infinite elastic body under uniformly distributed loads. The vertical additional force σ_z generated by the ground surcharge on the tunnel axis is expressed as follows:

$${\sigma }_z={\int }_{-{\displaystyle \frac{L}{2}}}^{{\displaystyle \frac{L}{2}}}{\int }_{-{\displaystyle \frac{B}{2}}}^{{\displaystyle \frac{B}{2}}}{\displaystyle \frac{3p{h}^3}{2\pi R^5}}d\xi d\eta$$

(1)

$$R=\sqrt{{\left(X-\xi \right)}^2+{\left(Y-\eta \right)}^2+h^2}$$

(2)

where (ξ, η) denote the coordinates within the ground surcharge in the global coordinate system, and (X, Y) denote the global coordinates of a point (x, y) situated on the tunnel.

The conversion expression for any point in the local coordinate system to the global coordinate system was derived and is shown below, and it is established through the geometric connection between the global and local coordinate systems:

$$\left\{\begin{array}{c}X=y\sin \alpha +x\cos \alpha +d\cos \beta \\ Y=y\cos \alpha -x\sin \alpha +d\sin \beta \end{array}\right.$$

(3)

2.3. The Total Potential Energy of the Model

The relative positional relationship between the ground surcharge and the existing shield tunnel is illustrated in Figure 1. When the load centerline is oblique to the tunnel axis, the displacement of the shield tunnel becomes asymmetrical. Therefore, when the displacement is expanded as a Fourier series, it can be decomposed into sine and cosine components.

$$\left\{\begin{array}{c}w\left(x\right)={\displaystyle \sum _{i=0}^{i=n}}{a}_i cos\left({\displaystyle \frac{n\pi x}{l}}\right)+{\displaystyle \sum _{i=n+1}^{i=2n}}{a}_i sin\left({\displaystyle \frac{n\pi x}{l}}\right)=\mathit{f}{\mathit{\alpha }}^{\mathrm{T}}\\ \phi \left(x\right)={\displaystyle \sum _{i=0}^{i=n}}{b}_i cos\left({\displaystyle \frac{n\pi x}{l}}\right)+{\displaystyle \sum _{i=n+1}^{i=2n}}{b}_i sin\left({\displaystyle \frac{n\pi x}{l}}\right)=\mathit{f}{\mathit{\beta }}^{\mathrm{T}}\end{array}\right.$$

(4)

$$\mathit{f}=\left(1,\cos {\displaystyle \frac{\pi x}{l}},\cos {\displaystyle \frac{2\pi x}{l}}\dots \cos {\displaystyle \frac{n\pi x}{l}},\sin {\displaystyle \frac{\pi x}{l}},\sin {\displaystyle \frac{2\pi x}{l}}\dots \sin {\displaystyle \frac{n\pi x}{l}}\right)$$

(5)

$$\mathit{\alpha }=\left({\alpha }_1,{\alpha }_2\dots {\alpha }_{2n+1}\right)$$

(6)

$$\mathit{\beta }=\left({\beta }_1,{\beta }_2\dots {\beta }_{2n+1}\right)$$

(7)

where $w\left(x\right)$ represents the vertical displacement of the tunnel and $\phi \left(x\right)$ represents the cross-sectional rotation angle. The matrices $\mathit{\alpha }$ and $\mathit{\beta }$ are the coefficients to be determined in the Fourier series, and n represents the order of expansion for the Fourier series.

The total potential energy of the model, caused by the ground surcharge, includes three main components. Given that the tunnel is subjected solely to the additional vertical force induced by surface loading and does not experience any bending moments, the work ($W_q)$ performed on the tunnel by external loads can be calculated using the following equations:

$$W_p={\displaystyle \frac{1}{2}}{\int }_{-l/2}^{l/2}w{\sigma }_z\mathrm{d}x={\mathit{K}}_p\mathit{\gamma }$$

(8)

$${\mathit{K}}_p={\displaystyle \frac{1}{2}}\left[\begin{array}{cc}\mathit{f}& \\ & \mathbf{0}\end{array}\right]\left[\begin{array}{c}{\mathit{\sigma }}_z\\ \mathbf{0}\end{array}\right]$$

(9)

$$\mathit{\gamma }=\left[\mathit{\alpha },\mathit{\beta }\right]$$

(10)

The existing tunnel is modeled using the Timoshenko beam theory; therefore, the deformation energy of the tunnel ${(W}_t)$ includes both bending deformation energy ($W_{tb})$ and shear deformation energy ($W_{ts})$.

$$W_{tb}={\displaystyle \frac{1}{2}}{\int }_{-l/2}^{l/2}{\left(EI\right)}_{eq}{\left({\displaystyle \frac{\mathrm{d}\phi }{\mathrm{d}x}}\right)}^2\mathrm{d}x={\mathit{\gamma }}^{\mathrm{T}}{\mathit{K}}_{tb}\mathit{\gamma }$$

(11)

$$W_{ts}={\displaystyle \frac{1}{2}}{\int }_{-l/2}^{l/2}{\left(GA\right)}_{eq}{\left({\displaystyle \frac{\mathrm{d}w}{\mathrm{d}x}}-\phi \right)}^2\mathrm{d}x={\mathit{\gamma }}^{\mathrm{T}}{\mathit{K}}_{ts}\mathit{\gamma }$$

(12)

$${\mathit{K}}_{tb}={\displaystyle \frac{1}{2}}{\left(EI\right)}_{eq}{\left[\begin{array}{cc}& \\ & {\displaystyle \frac{\mathrm{d}\mathit{f}}{\mathrm{d}x}}\end{array}\right]}^{\mathrm{T}}\left[\begin{array}{cc}& \\ & {\displaystyle \frac{\mathrm{d}\mathit{f}}{\mathrm{d}x}}\end{array}\right]$$

(13)

$${\mathit{K}}_{ts}={\displaystyle \frac{1}{2}}{\left(GA\right)}_{eq}{\left[\begin{array}{cc}{\displaystyle \frac{\mathrm{d}\mathit{f}}{\mathrm{d}\mathit{x}}}& -\mathit{f}\end{array}\right]}^{\mathrm{T}}\left[\begin{array}{cc}{\displaystyle \frac{\mathrm{d}\mathit{f}}{\mathrm{d}\mathit{x}}}& -\mathit{f}\end{array}\right]$$

(14)

$$W_t=W_{tb}+W_{ts}={\mathit{\gamma }}^{\mathrm{T}}{\mathit{K}}_t\mathit{\gamma }$$

(15)

$${\mathit{K}}_t={\mathit{K}}_{tb}+{\mathit{K}}_{ts}$$

(16)

where ${\left(EI\right)}_{eq}$ represents the equivalent longitudinal bending stiffness of the tunnel; ${\left(GA\right)}_{eq}$ denotes the equivalent longitudinal shear stiffness of the tunnel.

Assuming that the tunnel is in close contact with the surrounding foundation, the vertical deformation of the tunnel is equal to that of the foundation. The foundation is modeled using the Pasternak foundation model; consequently, the deformation energy of the foundation includes both compressive deformation energy $W_{fb}$ and shear deformation energy $W_{fs}$.

$$W_{fb}={\displaystyle \frac{1}{2}}{\int }_{-l}^{l}D{k}_f{w}^2\mathrm{d}x={\mathit{\gamma }}^{\mathrm{T}}{\mathit{K}}_{fb}\mathit{\gamma }$$

(17)

$$W_{fs}={\displaystyle \frac{1}{2}}{\int }_{-l}^{l}D{G}_f{\left({\displaystyle \frac{dw}{dx}}\right)}^2\mathrm{d}x={\mathit{\gamma }}^{\mathrm{T}}{\mathit{K}}_{fs}\mathit{\gamma }$$

(18)

$${\mathit{K}}_{fb}={\displaystyle \frac{1}{2}}D{k}_f{\left[\begin{array}{cc}\mathit{f}& \\ & \mathbf{0}\end{array}\right]}^T\left[\begin{array}{cc}\mathit{f}& \\ & \mathbf{0}\end{array}\right]$$

(19)

$${\mathit{K}}_{fs}={\displaystyle \frac{1}{2}}D{G}_f{\left[\begin{array}{cc}{\displaystyle \frac{d\mathit{f}}{d\mathit{x}}}& \\ & \mathbf{0}\end{array}\right]}^T\left[\begin{array}{cc}{\displaystyle \frac{d\mathit{f}}{d\mathit{x}}}& \\ & \mathbf{0}\end{array}\right]$$

(20)

$$W_f=W_{fb}+W_{fs}={\mathit{\gamma }}^{\mathrm{T}}{\mathit{K}}_f\mathit{\gamma }$$

(21)

$${\mathit{K}}_f={\mathit{K}}_{fb}+{\mathit{K}}_{fs}$$

(22)

where $k_f$ is the coefficient of the subgrade reaction and $G_f$ is the shear modulus of the foundation.

The total potential energy $\Pi$ of the model was thus calculated as follows:

$$\Pi =W_p+W_t+W_f$$

(23)

According to the principle of minimum potential energy, the undetermined coefficients were solved using the variational method.

$${\displaystyle \frac{\delta \Pi }{\delta \mathit{\gamma }}}={\displaystyle \frac{\delta W_p}{\delta \mathit{\gamma }}}+{\displaystyle \frac{\delta W_t}{\delta \mathit{\gamma }}}+{\displaystyle \frac{\delta W_f}{\delta \mathit{\gamma }}}=0$$

(24)

By substituting Equations (8), (16) and (21) into Equation (24), the following results were obtained:

$$\left({\mathit{K}}_t+{\mathit{K}}_f\right)\mathit{\gamma }={\mathit{K}}_p$$

(25)

Then, the undetermined coefficients were determined.

$$\mathit{\gamma }=\left({\mathit{K}}_t+{\mathit{K}}_f\right)\backslash {\mathit{K}}_p$$

(26)

The vertical displacement of the tunnel and cross-section rotation angle are obtained by substituting into Equation (4).

3. Recommendation for Parameters

3.1. The Coefficient of Subgrade Reaction

Given that the tunnel is buried at a certain depth within the stratum, Attewell et al. [30] recommend the following formula for estimating the coefficient of subgrade reaction:

$$k_f={\displaystyle \frac{1.3E_f}{\left(1-{\mu }^2\right)}}\sqrt[8]{{\displaystyle \frac{E_f{D}^4}{{\left(EI\right)}_{eq}}}}$$

(27)

where $E_f$ and v are the elastic modulus and Poisson’s ratio of soil, respectively.

This formula is derived under the assumption that the tunnel is buried at an infinite depth, which may lead to an overestimation of the soil stiffness [31]. However, in practical engineering applications, the tunnel is generally constructed at limited depths. To account for the influence of burial depth, Yu et al. [32] proposed a modified formula that incorporates a depth influence coefficient for tunnel burial:

$$k_f={\displaystyle \frac{3.08}{\eta }}{\displaystyle \frac{E_f}{\left(1-{\mu }^2\right)}}\sqrt[8]{{\displaystyle \frac{E_f{D}^4}{{\left(EI\right)}_{eq}}}}$$

(28)

where η represents the burial depth influence coefficient and is defined as follows:

$$\eta =\left\{\begin{array}{c}\hfill 2.18, {\displaystyle \frac{h}{D}}\le 0.5\\ \hfill 1+{\displaystyle \frac{D}{1.7h}}, {\displaystyle \frac{h}{D}}>0.5\end{array}\right.$$

(29)

The inclusion of the influence coefficient η facilitates a more accurate estimation of the subgrade reaction by considering the actual burial depth of the tunnel.

3.2. The Shear Modulus of the Foundation

According to Tanahashi [33], the shear modulus of the Pasternak foundation can be calculated using the following formula:

$$G_f={\displaystyle \frac{E_{f}t}{6\left(1+\mu \right)}}$$

(30)

where t is the shear layer thickness. According to Xu et al. [34], the shear layer thickness t can be assumed to be 2.5 times the tunnel diameter:

$$t=2.5D$$

(31)

3.3. The Equivalent Longitudinal Bending Stiffness

The shield tunnel is a barrel-shaped structure formed by connecting segments with bolts. Due to the presence of tunnel joints, its overall stiffness is lower than that of continuous concrete barrel structures. The longitudinal equivalent continuous model proposed by Shiba et al. [35] is a widely used method for determining the equivalent longitudinal bending stiffness of the shield tunnel. This model assumes that the tunnel is a uniform ring across the transverse direction. Furthermore, the shield tunnel is considered equivalent to a uniform continuous beam with identical stiffness, determined by the stiffness equivalence method. The equivalent longitudinal bending stiffness is expressed as follows:

$${\left(EI\right)}_{eq}={\displaystyle \frac{{\cos }^3\psi }{\cos \psi +\left(\psi +\frac{\pi }{2}\right)\sin \psi }}{E}_c{I}_c$$

(32)

where

$$\psi +\mathit{cot}\psi =\pi \left(0.5+{\displaystyle \frac{m{k}_b{l}_s}{E_c{A}_c}}\right)$$

(33)

$$k_b={\displaystyle \frac{E_b{A}_b}{l_b}}$$

(34)

In these formulas, $E_c$ is the Young’s modulus of the segment, $I_c$ is the moment of inertia of the tunnel’s cross-section, and $\psi$ is the position parameter of the central axis of the tunnel’s cross-section. The term $k_b$ is the average stiffness of the connecting bolt of the segment, where $E_b$ represents the Young’s modulus of the bolt, $A_b$ is the cross-sectional area of the bolt, and $l_b$ is the length of the bolt. m stands for the number of bolts, $l_s$ stands for the ring width, and $A_c$ stands for the cross-sectional area of the tunnel segment.

3.4. The Equivalent Longitudinal Shear Stiffness

The equivalent longitudinal shear stiffness of the shield tunnel was calculated according to the formula proposed by Wu [36].

$${\left(GA\right)}_{eq}=\zeta {\displaystyle \frac{l_{\mathrm{s}}}{\frac{l_{\mathrm{b}}}{n{\kappa }_{\mathrm{b}}{G}_{\mathrm{b}}{A}_{\mathrm{b}}}+\frac{l_{\mathrm{s}}-l_{\mathrm{b}}}{{\kappa }_{\mathrm{c}}{G}_{\mathrm{c}}{A}_{\mathrm{c}}}}}$$

(35)

where $\zeta$ is the correction coefficient, and in this study, its value is set to one. The variables ${\kappa }_{\mathrm{b}}$ and ${\kappa }_{\mathrm{c}}$ represent the Timoshenko shear coefficients for the bolts and segment rings, respectively. Specifically, ${\kappa }_{\mathrm{b}}$ is 0.9 for round section bolts, and ${\kappa }_{\mathrm{c}}$ is 0.5 for the ring tunnel segment structure. Additionally, $G_{\mathrm{b}}$ and $G_{\mathrm{c}}$ denote the shear stiffness of the bolt and the tunnel segment, respectively.

3.5. The Fourier Series Expansion Order n

In this paper, the order of the Fourier series expansion is a crucial parameter that affects the accuracy of the results. Under the premise of meeting the accuracy requirements, a reasonable order should be determined to reduce the computational burden as much as possible. To illustrate this, assume an engineering case where the load position and related parameters are shown in Figure 3. The equivalent flexural rigidity (EI)_eq is 7.8 × 10⁷ kN·m², the equivalent shear rigidity (GA)_eq is 2.1 × 10⁶ kN/m, the coefficient of subgrade reaction (k_f) is 2000 kN/m³, and the shear modulus of foundation $G_f$ is 1200 kN/m.

Figure 4 provides a comprehensive analysis of the convergence and accuracy of the proposed method. Figure 4a demonstrates the effects of different Fourier series expansion orders on the tunnel settlement at various points. As n increases (from n = 5 to n = 25), the settlement values gradually stabilize and converge with the results of the finite element method (FEM). When the Fourier series order reaches n = 15 or higher, the settlement values show almost no difference compared to the FEM results, thus verifying the accuracy of the proposed method. Figure 4b further examines the convergence of the method through the maximum settlement w_max; as the Fourier series expansion order increases, the maximum settlement stabilizes and reaches an asymptotic value when n ≥ 15. This indicates that once the Fourier series reaches a sufficient order, the method provides accurate and stable results. In conclusion, Figure 4a,b collectively demonstrate that the proposed Fourier series expansion method exhibits good convergence and high accuracy in the calculation of tunnel settlement. When the Fourier series expansion order reaches n = 15 or higher, the results are sufficiently precise and align with the FEM results. Considering the balance between the calculation accuracy and computational efficiency, the Fourier series expansion order is set to n = 20 for subsequent calculations.

4. Validations

In this section, the validity of the proposed method is established through a comparison with in situ measurement data from a ground surcharge case. The shield tunnel of the Zhongchun Road–Jiuting section of Shanghai Metro Line 9 was adversely affected by the backfill soil in the river course above, which led to the tunnel having significant uneven settlement. The relative positions of the tunnel and the backfill soil area are shown in Figure 5. The backfill soil area is 4.5 m high, 24 m wide, and 100 m long with a weight of 18 kN/m³. The distance between the load center and the tunnel is d = 0 m, and the tunnel axis is buried at a depth of 8.1 m. According to Fan et al. [7], the geological characteristics are presented in

Table 1. Physical and mechanical parameters of soils [7].

Soils	Thickness /(m)	Weight /(kN/m3)	Young’s Modulus /(MPa)	Poisson’s Ratio	c /(kPa)	$\mathit{\phi }$ /(°)
Backfill soil	4.5	17.0	-	-	-	-
②	3.0	18.7	15	0.31	22	18.5
③	11.0	17.2	9	0.33	13	15.5
④	9.0	16.7	6	0.33	14	10.5
⑥	6.0	19.4	21	0.30	45	15.5
⑦	14.0	19.0	36	0.29	3	31.5

, and the existing tunnel is situated within the third layer of soil. According to Wu et al. [22], the geometric and material parameters of the existing tunnel are shown in

Table 2. The geometric and material parameters of the shield tunnel [22].

Segmental Rings		Bolts
Outer diameter (mm)	6200	Number of longitudinal bolts	17
Thickness (mm)	350	Diameter (mm)	30
Length (mm)	1000	Length (mm)	400
Young’s modulus(kPa)	3.45 × 107	Young’s modulus (kPa)	2.06 × 108
Poisson’s ratio	0.2	Poisson’s ratio	0.3

.

Figure 6 shows the correlation between the vertical displacement w at different locations along the tunnel and the distance x from the tunnel center. The red line represents the method proposed in this paper, where the tunnel is modeled using the Timoshenko beam model and the foundation is modeled using the Pasternak foundation model; the blue line represents Liang’s predicted results, where the tunnel is modeled using the Euler beam model and the foundation is also modeled using the Pasternak foundation; and the black scatter points represent in situ measurement data. A comparison reveals that the maximum settlement calculated using this paper’s method is 28 mm, which is closer to the measured data, while the maximum settlement predicted using the Euler beam model is 24 mm, showing a larger deviation. The Timoshenko beam model considers not only the bending stiffness but also the deformation caused by shear forces, and it can more accurately predict the settlement characteristics of the shield tunnel under complex loading scenarios. In contrast, the Euler beam model only considers bending deformation and neglects the effects of shear forces, leading to smaller predicted vertical displacements and an inability to accurately reflect the tunnel’s actual stress and deformation under real working conditions. Thus, the Timoshenko beam model demonstrates significant advantages in handling both shear and bending deformations, particularly in the case of shield tunnels, which are large structures subjected to significant surcharge loading. The impact of shear deformation is particularly critical in such engineering projects. By taking shear deformation into account, the approach presented herein more accurately represents the actual settlement behavior of the tunnel. Furthermore, the high agreement between the predicted results and the in situ measurement data validates the method’s reliability in practical engineering applications.

5. Parametric Analysis

The relationship between the settlement of the shield tunnel and several parameters was analyzed based on the method proposed in this paper. These parameters include the interaction angle α between the load center line and the tunnel axis, the distance d between the load center O and the tunnel center o, the tunnel burial depth h, and the load q. A schematic diagram of the engineering case used for the parameter analysis is presented in Figure 3. The tunnel and foundation parameters are consistent with those described in Section 3.4.

5.1. The Influence of Load Position

Assuming that other parameters remain unchanged, the interaction angle ranges from 0° to 90°, and the distances are d = 0 m, 5 m, 10 m, and 15 m. Figure 7 illustrates the connection between the maximum settlement of the tunnel and both the interaction angle and the distance.

Overall, as the interaction angle increases, the maximum settlement tends to stabilize, indicating that a larger angle reduces the impact of the distance on tunnel settlement. When the distance is d = 0 m, the settlement significantly decreases with an increasing interaction angle, suggesting that when the load center is directly above the tunnel, changing the interaction angle can substantially reduce settlement. This trend is still noticeable at a smaller distance (e.g., (d = 5 m)), but as the distance increases (e.g., (d = 10 m) and (d = 15 m)), the effect diminishes. To mitigate the impact of the load on tunnel settlement, the load centerline should be positioned perpendicular to the tunnel axis when the load center is near the tunnel axis. Conversely, when the load center is positioned further from the tunnel axis, the load centerline should be aligned parallel to the tunnel axis.

5.2. The Influence of the Tunnel Burial Depth h

The burial depth varies from 5 m to 100 m, and the distances between the load center O and the tunnel center o are d = 0 m, 5 m, 10 m, 15 m, and 20 m. The relationship between the maximum settlement of the tunnel and burial depth and distance is shown in Figure 8. The maximum settlement gradually decreases as the burial depth increases for depths of 0 m, 5 m, and 10 m. A smaller burial depth results in a quicker reduction in the maximum settlement. This phenomenon occurs because as the tunnel approaches the load center, the additional force exerted by the load on the tunnel diminishes with the increased depth at which the tunnel is buried, leading to a decrease in the maximum settlement.

The maximum settlement initially increases and then gradually decreases with the increase in the burial depth for depths of 15 m and 20 m. The rates of both increases and decreases in the maximum settlement are higher when d = 15 m compared to d = 20 m. This behavior can be attributed to the fact that when the tunnel is situated far from the load center, the additional force induced by the load on the tunnel first rises and then declines with an increasing burial depth, resulting in the observed initial increase followed by a subsequent decrease in maximum tunnel settlement.

Figure 8 illustrates the trend in the maximum tunnel settlement as the burial depth increases while also examining the impact of different distances on the settlement. The burial depth ranges from 5 m to 100 m, and the distances between the load center O and the tunnel center o are d = 0 m, 5 m, 10 m, and 15 m. As the burial depth increases, the maximum settlement tends to stabilize. At smaller distances (e.g., d = 0 m and d = 5 m), the effect of the burial depth on the maximum settlement is more pronounced, showing a rapid increase in settlement at a lower burial depth followed by a tendency towards smaller negative values. This indicates that when the load center is closer to the tunnel center, an increase in the burial depth causes smaller settlement.

As the distance increases (for instance, d = 10 m and d = 15 m), the variation in the maximum settlement decreases and the settlement curves become more gradual. This signifies that when the point of load application is farther from the tunnel center, the impact of the burial depth on tunnel settlement is significantly reduced. This trend suggests that increasing the distance between the load and the tunnel center can effectively diminish settlement, thereby reducing potential risks to the tunnel structure.

5.3. The Influence of the Ground Surcharge Load q

The ground surcharge load varies from 10 kPa to 200 kPa, with distances set at 0 m, 5 m, 10 m, and 15 m, while keeping all other parameters constant. Figure 9 illustrates the correlation between the maximum tunnel settlement and both the ground surcharge load and the distance.

The maximum settlement decreases linearly with an increasing load regardless of the value of distance, indicating that a higher load results in greater settlement of the tunnel. This finding emphasizes the direct impact of the load on the tunnel structure. Simultaneously, as the distance increases from 0 m to 15 m, the maximum settlement under the same load decreases. This implies that when the load is applied farther from the tunnel center, the maximum settlement reduces. The greater the distance, the smaller the impact, suggesting that the effect of the load on tunnel settlement diminishes with an increasing distance. Thus, while an increase in load leads to greater tunnel settlement, increasing the distance from the load to the tunnel center can effectively reduce this effect. It is crucial to account for both the magnitude of the load and its position relative to the tunnel to minimize settlement risks.

6. Conclusions

This paper proposes a variational approach for assessing the settlement of the existing tunnel caused by ground surcharge to enhance existing analytical methods. The approach integrates the tunnel’s shear characteristics through the Timoshenko beam model and considers the continuous deformation of the foundation using the Pasternak foundation model. The Fourier series expansion order is set to n = 20 to balance both calculation accuracy and efficiency. Validation against in situ measurement data confirms the accuracy and reliability of the proposed method.

A detailed parametric analysis was carried out using the proposed method to assess the influence of various ground surcharge parameters on tunnel settlement. The maximum settlement of the tunnel is minimized when the interaction angle is α = 90° between the long side of the load and the tunnel axis and initially increases, then decreases, as the burial depth h increases, given that the distance d between the load center and the tunnel is less than the load width B. The maximum settlement of the tunnel is minimized when the angle is α = 0° and gradually decreases with the increase in the burial depth h, given that the distance is greater than the load width B. The tunnel’s maximum settlement increases linearly with the ground surcharge, and the rate of this increase decreases as the distance d increases.

It should be noted that the proposed method does not account for the nonlinear characteristics of the foundation reaction coefficient caused by foundation deformation. Further studies are needed to comprehensively evaluate the performance of the tunnel under ground surcharge.