[2]\fnmFu-quan \surChen

1]\orgdivFujian Provincial Key Laboratory of Advanced Technology and Informatization in Civil Engineering, \orgnameFujian University of Technology, \orgaddress\streetXueyuan Road, \cityFuzhou, \postcode350118, \stateFujian, \countryChina

[2]\orgdivCollege of Civil Engineering, \orgnameFuzhou University, \orgaddress\streetXueyuan Road, \cityFuzhou, \postcode350108, \stateFujian, \countryChina

Complex variable solution on asymmetrical sequential shallow tunnelling in gravitational geomaterial considering static equilibrium

\fnmLuo-bin \surLin luobin_lin@fjut.edu.cn phdchen@fzu.edu.cn \fnmChang-jie \surZheng zcj@fjut.edu.cn \fnmShang-shun \surLin linshangshun@fjut.edu.cn [ *

Abstract

Asymmetrical sequential excavation is common in shallow tunnel engineering, especially for large-span tunnels. Owing to the lack of necessary conformal mappings, existing complex variable solutions on shallow tunnelling are only suitable for symmetrical cavities, and can not deal with asymmetrical sequential tunnelling effectively. This paper proposes a new complex variable solution on asymmetrical sequential shallow tunnelling by incorporating a bidirectional conformal mapping scheme consisting of Charge Simulation Method and Complex Dipole Simulation Method. Moreover, to eliminate the far-field displacement singularity of present complex variable method, a rigid static equilibrium mechanical model is established by fixing the far-field ground surface to equilibriate the nonzero resultant along cavity boundary due to graviational shallow tunnelling. The corresponding mixed boundary conditions along ground surface are transformed into homogenerous Riemann-Hilbert problems with extra constraints of traction along cavity boundaries, which are solved in an iterative manner to obtain reasonable stress and displacement fields of asymmetrical sequential shallow tunnelling. The proposed solution is validated by sufficient comparisons with equivalent finite element solution with good agreements. The comparisons also suggest that the proposed solution should be more accurate than the finite element one. A parametric investigation is finally conducted to illustrate possible practical applications of the proposed solution with several engineering recommendations. Additionally, the theoretical improvements and defects of the proposed solution are discussed for objectivity.

keywords:

Sequential shallow tunnelling, Asymmetrical excavation, Static equilibrium, Gravitational gradient, Bidirectional conformal mapping

1 Introduction

Sequential excavation is a widely used and effective construction method in large-span tunnels, for example, the Niayesh Road Tunnel of Niayesh and Sadr highways in Tehran [1], the typical loess tunnel of the Zhengzhou-Xi’an high-speed railway [2], the Hejie Tunnel of the Guiyang-Guangzhou high-speed railway [3], the Yingpan Road Tunnel in Changsha [4], the Laodongshan Tunnel of the Guangtong-Kunming high-speed railway [5], the Lianchengshan Tunnel of the Yinchuan-Kunming National Expressway [6], the Laohushan Tunnel of the Jinan Ring Expressway [7]. According to excavation sequence, the construction methods can be classified into several catagories, such as side drift method, center diagram method, top heading and benche method, upper half vertical subdivision method, and so on. Fig. 1 shows typical side drift and center diagram schemes, and is cited from Ref [1].

Refer to caption — Figure 1: Typical sequential shallow excavation schemes with liners and temporary supports in Ref [1]

These above sequential excavation methods are generally conducted by asymmetrical excavation transversely, and can be investigated using numerical methods (finite element method and finite difference method for instance). Though the numerical methods would deliver satisfactory results of sequential shallow tunnelling, several potential defects exist: (1) Long runtime of numerical models, especially for a full parametric analysis; (2) Multiple modellings for relavant sequential shallow tunnellings of different but similar geometries; (3) Insufficient insight of the mechanism of shallow tunnellings; and (4) Possible license deficiency of necessary commercial numerical softwares (ABAQUS and FLAC for instance). By contrast, analytical methods, such as complex variable method, can be used to adapt the above-mentioned potential defects of the numerical methods in the preliminary design stage of sequential shallow tunnelling.

The most classical complex variable solution on shallow tunnelling is the Verruijt solution [8, 9], in which the Verruijt conformal mapping is proposed to bidirectionally map a lower half plane containing a circular cavity onto a unit annulus. Based on the Verruijt conformal mapping, several extending solutions are developed to estimate the traction along circular cavity boundary [10, 11, 12, 13], to simulate the displacement along circular cavity boundary [14, 15, 16], and to apply ground surface traction [17, 18, 19, 16]. However, no gravitational gradient of geomaterial is considered in the above extending solutions, which is a significant mechanical feature to distinguish shallow tunnelling from deep tunnelling. Strack [20] proposed the first complex variable solution on shallow tunnelling with consideration of gravitational gradient of geomaterial using Verruijt conformal mapping, and more relavant studies followed [21, 22, 23, 24, 25]. The complex variable solutions mentioned above all focus on circular shallow tunnelling, while noncircular tunnels are also widely used in real-world engineering.

To adapt noncircular shallow tunnelling, Zeng et al. [26] proposed an extension of Verruijt conformal mapping by adding finite items in formation of bilateral Laurent series to backwardly map a unit annulus onto a lower half plane containing a noncircular but symmetrical cavity, and more complex variable solutions on noncircular shallow tunnelling by consideration of gravitational gradient are are proposed [27, 28, 29, 30, 31]. The existing complex variable solutions on noncircular shallow tunnelling only focus on symmetrical and full cavity excavation, while asymmetrical and suquential excavation, which is commonly used in large-span tunnels, is rarely discussed. One possible reason is the lack of suitable conformal mapping. In this paper, we introduce a bidirectional conformal mapping scheme incorporating Charge Simulation Method [37] and Complex Dipole Simulation Method [38] to overcome such a mathematical obstable. Both simulation methods were originally put forward to solve electromagnetic problems, and were subsequently found efficient in constructing multiple types of conformal mappings.

Moreover, the above mentioned complex variable solutions considering gravitational gradient [26, 27, 30, 31] generally ignore the violation of the very fundamental static equilibrium owing to the excavation of gravitational geomaterial. The nonzero resultant along cavity boundaries caused by unloading (buoyancy) of excavated geomaerial can not be equilibriated by any counter-acting force, and a displacement singularity would be consequently raised in far-field geomaterial. To establish a rigid static equilibrium, several recent solutions [32, 33, 34, 35] have been put forward to construct new mechanical models by constraining far-field ground surface to artifically generate necessary counter-acting force to equilibriate the nonzero resultant caused by gravitational excavation. However, these recent solutions can not deal with asymmetrical sequential shallow tunnelling effectively.

In this paper, we seek a new complex variable solution to simultaneously solve both above mentioned problems: (i) Asymmetrical sequential shallow tunnelling with variety of cavity shapes, and (ii) Rigid static equilibrium to cancel far-field displacement singularity. With our new complex variable solution, reasonable stress and displacement fields of sequential shallow tunnelling with large-span cross section can be achieved, and the usage of complex variable method can further extended.

2 Typical sequential shallow tunnelling

2.1 Geometrical variations and notations

We may start from a gravitational geomaterial in Fig. 2a. The ground surface is denoted by ${\bm{C}}_{0}$ , the geomaterial region is denoted by ${\bm{\varOmega}}_{0}$ , and the closure $\overline{\bm{\varOmega}}_{0}={\bm{C}}_{0}\cup{\bm{\varOmega}}$ . The violet dash line denotes the possible final cavity boundary to be excavated. Two possible sequential shallow tunnelling schemes different from Fig. 1 are shown in Figs. 2b and 2c.

Fig. 2b shows one possible excavation scheme using top heading and bench method, in which the cavity is excavated from top to bottom stepwisely, and no temporary supports may be needed. However, the symmetrically noncircular cavity geometry of top-heading-bench method in Fig. 2b can already be mathematically studied by existing conformal mappings and complex variable solutions in above mentioned Refs [26, 27, 28, 29, 30, 31], and no obvious mathematical improvements of conformal mapping and complex variable method are thereby necessary.

By contrast, cavity geometry of the multi-stepwise upper half vertical subdivision method in Fig. 2c would no longer hold symmetry, and the exsiting conformal mappings and complex variable solutions in Refs [26, 27, 28, 29, 30, 31] would fail, and new conformal mapping scheme should be developed to form corresponding complex variable solution. The four-stage excavation scheme in Fig. 2c is only a possible example of sequential shallow tunnelling for visualized illustration. Practical excavation schemes may be more complicated in real-world tunnel engineering, and temporary supports are generally needed. In this paper, we would focus on the mechanical variation of gravitational geomaterial alone caused by multi-stepwise upper half vertical subdivision method, and temporary supports are not considered.

To adapt the real-world complicated excavation schemes, we should use abstract notations to present possible excavation schemes in the following notations. A sequential shallow tunnelling is decomposed into $N$ stages. The excavated region for $j^{\rm{th}}$ stage ( $j=1,2,3,\cdots,N$ ) are denoted by ${\bm{D}}_{j}$ , whose boundaries are denoted by by ${\bm{C}}_{j}$ . The above notations of regions and cavity boundaries are abstracted from Figs. 2b and 2c. After $j^{\rm{th}}$ excavation, the remaining geomaterial region is reduced as a closure $\overline{\bm{\varOmega}}_{j}=\overline{\bm{\varOmega}}_{0}-{\bm{D}}_{j}$ . As shown in Fig. 2c, for each stage, the excavated region after $j^{\rm{th}}$ excavation always remains simply-connected, so that the rest geoamterial region $\overline{\bm{\varOmega}}_{j}$ after $j^{\rm{th}}$ excavation should always remain doubly-connected to hold the consistent connectivity.

2.2 Initial stress field and mechanical properties of geomaterial

As shown in Fig. 2a, the geomaterial occupying the lower half plane $\overline{\bm{\varOmega}}_{0}$ is assumed to be of small deformation with elasticity $E$ , Poisson’s ratio $\nu$ , and shear modulus of $G=\frac{E}{2(1+\nu)}$ . A complex rectangular coordinate system $z=x+{\rm{i}}y$ locating at some point of the ground surface ${\bm{C}}_{0}$ . An initial stress field of gravitational gradient $\gamma$ and lateral coefficient $k_{x}$ is orthotropically subjected through the geomaterial as

\left\{\begin{aligned} &\sigma_{x}^{0}(z)=k_{x}\gamma y\\ &\sigma_{y}^{0}(z)=\gamma y\\ &\tau_{xy}^{0}(z)=0\\ \end{aligned}\right.,\quad z=x+{\rm{i}}y\in\overline{\bm{\varOmega}}_{0}

(2.1)

where $\sigma_{x}^{0}$ , $\sigma_{y}^{0}$ , and $\tau_{xy}^{0}$ denote horizontal, vertical, and shear components of the initial stress field, respectively. Before excavation, the displacement in geomaterial is artifically reset to be zero.

2.3 Sequential excavation and static equilibrium

Before any stage of sequential excavation, the traction along cavity boundary ${\bm{C}}_{j}$ can be expressed as

\left\{\begin{aligned} &X_{b}^{0,j}(S_{j})=\sigma_{x}^{0}(S_{j})\cdot\cos\langle\vec{n}_{j},\vec{x}\rangle+\tau_{xy}^{0}(S_{j})\cdot\cos\langle\vec{n}_{j},\vec{y}\rangle\\ &Y_{b}^{0,j}(S_{j})=\tau_{xy}^{0}(S_{j})\cdot\cos\langle\vec{n}_{j},\vec{x}\rangle+\sigma_{y}^{0}(S_{j})\cdot\cos\langle\vec{n}_{j},\vec{y}\rangle\\ \end{aligned}\right.,\quad S_{j}\in{\bm{C}}_{j}

(2.2)

where $X_{b}^{0,j}$ and $Y_{b}^{0,j}$ denote horizontal and vertical components of surface traction of arbitrary point $S_{j}$ along cavity boundary ${\bm{C}}_{j}$ under the initial stress field from $0^{\rm{th}}$ stage to $j^{\rm{th}}$ stage, respectively; $\langle\vec{n}_{j},\vec{x}\rangle$ denotes the angle between outward normal $\vec{n}_{j}$ and $x$ axis, and $\langle\vec{n}_{j},\vec{y}\rangle$ denotes the angle between outward normal $\vec{n}_{j}$ and $y$ axis.

The excavation after $j^{\rm{th}}$ stage is conducted by mechanically applying opposite surface traction of Eq. (2.2) along cavity boundary ${\bm{C}}_{j}$ in the integral formation as

\int_{\bm{C}_{j}}\left[X_{a}^{0,j}(S_{j})+{\rm{i}}Y_{a}^{0,j}(S_{j})\right]{\rm{d}}S_{j}=-\int_{\bm{C}_{j}}\left[X_{b}^{0,j}(S_{j})+{\rm{i}}Y_{b}^{0,j}(S_{j})\right]{\rm{d}}S_{j},\quad S_{j}\in{\bm{C}}_{j}

(2.3)

where $X_{a}^{0,j}$ and $Y_{a}^{0,j}$ denote horizontal and vertical components of opposite surface traction of arbitrary point $S_{j}$ along cavity boundary ${\bm{C}}_{j}$ under the initial stress field from $0^{\rm{th}}$ stage to $j^{\rm{th}}$ stage, respectively.

With the opposite surface traction in Eq. (2.3), the initial geomaterial $\overline{\bm{\varOmega}}_{0}$ is reduced to $\overline{\bm{\varOmega}}_{j}=\overline{\bm{\varOmega}}_{0}-{\bm{D}}_{j}$ . The nonzero resultant of surface traction in Eq. (2.3) can be given as

{\rm{i}}R_{y}^{0,j}=\varointclockwise_{\bm{C}_{j}}\left[X_{a}^{0,j}(S_{j})+{\rm{i}}Y_{a}^{0,j}(S_{j})\right]{\rm{d}}S_{j}={\rm{i}}\gamma\iint_{{\bm{D}}_{j}}{\rm{d}}x{\rm{d}}y

(2.4)

where the resultant is located as arbitrary point $z^{c}_{j}$ within cavity boundary ${\bm{C}}_{j}$ . It should be emphasized that the resultant ${\rm{i}}R_{y}^{0,j}$ is always upward for $j=1,2,3,\cdots,N$ for whatever cavity shape in the sequential excavation stages (see Fig. 2c-3), since the cavity boundary ${\bm{C}}_{j}$ is closed and surrounds the simply-connected region ${\bm{D}}_{j}$ , which is crossed through by the downward potential gravitational gradient $\gamma$ .

The nonzero resultant along cavity boundary ${\bm{C}}_{j}$ would cause the remaining geomaterial of fully free ground surface to be a non-static equilibrium system, and a displacement singularity at infinity would be raised. Detailed mathematical proof procedure can be found in our previous study [32], and no repitition should be conducted in this paper.

To equilibriate the nonzero upward resultant along cavity boundary ${\bm{C}}_{j}$ , the far-field ground surface ${\bm{C}}_{0c}$ is constrained to produce corresponding counter-acting force, and the remaining ground surface is left free and denoted by ${\bm{C}}_{0f}$ . In other words, ${\bm{C}}_{0}={\bm{C}}_{0c}\cup{\bm{C}}_{0f}$ , as shown in Fig. 3. As long as the width of segment ${\bm{C}}_{0f}$ is large enough, the finite segment ${\bm{C}}_{0f}$ can simulate an infinite ground surface. The following static equilibrium and mixed boundary conditions along ground surface can be established as

\int_{\bm{C}_{0c}}\left[X_{a}^{0,j}(T)+{\rm{i}}Y_{a}^{0,j}(T)\right]{\rm{d}}T=-{\rm{i}}R_{y}^{0,j},\quad T\in{\bm{C}}_{0c}

(2.5)

	$u_{c}^{0,j}(T)+{\rm{i}}v_{c}^{0,j}(T)=0,\quad T\in{\bm{C}}_{0c}$		(2.6a)
	$X_{f}^{0,j}(T)+{\rm{i}}Y_{f}^{0,j}(T)=0,\quad T\in{\bm{C}}_{0f}$		(2.6b)

where $u_{c}^{0,j}$ and $v_{c}^{0,j}$ denote horizontal and vertical components of displacement along ground surface segment ${\bm{C}}_{0c}$ due to the exacavation of geomaterial region ${\bm{D}}_{j}$ from $0^{\rm{th}}$ stage to $j^{\rm{th}}$ stage, respectively; $X_{f}^{0,j}$ and $Y_{f}^{0,j}$ dentote horizontal and vertical components of surface traction along ground surface segment ${\bm{C}}_{0f}$ due to the excavation of geomaterial region ${\bm{D}}_{j}$ from $0^{\rm{th}}$ stage to $j^{\rm{th}}$ stage, respectively. Eqs. (2.3) and (2.6) form the necessary boundary conditions for excavation of ${\bm{D}}_{j}$ geoamterial region from $0^{\rm{th}}$ stage to $j^{\rm{th}}$ stage.

Fig. 3 present a corresponding case for graphic illustration of the abstract boundary conditions of sequential excavation in this section using the example of Stage-3 excavation in Fig. 2c-3. In Fig. 3a, the geomterial region ${\bm{D}}_{3}$ is to be excavated, and the surface tractions along cavity boundary ${\bm{C}}_{3}$ before excavation caused by the initial stress field is denoted by $X_{b}^{0,3}$ and $Y_{b}^{0,3}$ . In Fig. 3b, the geomaterial region ${\bm{D}}_{3}$ is excavated, and the remaining geomaterial geometrically reduces from $\overline{\bm{\varOmega}}_{0}$ to $\overline{\bm{\varOmega}}_{3}=\overline{\bm{\varOmega}}_{0}-{\bm{D}}_{3}$ . Meanwhile, the opposite surface tractions $X_{a}^{0,3}$ and $Y_{a}^{0,3}$ are applied along cavity boundary ${\bm{C}}_{3}$ . The upward resultant $R_{y}^{0,3}$ due to excavation of geomaterial region ${\bm{D}}_{3}$ is located at point $z_{3}^{c}$ within cavity boundary ${\bm{C}}_{3}$ . The surface tractions $X_{b}^{0,3}$ and $Y_{b}^{0,3}$ in Fig. 3a and the opposite surface traction $X_{a}^{0,3}$ and $Y_{a}^{0,3}$ in Fig. 3b would cancel each other to free the boundary ${\bm{C}}_{3}$ .

With the mixed boundary conditions in Eqs. (2.3) and (2.6), the incremental stress and displacement fields in the remaining geomaterial can be expressed using the complex potentials as

	$\left\{\begin{aligned} &\sigma_{y}^{0,j}(z)+\sigma_{x}^{0,j}(z)=2\left[\varphi_{0,j}^{\prime}(z)+\overline{\varphi_{0,j}^{\prime}(z)}\right]\\ &\sigma_{y}^{0,j}(z)-\sigma_{x}^{0,j}(z)+2{\rm{i}}\tau_{xy}^{0,j}(z)=2\left[\overline{z}\varphi_{0,j}^{\prime\prime}(z)+\psi_{0,j}^{\prime}(z)\right]\\ \end{aligned}\right.,\quad z\in\overline{\bm{\varOmega}}_{j}$		(2.7a)
	$2G\left[u_{0,j}(z)+{\rm{i}}v_{0,j}(z)\right]=\kappa\varphi_{0,j}(z)-z\overline{\varphi_{0,j}^{\prime}(z)}-\overline{\psi_{0,j}(z)},\quad z\in\overline{\bm{\varOmega}}_{j}$		(2.7b)

where $\sigma_{x}^{0,j}$ , $\sigma_{y}^{0,j}$ , and $\tau_{xy}^{0,j}$ denote horizontal, vertical, and shear components of incremental stress field due to full excavation of ${\bm{D}}_{j}$ geomaterial region from $0^{\rm{th}}$ stage to $j^{\rm{th}}$ stage, respectively; $u_{0,j}$ and $v_{0,j}$ denote horizontal and vertical components of incremental displacement field due to full excavation of ${\bm{D}}_{j}$ geomaterial region from $0^{\rm{th}}$ stage to $j^{\rm{th}}$ stage, respectively; $\varphi_{0,j}^{\prime}(z)$ and $\psi_{0,j}^{\prime}(z)$ denote the complex potentials to be determined using the mixed boundary conditions in Eqs. (2.3) and (2.6); $\bullet^{\prime}$ and $\bullet^{\prime\prime}$ denote the first- and second-order derivatives of abstract function $\bullet$ . For simplicity, the range of $0\leq j\leq N$ is not repeated in the following deductions.

The total stress in geomaterial of $j^{\rm{th}}$ excavation stage can be obtained by accumulation of the incremental one in Eq. (2.7a) and the initial one in Eq. (2.1) as

\left\{\begin{aligned} &\sigma_{x}^{j}(z)=\sigma_{x}^{0}(z)+\sigma_{x}^{0,j}(z)\\ &\sigma_{y}^{j}(z)=\sigma_{y}^{0}(z)+\sigma_{y}^{0,j}(z)\\ &\tau_{xy}^{j}(z)=\tau_{xy}^{0}(z)+\tau_{xy}^{0,j}(z)\\ \end{aligned}\right.,\quad z\in\overline{\bm{\varOmega}}_{j},0\leq j\leq N

(2.7a’)

3 Riemann-Hilbert problem transformation

In this section, the complex variable method for sequential shallow tunnelling defined in Section 2.3 would be presented. To conduct the complex variable method, the physical geomaterial regions $\overline{\bm{\varOmega}}_{j}$ ( $j=1,2,3,\cdots,N$ ) should be bidirectionally mapped onto corresponding mapping unit annuli $\overline{\bm{\omega}}_{j}$ via suitable conformal mappings, so that the boundary conditions, complex potentials, stres and displacement within the physical geomaterial regions $\overline{\bm{\varOmega}}_{j}$ can be presented using the mathematically feasible Laurent series formations in the mapping unit annuli $\overline{\bm{\omega}}_{j}$ to facilitate further computation.

Our previous study [35] proposes a two-step conformal mapping scheme to bidirectionally map a lower half plane containing an arbitrary cavity onto a unit annulus. Each physical geomaterial region $\overline{\bm{\varOmega}}_{j}$ accords with the geometrical and topological requirements of the the two-step conformal mapping in our previous study [35], and can be thereby bidirectionally mapped onto corresponding mapping unit annuli $\overline{\bm{\omega}}_{j}$ as

	$\zeta_{j}=\zeta_{j}(z)=\zeta_{j}\left[w_{j}\left(z\right)\right],\quad z\in\overline{\bm{\varOmega}}_{j}\\$		(3.1a)
	$z=z_{j}(\zeta_{j})=z_{j}\left[w_{j}\left(\zeta_{j}\right)\right],\quad\zeta_{j}\in\overline{\bm{\omega}}_{j}\\$		(3.1b)

The mapping scheme in Eq. (3.1) is briefly explained in Appendix A. Via the bidirectional conformal mapping, the whole geomaterial in the lower half plane $\overline{\bm{\varOmega}}_{j}$ is bidirectionally mapped onto a bounded unit annulus $\overline{\bm{\omega}}_{j}$ ; the cavity boundary ${\bm{C}}_{j}$ is bidirectionally mapped onto the interior periphery of the unit annulus ${\bm{c}}_{j}$ ; the constrained ground surface ${\bm{C}}_{0c}$ and free one ${\bm{C}}_{0f}$ are bidirectionally mapped onto corresponding segments ${\bm{c}}_{0c,j}$ and ${\bm{c}}_{0f,j}$ of exterior periphery of the unit annulus, respectively; the joint points $T_{1}$ and $T_{2}$ connecting ${\bm{C}}_{0c}$ and ${\bm{C}}_{0f}$ are bidirectionally mapped onto joint points $t_{1,j}$ and $t_{2,j}$ connecting ${\bm{c}}_{0c,j}$ and ${\bm{c}}_{0f,j}$ , respectively. Joint points $t_{1j}$ and $t_{2,j}$ in mapping region $\overline{\bm{\omega}}_{j}$ should remain singularities after conformal mapping. We should emphasize that $\zeta_{j}=\rho_{j}\cdot\sigma_{j}=\rho_{j}\cdot{\rm{e}}^{{\rm{i}}\theta_{j}}$ holds in the following deductions.

With the conformal mapping scheme in Eq. (3.1), the complex potentials in Eq. (2.7) can be transformed as

\left\{\begin{aligned} &\varphi_{0,j}(z)=\varphi_{0,j}(\zeta_{j})\\ &\psi_{0,j}(z)=\psi_{0,j}(\zeta_{j})\\ \end{aligned}\right.,\quad\zeta_{j}\in\overline{\bm{\omega}}_{j}

(3.2)

Correspondingly, the incremental stress and displacement in Eq. (2.7) can be expressed using polar formation in mapping plane $\zeta_{j}=\rho_{j}\cdot{\rm{e}}^{{\rm{i}}\theta_{j}}$ as

	$\sigma_{\theta}^{0,j}(\zeta_{j})+\sigma_{\rho}^{0,j}(\zeta_{j})=2\left[\varPhi_{0,j}(\zeta_{j})+\overline{\varPhi_{0,j}(\zeta_{j})}\right],\quad\zeta_{j}\in\overline{\bm{\omega}}_{j}$	(3.3a)
	$\sigma_{\rho}^{0,j}(\zeta_{j})+{\rm{i}}\tau_{\rho\theta}^{0,j}(\zeta_{j})=\left[\varPhi_{0,j}(\zeta_{j})+\overline{\varPhi_{0,j}(\zeta_{j})}\right]-\frac{\overline{\zeta_{j}}}{\zeta_{j}}\left[\frac{z_{j}(\zeta_{j})}{z_{j}^{\prime}(\zeta_{j})}\overline{\varPhi_{0,j}^{\prime}(\zeta_{j})}+\frac{\overline{z_{j}^{\prime}(\zeta_{j})}}{z_{j}^{\prime}(\zeta_{j})}\overline{\varPsi_{0,j}(\zeta_{j})}\right],\quad\zeta_{j}\in\overline{\bm{\omega}}_{j}$	(3.3b)
	$g_{j}(\zeta_{j})=2G\left[u_{0,j}(\zeta_{j})+{\rm{i}}v_{0,j}(\zeta_{j})\right]=\kappa\varphi_{0,j}(\zeta_{j})-z_{j}(\zeta_{j})\overline{\varPhi_{0,j}(\zeta_{j})}-\overline{\psi_{0,j}(\zeta_{j})},\quad\zeta_{j}\in\overline{\bm{\omega}}_{j}$	(3.3c)
where $\sigma_{\theta}^{0,j}$ , $\sigma_{\rho}^{0,j}$ , and $\tau_{\rho\theta}^{0,j}$ denote hoop, radial, and tangential components of the curvlinear stress mapped onto the mapping unit annuli $\overline{\bm{\omega}}_{j}$ , respectively; and
	$\left\{\begin{aligned} &\varPhi_{0,j}(\zeta_{j})=\frac{\varphi_{0,j}^{\prime}(\zeta_{j})}{z_{j}^{\prime}(\zeta_{j})}\\ &\varPsi_{0,j}(\zeta_{j})=\frac{\psi_{0,j}^{\prime}(\zeta_{j})}{z_{j}^{\prime}(\zeta_{j})}\\ \end{aligned}\right.,\quad\zeta_{j}\in\overline{\bm{\omega}}_{j}$

The rectangular stress and displacement components in Eq. (2.7) can be computed using the curvilinear ones in Eq. (3.3) as

	$\left\{\begin{aligned} &\sigma_{y}^{0,j}(z)+\sigma_{x}^{0,j}(z)=\sigma_{\theta}^{0,j}\left[\zeta_{j}(z)\right]+\sigma_{\rho}^{0,j}\left[\zeta_{j}(z)\right]\\ &\sigma_{y}^{0,j}(z)-\sigma_{x}^{0,j}(z)+2{\rm{i}}\tau_{xy}^{0,j}(z)=\left\{\sigma_{\theta}^{0,j}\left[\zeta_{j}(z)\right]-\sigma_{\rho}^{0,j}\left[\zeta_{j}(z)\right]+2{\rm{i}}\tau_{\rho\theta}^{0,j}\left[\zeta_{j}(z)\right]\right\}\cdot\frac{\overline{\zeta_{j}(z)}}{\zeta_{j}(z)}\cdot\frac{\overline{z_{j}^{\prime}[\zeta_{j}(z)]}}{z_{j}^{\prime}[\zeta_{j}(z)]}\end{aligned}\right.,\quad z\in\overline{\bm{\varOmega}}_{j}$		(3.4a)
	$u_{0,j}(z)+{\rm{i}}v_{0,j}(z)=u_{0,j}\left[\zeta_{j}(z)\right]+{\rm{i}}v_{0,j}\left[\zeta_{j}(z)\right],\quad z\in\overline{\bm{\varOmega}}_{j}$		(3.4b)

With the conformal mapping scheme in Eq.(3.1), the mixed boundary conditions along ground surface in Eq. (2.6) can be mapped onto the mapping unit annuli $\overline{\bm{\omega}}_{j}$ as

	$u_{c}^{0,j}(t_{j})+{\rm{i}}v_{c}^{0,j}(t_{j})=0,\quad t_{j}\in{\bm{c}}_{0c,j}$		(3.5a)
	$\sigma_{\rho}^{0,j}(t_{j})+{\rm{i}}\tau_{\rho\theta}^{0,j}(t_{j})=0,\quad t_{j}\in{\bm{c}}_{0f,j}$		(3.5b)

where $t_{j}=\sigma_{j}={\rm{e}}^{{\rm{i}}\theta_{j}}$ .

The segmental traction free boundary condition in Eq. (3.5b) can be used to establish analytic continuation and mixed boundary condition along ground surface regarding first-order derivative of complex potential $\varPhi_{0,j}(\zeta_{j})$ . Substituting Eq. (3.3b) into Eq. (3.5b) yields

\overline{\varPhi_{0,j}(t_{j})}=-\varPhi_{0,j}(t_{j})+\frac{\overline{t}_{j}}{t_{j}}\frac{z_{j}(t_{j})}{z_{j}^{\prime}(t_{j})}\overline{\varPhi_{0,j}^{\prime}(t_{j})}+\frac{\overline{t}_{j}}{t_{j}}\frac{\overline{z_{j}^{\prime}(t_{j})}}{z_{j}^{\prime}(t_{j})}\overline{\varPsi_{0,j}(t_{j})},\quad t_{j}\in{\bm{c}}_{0f,j}

(3.6)

Partially substituting $t_{j}=\overline{t}_{j}^{-1}$ into the items of the right-hand side of Eq. (3.6) yields

\overline{\varPhi_{0,j}(t_{j})}=-\varPhi_{0,j}(\overline{t}_{j}^{-1})+\frac{\overline{t}_{j}}{\overline{t}_{j}^{-1}}\frac{z_{j}(\overline{t}_{j}^{-1})}{z_{j}^{\prime}(\overline{t}_{j}^{-1})}\overline{\varPhi_{0,j}^{\prime}(t_{j})}+\frac{\overline{t}_{j}}{\overline{t}_{j}^{-1}}\frac{\overline{z_{j}^{\prime}(t_{j})}}{z_{j}^{\prime}(\overline{t}_{j}^{-1})}\overline{\varPsi_{0,j}(t_{j})},\quad t_{j}\in{\bm{c}}_{0f,j}

(3.7)

Replacing $t_{j}=\sigma_{j}$ with $\zeta_{j}=\rho_{j}\cdot\sigma_{j}$ in Eq. (3.7) yields

\overline{\varPhi_{0,j}(\zeta_{j})}=-\varPhi_{0,j}(\overline{\zeta}_{j}^{-1})+\overline{\zeta}_{j}^{2}\frac{z_{j}(\overline{\zeta}_{j}^{-1})}{z_{j}^{\prime}(\overline{\zeta}_{j}^{-1})}\overline{\varPhi_{0,j}^{\prime}(\zeta_{j})}+\overline{\zeta}_{j}^{2}\frac{\overline{z_{j}^{\prime}(\zeta_{j})}}{z_{j}^{\prime}(\overline{\zeta}_{j}^{-1})}\overline{\varPsi_{0,j}(\zeta_{j})},\quad\zeta_{j}\in{\bm{\omega}}_{j}^{+}

(3.8)

We should note that all items on the right-hand side of Eq. (3.8) are analytic in region $\alpha_{j}\leq\rho_{j}<1$ , which is contained within the exterior region $\bm{\omega}_{j}^{-}$ . Thus, Eq. (3.8) shows that $\overline{\varPhi_{0,j}(\zeta_{j})}$ should be analytic in region $\alpha_{j}\leq\rho_{j}<1$ , indicating that analytic continuation has been conducted for $\overline{\varPhi_{0,j}(\zeta_{j})}$ , or simply its conjugate $\varPhi_{0,j}(\zeta_{j})$ .

Substituting Eq. (3.8) into Eq. (3.3b) yields

	$\displaystyle\sigma_{\rho}^{0,j}(\zeta_{j})+{\rm{i}}\tau_{\rho\theta}^{0,j}(\zeta_{j})=$	$\displaystyle\;\overline{\zeta}_{j}^{2}\left[\frac{z_{j}(\overline{\zeta}_{j}^{-1})}{z_{j}^{\prime}(\overline{\zeta}_{j}^{-1})}-\frac{1}{\zeta_{j}\overline{\zeta}_{j}}\frac{z_{j}(\zeta_{j})}{z_{j}^{\prime}(\zeta_{j})}\right]\overline{\varPhi_{0,j}^{\prime}(\zeta_{j})}+\overline{\zeta}_{j}^{2}\left[\frac{\overline{z_{j}^{\prime}(\zeta_{j})}}{z_{j}^{\prime}(\overline{\zeta}_{j}^{-1})}-\frac{1}{\zeta_{j}\overline{\zeta}_{j}}\frac{\overline{z_{j}^{\prime}(\zeta_{j})}}{z_{j}^{\prime}(\zeta_{j})}\right]\overline{\varPsi_{j}(\zeta_{j})}$		(3.3b’)
		$\displaystyle\;+\left[\varPhi_{0,j}(\zeta_{j})-\varPhi_{0,j}(\overline{\zeta}_{j}^{-1})\right],\quad\zeta_{j}\in{\bm{\omega}}_{j}^{+}$		(3.3b’)

The first-order of Eq. (3.3c) about $\zeta_{j}$ can be expressed as

\frac{{\rm{d}}g_{j}(\zeta_{j})}{{\rm{d}}\zeta_{j}}=\kappa z_{j}^{\prime}(\zeta_{j})\varPhi_{0,j}(\zeta_{j})-z_{j}^{\prime}(\zeta_{j})\overline{\varPhi_{0,j}(\zeta_{j})}+z_{j}(\zeta_{j})\overline{\varPhi_{0,j}^{\prime}(\zeta_{j})}+\frac{\overline{\zeta}_{j}}{\zeta_{j}}\overline{z_{j}^{\prime}(\zeta_{j})}\;\overline{\varPsi_{0,j}(\zeta_{j})}

(3.9)

Substituting Eq. (3.8) into Eq. (3.9) yields

	$\displaystyle\frac{{\rm{d}}g_{j}(\zeta_{j})}{{\rm{d}}\zeta_{j}}=$	$\displaystyle\;\frac{\overline{\zeta}_{j}}{\zeta_{j}}\left[z_{j}(\zeta_{j})-\zeta_{j}\overline{\zeta}_{j}\frac{z_{j}^{\prime}(\zeta_{j})}{z_{j}^{\prime}(\overline{\zeta}_{j}^{-1})}z(\overline{\zeta}_{j}^{-1})\right]\overline{\varPhi_{0,j}^{\prime}(\zeta_{j})}+\frac{\overline{\zeta}_{j}}{\zeta_{j}}\left[\overline{z_{j}^{\prime}(\zeta_{j})}-\zeta_{j}\overline{\zeta}_{j}\frac{z_{j}^{\prime}(\zeta_{j})}{z_{j}^{\prime}(\overline{\zeta}_{j}^{-1})}\overline{z_{j}^{\prime}(\zeta_{j})}\right]\overline{\varPsi_{0,j}(\zeta_{j})}$		(3.3c’)
		$\displaystyle\;+\left[\kappa z_{j}^{\prime}(\zeta_{j})\varPhi_{0,j}(\zeta_{j})+z_{j}^{\prime}(\zeta_{j})\varPhi_{0,j}(\overline{\zeta}_{j}^{-1})\right],\quad\zeta_{j}\in{\bm{\omega}}_{j}^{+}$		(3.3c’)

We should notice that when $\zeta_{j}\rightarrow\sigma_{j}$ (or $\rho_{j}\rightarrow 1$ ), the first lines of Eqs. (3.3b’) and (3.3c’) would turn to zero due to $\overline{\sigma}_{j}^{-1}=\sigma_{j}$ , and remaining second lines would form a homogenerous Riemann-Hilbert problem as [32, 34, 35]

	$\left.z_{j}^{\prime}(\zeta_{j})\left[\sigma_{\rho}^{0,j}(\zeta_{j})+{\rm{i}}\tau_{\rho\theta}^{0,j}(\zeta_{j})\right]\right\|_{\rho_{j}\rightarrow 1}=\varphi_{0,j}^{\prime+}(\sigma_{j})-\varphi_{0,j}^{\prime-}(\sigma_{j})=0,\quad\sigma_{j}\in{\bm{c}}_{0f,j}$		(3.10a)
	$\left.g_{j}^{\prime}(\zeta_{j})\right\|_{\rho_{j}\rightarrow 1}=\kappa\varphi_{0,j}^{\prime+}(\sigma_{j})+\varphi_{0,j}^{\prime-}(\sigma_{j})=0,\quad\sigma_{j}\in{\bm{c}}_{0c,j}$		(3.10b)

where $\varphi_{0,j}^{\prime+}(\sigma_{j})$ and $\varphi_{0,j}^{\prime-}(\sigma_{j})$ denote the boundary values of $\varphi_{0,j}^{\prime}(\zeta_{j})=z_{j}^{\prime}(\zeta_{j})\varPhi_{0,j}(\zeta_{j})$ approaching boundary ${\bm{c}}_{0,j}$ from ${\bm{\omega}}^{+}_{j}$ and ${\bm{\omega}}_{j}^{-}$ sides, respectively.

The integrand of the left-hand side of Eq. (2.3) can be mapped from the physical plane $z=x+{\rm{i}}y$ onto corresponding mapping plane $\zeta_{j}=\rho_{j}\cdot{\rm{e}}^{{\rm{i}}\theta_{j}}$ according to the backward conformal mapping in Eq. (3.1b) as

	$X_{a}^{0,j}(S_{j})+{\rm{i}}Y_{a}^{0,j}(S_{j})=\frac{z_{j}^{\prime}(s_{j})}{\|z_{j}^{\prime}(s_{j})\|}\frac{{\rm{d}}s_{j}}{\|{\rm{d}}s_{j}\|}\cdot\left[\sigma_{\rho}^{0,j}(s_{j})+{\rm{i}}\tau_{\rho\theta}^{0,j}(s_{j})\right],\quad s_{j}\in{\bm{c}}_{j}$	(3.11a)
where $s_{j}=\alpha_{j}\cdot\sigma_{j}=\alpha_{j}\cdot{\rm{e}}^{{\rm{i}}\theta_{j}}$ , denoting the point along the mapping contour ${\bm{c}}_{j}$ corresponding to the original point $S_{j}$ along the physical contour ${\bm{C}}_{j}$ . The integrand of the right-hand side of Eq. (2.3) can be expanded according to Eq. (2.2) as
	$X_{b}^{0,j}(S_{j})+{\rm{i}}Y_{b}^{0,j}(S_{j})=-k_{x}\gamma y(S_{j})\frac{{\rm{d}}y(S_{j})}{{\rm{d}}S_{j}}+{\rm{i}}\gamma y(S_{j})\frac{{\rm{d}}y(S_{j})}{{\rm{d}}S_{j}},\quad S_{j}\in{\bm{C}}_{j}$	(3.11b)
where $\cos\langle\vec{n}_{j},\vec{x}\rangle=\frac{{\rm{d}}y}{{\rm{d}}S_{j}}$ and $\cos\langle\vec{n}_{j},\vec{y}\rangle=-\frac{{\rm{d}}x}{{\rm{d}}S_{j}}$ owing to the clockwise integrating direction. The increment ${\rm{d}}S_{j}$ in Eq. (2.3) is clockwise length increment with ${\rm{d}}S_{j}=-\|{\rm{d}}S_{j}\|$ , which can be mapped from the physical plane $z=x+{\rm{i}}y$ onto mapping plane $\zeta_{j}=\rho\cdot{\rm{e}}^{{\rm{i}}\theta_{j}}$ as
	${\rm{d}}S_{j}=-\|{\rm{d}}S_{j}\|=-\|z_{j}^{\prime}(s_{j})\|\cdot\|{\rm{d}}s_{j}\|=-\|z_{j}^{\prime}(s_{j})\|\cdot\alpha_{j}\cdot\|\sigma_{j}\|\cdot\|{\rm{d}}\theta_{j}\|=\|z_{j}^{\prime}(s_{j})\|\cdot\alpha_{j}{\rm{d}}\theta_{j}$	(3.11c)
where $\|{\rm{d}}\theta_{j}\|=-{\rm{d}}\theta_{j}$ , since conformal mapping would not alter the clockwise integration direction.

With Eq. (3.11a) and (3.11c), the left-hand side of Eq. (2.3) can be rearranged as

	${\rm{i}}\int_{\bm{C}_{j}}\left[X_{a}^{0,j}(S_{j})+{\rm{i}}Y_{a}^{0,j}(S_{j})\right]{\rm{d}}S_{j}=\int_{\bm{c}_{j}}z_{j}^{\prime}(s_{j})\left[\sigma_{\rho}^{0,j}(s_{j})+{\rm{i}}\tau_{\rho\theta}^{0,j}(s_{j})\right]\cdot\frac{\|z_{j}^{\prime}(s_{j})\|}{\|z_{j}^{\prime}(s_{j})\|}\cdot\alpha_{j}\cdot{\rm{e}}^{{\rm{i}}\theta_{j}}\cdot{\rm{i}}\cdot{\rm{d}}\theta_{j}$	(3.12a)
where ${\rm{d}}s_{j}=\alpha_{j}\cdot{\rm{e}}^{\rm{i}\theta_{j}}\cdot{\rm{i}}\cdot{\rm{d}}\theta_{j}$ . With Eq. (3.11b), the right-hand side of Eq. (2.3) can be rearranged as
	$-{\rm{i}}\int_{\bm{C}_{j}}\left[X_{b}^{0,j}(S_{j})+{\rm{i}}Y_{b}^{0,j}(S_{j})\right]{\rm{d}}S_{j}=-\gamma\int_{\bm{C}_{j}}y(S_{j})\left[k_{x}{\rm{i}}{\rm{d}}y(S_{j})+{\rm{d}}x(S_{j})\right]$	(3.12b)

Via the backward conformal mapping in Eq. (3.1b), the right-hand side of Eq. (3.12b) can be mapped as

\left\{\begin{aligned} &{\rm{d}}y(S_{j})=\frac{{\rm{d}}y(s_{j})}{{\rm{d}}\sigma_{j}}{\rm{d}}\sigma_{j}\\ &{\rm{d}}x(S_{j})=\frac{{\rm{d}}x(s_{j})}{{\rm{d}}\sigma_{j}}{\rm{d}}\sigma_{j}\\ \end{aligned}\right.

The righ-hand sides of Eqs. (3.12a) and (3.12b) should be equal according to Eq. (2.3), and we should have

\int_{\bm{c}_{j}}z_{j}^{\prime}(s_{j})\left[\sigma_{\rho}^{0,j}(s_{j})+{\rm{i}}\tau_{\rho\theta}^{0,j}(s_{j})\right]{\rm{d}}s_{j}=-\gamma\int_{\bm{c}_{j}}y(s_{j})\left[k_{x}\frac{{\rm{i}}{\rm{d}}y(s_{j})}{{\rm{d}}\sigma_{j}}+\frac{{\rm{d}}x(s_{j})}{{\rm{d}}\sigma_{j}}\right]{\rm{d}}\sigma_{j}

(3.13)

where

\left\{\begin{aligned} &x(s_{j})=\frac{1}{2}\left[\overline{z_{j}(s_{j})}+z_{j}(s_{j})\right]\\ &y(s_{j})=\frac{{\rm{i}}}{2}\left[\overline{z_{j}(s_{j})}-z_{j}(s_{j})\right]\\ \end{aligned}\right.

Eqs. (3.10) and (3.13) are the mixed boundary conditions mapped from Eqs. (2.6) and (2.3) into the mapping plane $\zeta_{j}=\rho_{j}\cdot{\rm{e}}^{{\rm{i}}\theta_{j}}$ , respectively. Eq. (3.10) would form a homogenerous Riemann-Hilbert problem with extra constraint in Eq. (3.13), and the solution procedure would be briefly presented below.

4 Solution procedure of Riemann-Hilbert problem

Before any discussion of the solution, we should emphasize that all the discussion is conducted in the mapping plane $\zeta_{j}=\rho_{j}\cdot{\rm{e}}^{{\rm{i}}\theta_{j}}$ . The general solution of the homogenerous Riemann-Hilbert problem in Eq. (3.10) can be given below. Eq. (3.10a) would be simultaneously satisfied owing to the analytic continuation of Eq. (3.6), and no further discussion should be needed. As for Eq. (3.10b), we construct the following expression to simulate one of its potential solutions as

X_{j}(\zeta_{j})=(\zeta_{j}-t_{1,j})^{-\frac{1}{2}-{\rm{i}}\lambda}(\zeta_{j}-t_{2,j})^{-\frac{1}{2}+{\rm{i}}\lambda},\quad\lambda=\frac{\ln\kappa}{2\pi},\quad\zeta_{j}\in\overline{\bm{\omega}}_{j}\cup{\bm{\omega}}_{j}^{-}

(4.1)

where $t_{1,j}$ and $t_{2,j}$ are the singularities mapped from $T_{1}$ and $T_{2}$ in the lower half geomaterial region $\overline{\bm{\varOmega}}_{j}$ via the bidirectional conformal mapping in Eq. (3.1).

However, our problem consists of not only the mixed boundaries along ground surface in Eq. (3.10), but also the traction boundary condition in Eq. (3.13). The left-hand side of Eq. (3.13) contains both complex potentials $\varphi_{0,j}^{\prime}(\zeta_{j})$ and $\psi_{0,j}^{\prime}(\zeta_{j})$ to be determined, which should be analytic within the mapping region $\overline{\bm{\omega}}_{j}$ . To simulate the possible solution, we may assume the general solution of $\varphi_{0,j}^{\prime}(\zeta_{j})$ as a combination of $X_{j}(\zeta_{j})$ and an arbitrary analytic function, which should be defined within the mapping region $\overline{\bm{\omega}}_{j}$ with two poles of origin and complex infinity. Therefore, the general solution of $\varphi_{0,j}^{\prime}(\zeta_{j})$ can be expressed as

\varphi_{0,j}^{\prime}(\zeta_{j})=X_{j}(\zeta_{j})\sum\limits_{n=-\infty}^{\infty}f_{j,n}\zeta_{j}^{n},\quad\zeta_{j}\in\overline{\bm{\omega}}_{j}

(4.2)

where $f_{j,n}$ denote complex coefficients for $j^{\rm{th}}$ cavity to be determined. Comparing to Eq. (4.1), the definition domain of item $X_{j}(\zeta)$ in Eq. (4.2) is reduced from $\overline{\bm{\omega}}_{j}\cup{\bm{\omega}}_{j}^{-}$ to $\overline{\bm{\omega}}_{j}$ .

The complex items with complex fractional power $X_{j}(\zeta_{j})$ in Eq. (4.1) is difficult to reach computational results, and can be subsequently expanded into Taylor series near the poles origin and complex infinity, respectively, as

X_{j}(\zeta_{j})=\left\{\begin{aligned} &\sum\limits_{k=0}^{\infty}\alpha_{j,k}\zeta_{j}^{k},\quad\zeta_{j}\in{\bm{\omega}}_{j}^{+}\\ &\sum\limits_{k=0}^{\infty}\beta_{j,k}\zeta_{j}^{-k},\quad\zeta_{j}\in{\bm{\omega}}_{j}^{-}\\ \end{aligned}\right.

(4.3)

where

	$\left\{\begin{aligned} &\alpha_{j,0}=-t_{1,j}^{-\frac{1}{2}-{\rm{i}}\lambda}t_{2,j}^{-\frac{1}{2}+{\rm{i}}\lambda}\\ &\alpha_{j,k}=-t_{1,j}^{-\frac{1}{2}-{\rm{i}}\lambda}t_{2,j}^{-\frac{1}{2}+{\rm{i}}\lambda}\cdot(-1)^{k}\left(c_{k}t_{1,j}^{-k}+\overline{c}_{k}t_{2,j}^{-k}+\sum\limits_{l=1}^{k-1}c_{l}\overline{c}_{k-l}\cdot t_{1,j}^{-l}t_{2,j}^{-k+l}\right),k\geq 1\end{aligned}\right.$	(4.4a)
	$\left\{\begin{aligned} &\beta_{j,0}=0\\ &\beta_{j,1}=1\\ &\beta_{j,k}=(-1)^{k-1}\left(c_{k-1}t_{1,j}^{k-1}+\overline{c}_{k-1}t_{2,j}^{k-1}+\sum\limits_{l=1}^{k-2}c_{l}\overline{c}_{k-1-l}\cdot t_{1,j}^{l}t_{2,j}^{k-1-l}\right),\quad k\geq 2\end{aligned}\right.$	(4.4b)
with
	$c_{k}=\prod\limits_{l=1}^{k}\left(\frac{1}{2}-{\rm{i}}\lambda-l\right)/k!$	(4.4c)

Substituting Eq. (4.3) into Eq. (4.2) yields

	$\varphi_{0,j}^{\prime}(\zeta_{j})=\sum\limits_{k=-\infty}^{\infty}A_{j,k}\zeta_{j}^{k},\quad\zeta_{j}\in{\bm{\omega}}_{j}^{+}$		(4.5a)
	$\varphi_{0,j}^{\prime}(\zeta_{j})=\sum\limits_{k=-\infty}^{\infty}B_{j,k}\zeta_{j}^{k},\quad\zeta_{j}\in{\bm{\omega}}_{j}^{-}$		(4.5b)

where

	$A_{j,k}=\sum\limits_{n=-\infty}^{k}\alpha_{j,k-n}f_{j,n}$		(4.6a)
	$B_{j,k}=\sum\limits_{n=k}^{\infty}\beta_{j,-k+n}f_{j,n}$		(4.6b)

Eq. (4.5) provides expanding expressions of complex potential $\varphi_{0,j}^{\prime}(\zeta_{j})$ with undetermined coefficients $f_{j,n}$ in Eq. (4.2).

Now we may provide formation of the other complex potential $\psi_{0,j}^{\prime}(\zeta_{j})$ using Eqs. (4.5a) and (4.5b) below. Eq. (3.6) can be formulated as

\psi_{0,j}^{\prime}(t_{j})=\frac{\overline{t}_{j}}{t_{j}}\overline{\varphi_{0,j}^{\prime}(t_{j})}+\frac{\overline{t}_{j}}{t_{j}}\frac{\overline{z_{j}^{\prime}(t_{j})}}{z_{j}^{\prime}(t_{j})}\varphi_{0,j}^{\prime}(t_{j})+\frac{\overline{z_{j}(t_{j})}z_{j}^{\prime\prime}(t_{j})}{\left[z_{j}^{\prime}(t_{j})\right]^{2}}\varphi_{0,j}^{\prime}(t_{j})-\frac{\overline{z_{j}(t_{j})}}{z_{j}^{\prime}(t_{j})}\varphi_{0,j}^{\prime\prime}(t_{j}),\quad t_{j}\in{\bm{c}}_{0f,j}

(4.7)

Eq. (4.7) can be rewritten into a more compact manner as

\psi_{0,j}^{\prime}(t_{j})=\frac{\overline{t_{j}}}{t_{j}}\overline{\varphi_{0,j}^{\prime}(t_{j})}-{\rm{d}}\left[\frac{\overline{z_{j}(t_{j})}}{z_{j}^{\prime}(t_{j})}\varphi_{0,j}^{\prime}(t_{j})\right]/{\rm{d}}t_{j},\quad t_{j}\in{\bm{c}}_{0f,j}

(4.7’)

If we note that $\overline{z_{j}(t_{j})}=z_{j}(t_{j})$ denotes the $x$ axis in the physical plane $z=x+{\rm{i}}y$ and $\overline{t}_{j}=t_{j}^{-1}$ on the segment ${\bm{c}}_{0f,j}$ of the unit annuli $\overline{\bm{\omega}}_{j}$ , Eq. (4.7’) can be transformed as

\psi_{0,j}^{\prime}(t_{j})=\frac{1}{t_{j}^{2}}\overline{\varphi}_{0,j}^{\prime}(t_{j}^{-1})-{\rm{d}}\left[\frac{z_{j}(t_{j})}{z_{j}^{\prime}(t_{j})}\varphi_{0,j}^{\prime}(t_{j})\right]/{\rm{d}}t_{j},\quad t_{j}\in{\bm{c}}_{0f,j}

(4.8)

In above deductions, the closed-form expression of the other complex potential $\psi_{0,j}^{\prime}(\zeta_{j})$ is still not given, we may assume the following formation as

\psi_{0,j}^{\prime}(\zeta_{j})=\frac{1}{\zeta_{j}^{2}}\overline{\varphi}_{0,j}^{\prime}(\zeta_{j}^{-1})-{\rm{d}}\left[\frac{z_{j}(\zeta_{j})}{z_{j}^{\prime}(\zeta_{j})}\varphi_{0,j}^{\prime}(\zeta_{j})\right]/{\rm{d}}\zeta_{j},\quad\zeta_{j}\in\overline{\bm{\omega}}_{j}

(4.9)

Apparently, when $\zeta_{j}\rightarrow t_{j}\in{\bm{c}}_{0f,j}$ , Eq. (4.9) would turn into its boundary value in Eq. (4.8) inversely using $t_{j}^{-1}=\overline{t}_{j}$ in the first item and $z_{j}(t_{j})=\overline{z_{j}(t_{j})}$ in the second item, respectively. When $\zeta_{j}\in\overline{\bm{\omega}}_{j}\setminus\left\{t_{1,j},t_{2,j}\right\}$ , $\zeta_{j}^{-1}\in{\bm{\omega}}_{j}^{-}\cup{\bm{c}}_{0c,j}\cup{\bm{c}}_{0f,j}$ , and the first item of Eq. (4.9) can be expanded using Eq. (4.5b) as

\psi_{0,j}^{\prime}(\zeta_{j})=\sum\limits_{k=-\infty}^{\infty}\overline{B}_{j,-k-2}\zeta_{j}^{k}-{\rm{d}}\left[\frac{z_{j}(\zeta_{j})}{z_{j}^{\prime}(\zeta_{j})}\varphi_{0,j}^{\prime}(\zeta_{j})\right]/{\rm{d}}\zeta_{j},\quad\zeta_{j}\in\overline{\bm{\omega}}_{j}

(4.9’)

Eq. (4.9’) provides expressions of complex potential $\psi_{0,j}^{\prime}(\zeta_{j})$ with undetermined coefficients $f_{j,n}$ in Eq. (4.2) in the same definition domain. In the following deductions, we should determine coefficients $f_{j,n}$ in Eq. (4.2) using the boundary condition along cavity boundary in Eq. (3.13).

Substituting Eqs. (3.3b), (4.5a), and (4.9’) into the left-hand side of Eq. (3.13) yields

		$\displaystyle\sum\limits_{\begin{subarray}{c}k=-\infty\\ k\neq 0\end{subarray}}^{\infty}\left(A_{j,k-1}\frac{\alpha_{j}^{k}}{k}\sigma_{j}^{k}+B_{j,-k-1}\frac{\alpha_{j}^{k}}{k}\sigma_{j}^{-k}\right)+\frac{z_{j}(\alpha_{j}\sigma_{j})-\overline{z_{j}(\alpha_{j}\sigma_{j})}}{\overline{z_{j}^{\prime}(\alpha_{j}\sigma_{j})}}\sum\limits_{k=-\infty}^{\infty}\overline{A}_{j,k}\alpha_{j}^{k}\sigma_{j}^{-k}$		(4.10)
		$\displaystyle+(A_{j,-1}+B_{j,-1})\ln\alpha_{j}+C_{a,j}+(A_{j,-1}-B_{j,-1}){\rm{Ln}}\sigma_{j}=-\gamma\int_{{\bm{c}}_{j}}y(s_{j})\left[k_{x}\frac{{\rm{i}}{\rm{d}}y(s_{j})}{{\rm{d}}\sigma_{j}}+\frac{{\rm{d}}x(s_{j})}{{\rm{d}}\sigma_{j}}\right]{\rm{d}}\sigma_{j}$		(4.10)

where $C_{a,j}$ denotes an arbitrary complex integral constant, ${\rm{Ln}}$ denotes the possibly multi-valued natural logarithm. The mapping item and the integrand of the right-hand side of Eq. (3.13) can be expanded using the sample point method [34, 35] as

	$\frac{z_{j}(\alpha_{j}\sigma_{j})-\overline{z_{j}(\alpha_{j}\sigma_{j})}}{\overline{z_{j}^{\prime}(\alpha_{j}\sigma_{j})}}=\sum\limits_{k=-\infty}^{\infty}d_{j,k}\sigma_{j}^{k}$		(4.11a)
	$y(s_{j})\left[k_{x}\frac{{\rm{i}}{\rm{d}}y(s_{j})}{{\rm{d}}\sigma_{j}}+\frac{{\rm{d}}x(s_{j})}{{\rm{d}}\sigma_{j}}\right]=\sum\limits_{k=-\infty}^{\infty}h_{j,k}\sigma_{j}^{k}$		(4.11b)

Then the right-hand side of Eq. (3.13) can be integrated using Eq. (4.11b) as

-\gamma\int_{{\bm{c}}_{j}}\left[k_{x}\frac{{\rm{i}}{\rm{d}}y(s_{j})}{{\rm{d}}\sigma_{j}}+\frac{{\rm{d}}x(s_{j})}{{\rm{d}}\sigma_{j}}\right]{\rm{d}}\sigma_{j}=\sum\limits_{\begin{subarray}{c}k=-\infty\\ k\neq 0\end{subarray}}^{\infty}I_{j,k}\sigma_{j}^{k}+I_{j,0}{\rm{Ln}}\sigma_{j}

(4.12)

where

\left\{\begin{aligned} &I_{j,k}=-\gamma\frac{h_{j,k-1}}{k},\quad k\neq 0\\ &I_{j,0}=-\gamma h_{j,-1}\end{aligned}\right.

(4.13)

Substituting Eqs. (4.10), (4.11a), and (4.12) into Eq. (3.13) and comparing the coefficients of same power of $\sigma_{j}$ yields

	$A_{j,-k-1}=-k\alpha_{j}^{k}I_{j,-k}+\alpha_{j}^{2k}B_{j,-k-1}+k\alpha_{j}^{2k}\sum\limits_{l=-\infty}^{\infty}\alpha_{j}^{l}d_{j,k}\overline{A}_{j,l+k},\quad k\geq 1$		(4.14a)
	$B_{j,k-1}=-k\alpha_{j}^{k}I_{j,k}+\alpha_{j}^{2k}A_{j,k-1}+k\sum\limits_{l=-\infty}^{\infty}\alpha_{j}^{l}d_{j,k}\overline{A}_{j,l-k},\quad k\geq 1$		(4.14b)

A_{j,-1}-B_{j,-1}=I_{j,0}

(4.15)

The resultant equilibrium in Eq. (4.15) can be transformed using residual theorem as [32, 33, 34, 35]

A_{j,-1}-B_{j,-1}=-{\rm{i}}\frac{R_{y}^{0,j}}{2\pi}

(4.15’)

Eqs. (4.15) and (4.15’) reveal an implicit equilibrium as

I_{j,0}=-\frac{\rm{i}\gamma}{2\pi}\iint_{\bm{D}_{j}}{\rm{d}}x{\rm{d}}y

(4.16)

Eq. (4.16) should always be numerically examined. Furthermore, Eqs. (4.14) and (4.15) derived from Eqs. (3.10) and (3.13) only determines the first derivatives $\varphi_{0,j}^{\prime}(\zeta_{j})$ and $\psi_{0,j}^{\prime}(\zeta_{j})$ to be analytic and single-valued, however, the displacement in Eq. (3.3c), which contains the original functions $\varphi_{0,j}(\zeta_{j})$ and $\psi_{0,j}(\zeta_{j})$ , should be also analytic and single-valued. To guarantee the single-valuedness of displacement in Eq. (3.3c), the follow equilibrium should be established as [32, 33, 34, 35]

\kappa A_{j,-1}+B_{j,-1}=0

(4.17)

Eqs. (4.15) and (4.17) give the following invariants unaffected by conformal mappings as

	$A_{j,-1}=\frac{I_{j,0}}{1+\kappa}$		(4.18a)
	$B_{j,-1}=\frac{-\kappa I_{j,0}}{1+\kappa}$		(4.18b)

Substituting Eq. (4.6) into the left-hand sides of Eqs. (4.14) and (4.18) gives well-defined simultaneous complex linear system on $f_{n}$ ( $-\infty<n<\infty$ ) as

	$\left\{\begin{aligned} &\sum\limits_{n=1}^{\infty}\alpha_{j,-1+n}f_{j,n}=\frac{I_{j,0}}{1+\kappa}\\ &\sum\limits_{n=k+1}^{\infty}\alpha_{j,-k-1+n}f_{j,-n}=-k\alpha_{j}^{k}I_{j,-k}+\alpha_{j}^{2k}B_{j,-k-1}+k\alpha_{j}^{k}\sum\limits_{l=-\infty}^{\infty}\alpha_{j}^{l}d_{j,k}\overline{A}_{j,l+k},\quad k\geq 1\\ \end{aligned}\right.$		(4.19a)
	$\left\{\begin{aligned} &\sum\limits_{n=0}^{\infty}\beta_{j,1+n}f_{j,n}=\frac{-\kappa I_{j,0}}{1+\kappa}\\ &\sum\limits_{n=k}^{\infty}\beta_{j,-k+1+n}f_{j,n}=-k\alpha_{j}^{k}I_{j,k}+\alpha_{j}^{2k}A_{j,k-1}+k\sum\limits_{l=-\infty}^{\infty}\alpha_{j}^{l}d_{j,k}\overline{A}_{j,l-k},\quad k\geq 1\\ \end{aligned}\right.$		(4.19b)

Eq. (4.19) can be solved in the iterative method in Refs [32, 33, 34, 35] to uniquely determine $f_{n}$ in Eq. (4.2), so that the complex potentials in Eqs. (4.5a) and (4.9’) can be subsequently determined for further solution of stress and displacement of sequential shallow tunnelling.

5 Stress and displacement fields of sequential shallow tunnelling

Replacing $s_{j}=\alpha_{j}\sigma_{j}$ with $\zeta_{j}=\rho_{j}\sigma_{j}$ in Eqs. (3.13) and (4.10) yields the integration of the curvilinear traction components mapped onto the mapping annuli $\overline{\bm{\omega}}_{j}\setminus\left\{t_{1,j},t_{2,j}\right\}$ for a given polar radius $\rho_{j}$ as

	$\displaystyle\int z_{j}^{\prime}(\rho_{j}\sigma_{j})\left[\sigma_{\rho}^{0,j}(\rho_{j}\sigma_{j})+{\rm{i}}\tau_{\rho\theta}^{0,j}(\rho_{j}\sigma_{j})\right]{\rm{d}}(\rho_{j}\sigma_{j})$	(5.1)
$\displaystyle=$	$\displaystyle\sum\limits_{\begin{subarray}{c}k=-\infty\\ k\neq 0\end{subarray}}^{\infty}\left(A_{j,k-1}\frac{\rho_{j}^{k}}{k}-B_{j,k-1}\frac{\rho_{j}^{-k}}{k}\right)\sigma_{j}^{k}+\sum\limits_{k=-\infty}^{\infty}\sum\limits_{l=-\infty}^{\infty}g_{j,l}(\rho_{j})\overline{A}_{j,l-k}\rho_{j}^{l-k}\sigma_{j}^{k}$
	$\displaystyle+\left(A_{j,-1}+B_{j,-1}\right)\ln\rho_{j}+\left(A_{j,-1}-B_{j,-1}\right){\rm{Ln}}\sigma_{j}+C_{a,j},\quad\alpha_{j}\leq\rho_{j}\leq 1$

where

\sum\limits_{k=-\infty}^{\infty}g_{j,k}(\rho_{j})\sigma_{j}^{k}=\frac{z_{j}(\rho_{j}\sigma_{j})-\overline{z_{j}(\rho_{j}\sigma_{j})}}{\overline{z_{j}^{\prime}(\rho_{j}\sigma_{j})}}

(5.2)

When $\rho_{j}\rightarrow\alpha_{j}$ , Eq. (5.1) would be degenerated into Eq. (4.10). When $\rho_{j}\rightarrow 1$ , the right-hand side of Eq. (5.2) would be zero, since $\overline{z_{j}(\sigma_{j})}=z_{j}(\sigma_{j})$ denotes the $x$ axis in the physical plane $z=x+{\rm{i}}y$ .

The derivation of Eq. (5.1) gives the expansion of Eq. (3.3b) as

	$\sigma_{\rho}^{0,j}(\rho_{j}\sigma_{j})+{\rm{i}}\tau_{\rho\theta}^{0,j}(\rho_{j}\sigma_{j})=\frac{1}{z_{j}^{\prime}(\rho_{j}\sigma_{j})}\sum\limits_{k=-\infty}^{\infty}\left(\begin{aligned} &A_{j,k}\rho_{j}^{k}-B_{j,k}\rho_{j}^{-k-2}\\ +&(k+1)\sum\limits_{l=-\infty}^{\infty}g_{j,l}(\rho_{j})\overline{A}_{j,l-k-1}\rho_{j}^{l-k-2}\\ \end{aligned}\right)\sigma_{j}^{k},\quad\alpha_{j}\leq\rho_{j}\leq 1$		(5.3a)
Eqs. (3.3a) and (3.3c) can be respectively expanded as
	$\sigma_{\theta}^{0,j}(\rho_{j}\sigma_{j})+\sigma_{\rho}^{0,j}(\rho_{j}\sigma_{j})=4\Re\left[\frac{1}{z_{j}^{\prime}(\rho_{j}\sigma_{j})}\sum\limits_{k=-\infty}^{\infty}A_{j,k}\rho_{j}^{k}\sigma_{j}^{k}\right],\quad\alpha_{j}\leq\rho_{j}\leq 1$		(5.3b)

	$\displaystyle 2G[u_{0,j}(\rho_{j}\sigma_{j})+{\rm{i}}v_{0,j}(\rho_{j}\sigma_{j})]=$	$\displaystyle\sum\limits_{\begin{subarray}{c}k=-\infty\\ k\neq 0\end{subarray}}^{\infty}\frac{1}{k}\left(\kappa A_{j,k-1}\rho_{j}^{k}\sigma_{j}^{k}+B_{j,k-1}\rho_{j}^{-k}\sigma_{j}^{k}\right)+C_{0}+(\kappa A_{j,-1}-B_{j,-1})\ln\rho_{j}$	(5.3c)
		$\displaystyle-\sum\limits_{k=-\infty}^{\infty}\sum\limits_{l=-\infty}^{\infty}g_{j,l}(\rho_{j})\overline{A}_{j,l-k}\rho_{j}^{l-k}\sigma_{j}^{k},\quad\alpha_{j}\leq\rho_{j}\leq 1$
with
	$C_{0}=-\sum\limits_{\begin{subarray}{c}k=-\infty\\ k\neq 0\end{subarray}}^{k=\infty}\frac{1}{k}\left(\kappa A_{j,k-1}+B_{j,k-1}\right)$

The initial stress field in Eq. (2.1) can be mapped onto the mapping annuli $\overline{\bm{\omega}}_{j}$ as

\left\{\begin{aligned} &\sigma_{\theta}^{0}(\rho_{j}\sigma_{j})+\sigma_{\rho}^{0}(\rho_{j}\sigma_{j})=\sigma_{y}^{0}[z_{j}(\rho_{j}\sigma_{j})]+\sigma_{x}^{0}[z_{j}(\rho_{j}\sigma_{j})]\\ &\sigma_{\theta}^{0}(\rho_{j}\sigma_{j})-\sigma_{\rho}^{0}(\rho_{j}\sigma_{j})+2{\rm{i}}\tau_{\rho\theta}^{0}(\rho_{j}\sigma_{j})=\left\{\sigma_{y}^{0}[z_{j}(\rho_{j}\sigma_{j})]-\sigma_{x}^{0}[z_{j}(\rho_{j}\sigma_{j})]+2{\rm{i}}\tau_{xy}[z(\rho_{j}\sigma_{j})]\right\}\cdot\frac{z_{j}^{\prime}(\rho_{j}\sigma_{j})}{\overline{z_{j}^{\prime}(\rho_{j}\sigma_{j})}}\cdot\sigma_{j}^{2}\end{aligned}\right.,\quad\alpha_{j}\leq\rho_{j}\leq 1

(2.1’)

where $\sigma_{\theta}^{0}$ , $\sigma_{\rho}^{0}$ , and $\tau_{\rho\theta}^{0}$ denote hoop, radial, and shear components of the initial stress field mapped onto the mapping unit annuli $\overline{\bm{\omega}}_{j}$ , respectively.

According to Eqs. (5.3) and (2.1’), the total curvilinear stress field mapped onto mapping unit annuli $\overline{\bm{\omega}}_{j}$ can be obtained as

	$\left\{\begin{aligned} &\sigma_{\theta}^{j}(\rho_{j}\sigma_{j})=\sigma_{\theta}^{0}(\rho_{j}\sigma_{j})+\sigma_{\theta}^{0,j}(\rho_{j}\sigma_{j})\\ &\sigma_{\rho}^{j}(\rho_{j}\sigma_{j})=\sigma_{\rho}^{0}(\rho_{j}\sigma_{j})+\sigma_{\rho}^{0,j}(\rho_{j}\sigma_{j})\\ &\tau_{\rho\theta}^{j}(\rho_{j}\sigma_{j})=\tau_{\rho\theta}^{0}(\rho_{j}\sigma_{j})+\tau_{\rho\theta}^{0,j}(\rho_{j}\sigma_{j})\\ \end{aligned}\right.,\quad\alpha_{j}\leq\rho_{j}\leq 1$	(5.4a)
where $\sigma_{\theta}^{j}$ , $\sigma_{\rho}^{j}$ , and $\tau_{\rho\theta}^{j}$ denote hoop, radial, and tangential components of total curvilinear stress mapped onto the mapping unit annuli $\overline{\bm{\omega}}_{j}$ along the circles of radii $\|\zeta_{j}\|=\rho_{j}$ , respectively. The displacement field mapped onto mapping unit annuli $\overline{\bm{\omega}}_{j}$ can be obtained as
	$u_{j}(\rho_{j}\sigma_{j})+{\rm{i}}v_{j}(\rho_{j}\sigma_{j})=u_{j}[z_{j}(\rho_{j}\sigma_{j})]+{\rm{i}}v_{j}[z_{j}(\rho_{j}\sigma_{j})],\quad\alpha_{j}\leq\rho_{j}\leq 1,\sigma_{j}\neq t_{1,j},t_{2,j}$	(5.4b)

Significantly, when $\rho_{j}\rightarrow\alpha_{j}$ , Eq. (5.4a) gives the Mises stress (numerically equal to the absolute value of $\sigma_{\theta}^{j}(\alpha_{j}\sigma_{j})$ ) and the residual stresses (radial component $\sigma_{\rho}^{j}(\alpha_{j}\sigma_{j})$ and tangential component $\tau_{\rho\theta}^{j}(\alpha_{j}\sigma_{j})$ ) along cavity boundary ${\bm{C}}_{j}$ caused by $j^{\rm{th}}$ stage excavation; Eq. (5.4b) gives the horizontal and vertical deformation along cavity boundary ${\bm{C}}_{j}$ .

The curvilinear stress and displacement fields in Eq. (5.4) can be mapped onto rectangular ones as

	$\left\{\begin{aligned} &\sigma_{y}^{j}[z_{j}(\zeta_{j})]+\sigma_{x}^{j}[z_{j}(\zeta_{j})]=\sigma_{\theta}^{j}(\zeta_{j})+\sigma_{\rho}^{j}(\zeta_{j})\\ &\sigma_{y}^{j}[z_{j}(\zeta_{j})]-\sigma_{x}^{j}[z_{j}(\zeta_{j})]+2{\rm{i}}\tau_{xy}^{j}[z_{j}(\zeta_{j})]=\left[\sigma_{\theta}^{j}(\zeta_{j})-\sigma_{\rho}^{j}(\zeta_{j})+2{\rm{i}}\tau_{\rho\theta}^{j}(\zeta_{j})\right]\cdot\frac{\overline{\zeta}_{j}}{\zeta_{j}}\cdot\frac{\overline{z_{j}^{\prime}(\zeta_{j})}}{z_{j}^{\prime}(\zeta_{j})}\end{aligned}\right.,\quad\zeta_{j}=\rho_{j}\cdot\sigma_{j}$		(5.5a)
	$u_{j}[z_{j}(\zeta_{j})]+{\rm{i}}v_{j}[z_{j}(\zeta_{j})]=u_{j}(\zeta_{j})+v_{j}(\zeta_{j}),\quad\zeta_{j}=\rho_{j}\cdot\sigma_{j}$		(5.5b)

where $\sigma_{x}^{j}$ , $\sigma_{y}^{j}$ , and $\tau_{xy}^{j}$ denote horizontal, vertical, and shear components of the total stress field after $j^{\rm{th}}$ -stage excavation, respectively.

6 Numerical case and verification

The solution in Section 5 is analytical with infinite Laurent series of $f_{j,n}$ in Eq. (4.2), which should be bilaterally truncated into $2M_{j}+1$ items to reach numerical results. The coefficient series in Eq. (4.5) should be truncated as

A_{j,k}=\sum\limits_{n=-M_{j}}^{k}\alpha_{j,k+n}f_{j,-n}

(4.5a’)

B_{j,k}=\sum\limits_{n=k}^{M_{j}}\beta_{j,-k+n}f_{j,n}

(4.5b’)

where $j=1,2,3,\cdots,N$ . The infinite series to obtain solution of $f_{j,n}$ in Eqs. (4.2) should be truncated as well. Correspondingly, the stress and displacement fields in Section 5 should be truncated, which leads to numerical errors in oscillation form of Gibbs phenomena [32, 34]. To reduce the Gibbs phenomena, the Lanczos filtering is applied in Eq. (5.3) to replace $A_{j,k}$ and $B_{j,k}$ with $A_{j,k}\cdot L_{k}$ and $B_{j,k}\cdot L_{k}$ , respectively, where

L_{k}=\left\{\begin{aligned} &1,\quad k=0\\ &\sin\left(\frac{k}{M_{j}}\pi\right)/\left(\frac{k}{M_{j}}\pi\right),\quad{\rm{otherwise}}\end{aligned}\right.

(6.1)

with $-M_{j}\leq k\leq M_{j}$ .

In this section, a numerical case of 4 stage shallow tunnelling is performed to validate the present analytical solution by comparing with a corresponding finite element solution. The present analytical solution is conducted using the programming language Fortran of compiler GCC (version14.1.1). The linear systems are solved using the dgesv and zgesv packages of LAPACK (version 3.12.0). The figures are plotted using gnuplot (version 6.0).

6.1 Cavity geometries and bidirectional conformal mappings of sequential shallow tunnelling

The numerical case of sequential shallow tunnelling takes the geometry of the 4-stage excavation in Fig. 2b, while the sharp corner smoothening technique is applied to adapt the numerical schemes of bidirectional conformal mapping in A.2. Thus, the cavity boundaries of 4-stage excavation can be analytically expressed as

	${\rm{Stage\;1}}:\quad\left\{\begin{aligned} &x^{2}+(y+10)^{2}=5^{2},\quad-5\leq x\leq 0,-10\leq y\leq-5\\ &(x+4.5)^{2}+(y+10)^{2}=0.5^{2},\quad-5\leq x\leq-4.5,-10.5\leq y\leq-10\\ &y=-10.5,\quad-4.5\leq x\leq 0\\ &x^{2}+(y+10)^{2}=0.5^{2},\quad 0\leq x\leq 0.5,-10.5\leq y\leq-10\\ &x=0.5,\quad-10\leq y\leq-5.5\\ &x^{2}+(y+5.5)^{2}=0.5^{2},\quad 0\leq x\leq 0.5,-5.5\leq y\leq-5\\ \end{aligned}\right.$		(6.2a)
	${\rm{Stage\;2}}:\quad\left\{\begin{aligned} &x^{2}+(y+10)^{2}=5^{2},\quad-5\leq x\leq 5,-10\leq y\leq-5\\ &(x+4.5)^{2}+(y+10)^{2}=0.5^{2},\quad-5\leq x\leq-4.5,-10.5\leq y\leq-10\\ &y=-10.5,\quad-4.5\leq x\leq 4.5\\ &(x-4.5)^{2}+(y+10)^{2}=0.5^{2},\quad 4.5\leq x\leq 5,-10.5\leq y\leq-10\\ \end{aligned}\right.$		(6.2b)
	${\rm{Stage\;3}}:\quad\left\{\begin{aligned} &x^{2}+(y+10)^{2}=5^{2},\quad-5\leq x\leq 5,-10\leq y\leq-5\\ &x=-5,\quad-14.5\leq y\leq-10\\ &(x+4.5)^{2}+(y+14.5)^{2}=0.5^{2},\quad-5\leq x\leq-4.5,-15\leq y\leq-14.5\\ &y=-15,\quad-4.5\leq x\leq 0\\ &x^{2}+(y+14.5)^{2}=0.5^{2},\quad 0\leq x\leq 0.5,-15\leq y\leq-14.5\\ &x=0.5,\quad-14.5\leq y\leq-11\\ &(x-1)^{2}+(y+11)^{2}=0.5^{2},\quad 0.5\leq x\leq 1,-11\leq y\leq-10.5\\ &y=-10.5,\quad 1\leq x\leq 4.5\\ &(x-4.5)^{2}+(y+10)^{2}=0.5^{2},\quad 4.5\leq x\leq 5,-10.5\leq y\leq-10\\ \end{aligned}\right.$		(6.2c)
	${\rm{Stage\;4}}:\quad\left\{\begin{aligned} &x^{2}+(y+10)^{2}=5^{2},\quad-5\leq x\leq 5,-10\leq y\leq-5\\ &x=-5,\quad-14.5\leq y\leq-10\\ &(x+4.5)^{2}+(y+14.5)^{2}=0.5^{2},\quad-5\leq x\leq-4.5,-15\leq y\leq-14.5\\ &y=-15,\quad-4.5\leq x\leq 4.5\\ &(x-4.5)^{2}+(y+14.5)^{2}=0.5^{2},\quad 4.5\leq x\leq 5,-15\leq y\leq-14.5\\ &x=5,\quad-14.5\leq y\leq-10\\ \end{aligned}\right.$		(6.2d)

With the geometrical expressions in Eq. (6.2), the bidirectional conformal mappings of the 4-stage sequential shallow tunnelling can be obtained in Fig. 4, in which the collocation points are taken in the following pattern for good accuracy after trial computations: (1) 30 points are uniformly selected along the $90^{\circ}$ arc of $0.5{\rm{m}}$ radius; (2) 90 points are uniformly selected along the $90^{\circ}$ arc of $5{\rm{m}}$ radius; 60 points are uniformly selected along the rest straight lines. The collocation points along the cavity boundary in the physical plane $z$ are denoted by $z_{j,i}^{f}$ , where $j$ denotes the excavation stage, $i$ denotes the collocation point number, $f$ denotes the word "forward". After trial computations, the assignment factors take $K_{j}^{0}=k_{j}^{0}=2.2$ and $K_{j}^{1}=k_{j}^{1}=1.5$ , while $w_{j}^{c}=-2.5-7.5{\rm{i}}$ and $w_{j,0}=0.9\beta$ , where $\beta=5$ and $j=1,2,3,4$ .

To quantitatively determine the accuracy of the bidirectional conformal mapping, the following critieria are established as

\varepsilon_{j}=\max|z_{j,i}^{m,b}-z_{j,i}^{m,f}|,\quad i=1,2,3,\cdots,N_{j}

(6.3)

where $z_{j,i}^{m,f}$ denote the midpoints of the selected collocation points $z_{j,i}^{f}$ , and can be computed as

z_{j,i}^{m,f}=\frac{1}{2}(z_{j,i}^{f}+z_{j,i+1}^{f}),\quad i=1,2,3,\cdots,N_{j}

$z_{j,i}^{m,b}$ denote the coordinates of the corresponding midpoints $z_{j,i}^{m,f}$ after forward and backward conformal mappings, which can be computed in the following procedure.

Via Eq. (A.1a), the midpoints of the collocation points $z_{j,i}^{m,f}$ can be forwardly mapped onto corresponding mapping points $w_{j,i}^{m,f}$ in the interval mapping plane $\bm{\varOmega}_{j}^{w}$ . Then via Eq. (A.2a’), the midpoints in the mapping plane $\bm{\omega}_{j}$ can be computed as $\zeta_{j,i}^{m,f}$ . Subsequently, the midpoints for backward conformal mapping in the mapping plane $\bm{\omega}_{j}$ are given by only preserving the polar angles as

\zeta_{j,i}^{m,b}=\alpha_{j}\cdot\exp({\rm{i}}\arg\zeta_{j,i}^{m,f})

(6.4)

where $b$ in the superscript denotes "backward". Apparently, slight differences between $\zeta_{j,i}^{m,b}$ and $\zeta_{j,i}^{m,f}$ should exist. Whereafter, the backward midpoints in the interval mapping plane $\bm{\varOmega}_{j}^{w}$ can be computed as $w_{j,i}^{m,b}$ via Eq. (A.6). Finally, the backward midpoints $z_{j,i}^{m,b}$ in the physical plane $\bm{\varOmega}_{j}$ can be computed via Eq. (A.1b). Therefore, $\varepsilon_{j}$ in Eq. (4.3) should be able to present the maximum difference between the original midpoints and computed ones along cavity boundary for sequential excavation stages.

The orthogonal grids and corresponding values of $\varepsilon_{j}$ in Fig. 4 show that the bidirectional conformal mapping in A is very adaptive, and can be used to conformally map a lower half plane containing an asymmetrical cavity of very complicated shape (see Fig. 4c). Therefore, the conformal mapping can be used in further mechanical computation. However, we should also note in Fig. 4c that the orthogonal grid near concave boundary would be thin (as marked), indicating that the accuracy for a lower half plane containing a concave cavity would be slightly compromised.

6.2 Comparisons with finite element solution

To validate the present analytical solution in Section 5, a finite element solution using software ABAQUS 2020 is performed for comparisons. The mechanical parameters of geomaterial are listed in Table 1, where the coordinates of the joint points $T_{1}$ and $T_{2}$ are deliberately located near the shallow tunnel for better illustration of boundary conditions along ground surface (see Figs. 8 and 9).

Table 1: Mechanical parameters of numerical case

$\gamma{\rm{(kN/m^{3})}}$	$k_{x}$	$E{\rm{(MPa)}}$	$\nu$	$x_{0}$ (m)	$M_{j}$
$20$	$0.8$	$20$	$0.3$	$10$	$250$

The ${\rm{kN-m}}$ mechanical model of the finite element solution is shown in Fig. 5a, where the same 4-stage sequential shallow tunnelling is conducted. In Fig. 5a, the model geometry, geomaterial sketching, and seed distribution scheme outside of cavity are illustrated. In Fig. 5b, the 4-stage excavation procedure is marked by different colors, and more detailed seed distribution scheme of cavity boundaries is elaborated below: (1) 90 seeds are uniformly distributed along the $90^{\circ}$ arcs of $5{\rm{m}}$ ; (2) 30 seeds are uniformly distributed along the $90^{\circ}$ arcs of $0.5{\rm{m}}$ ; (3) 80 seeds are uniformly distributed along the rest straight lines. The meshing near cavity is shown in Fig. 5c, where the 4-stage excavation regions are marked with the same colors with the diagram of Fig. 5b. Both analytical and finite element solutions are in plane strain condition, thus, the 102499 finite elements use CPR4R type. The steps of finite element solution are listed in Table 2, where the 4 excavation stages are sequentially conducted, according to Fig. 5b.

Table 2: Steps of finite element solution

Step	Initial	Step 0	Step 1	Step 2	Step 3	Step 4
Procedure	(Initial)	Geostatic	Static,General	Static,General	Static,General	Static,General
Load	Applying constraints and geostress	Applying gravity	-	-	-	-
Interaction	-	-	Deactivating region I	Deactivating region II	Deactivating region III	Deactivating region IV

Substituting the mechanical parameters in Table 1 into the present solution and the finite element solution gives the results in Figs. 6-9. Fig. 6 shows the radial and tangential components of the residual stress along cavity boundaries of all four excavation stages, which are reduced by $10^{-2}$ for better illustration. In theory, the residual stress should be zero to accurately meet the zero boundary condition along cavity boundary (see Fig. 3). Fig. 6 indicates that both the present solution and the finite element solution agree to each other well for cavity boundaries with large curvature radii (a straight line has an infinite curvature radius). As for the rounded corners of small curvature radii, the present solution clearly surpasses the finite element solution. However, the residual stresses of $3^{\rm{rd}}$ stage excavation computed by the present solution contain significant errors near the concave boundary, indicating the defect of the present solution on concave cavity. Such a defect is caused by the bidirectional conformal mapping, and the curvilinear grid in Fig. 4c might provide certain clues. In Fig. 4c, the curvilinear grid density near the concave boundary is sparse and thin, which possibly causes loss of nearby mathematical information on the mechanical solution based on the curvilinear grid. Therefore, the present solution is more accurate than the finite element method for convex cavities.

Fig. 7 shows the Mises stress (reduced by $10^{-3}$ ) and deformation (magnified by $10^{2}$ ) along cavity boundaries for four excavation stages, and good agreements between these two solution are observed, except that the Mises stresses near the rounded corners obtained by the present solution are acuter and larger than those obtained by the finite element solution. The results of the present solution should be more accurate than the finite element solution, since the residual stresses near the rounded corners obtained by the finite element in Fig. 6 contains much larger errors.

Fig. 8 shows the comparisons of vertical and shear stress components along ground surface between the present solution and the finite element solution in all four excavation stages, and good agreements are observed. To be specific, the zero tractions in the range $x\in[-10,10]$ meet the requirement of the boundary condition in Eq. (2.6b), further indicating the correctness of the present solution.

Fig. 9 shows the comparisons of horizontal stress and ground deformation between the present solution and the finite element solution. The zero deformation in the range $x\in(\infty,-10]\cup[10,\infty)$ meets the requirement of the boundary condition in Eq. (2.6a), indicating that the present solution is capable of obtaining fixed far-field displacement as desired. The horizontal stress in the range $x\in(\infty,-10]\cup[10,\infty)$ between the present solution and the finite element solution is in good agreements. However, the horizontal stress and ground deformation in the range $x\in[-10,10]$ between present solution and finite element solution show clear discrepancy, especially for $1^{\rm{st}}$ , $2^{\rm{nd}}$ , and $3^{\rm{rd}}$ excavation stages. Noting that the overall indices $\varepsilon_{j}(j=1,2,3,4)$ are small with order of magnitutde $10^{-4}$ in Fig. 4, the bidirectional conformal mapping should be accurate, and such stress and ground deformation discrepancies are not caused by the mapping. The residual radial and tangential stress components in Fig. 6 (except for 6c), the vertical and shear stress components in the range $[-10,10]$ in Fig. 8, and the ground deformation in the range $(-\infty,-10]\cup[10,\infty)$ in Fig. 9 match the corresponding mixed boundary conditions in Eqs. (2.3), (2.6b), and (2.6a), respectively. By contrast, the residual radial and tangential stress components of the finite element solution in Fig. 6 contain significant errors, especially around the rounded corners. Therefore, it would be reasonable to presume that the results computed by the present solution should be more accurate than those computed by the finite element method. In other words, the present solution should be more accurate than the finite element solution. However, we should also notice that the obvious errors of residual radial and tangential stress components computed by the present solution in Fig. 6c, which are probably caused by the thin orthogonal grid near the concave. Such errors indicate that the present solution may be inaccurate for excavation with concave cavity.

In summary, Figs. 6-9 show the validation of the boundary conditions of the present solution, and the good agreements between the present solution and corresponding finite element solution. Additionally, detailed discussions reveal that the present solution should be more accurate than the finite element solution. The present solution can be used for further application.

7 Parametric investigation

In this section, several practical applications of present solution are illustrated to show its possible usage in sequential shallow tunnelling. To be consistent and simple, the benchmark geometry and sequential excavation stages take the same ones in Section 4(see Figs. 4 and 5b), and the benchmark parameters take the ones in Table 1.

7.1 Deformation along cavity boundary and solution convergence

Deformation along cavity boundary is significant in shallow tunnelling to estimate possible displacement around tunnel for further design optimization. In the validation of present solution in Section 5, the free ground surface is kept within the range of $x\in[-x_{0},x_{0}]$ of $x_{0}=10{\rm{m}}$ for good visualization and comparisons with finite element solution. However, the free ground surface should range wider in practical shallow tunnelling, thus, we select the joint coordinate as $x_{0}=10^{1},10^{2},10^{3},10^{4}{\rm{m}}$ to investigate the possible deformation along cavity boundary, while the rest parameters remain the same to the ones in Table 1.

The results of deformation along cavity boundary against joint point coordinate $x_{0}$ for four excavation stages are shown in Fig. 10 with a magnification of $10^{2}$ times for better illustration. Fig. 10 shows overall upheavals of geomaterial around tunnels and clear inward contractions against cavity boundaries for all four excavation stages, which are coincident to engineering expectations. Moreover, the deformed cavity boundaries of $x_{0}=10^{3}{\rm{m}}$ and $x_{0}=10^{4}{\rm{m}}$ overlap to each other, which indicates the same deformation for large free ground segment ${\bm{C}}_{0f}$ above tunnel, as well as convergence of present solution.

7.2 Mises stress along cavity boundary

Mises stress along cavity boundary can be mechanically interpreted as the hoop stress along cavity boundary, which is significant in shallow tunnelling to estimate possible damage zones around cavity for consideration of necessary reinforcement measures. In the validation, the lateral coefficient is set to be $k_{x}=0.8$ , however, the shallow strata to excavate are complicated with different lateral coefficients. Thus, we select lateral coefficient as $k_{x}=0.6,0.8,1.0,1.2,1.4,1.6$ to investigate the variation of Mises stress along cavity boundaries of sequential shallow tunnelling, and the rest parameters remain the same to the ones in Table 1.

The results of Mises stress along cavity boundaries against lateral coefficient $k_{x}$ for for excavation stages are shown in Fig. 11 with reduction of $10^{-3}$ times for better illustration. Fig. 11 shows that a larger lateral coefficient would cause higher stress concentrations near the rounded corners during sequential excavation, since the overall stress of the initial stress field increases. Therefore, tunnel corners may need more reinforcement measures in sequential excavation to avoid possible damage of geomaterial due to sress concentration. Moreover, Fig. 11c shows that the Mises stress near the right bottom geomaterial is very approximate to zero, indicating possible negative hoop stress. In other words, possible tensile stress may exist near the to-be-excavated region of $3^{\rm{rd}}$ -stage excavation, which may be hazardous in construction safety, and necessary monitoring of nearby geomaterial should be conducted.

7.3 Stress concentration of rounded corners

Fig. 11 shows remarkable stress concentrations near the corners during sequential shallow tunnelling, thus, the reduction of stress concentration near corners should be discussed due to its high value in shallow tunnelling safety. Among the four excavation stages in Fig. 11, the geometry of the $3^{\rm{th}}$ excavation stage is the most complicated, and worthy of discussion. The corners are the geometrical cause of stress concentration in Fig. 11, thus, we select different radii of the rounded corners (see the excavation scheme in Fig. 5b) as $0.3,0.4,0.5,0.6,0.7,0.8{\rm{m}}$ , and the rest cavity geometry is the same to that in Fig. 4c. The mechanical parameters take the same values in Table 1.

The Mises stresses along cavity boundary of different radii of rounded corners for the $3^{\rm{rd}}$ -stage excavation are shown in Fig. 12 with reduction of $10^{-3}$ times for better illustration. Fig. 12 shows that as the radius of rounded corners increases from $0.3{\rm{m}}$ to $0.8{\rm{m}}$ , the stress concentration decreases from $1960.11{\rm{kPa}}$ to $1247.92{\rm{kPa}}$ by a remarkable $36.33\%$ . Therefore, relatively large radius of rounded corners might be considered during sequential excavation to avoid overlarge stress concentration.

8 Further discussions

In exsiting complex variable solutions of shallow tunnelling in gravitational geomaterial, the cavity shapes are always fully circular [24, 25] or symmetrically noncircular [26, 31, 30], and no sequential shallow tunnelling is considered. The most fundamental reason is the mathematical difficulty to seek a suitable conformal mapping scheme that could bidirectionally map a lower half plane containing a noncircular and possibly asymmetrical cavity with shape changing due to sequential excavation. The conformal mapping scheme in A extends the bidirectional stepwise conformal mapping in our previous study [35], which incorporates Charge Simulation Method [37] and Complex Dipole Simulation Method [38], to a lower half plane containing more complicated cavity shapes, as shown in Fig. 4 (especially in Fig. 4c). Furthermore, the validation in Figs. 6-9 indicate the possibilities to extend the mathematical usage of the complex variable method from the full shallow tunnelling schemes with circular or symmetrically noncircular cavity shapes in the exsiting complex variable solutions [24, 25, 26, 31, 30] to sequential ones with complicated interval cavity shapes (see Figs. 6c and 7c).

Despite the theoretical improvements mentioned above, several defects of the present solution exist, and should be disclosed as:

(1) The most obvious defect of the present solution is that the present solution requires the rest geomaterial after sequential excavation to always remain doubly-connected to hold topological consistence of geomaterial and ensure solvability of the solution, while temporary supports are not considered for seeking mechanical variation within gravitational geomaterial alone. Temporary supports are important mechanical structure in sequential shallow tunnelling, and would alter the stress and displacement fields within rest geomaterial. Ignoring temporary supports would compromise the present solution to a certain extent. In our future studies, the complex variable method should be further improved and modified to be able to consider liners and temporary supports.

(1^′) On the other hand, however, the numerical cases in Sections 4 and 5, which simulate the multi-stepwise upper half vertical subdivision method in Fig. 2c without consideration of temporary supports, show that the present solution is capable of dealing with complicated interval cavity shapes of sequential shallow tunnelling. Then the present solution can thereby be degenerated to adapt more simple cavity shapes, for example, the top heading and bench method in Fig. 2b, where temporary supports may be not necessary. Therefore, the present solution is suitable for sequential shallow tunnelling methods that need no temporary supports without solution modification. The numerical cases in Sections 4 and 5, which certainly deviate from tunnel engineering fact without consideration of temporary supports, are deliberately conducted to demonstrate the significant theoretical improvements of present solution.

(2) Sequential excavation generally occurs with three-dimensional effect of tunnel face, which is not considered in the present solution. Thus, the stress and displacement fields computed by the present solution would deviate from real-world stress and displacement fields to a certain extent. To be conservative and to ensure construction safety, the present solution can be only used as very preliminary estimation of sequential shallow tunnelling.

(3) The thin orthogonal grid in Fig. 4c and residual stress errors in Fig. 6c indicate a possible weakness of the present solution that the stress and displacement for sequential excavation stage of shallow tunnels with concave cavity may be inaccurate. For computation and design of concave cavity, the present solution may be not suitable.

9 Conclusions

This paper proposes a new complex variable solution on asymmetrical sequential excavation of large-span shallow tunnels in gravitational geomaterial with consideration of rigid static equilibrium. The asymmetrical cavities of sequential tunnelling are conformally mapped using a new bidirectional stepwise mapping scheme consiting of Charge Simulation Method and Complex Dipole Simulation Method. The sequential shallow tunnelling procedure is mathematically disassembled into parallel and indepedent shallow tunnelling problems to seek similar pattern of mixed boundary conditions and solution procedure, which are subsequently transformed into a homogenerous Riemann-Hilbert problem to obtain reasonable stress and displacement field. Via comparisons with corresponding finite element solution, the bidirectional stepwise conformal mapping scheme and the proposed solution are both sufficiently validated. A final parametric investigation shows several possible applications of the proposed solution in asymmetrical sequential shallow tunnelling, and certain engineering recommendations are proposed. This new solution adapts lower half geomaterial containing asymmetrical cavity of complicated shape, and extends the usage of complex variable method in shallow tunnelling.

Acknowlegements

This study is financially supported by the Fujian Provincial Natural Science Foundation of China (Grant No. 2022J05190), the Scientific Research Foundation of Fujian University of Technology (Grant No. GY-H-22050), and the National Natural Science Foundation of China (Grant No. 52178318). The authors would like to thank Ph.D. Yiqun Huang for his suggestion on this study.

Author Contributions

Conceptualization: Luo-bin Lin, Fu-quan Chen; Methodology: Luo-bin Lin; Validation: Luo-bin Lin; Writing - original draft preparation: Luo-bin Lin; Writing - review and editing: Fu-quan Chen, Chang-jie Zheng, Shan-shun Lin; Funding acquisition: Chang-jie Zheng, Shang-shun Lin; Resources: Fu-quan Chen, Chang-jie Zheng; Supervision: Fu-quan Chen, Chang-jie Zheng, Shang-shun Lin

Declarations

The authors have no relevant financial or non-financial interests to disclose.

Data Availability

The research codes can be found via link github.com/luobinlin987/asymmetrical-sequential-static-equilibrium.

Appendix A Bidirectional conformal mapping

The geomaterials $\overline{\bm{\varOmega}}_{j}$ ( $j=1,2,3,\cdots,n$ ) after excavation are lower half planes containing corresponding cavities of arbitrary shapes ${\bm{D}}_{j}$ , which are doubly-connected regions and can be bidirectionally mapped onto unit annuli $\overline{\bm{\omega}}_{j}$ using the two-step mapping scheme in our previous study [35], respectively. Since the theory of such bidirectinal conformal mapping has been elaborated in detail in Ref [35], the mapping scheme below only illustrates necessary numerical procedure for conciseness. Comparing to existing mapping schemes [39, 26, 40], the proposed bidirectional conformal mapping has the following advantages: (1) Good adaptivity for a lower half plane containing an irregular and asymmetrical cavity; (2) Fast computation of solving a pair of linear systems below; (3) Closed logic to provide both forward and backward conformal mappings for mathematical rigor.

A.1 Mapping scheme

The first step is to bidirectionally map the geomaterial $\overline{\bm{\varOmega}}_{j}$ in the physical plane $z=x+{\rm{i}}y$ onto interval annuli $\overline{\bm{\varOmega}}_{j}^{w}$ in interval mapping plane $w_{j}=\Re w_{j}+{\rm{i}}\Im w_{j}$ using the following mapping schemes as

	$w_{j}=w_{j}(z)=\beta_{j}\frac{z-z_{c}^{j}}{z-\overline{z}_{c}^{j}},\quad z\in\overline{\bm{\varOmega}}_{j}$		(A.1a)
	$z=z_{j}(w_{j})=\frac{w_{j}\overline{z}_{c}^{j}-\beta_{j}z_{c}^{j}}{w_{j}-\beta_{j}},\quad w_{j}\in\overline{\bm{\varOmega}}_{j}^{w}$		(A.1b)

where $\beta_{j}$ denotes the exterior radius of the interval mapping annulus $\overline{\bm{\varOmega}}_{j}^{w}$ , which should be given before-hand.

Via Eq. (A.1), the finite free segment ${\bm{C}}_{0f,j}$ and far-field segment ${\bm{C}}_{0c,j}$ of the ground surface ${\bm{C}}_{0,j}$ in the physical plane $z=x+{\rm{i}}y$ are bidirectionally mapped onto free arc segment ${\bm{c}}_{0f,j}^{w}$ and constrained arc segment ${\bm{c}}_{0c,j}^{w}$ of the exterior boundary ${\bm{c}}_{0,j}^{w}={\bm{c}}_{0f,j}^{w}\cup{\bm{c}}_{0c,j}^{w}$ in the interval mapping plane $w_{j}=\Re w_{j}+{\rm{i}}\Im w_{j}$ , respectively. The joint points $T_{1}$ and $T_{2}$ connecting segments ${\bm{C}}_{0f,j}$ and ${\bm{C}}_{0c,j}$ in the physical plane $z=x+{\rm{i}}y$ are bidirectionally mapped onto corresponding joint points $t_{1,j}^{w}$ and $t_{2,j}^{w}$ connecting arc segments ${\bm{c}}_{0f,j}^{w}$ and ${\bm{c}}_{0c,j}^{w}$ in the interval mapping plane $w_{j}=\Re w_{j}+{\rm{i}}\Im w_{j}$ , respectively. The cavity boundary ${\bm{C}}_{j}$ in the physical plane $z=x+{\rm{i}}y$ is bidirectionally mapped onto ${\bm{c}}_{j}^{w}$ in the interval mapping plane $w_{j}=\Re w_{j}+{\rm{i}}\Im w_{j}$ . Figs. A-1a and A-1b graphically and schematically show the bidirectional conformal mapping in Eq. (A.1) for the example in Fig. 3b.

The second step is to bidirectionally map the interval annulus $\overline{\bm{\varOmega}}_{j}^{w}$ in the interval mapping plane $w_{j}=\Re w_{j}+{\rm{i}}\Im w_{j}$ onto corresponding unit mapping annulus $\overline{\bm{\omega}}_{j}$ of interior radius $\alpha_{j}$ in the mapping plane $\zeta_{j}=\rho_{j}\cdot{\rm{e}}^{{\rm{i}}\theta_{j}}$ . The forward and backward conformal mappings of second-step bidirectional conformal mapping can be respectively expressed using Charge Simulation Method [37, 41, 35] and Complex Dipole Simulation Method [38, 35] as

	$\zeta_{j}=\zeta_{j}(w_{j})=\frac{w_{j}-w_{j}^{c}}{w_{j}^{\beta}-w_{j}^{c}}\cdot\exp\left[\sum\limits_{k=1}^{N_{0}}P_{j,k}\ln\frac{w_{j}-U_{j,k}}{w_{j}^{\beta}-U_{j,k}}+\sum\limits_{k=1}^{N_{j}}Q_{j,k}\ln\frac{w_{j}-V_{j,k}}{w_{j}^{\beta}-V_{j,k}}\right]$		(A.2a)
	$w_{j}=w_{j}(\zeta_{j})=\sum\limits_{k=1}^{N_{0}}\frac{p_{j,k}}{\zeta_{j}-\eta_{j,k}}+\sum\limits_{k=1}^{N_{j}}\frac{q_{j,k}}{\zeta_{j}-\mu_{j,k}}$		(A.2b)

where $w_{j}^{c}$ denotes arbitrary point within interior boundary ${\bm{c}}_{j}^{w}$ of the interval mapping annulus $\overline{\bm{\varOmega}}_{j}^{w}$ for fixed point normalization; $w_{j}^{\beta}$ denotes arbitrary point along exterior bounadry ${\bm{c}}_{0,j}$ (the radius is $\beta_{j}$ ) for rotation normalization, and generally can take the coordinate $w_{j}^{\beta}=\beta$ for simplicity; $P_{j,k}$ and $Q_{j,k}$ denote charges of Charge Simulation Method, and are real variables to be determined by Eq. (A.5); $U_{j,k}$ and $V_{j,k}$ denote charge points of Charge Simulation Method, and are complex coefficients determined by Eq. (A.3); $p_{j,k}$ and $q_{j,k}$ are charges of Complex Dipole Simulation Method, and are complex variables to be determined by Eq. (A.6); $\eta_{j,k}$ and $\mu_{j,k}$ are charge points of Complex Dipole Simulation Method, and are complex coefficients determined by Eq. (A.4); $N_{0}$ and $N_{j}$ denote quantities of charge points along exterior boundary ${\bm{c}}_{0,j}^{w}$ (or ${\bm{c}}_{0,j}$ ) and interior boundary ${\bm{c}}_{j}^{w}$ (or ${\bm{c}}_{j}$ ), respectively.

Via Eq. (A.2), the free arc segment ${\bm{c}}_{0f,j}^{w}$ and constrained arc segment ${\bm{c}}_{0c,j}^{w}$ of exterior boundary ${\bm{c}}_{0f}^{w}$ in the interval mapping plane $w_{j}=\Re w_{j}+{\rm{i}}\Im w_{j}$ are bidirectinally mapped onto free segment ${\bm{c}}_{0f,j}$ and constrained arc segment ${\bm{c}}_{0c,j}$ of unit exterior boundary ${\bm{c}}_{0,j}={\bm{c}}_{0f,j}\cup{\bm{c}}_{0c,j}$ in the mapping plane $\zeta_{j}=\rho_{j}\cdot{\rm{e}}^{{\rm{i}}\theta_{j}}$ , respectively. The joint points $t_{1,j}^{w}$ and $t_{2,j}^{w}$ connecting arc segments ${\bm{c}}_{0f,j}^{w}$ and ${\bm{c}}_{0c,j}^{w}$ in the interval mapping plane $w=\Re w+{\rm{i}}\Im w$ are bidirectionally mapped onto corresponding joint points $t_{1,j}$ and $t_{2,j}$ connecting unit arc segments ${\bm{c}}_{0f,j}$ and ${\bm{c}}_{0c,j}$ in the mapping plane $\zeta_{j}=\rho_{j}\cdot{\rm{e}}^{{\rm{i}}\theta_{j}}$ , respectively. The possibly noncircular interior boundary ${\bm{c}}_{j}^{w}$ in the interval mapping plane $w=\Re w+{\rm{i}}\Im w$ is bidirectionally mapped onto the circular interior boundary ${\bm{c}}_{j}$ of radius $\alpha_{j}$ . The interior and exterior regions of boundary ${\bm{c}}_{0,j}$ are denoted by ${\bm{\omega}}_{0,j}^{+}$ and ${\bm{\omega}}_{0,j}^{-}$ , respectively. Figs. A-1b and A-1c graphically and schematically show the subsequent bidirectional conformal mapping in Eq. (A.2) for the example in Fig. A-1a.

The charge points of Charge Simulation Method can be given as

	$\left\{\begin{aligned} &U_{j,k}=w_{j,k}^{0}+K_{j}^{0}\cdot H_{j,k}^{0}\cdot\exp\left({\rm{i}}\varTheta_{j,k}^{0}\right)\\ &H_{j,k}^{0}=\frac{1}{2}\left(\left\|w_{j,k+1}^{0}-w_{j,k}^{0}\right\|+\left\|w_{j,k}^{0}-w_{j,k-1}^{0}\right\|\right)\\ &\varTheta_{j,k}^{0}=\arg\left(w_{j,k+1}^{0}-w_{j,k-1}^{0}\right)-\frac{\pi}{2}\end{aligned}\right.,\quad k=1,2,3,\cdots,N_{0}$		(A.3a)
	$\left\{\begin{aligned} &V_{j,k}=w_{j,k}^{1}+K_{j}^{1}\cdot H_{j,k}^{1}\cdot\exp\left({\rm{i}}\varTheta_{j,k}^{1}\right)\\ &H_{j,k}^{1}=\frac{1}{2}\left(\left\|w_{j,k+1}^{1}-w_{j,k}^{1}\right\|+\left\|w_{j,k}^{1}-w_{j,k-1}^{1}\right\|\right)\\ &\varTheta_{j,k}^{1}=\arg\left(w_{j,k+1}^{1}-w_{j,k-1}^{1}\right)-\frac{\pi}{2}\end{aligned}\right.,\quad k=1,2,3,\cdots,N_{j}$		(A.3b)

where $w_{j,k}^{0}$ and $w_{j,k}^{1}$ denote collocation points of Charge Simulation Method distributed along exterior boundary ${\bm{c}}_{0,j}^{w}$ and interior boundary ${\bm{c}}_{j}^{w}$ of the interval mapping region $\overline{\bm{\varOmega}}_{j}^{w}$ , respectively, and the selection schemes can be seen in A.2; $K_{j}^{0}$ and $K_{j}^{1}$ denote assignment factors of Charge Simulation Method for collocation points along exterior boundary ${\bm{c}}_{0,j}^{w}$ and interior boundary ${\bm{c}}_{j}^{w}$ , respectively; $H_{j,k}^{0}$ and $H_{j,k}^{1}$ denote distances between charge points and corresponding collocation points along exterior boundary ${\bm{c}}_{0,j}^{w}$ and interior boundary ${\bm{c}}_{j}^{w}$ , respectively; $\varTheta_{j,k}^{0}$ and $\varTheta_{j,k}^{1}$ denote outward normal angles pointing from collocation points to corresponding charge points along exterior boundary ${\bm{c}}_{0,j}^{w}$ and interior boundary ${\bm{c}}_{j}^{w}$ , respectively.

The charge points of Complex Dipole Simulation Method can be given as

	$\left\{\begin{aligned} &\eta_{j,k}=\zeta_{j,k}^{0}+k_{j}^{0}\cdot h_{j,k}^{0}\cdot\exp\left({\rm{i}}\vartheta_{j,k}^{0}\right)\\ &h_{j,k}^{0}=\frac{1}{2}\left(\left\|\zeta_{j,k+1}^{0}-\zeta_{j,k}^{0}\right\|+\left\|\zeta_{j,k}^{0}-\zeta_{j,k-1}^{0}\right\|\right)\\ &\vartheta_{j,k}^{0}=\arg\left(\zeta_{j,k+1}^{0}-\zeta_{j,k-1}^{0}\right)-\frac{\pi}{2}\end{aligned}\right.,\quad k=1,2,3,\cdots,N_{0}$		(A.4a)
	$\left\{\begin{aligned} &\mu_{j,k}=\zeta_{j,k}^{1}+k_{j}^{1}\cdot h_{j,k}^{1}\cdot\exp\left({\rm{i}}\vartheta_{j,k}^{1}\right)\\ &h_{j,k}^{1}=\frac{1}{2}\left(\left\|\zeta_{j,k+1}^{1}-\zeta_{j,k}^{1}\right\|+\left\|\zeta_{j,k}^{1}-\zeta_{j,k-1}^{1}\right\|\right)\\ &\vartheta_{j,k}^{1}=\arg\left(\zeta_{j,k+1}^{1}-\zeta_{j,k-1}^{1}\right)-\frac{\pi}{2}\end{aligned}\right.,\quad k=1,2,3,\cdots,N_{j}$		(A.4b)

where $\zeta_{j,k}^{0}$ and $\zeta_{j,k}^{1}$ denote collocation points of Complex Dipole Simulation Method distributed along exterior boundary ${\bm{c}}_{0,j}$ and interior boundary ${\bm{c}}_{j}$ of the interval mapping region $\overline{\bm{\omega}}_{j}$ , respectively, and are determined by Eq. (A.7); $k_{j}^{0}$ and $k_{j}^{1}$ denote assignment factors of Complex Dipole Simulation Method for collocation points along exterior boundary ${\bm{c}}_{0,j}$ and interior boundary ${\bm{c}}_{j}$ , respectively; $h_{j,k}^{0}$ and $h_{j,k}^{1}$ denote distances between charge points and corresponding collocation points along exterior boundary ${\bm{c}}_{0,j}$ and interior boundary ${\bm{c}}_{j}$ , respectively; $\vartheta_{j,k}^{0}$ and $\vartheta_{j,k}^{1}$ denote outward normal angles pointing from collocation points to corresponding charge points along exterior boundary ${\bm{c}}_{0,j}$ and interior boundary ${\bm{c}}_{j}$ , respectively.

With the charge points in Eq. (A.3), the charges in Eq. (A.2a) can be determined as

	$\sum\limits_{k=1}^{N_{0}}\ln\left\|\frac{w_{j,i}^{0}-U_{j,k}}{w_{j}^{\beta}-U_{j,k}}\right\|\cdot P_{j,k}+\sum\limits_{k=1}^{N_{j}}\ln\left\|\frac{w_{j,i}^{0}-V_{j,k}}{w_{j}^{\beta}-V_{j,k}}\right\|\cdot Q_{j,k}-\ln r_{e}=-\ln\left\|\frac{w_{j,k}^{0}-w_{j}^{c}}{w_{j}^{\beta}-w_{j}^{c}}\right\|,\quad i=1,2,3,\cdots,N_{0}$		(A.5a)
	$\sum\limits_{k=1}^{N_{0}}\ln\left\|\frac{w_{j,i}^{1}-U_{j,k}}{w_{j}^{\beta}-U_{j,k}}\right\|\cdot P_{j,k}+\sum\limits_{k=1}^{N_{j}}\ln\left\|\frac{w_{j,i}^{1}-V_{j,k}}{w_{j}^{\beta}-V_{j,k}}\right\|\cdot Q_{j,k}-\ln\alpha_{j}=-\ln\left\|\frac{w_{j,k}^{1}-w_{j}^{c}}{w_{j}^{\beta}-w_{j}^{c}}\right\|,\quad i=1,2,3,\cdots,N_{j}$		(A.5b)
	$\sum\limits_{k=1}^{N_{0}}P_{j,k}=-1$		(A.5c)
	$\sum\limits_{k=1}^{N_{j}}Q_{j,k}=0$		(A.5d)

Eq. (A.5) contains $N_{0}+N_{j}+2$ to-be-determined real variables and real linear equations to form simultaneous real linear system, and unique solution can be obtained.

With the charge points in Eq. (A.4), the charges in Eq. (A.2b) can be determined as

	$\sum\limits_{k=1}^{N_{0}}\frac{p_{j,k}}{\zeta_{j,i}^{0}-\eta_{j,k}}+\sum\limits_{k=1}^{N_{j}}\frac{q_{j,k}}{\zeta_{j,i}^{0}-\mu_{j,k}}=w_{j,i}^{0},\quad i=1,2,3,\cdots,N_{0}$		(A.6a)
	$\sum\limits_{k=1}^{N_{0}}\frac{p_{j,k}}{\zeta_{j,i}^{1}-\eta_{j,k}}+\sum\limits_{k=1}^{N_{j}}\frac{q_{j,k}}{\zeta_{j,i}^{1}-\mu_{j,k}}=w_{j,i}^{1},\quad i=1,2,3,\cdots,N_{j}$		(A.6b)

where

	$\zeta_{j,i}^{0}=\zeta_{j}(w_{j,i}^{0}),\quad i=1,2,3,\cdots,N_{0}$		(A.7a)
	$\zeta_{j,i}^{1}=\zeta_{j}(w_{j,i}^{1}),\quad i=1,2,3,\cdots,N_{j}$		(A.7b)

Eq. (A.7) is computed using Eq. (A.2a’). Eq. (A.6) contains $N_{0}+N_{j}$ to-be-determined complex variables and complex linear equations to form simultaneous complex linear system, and unique solution can be obtained.

In summary, the bidirectional conformal mapping in Eq. (A.2) can be determined by the linear systems in Eqs. (A.5) and (A.6) with auxiliary supports of Eqs. (A.3), (A.4), and (A.7), except for the before-hand determination of collocation points $w_{j,k}^{0}$ and $w_{j,k}^{1}$ , which are given in A.2 with certain numerical principles.

A.2 Numerical schemes

The collocation point selection of $w_{j,k}^{0}$ and $w_{j,k}^{1}$ in Eq. (A.3) should be suitable to guarantee the analyticity and accuracy of the backward conformal mapping in Eq. (A.2b). The collocation points along boundaries ${\bm{c}}_{0,j}^{w}$ and ${\bm{c}}_{j}^{w}$ of the interval mapping annulus $\overline{\bm{\varOmega}}_{j}^{w}$ should be distributed as uniformly as possible, and sharp corners should be eliminated to avoid boundary singularities [42, 43]. Therefore, we should first ensure that no sharp corners exist along the cavity boundary ${\bm{C}}_{j}$ in the physical geomaterial region $\overline{\bm{\varOmega}}_{j}$ . In Figs. 2 and 3 and A-1a, the sharp corners of cavity boundaries for sequential shallow tunnelling are kept merely due to convenience of schematic plotting. In practical tunnel engineering, sharp corners are generally rounded on purpose to avoid stress concentration. Therefore, the sharp corners in this paper do not exist, and Fig. A-2 graphically shows the replacement of rounded corners with sharp corners for the example of Fig. A-1a.

Therefore, the rounded corners guarantee the continuity and analyticity of the cavity boundary ${\bm{C}}_{j}$ in the physical plane $z=x+{\rm{i}}y$ . Subsequently, via the first-step conformal mapping in Eq. (A.1), the exterior circular boundary ${\bm{c}}_{0,j}^{w}$ (mapped from the straight ground surface ${\bm{C}}_{0}$ ) and interior boundary ${\bm{c}}_{j}^{w}$ in the inteval mapping annulus $\overline{\bm{\varOmega}}_{j}^{w}$ (possibly not circular) would contain no sharp corners, so that the continuity of the outward normal angles $\varTheta_{j,k}^{0}$ and $\varTheta_{j,k}^{1}$ in Eq. (A.3) would be potentially guaranteed by mandatorily setting the following numerical principle as

	$\max\left[\max\limits_{1\leq k\leq N_{0}}\left\|\arg\left(w_{j,k+1}^{0}-w_{j,k}^{0}\right)\right\|,\max\limits_{1\leq k\leq N_{j}}\left\|\arg\left(w_{j,k+1}^{1}-w_{j,k}^{1}\right)\right\|\right]\leq 10^{\circ}$	(A.8a)
Moreover, the quantity of collocation points should be large enough, thus, the following numerical principle is mandatorily set as
	$\max\left[\max\limits_{1\leq k\leq N_{0}}\frac{\|w_{j,k+1}^{0}-w_{j,k}^{0}\|}{\sum\limits_{l=1}^{N_{0}}\|w_{j,l}^{0}-w_{j,l}^{0}\|},\max\limits_{1\leq k\leq N_{j}}\frac{\|w_{j,k+1}^{1}-w_{j,k}^{1}\|}{\sum\limits_{l=1}^{N_{j}}\|w_{j,l}^{1}-w_{j,l}^{1}\|}\right]\leq 10^{-2}$	(A.8b)

Eq. (A.8) is based on the numerical principles proposed in our previous study [35], and controls the varying angle and distance between adjacent collocation points along boundaries ${\bm{c}}_{0,j}^{w}$ and ${\bm{c}}_{j}^{w}$ .

Eq. (A.2a) should be slightly modified after solving the charges by Eq. (A.5) to guarantee the single-valuedness of arguments without changing the principal values of the computation results [41, 35] as

\zeta_{j}=\zeta_{j}(w_{j})=\frac{w_{j}-w_{j}^{c}}{w_{j}^{\beta}-w_{j}^{c}}\cdot\exp\left[\begin{aligned} &\sum\limits_{k=1}^{N_{0}}P_{j,k}\left(\ln\frac{w_{j}-U_{j,k}}{w_{j,0}-U_{j,k}}-\ln\frac{w_{j}^{\beta}-U_{j,k}}{w_{j,0}-U_{j,k}}\right)\\ +&\sum\limits_{k=1}^{N_{j}}Q_{j,k}\left(\ln\frac{w_{j}-V_{j,k}}{w_{j}-w_{j}^{c}}-\ln\frac{w_{j}^{\beta}-V_{j,k}}{w_{j}^{\beta}-w_{j}^{c}}\right)\end{aligned}\right]

(A.2a’)

where $w_{j,0}$ denote an arbitrary point within interval mapping region $\overline{\bm{\varOmega}}_{j}^{w}$ .

References

\bibcommenthead
Sharifzadeh et al. [2013] Sharifzadeh, M., Kolivand, F., Ghorbani, M., Yasrobi, S.: Design of sequential excavation method for large span urban tunnels in soft ground–niayesh tunnel. Tunn. Undergr. Space Technol. 35, 178–188 (2013)
Li et al. [2016] Li, P., Zhao, Y., Zhou, X.: Displacement characteristics of high-speed railway tunnel construction in loess ground by using multi-step excavation method. Tunn. Undergr. Space Technol. 51, 41–55 (2016)
Fang et al. [2017] Fang, Q., Liu, X., Zhang, D., Lou, H.: Shallow tunnel construction with irregular surface topography using cross diaphragm method. Tunn. Undergr. Space Technol. 68, 11–21 (2017)
Shi et al. [2017] Shi, C., Cao, C., Lei, M.: Construction technology for a shallow-buried underwater interchange tunnel with a large span. Tunn. Undergr. Space Technol. 70, 317–329 (2017)
Cao et al. [2018] Cao, C., Shi, C., Lei, M., Yang, W., Liu, J.: Squeezing failure of tunnels: a case study. Tunn. Undergr. Space Technol. 77, 188–203 (2018)
Liu et al. [2021] Liu, W., Chen, J., Luo, Y., Chen, L., Shi, Z., Wu, Y.: Deformation behaviors and mechanical mechanisms of double primary linings for large-span tunnels in squeezing rock: a case study. Rock Mech. Rock Eng. 54, 2291–2310 (2021)
Zhou et al. [2021] Zhou, S., Li, L., An, Z., Liu, H., Yang, G., Zhou, P.: Stress-release law and deformation characteristics of large-span tunnel excavated with semi central diaphragm method. KSCE J. Civil Eng. 25, 2275–2284 (2021)
Verruijt [1997a] Verruijt, A.: Deformations of an elastic plane with a circular cavity. Int. J. Solids Struct. 35(21), 2795–2804 (1997)
Verruijt [1997b] Verruijt, A.: A complex variable solution for a deforming circular tunnel in an elastic half-plane. Int. J. Numer. Anal. Methods Geomech. 21(2), 77–89 (1997)
Zhang et al. [2018] Zhang, Z., Huang, M., Xi, X., Yang, X.: Complex variable solutions for soil and liner deformation due to tunneling in clays. Int. J. Geomech. 18(7), 04018074 (2018)
Wang et al. [2020] Wang, H., Gao, X., Wu, L., Jiang, M.: Analytical study on interaction between existing and new tunnels parallel excavated in semi-infinite viscoelastic ground. Comput. Geotech. 120, 103385 (2020)
Zhang et al. [2020] Zhang, Z., Pan, Y., Zhang, M., Lv, X., Jiang, K., Li, S.: Complex variable analytical prediction for ground deformation and lining responses due to shield tunneling considering groundwater level variation in clays. Comput. Geotech. 120, 103443 (2020)
Zhang et al. [2021] 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. Comput. Geotech. 129, 103858 (2021)
Kong et al. [2019] Kong, F., Lu, D., Du, X., Shen, C.: Displacement analytical prediction of shallow tunnel based on unified displacement function under slope boundary. Int. J. Numer. Anal. Methods Geomech. 43(1), 183–211 (2019)
Lu et al. [2019] Lu, D., Kong, F., Du, X., Shen, C., Gong, Q., Li, P.: A unified displacement function to analytically predict ground deformation of shallow tunnel. Tunn. Undergr. Space Technol. 88, 129–143 (2019)
Kong et al. [2021] Kong, F., Lu, D., Du, X., Li, X., Su, C.: Analytical solution of stress and displacement for a circular underwater shallow tunnel based on a unified stress function. Ocean Eng. 219, 108352 (2021)
Wang et al. [2018a] Wang, H., Wu, L., Jiang, M., Song, F.: Analytical stress and displacement due to twin tunneling in an elastic semi-infinite ground subjected to surcharge loads. Int. J. Numer. Anal. Methods Geomech. 42(6), 809–828 (2018)
Wang et al. [2018b] Wang, H., Chen, X.P., Jiang, M., Song, F., Wu, L.: The analytical predictions on displacement and stress around shallow tunnels subjected to surcharge loadings. Tunn. Undergr. Space Technol. 71, 403–427 (2018)
Lin et al. [2021] Lin, L., Lu, Y., Chen, F., Li, D., Zhang, B.: Analytic study of stress and displacement for shallow twin tunnels subjected to surcharge loads. Appl. Math. Model. 93, 485–508 (2021)
Strack [2002] Strack, O.E.: Analytic solutions of elastic tunneling problems. PhD thesis, Delft University of Technology, Amsterdam (2002)
Strack and Verruijt [2002] Strack, O.E., Verruijt, A.: A complex variable solution for a deforming buoyant tunnel in a heavy elastic half-plane. Int. J. Numer. Anal. Methods Geomech. 26(12), 1235–1252 (2002)
Verruijt and Strack [2008] Verruijt, A., Strack, O.E.: Buoyancy of tunnels in soft soils. Géotech. 58(6), 513–515 (2008)
Fang et al. [2015] Fang, Q., Song, H., Zhang, D.: Complex variable analysis for stress distribution of an underwater tunnel in an elastic half plane. Int. J. Numer. Anal. Methods Geomech. 39(16), 1821–1835 (2015)
Lu et al. [2016] Lu, A., Zeng, X., Xu, Z.: Solution for a circular cavity in an elastic half plane under gravity and arbitrary lateral stress. Int. J. Rock Mech. Min. Sci. 89, 34–42 (2016)
Lu et al. [2019] Lu, A., Cai, H., Wang, S.: A new analytical approach for a shallow circular hydraulic tunnel. Mecc. 54(1-2), 223–238 (2019)
Zeng et al. [2019] Zeng, G., Cai, H., Lu, A.: An analytical solution for an arbitrary cavity in an elastic half-plane. Rock Mech. Rock Eng. 52, 4509–4526 (2019)
Lu et al. [2021] Lu, A., Zeng, G., Zhang, N.: A complex variable solution for a non-circular tunnel in an elastic half-plane. Int. J. Numer. Anal. Methods Geomech. 45(12), 1833–1853 (2021)
Zeng et al. [2022] Zeng, G., Wang, H., Jiang, M.: Analytical solutions of noncircular tunnels in viscoelastic semi-infinite ground subjected to surcharge loadings. Appl. Math. Model. 102, 492–510 (2022)
Zeng et al. [2023] Zeng, G., Wang, H., Jiang, M.: Analytical stress and displacement of twin noncircular tunnels in elastic semi-infinite ground. Comput. Geotech. 160, 105520 (2023)
Fan et al. [2024] Fan, T., Wang, J., Wang, Z.: Analytical solutions of elastic complex variables for tunnels with complicated shapes under geostress field. Rock Mech. Rock Eng., 1–21 (2024)
Zhou et al. [2024] Zhou, Y., Lu, A., Zhang, N.: An analytical solution for the orthotropic semi-infinite plane with an arbitrary-shaped hole. Math. Mech. Solids, 10812865231225131 (2024)
Lin et al. [2024a] Lin, L., Chen, F., Huang, X.: Reasonable mechanical model on shallow tunnel excavation to eliminate displacement singularity caused by unbalanced resultant. Appl. Math. Model. 127, 396–431 (2024)
Lin et al. [2024b] Lin, L., Chen, F., Lin, S.: A new complex variable solution on noncircular shallow tunnelling with reasonable far-field displacement. Comput. Geotech. 166, 106008 (2024)
Lin et al. [2024c] Lin, L., Chen, F., Zhuang, J.: Complex variable solution on over-/under-break shallow tunnelling in gravitational geomaterial with reasonable far-field displacement. Comput. Geotech. 174, 106595 (2024)
Lin et al. [2024d] Lin, L., Chen, F., Zheng, C., Lin, S.: Bidirectional conformal mapping for over-break and under-break tunnelling and its application in complex variable method. arXiv preprint arXiv:2406.12148 (2024)
Muskhelishvili [1966] Muskhelishvili, N.I.: Some Basic Problems of the Mathematical Theory of Elasticity, 4th edn. Cambridge University Press, Cambridge (1966)
Amano [1994] Amano, K.: A charge simulation method for the numerical conformal mapping of interior, exterior and doubly-connected domains. J. Comput. Appl. Math. 53(3), 353–370 (1994)
Sakakibara [2020] Sakakibara, K.: Bidirectional numerical conformal mapping based on the dipole simulation method. Eng. Anal. Bound. Elem. 114, 45–57 (2020)
Exadaktylos and Stavropoulou [2002] Exadaktylos, G., Stavropoulou, M.: A closed-form elastic solution for stresses and displacements around tunnels. Int. J. Rock Mech. Min. Sci. 39(7), 905–916 (2002)
Ye and Ai [2023] Ye, Z., Ai, Z.: A novel complex variable solution for the stress and displacement fields around a shallow non-circular tunnel. Comput. Geotech. 164, 105797 (2023)
Okano et al. [2003] Okano, D., Ogata, H., Amano, K., Sugihara, M.: Numerical conformal mappings of bounded multiply connected domains by the charge simulation method. J. Comput. Appl. Math. 159(1), 109–117 (2003)
Hough and Papamichael [1981] Hough, D.M., Papamichael, N.: The use of splines and singular functions in an integral equation method for conformal mapping. Numer. Math. 37, 133–147 (1981)
Hough and Papamichael [1983] Hough, D., Papamichael, N.: An integral equation method for the numerical conformal mapping of interior, exterior and doubly-connected domains. Numer. Math. 41, 287–307 (1983)