采空区长期变形对桑掌隧道影响的数值模拟研究

刘晋华 朱义嘉 郑 刚 刘伟博 夏江都

（1.山西省公路局阳泉分局，山西 阳泉 045000；2.中国建设基础设施有限公司，北京 100029；3.山西中建基础设施建设运营有限公司，山西 阳泉 045000；4.中国地质大学（北京）工程技术学院，北京 100083）

0 引言

我国拥有丰富的煤炭资源，同时也是煤炭开采和消费大国。经过几十年的大规模开采，山西等采煤大省形成了大量采空区［1］。采空区变形经历一个相对较长的过程，对于桥梁、隧道等对变形要求严格的结构类公路工程，采空区的长期变形影响会对其造成较大的危害［2］。当建设场地位于采空区影响范围内时，采空区长期变形对隧道沿线稳定性影响的问题值得深入研究［3］。

众多学者针对采空区影响下隧道建设进行了深入研究。时均伟［4］总结分析了华蓥山隧道采空区段及国内外相关工程，并进行数值模拟得到隧道采空区段稳定性影响的情况，为相关工程提供宝贵的施工经验。黄明等［5］利用数值模拟软件对沪昆高速公路灯草塘隧道的稳定性进行了分析，得到靠近采空区一侧与另一侧的应力、位移有明显差异。姚宣成等［6］对隧道下伏采空区稳定性进行了数值模拟，指出采空区加固前后，隧道出现不同的现象，加固采空区后，隧道稳定性得到大幅提升。符亚鹏等［7］研究了不同倾角采空区对邻近公路隧道的稳定性的影响，提出了下伏采空区会增大隧道沉降值、塑性区范围、应力值，并指出倾角越小对隧道的影响越大。蒋永亨［8］以D-InSAR 监测为验证方式，通过正交试验得到采空区上覆岩层的力学参数，最后通过数值模拟得到隧道采空区段暂时处于稳定状态。

桑掌隧道穿越阳泉二矿采空区，煤层停采时间与隧道建设时间相距较近，本文以桑掌隧道高速公路采空区段工程为依托，采用数值模拟方法，针对受下伏采空区长期变形影响下的隧道围岩、支护结构在隧道建设不同时期的变形破坏情况展开研究，分析经历长期变形作用后隧道围岩、支护结构在不同时期的应力响应特征，研究成果对隧道稳定性评价具有指导意义。

1 研究区工程概况

桑掌隧道位于山西省东部阳泉二矿范围内，阳泉二矿井田面积约60.06 km2，开采深度463.3～713.5 m，开采煤层分别为山西组（C3s）地层中的3#煤层和太原组（C3t）地层中的8#、15#煤层，平均厚度分别为1.96 m、2.26 m、6.42 m，煤层发育近似水平。新建桑掌隧道为左右线分离式洞体，全长925 m，平均埋深80 m（最大埋深113 m），修建时间为2019—2020 年。隧道修建时下伏煤层经过开采已形成采空区。隧道从二叠系上统上石河子组下段地层穿过，地层岩性为砂质页岩和砂岩。由于桑掌隧道修建距15#煤层开采的时间最近，并且15#煤层埋深、采厚均大于3#、8#煤层，故桑掌隧道主要受15#煤层开采后的长期变形影响。15#煤层工作面停采年份及隧道穿越采空区位置如图1所示。

图1 隧道与15#煤层采空区相对位置示意

2 模拟方法

拟建桑掌隧道的所在位置受到下伏多煤层采空区长期变形影响，若将隧道与采空区作为一个整体进行研究，需要建立范围较大的数值模型。在此大模型中，隧道所处区域附近地层网格尺寸远远超出了隧道支护结构的尺寸大小，不利于隧道结构建模及后续计算。基于上述问题，本研究借鉴位移反分析法，建立隧道及围岩的小尺度模型。在模拟开挖后，根据研究区变形位移场计算小尺度模型中各节点在不同时间的位移速率，最后以不同时间段各节点位移速率为边界条件，耦合计算得到隧道及围岩的变形情况。

具体模拟方法如下：首先，整理研究区域内变形位移场；其次，建立隧道及围岩模型并进行开挖获得初始应力场；再次，根据编写的三维克里金代码插值还原隧道数值模型残余变形位移量，得到隧道所处区域的采空区围岩的位移变形量；最后，通过FLAC3D还原精细化隧道计算模型残余变形位移场，研究采空区长期变形对桑掌隧道的影响。

2.1 桑掌隧道场地长期变形位移场的建立

为探究采空区长期变形对桑掌隧道的影响，在进行本次研究前先对桑掌隧道周围采空区场地进行了长期变形计算，计算结果如图2 至图4 所示。由此可知，15#煤层的80804、80806 工作面是导致场地残余变形的主要因素。桑掌隧道处于沉陷中心，隧道中心区域沉降大，变化量小，隧道两侧沉降小，变化量大；场地的X方向水平变形以桑掌隧道为中心向桑掌隧道所处区域移动；受到场地Y方向水平位移影响的主要为北侧隧道洞口。

图2 2019—2109年地表X方向水平位移等值线

图3 2019—2109年地表Y方向水平位移等值线

图4 2019—2109年地表残余沉降等值线

2.2 隧道计算模型的建立及开挖模拟

相较于采空区影响场地，桑掌隧道建设场地范围小，需要建立节点间距与网格尺寸更小的精细化计算模型以便后续计算。根据前人的研究，由洞室开挖引起的围岩变形一般不超过三倍的洞室跨度［9］，桑掌隧道单洞跨度12 m、高度10 m。为保证隧道及围岩高精度计算模型边界不受隧道开挖扰动，模型左右侧边界距离隧道边界取45 m。隧道全长986 m，前后边界与隧道洞口的距离取20 m，模型下边界距离隧道底板的距离取75 m，模型上边界为地表。隧道及围岩高精度模型建模范围如图5 所示，最终模型大小为长1 006 m、宽140 m、高85～193 m。

图5 隧道及围岩高精度模型建模范围

为了更加准确地模拟不同的隧道支护结构在长期变形作用下的变形、受力、破坏情况，本研究采用不同的单元来模拟不同的构件。根据深部钻孔岩芯力学试验数据，并结合围岩等级对岩体进行取参，岩体采用四面体单元进行剖分，使用摩尔-库伦本构模型计算；初期支护和二次衬砌的计算参数分别按照C25 混凝土和C30 混凝土进行取值，初期支护和二次衬砌采用实体单元进行模拟，使用弹性本构模型计算。模型中各部分的计算参数见表1。

表1 数值模拟物理力学参数

模型建立后，采用置空的方式模拟隧道开挖。沿隧道开挖方向按10 m 一段对被开挖岩体网格进行置空，在开挖形成的临空面上设置Cable 单元以模拟锚杆支护过程，同时以赋予置空处单元弹性模型及混凝土参数的方式模拟喷射混凝土的过程，并对需要进行二次衬砌支护的位置采用同样的赋参方式进行模拟。按照上述流程依次完成开挖—初期支护—二次衬砌—开挖的循环，直至计算完成。

隧道开挖后围岩横截面的应力分布如图6、图7所示，隧道开挖后围岩在隧道周边约1倍洞径的范围内应力等值线有弯曲现象，说明该区域为隧道开挖的主要扰动区域。从总体来看，洞室及周边围岩处于相对稳定和可控的状态，可以进行后续的模拟计算。

图6 隧道围岩横断面（1 433 m）最大主应力

图7 隧道围岩横断面（1 433 m）最小主应力

2.3 隧道围岩长期变形位移场的施加

将隧道开挖后产生的围岩应力作为初始应力，在此基础上使用场量插值方法对隧道围岩网格节点进行位移插值。插值后可得到隧道区域围岩各网格节点的残余变形位移量，将残余变形量按照边界条件加载到隧道围岩的三维数值模型上，如图8所示。

图8 节点位移边界示意

在获得隧道围岩各个节点在运营期的残余变形量后，为计算围岩残余变形在较长时期对隧道及支护结构的影响，需要采用对围岩四面体网格节点按不同时期添加变形速率的方式进行模拟。如图9所示，对围岩中任意四面体网格节点，按照长期变形各时间段对网格节点的位移，结合计算时间添加变形速率vx1、vy1、vz1。以此方法可建立隧道围岩模型在采空区长期变形影响下不同时期的速率/位移场。

图9 隧道单元节点速度—位移加载示意

3 计算结果

模型建立位移场后，计算得出隧道围岩在采空区长期变形作用下塑性区的变化情况，如图10 所示。由图10 可知，桑掌隧道从2019 年隧道开挖完成到2109年的90年间，围岩塑性区将逐渐增大，隧道围岩的剪切破坏、拉伸破坏区域逐渐向隧道中心扩大。前20 年隧道围岩主要以拉伸破坏为主；后20年拉伸破坏范围不断扩大，同时开始出现剪切破坏；之后的50 年里，隧道拉伸破坏、剪切破坏范围不断扩大。隧道围岩塑性区最先出现在隧道两端洞口处，首先发生沿隧道方向的拉伸破坏，这是因为隧道中间段处于沉降中心，两侧洞口端相较于中心段沉降较小，但此处倾斜变形较大，引起围岩的拉伸、剪切破坏；随着时间的推移，隧道洞口端也会发生剪切破坏，同时隧道围岩塑性区逐渐增大，围岩破坏区域逐渐由两端向隧道中心扩大；但是桑掌隧道围岩塑性区的扩张范围较小，主要集中在隧道两端洞口处，隧道洞口处的围岩处于欠稳定状态，其余部分围岩基本没有发生破坏，处于稳定状态。需要着重检查隧道洞口段洞室后期的维修工作，确保隧道运营期间的安全。

图10 长期变形作用下（2109 年）隧道围岩塑性区

残余变形作用下隧道衬砌支护应力分布情况如图11 所示。由图11 可知，桑掌隧道衬砌应力由对称分布沿着隧道走向向右（东侧）发生偏转现象。隧道斜穿采空区下沉盆地，导致隧道与沉降等值线斜交，且隧道东侧更靠近下沉中心，隧道东侧竖向位移最大，隧道横截面的下沉中心在底板东侧，因此隧道衬砌应力在采空区残余变形作用下由原来的对称分布逐渐向右侧偏转，在右侧墙、左侧墙墙角出现拉应力集中、剪切应力集中现象；而在左侧墙、右侧墙墙角出现压应力集中现象。隧道衬砌yz方向剪应力在左右侧墙均出现集中现象，这是因为在采空区残余变形作用下，隧道沿走向的残余沉降值不断变化，发生不均匀沉降，导致隧道侧墙横截面受到来自围岩z方向剪应力发生上下错动。

图11 残余变形（2109年）作用下隧道衬砌支护应力云图

在桑掌隧道沿走向线上应力集中区的右侧墙、墙脚的单元进行应力提取并绘制隧道衬砌支护残余应力曲线，如图12所示。由图12可知，隧道各方向应力沿隧道走向会出现两个峰值，即距隧道洞口距离30～120 m、750～850 m 区域处，其中30～120 m段隧道的zz方向拉应力最大值在隧道建设20 年后（即2040 年）衬砌C30 混凝土抗拉强度达到2.01 MPa，并于2119 年达到2.8 MPa。原因是桑掌隧道在建设后20 年至40 年间累积破坏程度较大，隧道围岩率先发生拉伸破坏，主要集中在隧道两侧洞口段，在后期运营检查中应重点关注。

图12 隧道衬砌支护六项残余应力曲线

桑掌隧道衬砌结构位移情况如图13 所示。由图13（a）可知，桑掌隧道在建设后100 年内，隧道沿线沉降速率逐年减小，前20 年的下沉量达到了总下沉量的80%。隧道沿线残余沉降曲线呈U 形，隧道南侧洞口处最终沉降量约为80 mm；在隧道中间区域，即距离南侧洞口400～600 m 范围内，隧道残余沉降维持在230～260 mm 之间；至隧道北侧洞口处残余沉降减小20 mm 左右。这是由于隧道中心位于15#煤层80804 工作面中心的正上方区域，而残余沉降中心与80804 工作面中心点几乎重合，隧道斜穿采空区下沉盆地中心区域，导致隧道地基竖向位移呈现中间大、两侧洞口小的状况。

图13 桑掌隧道衬砌结构位移

由图13（b）可知，在采空区残余变形作用下，桑掌隧道大部分区域由东向西侧移动。隧道沿线x方向的水平移动曲线呈V 形。南北两侧隧道洞口附近水平移动为正值，向东侧移动；在距离南侧洞口100～850 m 处水平移动为负值，向西侧移动，在距离南侧洞口约620 m 处水平移动达到最大值，约为48 mm。

由图13（c）可知，在采空区残余变形的作用下，隧道y方向水平位移在距隧道南侧洞口0～400 m 段水平移动为正，在距离洞口320 m 达到峰值3 mm；400～980 m 段（北侧）为负，在距离南侧洞口750 m处达到峰值约-14 mm；隧道南北两端的围岩朝隧道中部移动，隧道围岩在y方向受到挤压作用。

4 结论

通过采用数值模拟方法，研究阳泉二矿采空区长期变形对桑掌隧道围岩及衬砌的影响，得出以下结论。

①采空区长期变形将对隧道围岩的变形及稳定性产生影响，围岩的竖向位移表现为整体下沉。隧道围岩塑性区逐渐增大，隧道围岩的破坏主要集中在隧道两端洞口处，隧道洞口处的围岩处于欠稳定状态，其余部分围岩基本没有发生破坏，处于稳定状态。

②隧道衬砌最大应力沿着隧道走向向右（东侧）产生偏移、旋转现象。在桑掌隧道南北侧洞口段区域应力最大并出现峰值，其中距南侧洞口30～120 m 段隧道的zz方向拉应力最大值2.8 MPa，超过了C30 混凝土抗拉强度。隧道围岩率先发生拉伸破坏，且主要集中在隧道两侧洞口段，在后期运营检查中应重点关注。

③隧道位移量不断增加，每一年位移增量逐渐减小，前20 年的位移量达到了总位移量的80%。隧道残余沉降曲线呈U 形曲线，最大值为257.79 mm。隧道沿线x方向的水平移动曲线呈V形，最大值约为48 mm。隧道沿线y方向水平位移由南北两侧向距隧道南侧洞口400 m 区域移动，在距离洞口750 m达到最大值14 mm。