Mechanical responses of multi-layered ground due to shallow tunneling with arbitrary ground surface load

Abstract

An analytical model based on complex variable theory is proposed to investigate ground responses due to shallow tunneling in multi-layered ground with an arbitrary ground surface load. The ground layers are assumed to be linear-elastic with full-stick contact between them. To solve the proposed multi-boundary problem, a series of analytic functions is introduced to accurately express the stresses and displacements contributed by different boundaries. Based on the principle of linear-elastic superposition, the multi-boundary problem is converted into a superposition of multiple single-boundary problems. The conformal mappings of different boundaries are independent of each other, which allows the stress and displacement fields to be obtained by the sum of components from each boundary. The analytical results are validated based on numerical and in situ monitoring results. The present model is superior to the classical model for analyzing ground responses of shallow tunneling in multi-layered ground; thus, it can be used with assurance to estimate the ground movement and surface building safety of shallow tunnel constructions beneath surface buildings. Moreover, the solution for the ground stress distribution can be used to estimate the safety of a single-layer composite ground.

Graphical abstract

Keywords

analytical model / mechanical response / multi-layered ground / shallow tunneling / ground surface load / complex variable solution

1 Introduction

Owing to rapid urbanization in recent years, shallow tunnel constructions in multi-layered ground have become increasingly common in many cities worldwide. During the construction of shallow urban tunnels, induced ground movements may cause adjacent surface or underground structures to deform or become damaged [1–4]. To ensure safety during shallow tunnel constructions, tunneling-induced ground responses have been extensively investigated, most of which have been performed for single-layer homogeneous ground. Previous studies have provided insights into the ground responses of shallow tunneling. However, owing to the differences in the mechanical properties of each ground layer, the ground responses of tunneling in multi-layered ground differ significantly from those in single-layer homogeneous ground [5–9]. Therefore, the ground responses of shallow tunneling in multi-layered ground must be elucidated to achieve desirable tunnel designs and adjacent structure protection.

To solve multi-layer problems, three main approaches are currently used: empirical, numerical, and analytical [10–14]. The empirical approach is characterized by a set of empirical equations from field monitoring results, such as the modified formulae based on the Peck formula [15] and the fitting formula presented by Selby [16]. However, the applicability of the empirical approach is ambiguous and should be evaluated closely when applied to complex and various ground conditions. The numerical approach, which considers the nonlinear behavior of the ground [17,18], is typically used to predict the tunneling-induced responses of multi-layered ground. However, the calculation is time-consuming when using commercial software. In particular, remodeling is typically required when performing parametric studies, which incurs a significant amount of time to obtain the results [19–21]. The analytical approach combines the first two approaches and can yield solutions more efficiently based on scientific theory [22–24]. In general, analytical solutions for the responses of a multi-layered ground are limited. Zhang et al. [25] presented a solution for the tunneling-induced displacement of multi-layered ground using a deformation-controlled boundary element method. Their model adopts the fundamental solutions for multi-layered ground to consider the stratification effect on ground movements. However, the calculation is relatively complicated and requires considerable effort to achieve satisfactory results. Zymnis et al. [26] presented ground displacement solutions for cross-anisotropic soils based on the superposition of fundamental singularity solutions for elastic materials. However, their model requires some empirical parameters, which are relatively difficult to determine for a specific project, thus limiting its application. Cao et al. [27] proposed a semi-analytical solution for multi-layered clayey soils based on the elastic equivalent theory with relatively simple formulae. However, the elastic equivalent method is based on an empirical formula that cannot consider the actual stratigraphic variations. When the positions of any ground layers are exchanged, the surface displacement remains unchanged, which is not consistent with reality. Moreover, all the abovementioned methods focus on displacement solutions without including tunneling-induced ground stress. Furthermore, owing to differences in the mechanical parameters of different ground layers, the stress distribution in those ground layers differs from that in homogenous ground. Unlike the stress distribution in homogenous ground, characterized as smoothly differentiable and gradually varying monotonically from the tunnel boundary, the stresses at the boundaries of multi-layered ground may exhibit significant fluctuations. However, the ground stress distribution is important for estimating the safety of a single-layer for the composite stratified ground.

The complex variable method, which is effective for analyzing the construction effects of tunnels with complex boundaries, has been widely used to solve shallow tunneling problems [28–33]. Through the conformal mapping technique, a semi-infinite medium with a shallow tunnel, characterized as a complex half-plane problem, can be converted into a simpler circular domain problem. Notably, the analytical solutions based on displacement control from the abovementioned studies provide important insight for this study, particularly the non-iterative analytical method for multiconnected domains proposed by Fang et al. [34]. However, the currently available complex variable solutions for the ground responses of shallow tunneling aimed at homogeneous ground conditions and the mechanical responses of multi-layered ground have not been solved. In fact, using the complex variable method to solve multi-layered ground responses would be faced with the challenge to satisfy all the boundary conditions of the ground boundaries, which cannot be solved using previous methods. Except for the Schwarz alternating method, most methods are unsuitable for handling multi-boundary conditions; however, the accuracy and efficiency of the alternating method are limited by the number of iterations, and a uniform conformal transformation does not exist for multiple boundaries. In this regard, all previous methods are inadequate for solving the proposed problem.

In addition, for shallow tunnels, the ground surface loading, such as the surface building weight, wind load on the building, crane load, and landfill load, significantly affects the ground response [35–40]. However, most previous analytical approaches rely significantly on a free surface or a specified loading pattern imposed on a homogeneous ground. In this regard, analytical solutions for the ground responses of multi-layered ground with arbitrary surface-loading effects have not been reported.

This study proposes a non-iterative analytical method based on the complex variable method to predict multi-layered ground responses due to shallow tunneling with arbitrary surface loading. First, the multi-boundary problem is transformed into multiple single-boundary problems that are independent of each other, by introducing a series of analytic functions for the stress and displacement of each boundary. Thereafter, the contributions of the tunnel and each ground layer are analyzed independently and superposed to obtain accurate solutions for the ground stress and displacement fields of multi-layered ground. Finally, the proposed method is verified based on comparisons with numerical simulations, previous analytical methods, and in situ monitoring results. The proposed model is particularly suitable for investigating multi-layered ground responses of shallow tunneling beneath surface buildings.

2 Problem definition

2.1 Analytical model and basic assumptions

A shallow circular tunnel is considered to be constructed in a multi-layered ground with ground surface loads, as shown in Fig.1. Each ground layer is assumed to be an elastic, homogeneous, and isotropic medium, characterized by an elastic modulus E_i, the Poisson’s ratio {\mu }_i, a unit weight {\gamma }_i, and an upper boundary depth h_i. Considering the ground surface structures, an arbitrarily distributed vertical load {\widetilde{\mathit{\text{σ}}}}_y (e.g., structure weight) and a horizontal load {\widetilde{\mathit{\text{τ}}}}_{xy} (e.g., wind load on buildings or crane load) exist on the ground surface. Regarding the contact conditions, two neighboring layers are assumed to be fully adhered, such that the ground displacement is continuous over the different layers. The tunnel is located within the tth layer with radius R and central buried depth h₀. A two-dimensional plane strain condition is adopted in the analysis, and the excavation effect of the tunnel is modeled by the double-circle convergence pattern suggested by Park [41] (see Fig.1).

The global coordinate system Z is established using the intersection of the geometrically symmetric vertical axis and the surface as the origin. Moreover, to facilitate the derivation, the geometric boundaries, including the layer and tunnel boundaries, are numbered based on layer number i. Next, B_{ij} is defined as the jth (j = 1,2,...,m_i) boundary of the ith layer with m_i boundaries, where B_{t{m}_t} denotes the tunnel boundary. For each boundary, a local coordinate system, Z_{ij}, is introduced. Notably, the interfaces that coincide in the physical plane belong to different computational domains and must be numbered separately. Hence, the transformation between the global and local coordinate systems can be expressed as

(1){\mathit{\text{z}}}_{ij}=\mathit{\text{z}}-v_{ij},

where z_ij is the local coordinate system, Z is the global coordinate system, and v_{ij} is the origin of the local coordinate system in the global coordinate system. The process for solving the problem can be summarized as follows.

(1) First, the original multi-boundary problem is converted into multiple single-boundary problems. Next, the conformal transformation for each independent boundary and the corresponding ground are identified to independently transform the boundary into a unit circle in the mapping plane with the ground transformed into inside/outside the circle. Subsequently, the relationship between the local physical coordinate z_ij and local mapping coordinate ζ_ij is established.

(2) Based on the superposition principle of the linear elasticity problem, the expressions of the ground stresses and displacements are determined using the analytic functions {\phi }_{ij}(ζ_ij) and ψ_ij(ζ_ij) in the mapping plane. Subsequently, the determination of the ground mechanical responses can be converted into the problem of solving unknown complex coefficients from the analytic functions.

(3) The calculated boundary points are selected based on the symmetric equalization of the argument in the reference local mapping coordinate system.

(4) All linear equations for the unknown complex coefficients at the aforementioned boundary points are established using the boundary conditions. Subsequently, the linear equations are solved to determine the aforementioned complex coefficients and to obtain the solution to the problem.

The derivation process of the problem is presented in detail below.

2.2 Conformal transformation for multi-layered ground

The conformal mapping function must be determined when performing a complex variable analysis on a specific boundary and the corresponding ground. The conformal transformations of the ground layer and tunnel boundaries differ owing to the geometry, as shown in Fig.2. The boundaries of the ground layers are transformed into a unit circle in the mapping plane, with the ground transformed within the circle. The mapping function of the ith ground layer is shown next. For the upper boundary and the ground underneath it,

(2){\mathit{\text{z}}}_{ij}=\mathit{\text{z}}+\text{i}{h}_i={\omega }_{ij}({\zeta }_{ij})=-\text{i}{\lambda }_{ij}\frac{1+{\zeta }_{ij}}{1-{\zeta }_{ij}},

where {\omega }_{ij}({\zeta }_{ij}) is the mapping function, and for the lower boundary and the ground above it,

(3){\mathit{\text{z}}}_{ij}=\mathit{\text{z}}+\text{i}{h}_{i+1}={\omega }_{ij}({\zeta }_{ij})=\text{i}{\lambda }_{ij}\frac{1+{\zeta }_{ij}}{1-{\zeta }_{ij}},

where {\zeta }_{ij} is the polar coordinate of {\mathit{\text{z}}}_{ij} in the mapping plane, expressed as {\zeta }_{ij}={\rho }_{ij}{\text{e}}^{\text{i}{\theta }_{ij}}; and {\lambda }_{ij} is the mapping scaling constant.

The tunnel boundary is transformed into a unit circle, with the outer ground transformed outside the circle. The mapping function is expressed as

(4){\mathit{\text{z}}}_{t{m}_t}=\mathit{\text{z}}+\mathrm{i}{h}_0={\omega }_{t{m}_t}({\zeta }_{t{m}_t})=R{\zeta }_{t{m}_t}.

2.3 Analytical expressions of ground responses

A set of single analytic functions in the conventional complex variable method expresses the stress and displacement fields. Because the boundary conditions cannot be satisfied simultaneously, the Schwarz alternating method is required to solve the multi-boundary problem. In this method, the multi-boundary problem is converted into multiple single-boundary problems by introducing a series of analytic functions for all boundaries, as shown in Fig.3. The contribution of each boundary to the mechanical responses in the computational domain is independent, and the single-boundary problem can be solved conveniently using complex function theory. According to the superposition principle of the linear elasticity problem, the mechanical response vector of any point in the computational domain is equal to the sum of the mechanical response vectors contributed by all the boundaries. In particular, in the Cartesian coordinate system, the vector sums of the aforementioned mechanical responses can be transformed into their algebraic sums in a uniform direction.

The non-iterative analytical process for multi-layered ground is illustrated based on a ground layer comprising a tunnel, as shown in Fig.3. The ground contains three boundaries, i.e., B_k, B_l1, and B_l2, which are in the global Cartesian coordinate system, Z. The local Cartesian coordinate systems of all boundaries are established, and the three-boundary problem is transformed into three independent single-boundary problems. In local Cartesian coordinate systems, complex function theory is used to solve the Cartesian coordinate component of the mechanical response for each boundary. For a point z₀, the mechanical responses are the algebraic sum of the contributions from the boundaries, B_k, B_l1, and B_l2. For point z₀ on the outer boundary, the mechanical responses toward the ground surface and tunnel opening should satisfy the corresponding boundary conditions. Similarly, when point z₀ is located on the inner boundary (the ground interface), the mechanical responses derived from two adjacent layers, except for the horizontal normal stress, should be equal at the interface. By satisfying these conditions, the quantitative relationships between the contributions from each boundary can be successfully established.

Based on the concept illustrated in Fig.3, a set of analytic functions Φ_ij(z_ij) and Ψ_ij(z_ij) is introduced to reformulate the contribution of the arbitrary boundary B_ij to the additional stress and displacement of a specific point in the ground. The analytic functions in the Cartesian coordinate system can be expressed as [34]

(5)\{\begin{array}{c}{\mathrm{\Phi }}_{ij}({\mathit{\text{z}}}_{ij})=\sum _{k=1}^{\mathrm{\infty }}{a}_{ij,k}^{\prime }{\mathit{\text{z}}}_{ij}^{-k}+a_{ij,0}^{\prime }\ln ({\mathit{\text{z}}}_{ij}),\\ {\mathrm{\Psi }}_{ij}({\mathit{\text{z}}}_{ij})=\sum _{k=1}^{\mathrm{\infty }}{b}_{ij,k}^{\prime }{\mathit{\text{z}}}_{ij}^{-k}-{\kappa }_i{\overline{a}}_{ij,0}^{\prime }\ln ({\mathit{\text{z}}}_{ij}),\end{array}

where a_{ij,0}^{\prime },a_{ij,k}^{\prime }, and b_{ij,k}^{\prime } are the undetermined complex coefficients; {\kappa }_i is correlated with the Poisson’s ratio, i.e., {\kappa }_i=3-4{\mu }_i; and the superscript “–” represents the conjugation of a complex number.

The stress components are correlated as follows:

(6)\{\begin{array}{c}{\sigma }_x+{\sigma }_y=4\mathrm{R}\mathrm{e}\{{\phi }^{\prime }(\mathit{\text{z}})\},\\ {\sigma }_y-{\sigma }_x+2\mathrm{i}{\tau }_{xy}=2(\overline{\mathit{\text{z}}}{\phi }^{″}(\mathit{\text{z}})+{\psi }^{\prime }(\mathit{\text{z}})),\end{array}

and the displacement components are correlated as follows:

(7)2G(u_x+\mathrm{i}{u}_y)=\kappa \phi (\mathit{\text{z}})-\mathit{\text{z}}\overline{{\phi }^{\prime }(\mathit{\text{z}})}-\overline{\psi (\mathit{\text{z}})}+A,

where G is the ground shear modulus, A is the complex constant associated with a rigid displacement, \phi (\mathit{\text{z}}) and \psi (\mathit{\text{z}}) are the analytic functions.

Therefore, the additional stress field contributed by the boundary B_ij can be expressed as

(8)\{\begin{array}{c}{\sigma }_{x,ij}^1(x,y)=\mathrm{R}\mathrm{e}\left\{2{\mathrm{\Phi }}_{ij}^{\prime }({\mathit{\text{z}}}_{ij})-{\overline{\mathit{\text{z}}}}_{ij}{\mathrm{\Phi }}_{ij}^{″}({\mathit{\text{z}}}_{ij})-{\mathrm{\Psi }}_{ij}^{\prime }({\mathit{\text{z}}}_{ij})\right\},\\ {\sigma }_{y,ij}^1(x,y)=\mathrm{R}\mathrm{e}\left\{2{\mathrm{\Phi }}_{ij}^{\prime }({\mathit{\text{z}}}_{ij})+{\overline{\mathit{\text{z}}}}_{ij}{\mathrm{\Phi }}_{ij}^{″}({\mathit{\text{z}}}_{ij})+{\mathrm{\Psi }}_{ij}^{\prime }({\mathit{\text{z}}}_{ij})\right\},\\ {\tau }_{xy,ij}^1(x,y)=\mathrm{I}\mathrm{m}\left\{{\overline{\mathit{\text{z}}}}_{ij}{\mathrm{\Phi }}_{ij}^{″}({\mathit{\text{z}}}_{ij})+{\mathrm{\Psi }}_{ij}^{\prime }({\mathit{\text{z}}}_{ij})\right\}.\end{array}

Similarly, the additional displacement field contributed by the boundary B_ij can be expressed as

(9)\{\begin{array}{c}{u}_{x,ij}^1(x,y)=\frac{1}{2G_i}\mathrm{R}\mathrm{e}\left\{{\kappa }_i{\mathrm{\Phi }}_{ij}({\mathit{\text{z}}}_{ij})-{\mathit{\text{z}}}_{ij}\overline{{\mathrm{\Phi }}_{ij}^{\prime }({\mathit{\text{z}}}_{ij})}-\overline{{\mathrm{\Psi }}_{ij}({\mathit{\text{z}}}_{ij})}+A_{ij}\right\},\\ u_{y,ij}^1(x,y)=\frac{1}{2G_i}\mathrm{I}\mathrm{m}\left\{{\kappa }_i{\mathrm{\Phi }}_{ij}({\mathit{\text{z}}}_{ij})-{\mathit{\text{z}}}_{ij}\overline{{\mathrm{\Phi }}_{ij}^{\prime }({\mathit{\text{z}}}_{ij})}-\overline{{\mathrm{\Psi }}_{ij}({\mathit{\text{z}}}_{ij})}+A_{ij}\right\},\end{array}

where G_i denotes the shear modulus of the ith ground layer and A_{ij} is a complex constant associated with the rigid displacement of boundary B_ij.

In tunnel engineering, the evaluation of the stress field should include the initial ground stress. In contrast, the initial displacement should be removed to illustrate the additional displacement field caused by the surface distribution load and tunnel excavation. The vertical and horizontal stresses are the main stresses when gravitational stress is considered. By defining the lateral pressure coefficient as {\lambda }_i={\mu }_i/\phantom{{\mu }_i\left(1-{\mu }_i\right)}\left(1-{\mu }_i\right), the initial ground stress components can be expressed as

(10)\{\begin{array}{c}{\sigma }_{x,i}^0(x,y)={\lambda }_i{\sigma }_{y,i}^0(x,y),\\ {\sigma }_{y,i}^0(x,y)=-\sum _{k=1}^{i-1}{\gamma }_k(h_k-h_{k-1})+{\gamma }_i(y+h_i),\\ {\tau }_{xy,i}^0(x,y)=0.\end{array}

Notably, the additional normal stresses, {\sigma }_{x,i}^1(x,y) and {\sigma }_{y,i}^1(x,y), as well as the shear stress {\tau }_{xy,i}^1(x,y) are the sum of the contributions of all the boundaries, which can be derived using Eq. (8) as follows:

(11)\{\begin{array}{rl}{\sigma }_{x,i}^1(x,y)& =\sum _{j=1}^{m_i}{\sigma }_{x,ij}^1(x,y)\\ & =\mathrm{R}\mathrm{e}\left\{\sum _{j=1}^{m_i}[2{\mathrm{\Phi }}_{ij}^{\prime }({\mathit{\text{z}}}_{ij})-{\overline{\mathit{\text{z}}}}_{ij}{\mathrm{\Phi }}_{ij}^{″}({\mathit{\text{z}}}_{ij})-{\mathrm{\Psi }}_{ij}^{\prime }({\mathit{\text{z}}}_{ij})]\right\},\\ {\sigma }_{y,i}^1(x,y)& =\sum _{j=1}^{m_i}{\sigma }_{y,ij}^1(x,y)\\ & =\mathrm{R}\mathrm{e}\left\{\sum _{j=1}^{m_i}[2{\mathrm{\Phi }}_{ij}^{\prime }({\mathit{\text{z}}}_{ij})+{\overline{\mathit{\text{z}}}}_{ij}{\mathrm{\Phi }}_{ij}^{″}({\mathit{\text{z}}}_{ij})+{\mathrm{\Psi }}_{ij}^{\prime }({\mathit{\text{z}}}_{ij})]\right\},\\ {\tau }_{xy,i}^1(x,y)& =\sum _{j=1}^{m_i}{\tau }_{xy,ij}^1(x,y)\\ & =\mathrm{I}\mathrm{m}\left\{\sum _{j=1}^{m_i}[{\overline{\mathit{\text{z}}}}_{ij}{\mathrm{\Phi }}_{ij}^{″}({\mathit{\text{z}}}_{ij})+{\mathrm{\Psi }}_{ij}^{\prime }({\mathit{\text{z}}}_{ij})]\right\}.\end{array}

Subsequently, the total stress field can be expressed by

(12)\{\begin{array}{c}{\sigma }_{x,i}(x,y)={\sigma }_{x,i}^0(x,y)+{\sigma }_{x,i}^1(x,y),\\ {\sigma }_{y,i}(x,y)={\sigma }_{y,i}^0(x,y)+{\sigma }_{y,i}^1(x,y),\\ {\tau }_{xy,i}(x,y)={\tau }_{xy,i}^0(x,y)+{\tau }_{xy,i}^1(x,y).\end{array}

Similarly, the displacement field due to tunnel excavation can be expressed as

(13)\{\begin{array}{rl}{u}_{x,i}^1(x,y)& =\sum _{j=1}^{m_i}{u}_{x,ij}^1(x,y)\\ & =\frac{1}{2G_i}\mathrm{R}\mathrm{e}\left\{\sum _{j=1}^{m_i}[{\kappa }_i{\mathrm{\Phi }}_{ij}({\mathit{\text{z}}}_{ij})-{\mathit{\text{z}}}_{ij}\overline{{\mathrm{\Phi }}_{ij}^{\prime }({\mathit{\text{z}}}_{ij})}-\overline{{\mathrm{\Psi }}_{ij}({\mathit{\text{z}}}_{ij})}]+A_i\right\},\\ u_{y,i}^1(x,y)& =\sum _{j=1}^{m_i}{u}_{y,ij}^1(x,y)\\ & =\frac{1}{2G_i}\mathrm{I}\mathrm{m}\left\{\sum _{j=1}^{m_i}[{\kappa }_i{\mathrm{\Phi }}_{ij}({\mathit{\text{z}}}_{ij})-{\mathit{\text{z}}}_{ij}\overline{{\mathrm{\Phi }}_{ij}^{\prime }({\mathit{\text{z}}}_{ij})}-\overline{{\mathrm{\Psi }}_{ij}({\mathit{\text{z}}}_{ij})}]+A_i\right\},\end{array}

where A_i={\displaystyle \sum _{j=1}^{m_i}{A}_{ij}}, which is the rigid displacement constant of the ith ground layer, should be determined by combining the continuous displacement condition of each layer and the vertical displacements of two symmetric points at the ground surface.

3 Determination of analytic functions

3.1 Transformation of analytic functions into mapping plane

For the ground layer boundaries (ij ≠ tm_t), the analytic functions Φ_ij(z_ij) and Ψ_ij(z_ij) are single-valued in the ith-layer computational domain, and the Taylor expansion can be performed directly in the mapping plane. Using the Taylor expansion in the unit circle obtained using Eq. (2) or (3), Eq. (5) for the ground boundary can be transformed into the mapping plane as follows:

(14)\{\begin{array}{c}{\mathrm{\Phi }}_{ij}({\mathit{\text{z}}}_{ij})={\mathrm{\Phi }}_{ij}({\zeta }_{ij})=\sum _{k=1}^{\mathrm{\infty }}{a}_{ij,k}{\zeta }_{ij}^k,\\ {\mathrm{\Psi }}_{ij}({\mathit{\text{z}}}_{ij})={\mathrm{\Psi }}_{ij}({\zeta }_{ij})=\sum _{k=1}^{\mathrm{\infty }}{b}_{ij,k}{\zeta }_{ij}^k,\end{array}

where a_{ij,k} and b_{ij,k} are the undetermined complex coefficients.

For the tunnel boundary (ij = tm_t), the logarithmic term \ln ({\mathit{\text{z}}}_{ij}) in {\mathrm{\Phi }}_{ij}({\mathit{\text{z}}}_{ij}) and {\mathrm{\Psi }}_{ij}({\mathit{\text{z}}}_{ij}) is a multiple-valued function. Therefore, the single-valued component of {\mathrm{\Phi }}_{ij}({\mathit{\text{z}}}_{ij}) and {\mathrm{\Psi }}_{ij}({\mathit{\text{z}}}_{ij}) can be obtained directly via Taylor expansion in the mapping domain, whereas the multivalued logarithmic term \ln ({\mathit{\text{z}}}_{ij}) must be considered separately via conformal transformation. Using Eq. (4), {\mathrm{\Phi }}_{ij}({\mathit{\text{z}}}_{ij}) and {\mathrm{\Psi }}_{ij}({\mathit{\text{z}}}_{ij}) for the tunnel boundary can be transformed into the mapping plane as follows:

(15)\{\begin{array}{c}{\mathrm{\Phi }}_{ij}({\mathit{\text{z}}}_{ij})={\mathrm{\Phi }}_{ij}({\zeta }_{ij})=\sum _{k=1}^{\mathrm{\infty }}{a}_{ij,k}{\zeta }_{ij}^{-k}+a_{ij,0}\ln ({\zeta }_{ij}),\\ {\mathrm{\Psi }}_{ij}({\mathit{\text{z}}}_{ij})={\mathrm{\Psi }}_{ij}({\zeta }_{ij})=\sum _{k=1}^{\mathrm{\infty }}{b}_{ij,k}{\zeta }_{ij}^{-k}-{\kappa }_t{\overline{a}}_{ij,0}\ln ({\zeta }_{ij}),\end{array}

where a_{ij,k}, b_{ij,k}, and a_{ij,0} are the undetermined complex coefficients; and {\kappa }_t=3-4{\mu }_t, where {\mu }_t is Poisson’s ratio of the ground in which the tunnel is located.

Using conformal transformation, the derivatives of {\mathrm{\Phi }}_{ij}({\mathit{\text{z}}}_{ij}) and {\mathrm{\Psi }}_{ij}({\mathit{\text{z}}}_{ij}) can be expressed in terms of the analytic functions {\phi }_{ij}({\zeta }_{ij}), {\psi }_{ij}({\zeta }_{ij}) and the mapping function {\omega }_{ij}({\zeta }_{ij}) as follows:

(16)\{\begin{array}{rl}{\mathrm{\Phi }}_{ij}^{\prime }({\mathit{\text{z}}}_{ij})=& \frac{\text{d}{\mathrm{\Phi }}_{ij}({\zeta }_{ij})}{\text{d}{\zeta }_{ij}}\frac{\text{d}{\zeta }_{ij}}{\text{d}{\omega }_{ij}({\zeta }_{ij})}=\frac{{\mathrm{\Phi }}_{ij}^{\prime }({\zeta }_{ij})}{{\omega }_{ij}^{\prime }({\zeta }_{ij})},\\ {\mathrm{\Psi }}_{ij}^{\prime }({\mathit{\text{z}}}_{ij})=& \frac{\text{d}{\mathrm{\Psi }}_{ij}({\zeta }_{ij})}{\text{d}{\zeta }_{ij}}\frac{\text{d}{\zeta }_{ij}}{\text{d}{\omega }_{ij}({\zeta }_{ij})}=\frac{{\mathrm{\Psi }}_{ij}^{\prime }({\zeta }_{ij})}{{\omega }_{ij}^{\prime }({\zeta }_{ij})},\\ {\mathrm{\Phi }}_{ij}^{″}({\mathit{\text{z}}}_{ij})=& \frac{\text{d[}{\mathrm{\Phi }}_{ij}^{\prime }({\zeta }_{ij})/{\omega }_{ij}^{\prime }({\zeta }_{ij})\text{]}}{\text{d}{\zeta }_{ij}}\frac{\text{d}{\zeta }_{ij}}{\text{d}{\omega }_{ij}({\zeta }_{ij})}\\ =& \frac{{\mathrm{\Phi }}_{ij}^{″}({\zeta }_{ij}){\omega }_{ij}^{\prime }({\zeta }_{ij})-{\mathrm{\Phi }}_{ij}^{\prime }({\zeta }_{ij}){\omega }_{ij}^{″}({\zeta }_{ij})}{{[{\omega }_{ij}^{\prime }({\zeta }_{ij})]}^3}.\end{array}

By substituting Eqs. (14)–(16) into Eqs. (11) and (13), the additional stress and displacement field due to excavation can be transformed into the mapping plane as follows:

(17)\{\begin{array}{c}{\sigma }_{x,i}^1(x,y)=\mathrm{R}\mathrm{e}\left\{\sum _{j=1}^{m_i}[{\mathrm{\Delta }}_{ij}({\zeta }_{ij})-T_{ij}({\zeta }_{ij})]\right\},\\ {\sigma }_{y,i}^1(x,y)=\mathrm{R}\mathrm{e}\left\{\sum _{j=1}^{m_i}[{\mathrm{\Gamma }}_{ij}({\zeta }_{ij})+T_{ij}({\zeta }_{ij})]\right\},\\ {\tau }_{xy,i}^1(x,y)=\mathrm{I}\mathrm{m}\left\{\sum _{j=1}^{m_i}[{\mathrm{\Theta }}_{ij}({\zeta }_{ij})+T_{ij}({\zeta }_{ij})]\right\},\end{array}

and

(18)\{\begin{array}{c}{u}_{x,i}^1(x,y)=\mathrm{R}\mathrm{e}\left\{\sum _{j=1}^{m_i}[{\mathrm{{\rm Y}}}_{ij}({\zeta }_{ij})-{\mathrm{\Xi }}_{ij}({\zeta }_{ij})]+A_i\right\},\\ u_{y,i}^1(x,y)=\mathrm{I}\mathrm{m}\left\{\sum _{j=1}^{m_i}[{\mathrm{{\rm Y}}}_{ij}({\zeta }_{ij})-{\mathrm{\Xi }}_{ij}({\zeta }_{ij})]+A_i\right\},\end{array}

where

(19)\{\begin{array}{c}\begin{array}{rl}{\mathrm{\Delta }}_{ij}({\zeta }_{ij})=& 2\frac{{\phi }_{ij}^{\prime }({\zeta }_{ij})}{{\omega }_{ij}^{\prime }({\zeta }_{ij})}-\frac{\overline{{\omega }_{ij}({\zeta }_{ij})}{\phi }_{ij}^{″}({\zeta }_{ij})}{{[{\omega }_{ij}^{\prime }({\zeta }_{ij})]}^2}\\ & +\frac{\overline{{\omega }_{ij}({\zeta }_{ij})}{\omega }_{ij}^{″}({\zeta }_{ij}){\phi }_{ij}^{\prime }({\zeta }_{ij})}{{[{\omega }_{ij}^{\prime }({\zeta }_{ij})]}^3},\end{array}\\ \begin{array}{rl}{\mathrm{\Gamma }}_{ij}({\zeta }_{ij})=& 2\frac{{\phi }_{ij}^{\prime }({\zeta }_{ij})}{{\omega }_{ij}^{\prime }({\zeta }_{ij})}+\frac{\overline{{\omega }_{ij}({\zeta }_{ij})}{\phi }_{ij}^{″}({\zeta }_{ij})}{{[{\omega }_{ij}^{\prime }({\zeta }_{ij})]}^2}\\ & -\frac{\overline{{\omega }_{ij}({\zeta }_{ij})}{\omega }_{ij}^{″}({\zeta }_{ij}){\phi }_{ij}^{\prime }({\zeta }_{ij})}{{[{\omega }_{ij}^{\prime }({\zeta }_{ij})]}^3},\end{array}\\ {\mathrm{\Theta }}_{ij}({\zeta }_{ij})=\overline{{\omega }_{ij}({\zeta }_{ij})}\frac{{\phi }_{ij}^{″}({\zeta }_{ij}){\omega }_{ij}^{\prime }({\zeta }_{ij})-{\phi }_{ij}^{\prime }({\zeta }_{ij}){\omega }_{ij}^{″}({\zeta }_{ij})}{{[{\omega }_{ij}^{\prime }({\zeta }_{ij})]}^3},\\ {\mathrm{{\rm Y}}}_{ij}({\zeta }_{ij})=\frac{1}{2G_i}[{\kappa }_i{\phi }_{ij}({\zeta }_{ij})-{\omega }_{ij}({\zeta }_{ij})\frac{\overline{{\phi }_{ij}^{\prime }({\zeta }_{ij})}}{\overline{{\omega }_{ij}^{\prime }({\zeta }_{ij})}}],\\ T_{ij}({\zeta }_{ij})=\frac{{\psi }_{ij}^{\prime }({\zeta }_{ij})}{{\omega }_{ij}^{\prime }({\zeta }_{ij})},{\mathrm{\Xi }}_{ij}({\zeta }_{ij})=\frac{\overline{{\psi }_{ij}({\zeta }_{ij})}}{2G_i}.\end{array}

3.2 Boundary conditions

3.2.1 Stress boundary conditions of ground surface B₁₁

Considering the surface loads {\widetilde{\mathit{\text{σ}}}}_y and {\widetilde{\mathit{\text{τ}}}}_{xy}, the stress boundary conditions of the ground surface can be written as

(20)\{\begin{array}{l}{\mathrm{R}\mathrm{e}\left\{{\displaystyle \sum _{j=1}^{m_1}\left[{\text{Γ}}_{1j}\left({\tilde{\zeta }}_{1j}\right)+T_{1j}\left({\tilde{\zeta }}_{1j}\right)\right]}\right\}|}_{\left|{\tilde{\xi }}_{11}\right|=1}={\tilde{\sigma }}_y\left({\tilde{\zeta }}_{11}\right),\\ {\mathrm{I}\mathrm{m}\left\{{\displaystyle \sum _{j=1}^{m_1}\left[{\text{Θ}}_{1j}\left({\tilde{\zeta }}_{1j}\right)+T_{1j}\left({\tilde{\zeta }}_{1j}\right)\right]}\right\}|}_{\left|{\tilde{\xi }}_{11}\right|=1}={\tilde{\tau }}_{xy}\left({\tilde{\zeta }}_{11}\right),\end{array}

where {\tilde{\zeta }}_{\text{1}j} corresponds to the local coordinate system {\zeta }_{1j} in the mapping plane, if j=1, \left|{\tilde{\zeta }}_{\text{1}j}\right|=1; and the superscript “\stackrel{~}{}” represents the positions located at the boundaries. Considering the mapping of z into two associated local coordinate systems {\zeta }_{ij} and {\zeta }_{kl} in the mapping plane, the relationship between {\tilde{\zeta }}_{ij} and {\tilde{\zeta }}_{kl} can be expressed as

(21)\mathit{\text{z}}={\omega }_{ij}\left({\tilde{\zeta }}_{ij}\right)+v_{ij}={\omega }_{kl}\left({\tilde{\zeta }}_{kl}\right)+v_{kl}.

The physical meaning of Eq. (21) is illustrated in Fig.2. For ground surface B₁₁, substituting i=1 into Eq. (21) and replacing {\tilde{\zeta }}_{kl} with {\tilde{\zeta }}_{\text{11}}, {\tilde{\zeta }}_{\text{1}j} can be obtained for a specified {\tilde{\zeta }}_{\text{11}}.

3.2.2 Stress and displacement continuity conditions at ground interfaces

Because the full-adherence condition is adopted between the ground layers, the vertical and horizontal displacements as well as the vertical normal and horizontal shear stresses derived from two adjacent layers, are equal at the interface. By considering the ground interface B_{(i + 1)l} between the ith and (i + 1)th layers, the stress and displacement continuity equations for all ground interfaces are established as follows: the stress continuity can be expressed as

(22)\{\begin{array}{rl}& \mathrm{R}\mathrm{e}\{\sum _{j=1}^{m_i}\left[{\mathrm{\Gamma }}_{ij}\left({\tilde{\zeta }}_{ij}\right)+T_{ij}\left({\tilde{\zeta }}_{ij}\right)\right]\\ -& {\sum _{j=1}^{m_{i+1}}\left[{\mathrm{\Gamma }}_{(i+1)j}\left({\tilde{\zeta }}_{(i+1)j}\right)+T_{(i+1)j}\left({\tilde{\zeta }}_{(i+1)j}\right)\right]\}|}_{\left|{\tilde{\zeta }}_{(i+1)1}\right|=1}=0,\\ & \mathrm{I}\mathrm{m}\{\sum _{j=1}^{m_i}\left[{\mathrm{\Theta }}_{ij}\left({\tilde{\zeta }}_{ij}\right)+T_{ij}\left({\tilde{\zeta }}_{ij}\right)\right]\\ -& {\sum _{j=1}^{m_{i+1}}\left[{\mathrm{\Theta }}_{(i+1)j}\left({\tilde{\zeta }}_{(i+1)j}\right)+T_{(i+1)j}\left({\tilde{\zeta }}_{(i+1)j}\right)\right]\}|}_{\left|{\tilde{\zeta }}_{(i+1)1}\right|=1}=0.\end{array}

The displacement continuity can be expressed as

(23)\{\begin{array}{rl}& \mathrm{R}\mathrm{e}\{\sum _{j=1}^{m_i}\left[{\mathrm{{\rm Y}}}_{ij}\left({\tilde{\zeta }}_{ij}\right)-{\mathrm{\Xi }}_{ij}\left({\tilde{\zeta }}_{ij}\right)\right]+A_i-A_{i+1}\\ -& {\sum _{j=1}^{m_{i+1}}\left[{\mathrm{{\rm Y}}}_{(i+1)j}\left({\tilde{\zeta }}_{(i+1)j}\right)-{\mathrm{\Xi }}_{(i+1)j}\left({\tilde{\zeta }}_{(i+1)j}\right)\right]\}|}_{\left|{\tilde{\zeta }}_{(i+1)1}\right|=1}=0,\\ & \mathrm{I}\mathrm{m}\{\sum _{j=1}^{m_i}\left[{\mathrm{{\rm Y}}}_{ij}\left({\tilde{\zeta }}_{ij}\right)-{\mathrm{\Xi }}_{ij}\left({\tilde{\zeta }}_{ij}\right)\right]+A_i-A_{i+1}\\ -& {\sum _{j=1}^{m_{i+1}}\left[{\mathrm{{\rm Y}}}_{(i+1)j}\left({\tilde{\zeta }}_{(i+1)j}\right)-{\mathrm{\Xi }}_{(i+1)j}\left({\tilde{\zeta }}_{(i+1)j}\right)\right]\}|}_{\left|{\tilde{\zeta }}_{(i+1)1}\right|=1}=0,\end{array}

where {\tilde{\zeta }}_{ij} corresponds to the local coordinate system {\zeta }_{ij} in the mapping plane and \left|{\tilde{\zeta }}_{(i+1)1}\right|\text{= 1}. When {\tilde{\zeta }}_{(i+1)1} is provided, {\tilde{\zeta }}_{ij} can be obtained using Eq. (21).

3.2.3 Displacement boundary condition at tunnel boundary B_{t{m}_t}

Because the double-circle model is used to characterize tunnel convergence, the displacement boundary equation at the tunnel boundary B_{t{m}_t} can be expressed as

(24)\{\begin{array}{rl}& {\mathrm{R}\mathrm{e}\left\{\sum _{j=1}^{m_t}\left[{\mathrm{{\rm Y}}}_{tj}\left({\tilde{\zeta }}_{tj}\right)-{\mathrm{\Xi }}_{tj}\left({\tilde{\zeta }}_{tj}\right)\right]+A_t\right\}|}_{{\tilde{\xi }}_{t{m}_t}=1}\\ =& {\tilde{u}}_x\left({\tilde{\zeta }}_{t{m}_t}\right)=\frac{Gap}{2}(1+\sin \alpha )\cos \alpha ,\\ & {\mathrm{I}\mathrm{m}\left\{\sum _{j=1}^{m_t}\left[{\mathrm{{\rm Y}}}_{tj}\left({\tilde{\zeta }}_{tj}\right)-{\mathrm{\Xi }}_{tj}\left({\tilde{\zeta }}_{tj}\right)\right]+A_t\right\}|}_{{\xi }_{t{m}_t}\mid =1}\\ =& {\tilde{u}}_y\left({\tilde{\zeta }}_{t{m}_t}\right)=\frac{Gap}{2}(1+\sin \alpha )\sin \alpha ,\end{array}

where the angle α denotes the position at the tunnel boundary; \sin \alpha =\left(y_t+h_0\right)/\phantom{\left(y_t+h_0\right)R}R and \cos \alpha =x_t/\phantom{x_{t}R}R; x_t+\text{i}{y}_t= {\mathit{\text{z}}}_t={\omega }_{t{m}_t}({\tilde{\zeta }}_{{tm}_t})-\text{i}{h}_0 defines the global coordinate value of any point at B_{t{m}_t}; \left|{\tilde{\zeta }}_{{tm}_t}\right|=1, {\tilde{\zeta }}_{tj} is the local coordinate of the tunnel boundary points in the mapping plane ζ_tj, which can be obtained based on a specified {\tilde{\zeta }}_{{tm}_t} using Eq. (21); and Gap is the gap parameter, which can be evaluated by referring to Lee et al. [42].

3.3 Solutions for complex coefficients

In the proposed method, the stress field of all the ground layers is reformulated using a series of analytic functions to satisfy the boundary conditions simultaneously without using the Schwarz alternative method. In this subsection, the general method of power series is used to establish the linear equations of the undetermined complex coefficients of all analytic functions for each ground layer to obtain the analytic functions directly.

3.3.1 Determination of all undetermined complex coefficient vectors

To avoid excessive calculations and arrive at accurate results, the highest and lowest powers of {\phi }_{ij}({\zeta }_{ij}) and {\psi }_{ij}({\zeta }_{ij}) corresponding to the boundary B_ij are assumed to be s_ij. Therefore, {\phi }_{ij}({\zeta }_{ij}) and {\psi }_{ij}({\zeta }_{ij}) in Eqs. (14) and (15) can be rewritten for ij\ne t{m}_t as follows:

\{\begin{array}{c}{\phi }_{ij}({\zeta }_{ij})=\sum _{k=1}^{s_{ij}}{a}_{ij,k}{\zeta }_{ij}^k,\\ {\psi }_{ij}({\zeta }_{ij})=\sum _{k=1}^{s_{ij}}{b}_{ij,k}{\zeta }_{ij}^k,\end{array}

and for ij=t{m}_t, these can be rewritten as

\{\begin{array}{c}{\phi }_{ij}({\zeta }_{ij})=\sum _{k=1}^{s_{ij}}{a}_{ij,k}{\zeta }_{ij}^{-k}+a_{ij,0}\ln ({\zeta }_{ij}),\\ {\psi }_{ij}({\zeta }_{ij})=\sum _{k=1}^{s_{ij}}{b}_{ij,k}{\zeta }_{ij}^{-k}-{\kappa }_t{\overline{a}}_{ij,0}\ln ({\zeta }_{ij}).\end{array}

The undetermined complex coefficient vectors in Eq. (25) are defined as Coef_ij, which is expressed in the following form:

(26){Coef}_{ij}=\{\begin{array}{c}{\left(\dots ,a_{ij,k},\dots ,b_{ij,k},\dots \right)}_{1\times 4s_{ij}},ij\ne t{m}_t,\\ {\left(a_{ij,0},\dots ,a_{ij,k},\dots ,b_{ij,k},\dots \right)}_{1\times (4s_{ij}+2)},ij=t{m}_t,\end{array}

where a_{ij,k}=\left(\text{Re}\{a_{ij,k}\},\text{Im}\{a_{ij,k}\}\right) and b_{ij,k}=\left(\mathrm{R}\mathrm{e}\{b_{ij,k}\},\text{Im}\{b_{ij,k}\}\right), with k = 1,2,...,s_ij.

Considering that each ground layer contains an undetermined complex constant related to the rigid body displacement, the undetermined complex coefficient vector set for the ith ground layer can be expressed as

(27)Coe{f}_i=\{\begin{array}{c}{\left(Coe{f}_{i1},\dots ,Coe{f}_{ij},\dots ,Coe{f}_{i{m}_i},A_i\right)}_{1\times \left(4\sum _{j=1}^{m_i}{s}_{ij}+2\right)},i\ne t,\\ {\left(Coe{f}_{i1},\dots ,Coe{f}_{ij},\dots ,Coe{f}_{i{m}_i},A_i\right)}_{1\times \left(4\sum _{j=1}^{m_i}{s}_{ij}+4\right)},i=t.\end{array}

Therefore, the entire set of undetermined complex coefficient vectors for the proposed model can be expressed as

(28)Coef={\left(Coe{f}_1,\dots ,Coe{f}_i,\dots ,Coe{f}_n\right)}_{1\times (4\sum _{i=1}^n\sum _{j=1}^{m_i}{s}_{ij}+2n+2)}.

Furthermore, the values of parameter s_ij are set as a constant s for the convenience of establishing linear equations.

3.3.2 Selection of boundary points for calculation

The boundary points should be widely distributed and concentrated in areas with significant mechanical responses. Furthermore, because the establishment of linear equations is based on the local mapping coordinate of each calculated boundary point, the images of the selected boundary points in the reference local coordinate system denoted as reference images must be distributed as evenly spaced as possible in the unit circle. Thus, the possible coincidence of the reference images can be avoided to ensure the independence of the linear equations. Therefore, the boundary points for the calculation are selected using the following method based on the symmetric equalization of the argument in the reference local mapping coordinate system.

(1) First, the reference boundary B_kl where the calculated boundary points are selected and the corresponding reference local mapping coordinate system ζ_kl are determined for a specific boundary. Three cases are established (Fig.4), as follows.

Case 1: The ground surface is selected as the reference boundary B₁₁, and the corresponding reference local mapping coordinate system is ζ₁₁.

Case 2: For the interface between the ith and (i + 1)th ground layers, the upper boundary of the (i + 1)th layer is selected as the reference boundary B_(i+1)1, and the corresponding reference local mapping coordinate system is ζ_(i+1)1.

Case 3: The tunnel boundary is selected as the reference boundary B_{t{m}_t}, and the corresponding reference local mapping coordinate system is {\zeta }_{t{m}_t}.

(2) A total of N = 2^m points are symmetrically selected anticlockwise on the unit circle mapped by B_kl, with m\in N^{+} and N>2s. The start angle of the reference images is set as {\theta }_0=\text{π}/\phantom{\text{π}N}N, and the angular interval is set as \mathrm{\Delta }\theta =2\text{π}/\phantom{2\text{π}N}N. Subsequently, the reference local mapping coordinate vectors of the calculated boundary points (reference mapping vector) are obtained as {\tilde{\xi }}_{kl}=e^{l{e}_2}, with {\mathit{\text{θ}}}_{kl}={\displaystyle \frac{\text{π}}{N}}+{\displaystyle \frac{2\text{π}}{N}}{\left(0,1,\dots ,N-1\right)}_{1\times N}^T.

(3) After determining {\tilde{\xi }}_{kl}, the correlated local mapping coordinate vector {\tilde{\xi }}_{ij} can be obtained using Eq. (21). To ensure that 50% of the calculated boundary points are located within the tunnel radius, the mapping scaling constant λ_ij in Eq. (3) is unified as R.

3.3.3 Linear equations for ground layers

To simplify the derivation, the dot product of two complex matrices with the same order, i.e., A = a_ij and B = b_ij (A, B \in C^m×n), is defined as the Hadamard product \mathit{C}=\mathit{A}\circ \mathit{B}, with C = c_ij = a_ij∙b_ij, and their dot division is defined as \mathit{D}={\displaystyle \frac{\mathit{A}}{\mathit{B}}}, with D = d_ij = a_ij/b_ij. Moreover, O_m×n and P_m×n refer to m × n matrices, whose elements are 0 and 1, respectively. For clarity and completeness, the general solutions for the additional stresses and displacements based on the boundary conditions are provided in Appendix A and Appendix B (Electronic Supplementary Materials), respectively.

For the ground surface, by considering the stress boundary expressed in Eq. (20), Eq. (A.1) can be rewritten as follows:

(29)\{\begin{array}{rl}& \mathrm{Re}\left\{\sum _{j=1}^{m_1}\sum _{k=1}^s\left[a_{1j,k}{\mathit{f}}_{1j,k}^{B_{11},{\sigma }_y}\left({\tilde{\mathit{\zeta }}}_{1j}\right)+b_{1j,k}{\mathit{g}}_{1j,k}^{B_{11},{\sigma }_y}\left({\tilde{\mathit{\zeta }}}_{1j}\right)\right]\right\}={\tilde{\mathit{\sigma }}}_y\left({\tilde{\mathit{\zeta }}}_{11}\right),\\ & \mathrm{Im}\left\{\sum _{j=1}^{m_1}\sum _{k=1}^s\left[a_{1j,k}{\mathit{f}}_{1j,k}^{B_{11},{\tau }_{yy}}\left({\tilde{\mathit{\zeta }}}_{1j}\right)+b_{1j,k}{\mathit{g}}_{1j,k}^{B_{11},{\tau }_y}\left({\tilde{\mathit{\zeta }}}_{1j}\right)\right]\right\}={\tilde{\mathit{\tau }}}_{xy}\left({\tilde{\mathit{\zeta }}}_{11}\right).\end{array}

Substituting Eq. (A.2) into Eq. (29) by considering Eq. (A.3) yields

(30){\mathit{K}}_1^{B_{11},S}{\mathit{C}\mathit{o}\mathit{e}\mathit{f}}_1^{\mathrm{T}}={\left(\begin{array}{c}{\tilde{\mathit{\sigma }}}_y\left({\tilde{\mathit{\zeta }}}_{11}\right)\\ {\tilde{\mathit{\tau }}}_{xy}\left({\tilde{\mathit{\zeta }}}_{11}\right)\end{array}\right)}_{2N\times 1}.

For the ground interfaces, by considering the stress continuity condition in Eq. (22), Eq. (A.1) can be rewritten as (for i\ne t-1\mathrm{o}\mathrm{r}t)

(31)\{\begin{array}{l}\begin{array}{rl}\mathrm{Re}& \{\sum _{j=1}^{m_i}\sum _{k=1}^s\left[a_{ij,k}{\mathit{f}}_{ij,k}^{R_{i(i+1)},{\sigma }_y}\left({\tilde{\zeta }}_{ij}\right)+b_{ij,k}{\mathit{g}}_{ij,k}^{B_{i+1},{\sigma }_y}\left({\tilde{\zeta }}_{ij}\right)\right]-\\ & \sum _{j=1}^{m_{i+1}}\sum _{k=1}^s\left[a_{(i+1)j,k}{\mathit{f}}_{(i+1)j,k}^{B_{i+1,1},{\sigma }_y}\left({\tilde{\mathit{\zeta }}}_{(i+1)j}\right)+b_{(i+1)j,k}{\mathit{g}}_{(i+1)j,k}^{B_{i+1},{\sigma }_y}\left({\tilde{\zeta }}_{(i+1)j}\right)\right]\}=0,\end{array}\\ \begin{array}{rl}\mathrm{I}\mathrm{M}& \{\sum _{j=1}^{m_i}\sum _{k=1}^s\left[a_{ij,k}{\mathit{f}}_{ij,k}^{B_{ij+1},{\tau }_{xy}}\left({\tilde{\zeta }}_{ij}\right)+b_{ij,k}{\mathit{g}}_{ij,k}^{B_{(i+1)1},{\tau }_{ij}}\left({\tilde{\zeta }}_{ij}\right)\right]-\\ & \sum _{j=1}^{m_{i+1}}\sum _{k=1}^s\left[a_{(i+1)j,k}{\mathit{f}}_{(i+1)j,k}^{B_{i+1},{\tau }_{iy}}\left({\tilde{\mathit{\zeta }}}_{(i+1)j}\right)+b_{(i+1)j,k}{\mathit{g}}_{(i+1)j,k}^{B_{i+1},{\tau }_{ij}}\left({\tilde{\mathit{\zeta }}}_{(i+1)j}\right)\right]\}=0,\end{array}\end{array}

when i=t-1\mathrm{o}\mathrm{r}t, multivalued logarithmic terms exist in the analytic functions {\phi }_{ij}({\zeta }_{ij}) and {\psi }_{ij}({\zeta }_{ij}), the term (-1{)}^{t-i}\left[a_{t{m}_t,0}{f}_{t{m}_t,0}^{B_{t1},{\sigma }_y}\left({\tilde{\zeta }}_{t{m}_t}\right)+{\overline{a}}_{t{m}_t,0}{\overline{f}}_{t{m}_t,0}^{B_{t1},{\sigma }_y}\left({\tilde{\zeta }}_{t{m}_t}\right)\right] must be added to the real component of Eq. (31), and the term (-1{)}^{t-i}\left[a_{t{m}_t,0}{f}_{t{m}_t,0}^{B_{t1},{\tau }_{xy}}\left({\tilde{\zeta }}_{t{m}_t}\right)+{\overline{a}}_{t{m}_t,0}{\overline{f}}_{t{m}_t,0}^{B_{t1},{\tau }_{xy}}\left({\tilde{\zeta }}_{t{m}_t}\right)\right] must be added to the imaginary component.

Substituting Eqs. (A.2) and (A.7) into Eq. (31) by considering Eq. (A.3) yields

(32)\left({\mathit{K}}_i^{B_{i+1},S},-{\mathit{K}}_{i+1}^{B_{i+1},S}\right)\left(\begin{array}{l}{\mathit{C}\mathit{o}\mathit{e}\mathit{f}}_i^{\mathrm{T}}\\ {\mathit{C}\mathit{o}\mathit{e}\mathit{f}}_{i+1}^{\mathrm{T}}\end{array}\right)={\mathit{O}}_{2N\times 1}.

For the ground interfaces, by considering the displacement continuity condition in Eq. (23), Eq. (B.1) can be rewritten as (for i\ne t-1\mathrm{o}\mathrm{r}t)

(33)\{\begin{array}{l}\begin{array}{rl}\mathrm{R}\mathrm{e}& \{\sum _{j=1}^{m_i}\sum _{k=1}^s\left[a_{ij,k}{\mathit{f}}_{ij,k}^{B_{(i+1)1},u_x}\left({\tilde{\mathit{\zeta }}}_{ij}\right)+{\overline{a}}_{ij,k}{\overline{\mathit{f}}}_{ij,k}^{B_{(i+1)1},u_x}\left({\tilde{\mathit{\zeta }}}_{ij}\right)+{\overline{b}}_{ij,k}{\overline{\mathit{g}}}_{ij,k}^{B_{(i+1)1},u_x}\left({\tilde{\mathit{\zeta }}}_{ij}\right)\right]+\left(A_i-A_{i+1}\right)\mathit{P}\\ & -\sum _{j=1}^{m_{i+1}}\sum _{k=1}^s\left[a_{(i+1)j,k}{\mathit{f}}_{(i+1)j,k}^{B_{(i+1)1},u_x}\left({\tilde{\mathit{\zeta }}}_{(i+1)j}\right)+{\overline{a}}_{(i+1)j,k}{\overline{\mathit{f}}}_{(i+1)j,k}^{B_{(i+1)1},u_x}\left({\tilde{\mathit{\zeta }}}_{(i+1)j}\right)+{\overline{b}}_{(i+1)j,k}{\overline{\mathit{g}}}_{(i+1)j,k}^{B_{(i+1)1},u_x}\left({\tilde{\mathit{\zeta }}}_{(i+1)j}\right)\right]\}=0,\end{array}\\ \begin{array}{rl}\mathrm{I}\mathrm{m}& \{\sum _{j=1}^{m_i}\sum _{k=1}^s\left[a_{ij,k}{\mathit{f}}_{ij,k}^{B_{(i+1)1},u_y}\left({\tilde{\mathit{\zeta }}}_{ij}\right)+{\overline{a}}_{ij,k}{\overline{\mathit{f}}}_{ij,k}^{B_{(i+1)1},u_y}\left({\tilde{\mathit{\zeta }}}_{ij}\right)+{\overline{b}}_{ij,k}{\overline{\mathit{g}}}_{ij,k}^{B_{(i+1)1},u_y}\left({\tilde{\mathit{\zeta }}}_{ij}\right)\right]+\left(A_i-A_{i+1}\right)\mathit{P}\\ & -\sum _{i+1}\sum _{j=1}^s\left[a_{(i+1)j,k}{\mathit{f}}_{(i+1)j,k}^{B_{(i+1)1},u_y}\left({\tilde{\mathit{\zeta }}}_{(i+1)j}\right)+{\overline{a}}_{(i+1)j,k}{\overline{\mathit{f}}}_{(i+1)j,k}^{B_{(i+1)1},u_y}\left({\tilde{\mathit{\zeta }}}_{(i+1)j}\right)+{\overline{b}}_{(i+1)j,k}{\overline{\mathit{g}}}_{(i+1)j,k}^{B_{(i+1)1},u_y}\left({\tilde{\mathit{\zeta }}}_{(i+1)j}\right)\right]\}=0.\end{array}\end{array}

Note that when (for i\ne t-1\mathrm{o}\mathrm{r}t), another term for (-1{)}^{t-i}\left[a_{t{m}_t,0}{f}_{t{m}_t,0}^{B_{t1},u_x}\left({\tilde{\zeta }}_{t{m}_t}\right)+{\overline{a}}_{t{m}_t,0}{\overline{f}}_{t{m}_t,0}^{B_{t1},u_x}\left({\tilde{\zeta }}_{t{m}_t}\right)\right] must be added to the real component of Eq. (31), and (-1{)}^{t-i}\left[a_{t{m}_t,0}{f}_{t{m}_t,0}^{B_{t1},u_y}\left({\tilde{\zeta }}_{t{m}_t}\right)+{\overline{a}}_{t{m}_t,0}{\overline{f}}_{t{m}_t,0}^{B_{t1},u_y}\left({\tilde{\zeta }}_{t{m}_t}\right)\right] must be added to the imaginary component.

Substituting Eqs. (B.2) and (B.7) into Eq. (31) by considering Eq. (B.3) yields

(34)\left({\mathit{K}}_i^{B_{(i+1)1},D},-{\mathit{K}}_{i+1}^{B_{(i+1)1},D}\right)\left(\begin{array}{c}{\mathit{C}\mathit{o}\mathit{e}\mathit{f}}_i^{\mathrm{T}}\\ {\mathit{C}\mathit{o}\mathit{e}\mathit{f}}_{i+1}^{\mathrm{T}}\end{array}\right)={\mathit{O}}_{2N\times 1}.

Combining Eqs. (32) and (34) yields

(35)\{\begin{array}{c}\left({\mathit{K}}_i^{B_{(i+1)1}},-{\mathit{K}}_{i+1}^{B_{(i+1)1}}\right)\left(\begin{array}{c}{\mathit{C}\mathit{o}\mathit{e}\mathit{f}}_i^{\mathrm{T}}\\ {\mathit{C}\mathit{o}\mathit{e}\mathit{f}}_{i+1}^{\mathrm{T}}\end{array}\right)={\mathit{O}}_{4N\times 1},\\ \begin{array}{cc}{\mathit{K}}_i^{B_{(i+1)1}}=\left(\begin{array}{c}{\mathit{K}}_i^{B_{(i+1)1},S}\\ {\mathit{K}}_i^{B_{(i+1)1},D}\end{array}\right),& {\mathit{K}}_{i+1}^{B_{(i+1)1}}=\left(\begin{array}{c}{\mathit{K}}_{i+1}^{B_{(i+1)1},S}\\ {\mathit{K}}_{i+1}^{B_{(i+1)1},D}\end{array}\right)\end{array}.\end{array}

3.3.4 Linear equations for tunnel boundary

For the tunnel boundary, considering the displacement boundary expressed in Eq. (24), Eq. (B.6) can be rewritten as

(36)\{\begin{array}{rl}& \mathrm{R}\mathrm{e}\left\{\sum _{j=1}^{m_t}\sum _{k=1}^s\left[a_{tj,k}{\mathit{f}}_{tj,k}^{B_{t{m}_t},u_x}\left({\tilde{\mathit{\zeta }}}_{tj}\right)+{\overline{a}}_{tj,k}{\overline{\mathit{f}}}_{tj,k}^{B_{t{m}_t},u_x}\left({\tilde{\mathit{\zeta }}}_{tj}\right)+{\overline{b}}_{tj,k}{\overline{\mathit{g}}}_{tj,k}^{B_{t{m}_t},u_x}\left({\tilde{\mathit{\zeta }}}_{tj}\right)\right]+a_{t{m}_t,0}{\mathit{f}}_{t{m}_t,0}^{B_{t{m}_t},u_x}\left({\tilde{\mathit{\zeta }}}_{t{m}_t}\right)+{\overline{a}}_{t{m}_t,0}{\overline{\mathit{f}}}_{t{m}_t,0}^{B_{t{m}_t},u_x}\left({\tilde{\mathit{\zeta }}}_{t{m}_t}\right)+A_t\mathit{P}\right\}={\tilde{\mathit{u}}}_x\left({\tilde{\mathit{\zeta }}}_{t{m}_t}\right),\\ & \mathrm{I}\mathrm{m}\left\{\sum _{j=1}^{m_t}\sum _{k=1}^s\left[a_{tj,k}{\mathit{f}}_{tj,k}^{B_{t{m}_t},u_y}\left({\tilde{\mathit{\zeta }}}_{tj}\right)+{\overline{a}}_{tj,k}{\overline{\mathit{f}}}_{tj,k}^{B_{t{m}_t},u_y}\left({\tilde{\mathit{\zeta }}}_{tj}\right)+{\overline{b}}_{tj,k}{\overline{\mathit{g}}}_{tj,k}^{B_{t{m}_t},u_y}\left({\tilde{\mathit{\zeta }}}_{tj}\right)\right]+a_{t{m}_t,0}{\mathit{f}}_{t{m}_t,0}^{B_{t{m}_t},u_y}\left({\tilde{\mathit{\zeta }}}_{t{m}_t}\right)+{\overline{a}}_{t{m}_t,0}{\overline{\mathit{f}}}_{t{m}_t,0}^{B_{t{m}_t},u_y}\left({\tilde{\mathit{\zeta }}}_{t{m}_t}\right)+A_t\mathit{P}\right\}={\tilde{\mathit{u}}}_y\left({\tilde{\mathit{\zeta }}}_{t{m}_t}\right).\end{array}

Substituting Eqs. (B.2) and (B.7) into Eq. (36) by considering Eq. (B.3) yields

(37){\mathit{K}}_t^{B_{t{m}_t},D}{\mathit{C}\mathit{o}\mathit{e}\mathit{f}}_t^{\mathrm{T}}={\left(\begin{array}{l}{\tilde{\mathit{u}}}_x\left({\tilde{\mathit{\zeta }}}_{t{m}_t}\right)\\ {\tilde{\mathit{u}}}_y\left({\tilde{\mathit{\zeta }}}_{t{m}_t}\right)\end{array}\right)}_{2N\times 1}.

3.3.5 Establishment of whole linear equations

Two symmetric points on the ground surface with a vertical displacement of 0 should be assumed to obtain the undetermined complex constant A_i for each ground layer. Assuming that the two points located at a distance x₀ from the coordinate origin are additional calculation points, their global coordinate vector is {\tilde{\mathit{z}}}^{\text{add}}={\left(-x_0,x_0\right)}^{\text{T}}. Therefore, the local mapping vector {\tilde{\mathit{\zeta }}}_{1j}^{\text{add}} can be obtained as follows:

(38){\mathit{\omega }}_{1j}\left({\tilde{\mathit{\zeta }}}_{1j}^{\mathrm{a}\mathrm{d}\mathrm{d}}\right)+v_{1j}\mathit{P}={\tilde{\mathit{\xi }}}^{\mathrm{a}\mathrm{d}\mathrm{d}}.

Using the linear equations for the tunnel boundary, the matrix for the additional vertical displacement condition at the ground surface can be expressed as

(39){\mathit{K}}_1^{{\mathrm{B}}_{11},\mathrm{D}}{\mathit{C}\mathit{o}\mathit{e}\mathit{f}}_1^{\mathrm{T}}={\mathbf{O}}_{2\times 1}.

Subsequently, by combining Eqs. (29), (35), (37), and (39), the linear equations can be expressed as

(40){\mathit{K}}_{(4nN+2)\times (8n\cdot s+2n+2)}{{\mathit{C}\mathit{o}\mathit{e}\mathit{f}}^{\mathrm{T}}}_{(8n\cdot s+2n+2)\times 1}=(\mathit{U}{)}_{(4nN+2)\times 1}.

The explicit expression of Eq. (40) is

(41)\left(\begin{array}{ccccccc}{\mathit{K}}_1^{B_{11},S}& \mathbf{O}& \cdots & \mathbf{O}& \cdots & \mathbf{O}& \mathbf{O}\\ {\mathit{K}}_1^{B_{21}}& -{\mathit{K}}_2^{B_{21}}& \cdots & \mathbf{O}& \cdots & \mathbf{O}& \mathbf{O}\\ ⋮& ⋮& & ⋮& & ⋮& ⋮\\ \mathbf{O}& \mathbf{O}& \cdots & \mathbf{O}& \cdots & {\mathit{K}}_{n-1}^{B_{n1}}& -{\mathit{K}}_n^{B_{n1}}\\ \mathbf{O}& \mathbf{O}& \cdots & {\mathit{K}}_t^{B_{t{m}_t},D}& \cdots & \mathbf{O}& \mathbf{O}\\ {\mathit{K}}_1^{B_{11},D}& \mathbf{O}& \cdots & \mathbf{O}& \cdots & \mathbf{O}& \mathbf{O}\end{array}\right)\left(\begin{array}{c}{\mathit{C}\mathit{o}\mathit{e}\mathit{f}}_1^{\mathrm{T}}\\ {\mathit{C}\mathit{o}\mathit{e}\mathit{f}}_2^{\mathrm{T}}\\ ⋮\\ {\mathit{C}\mathit{o}\mathit{e}\mathit{f}}_t^{\mathrm{T}}\\ ⋮\\ {\mathit{C}\mathit{o}\mathit{e}\mathit{f}}_{n-1}^{\mathrm{T}}\\ {\mathit{C}\mathit{o}\mathit{e}\mathit{f}}_n^{\mathrm{T}}\end{array}\right)=\left(\begin{array}{c}{\stackrel{\mathbf{~}}{\mathit{\sigma }}}_y({\stackrel{\mathbf{~}}{\mathit{\zeta }}}_{11})\\ {\stackrel{\mathbf{~}}{\mathit{\tau }}}_{xy}({\stackrel{\mathbf{~}}{\mathit{\zeta }}}_{11})\\ {\mathbf{O}}_{4(n-1)N\times 1}\\ {\stackrel{\mathbf{~}}{\mathit{u}}}_x({\stackrel{\mathbf{~}}{\mathit{\zeta }}}_{t{m}_t})\\ {\stackrel{\mathbf{~}}{\mathit{u}}}_y({\stackrel{\mathbf{~}}{\mathit{\zeta }}}_{t{m}_t})\\ {\mathbf{O}}_{2\times 1}\end{array}\right).

Based on Eq. (41), (8ns + 2n + 2) undetermined real coefficients exist in the objective problem. Therefore, the determination of these coefficients requires at least (8ns + 2n + 2) linear equations. However, the number of linear equations present is (4nN + 2). By performing the abovementioned selection steps for the calculated points, the number of calculated points on each boundary is N ≥ 2s + 1. This indicates that the number of equations is greater than the unknown number, which satisfies the basic requirements for a definite solution of unknown coefficients. Moreover, under certain conditions, increasing the power term s of the local mapping analytic function can improve the accuracy of the analytical solution, and a reasonable value of s can be determined by considering the calculation accuracy and efficiency comprehensively, which will be discussed later. Subsequently, when all the complex coefficient vectors are determined using Eq. (41), the mechanical responses can be obtained using Eqs. (12) and (13). Notably, the convergence pattern of the tunnel is not necessarily limited to that shown in Fig.1; in fact, the uniform radial convergence presented by Verruijt and Booker [28] and the non-uniform convergence presented by Loganathan and Poulos [43] can be used in the proposed solutions.

4 Validation of proposed method

4.1 Boundary and continuity verification

The additional mechanical responses of multi-layered ground due to tunnel excavation are examined since the in situ stress field is relatively simple. However, the boundary and continuity conditions of stress and displacement at the interface should be validated prior to determining the appropriate term. Because of the algebraic complexity involved, a computer program was prepared in MATLAB to perform computations for the analytical solution. A three-layered ground with a tunnel located in the bottom layer was considered, where the ground surface load was set to zero. The tunnel radius, tunnel center depth, and gap parameter were set as R=6\text{m}, h_0=25\text{m}, and Gap=-0.06\text{m}, respectively. The ground parameters are listed in Tab.1.

First, to select a reasonable term number for the Taylor series, the effects of the term number on the accuracy of the analytical solutions and the computational time were investigated, as shown in Fig.5. The maximum error stresses at the ground surface were low when the number of terms was 2^m−1. However, as the number of terms increased, the computational time increased with the emergence of impulse jumps at 2^m−1. When m = 7, the maximum error stresses were {\sigma }_y=1.5\text{Pa} and {\tau }_{xy}= 2.1 Pa, based on a computational time of 6.0 s. Therefore, to balance accuracy and computational time, we set the term number of the Taylor series to 2⁷−1 = 127.

To validate the self-consistency of the analytical solutions, we compared the input boundary conditions and analytical results (see Fig.6 and Fig.7). The analytical results satisfy the boundary conditions, indicating that the proposed non-iterative analytical method and derivation process are reliable.

4.2 Comparison with in situ monitoring results

Based on two cases, i.e., the Green Park tunnel [44,45] and Bangkok Sewer tunnel [46], the proposed method was compared with other methods and the in situ data to verify the applicability of the proposed solutions to actual situations. The ground parameters are listed in Tab.2, and the gap parameters used by Loganathan and Poulos [43] and Park [41] were adopted for the analysis.

The results of ground settlement and horizontal displacement were compared (see Fig.8 and Fig.9). For both cases characterized by soft-over-stiff ground conditions, the maximum surface settlement predicted by the proposed method was slightly higher than that predicted using the Park formula; however, the width of the settlement trough was narrower, and the settlement range was more similar to the observed data. This observation indicates that the proposed method can reflect the effects of the ground parameters on the ground settlement. In fact, the proposed method can obtain reasonable prediction results under significantly different ground stiffness conditions. For subsurface settlement, the results of both the proposed method and Park formula are similar to the observed data. In contrast, the values obtained using the Loganathan–Poulos (L&P) formula [43] are lower than the measured values, particularly for the Bangkok sewer project. Regarding the horizontal displacement, the proposed method can reflect the effect of the ground stiffness, i.e., a step at the ground interface. In the hard lower layer, the horizontal displacement yielded by the proposed method is more similar to the measured value compared with those yielded by the other two methods. Meanwhile, in the soft upper layer, the values yielded by the proposed method were slightly higher than those yielded by the L&P and Park formulae. The horizontal displacement trend predicted using the proposed method was similar to those predicted using the L&P and Park formulae.

In both cases, the L&P formula predicted a narrower settlement trough than the Park formula. The proposed method can be explained by the boundary condition applied by L&P [43], where the surface ground loss at the horizontal distance (h₀ + R) is 25% of the cumulative equivalent ground loss. The results obtained using the Cao and Zeng formulas [27,47] and the proposed method were approximate but similar to the observed data. Generally, the maximum values of the ground surface settlement, subsurface settlement, and horizontal displacement yielded by the proposed method are consistent with field observations and validates the proposed method.

Furthermore, if the mechanical parameters of different layers are identical, the proposed model will degenerate to the classical problem of “deforming a tunnel into an elastic half-plane”. When using uniform radial tunnel convergence under the unified parameters of different layers, the proposed solutions can be transformed into the classical Verruijt solution [29].

4.3 Comparison with numerical simulation

The analytical results were compared with the numerical results to validate the analytical model further. Numerical simulations were conducted using the finite element code ABAQUS under plane-strain conditions and assuming full adherence between layers, which is consistent with the analytical model. Two ground loading patterns were considered for model symmetry, as shown in Fig.10. Models A and B refer to asymmetric and symmetric surface loading about the tunnel axis, respectively. The monitoring path was set as a concentric circle with a radius of 20 m in the tunnel such that the responses of the typical points at layer interfaces and within the ground layer were included in the comparison. The physical and mechanical parameters of the ground are listed in Tab.1, and the other tunnel-related parameters are the same as those presented in Subsection 4.1.

Fig.11 and Fig.12 show a comparison of the ground responses between the analytical solution and numerical results. For Models A and B, the analytical solution agreed well with the numerical results. Both the analytical model and numerical simulations show that the horizontal stress {\sigma }_x at points {\alpha }_1, {\alpha }_2, {\alpha }_3, and {\alpha }_4 indicated the step phenomenon at different degrees, which implies that {\sigma }_x is not continuous at the ground interfaces. In fact, for multi-layered ground, the step phenomenon of {\sigma }_x is caused by the mechanical difference of the ground on both sides of the interface, which is different from the case for a continuous medium. Under symmetrically distributed surface loads, the normal stresses {\sigma }_x and {\sigma }_y are symmetric around the tunnel centerline, whereas the shear stress {\tau }_{xy} is antisymmetric. When the surface loads are asymmetrically distributed, the symmetry properties of all the stresses are broken. Because the magnitude of {\tilde{\sigma }}_y is greater than that of {\tilde{\tau }}_{xy}, the effect of surface loads on the vertical stress {\sigma }_y is more significant than that on {\sigma }_x. Because the surface loads are distributed over the right side of the tunnel, the mechanical responses of the right side are affected more significantly. This indicates that when tunneling beneath a surface building, the ground responses near the building should be considered to evaluate the necessity for ground reinforcement [48].

4.4 Comparison with existing analytical method

The proposed analytical model was compared with the classical solution of Verruijt [29]. A two-layered ground system with a tunnel in the second layer was considered, as shown in Fig.13. Because Verruijt’s solution can only be used for homogeneous ground, the surface settlement troughs yielded by the classical model were based on the parameters of either the first or second layer. A comparison of the results is shown in Fig.14. The settlement troughs yielded by the proposed method were consistent with Verruijt’s solution. In fact, the settlement troughs based on Verruijt’s solution were the lower and upper bounds, and the results yielded by the proposed model differed significantly from those of the classical model. Moreover, the settlement in Case 1 was greater than that in Case 2, which cannot be distinguished via the elastic equivalent method. Therefore, the proposed method outperformed the classical model in solving the aforementioned problem.

To illustrate the differences between the proposed method and the classical model more clearly, a comparison of the subsurface settlement results is presented in Fig.14(b). For the second ground layer (y⩾10\text{m}), the Case 1 results based on the proposed method were similar to those obtained from the classical model using the hard parameters. Based on the specified displacement boundary of tunnel convergence, the lower layer deformed with the constraint effect from the upper ground layer. Because the constraint effect from the soft upper layer in Case 1 was weaker than that from the hard upper layer in the classical model, the settlement of the lower ground layer in the proposed method was greater than that in the classical model. However, the attenuation of the displacement propagation in soft ground was more significant than that in hard ground. Therefore, an inflection point appeared at the ground interface when the proposed model was used. For Case 1, as the depth decreased, the settlement decreased more significantly, and the surface settlement was significantly lower compared with the results yielded by the classical model. Therefore, the comparison between Case 2 and the classical model can be similarly interpreted with reference to Case 1.

5 Conclusions

This paper presents an analytical model for predicting the mechanical responses of multi-layered ground due to shallow tunneling with arbitrary ground surface loading. The proposed method is based on the complex variable method, which can yield accurate solutions for the stress and displacement fields of multi-layered ground. Several examples are presented to demonstrate the efficiency and advantages of the proposed analytical solutions. The evaluation of ground stresses and displacement can facilitate the assessment of tunnel stability and surface building safety.

To solve a multi-boundary problem, a series of analytic functions was introduced to express the mechanical responses contributed by different boundaries. Because the conformal mappings of different boundaries are independent, the multi-boundary problem can be addressed by superpositioning multiple single-boundary problems. Therefore, the mechanical responses of any point in the ground can be obtained directly from the sum of the components contributed by different boundaries. Thus, accurate solutions for the ground responses can be obtained without any iteration.

The results showed that the number of terms in the Taylor series significantly affected the accuracy of the proposed method. The reasonable term number was obtained in the series of 2^m−1 (m = 6,7,…), and as m increased, the boundary conditions were satisfied to a greater extent, although the computational time increased. By balancing accuracy and computational time, the number of terms in this study was set to 127. Subsequently, the results yielded by the proposed method were compared with the numerical and in situ monitoring results for validation. Furthermore, another comparison revealed that the proposed model outperformed the classical model in solving the aforementioned problem.

Because the proposed method can predict ground displacement easily, it can be adopted in the preliminary design stage to estimate the effect of shallow tunneling in multi-layered ground on the surrounding urban infrastructures. The prediction accuracy of the proposed method was validated based on field monitoring results obtained from two actual cases, and the results confirmed the applicability of the proposed method. In fact, the proposed method can be used to predict ground movement when field monitoring data are insufficient, although the different ground qualities must be further discussed.

References

[1] Jenck O, Dias D. 3D-finite difference analysis of the interaction between concrete building and shallow tunnelling. Geotechnique, 2004, 54(8): 519–528

[2] Fang Q, Zhang D, Wong L N Y. Shallow tunnelling method (STM) for subway station construction in soft ground. Tunnelling and Underground Space Technology, 2012, 29: 10–30

[3] Zhang D, Fang Q, Hou Y, Li P, Yuen Wong L N. Protection of buildings against damages as a result of adjacent large-span tunneling in shallowly buried soft ground. Journal of Geotechnical and Geoenvironmental Engineering, 2013, 139(6): 903–913

[4] Fang Q, Liu X, Zeng K, Zhang X, Zhou M, Du J. Centrifuge modelling of tunnelling below existing twin tunnels with different types of support. Underground Space, 2022, 7(6): 1125–1138

[5] Shou K J, Napier J A L. A two-dimensional linear variation displacement discontinuity method for three-layered elastic media. International Journal of Rock Mechanics and Mining Sciences, 1999, 36(6): 719–729

[6] Zhang D, Huang H, Hu Q, Jiang F. Influence of multi-layered soil formation on shield tunnel lining behavior. Tunnelling and Underground Space Technology, 2015, 47: 123–135

[7] Yuan Z, Cao Z, Tang H, Xu Y, Wu T. Analytical layer element with a circular cavity and its application in predicting ground vibrations from surface and underground moving sources. Computers and Geotechnics, 2021, 137: 104262

[8] Zhang J, Gao Y, Liu X, Zhang Z, Yuan Y, Mang H. A shield tunneling method for enlarging the diameter of existing tunnels: Experimental investigations. Tunnelling and Underground Space Technology, 2022, 128: 104605

[9] Do N A, Dias D. Tunnel lining design in multi-layered grounds. Tunnelling and Underground Space Technology, 2018, 81: 103–111

[10] Ieronymaki E, Whittle A J, Einstein H H. Comparative study of the effects of three tunneling methods on ground movements in stiff clay. Tunnelling and Underground Space Technology, 2018, 74: 167–177

[11] Chen Z, He C, Xu G, Ma G, Wu D. A case study on the asymmetric deformation characteristics and mechanical behavior of deep-buried tunnel in phyllite. Rock Mechanics and Rock Engineering, 2019, 52(11): 4527–4545

[12] Zheng H, Li P, Ma G, Zhang Q. Experimental investigation of mechanical characteristics for linings of twins tunnels with asymmetric cross-section. Tunnelling and Underground Space Technology, 2022, 119: 104209

[13] Li P, Wang F, Zhang C, Li Z. Face stability analysis of a shallow tunnel in the saturated and multilayered soils in short-term condition. Computers and Geotechnics, 2019, 107: 25–35

[14] Zhang D, Chen S, Wang R, Zhang D, Li B. Behaviour of a large-diameter shield tunnel through multi-layered strata. Tunnelling and Underground Space Technology, 2021, 116: 104062

[15] PeckR B. Deep excavations and tunneling in soft ground. In: Proceedings of 7th International Conference on Soil Mechanics and Foundation Engineering. Mexico City: Sociedad Mexicana de Mecanica de Suelo, 1969, 225–290

[16] Selby A R. Surface movements caused by tunnelling in two-layer soil. Geological Society Engineering Geology Special Publication, 1988, 5(1): 71–77

[17] Sun Z, Zhang D, Li A, Lu S, Tai Q, Chu Z. Model test and numerical analysis for the face failure mechanism of large cross-section tunnels under different ground conditions. Tunnelling and Underground Space Technology, 2022, 130: 104735

[18] Zhang J, Liu X, Ren T, Shi Y, Yuan Y. Numerical analysis of tunnel segments strengthened by steel-concrete composites. Underground Space, 2022, 7(6): 1115–1124

[19] Di Q, Li P, Zhang M, Wu J. Influence of permeability anisotropy of seepage flow on the tunnel face stability. Underground Space, 2022, 8: 1–14

[20] Di Q, Li P, Zhang M, Cui X. Investigation of progressive settlement of sandy cobble strata for shield tunnels with different burial depths. Engineering Failure Analysis, 2022, 141: 106708

[21] Di Q, Li P, Zhang M, Cui X. Influence of relative density on deformation and failure characteristics induced by tunnel face instability in sandy cobble strata. Engineering Failure Analysis, 2022, 141: 106641

[22] Sun Z, Zhang D, Fang Q, Dui G, Tai Q, Sun F. Analysis of the interaction between tunnel support and surrounding rock considering pre-reinforcement. Tunnelling and Underground Space Technology, 2021, 115: 104074

[23] Sun Z, Zhang D, Fang Q, Liu D, Dui G. Displacement process analysis of deep tunnels with grouted rockbolts considering bolt installation time and bolt length. Computers and Geotechnics, 2021, 140: 104437

[24] Sun Z, Zhang D, Fang Q, Dui G, Chu Z. Analytical solutions for deep tunnels in strain-softening rocks modeled by different elastic strain definitions with the unified strength theory. Science China. Technological Sciences, 2022, 65(10): 2503–2519

[25] Zhang Z, Huang M, Zhang M. Theoretical prediction of ground movements induced by tunnelling in multi-layered soils. Tunnelling and Underground Space Technology, 2011, 26(2): 345–355

[26] Zymnis D M, Chatzigiannelis I, Whittle A J. Effect of anisotropy in ground movements caused by tunnelling. Geotechnique, 2013, 63(13): 1083–1102

[27] Cao L, Zhang D, Fang Q. Semi-analytical prediction for tunnelling-induced ground movements in multi-layered clayey soils. Tunnelling and Underground Space Technology, 2020, 102: 103446

[28] Verruijt A, Booker J R. Surface settlements due to deformation of a tunnel in an elastic half plane. Geotechnique, 1996, 46(4): 753–756

[29] Verruijt A. A complex variable solution for a deforming circular tunnel in an elastic half-plane. International Journal for Numerical and Analytical Methods in Geomechanics, 1997, 21(2): 77–89

[30] Fang Q, Song H, Zhang D. Complex variable analysis for stress distribution of an underwater tunnel in an elastic half plane. International Journal for Numerical and Analytical Methods in Geomechanics, 2015, 39(16): 1821–1835

[31] Zhang D, Xu T, Fang H, Fang Q, Cao L, Wen M. Analytical modeling of complex contact behavior between rock mass and lining structure. Journal of Rock Mechanics and Geotechnical Engineering, 2022, 14(3): 813–824

[32] Fang H, Zhang D, Fang Q. A semi-analytical method for frictional contact analysis between rock mass and concrete linings. Applied Mathematical Modelling, 2022, 105: 17–28

[33] ZhangZHuangMPanYLiZMaSZhangY. Time-dependent analyses for ground movement and stress field induced by tunnelling considering rainfall infiltration mechanics. Tunnelling and Underground Space Technology, 2022, 122, 104378

[34] Fang H, Zhang D, Fang Q, Wen M. A generalized complex variable method for multiple tunnels at great depth considering the interaction between linings and surrounding rock. Computers and Geotechnics, 2021, 129: 103891

[35] Katebi H, Rezaei A H, Hajialilue-Bonab M, Tarifard A. Assessment the influence of ground stratification, tunnel and surface buildings specifications on shield tunnel lining loads (by FEM). Tunnelling and Underground Space Technology, 2015, 49: 67–78

[36] Huang H, Zhang D. Resilience analysis of shield tunnel lining under extreme surcharge: Characterization and field application. Tunnelling and Underground Space Technology, 2016, 51: 301–312

[37] Simon H P, Marte G. Effect of surface loading on the hydro-mechanical response of a tunnel in saturated ground. Underground Space, 2016, 1(1): 1–19

[38] Wang H, Chen X, Jiang M, Song F, Wu L. The analytical predictions on displacement and stress around shallow tunnels subjected to surcharge loadings. Tunnelling and Underground Space Technology, 2018, 71: 403–427

[39] Gao X, Wang H, Jiang M. Analytical solutions for the displacement and stress of lined circular tunnel subjected to surcharge loadings in semi-infinite ground. Applied Mathematical Modelling, 2021, 89: 771–791

[40] Zhang Z, Huang M, Pan Y, Jiang K, Li Z, Ma S, Zhang Y. Analytical prediction of time-dependent behavior for tunneling-induced ground movements and stresses subjected to surcharge loading based on rheological mechanics. Computers and Geotechnics, 2021, 129: 103858

[41] Park K H. Analytical solution for tunnelling-induced ground movement in clays. Tunnelling and Underground Space Technology, 2005, 20(3): 249–261

[42] Lee K M, Rowe R K, Lo K Y. Subsidence owing to tunnelling. I. Estimating the gap parameter. Canadian Geotechnical Journal, 1992, 29(6): 929–940

[43] Loganathan N, Poulos H G. Analytical prediction for tunneling-induced ground movements in clays. Journal of Geotechnical and Geoenvironmental Engineering, 1998, 124(9): 846–856

[44] Attewell P B, Farmer I W. Ground deformations resulting from shield tunnelling in London clay. Canadian Geotechnical Journal, 1974, 11(3): 380–395

[45] Rowe R K, Kack G J. Theoretical examination of the settlements induced by tunneling: Four case histories. Canadian Geotechnical Journal, 1983, 20(2): 299–314

[46] PhienwejN. Ground movements in shield tunnelling in Bangkok soils. In: Proceedings of 14th International Conference on Soil Mechanics and Foundation Engineering. Hamburg: International Society for Soil Mechanics and Foundation Engineering, 1997, 1469–1472

[47] ZengBHuangDPengNChenF. Analogous stochastic medium theory method (ASMTM) for predicting soil displacement induced by general and special-section shield tunnel construction. Chinses Journal of Rock Mechanics and Engineering, 2018, 37: 4356–4366 (in Chinese)

[48] SunZZhangDFangQHuangfuNChuZ. Convergenceconfinement analysis for tunnels with combined bolt–cable system considering the effects of intermediate principal stress. Acta Geotechnica, 2022 (in press)

RIGHTS & PERMISSIONS

Higher Education Press