考虑接触非线性的盾构隧道纵向分析模型及其应用

摘要

盾构隧道环缝接触状态与传递的内力有关，但传统盾构隧道结构分析模型不能考虑轴力、剪力、弯矩共同作用下的环缝接触非线性。提出了一种考虑3种内力共同作用下接触非线性的盾构隧道纵向分析模型，并据此对汕头湾盾构隧道穿越蠕滑断层的力学行为进行分析。研究发现：①当考虑纵向接头同时传递轴力、剪力、弯矩时，其传力作用可等效为考虑轴向和弯曲刚度折减的梁，其折减系数与位移偏心距有关；②汕头湾盾构隧道在蠕滑作用下的力学分析表明，纵向接头接触非线性对结构内力影响较大，忽视非线性的影响可能带来1.4~3.3倍的误差。研究成果可用于蠕滑断层位错等复杂荷载作用下的盾构隧道纵向力学行为分析。

关键词

盾构隧道; 纵向非线性分析; 断层蠕滑

断层经常造成岩层之间的不整合面，可能对隧道、地基等的力学性质造成影响。活动断层发生日积月累的蠕滑则会对隧道、地面建筑等造成不可逆的影响，因此各国建筑规范对在活动断层区域进行工程建造有严格的要求。然而，对于隧道等类型的线状工程，有时不可避免需要穿越蠕滑断层，因此，在长期蠕滑作用下，结构力学性能预测尤为重要。

因为不连续的结构构造，盾构隧道力学分析为各类穿越蠕滑断层隧道分析中最为复杂的分析之一，其面临的主要困难是盾构隧道纵向接头的非线性问题。目前盾构隧道纵向接头分析模型主要有2类：管片-环缝模型和接缝界面模型。其中，管片-环缝模型中接头力学分析考虑管片和环缝参与受力，基于盾构隧道分别在轴拉力、轴压力、纯弯矩作用下进行分析，环缝分别处于全张开、全闭合、半张开状态，常采用考虑刚度折减的连续梁等效。文献［1］最早研究了该类模型，通过引入平截面假定，推导了盾构隧道纵向接头在纯轴力和纯弯矩作用下的轴向刚度和弯曲刚度解析解。接缝界面模型中接头力学分析不考虑管片参与受力，以接缝在外力作用下的力学行为为分析对象，常采用轴向、剪切、弯曲弹簧模拟。文献［2-3］对隧道纵断面方向上环间接头转动弹簧系数进行了理论推导。2类方法中管片-环缝模型应用最广泛，徐凌^［4］、黄正荣^［5］、张文杰等^［6］、叶飞等^［7-8］、李翔宇等^［9］通过理论分析、模型试验等手段研究了横向变形、环缝作用范围和螺栓预紧力等因素对盾构隧道纵向等效抗弯刚度的影响，改进了该模型；黄亮等^［10］采用与志波模型相同的思路，建立了类矩形盾构隧道纵向弯曲刚度的计算方法；文献［11-14］则重点关注了纵向弯曲和剪切效应；卢慈荣^［15］、耿萍等^［16］、张景等^［17］针对弯矩和轴力共同作用下盾构隧道纵向抗弯刚度非线性变化进行了研究，建立了轴力、弯矩共同作用下的纵向接头力学分析方法。

综合来看，盾构隧道纵向结构为传递轴力、弯矩、剪力作用的构造，具有轴向、转动、切向3个自由度。现有纵向接头分析模型进行力学分析时，一般将3个自由度独立考虑，然而，在一些特殊工况下盾构隧道同时传递着轴力、弯矩和剪力，如盾构隧道在蠕滑位错荷载作用下的力学响应、在斜入射地震动作用下的地震响应，传统纵向接头分析模型不能考虑轴力、弯矩、剪力共同作用下的接头力学特性，使得传统分析方法在特殊工况下不适用。

笔者曾对轴力和弯矩共同作用下的盾构隧道纵向抗弯刚度进行了解析^［16］，然而该模型亦存在如下缺点：模型由多种接触状态分别表达，需要分状态计算，且受压区计算高度复杂；仅考虑轴力和弯矩的共同作用，未考虑剪力作用。本文通过引入基于位移的接头力学分析方法，建立轴力、剪力、弯矩共同作用下考虑接触非线性的盾构隧道纵向分析模型；然后以汕头湾盾构隧道为背景，分析蠕滑作用下盾构隧道结构内力响应。与轴力和弯矩共同作用模型相比，2个模型采用了相同的基本假定，但轴力和弯矩共同作用模型以轴力和弯矩为已知量。本文模型以位移为已知量，主要有如下改进：在相同假定的条件下，能考虑轴向、切向、转动3个自由度的相互影响；采用基于位移的推导，将环间接缝多种接触状态使用统一力学方程表达。

1 盾构隧道纵向分析模型推导

1.1　盾构隧道纵向力学分析的基本问题

当考虑轴力、剪力、弯矩共同作用时，取分析区域为相邻管片环中心截面间的区域，分别为截面1和截面2，分析区域长度为$l_{\mathrm{s}}$，如图 1 所示，其中作用于截面1和截面2上的位移和力的大小为$(u_i,v_i,{\theta }_i)$和$(N_i,V_i,M_i)$（$i=\mathrm{1,2}$）。

图1 盾构隧道纵向力学分析

Fig. 1 Basic problems in longitudinal mechanical analysis of shield tunnel

分析区域控制方程可用面2与面1的相对位移表示，如式（1）所示。

$$\left[\begin{array}{c}N\\ V\\ M\end{array}\right]=\left[\begin{array}{ccc}{k}_{\mathrm{11}}& k_{\mathrm{12}}& k_{\mathrm{13}}\\ k_{\mathrm{21}}& k_{\mathrm{22}}& k_{\mathrm{23}}\\ k_{\mathrm{31}}& k_{\mathrm{32}}& k_{\mathrm{33}}\end{array}\right]\left[\begin{array}{c}u\\ v\\ \theta \end{array}\right]$$

（1）

式中：M、V、N为截面2上的作用力；u、v、θ为截面2与截面1的相对位移，见式（2）~（4）。

$$u=u_{\mathrm{2}}-u_{\mathrm{1}}$$

（2）

$$v=v_{\mathrm{2}}-v_{\mathrm{1}}+l_{\mathrm{s}}{\theta }_{\mathrm{1}}$$

（3）

$$\theta ={\theta }_{\mathrm{2}}-{\theta }_{\mathrm{1}}$$

（4）

因此，盾构隧道纵向分析的基本问题可表述为：在考虑隧道结构存在环间接缝的情况下，确定如式（1）所示的3个自由度上力与位移的关系，并计算隧道结构的环间接缝张开量、螺栓拉力、管片应力的问题。

引入如下基本假设：①假定管片环中心截面在同一平面内变形；②相邻管片环中心截面范围内的平均变形服从平截面假定；③环间接缝受压侧仅管片受压，螺栓不参与受力，受拉侧管片和螺栓共同受拉；④连接件（如螺栓）在环向的离散分布转化为等效的连续均匀分布；⑤不考虑环间接缝的滑移。

1.2　基于位移的力学分析方法

通过尝试基于力的方法和基于位移的方法，发现：当采用基于力的方法时需要建立力平衡方程和变形协调方程，联立求解常得到以受压区高度为未知量的隐式超越方程，使得内力和变形的求解十分复杂；当采用基于位移的推导方法时，接头分析区域变形分布容易求得，引入考虑了接头不连续条件的物理方程后，直接积分可得内力，使得问题的分析得到简化。因此，本文采用基于位移的接头力学分析方法，如图 2 所示。该方法以位移作为已知量，利用平截面假定得到分析区域变形函数，然后建立盾构隧道纵向力学方程。

图2 基于位移的力学方程推导

Fig. 2 Derivation of mechanical equation based on displacement

1.3　盾构隧道纵向力学方程推导

1.3.1　分析区域变形函数

由平截面假定可知分析区域的挠曲变形符合梁的挠曲变形函数，轴向变形符合轴心受力构件变形函数。若将分析区域截面1固定，在截面2上施加位移$\mathit{U}={\left(u,v,\theta \right)}^{\mathrm{T}}$，则分析区域的挠曲函数和轴向变形函数分别如式（5）、式（6）^［18］所示。

$$\lambda \left(\xi \right)=\left(\mathrm{3}{\xi }^{\mathrm{2}}-\mathrm{2}{\xi }^{\mathrm{3}}\right)v+\left({\xi }^{\mathrm{3}}-{\xi }^{\mathrm{2}}\right)l_{\mathrm{s}}\theta$$

（5）

$$u\left(\xi \right)=\xi u$$

（6）

式中：$\xi =x/l_{\mathrm{s}}$，x为以截面1中心点为原点、指向截面2中心点坐标轴上的坐标，如图 3 所示；$l_{\mathrm{s}}$为管片环宽度。

图3 纵向受力分析示意

Fig. 3 Diagram of longitudinal mechanical analysis

由式（5）可得，各截面相对于截面1的转角为

$$\theta \left(\xi \right)=\frac{\mathrm{d}\lambda \left(\xi \right)}{\mathrm{d}\xi }\frac{\mathrm{d}\xi }{\mathrm{d}x}=\frac{\mathrm{1}}{l_{\mathrm{s}}}\left(\mathrm{6}\xi -\mathrm{6}{\xi }^{\mathrm{2}}\right)v+\left(\mathrm{3}{\xi }^{\mathrm{2}}-\mathrm{2}\xi \right)\theta$$

（7）

则各截面处的曲率$\kappa \left(\xi \right)$和轴向变形率$e\left(\xi \right)$分别如式（8）、式（9）所示。

$$\kappa \left(\xi \right)=\frac{\mathrm{d}\theta \left(\xi \right)}{\mathrm{d}\xi }\frac{\mathrm{d}\xi }{\mathrm{d}x}=\frac{\mathrm{1}}{l_{\mathrm{s}}^{\mathrm{2}}}\left(\mathrm{6}-\mathrm{12}\xi \right)v+\frac{\mathrm{1}}{l_{\mathrm{s}}}\left(\mathrm{6}\xi -\mathrm{2}\right)\theta$$

（8）

$$e\left(\xi \right)=\frac{\mathrm{d}u\left(\xi \right)}{\mathrm{d}\xi }\frac{\mathrm{d}\xi }{\mathrm{d}x}=\frac{u}{l_{\mathrm{s}}}$$

（9）

将$\xi =x/l_{\mathrm{s}}$代入得

$$\kappa \left(x\right)=\frac{\mathrm{1}}{l_{\mathrm{s}}^{\mathrm{2}}}\left(\mathrm{6}-\mathrm{12}\frac{x}{l_{\mathrm{s}}}\right)v+\frac{\mathrm{1}}{l_{\mathrm{s}}}\left(\mathrm{6}\frac{x}{l_{\mathrm{s}}}-\mathrm{2}\right)\theta$$

（10）

$$e\left(x\right)=\frac{u}{l_{\mathrm{s}}}$$

（11）

将u、v、θ作用下的分析区域变形分解为常曲率变形和变曲率变形两部分，分解后常曲率变形由常曲率的转动变形和轴向变形组成，变曲率变形为线性变化曲率的转动变形，如图 3 所示。

对于常曲率变形部分，曲率变化率应为零，即

$$\frac{\mathrm{d}\kappa \left(x\right)}{\mathrm{d}x}=\mathrm{0}$$

（12）

可得

$$v=\frac{l_{\mathrm{s}}\theta }{\mathrm{2}}$$

（13）

则常曲率变形和变曲率变形对应的位移量为

$$U_{\mathrm{p}}={\left[u,\frac{l_{\mathrm{s}}\theta }{\mathrm{2}},\theta \right]}^{\mathrm{T}}$$

（14）

$$U_{\mathrm{s}}={\left[\mathrm{0},v-\frac{l_s\theta }{\mathrm{2}},\mathrm{0}\right]}^{\mathrm{T}}$$

（15）

$$U=U_{\mathrm{p}}+U_{\mathrm{s}}$$

（16）

式中：$U_{\mathrm{p}}$为常曲率变形部分的位移量；$U_{\mathrm{s}}\mathrm{为}$变曲率变形部分的位移量；$U$为分析区域总的相对位移量。

将式（14）代入式（10），则常曲率变形的曲率函数为

$${\kappa }_{\mathrm{p}}(x)=\frac{\theta }{l_{\mathrm{s}}}$$

（17）

将式（15）代入式（10），则变曲率变形的曲率函数为

$${\kappa }_{\mathrm{s}}\left(x\right)=\frac{\mathrm{1}}{l_{\mathrm{s}}^{\mathrm{2}}}\left(\mathrm{6}-\mathrm{12}\frac{x}{l_{\mathrm{s}}}\right)(v-\frac{l_{\mathrm{s}}\theta }{\mathrm{2}})$$

（18）

1.3.2　物理方程及接头不连续条件

根据环间接缝受压侧仅管片受压、螺栓不参与受力的假定，可知长为$l_{\mathrm{s}}$的分析区内其抗压刚度可表示为

$$k_{\mathrm{c}}=E_{\mathrm{c}}t$$

（19）

受拉侧由管片混凝土和螺栓共同受拉，则在长为$l_{\mathrm{s}}$的分析区内其抗拉刚度可表示为

$$k_t=\zeta E_{\mathrm{c}}t$$

（20）

$$\zeta =\frac{k_j}{E_{\mathrm{c}}t+l_{\mathrm{s}}{k}_j}$$

（21）

式中：$E_{\mathrm{c}}$为管片弹性模量；t为管片厚度；$\zeta$为分析区拉压刚度比；$k_j=n{K}_j/\left(\mathrm{2}\mathrm{\pi }r\right)$，$K_j$为单根螺栓沿隧道纵向的抗拉刚度，$n$为环间接头螺栓总个数。从上述抗压和抗拉刚度定义可以发现，该刚度为沿盾构隧道横截面圆周方向的线性刚度。

由于环间接缝不连续，导致隧道沿圆周方向上的受拉刚度和抗压刚度不同。令某截面y处应变为$\epsilon \left(y\right)$，则考虑接头不连续性后的截面分布力为

$$\sigma \left(\phi \right)=\epsilon \left(y\right)\left\{\begin{array}{c}\begin{array}{cc}{k}_{\mathrm{t}}& \epsilon \left(y\right)\ge \mathrm{0}\end{array}\\ \begin{array}{cc}{k}_{\mathrm{c}}& \epsilon \left(y\right)<\mathrm{0}\end{array}\end{array}\right.$$

（22）

显然，考虑了接头不连续性后的应力函数为分段函数，这使得内力计算需要分状态进行，计算较为繁琐，因此将其改写为式（23）。在该式中，当$\epsilon \left(y\right)$为正时$k_{\mathrm{c}}$项被约去，为负值时$k_{\mathrm{t}}$项被约去，因此该式与式（22）是等效的。

$$\sigma \left(y\right)=\frac{\epsilon \left(y\right)}{\mathrm{2}}\left(k_{\mathrm{t}}+k_{\mathrm{c}}\right)+\frac{\left|\epsilon \left(y\right)\right|}{\mathrm{2}}\left(k_{\mathrm{t}}-k_{\mathrm{c}}\right)$$

（23）

1.3.3　盾构隧道纵向力学方程推导

由1.3.1节分析区域变形函数的分解，分别推导常曲率变形和变曲率变形的力学方程，然后将其组合得到盾构隧道纵向力学方程。

1.3.3.1　常曲率变形力学方程

由式（11）和式（17），可得截面上各点的平均应变函数为

$$\epsilon \left(y\right)=e(x)-{\kappa }_{\mathrm{p}}(x)y=\frac{u}{l_{\mathrm{s}}}-\frac{\theta }{l_{\mathrm{s}}}y$$

（24）

则考虑了接头不连续条件后的分布力函数为

$$\sigma \left(\phi \right)=\frac{\epsilon \left(\phi \right)}{\mathrm{2}}\left(k_{\mathrm{t}}+k_{\mathrm{c}}\right)+\frac{\left|\epsilon \left(\phi \right)\right|}{\mathrm{2}}\left(k_{\mathrm{t}}-k_{\mathrm{c}}\right)$$

（25）

$$\epsilon \left(\phi \right)=\frac{u}{l_{\mathrm{s}}}-r\mathrm{s}\mathrm{i}\mathrm{n}\phi \frac{\theta }{l_{\mathrm{s}}}$$

（26）

式（26）采用了极坐标表达方式，其中φ的含义如图 4 所示。

图4 积分变量φ

Fig. 4 Integral variable φ

$\mathrm{如}\mathrm{图}\mathrm{4}\mathrm{所}\mathrm{示},\mathrm{在}{U}_{\mathrm{p}}$作用下，面2上作用u、θ位移时的作用轴力和弯矩分别为

$$N_{\mathrm{p}}=\mathrm{2}{\int }_{-\frac{\mathrm{\pi }}{\mathrm{2}}}^{\frac{\mathrm{\pi }}{\mathrm{2}}}\sigma \left(\phi \right)·r\mathrm{d}\phi$$

（27）

$$M_{\mathrm{p}}=-\mathrm{2}{\int }_{-\frac{\mathrm{\pi }}{\mathrm{2}}}^{\frac{\mathrm{\pi }}{\mathrm{2}}}\sigma \left(\phi \right)r\mathrm{s}\mathrm{i}\mathrm{n}\phi ·r\mathrm{d}\phi$$

（28）

其中，下标p表示该力为$U_{\mathrm{p}}$作用下的力。将式（25）、（26）代入式（27）、（28）得

$$N_{\mathrm{p}}=N_{\mathrm{p}\mathrm{1}}+N_{\mathrm{p}\mathrm{2}}$$

（29）

$$M_{\mathrm{p}}={\mathrm{M}}_{\mathrm{p}\mathrm{1}}+M_{\mathrm{p}\mathrm{2}}$$

（30）

$$N_{\mathrm{p}\mathrm{1}}=\frac{r\left(k_{\mathrm{t}}+k_{\mathrm{c}}\right)}{l_{\mathrm{s}}}{\int }_{-\frac{\mathrm{\pi }}{\mathrm{2}}}^{\frac{\mathrm{\pi }}{\mathrm{2}}}\left(u-\theta r\mathrm{s}\mathrm{i}\mathrm{n}\phi \right)\mathrm{d}\phi$$

（31）

$$M_{\mathrm{p}\mathrm{1}}=-\frac{r^{\mathrm{2}}\left(k_{\mathrm{t}}+k_{\mathrm{c}}\right)}{l_{\mathrm{s}}}{\int }_{-\frac{\mathrm{\pi }}{\mathrm{2}}}^{\frac{\mathrm{\pi }}{\mathrm{2}}}(u-\theta r\mathrm{s}\mathrm{i}\mathrm{n}\phi )\mathrm{s}\mathrm{i}\mathrm{n}\phi \mathrm{d}\phi$$

（32）

$$N_{\mathrm{p}\mathrm{2}}=\frac{r\left(k_{\mathrm{t}}-k_{\mathrm{c}}\right)}{l_{\mathrm{s}}}{\int }_{-\frac{\mathrm{\pi }}{\mathrm{2}}}^{\frac{\mathrm{\pi }}{\mathrm{2}}}\left|u-\theta r\mathrm{s}\mathrm{i}\mathrm{n}\phi \right|\mathrm{d}\phi$$

（33）

$$M_{\mathrm{p}\mathrm{2}}=-\frac{r^{\mathrm{2}}\left(k_{\mathrm{t}}-k_{\mathrm{c}}\right)}{l_{\mathrm{s}}}{\int }_{-\frac{\mathrm{\pi }}{\mathrm{2}}}^{\frac{\mathrm{\pi }}{\mathrm{2}}}\left|u-\theta r\mathrm{s}\mathrm{i}\mathrm{n}\phi \right|\mathrm{s}\mathrm{i}\mathrm{n}\phi \mathrm{d}\phi$$

（34）

其中，式（31）、（32）积分有解析解，式（33）、（34）为超越积分函数。因此，式（33）、（34）引入辛普森数值积分方法，可得$U_{\mathrm{p}}$作用下的力为

$$N_{\mathrm{p}}=\frac{\mathrm{\pi }r\left(k_{\mathrm{t}}+k_{\mathrm{c}}\right)}{l_{\mathrm{s}}}u+\frac{r\left(k_{\mathrm{t}}-k_{\mathrm{c}}\right)}{l_{\mathrm{s}}}\frac{\mathrm{\pi }}{\mathrm{6}}(\left|u-\theta r\right|+\left|u+\theta r\right|+\mathrm{4}\left|u\right|)$$

（35）

$$M_{\mathrm{p}}=r^{\mathrm{3}}\frac{k_{\mathrm{t}}+k_{\mathrm{c}}}{l_{\mathrm{s}}}\frac{\mathrm{\pi }}{\mathrm{2}}\theta -\frac{r^{\mathrm{2}}\left(k_{\mathrm{t}}-k_{\mathrm{c}}\right)}{l_{\mathrm{s}}}\frac{\mathrm{\pi }}{\mathrm{6}}\left(\left|u-\theta r\right|-\left|u+\theta r\right|\right)$$

（36）

令

$$e_u=\frac{\left|\theta \right|r}{u}$$

（37）

其物理含义为，转动θ与轴向u在分析区上缘点处引起的位移之比，称为位移偏心距。

再将式（19）、（21）代入式（35）、（36），得

$$N_{\mathrm{p}}=\frac{E_{\mathrm{c}}A}{l_{\mathrm{s}}}{\lambda }_{\mathrm{1}}u$$

（38）

$$M_{\mathrm{p}}=\frac{E_{\mathrm{c}}I}{l_{\mathrm{s}}}{\lambda }_{\mathrm{2}}\theta$$

（39）

$$\begin{array}{l}{\lambda }_{\mathrm{1}}=\frac{\zeta +\mathrm{1}}{\mathrm{2}}+e_u\frac{\zeta -\mathrm{1}}{\mathrm{12}}(\left|\frac{\mathrm{1}}{e_u}-\mathrm{1}\right|+\\ \left|\frac{\mathrm{1}}{e_u}+\mathrm{1}\right|+\mathrm{4}\left|\frac{\mathrm{1}}{e_u}\right|)\end{array}$$

（40）

$${\lambda }_{\mathrm{2}}=\frac{\zeta +\mathrm{1}}{\mathrm{2}}-\frac{\zeta -\mathrm{1}}{\mathrm{6}}(\left|\frac{\mathrm{1}}{e_u}-\mathrm{1}\right|-\left|\frac{\mathrm{1}}{e_u}+\mathrm{1}\right|)$$

（41）

式中：${\lambda }_{\mathrm{1}}$、${\lambda }_{\mathrm{2}}$为与位移偏心距和拉压刚度比相关的系数，称轴向刚度折减系数和弯曲刚度折减系数；A为管片环横截面的面积，A=$\mathrm{2}\mathrm{\pi }rt$；I为管片环的近似惯性矩$,I=\mathrm{\pi }{r}^{\mathrm{3}}t$。

由式（38）、（39）即可求得$U_{\mathrm{p}}$作用时的轴力和弯矩项。此外，盾构隧道纵向力学分析需要分析纵向接头受压区高度、环间接缝张开量、接触压力等力学响应。

若令环间接缝的上部和下部的变形量分别为$d_{\mathrm{1}}$和$d_{\mathrm{2}}$，则其表达式为

$$d_{\mathrm{1}}=\frac{u-\theta r}{l_{\mathrm{s}}}$$， $$d_{\mathrm{2}}=\frac{u+\theta r}{l_{\mathrm{s}}}$$

（42）

由$d_{\mathrm{1}}$、$d_{\mathrm{2}}$可以判断环间接缝的接触状态：①当$d_{\mathrm{1}}>\mathrm{0}$，$d_{\mathrm{2}}>\mathrm{0}$时，为无接触；②当$d_{\mathrm{1}}\le \mathrm{0}$，$d_{\mathrm{2}}\le \mathrm{0}$时，为全接触；③当$d_{\mathrm{1}}\le \mathrm{0}$，$d_{\mathrm{2}}>\mathrm{0}$时，为上部接触；④当$d_{\mathrm{1}}>\mathrm{0}$，$d_{\mathrm{2}}\le \mathrm{0}$时，为下部接触。

对于状态③上部接触和④下部接触，曲率${\kappa }_{\mathrm{p}}(x)$表示该截面处的转角变化率，$e(x)$为该截面的轴向变形率。图 5 所示为接缝下部接触示意图，所以接缝截面的受压区和受拉区的临界坐标可表示为

$$h_{\mathrm{c}}=\frac{e(x)}{{\kappa }_{\mathrm{p}}(x)}=\frac{u}{\theta }$$

（43）

图5 盾构隧道环间接缝接触示意

Fig. 5 Schematic diagram of contact in longitudinal shield tunnel

相应的极坐标圆心角为

$${\phi }_{\mathrm{c}}=\mathrm{s}\mathrm{i}{\mathrm{n}}^{-\mathrm{1}}\left(\frac{h_{\mathrm{c}}}{r}\right)$$

（44）

式中：$h_{\mathrm{c}}$以向y轴正方向为正；${\phi }_{\mathrm{c}}$以自z轴向y轴方向为正。

环间接缝的上沿和下沿纵向应力由式（45）、（46）求得

$$p_{\mathrm{1}}=\frac{\sigma \left(r\right)}{t}=$$

$$\frac{\left(u-\theta r\right)\left(k_t+k_{\mathrm{c}}\right)+\left|u-\theta r\right|\left(k_{\mathrm{t}}-k_{\mathrm{c}}\right)}{\mathrm{2}{l}_{\mathrm{s}}}$$

（45）

$$p_{\mathrm{2}}=\frac{\sigma \left(r\right)}{t}=$$$$\frac{\left(u+\theta r\right)\left(k_{\mathrm{t}}+k_{\mathrm{c}}\right)+\left|u+\theta r\right|\left(k_{\mathrm{t}}-k_{\mathrm{c}}\right)}{\mathrm{2}{l}_{\mathrm{s}}}$$

（46）

式中：$p_{\mathrm{1}}$为上沿纵向应力；$p_{\mathrm{2}}$为下沿纵向应力。

由式（40）、（41）可得到环间接缝上沿和下沿的张开量为

$$o_{\mathrm{1}}=\mathrm{m}\mathrm{a}\mathrm{x}(\frac{p_{\mathrm{1}}t}{k_j},\mathrm{0})$$

（47）

$$o_{\mathrm{2}}=\mathrm{m}\mathrm{a}\mathrm{x}(\frac{p_{\mathrm{1}}t}{k_j},\mathrm{0})$$

（48）

1.3.3.2　变曲率变形力学方程及组合方程

对于变曲率变形部分，轴向变形量为零，又由环缝对应的x坐标为$l_{\mathrm{s}}/\mathrm{2}$，代入式（18）得环缝截面曲率为

$${\kappa }_{\mathrm{s}}\left(x=\frac{l_{\mathrm{s}}}{\mathrm{2}}\right)=\mathrm{0}$$

（49）

可以发现，按常曲率变形和变曲率变形分解后，常曲率部分在环缝面位置处曲率和轴向变形率不为零，由此决定了接缝面的接触状态。变曲率变形在环缝位置处的曲率函数为零，轴向变形率亦为零，环缝面不产生相对转动和轴向变形，不影响接缝面的法向接触状态，因此变曲率变形的刚度性质应与常曲率变形相同。

由此，变曲率变形部分剪力项和弯矩项为

$$\begin{array}{l}{V}_{\mathrm{s}}=-E_{\mathrm{c}}I{\lambda }_{\mathrm{2}}{\left.\frac{{\kappa }_{\mathrm{s}}\left(x\right)}{\mathrm{d}x}\right|}_{x=l_{\mathrm{s}}}=\\ \mathrm{12}\frac{E_{\mathrm{c}}I{\lambda }_{\mathrm{2}}}{l_{\mathrm{s}}^{\mathrm{3}}}v-\mathrm{6}\frac{E_{\mathrm{c}}I{\lambda }_{\mathrm{2}}}{l_{\mathrm{s}}^{\mathrm{2}}}\theta \end{array}$$

（50）

$$\begin{array}{l}{M}_{\mathrm{s}}=E_{\mathrm{c}}I{\lambda }_{\mathrm{2}}{\left.{\kappa }_{\mathrm{s}}(x)\right|}_{x=l_{\mathrm{s}}}=\\ -\mathrm{6}\frac{E_{\mathrm{c}}I{\lambda }_{\mathrm{2}}}{l_{\mathrm{s}}^{\mathrm{2}}}v+\mathrm{3}\frac{E_{\mathrm{c}}I{\lambda }_{\mathrm{2}}}{l_{\mathrm{s}}}\theta \end{array}$$

（51）

式中，$E_{\mathrm{c}}I{\lambda }_{\mathrm{2}}$为常曲率变形下的盾构隧道纵向的等效弯曲刚度。

将$U_{\mathrm{p}}$与$U_{\mathrm{s}}$作用下产生的力对应求和，可得截面2施加$U=(u,v,\theta )$位移时的力为

$$N=N_{\mathrm{p}}=\frac{E_{\mathrm{c}}A}{l_{\mathrm{s}}}{\lambda }_{\mathrm{1}}u$$

（52）

$$V=V_{\mathrm{s}}=\mathrm{12}\frac{E_{\mathrm{c}}I{\lambda }_{\mathrm{2}}}{l_{\mathrm{s}}^{\mathrm{3}}}v-\mathrm{6}\frac{E_{\mathrm{c}}I{\lambda }_{\mathrm{2}}}{l_{\mathrm{s}}^{\mathrm{2}}}\theta$$

（53）

$$M=M_{\mathrm{p}}+M_{\mathrm{s}}=-\mathrm{6}\frac{E_{\mathrm{c}}I{\lambda }_{\mathrm{2}}}{l_{\mathrm{s}}^{\mathrm{2}}}v+\mathrm{4}\frac{E_{\mathrm{c}}I{\lambda }_{\mathrm{2}}}{l_{\mathrm{s}}}\theta$$

（54）

写为矩阵形式得

$$\mathit{F}=\mathit{K}\mathit{U}$$

（55）

其中

$$\mathit{K}=\left[\begin{array}{ccc}\frac{E_{\mathrm{c}}A}{l_{\mathrm{s}}}{\lambda }_{\mathrm{1}}& \mathrm{0}& \mathrm{0}\\ \mathrm{0}& \mathrm{12}\frac{E_{\mathrm{c}}I}{l_{\mathrm{s}}^{\mathrm{3}}}{\lambda }_{\mathrm{2}}& -\mathrm{6}\frac{E_{\mathrm{c}}I}{l_{\mathrm{s}}^{\mathrm{2}}}{\lambda }_{\mathrm{2}}\\ \mathrm{0}& -\mathrm{6}\frac{E_{\mathrm{c}}I}{l_{\mathrm{s}}^{\mathrm{2}}}{\lambda }_{\mathrm{2}}& \mathrm{4}\frac{E_{\mathrm{c}}I}{l_{\mathrm{s}}}{\lambda }_{\mathrm{2}}\end{array}\right]$$

（56）

式（55）、（56）即为盾构隧道纵向的力学控制方程。可以发现，该式实际上与轴向刚度为$E_{\mathrm{c}}A{\lambda }_{\mathrm{1}}$、弯曲刚度为$E_{\mathrm{c}}I{\lambda }_{\mathrm{2}}$的梁构件控制方程相同，因此可采用等效梁模拟其力学性质，其轴向和弯曲刚度为考虑折减系数${\lambda }_{\mathrm{1}}$和${\lambda }_{\mathrm{2}}$的刚度。

2 蠕滑作用下盾构隧道纵向力学行为分析

2.1　工程背景

汕头至汕尾铁路途经的汕头湾海底隧道位于广东省汕头市濠江区和龙湖区境内，线路呈东北—西南走向，属控制性重点隧道工程。工程场区主要有丘陵地貌、海底地貌及平原地貌，隧道全长9 771 m，其中越岭段长约5 741 m，穿汕头湾海底段约2 000 m，穿三角洲平原段约2 030 m，最大埋深约180 m，汕头湾海底隧道纵断面如图 6 所示。

图6 汕头湾海底隧道纵断面

Fig. 6 Transverse section of tunnel in Shantou Bay

盾构隧道主体结构拟采用单层装配式管片衬砌，内径12.8 m，外径14.0 m，管片环采用“9+1、1/3封顶块”型式，由1个封顶块（F）、2个邻接块（L1、L2）和7个标准块（B1~B7）组成，平均幅宽2.0 m。混凝土强度等级为C60，抗渗等级为P12，块与块之间每块设置3颗M36环向斜螺栓，螺栓机械等级10.9级。环与环之间每环布置56颗M36纵向斜螺栓，螺栓机械等级10.9级。螺栓中心线与管片接缝交点位于0.5倍厚度处，螺栓中心线与管片径向线近似呈60°，盾构隧道横断面如图 7 所示。

图7 汕头湾海底盾构隧道横断面

Fig.7 Cross-section of tunnel in Shantou Bay

盾构段隧道横穿桑浦山活动断裂带，其平面关系如图 8 所示。该断层为正断层，宽约20 m。根据地勘资料及活动断裂带工程活动性专题研究报告，针对汕头湾海底隧道近场新建GNSS（全球导航卫星系统）的观测资料，通过解算得到的速度场分布结果知断层主要产生垂直运动，推断断层垂直错动速率为（0.15±0.44） mm·a^-1。

图8 汕头湾海底隧道与断层关系平面

Fig. 8 Relationship between tunnel and fault in Shantou Bay

2.2　蠕滑作用下盾构隧道力学分析方法

采用图 9 所示流程分析蠕滑作用下盾构隧道的力学行为，主要流程如下：

图9 蠕滑作用下盾构隧道纵向力学分析流程

Fig. 9 Process of longitudinal mechanical analysis of shield tunnel under creep action

（1）建立蠕滑位错作用下的地层宏观计算模型。汕头湾隧道地层计算模型采用连续介质地层计算模型，如图 10 所示，地层物理力学参数见表 1。通过计算得到隧道所在位置地层变形。

图10 蠕滑作用下地层计算模型

Fig. 10 Soil calculation model under creep action

表1 地层物理力学参数

Tab. 1 Mechanical parameters of soil

地层编号	岩土名称	天然重度/ （kN·m-3）	泊松比	变形模量/MPa	黏聚力/kPa	内摩擦角/（°）	层厚/m
①	粉砂	18.5	0.30	10		23.00	2.5
②	淤泥质土	16.5	0.45	2	6.23	7.14	4.0
③	全风化花岗岩	18.8	0.37	50	29.17	20.99	9.5
④	强风化花岗岩	21.0	0.35	500	100	21.00	49.0
⑤	弱风化花岗岩	26.7	0.25	83 700	13 300	42.90	55.0
⑥	断碎带区域	16.4	0.42	100	50	21.50

（2）建立盾构隧道纵向节点-接头结构计算模型。将隧道所在位置地层变形作为荷载输入，如图 11 所示，为考虑接头非线性，采用如图 12 所示流程计算。

图11 盾构隧道纵向节点-接头结构模型

Fig. 11 Calculation model of longitudinal node-joint of shield tunnel

图12 盾构隧道纵向节点-接头结构模型计算流程

Fig. 12 Process of longitudinal structure analysis of shield tunnel under creep action

（3）提取结构内力和节点位移，由节点位移计算纵向接头力学响应。

汕头湾盾构隧道结构计算采用如图 11 所示的节点-接头结构计算模型，其中接头采用轴向刚度为$E_{\mathrm{c}}A{\lambda }_{\mathrm{1}}$、弯曲刚度为$E_{\mathrm{c}}I{\lambda }_{\mathrm{2}}$的等效梁元件模拟，梁长与管片环幅宽相同。计算参数列于表 2，其中r为隧道计算半径，t为衬砌厚度，$l_{\mathrm{s}}$为管片环幅宽，$K_{\mathrm{j}}$为环间接头单根螺栓沿隧道纵向的抗拉刚度，$k_{\mathrm{t}}$为地基弹簧纵向刚度，$k_{\mathrm{n}}$为地基弹簧横向刚度。

表 2 盾构隧道结构计算参数

Tab. 2 Structural calculation parameters of shield tunnel

r/m	t/m	ls/m	Ec/GPa	$K_{\mathrm{j}}$/（MN·m-1）	Kt/（MN·m-1）	Kn/（MN·m-1）
6.7	0.6	2.0	36	486	607	1 820

2.3　计算工况

根据汕头湾隧道地勘资料及活动断裂带工程活动性专题研究报告，桑浦山断裂带主要表现为垂直错动，垂直错动速率为（0.15±0.44）mm·a^-1，因此计算工况选取了蠕滑速率介于（0.15~0.59）mm·a^-1的3种情况，蠕滑时间取介于20~120 a的5种情况，具体工况见表 3。

表 3 计算工况

Tab. 3 Calculation condition

编号	蠕滑速率/ （mm·a-1）	时间/a
1	0.59	20、40、60 80、100
2	0.30
3	0.15

将蠕滑位移作为荷载，输入图 10 的地层计算模型，得到隧道轴向位置处的竖向（Y轴向）和水平向（X轴向）地层位移。图13 为蠕滑速率为0.59 mm·a^-1 （工况1）时隧道位置处的地层位移，其中隧道与断裂带交叉区段的隧道纵轴向坐标为191~212m。可以看出，由于断裂带区域地层刚度相对较小，当断层上下盘发生相对位移时，断裂带区域变形最大。

图13 蠕滑速率为0.59 mm·a^-1隧道位置处地层位移

Fig. 13 Ground displacement at tunnel location at a creep rate of 0.59 mm·a^-1

2.4　盾构隧道结构力学响应

图 14 为计算所得的盾构隧道150~250 m段的结构内力图，可以看出：①桑浦山断裂为正断层，因此隧道结构表现为轴力为拉力，远离断层带时轴力先增大后减小，表现出上盘隧道底部受拉、下盘隧道顶部受拉；②不同蠕滑速率下隧道结构内力分布形式基本不变，仅量值大小发生变化，以0.15 mm·a^-1和0.59 mm·a^-1蠕滑为例，设计使用100 a时内力增大约3倍。

图14 考虑接头非线性时盾构隧道纵向内力响应

Fig. 14 Longitudinal internal force of shield tunnel considering nonlinear stiffness

图 15 比较了采用常数刚度和非线性刚度情况

图15 采用非线性和常数刚度时盾构隧道纵向内力响应（0.30 mm·a^-1， 100 a）

Fig. 15 Longitudinal internal force of shield tunnel with nonlinear and constant stiffness (0.30 mm·a^-1, 100 a)

下结构内力响应的差异，图中${\lambda }_{\mathrm{1}}$、${\lambda }_{\mathrm{2}}$分别为采用常数轴向刚度和常数弯曲刚度分析时的刚度折减系数，可以看出：①轴向刚度对轴力和剪力影响大，对弯矩影响小，采用常数轴向刚度分析时，轴力分布影响大，弯矩分布基本无影响，轴力分布由非线性变为接近线性，相差最大时，非线性刚度下的轴力是常数轴向刚度下的3.3倍；②弯曲刚度对弯矩影响较大，对轴力和剪力影响较小，弯曲刚度对弯矩分布形式基本无影响，相差最大时，非线性刚度下的弯矩是常数弯曲刚度下的1.4倍；③蠕滑作用下，采用常数刚度与非线性刚度时的盾构隧道结构内力差异较大，考虑内力共同作用下的接触非线性是有必要的。

3 结论

盾构隧道因由管片环拼接而成，其纵向力学性能表现出显著的接触非线性现象，如在活动断裂带蠕滑作用下，其环间接缝变形与传统的将轴向、弯矩性能独立考虑的情况有较大差异，因此本文主要针对考虑轴力、剪力、弯矩共同作用下环间接缝接触非线性的盾构隧道的纵向力学分析方法做了研究，在此基础上分析了蠕滑作用下盾构隧道结构内力响应，主要有如下结论：

（1）针对盾构隧道纵向接头的非线性接触变形问题，提出了一种盾构隧道纵向分析模型，该模型能考虑盾构隧道纵向接头同时传递轴力、剪力、弯矩时的接触非线性，可以更精确地分析盾构隧道在复杂载荷作用下的力学响应。

（2）当考虑盾构隧道纵向同时传递轴力、剪力、弯矩时，其传力作用可等效为考虑轴向和转动刚度折减的梁单元，轴向和转动刚度折减系数与位移偏心距有关，可由位移偏心距计算梁的等效刚度。

（3）汕头湾盾构隧道在蠕滑作用下的力学分析表明，在不同蠕滑量时，隧道纵向内力分布规律相同，纵向接头非线性对隧道结构内力响应影响较大，考虑接头非线性时轴力可达到常数刚度下的3.3倍，弯矩可达到常数刚度下的1.4倍，考虑接头的非线性刚度是有必要的。

作者贡献声明

耿 萍：制定研究方案和论文撰写，全文审阅。

陈枰良：理论推导，研究框架，论文撰写。

陈昌健：非线性算法实现，论文撰写。

杨 琪：数值模拟，论文审阅及修改。

郭翔宇：提供研究思路和技术指导。

何 川：提供研究支持，论文审阅。