全断面隧道掘进机(tunnel boring machine, TBM)凭借其安全、高效的特点在各种地下工程施工中得到了广泛的应用，尤其对于深埋长隧洞，TBM的优势更加明显^[1]。TBM开挖岩体后，围岩与护盾之间存在一个较小的间隙，当围岩的收敛变形量过大时，便会挤压护盾而产生摩阻力，当摩阻力大于TBM的最大推进力时，就会引发卡机的风险。对于深部高地应力岩体，围岩的时效变形明显，卡机事故频发。卡机不仅影响施工进程，也会增加工程的成本^[2-4]。因此，准确分析围岩与TBM的相互作用，评价卡机风险，并研究卡机控制措施的影响具有重要意义。

卡机风险分析的本质是分析围岩与TBM的相互作用机制，方法包括经验性分析模型、理论分析模型和数值仿真分析模型3种。早期经验性分析模型主要基于实际工程经验，结合图表对卡机进行定性分析，或通过构建一些指标来判断围岩的挤压变形程度，从而间接评价卡机风险。如：Ramoni等^[5]根据N²关系图表法定性分析了TBM掘进速度、停机时间对卡机的影响。Aydan等^[6]根据日本某隧洞工程的经验，提出通过完整岩石强度与上覆围岩压力的比值来评价隧道的挤压风险。理论分析模型主要基于弹塑性力学理论对隧洞围岩的应力状态和收敛变形进行解析求解。例如，Carranza-Torres等^[7]根据收敛-约束法提出了挤压变形的评价方法；温森等^[8]基于Hoek-Brown准则和收敛-约束法建立了停机和连续掘进2种工况下卡机计算模型；Zhang等^[9]着重考虑了地层变形的时间效应，推导了挤压性地层卡机计算模型。由于TBM卡机的影响因素众多，包括围岩特性(抗压强度、弹性模量)、TBM掘进速度、超挖量、地应力、隧洞直径等。经验性分析模型和理论分析模型的考虑因素均较少，存在一定的局限性。近年来，数值仿真技术成为分析岩土工程问题的重要手段，并取得了很多成果。Hasanpour等^[10]基于FLAC^3D(fast Lagrangian analysis of continua in 3 dimensions)提出了不良地质条件下的双护盾TBM掘进分析模型，并对TBM掘进过程中的卡机风险进行了分析；黄兴等^[11]提出了可用于模拟TBM开挖围岩挤压大变形的本构模型，对引大济湟隧道工程卡机段进行了数值计算和分析；程建龙等^[12]采用FLAC^3D建立了复合地层三维双护盾TBM开挖模型，研究了不同地层TBM与围岩的作用机制。

本文以某公路隧洞工程为背景，建立考虑双护盾TBM施工过程的数值仿真模型，通过基于内变量热力学的蠕变损伤模型模拟围岩的时效变形，通过FLAC^3D中的接触面单元，对围岩-护盾相互作用进行模拟。在模型中考虑不同的超挖量、TBM推进速度2个卡机控制参数，通过设置不同的参数组合研究了不同地应力和不同围岩弹性模量下，2个卡机控制参数的取值对围岩-护盾接触压力及TBM卡机的影响。

1 工程概况与卡机情况 1.1 工程概况

隧洞总长度4 789 m，隧洞施工主要采用德国海瑞克生产的双护盾TBM掘进，开挖洞径9.13 m，衬砌采用“6+1”型四边形预制混凝土管片，管片厚度35 cm，衬砌后内径8.1 m。隧道沿线大多基岩裸露，穿越地层为元古界南迦巴瓦岩群多雄拉岩组(Pt2-3d)，属喜马拉雅地层区，岩性单一，主要由花岗片麻岩、混合片麻岩等组成，片麻理发育。

该隧洞为越岭深埋隧道，其中垂直埋深大于600 m的洞段超过2 km，最大埋深达830 m，通过现场深钻孔水压致裂法地应力测试可推算出隧洞中心最大水平应力分量约为30 MPa，最大竖向应力分量约为20 MPa，竖向与横向地应力比约为2∶3。现场围岩以Ⅲ级为主，部分段落为Ⅳ级围岩，局部存在Ⅴ级围岩。根据地质勘探，隧洞施工区域岩体的物理力学参数如表 1所示。

1.2 卡机情况介绍

卡机分为卡刀盘和卡护盾2种^[13]，其中刀盘被卡主要发生在掌子面处，当掌子面处围岩破碎、自稳能力差时，开挖扰动后在高地应力作用下，掌子面处围岩破碎坍塌而导致刀盘被卡^[14]，通常卡刀盘可采用增大刀盘扭矩和增大超前注浆等方法进行脱困处理。TBM卡护盾在深埋高地应力等复杂地质条件区域出现较多，具有更大的危害。该隧洞双护盾TBM掘进施工过程中，先后出现了4次卡机事故，均为护盾被卡。图 1为工程现场2次卡机护盾不同部位围岩挤压情况。

图中2次卡机发生区域均存在一定的不良地质情况，如K10+145桩号位置，发现多条沿裂隙发育的挤压带；K10+209位置岩体整体性较好，但处于围岩由差变好的交界部位。另外，2次卡机段均处于深埋高地应力区域，通过对观测窗部位观察，围岩收敛变形大，岩体与盾体基本接触，初步分析认为，高地应力作用造成的围岩收敛变形仍是引发卡机的主要原因。

2 TBM掘进过程及围岩-护盾相互作用模拟 2.1 计算模型

建立了考虑TBM掘进施工过程的数值仿真模型，模型网格如图 2所示。为了减小边界效应，模型的横向模拟范围取不小于15倍洞径，纵向模拟范围取不小于10倍洞径^[15]，计算模型范围为横向150 m×150 m，纵向100 m。同时为了使模拟更加精细，对隧洞开挖区域半径15 m范围内的网格进行加密处理，以获得更高的应力应变梯度。

模型对双护盾TBM的主要部件进行了模拟，包括刀盘、前护盾、后护盾、衬砌管片和豆砾石回填灌浆层等。其中刀盘采用实体单元，前后护盾、衬砌管片、回填灌浆层采用cshell单元进行模拟；本文主要研究围岩与护盾的相互作用，围岩采用基于内变量热力学的蠕变本构，其余TBM部件均采用线弹性模型进行模拟。TBM主要部件的力学参数如表 2所示。采用等效密度法^[16]对TBM自重进行模拟，根据TBM总重量和各部件的体积，将TBM自重换成等效密度(即TBM总重量除以模型中各部件的总体积)。管片的参数根据工程实际情况选取；豆砾石回填灌浆层硬化后的弹性模量一般设置为1.0 GPa^[17-18]，Poisson比设置为0.3。模型中双护盾TBM几何尺寸依据工程实际情况而定，如表 3所示。

对于TBM的掘进施工，采用逐步开挖法的方式进行模拟^[19]。模型纵向的网格尺寸为1 m，每个开挖步的长度为1 m，通过调整每个开挖步的蠕变时间，来模拟不同的掘进速度。模型中，对TBM岩体开挖、TBM掘进、管片安装及豆砾石回填灌浆层的回填等过程进行了模拟。每次开挖岩体后，在掌子面施加0.255 MPa的面力模拟刀盘推力，并在进入下一个开挖步前将该掌子面的面力去除。考虑到边界效应的影响，不能将整个模型挖穿至100 m，掌子面和纵向边界应有一定距离，数值模型中TBM均计算至掘进60 m。

2.2 基于内变量热力学的蠕变本构

深部岩体开挖后，围岩变形具有明显的时间效应。本文引入了基于内变量热力学的蠕变损伤模型^[20]，该蠕变本构模型能够较好地描述岩体结构的黏塑性特性，可较好地反映蠕变的3个阶段，并基于FLAC^3D的二次开发程序接口实现了^[21]。

与黏塑性内变量λ(λ₁, λ₂)和χ共轭的热力学力为

$ f_{1}^{v p}=\sqrt{J_{2}} \text { , } $

(1)

$ f_{2}^{v p}=(1+b \chi)\left(a I_{1}+\sqrt{J_{2}}\right), $

(2)

$ f_{s}=b \lambda_{2}\left(a I_{1}+\sqrt{J_{2}}\right). $

(3)

黏塑性应变率方程为

$ \dot{\varepsilon}_{m}^{v p}=a\left[(1+b \chi) \dot{\lambda}_{2}+b \lambda_{2} \dot{\chi}\right], $

(4)

$ \dot{e}_{i j}^{\tau p}=\left[\dot{\lambda}_{1}+(1+b \chi) \dot{\lambda}_{2}+b \lambda_{2} \dot{\chi}\right] \frac{s_{i j}}{2 \sqrt{J_{2}}}. $

(5)

内变量演化方程为

$ \dot{\lambda}_{1}=\frac{1}{\eta_{p 1}}\left\langle f_{1}^{v p}-h \lambda_{1}\right\rangle ; $

(6)

$ \dot{\lambda}_{2}=\kappa_{p 2}\left\langle\frac{f_{2}^{v p}-R}{R}\right\rangle^{p} ; $

(7)

$ \dot{\chi}=\kappa_{p 3} \exp (m \chi)\left(\frac{f_{s}}{R}\right)^{2} \operatorname{sign}\left(\dot{\lambda}_{1}\right) ; $

(8)

$ \langle x\rangle=\left\{\begin{array}{ll} x, & x>0 ;\\ 0, & x \leqslant 0 ; \end{array}\right. $

(9)

$ \operatorname{sign}(x)=\left\{\begin{array}{ll} 1, & x>0 ; \\ 0, & x=0 ; \\ -1, & x<0. \end{array}\right. $

(10)

其中：λ(λ₁, λ₂)、χ为黏塑性内变量，f₁^vp、f₂^vp为与λ(λ₁, λ₂)共轭的热力学力，f_s为与内变量χ共轭的热力学力，a、b为材料参数，I₁为应力张量第一不变量，J₂为应力偏量第二不变量，ε_m^vp为黏塑性体应变，e_ij^vp为黏塑性偏应变，s_ij为应力偏量，η_p1、κ_p2和κ_p3均为黏滞系数，m为方程参数，p、h和R是材料常数，sign(x)为符号函数。

蠕变本构模型的参数中，基本力学参数通过室内试验获取，其余参数结合现场超前孔变形监测结果进行反演分析获得，超前孔及监测点(M1~M8)的布置详见文^[15]。反演分析得到围岩蠕变本构参数如表 4所示。基于该蠕变参数进行数值模拟，围岩测点M2变形实测值与计算值的对比如图 3所示。可以看出：数值模拟的围岩变形与实测围岩变形具有较好的吻合程度，说明建立的模型和设置的本构参数合理有效。本文研究中，主要考虑对弹性模量E的影响进行研究，后文中对E设置不同的值进行调整，其余参数值保持不变。

2.3 围岩-护盾相互作用模拟

TBM卡机风险分析的关键在于判断围岩与护盾接触压力的大小。在FLAC^3D数值模拟中，可通过在前后护盾的外表面施加接触面(interface)单元，计算围岩与护盾间的接触压力^[22-23]。Interface单元在岩土工程中具有较为广泛的应用，其是由一系列三节点三角形单元组成的无厚度单元，每个节点都有一个与之相关的特有面积，如图 4所示。接触面单元一般建立于实体单元(目标单元)表面，相应的实体单元四边形区域(目标面)会定义2个三角形接触面单元，并在三角形角点自动生成接触面节点，当有另外一组网格表面与接触面单元接触时，接触面节点能够识别该接触。

接触面单元的本构模型如图 5所示，图中S表示滑块元件，S_s表示抗剪强度，T_s表示抗拉强度，k_s表示切向刚度，k_n表示法向刚度，D表示剪胀角，P表示本构模型作用的接触面节点。

围岩对护盾挤压力的计算主要是计算接触面单元的法向压力，计算公式如下：

$ F_{\mathrm{n}}^{(t+\Delta t)}=k_{\mathrm{n}} u_{\mathrm{n}} A+\sigma_{\mathrm{n}} A . $

(11)

其中：F_n^(t+Δt)为t+Δt时刻接触面节点所测得的法向应力，u_n为接触面节点贯入目标面的绝对法向位移，k_n为接触面的法向刚度，A为接触面节点的代表面积，σ_n为接触面应力初始化导致的附加法向应力。

则TBM掘进所受的摩擦力可计算如下：

$ F_{\mathrm{f}}=\mu \sum\limits_{i=1}^{m} F_{\mathrm{n}}^{(i)} . $

(12)

其中：F_f为护盾上所受的总的摩擦力，F_n⁽ⁱ⁾为第i个接触面节点所测得的法向压力，m为接触面节点的个数，μ为围岩与护盾间的摩擦因数。TBM连续掘进过程中，围岩与护盾间的动摩擦因数一般为0.3~0.5；用于TBM静止后重启动的静摩擦因数一般为0.25~0.45 ^{[5, 17]}。本文中取μ的值为0.3。

图 6为前后护盾上施加接触面单元的示意图，考虑到双护盾TBM前后护盾有重叠，而重叠部分后护盾外表面不与围岩相接触，故仅在可以接触到围岩的部分施加接触面单元。

接触面单元的参数设置会影响接触压力模拟的精度，接触面的法向刚度k_n和切向刚度k_s不能过高，应低于100倍相邻单元的刚度^[25]，否则会导致FLAC^3D计算程序收敛缓慢甚至发生错误。FLAC^3D用户手册中给出了较好的经验取值方法：k_n和k_s取值为周围“最硬”相邻区域等效刚度的10倍^[24]。

$ k_{\mathrm{n}}=k_{\mathrm{s}}=10 \max \left[\frac{K+\frac{4}{3} G}{\Delta Z_{\min }}\right] . $

(13)

其中：K和G分别为体积模量和剪切模量，ΔZ_min为相邻区域法向方向单元最小尺寸。

通过一个小算例对接触面单元对接触压力的模拟效果进行验证，如图 7所示。在围岩顶部施加总大小为2 500 MN的面力，接触面单元的参数根据式(13)进行设置。可以看出，围岩与护盾间的间隙逐渐减小，直至围岩与护盾相互接触，图 8为计算得到的围岩与护盾接触面的法向应力云图，则可求得围岩与护盾间总法向压力也为2 500 MN。小算例的研究说明，当接触面单元参数取值合适时，可以精确地模拟围岩与护盾间的相互作用力。

2.4 TBM卡机状态判断

围岩开挖后，TBM掘进的阻力包括2部分：1) 由围岩与护盾挤压作用产生的摩擦力(含自重影响)；2) TBM正常连续掘进破岩所需的推进力。当上述两部分力之和大于等于TBM的推进力(额定推力)时，则会发生TBM护盾被卡的现象。故TBM卡机状态的判断条件如下：

$ \left\{\begin{array}{ll} F_{\mathrm{f}}+F_{\mathrm{e}} \geqslant F_{\mathrm{I}}, & \text { 护盾卡机 } ; \\ F_{\mathrm{f}}+F_{\mathrm{e}}<F_{\mathrm{I}}, & \text { 不发生卡机. } \end{array}\right. $

(14)

其中：F_e为TBM正常连续掘进破岩所需的推进力，F_I为TBM的最大推力(额定推力)。本文中F_e和F_I分别取值为16和61.575 MN。

3 卡机控制因素的影响分析

根据节2.4，除正常连续掘进破岩所需的推进力外，TBM卡机的摩擦力主要是由于卸荷岩体大变形对护盾产生的接触压力引起的。在TBM施工过程中，已经开挖的岩体会向洞内产生收敛变形，且越先开挖的岩体变形往往越大。推进速度对接触压力的影响主要表现为：当推进速度增大时，可在围岩收敛变形还较小时，使TBM离开较早的已开挖洞段，避免围岩与护盾的接触或减小接触面积，从而减小接触压力。超挖是减小TBM卡机风险的重要手段，通过增大TBM开挖后围岩与护盾间隙大小的方式，减小围岩与护盾的接触面积和接触压力。

3.1 TBM掘进速度的影响分析

这里不考虑TBM的停机检修等情况，认为TBM是连续掘进的，即考虑TBM的净掘进速度。根据工程经验，为研究TBM推进速度对卡机的影响，设置推进速度值分别为0.5、1.0、1.5、2.0、2.5 m/h。超挖量设置为0 cm，竖向地应力设置分别设置为20、30、40 MPa，竖向与横向地应力比与实际工程保持一致，围岩的弹性模量分别设置为2、6、10、14、18 GPa。

图 9为不同竖向地应力和围岩弹性模量下TBM推进速度对围岩-护盾接触压力的影响曲线。可以看出：

1) 随着TBM推进速度的增大，围岩与护盾的接触压力逐渐减小；且随着推进速度的增大，曲线有逐渐变缓的趋势。

2) 当围岩较软，即围岩弹性模量较小时，相邻2个数据点(相邻推进速度)连线的斜率越大。说明围岩弹性模量越小时，增加推进速度对减小接触压力(卡机风险)的收益越大。

3) 地应力越高时，相邻TBM推进速度的数据点接触压力的减小值越大；说明地应力越高时，增加推进速度对减小接触压力(卡机风险)的收益也越大。

4) 高地应力和低弹性模量下，围岩与护盾的接触压力会较大。如图 7所示，弹性模量E=2 GPa情况下：竖向地应力为20 MPa时，推进速度应达到1.5 m/h则不发生卡机；而竖向地应力为30 MPa时，推进速度达到2.5 m/h时仍会发生卡机。因此，为避免卡机或减小卡机风险，地应力越高和弹性模量越低时，则推进速度应取较大值。然而由于工程中其他施工条件的制约，TBM净推进速度通常为0.5~4 m/h，较难无限制增加。故对于高地应力和围岩变形收敛变形过大的施工条件，应尽可能增加推进速度，并同时采取其他卡机防控措施。

3.2 超挖量的影响分析

为研究超挖量对卡机的影响，分别设置超挖量为0、5、10 cm。设置推进速度值分别为0.5 m/h。竖向地应力设置分别设置为20、30、40 MPa，竖向与横向地应力比与实际工程保持一致，围岩的弹性模量分别设置为2、6、10、14、18 GPa。图 10所示为不同竖向地应力和围岩弹性模量下TBM超挖量对接触压力的影响柱状图。从图中可以看出：

1) 超挖量对减小围岩-护盾接触压力有较为显著作用，两者呈负相关关系，随着超挖量增大，围岩-护盾的接触压力减小。对于本文考虑的计算工况，当超挖量达到10 cm时，围岩接触压力均降低到较小的值或零值，不会发生卡机。

2) 围岩弹性模量越小时，增加相同的超挖量，接触压力的减小值越大。如在地应力为20 MPa情况下：弹性模量E=2 GPa时，超挖量由0增加到5 cm的接触压力减少量为86 360 kN(约41.1%)，TBM由发生卡机状态变为不发生卡机状态；弹性模量E=6 GPa时，超挖量由0增加到5 cm的接触压力减少量为33 297 kN(约48.3%)。虽然E=6 GPa相对于E=2 GPa的接触压力减少率更大，但是E=2 GPa时超挖带来的收益更大，故对于超挖对卡机的影响，应更关注对接触压力的减少量。上述分析说明，相对于硬岩(弹性模量较大)隧洞，在软岩(弹性模量较小)隧洞中增加超挖量对于减小卡机风险的收益更大。

3) 对于不同的地应力，统计E=2 GPa时，不同地应力下相邻超挖量的接触压力减小量如表 5所示。竖向地应力为20 MPa时，相邻超挖量的接触压力减小量为86 360和88 221 kN；竖向地应力为30 MPa时，相邻超挖量的接触压力减小量为179 221和180 691 kN。可以看出：同一地应力下，超挖量增加相同的增量时，接触压力的减少量近似相同，但超挖量更大时的接触压力减小量均略大；分析是由于本文计算工况围岩变形尚未达到收敛阶段，仍处于变形增长较快的阶段。对于不同地应力，高地应力下增加超挖量对减小接触压力的作用更加明显。

超挖是避免TBM卡机的有效措施，但是过分的超挖会增加刀具的磨损和施工成本。由于隧洞的自稳特性，围岩往往有一个收敛稳定值，当超挖量越接近围岩收敛变形稳定值时，增加超挖量对减小卡机风险的收益将越来越小；当超挖量大于围岩收敛变形稳定值时，增加超挖量只会增加工程成本。实际工程中，应结合具体的地应力、围岩特性、TBM推进速度等条件进行数值仿真，综合比较超挖对护盾所受摩阻力的减小效果和增加的施工成本，选取最优的超挖量值。

4 结论

本文提出了考虑双护盾TBM掘进施工过程和围岩蠕变特性的数值模型，在模型中采用基于内变量热力学的蠕变损伤模型来模拟卸荷岩体的时效变形行为，通过FLAC^3D中的接触面单元，对围岩与护盾的相互作用进行模拟，并对超挖量和TBM推进速度2个卡机控制因素的影响进行分析。主要结论如下：

1) 建立的数值模型可较为精确地模拟双护盾TBM的施工过程和围岩-护盾间的相互作用力，从而可有效地分析TBM掘进过程中的卡机风险。

2) 增大TBM推进速度能在一定程度上减小围岩与护盾的接触压力。围岩弹性模量越小、地应力越高时，增加推进速度对减小接触压力(卡机风险)的收益也越大。受其他施工条件影响，推进速度无法一直增加，通常为0.5~4 m/h，从避免卡机的角度出发，应尽可能采取较大的推进速度进行掘进施工。

3) 超挖量对减小围岩与护盾的接触压力有较为显著的作用，当超挖量达到10 cm，本文计算工况均不会发生卡机。对于高地应力、软岩(即弹性模量较小)隧洞，增加超挖量对于减小卡机风险的收益更大。但超挖量的增加会增加刀盘磨损和施工成本，且超挖量接近或超过围岩的收敛变形值后，增加超挖量的收益将变得很小，合理超挖量需要结合具体工程进行分析。