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

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

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

隧洞总长度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所示。

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

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

建立了考虑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。

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

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

黏塑性应变率方程为

内变量演化方程为