孙 博

(中国电建集团成都勘测设计研究院有限公司，四川 成都 610072)

深埋隧洞双护盾TBM施工仿真分析的关键技术初探

孙 博

(中国电建集团成都勘测设计研究院有限公司，四川 成都 610072)

长距离深埋隧洞采用双护盾TBM施工是一种主要施工方式，对其施工过程的数值仿真分析一直是岩土工程界研究的重点和难点。收集已有研究文献，归纳、总结并初步探讨数值模拟深埋隧洞双护盾TBM施工过程的关键技术。基于已有研究成果，采用FLAC3D软件及其自带FISH语言，在局部改进的基础上重新编写仿真分析程序。以某深埋隧洞为例，介绍了该仿真分析程序的计算流程，同时验证了其可行性。

深埋隧洞； 双护盾TBM； 施工过程； 数值仿真分析； 关键技术

0 前 言

目前，我国越来越多的水工隧洞和交通隧洞位于高山峡谷地区，高埋深、长距离隧洞的设计与施工已成为国内外岩土工程界的主要难题之一。针对这种深埋长隧洞，目前常用的施工方法主要分为三种：钻爆法、开敞式TBM法和双护盾TBM法。由于钻爆法施工速率慢、施工条件差、安全性差等明显缺点，近年来有逐渐被TBM法取代的趋势。虽然开敞式TBM和双护盾TBM都是采用全断面掘进机掘进施工，但两者在主机构造和支护形式上还是有较大差别的[1-3]。一般来说，开敞式TBM在应对岩爆方面较弱，对施工人员的威胁较大；而双护盾TBM在软岩大变形洞段极易出现卡机现象，延误工期。总体来看，两种施工方式各有优缺点，一般情况下采用何种施工方式主要是由地质条件决定的。根据已有的工程经验可知，当隧洞埋深超过1 000 m、甚至2 000 m时，施工过程中隧洞发生强岩爆的可能性急剧增加[4-5]。从施工人员安全角度考虑，双护盾TBM施工是一种较优的选择。因此，本文针对深埋隧洞双护盾TBM施工过程数值仿真分析的关键技术进行初步探讨，以期从数值分析上预判双护盾TBM施工过程中可能出现的各种问题，为隧洞的施工过程及支护设计提供参考。

1 研究现状

双护盾TBM作为一种主要施工方式，以其高效、安全的特性受到了广大岩土工程师的青睐。对双护盾TBM来说，其在施工中面临的最大问题是围岩大变形导致的卡机问题[6-7]。针对该问题，国内外众多学者从双护盾TBM机械构造、掘进过程及处理措施等方面进行了相关研究。从数值分析来看，现阶段主要研究成果如下：

刘泉声等[8]基于TBM卡机孕育致灾机理的理论分析，通过FLAC3D软件计算护盾区域围岩的三维应力场和位移场，从而计算出护盾前移所受的摩擦阻力，并根据卡机状态判据判断护盾是否被卡。杜立杰等[9]在进行贯入度指数FPI与主要地质因素间相关性分析的基础上，建立了TBM掘进贯入度P与FPI之间的拟合关系式，以及推进力F与FPI之间的拟合关系式。Kai Zhao等[10]在考虑复杂的围岩与护盾的相互作用、管片和豆砾石支护的基础上，采用Midas GTS软件建立了模拟双护盾TBM施工过程的三维计算模型。Rohola Hasanpour等[11]在考虑TBM掘进速率的基础上，采用有限差分软件FLAC3D模拟双护盾TBM在高挤压区的施工过程。Kai Zhao等[12]采用三维模型分析了双护盾TBM穿越断层区时围岩失稳机理。程建龙等[13]采用FLAC3D软件建立完整的TBM模型，研究复合地层中双护盾TBM与围岩相互作用机制。

从现有的研究成果来看，国内外已经在采用数值方法分析双护盾TBM的施工掘进过程方面取得了长足的进步，由早期的二维计算分析[14-16]逐步发展到三维计算分析，且考虑因素越来越接近实际施工过程。下文将对双护盾TBM施工掘进过程的数值模拟关键技术进行深入探讨。

2 关键技术的探讨

根据目前双护盾TBM施工仿真分析的研究成果，同时结合现阶段深埋隧洞采用双护盾TBM遇到的各种问题，现将数值模拟深埋隧洞双护盾TBM施工过程的关键技术汇总如下：

(1)围岩本构模型的选择：对于深埋隧洞来说，其围岩力学特征与浅埋隧洞差别较大，表现出高度的非线性特征。对于硬岩，其峰后应力跌落特征非常明显，在数值模拟中常采用弹脆塑性本构模型，如文献[10]。对于软岩，其蠕变特征十分明显，但考虑到蠕变模拟的复杂性，现有研究成果大多仍采用常用的理想弹塑性Mohr-Coulomb模型，如文献[10-12]。在今后的研究中，应进一步考虑硬岩本构模型的选择和软岩蠕变效应对计算成果准确性的影响。

(2)外水压力的影响：在深埋隧洞施工过程中，高外水是经常遇到的一大难题，对围岩稳定及施工造成极大影响。但是在现有的研究成果中，大都未考虑外水压力的影响。一方面，考虑外水会大大增加计算工作量；另一方面，外水在深埋隧洞中是按裂隙水还是孔隙水考虑尚有争议。在后续的研究工作中，将进一步分析外水压力的影响。

(3)开挖速率(即:时效性)的影响：在实际施工中，双护盾TBM的掘进速率不是固定不变的，由于围岩类别、管片拼装和机械检修等因素的影响，掘进速率是不断变化的。在已有的研究中，只有文献[11]通过围岩应力释放率来粗略考虑掘进速率，但该方法只能应用于有限差分法，而无法应用于有限单元法。此外，该方法模拟的准确性严重依赖于应力释放率的取值，对计算人员要求较高。同时，文献[11]中围岩应力释放率采用的计算步数是固定的；而在实际计算过程中，不同开挖步达到同一应力释放率需要的计算步数是不同的，这点在后续工作中需要进一步改进。

(4)围岩与护盾的相互作用：在数值分析中，模拟围岩与护盾的相互作用是判断护盾卡机的关键。现有研究成果中，大都是在围岩与护盾之间建立接触单元[10-11,13]，再通过接触单元分析护盾与围岩之间的摩擦力，进而判断护盾是否发生卡机。此技术的关键在于真实考虑围岩与前护盾、尾盾之间的不均匀空隙。

(5)豆砾石和管片的支护模拟：对于双护盾TBM来说，其支护措施主要是管片和豆砾石回填灌浆；且在实际施工过程中，管片是紧跟尾盾拼装的，而豆砾石回填与灌浆距尾盾的距离较远，一般大于50 m。此外，豆砾石灌浆层是逐渐硬化的，其力学特性是与灌浆材料及时间密切相关的。但在数值模拟中，如果按实际施工过程考虑，则一方面会导致计算模型过大而无法计算，另一方面会导致计算过程十分复杂，不便于分析。鉴于此，文献[10]提出管片和豆砾石滞后尾盾一定距离再支护的方法，且对于软岩洞段，将豆砾石灌浆层分为软化和硬化两阶段来计算。该简化方法是根据已有深埋隧洞监测结果进行反演分析得到的，其普适性还需要进一步验证。

(6)围岩开挖的模拟：在数值分析中，模拟围岩开挖的方法主要有两种：一种是直接“杀死”(或“赋空”)围岩单元；另一种是先逐步折减要挖除围岩单元的变形模量，然后再“杀死”(或“赋空”)围岩单元。这两种方法从理论上都是合理可行的，但现有的研究成果[17]表明在有限差分法中，采用直接“杀死”单元法会导致围岩损伤塑性区过大，与钻爆法的实测成果更相符；而先折减围岩变形模量再“杀死”单元的方法更适用于TBM法开挖，因其对围岩的扰动较小。但现有文献[10-12]大都未考虑此因素，在后续研究中需考虑该因素对计算成果的影响。

基于已有研究成果，笔者基于FLAC3D软件及其自带FISH语言，初步实现了双护盾TBM施工过程的数值仿真程序，计算流程见图1。

图1 计算流程

与已有研究成果相比，本仿真分析程序改进之处在于：

(1)采用弱化围岩变形模量的方法实现围岩的逐步开挖，该过程是通过FISH语言实现的。此步骤的改进，可使围岩变形及塑性损伤区等结果更符合双护盾TBM施工掘进的实测成果;

(2)围岩应力释放率的考虑：在本分析程序中，不考虑围岩变形模量弱化过程中的应力释放率，仅考虑围岩开挖后的应力释放。在计算过程中，以体系开挖前初始的最大不平衡力作为基准，以体系剩余最大不平衡力与其商值作为评判标准。在本次计算分析时，暂定该标准为13%(即:围岩应力释放率为87%，与文献[11]相同)。该步骤也是通过软件自带的FISH语言实现的。

3 实例分析

基于现有研究成果，以某深埋引水隧洞双护盾TBM施工掘进过程为例，下文将详细介绍本仿真分析程序的计算过程，验证现有计算程序的可行性及可能存在的问题。

3.1 计算模型

由本程序的计算流程(见图1)可知，该仿真分析程序的基础是计算模型的建立。本文针对计算模型主要讨论两方面:计算模型的包含对象和计算模型的范围。

3.1.1 计算模型的包含对象

本仿真分析程序的目的主要有三点：一是围岩稳定的分析，二是TBM卡机的判别，三是管片衬砌等支护措施的可行性。根据这三个目的，且考虑计算的简化，计算模型应包含的对象有围岩、管片、豆砾石、前护盾、尾盾和刀盘等。在这些对象中，一般围岩、管片、豆砾石采用实体单元；而前护盾、尾盾和刀盘既可以采用实体单元，也可以采用壳单元，可根据计算分析的内容、规模及计算软件确定。本次计算时均采用了实体单元，图2为本次计算模型示意(未给出刀盘)。

图2 计算模型

3.1.2 计算模型的范围

根据前人的研究成果[18-19]，对于硬岩，在横断面方向，计算模型的尺寸应不小于隧洞直径的10倍；在隧洞轴线方向，计算模型的尺寸应不小于隧洞洞径的10倍。对于软岩，在横断面方向，计算模型的尺寸应不小于隧洞直径的15倍；在隧洞轴线方向，计算模型的尺寸应不小于隧洞洞径的20倍。在实际应用时，对于Ⅱ、Ⅲ类围岩，一般按硬岩考虑；对于Ⅴ类围岩，一般按软岩考虑；对于Ⅳ类围岩，可根据初始地应力及围岩变形情况合理选择。本次计算主要用于介绍仿真程序的计算流程及其可行性，故计算模型的范围是按硬岩考虑的。

3.2 计算过程

本仿真程序的主要计算过程如下：

(1)初始地应力平衡。在进行初始地应力平衡之前，首先赋空非围岩单元(刀盘、前护盾、尾盾、管片和豆砾石)，然后按照常规方法进行初始地应力平衡的计算。

(2)TBM开挖掘进。模拟双护盾TBM的开挖掘进，主要包括两方面:围岩的开挖和前护盾、尾盾的跟进。在本程序中，首先对要开挖的围岩进行变形模量的弱化。当围岩的变形模量弱化完成后，赋空相应围岩；同时前护盾和尾盾均前进一步，且相应的末端单元赋空，以保证前护盾和尾盾的长度不变。与此同时，前护盾、尾盾和相应的围岩之间建立接触单元，且删除三者之间无效的接触单元。最后再进行计算，且计算时通过围岩应力释放率控制计算何时终止。

(3)豆砾石和管片支护。当TBM开挖掘进计算终止后，激活滞后尾盾的豆砾石和管片。豆砾石和管片滞后尾盾的距离由围岩是硬岩还是软岩来确定[10-11]。此外，还需根据围岩是硬岩还是软岩确定豆砾石如何考虑。当围岩为硬岩时，豆砾石直接按硬化状态考虑；当围岩为软岩时，豆砾石先按软化状态考虑，随后按硬化状态考虑。根据笔者的计算分析，建议当围岩类别为Ⅱ、Ⅲ类时，豆砾石按硬化状态考虑；当围岩类别为Ⅳ、Ⅴ类时，豆砾石按软化和硬化两种状态考虑。

(4)计算结束的判别。当开挖和支护均未结束时，仍按步骤(2)、(3)计算；当开挖结束、支护未结束时，按步骤(3)计算；当开挖和支护均结束时，整个计算过程结束。

3.3 结果分析内容

根据双护盾TBM施工仿真分析程序的目的，需要分析施工掘进过程结束后围岩稳定性、TBM是否卡机和支护受力，三方面分析内容简述如下:

(1)围岩稳定分析。为了准确分析双护盾TBM施工完成后的围岩稳定性，一般选择模型轴向中间部分的围岩。一般来说，对围岩稳定分析，主要是分析开挖结束后围岩的压应力分布、水平向和竖向位移、围岩塑性区等内容。对于双护盾TBM施工方式来说，采用这几个指标有时并不能准确确定围岩稳定性，还需要结合卡机分析和管片受力分析来综合评价围岩稳定性。

(2)双护盾TBM卡机分析。分析双护盾TBM是否卡机，主要是计算前护盾和尾盾与围岩的接触摩擦力是否大于机械提供的最大推力。在计算接触摩擦力之前，首先需判断围岩与前护盾或尾盾是否接触，这可以从前护盾或尾盾处的围岩变形来初步判断。当初步判断围岩与前护盾或尾盾接触后，再查看接触面的法向应力，以确定可能卡机部位。确定可能卡机部位后，再计算接触面上的法向力(即:围岩传递给护盾的压力)，可通过FISH语言对接触面节点法向应力与其控制面积的乘积求和获得。最后，将围岩与护盾之间的摩擦系数乘以接触面的法向力得到围岩与护盾之间的摩擦力。若摩擦力小于TBM的最大推力，则双护盾TBM不会卡机；否则，双护盾TBM将会卡机。

(3)管片受力分析。当分析完围岩稳定性和双护盾TBM卡机后，还需要进一步分析管片受力。这是因为一方面管片自身承载力是有限的，另一方面从管片受力可间接评价围岩稳定性。当双护盾TBM卡机时，管片受力可不做分析；但当双护盾TBM未卡机时，必须要分析管片受力，以防管片自身承载力不足而导致支护措施不够。

通过上文的介绍，以埋深1 000 m的Ⅳ类围岩为例，对其施工掘进的结果进行分析。图3～4分别为围岩竖向位移云图及其塑性区分布(图中位移单位以m计；塑性区中深色区表示张拉-剪切屈服，浅色区表示剪切屈服)。从位移来看，顶拱沉降为-123.7 mm，底拱向上鼓起129.2 mm。从围岩塑性区来看，洞周0.6 m范围内为拉-剪破坏，向内2.8 m范围内基本为剪切破坏。仅从围岩位移和塑性区来看，围岩可初步认为是稳定的。

图5～8分别为前护盾和尾盾位置的围岩变形和接触面法向应力云图。

图3 围岩竖向位移云图 图4 围岩塑性区分布

图7 尾盾处围岩水平向变形 图8 尾盾接触面法向应力

由图5～8可知，前护盾处的围岩最大变形量为55.4 mm，大于前护盾与围岩之间的空隙(50 mm)，这表明围岩与前护盾发生接触。接触面法向应力最大值为37.61 MPa，较大法向应力主要位于后部一环。尾盾处的围岩最大变形量为100.4 mm，大于围岩与尾盾之间的空隙(100 mm)，这表明围岩与尾盾也发生了接触。接触面法向应力最大值为21.64 MPa，较大法向应力主要分布在两腰。

由前文分析可知，前护盾、尾盾与围岩之间均发生接触，且接触摩擦力分别为373.2 MPa和281.7 MPa，远大于双护盾TBM的最大主推力(61.6 MN)，故TBM会发生卡机现象。

由于双护盾TBM会发生卡机现象，故可不必再分析管片受力。即使分析管片受力，也可知管片受力会超过其抗压强度，如图9所示。

图9 管片第三主应力云图

4 结 论

基于已有研究成果，本文采用FLAC3D软件及其自带FISH语言重新编写了双护盾TBM施工过程的数值仿真分析程序，并以某深埋引水隧洞为例，介绍了该仿真分析程序的计算流程，并验证了其可行性。

与原有研究成果相比，本仿真程序主要在围岩开挖模拟和围岩应力释放率两方面进行了改进。总体来看，现有程序还不完善，还无法考虑外水、超前处理措施及时间效应等因素的影响，在后续工作中将进一步完善该仿真程序，以使其更具有普遍适用性及工程应用价值。

[1] 编纂委员会.岩石隧道掘进机(TBM)施工及工程实例[M].北京:中国铁道出版社，2004.

[2] 山西省万家寨引黄工程管理局.双护盾TBM的应用与研究[M].北京:中国水利水电出版社,2011.

[3] 杜士斌.开敞式TBM的应用[M].北京:中国水利水电出版社，2011.

[4] 吴世勇，周济芳，陈炳瑞，等.锦屏二级水电站引水隧洞TBM开挖方案对岩爆风险影响研究[J].岩石力学与工程学报,2015,34(4)：1-7.

[5] Q.M. GONG, L.J. YIN, S.Y. WU, et al. Rock burst and slabbing failure and its influence on TBM excavation at headrace tunnels in Jinping Ⅱ hydropower station[J]. Engineering Geology, 2012, 124(4): 98-108.

[6] 刘泉声，黄兴，时凯，等. 超千米深部全断面岩石掘进机卡机机理[J]. 煤炭学报, 2013, 38(1)：78-84.

[7] 黄兴. 深部软弱地层TBM掘进围岩挤压大变形与卡机致灾机制[D]. 武汉：中国科学院武汉岩土力学研究所，2014.

[8] 刘泉声，黄兴，时凯，等. 深部挤压性地层TBM掘进卡机孕育致灾机理[J]. 煤炭学报，2014,39(增1)：75-82.

[9] 杜立杰，齐志冲，韩小亮，等. 基于现场数据的TBM可掘进性和掘进性能预测方法[J]. 煤炭学报，2015,40(6)：1248-1289.

[10]Kai Zhao, Michele Janutolo, Giovanni Barla. A completely 3D model for the simulation of mechanized tunnel excavation[J]. Rock Mech Rock Eng, 2012,45(4): 475-497.

[11]Rohola Hasanpour, Jamal Rostami, Bahtiyar Yunver. 3D finite difference model for simulation of double shield TBM tunneling in squeezing grounds[J]. Tunnelling and Underground Space Technology, 2014, 40(2): 109-126.

[12]Kai Zhao, Michele Janutolo, Giovanni Barla, et al. 3D simulation of TBM excavation in brittle rock associated with fault zones: the Brenner Exploratory Tunnel case[J]. Engineering Geology, 2014, 181: 93-111.

[13]程建龙，杨圣奇，杜立坤，等. 复合地层中双护盾TBM与围岩相互作用机制三维数值模拟研究[J]. 岩石力学与工程学报，2016,35(3)：511-523.

[14]Ramoni M, Anagnostou G. On the feasibility of TBM drives in squeezing rock conditions[J]. Tunnelling and Underground Space Technology, 2006, 21(3-4): 262-268.

[15]Ramoni M, Anagnostou G. Tunnel boring machines under squeezing conditions[J]. Tunnelling and Underground Space Technology, 2010, 25(2): 139-157.

[16]Ramoni M, Anagnostou G. The interaction between shield, ground and tunnel support in TBM tunneling through squeezing conditions[J]. Rock Mech Rock Eng, 2011, 44(1): 37-61.

[17]Itasca Consulting Group. FLAC 3D, fast lagrangian analysis of continua in 3dimensions, version 5.0[M]. Itasca Consulting Group, Minneapolis, Minnesota, USA, 2012.

[18]Eberhardt E. Numerical modelling of three-dimension stress rotation ahead of an advancing tunnel face[J]. Int J Rock Mech Min Sci, 2001, 38(4):499-518.

[19]Graziani A, Ribacchi R, Capata A. 3D modelling of TBM excavation in squeezing rock masses[C]. Brenner Basistunnel und Zulaufstrecken, Internationales Symposium BBT 2007, Innsbruck, Austria. Innsbruck University Press, Innsbruck, 2007.

2017-01-12

孙博(1984-)，男，河北三河人，博士，高级工程师，从事水工结构、隧洞及边坡的设计与研究工作。

TV554

：B

：1003-9805(2017)03-0065-05