Introduction

The construction of large-section tunnels inevitably encounters shallow cover, unsymmetrical pressure, and adverse geological conditions at the entrances and exits. It poses a great challenge to construction safety, which requires appropriate excavation methods1,2,3,4,5. The Sequential Excavation Method (SEM) is the tunnel construction method that alleviates the pressure of surroundings through staged excavation, using permanent and temporary support structures to achieve self-support and balance stability6,8.The Center Diaphragm Method (CDM), as a form within SEM, divides the excavation cross section into smaller parts through the center diaphragm, improving the safety and stability of construction under poor geological conditions. Global underground engineering practices and research have confirmed that CDM is the way to reduce the construction risk for shallow large-span tunnels8,9,10,11.

Currently, research on the constructional behavior of CDM is conducted primarily through (1) finite element simulation12,13,14,15,16,17,18,19,20,21, (2) laboratory and field monitoring12,13,14,15,16,22, and (3) mechanical model23. Zhou et al.12 used finite element simulation and field monitoring to analyze deformation and stress relief characteristics in large-section tunnel excavation using the upper bench CDM. GUO et al.22 analyzed field monitoring data to conclude that full-section grouting assisted CDM construction significantly reduces surface settlement and enhances construction quality. Luo et al.23 based on upper bench CDM construction of shallow, large-span tunnels, developed mechanical models for different construction stages to back-calculate rock pressure, revealing the evolution pattern of load release. Remarkably, the most crucial component in CDM is the center diaphragm (Fig. 1), which serves as a temporary support. Its significance has been increasingly acknowledged and extensive studies have been conducted. Luo et al. combining field monitoring and numerical simulation results, identified deformation patterns in the horizontal and vertical directions of the center diaphragm and recommended using longitudinal reinforcement to connect steel frames and reinforcement nets based on forces acting on the center diaphragm24. For large-section shallowly buried tunnels excavated using the Center Cross Diaphragm Method(CRDM), the center diaphragm was considered a small-curvature beam, and a mechanical model under the horizontal pressure of surroundings and upper structure loads was established25. Based on the upper bench CDM excavation of large-span loess tunnels, a mechanical model was established in which the primary lining steel frames and the center diaphragm structure jointly bear the load and coordinate deformation26. Zhang et al.27 proposed the method of sectional construction with the combined support of steel frames and rock walls. It was discovered that during the excavation stage of the upper surroundings of the tunnel, the tunnel support structure is at the highest risk of being damaged. Among them, the center diaphragm is the most vulnerable, while the forces that act on the remaining parts are hardly influenced by the construction of the tunnel. Zhao28 analyzed field monitoring data to uncover the mechanical properties and displacement trends of support structures in large-section highway tunnels constructed using half-bench CDM, utilizing a partial excavation support structure mechanical response model for analysis. Hou et al.29 conducted dynamic damage experiments on center diaphragms under blast loads, identifying failure modes and types of center diaphragms. Wang et al.30 utilized numerical simulation software to determine the impact and deformation characteristics of temporary steel supports in shallow tunnel excavation with unsymmetrical pressure using CDM. However, despite being a key temporary support structure, center diaphragms are often designed empirically, and scholarly research on the impact of different center diaphragm forms within CDM is still limited. Although Zhou et al.31 based their study on the traditional CDM and modified the temporary support form into vertical diaphragms for comparative analysis with conventional CDM, the construction impacts of center diaphragms with different curvatures are not analyzed and discussed sufficiently.

Fig. 1
figure 1

Tunnel construction site adopting CDM.

Full size image

This study proposed a mechanical model of CDM. By establishing a double-arch support structure model of the primary lining-upper center diaphragm and its coordinate system, with the covering depth, unsymmetrical pressure angle, and temporary support position as variable conditions, the model can quickly calculate the internal forces of the support structure corresponding to different curvatures of the center diaphragm. Based on practical engineering, the three-dimensional numerical models of the construction of the large diameter tunnel with CDM were established, and the surrounding deformed characteristics and the mechanical response of the support, which are influenced by the varying curvatures of the center diaphragm, were investigated, respectively.

Theoretical analysis

In the tunnel construction process with the center diaphragm, the internal force of the support structure increases rapidly in the upper bench excavation in leading pilot tunnel, and the upper center diaphragm is the most vulnerable to damage27,32. To explore the influence of the curvature of the center diaphragm on construction at this stage, the center diaphragm curvature \(\kappa\) is set as the core variable, and a mechanical model is developed to quickly calculate the internal force of the support structure, which provides a reference for the design of the curvature of the center diaphragm in CDM. The mechanical model set the support structure as a statically indeterminate three-time double-arch structure. By establishing a coordinate system to determine the geometric relationship of the structure and taking the tunnel covering depth z, unsymmetrical pressure angle \(\eta\), and the intersection positions (x, y) between the center diaphragm and the primary lining as variable parameters, it can be widely applied to different operating conditions of shallow buried tunnels.

Mechanical model

The primary lining-upper center diaphragm double-arch support structure shown in Fig. 2a . The top of the primary lining (left arch) and the upper center diaphragm (right arch) are rigidly connected, and the bottom is fixed by feet-locking bolts. The mechanical model in this paper hypothesized that the support structure materials (including primary lining and upper center diaphragm) were homogeneous with a rectangular cross-section. \(h_1\) represents the height of the upper section; \(l_L\) and \(l_R\) correspond to the projection lengths of the left and right arch structures on the springing line, respectively ; similarly, \(\varphi _L\) and \(\varphi _R\) are the central angles of the left and right structures, respectively; and \(r_L\) and \(r_R\) are their respective radii; \(\alpha _L\) and \(\alpha _R\) represent the angles between the chords of the left and right arcs and the springing line, respectively. The geometric parameters \(\alpha _R\), \(l_R\) and \(\varphi _R\) change with the variation of the center diaphragm curvature \(\kappa\) (e.g., \(\alpha _R'\) and \(l_R'\) in Fig. 2a ), thereby affecting the internal force of the structure. To quantitatively analyze this effect, a coordinate system was established in the following section to relate the aforementioned variables to the curvature \(\kappa\).

Fig. 2
figure 2

Support structure and coordinate system.

Full size image

As shown in Fig. 2b , point O is the intersection of the vertical bisector OB of the line segment connecting the two endpoints A and C of the center diaphragm and the springing line OD. The coordinate system of the support structure takes point O as the origin, the springing line as the x-axis, and constructs the y-axis perpendicular to it.

From the geometric relationships, the equation of circle \(O^{\prime }\), corresponding to the arc of the center diaphragm, can be derived as Eq. (1). And the relation between \(\alpha _R\), \(\varphi _R\), \(l_R\) and \(\kappa\) can be obtained as Eq. (2), where \(k_1\) and \(k_2\) are the slopes of the lines \(O'A\) and \(O'D\), respectively.

$$\begin{aligned} & \begin{aligned}&(x-x_0)^2+(y-y_0)^2=1/\kappa ^2\\&\left\{ \begin{aligned}&x_0=-|OO'|h/\sqrt{b^2+h^2}\\&y_0=-|OO'|b/\sqrt{b^2+h^2}\\&|OO'|=\sqrt{1/\kappa ^2-(b^2+h^2)/4}-\frac{(h_1-h_2)\sqrt{b^2+h^2}}{2b} \end{aligned} \right. \end{aligned} \end{aligned}$$
(1)
$$\begin{aligned} & \begin{aligned}&\left\{ \begin{aligned}&\alpha _R=\tan ^{-1}(\frac{h_1}{\sqrt{ 1/\kappa ^2-{y_0}^2}+x_0 -\frac{{h_1}^2-{h_2}^2-b^2}{2b}})\\&l_R = h_1/\tan \alpha _R\\&\varphi _R = \tan ^{-1}\left( \frac{k_1 - k_2}{1 + k_1 k_2} \right) \\ \end{aligned} \right. \\&\left\{ \begin{aligned}&k_1 = \frac{h_1 - y_0}{(h_1^2 - h_2^2 - b^2)/2b - x_0}\\&k_2 = \frac{-y_0}{\sqrt{1/\kappa ^2 - y_0^2} - x_0} \end{aligned} \right. \end{aligned} \end{aligned}$$
(2)

Under load, the support structure is a three-times statically indeterminate double-arch structure without hinges. In Fig. 3, the covering depths at various points are \(z_i\), \((i = 1, 2, 3, 1', 2')\) and the unsymmetrical pressure angle is \(\eta\). The vertical and horizontal loads acting on the structure are denoted as \(q_i\) and \(p_i\), respectively. Their relationship with \(z_i\) can be expressed by Eq. (3), where \(\gamma _{s}\) represents the unit weight of the surroundings, and \(\lambda\) represents the coefficient of lateral pressure.

Fig. 3
figure 3

Basic system for force method calculation of support structure.

Full size image
$$\begin{aligned} \left\{ \begin{aligned}&q_1=\gamma _{s}z_1\\&q_2=\gamma _{s}z_2=\gamma _{s}(z_1+l_L\tan \eta )\\&q_3=\gamma _{s}z_3=\gamma _{s}(z_1+(l_L\tan \eta +l_R\tan \eta )) \end{aligned} \right. \left\{ \begin{aligned} p_1&=\gamma _{s}z_1\lambda \\ p_2&=\gamma _{s}z_1'\lambda \\&=\gamma _{s}(z_1+h_1)\lambda \end{aligned} \right. \left\{ \begin{aligned} p'_1&=\gamma _{s}z_3\lambda \\ p'_2&=\gamma _{s}z_3'\lambda \\&=\gamma _{s}(z_3+h_1)\lambda \end{aligned} \right. \end{aligned}$$
(3)

In Fig. 3b , the constraints at the intersection of the double-arch structure were released. \(X_1\), \(X_2\), and \(X_3\) represent the horizontal force, vertical force, and bending moment at the top node of the steel frame, respectively. The equation of the force method was established as shown in Eq. (4). Here, \(\delta _{ij}\) denotes the displacement in the \(X_i\) direction at the top node caused by a unit force or moment applied in the \(X_j\) direction. Furthermore, \(\Delta _{ic}\) and \(\Delta _{ip}\) \((i = 1, 2, 3)\) represent displacements at the top node in the Xi direction due to the movement of the support and external forces.

$$\begin{aligned} \left\{ \begin{aligned}&X_1\delta _{11}+X_2\delta _{12}+X_3\delta _{13}+\Delta _{1c}+\Delta _{1p}=0\\&X_1\delta _{21}+X_2\delta _{22}+X_3\delta _{23}+\Delta _{2c}+\Delta _{2p}=0\\&X_1\delta _{31}+X_2\delta _{32}+X_3\delta _{33}+\Delta _{3c}+\Delta _{3p}=0\\ \end{aligned} \right. \end{aligned}$$
(4)

\(\Delta _{ic}\), \((i = 1, 2, 3)\) is obtained from Eq. (5), where \(\beta _L\), \(\beta _R\), \(\mu _L\), \(\mu _R\), \(v_L\) and \(v_L\) are the rotation angles, horizontal displacements and vertical displacements at the supports on the left and right sides, respectively (positive directions are shown in Fig. 3b ).

$$\begin{aligned} \begin{aligned}&\left\{ \begin{aligned} \Delta _{1c}&= - \mu _L - \mu _R - h_1 \beta _L - h_1 \beta _R \\ \Delta _{2c}&= - v_L + v_R + l_L \beta _L - l_R \beta _R \\ \Delta _{3c}&= \beta _L + \beta _R \end{aligned} \right. \\&\left\{ \begin{aligned} \beta _L&= - X_1 (\beta _{1L} + h_1 \beta _{3L}) + X_2 (\beta _{2L} + l_L \beta _{3L}) + X_3 \beta _{3L} + \beta _{pL} \\ \beta _R&= - X_1 (\beta _{1R} + h_1 \beta _{3R}) - X_2 (\beta _{2R} + l_L \beta _{3R}) + X_3 \beta _{3R} + \beta _{pR} \\ \mu _L&= - X_1 (\mu _{1L} + h_1 \mu _{3L}) + X_2 (\mu _{2L} + l_L \mu _{3L}) + X_3 \mu _{3L} + \mu _{pL} \\ \mu _R&= - X_1 (\mu _{1R} + h_1 \mu _{3R}) - X_2 (\mu _{2R} + l_L \mu _{3R}) + X_3 \mu _{3R} + \mu _{pR} \\ v_L&= X_1 (v_{1L} + h_1 v_{3L}) - X_2 (v_{2L} + l_L v_{3L}) - X_3 v_{3L} + v_{pL} \\ v_R&= X_1 (v_{1R} + h_1 v_{3R}) + X_2 (v_{2R} + l_L v_{3R}) - X_3 v_{3R} + v_{pR} \end{aligned} \right. \end{aligned} \end{aligned}$$
(5)

In Eq. (5), \(\beta _{iL}\), \(\beta _{iR}\), \(\mu _{iL}\), \(\mu _{iR}\), \(v_{iL}\) and \(v_{iR }\), \((i=1, 2, 3)\), respectively, represent the rotation angle, horizontal displacement, and vertical displacement of the left and right supports under the action of the unit bending moment, horizontal unit force, and vertical unit force. \(\beta _{iL}\), \(\beta _{iR}\), \(\mu _{iL}\), \(\mu _{iR}\), \(v_{iL}\) and \(v_{iR }\), \((i=1, 2, 3)\)) are related to the setting angle of lock-foot bolt, the deformation coefficient of lock-foot bolt, the resistance coefficient of the surroundings, and other parameters, which can be obtained by using the Winkler foundation beam model33,34,35 and the elastic foundation beam theory36. \(\beta _{pL}\), \(\beta _{pR}\), \(\mu _{pL}\), \(\mu _{pR}\), \(v_{pL}\) and \(v_{pR }\) represent the rotation angle, horizontal displacement, and vertical displacement of the left and right supports under the load of the surroundings, respectively.

\(\Delta _{ip}\), \((i = 1, 2, 3)\) can be solved by Eqs. (6–8), where the elastic modulus, shear modulus, inertial section moment and cross-sectional area of the primary lining are \(E_L\), \(G_L\), \(I_L\), \(A_L\) respectively, and those of the center diagram are \(E_R\), \(G_R\), \(I_R\), \(A_R\) respectively. The definitions of \(x_L\), \(x_R\), \(y_L\), \(y_R\), \(\omega _L\) and \(\omega _R\) are shown in Eq. (9):

$$\begin{aligned} & \begin{aligned} \Delta _{1p}&= \int \frac{M_p \overline{M_1}}{EI}ds+\int \frac{N_p \overline{N_1}}{EA}ds+\int \frac{n V_p \overline{V_1}}{GA}ds\\&=\int _{0}^{\varphi _L}\left( \frac{q_2 {x_L}^2y_L}{2E_LI_L}-\frac{(q_2 - q_1){x_L}^3y_L}{6l_L E_LI_L}+\frac{p_1 {y_L}^3}{2E_LI_L}+\frac{(p_2 - p_1){ y_L}^4}{6h_1 E_LI_L}\right) r_L d\theta \\&\quad +\int _{0}^{\varphi _R}\left( \frac{q_2 {x_R}^2y_R}{2E_RI_R}+\frac{(q_3 - q_2){x_R}^3y_R}{6l_R E_RI_R}+\frac{p'_1 {y_R}^3}{2E_RI_R}+\frac{(p'_2 - p'_1){ y_R}^4}{6h_1 E_RI_R}\right) r_R d\theta \\&\quad +\int _{0}^{\varphi _L}\left( -\frac{q_2 x_L-(q_2 - q_1){x_L}^2/2l_L}{2E_LA_L}\sin 2\omega _L+ \frac{p_1 y_L+(p_2 - p_1)y_L^2/2h_1}{E_LA_L}\cos ^2\omega _L\right) r_L d\theta \\&\quad +\int _{0}^{\varphi _R}\left( -\frac{q_2 x_R+(q_3 - q_2){x_R}^2/2l_R}{2E_RA_R}\sin 2\omega _R+ \frac{p'_1 y_R+(p'_2 - p'_1)y_R^2/2h_1}{E_RA_R}\cos ^2\omega _R\right) r_R d\theta \\&\quad +\int _{0}^{\varphi _L}\left( \frac{n(q_2 x_L-(q_2 - q_1){x_L}^2/2l_L)}{2G_LA_L}\sin 2\omega _L+ \frac{n(p_1 y_L+(p_2 - p_1)y_L^2/2h_1)}{G_LA_L}\sin ^2\omega _L\right) r_L d\theta \\&\quad +\int _{0}^{\varphi _R}\left( \frac{n(q_2 x_R+(q_3 - q_2){x_R}^2/2l_R)}{2G_RA_R}\sin 2\omega _R+ \frac{n(p'_1 y_R+(p'_2 - p'_1)y_R^2/2h_1)}{G_RA_R}\sin ^2\omega _R\right) r_Rd\theta \\ \end{aligned} \end{aligned}$$
(6)
$$\begin{aligned} & \begin{aligned} \Delta _{2p}&= \int \frac{M_p \overline{M_2}}{EI}ds+\int \frac{N_p \overline{N_2}}{EA}ds+\int \frac{n V_p \overline{V_2}}{GA}ds\\&=\int _{0}^{\varphi _L}\left( -\frac{q_{2}x_{L}^{3}}{2E_LI_L}+\frac{(q_{2}-q_{1})x_{L}^{4}}{6l_{L}E_LI_L}-\frac{p_{1}x_{L}y_{L}^2}{2E_LI_L}-\frac{(p_{2}-p_{1})x_{L}y_{L}^{3}}{6h_{1}E_LI_L}\right) r_{L}d\theta \\&\quad +\int _{0}^{\varphi _R}\left( \frac{q_{2}x_{R}^{3}}{2E_RI_R}+\frac{(q_{3}-q_{2})x_{R}^{4}}{6l_{R}E_RI_R}+\frac{p'_{1}x_{R}y_{R}^2}{2E_RI_R}+\frac{(p'_{2}-p'_{1})x_{R}y_{R}^3}{6h_{1}E_RI_R}\right) r_{R}d\theta \\&\quad +\int _{0}^{\varphi _L}\left( -\frac{q_{2}x_{L}-(q_{2}-q_{1})x_{L}^{2}/2l_{L}}{E_LA_L}\sin ^{2}\omega _{L}+\frac{p_{1}y_{L}+(p_{2}-p_{1})y_{L}^{2}/2h_{1}}{2E_LA_L}\sin 2\omega _{L}\right) r_{L}d\theta \\&\quad +\int _{0}^{\varphi _R}\left( \frac{q_{2}x_{R}+(q_{3}-q_{2})x_{R}^{2}/2l_{R}}{E_RA_R}\sin ^{2}\omega _{R}-\frac{p'_{1}y_{R}+(p'_{2}-p'_{1})y_{R}^{2}/2h_{1}}{2E_RA_R}\sin 2\omega _{R}\right) r_{R}d\theta \\&\quad +\int _{0}^{\varphi _L}\left( -\frac{n(q_{2}x_{L}-(q_{2}-q_{1})x_{L}^{2}/2l_{L})}{G_LA_L}\cos ^{2}\omega _{L}-\frac{n(p_{1}y_{L}+(p_{2}-p_{1})y_{L}^{2}/2h_{1})}{2G_LA_L}\sin 2\omega _{L}\right) r_{L}d\theta \\&\quad +\int _{0}^{\varphi _R}\left( \frac{n(q_{2}x_{R}+(q_{3}-q_{2})x_{R}^{2}/2l_{R})}{G_RA_R}\cos ^{2}\omega _{R}+\frac{n(p'_{1}y_{R}+(p'_{2}-p'_{1})y_{R}^{2}/2h_{1})}{2G_RA_R}\sin 2\omega _{R}\right) r_{R}d\theta \end{aligned} \end{aligned}$$
(7)
$$\begin{aligned} & \begin{aligned} \Delta _{3p}&=\int \frac{M_{p}\overline{M}_{3}}{EI}ds\\&=\int _{0}^{\varphi _L}\left( -\frac{q_{2}x_{L}^{2}}{2E_LI_L}+\frac{(q_{2}-q_{1})x_{L}^{3}}{6l_{L}E_LI_L}-\frac{p_{1}y_{L}^2}{2E_LI_L}-\frac{(p_{2}-p_{1})y_{L}^{3}}{6h_{1}E_LI_L}\right) r_{L}d\theta \\&\quad +\int _{0}^{\varphi _R}\left( -\frac{q_{2}x_{R}^{2}}{2E_RI_R}-\frac{(q_{3}-q_{2})x_{R}^{3}}{6l_{R}E_RI_R}-\frac{p'_{1}y_{R}^{2}}{2E_RI_R}-\frac{(p'_{2}-p'_{1})y_{R}^{3}}{6h_{1}E_RI_R}\right) r_{R}d\theta \end{aligned} \end{aligned}$$
(8)
$$\begin{aligned} & \left\{ \begin{aligned} x_{R}&=l_{R}-2r_{R}\sin \left( \left( \varphi _{R}-\theta \right) /2\right) \cos \left( \alpha _{R}+\theta /2\right) \\ y_{R}&=h_{1}-2r_{R}\sin \left( \left( \varphi _{R}-\theta \right) /2\right) \sin \left( \alpha _{R}+\theta /2\right) \\ x_{L}&=l_{L}-2r_{L}\sin \left( \left( \varphi _{L}-\theta \right) /2\right) \cos \left( \alpha _{L}+\theta /2\right) \\ y_{L}&=h_{1}-2r_{L}\sin \left( \left( \varphi _{L}-\theta \right) /2\right) \sin \left( \alpha _{L}+\theta /2\right) \\ \omega _{L}&=\theta +\alpha _{L}-\varphi _{L}/2\\ \omega _{R}&=\theta +\alpha _{R}-\varphi _{R}/2 \end{aligned} \right. \end{aligned}$$
(9)

The compliance coefficient \(\delta _{ij}\), \((i, j = 1, 2, 3)\) can be determined by Eq. (10):

$$\begin{aligned} \left\{ \begin{aligned} \delta _{11}&= \int _{0}^{\varphi _L} \left( \frac{y_L^2}{E_LI_L} + \frac{\cos ^2 \omega _L}{E_LA_L} + \frac{\sin ^2 \omega _L}{G_LA_L} \right) r_L d\theta + \int _{0}^{\varphi _R} \left( \frac{y_R^2}{E_RI_R} + \frac{\cos ^2 \omega _R}{E_RA_R} + \frac{\sin ^2 \omega _R}{G_RA_R} \right) r_R d\theta \\ \delta _{22}&= \int _{0}^{\varphi _L} \left( \frac{x_L^2}{E_LI_L} + \frac{\sin ^2 \omega _L}{E_LA_L} + \frac{\cos ^2 \omega _L}{G_LA_L} \right) r_L d\theta + \int _{0}^{\varphi _R} \left( \frac{x_R^2}{E_RI_R} + \frac{\sin ^2 \omega _R}{E_RA_R} + \frac{\cos ^2 \omega _R}{G_RA_R} \right) r_R d\theta \\ \delta _{33}&= \frac{\varphi _L r_L}{E_LI_L} +\frac{\varphi _R r_R}{E_RI_R}\\ \delta _{12}&= \delta _{21} = \int _{0}^{\varphi _L} \left( -\frac{x_L y_L}{E_LI_L} + \frac{\sin 2\omega _L}{2E_LA_L} - \frac{\sin 2\omega _L}{2G_LA_L} \right) r_L d\theta + \int _{0}^{\varphi _R} \left( \frac{x_R y_R}{E_RI_R} - \frac{\sin 2\omega _R}{2E_RA_R} +\frac{\sin 2\omega _R}{2G_RA_R} \right) r_R d\theta \\ \delta _{13}&= \delta _{31} = - \int _{0}^{\varphi _L} \frac{y_L r_L}{E_LI_L} d\theta - \int _{0}^{\varphi _R} \frac{y_R r_R}{E_RI_R}d\theta \\ \delta _{23}&= \delta _{32} = \int _{0}^{\varphi _L} \frac{x_L r_L}{E_LI_L} d\theta - \int _{0}^{\varphi _R} \frac{x_R r_R}{E_RI_R} d\theta \end{aligned} \right. \end{aligned}$$
(10)

At this juncture, the constraint force \(X_i\), \((i = 1, 2, 3)\) at the intersection of the center diaphragm and the primary lining can be obtained by solving the Eqs. (6–10) . Stipulating the bending moment is positive when the inner side of the structure is in tension, and the shear force is positive when it points to the center of the circle. The bending moment \(M_L(\theta _0)\), \(M_R(\theta _0)\) and the shear force \(V_L(\theta _0)\), \(V_R(\theta _0)\) at any cross section on the primary lining and the center diaphragm can be calculated with Eq. (11).

$$\begin{aligned} \left\{ \begin{aligned} M_L(\theta _0)&= -\left( \frac{q_2 x_L^2}{2} - \frac{(q_2 - q_1)x_L^3}{6l_L} + \frac{p_1 y_L^2}{2} + \frac{(p_2 - p_1)y_L^3}{6h_1} \right) - X_1 y_L + X_2 x_L + X_3 \\ M_R(\theta _0)&= -\left( \frac{q_2 x_R^2}{2} + \frac{(q_3 - q_2)x_R^3}{6l_R} + \frac{p'_1 y_R^2}{2} + \frac{(p'_2 - p'_1)y_R^3}{6h_1} \right) - X_1 y_R - X_2 x_R + X_3 \\ V_L(\theta _0)&= -\left( p_1 y_L + \frac{ (p_2 - p_1)y_L^2}{2h_1}\right) \sin \omega _L +\left( q_2 x_L - \frac{(q_2 - q_1)x_L^2}{2l_L} \right) \cos \omega _L - X_1 \sin \omega _L + X_2 \cos \omega _L\\ V_R(\theta _0)&= -\left( p'_1 y_R + \frac{ (p'_2 - p'_1)y_R^2}{2h_1}\right) \sin \omega _R -\left( q_2 x_R + \frac{(q_3 - q_2)x_R^2}{2l_R} \right) \cos \omega _R - X_1 \sin \omega _R - X_2 \cos \omega _R \end{aligned} \right. \end{aligned}$$
(11)

Application of mechanical model

For practical engineering, the proposed model can be applied following the flowchart in Fig. 4. Notably, for preliminary or approximate calculations, the displacements and unit displacements from axial and shear forces in Eqs. (6, 7 and 10) can be ignored.

Fig. 4
figure 4

The application process of the mechanical model.

Full size image

The support structure (two-dimensional model) was built in Midas GTS NX using the load structure method. The internal force outcomes of the finite element and mechanical models matched well, validating the accuracy of the theoretical analysis.

Numerical analysis

The influence of the change in the curvature of the center diaphragm on the whole tunnel construction process can be more truly and comprehensively reflected in the numerical model of the three-dimensional tunnel. Based on practical engineering of the Panlong Tunnel in China, a set of numerical simulations are conducted with the curvature of the center diaphragm as a variable for an in-depth exploration of the tunneling response.

Project overview

The Panlong Tunnel serves as a crucial node of the Yan’an East Bypass Expressway, which is an essential link in promoting rapid economic development of the northern Shaanxi region, China. There are four lanes in both directions on the traffic lanes of the tunnel. The terrain at the entrance and exit sections of it is relatively steep, exposing geological formations including Triassic sandstone and mudstone (\(T_{3y}\)), Middle Pleistocene Lishi Loess (\(Q_{2}^{eol}\)), Upper Pleistocene Malan Loess (\(Q_{3}^{eol}\)), and quaternary landslide deposits and slope accumulations (\(Q_{4}^{del}\), \(Q_{4}^{dl}\))(Fig. 5). The study was based on a selected segment of the left line of the Panlong Tunnel, specifically from ZK28+626 to ZK28+656. This segment has an average overburden depth of 0-30 m, which is categorize as a shallow-buried tunnel with significant unsymmetrical pressure. The surroundings consists mainly of Lishi Loess (\(Q_{2}^{eol}\)), which is uniform in composition and compact in structure. However, without adequate support, the arch section is prone to collapse after excavation, indicating poor rock mass stability37,38. The comprehensive evaluation classifies the surroundings as Grade V. A composite support system was adopted for the tunnel, with the support of fore-poling pipe roof and construction employing the CDM. The tunnel support system is shown in Fig. 6.

Fig. 5
figure 5

Project overview (China Map: Map Review No. GS(2019)1659 (Revised); Yan’an East Bypass Expressway Map: Source: Windows Maps (V. 11.2411.7.0, URL: https://apps.microsoft.com/detail/9wzdncrdtbvb)).

Full size image
Fig. 6
figure 6

Tunnel support system.

Full size image

Numerical model

The three-dimensional tunnel model was constructed using the Midas GTS NX finite element simulation software to dynamically simulate the construction processes for the Panlong Tunnel. To eliminate the influence of boundaries, the dimensions of the model are 93 m × 3 0 m × 83.7 m. Taking into account the actual conditions of the project, the distance from the tunnel crown to the ground surface was set to 24.0 m. The upper surface of the model was designed with a slope from left-low to right-high, featuring a gradient angle of \(21.6^\circ\). The upper surface of the model was set as unrestrained, while the remaining external surfaces were constrained in the normal displacement direction. The permanent support system of the model consists of bolts, grouting reinforcement zones, pipe roofs, primary lining, and secondary lining. Taking the four curvatures of the temporary support (center diaphragm) as the research variable, the situations corresponding to the four working conditions were analyzed. The schematic diagram of the model is shown in Fig. 7.

Fig. 7
figure 7

Tunnel numerical model.

Full size image

Constitutive model and material parameters

The Mohr-Coulomb constitutive model was applied for the surroundings and advanced support grouting reinforcement zone of the tunnel. And the effect of fore-poling support grouting reinforcement was simulated by adjusting the parameters of the reinforced surroundings. The fore-poling pipe roof, bolts, primary lining, temporary support, and secondary lining were all modeled using the elastic constitutive model. To simplify the model, the equivalent elastic moduli E of the primary lining, the temporary support (center diaphragm) and the secondary lining were calculated according to Eq. (12), where \(E_C\), \(A_C\), \(E_S\) and \(A_S\) respectively represent the elastic modulus of concrete, its cross-sectional area, the elastic modulus of the steel frame and that of the steel frame. The elastic modulus \(E'\) of the pipe roof and the equivalent unit weight \(\gamma '\) of grouting inside the pipes were calculated with reference to Eq. (13), where \(E_1\), \(I_1\), \(\gamma _1\), \(A_1\) are the elastic modulus, moment of inertia, unit weight, and cross-sectional area of the steel pipe, respectively, and \(E_2\), \(I_2\), \(\gamma _2\), \(A_2\) represent the corresponding parameters for the grout (Table 1).

$$\begin{aligned} E = E_C + (A_S + E_S)/A_C \end{aligned}$$
(12)
$$\begin{aligned} \left\{ \begin{aligned}&E'=(E_1I_1+E_2I_2)/(I_1+I_2)\\&\gamma '=(A_1\gamma _1+A_2\gamma _2)/(A_1+A_2) \end{aligned} \right. \end{aligned}$$
(13)
Table 1 Mechanical parameters of materials39,40.
Full size table

Construction simulation

The process of simulation construction and mounting of the main support structures are shown in Fig. 8. CDM was carried out, with a cyclic excavation advance of 1.5 m. The excavation spacing for the upper and lower benches is 6 m, and the excavation spacing of the leading and trailing pilot tunnel is 21 m.

Fig. 8
figure 8

Process of simulation construction.

Full size image

Results

The middle cross-section along the excavation direction was selected as the study section to minimize the effects of the boundary conditions on the simulation results. The research focused on the construction responses caused by the change in the curvature of the center diaphragm. The intersection of the center diaphragm and the primary lining was set as fixed nodes, and the curvatures of the center diaphragms, \(\kappa\), were set to \(0, 0.05, 0.10\) and \(0.15\). Among them, the case with \(\kappa = 0\) corresponds to a straight center diaphragm, while the cases with \(\kappa = 0.05\), \(\kappa = 0.1\), and \(\kappa = 0.15\) correspond to curved center diaphragms. The corresponding radii \(r_R\) were \(\infty\) m, \(20.00\) m, \(10.00\) m and \(6.67\) m respectively, where \(\kappa =\frac{1}{r_R}\). This resulted in four configurations of the center diaphragm, as shown in Fig. 9.

Fig. 9
figure 9

Center diaphragms of different curvatures within tunnel numerical models.

Full size image

Displacement analysis

Surface settlement

Figure 10 illustrates the surface settlement curves for models with different curvatures of the center diaphragms at the end of the construction stage. The results showed that the widths of the settlement troughs were nearly consistent among all models, measuring approximately 80 m. A slightly greater settlement was observed on the right side compared to the left, which can be attributed to unsymmetrical pressure conditions. When the curvature of the center diaphragm increased from 0 to 0.15, the maximum surface settlement increased slightly, from 22.94 mm to 24.27 mm, with an increase rate of 0.06%. However, the differences between the models remained minimal, indicating that the curvature of the center diaphragm had a negligible impact on the surface settlement results.

Fig. 10
figure 10

Surface settlement curves.

Full size image

Crown settlement

Taking the model with a curvature of the center diaphragm of 0.1 as an example, Fig. 11 illustrates the variation of vertical displacement at the crown of the leading and trailing pilot tunnels41. Prior to and after the upper bench excavation of the leading (step 12) and trailing pilot tunnels (step 26), a common pattern was observed for displacement measurement points. Both points 1 and 2 typically experienced a settlement process that transitioned from slow to fast and back to slow. However, during excavation of the upper bench of the leading pilot tunnel, the range of settlement changes at point 1 surpassed that at point 2. In contrast, during the excavation of the upper bench of the trailing pilot tunnel, the situation reversed. After installing the primary lining for the leading and trailing pilot tunnels, the crown settlement values \(\Delta\) \(h_1\) and \(\Delta\) \(h_2\) were 17.61 mm and 19.25 mm, respectively. The discrepancy in the settlement values could be attributed to the installation of the primary lining of the upper bench in the leading pilot tunnel independently, without sharing the load with the primary lining of the trailing pilot tunnel. Furthermore, as excavation progressed in the trailing pilot tunnel, the crown of the leading tunnel continued to experience additional settlement. The cumulative crown settlements \(\Delta\) \(H_1\) and \(\Delta\) \(H_2\) were 23.65 mm and 42.69 mm, respectively. The larger settlement in the trailing pilot tunnel was attributed to the unsymmetrical pressure and the excavation sequence.

Fig. 11
figure 11

Crown settlement changes of surroundings at displacement measurement points 1 and 2.

Full size image

Histograms of post-support crown settlement \(\Delta\) \(h\), at the crowns of the leading and trailing pilot tunnels in models withdifferent curvatures of the center diaphragm \(\kappa\) (Fig. 12a), and histograms of cumulative crown settlement \(\Delta\) \(H\) (Fig. 12b) are presented in Fig. 12. The crown settlement of the leading pilot tunnel was significantly affected by the curvature of the center diaphragm, with a clear pattern. Both post-support and cumulative crown settlement of measure point 1 increased as the curvature increased. As \(\kappa\) increased from 0 to 0.15, \(\Delta h_1\) increased by \(59.5\%\) and \(\Delta H_1\) increased by \(43.1\%\).While the settlement values \(\Delta h_2\) and \(\Delta H_2\) at the corresponding measuring point 2 of the trailing pilot tunnel fluctuated within the ranges of \(10\%\) and \(5\%\), respectively. This indicated that the curvature of the center diaphragm significantly affected the crown settlement of the leading pilot tunnel, but had little impact on that of the trailing pilot tunnel.

Fig. 12
figure 12

Crown settlement at measurement points 1 and 2.

Full size image

Stress analysis

Primary lining stress

To analyze the stress distribution of the primary lining, ten measurement points were selected. Figure 13 shows the maximum and minimum principal stresses at these measurement points before dismantling the center diaphragm and in the final construction step, respectively. The analysis focused on the absolute values of stress, with positive values indicating tensile stress and negative values representing compressive stress.

Fig. 13
figure 13

Distribution curves of principal stress for primary lining.

Full size image

The relationship between the distribution of maximum principal stress in the primary lining and the curvature of the center diaphragm is shown in Fig. 13a and b. Before dismantling the center diaphragm (Fig. 13a), the crown and shoulders were the most affected by curvature. As the curvature increased from 0 to 0.15, the maximum principal stresses at the corresponding measuring points 1, 2, and 3 decreased, with the decreases being 2.79 MPa, 2.17 MPa, and 1.98 MPa, respectively. Tensile stress at the crown decreased and shoulder stresses shifted from tensile to compressive. The curvature change had less impact on other parts. In the final step of construction (Fig. 13b ), the shoulders and invert were significantly affected by curvature. As the curvature of the center diaphragm increased, the maximum principal stresses at the left shoulder and invert increased, while the stress at the right shoulder decreased. All of these parts were under compression. When the curvature increased from 0 to 0.15, the compressive stresses at the left shoulder (point 2) and the invert (point 10) decreased by 2.81 MPa and 1.93 MPa, respectively, while the stress at the right shoulder (point 3) increased by 2.51 MPa. The curvature change had less impact on other parts. When comparing the stress distributions before and after center diaphragm dismantling, the maximum principal stress at the left shoulder and the invert decreased significantly after dismantling, which meant that the compressive stress increased. This indicated that the compressive stress near the center diaphragm of the primary lining increased after removal.

The connection between the minimum principal stress distribution of the primary lining and the curvature is presented in Fig. 13c and d . The primary lining was compressed in the direction of the minimum principal stress. Before dismantling the center diaphragm (Fig. 13c ), for the curved center diaphragms (\(\kappa = 0.05, 0.10\) and 0.15), the compressive stress of all parts of the primary lining increased with curvature. Specifically, for the crown, shoulders, waists, left side wall and invert, when the curvature increased from 0.05 to 0.15, the increase in stress at the corresponding measuring points 1-6 and 10 was 29.0%, 27.8%, 6.3%, 11.8%, 10.8%, 8.9%, and 49.3%, respectively. For other parts, the increase in stress due to curvature was within 5%.Although the straight center diaphragm had the minimum curvature (\(\kappa = 0\)), the stress of all parts was not minimum, reflecting differences between the straight and curved center diaphragms. In the final construction step (Fig. 13d ), for the curved center diaphragms, the compressive stress at the crown, right shoulder, waists, side walls, and feet increased with curvature. When the curvature increased from 0.05 to 0.15, the stress increased at the corresponding measuring points 1, 3-9, which were 13.0%, 18.3%, 37.0%, 15.0%, 47.2%, 7.3%, 294.8%, and 17.0%, respectively. However, the compressive stress at the left shoulder and invert decreased with increasing curvature. At these two positions, the curved and straight center diaphragms followed a unified pattern. When the curvature increased from 0 to 0.15, the stress at the left shoulder decreased by 32.1% and at the invert by 12.9%. The stress distribution for the straight center diaphragm was special, not forming a common pattern with the curved center diaphragms at the crown, right shoulder, left waist, both side walls, and left foot, reflecting differences caused by different curvatures. When comparing the stress distributions of the primary lining before and after the dismantling of the center diaphragm, compressive stress at the crown, left shoulder, and invert increased significantly post-dismantling, confirming the increased compressive stress near the original position of the center diaphragm in the primary lining after its removal.

Center diaphragm stress

Stress concentration tends to occur at the upper end, middle part, and lower end of the center diaphragm, making these areas the most vulnerable29. The stress measurement points were set at these locations, and the maximum values of the maximum and minimum principal stress during the simulated construction process were recorded in the form of histograms (Fig. 14). Figure 14a shows that the maximum principal stress of each part of the curved center diaphragms was significantly higher than those of the straight center diaphragm, implying a higher susceptibility to tensile failure. For the curved center diaphragms, the maximum value of the maximum principal stress was at the lower end, indicating a risk of tensile failure39. The stress at the lower end decreased as the curvature increased. When the curvature increased from 0.05 to 0.15, the stress at the corresponding measurement point 3 decreased by 31.2%. In contrast, for the straight center diaphragm, the maximum value was at the upper end. Figure 14b shows that the maximum value of the minimum principal stress for the center diaphragms with various curvatures was at the middle part, indicating that this part bore the maximum compressive stress, illustrating that this part bore the maximum compressive stress. As the curvature increased from 0 to 0.15, the maximum value at the corresponding measuring point 2 increased by 80.9%. Only considering the curved center diaphragms, as the curvature increased from 0.05 to 0.15, the maximum stress increased by 40.9%.

Fig. 14
figure 14

Maximum values of principal stress at each stress measurement point on the center diaphragms.

Full size image

To explore how the curvature of the center diaphragm affects the stress changes at locations with maximum values of principal stress during construction, Fig. 15 presents the variation curves of the maximum principal stress at measurement point 3 (Fig. 15a ) at the lower end and the minimum principal stress at measurement point 2 (Fig. 15b ) at the middle part of the center diaphragm as the simulation progresses. Note that the maximum value of the maximum principal stress in the straight center diaphragm appeared at measuring point 1 located at its upper end. Given the small stress value and the absence of failure risk, it was not discussed further.

Fig. 15
figure 15

Maximum values of principal stress at each stress measurement point on the center diaphragms.

Full size image

Figure 15a shows that as construction progressed, the stress of the curved center diaphragms initially increased, then decreased, reaching its peak before the primary lining was fully installed. The possible reason was that, once the primary lining formed a ring, the stress was transferred and dispersed from the center diaphragm. Throughout construction, the curved center diaphragms bore a lower stress with a greater curvature. However, the maximum principal stress at the lower end of the straight center diaphragm varied little, ranging from −1.5 to 0.5 MPa. These differences further reflect the influence of curvature on construction stress.

Figure 15b demonstrates that the minimum principal stress in the middle part of the curved center diaphragms with different curvatures first increased and then decreased. After the center diaphragm and the primary lining of the leading pilot tunnel formed a ring, stress stabilized within the support system and its increase periodically slowed. The maximum stress was reached before the construction of the primary lining of the upper bench of the trailing pilot tunnel. The possible reason was that the primary lining of the upper bench in the trailing section, combined with the primary lining of the leading pilot tunnel, formed a stable arched support system, which dispersed the stress on the center diaphragm. Throughout construction, the stress of the curved center diaphragms increased with increasing curvature. The stress of the straight center diaphragm rose continuously, with the primary lining construction not significantly dispersing the stress, highlighting the differences between the curved and straight center diaphragms and demonstrating the impact of curvature on the construction stress.

Conclusions

This study developed a mechanical model for the support structure and three-dimensional numerical tunneling models to investigate the effects of the curvature of the center diaphragm on the shallow tunnel with large-span on slope topography. The conclusions of this work can be summarized as follows:

  1. (1)

    A method applicable to CDM was proposed for quickly determining the internal forces of the support structure. This method established a double-arch model of primary lining-upper center diaphragm double-arch model and its coordinate system, linking the curvature of the center diaphragm with other geometric parameters of the structure, while taking into account factors of tunnel covering depth, unsymmetrical pressure angle, and temporary support ___location.

  2. (2)

    The influence of the curvature of the center diaphragm on surface settlement was negligible. When the curvature of the center diaphragm increased, the maximum surface settlement increased, but the increase rate was within 0.1%.

  3. (3)

    The curvature had a significant impact on the crown settlement of the leading pilot tunnel, while having a negligible impact on that of the trailing pilot tunnel. As the curvature increased from 0 to 0.15, the post-support and cumulative crown settlements of the leading pilot tunnel increased by 59.5% and 43.1%, respectively, while those of the trailing pilot tunnel changed within a range of 10%.

  4. (4)

    The dismantling of the center diaphragm increased the compressive stress at the positions adjacent to the center diaphragm (left shoulder and invert) of the primary lining. When the curvature of the center diaphragm increased from 0 to 0.15, the compressive stress at these two positions decreased by 32.1% and 12.9%, respectively.

  5. (5)

    During the construction process, for the curved center diaphragms, the maximum principal stress decreased as the curvature increased, while the minimum principal stress (in absolute value) increased with increasing curvature. When the curvature increased from 0.05 to 0.15, the maximum values of the two principal stresses decreased by 31.2% and increased by 40.9%, respectively. The maximum tensile stress and compressive stress of the straight center diaphragm were both smaller than those of the curved center diaphragms, and the straight center diaphragm exhibited a higher level of safety performance.