基于GDEM的隐伏岩溶隧道隔水岩体水压致裂安全厚度及破裂演化规律分析

李 涛, 张 丽, 蒋 庆, 冯 春, 赵 然

(1. 山东高速济莱城际公路有限公司， 山东 济南 250014； 2. 中国科学院力学研究所， 北京 100190)

0 引言

在富水岩溶地区，隧道与溶洞间的岩层由于开挖引起的扰动而极易失稳，引发突水、涌泥灾害，是隧道安全施工的重大问题，也是目前工程建设研究的重点之一。目前研究岩溶隧道隔水岩层稳定性的方法主要有模型试验、理论研究和数值模拟等。

模型试验需根据原型确定试验模型的几何相似比，然后根据力学方程和边界条件确定应力相似比，并采用相似材料模拟一定工况下的突水现象。魏星等[1]以北岗隧道DK462+527～+537段为原型，采用1∶100几何相似比模型，模拟了围岩出现位移剧增、突泥涌水和破坏的现象，但材料的力学指标未能达到相似性要求。赵瑜等[2]建立了应力相似比为40、几何相似比为55的某在建隧道模型，并与FLAC3D数值模拟结果进行对比，发现在各种趋势上有较好的一致性，且物理模型的最终破坏及影响范围与数值模拟有较好的一致性，但局限于材料特性和测量上的误差，未能体现数值模拟的全部现象。李浪等[3]研制出一套深长隧道突水地质灾害三维模型试验系统，该系统不仅能够模拟地应力场、水压等初始环境，还可精确模拟隧洞开挖，但该模型试验是在隧道模型简化条件下进行的。Yang等[4]研发了高地应力高渗隧道突水模型试验系统,并以重庆西马隧道为研究对象，建立了三维流固耦合模型，通过监测应力、位移、渗流压力等的变化趋势，较为真实地模拟了突水发生的过程，但采用的是相似比材料和简化模型。鉴于地层结构的复杂性和研究条件限制，模型试验只能一定程度地模拟突水现象，对工程的指导性有限。大多数理论研究是通过对现场地质模型进行一系列简化而得到相应的理论模型，对于简化后的理论模型可以用结构力学、弹性力学和断裂力学中的相关公式加以计算，从而得到临界安全厚度[5]。李利平等[6]以断裂力学、弹塑性力学理论为基础，推导出隧道最小岩石保护厚度的半解析解表达式，并认为岩溶隧道裂隙突水具有明显的时空效应特征。郭佳奇等[7-8]针对侧壁溶腔，以塑性区和高渗透带贯通与否作为判定中间岩柱最小安全厚度的标准，建立最小防突厚度计算公式；基于弹性厚板理论，推导了边界条件分别为固支和简支2种模式的岩溶隧道掌子面岩墙安全防突厚度及临界水压计算公式。Yang等[9]基于Hoek-Brown非线性破坏准则，利用上界定理和变分原理推导获得了隔水岩柱的安全厚度，并研究了各参数的影响，绘制出了突水破坏区域。地质工程的复杂性使各种理论研究和模型试验缺乏有力的支持，相比之下，数值模拟方法具有较广泛的适用性,成为解决岩土工程问题的有效工具。随着计算机技术的迅猛发展, 各种数值计算方法越来越普遍地被应用到围岩稳定性分析中[10-11]。

目前，隧道围岩采用的数值计算方法主要为有限元法、有限差分法和离散元法等，ABAQUS、MIDASGTS等为有限元软件；FLAC3D为有限差分软件；UDEC、3DEC等为块体离散元软件。佘健等[12]针对笔架山隧道建立了三维有限元模型，采用弹塑性模型进行计算，通过监测地表、拱顶、拱腰等的位移变化，分析开挖的影响。谢海文等[13]利用二维离散元数值模拟软件UDEC，分析了贵州省德江隧道围岩中裂隙的响应特征、流固耦合效应和渗流场分布特征。黄明利等[14]利用岩石破裂过程分析程序(RFPA), 建立了二维平面应变模型, 对隧道施工诱发隐伏溶洞破裂突水过程中的应力场、位移场和声发射等特征进行了系统研究。徐长金等[15]基于MIDAS/GTS有限元分析软件，分析了溶洞半径以及溶洞与隧道底板距离的变化对围岩应力分布和塑性区范围的影响。雷霆等[16]利用FLAC3D有限差分软件进行三维隧道开挖模拟，分析了顶部溶洞与隧道掌子面附近围岩塑性区的发展趋势，以及不同溶洞水压下隧道的安全厚度。李红卫[17]采用数值分析软件ABAQUS对贵阳轨道交通建设中的隧道底部溶洞围岩稳定性进行三维数值模拟，分析了距隧道断面不同间距下溶洞对隧道底部、顶部和地表位移的影响。Shan等[18]利用FLAC3D建立了三维隧道模型，对3种不同位置的溶洞开展了数值模拟，并采用经验公式分析了各因素对安全厚度的影响程度，认为数值方法更具合理性。岩体的破裂包含小变形、损伤演化、裂纹萌生、裂纹扩展以及大位移、大转动等阶段，是一个连续到非连续的破坏过程。因此，隔水岩体的破坏是开挖扰动引起的裂隙扩展、贯通直至破裂的渐进过程[6]。有限元法及有限差分法能够较好地模拟材料在连续状态下的特性，但不能模拟材料从连续到非连续的过程及在非连续状态下的运动特性；块体离散元在模拟非连续体的运动特性方面具有一定的优势，但较难模拟材料的连续变形过程。目前大多数研究以溶洞的空间二维分布形态为主，三维模型受软件本身计算能力的影响，网格较大，精度较差，计算耗时多。因此，大多数值模拟软件无法真实模拟三维隧道的开挖突水破坏过程。

从上述分析可见，模型试验、理论研究对于描述突水这种复杂地质体的力学行为具有一定局限性，而传统数值模拟方法在描述岩体渐进破坏过程和计算精度方面均有不足。本文选取GDEM系列中的DAS软件，建立百万单元三维溶洞隧道模型，分析在高地应力、高溶洞水压下，开挖扰动引起的岩体破坏演化规律；为了分析隔水岩体的渐进破坏过程，引入无量纲指标破裂度来定量化描述岩体的破裂状态；通过监测掌子面位移变化特征，获得隔水岩体的安全厚度，并分析在这一过程中隔水岩柱的破裂发展趋势与位移之间的对应关系。

1 基于GPU加速的CDEM算法

GDEM应力分析系统(GDEM-DAS)是基于连续-非连续介质力学的离散元方法CDEM(continuum-based discrete element method)[19-21]的高性能有限元-离散元计算软件。该软件以CDEM理论为基础，利用GPU(显卡)进行加速计算。CDEM将有限元与离散元进行耦合，在块体内部进行有限元计算，在块体边界进行离散元计算，通过块体内部及块体边界的断裂，不仅可以模拟材料在连续状态下及非连续状态下的变形、运动特性，更可以实现材料由连续体到非连续体的渐进破坏过程，目前已广泛应用于滑坡、爆破、地下工程等行业[22-24]。

1.1 CDEM算法的基本概念

CDEM中的数值模型由块体及界面2部分构成，如图1所示。块体由1个或多个有限元单元组成，用于表征材料的弹性、塑性、损伤等连续特征；2个块体间的公共边界即为界面，用于表征材料的断裂、滑移、碰撞等非连续特征。界面包含真实界面及虚拟界面2个概念，真实界面用于表征材料的交界面、断层、节理等真实的不连续面，其强度参数与真实界面的参数一致；虚拟界面主要有2个作用，一是连接2个块体，用于传递力学信息，二是为显式裂纹的扩展提供潜在的通道(即裂纹可沿着任意一个虚拟界面进行扩展)。

(a) 数值模型 (b) 块体 (c) 界面

1.2 离散特征

采用CDEM方法描述连续-非连续问题时，连续部分采用有限元离散，非连续部分添加界面接触表征其力学特征，相邻接触面通过空间拓扑信息进行查找，在公共接触面两侧块体的相应节点间建立弹簧(包括法向弹簧和切向弹簧)，接触面之间的相互作用力由弹簧力和弹簧的特征面积来表征。假设计算区域被结构面T1、T2和T3切割后形成3个区域，各个区域内需要进行有限元离散，接触边界也需要进行离散化并建立界面接触来表征结构面，如图2所示。

(a) 块体单元 (b) 接触弹簧

1.3 破裂判断准则

CDEM采用基于时程的动态松弛技术进行显式迭代计算，因此可获知每时步的单元应力和节点弹簧力，并根据破裂准则进行强度判断，将发生破裂的节点弹簧力置为0，并记录破裂弹簧的特征面积。本文拟在块体部分采用线弹性模型，界面接触弹簧采用张拉-压剪复合准则判断其破坏状态。张拉模式的弹脆性破坏采用最大张力准则：

(1)

式中：σn为接触法向应力；kn为法向弹簧刚度；u为弹簧的法向位移；σt为接触抗拉强度。

当界面发生张拉破坏时，界面上的剪切应力也相应地置为0。剪切模式的弹脆性破裂采用摩尔-库仑准则 ：

(2)

式中：τ为接触切向应力；kτ为切向弹簧刚度；ν为弹簧的切向位移；c为接触黏聚力；φ为内摩擦角。

当界面接触发生破裂之后，接触所连接的两侧块体发生张开和滑动。

2 计算模型及材料参数

2.1 模型尺寸

根据赵明阶等[25]的研究成果和弹塑性理论，隧道周边影响范围为3～5倍的洞径，因此，建立三维隧道计算模型为100 m×100 m×100 m，三心圆隧道位于模型中间，开挖跨度为19 m，高度为11 m，衬砌厚度为55 cm，如图3(a)所示。溶洞简化为椭球模型，长轴为30 m，短轴为14 m，小于隧道跨度。考虑突水最不利位置，将溶洞置于隧道掘进的正前方，溶腔垂直隧道轴线方向(z向)的尺寸约为隧道高度的2倍，如图3(b)所示。采用前处理功能较强的ANSYS建立几何模型后，共划分105万个四面体单元、420万个节点，然后通过数据转换程序将节点单元导入到GDEM-DAS中，再进一步进行岩层属性设定、开挖衬砌模拟和计算分析等。

(a) 三维溶洞隧道网格

2.2 边界条件

计算模型设定为埋深130 m的深埋岩溶隧道，符合平面应变模型，因此，计算模型左右侧边界、前后边界和底部均采用法向位移约束；模型上边界(z=100 m)施加上端岩体自重；根据京沪高速某公路段地应力资料，围岩侧压力系数取1.0，转换为线性分布荷载分别施加于计算模型z=0 m、z=100 m、x=0 m和x=100 m 4个平面上。

鉴于溶腔轴的尺寸相对较小，可忽略溶腔内部上下水压差别，水压简化为施加于椭球内表面的均布压力。李利平等[6]研究发现瞬时突水突泥型溶洞水压一般为1～3 MPa，本文拟分别计算水压1、2、3 MPa时的隔水岩层破坏状态。

2.3 材料参数

因地质情况复杂，各地岩石参数差异较大，本文围岩参数选自京沪高速某公路段现场测试结果。该路段采用初期支护和二次衬砌支护，初期支护包含钢筋网、钢拱架、锚杆和混凝土等材料，二次衬砌包含C30模筑混凝土和钢筋，各材料属性均不同，故采用等效参数进行计算。围岩和支护参数取值见表1。

表1 围岩与支护参数

2.4计算过程

计算过程为先采用弹性模型计算获得岩体初始应力状态，再采用强度准则模拟出围岩的初始破裂状态，然后进行开挖过程计算： 采用全断面开挖模型，开挖进尺2 m，在开挖下一进尺的同时完成上一进尺的衬砌支护，每次开挖计算均采用破坏准则迭代至40 000步。

选取掌子面为监测面，分别在每个监测面的拱顶、拱底、拱腰和中心位置设置监测点，如图4所示。通过监测开挖过程中监测点的位移变化规律，分析开挖扰动对隔水岩层的影响。

(a) 监测面 (b) 监测点

2.5 隔水岩体安全厚度

从数值模拟的角度来讲，位移能够直观地表征岩体的渐进破坏过程。在同一计算工况下，随着掌子面向溶腔方向的推进，开挖引起隔水岩体出现裂缝，随后逐渐扩展，此时掌子面的位移缓慢增加；当完全贯通时，掌子面的位移会迅速增加，存在突变现象，说明此时隔水岩体的厚度已不具备足够的安全储备。因此，可将位移值突变发生的前一开挖步视为安全厚度，在此时采取工程措施进行防治。

2.6 隔水岩体的破裂度

2.6.1 破裂度的定义

破裂是表征岩体灾变过程的重要现象，隔水岩体的破裂程度是影响突水发生的重要因素。现有的测量手段可直接测得暴露的裂缝长度，但无法由此推测出岩体内部的损伤程度以及渐进破坏过程。

为了建立隔水岩体渐进破裂程度与开挖进程之间的关系，本文借助于连续-非连续数值计算方法，模型中块体之间的界面即为潜在的破坏面，利用强度准则判断界面状态，通过弹簧断裂模拟裂缝扩展，进而通过统计弹簧力为0的特征面积获得当前状态下岩体的破裂面积。将当前状态下的破裂面积与界面总面积之比定义为破裂度，即破裂度

D=S/St×100%。

(3)

式中：S为研究区域内当前断裂弹簧特征总面积；St为研究区域内界面总面积。

2.6.2 单轴压缩破裂度计算和网格依赖性验证

下面以单轴压缩为例说明破裂度的计算方法，并进行网格依赖性验证。计算模型为直径50 cm、长100 cm的圆柱，上下端同时施加相同的加载速率0.5×10-9m/s，计算模型参数见表2，界面弹簧参数通过材料参数计算获得。

表2 计算模型参数

计算模型采用张拉-压剪复合准则进行计算，由于破裂发生在界面，通常选用随机性较好的四面体网格进行计算。不同网格数量的计算模型如图5所示。为了验证计算结果的可靠性，对四面体网格进行加密，计算模型网格数量见表3。

(a) 模型1

表3 计算模型网格、接触面和弹簧组数量

研究区域为整个圆柱，单元之间均采用界面接触模型进行计算，通过弹簧连接，接触面数量和弹簧组数见表3。采用1.3节的破裂准则进行强度计算，通过统计各监测时步弹簧力为0的特征面获得当前时步的破裂总面积，弹簧力不为0的特征面积为当前时步的界面总面积，二者之比即为当前时步的破裂度。破裂度随时步变化曲线如图6所示。计算结果显示，在同一计算工况下，不同密度网格的破裂度随加载时步的变化趋势基本一致，都存在突变点，在突变点后破裂度逐渐趋于稳定；相邻网格数量的破裂度误差小于4.3%。由此可见，破裂度对网格的依赖性较小。

图6 破裂度随时步变化曲线

3 隔水岩层破坏特征分析

3.1 隔水岩体监测点位移变化分析

以2 MPa溶洞水压为例，计算模型在重力和侧压力作用下，弹性稳定时的y向位移云图如图7所示，开挖过程中的z向位移云图如图8—10所示，开挖进程中各监测点位移见表4，监测点处位移随开挖进程的变化曲线如图11所示。

计算结果说明：

1)随着掌子面的掘进，轴向位移扰动区呈漏斗形，掌子面中心处位移值最大，破坏模式比较符合弹性圆板模型，以张拉破坏和剪切破坏为主。

2)掌子面中心点处的位移值最大，说明中心处最容易发生破坏，符合隔水岩体一侧受溶洞水压、另一侧临空的受力状态。

3)侧边位移较小且变化不大，说明溶洞水压对拱腰影响很小。

4)拱顶、拱底和中心处监测点位移均在距溶洞7 m处存在突变点。突变点左边位移值较低，最大值在中心处(约0.003 m)，为溶洞水压弱影响区，隔水岩体处于稳定阶段；突变点右边，位移先缓慢增加后迅速增加，为溶洞水压强影响区，隔水岩体失稳。

5)掌子面距溶洞5 m处，掌子面中心处位移突然大幅增加至0.01 m左右，且有快速增加的趋势，说明此时隔水岩体已开始失稳，岩体强度已不具备足够的安全储备。因此，可以认为该工况下的安全厚度为7 m。

图7 弹性稳定时y向位移云图

图8 掌子面距溶洞17 m时z向位移云图

图9 掌子面距溶洞13 m时z向位移云图

图10 掌子面距溶洞7 m时z向位移云图

表4 监测点位移

图11 监测点处位移随开挖进程的变化曲线

3.2 破裂状态分析

由监测点位移分析结果显示，掌子面距溶洞7 m是位移的突变点，为隔水岩体的安全厚度，下面将重点研究该部分隔水岩体的破裂状态。隔水岩体随掌子面推进的破裂面状态如图12所示。未开挖时，隔水岩体相对完整，无明显破裂面；在距溶洞17 m时，开挖引起的应力重分布逐渐对岩体造成损伤，出现少量分散型破裂面，且随着掌子面的推进逐渐增多；至距溶洞7 m时，隔水岩体已产生大量破裂面，呈放射状分布。

根据上述分析结果，建立破裂度统计研究区域30 m×20 m×7 m，如图12(a)的红框所示。为分析开挖对研究区域的位移扰动，设立监测点于距溶洞3 m的掌子面中心。研究区域破裂度和监测点的z向位移随开挖进程的变化曲线如图13所示。

(a) 未开挖时

图13 破裂度和监测点z向位移随开挖进程变化曲线

计算结果显示：

1)破裂度和位移曲线均有突变点，分为2阶段。突变点左边为开挖影响较弱区，二者变化均很小；突变点右边为开挖影响强烈区，破裂度和位移均迅速增加。

2)从突变点开始，破裂度和位移均呈先缓慢增加后迅速增加，说明开挖引起的应力重分布对隔水岩体的损伤是逐渐加剧的，灾变过程发生前有孕育阶段。

3)破裂度突变点左边，破裂面占比小于1%，呈零星分布，此时隔水岩体处于微损伤阶段，如图12(b)和图12(c)所示；突变点右边，破裂面占比快速增加，呈放射状，掌子面处破裂面数量大于溶洞，失稳始发于掌子面。

4)破裂度突变点早于位移突变点6 m，说明隔水岩体破裂度达到一定程度后才会引起灾变，位移突变是灾变的表征现象。

5)距溶洞7 m处的破裂度约为8%，5 m处的破裂度约为14%，结合位移监测曲线可知，破裂度为8%时岩体损伤已处于临界状态。

3.3 不同水压下隔水岩层的破裂特征分析

当溶洞水压分别为1、2、3 MPa时，不同水压下监测点z向位移和破裂度随开挖进程变化曲线如图14和图15所示。计算结果显示：

1)不同水压作用下的破裂度突变点都在距离溶洞13 m处。突变点左边破裂度基本无变化；突变点右边，水压越大，破裂度增大，但增幅较小，说明此种工况下水压对隔水岩层破裂状态的影响不大。

2)不同水压作用下监测点位移的突变点都在距溶洞7 m处，突变点左边位移相对稳定，水压对远端位移的影响很小；突变点右边，水压越大，位移增幅越大，因此，可认为3种水压下的安全厚度均为7 m。

3)不同水压下破裂度突变点均早于位移突变点6 m；在突变点右边，二者均先缓慢增加后迅速增加。

图14 不同水压下监测点z向位移随开挖进程变化曲线

图15 不同水压下破裂度随开挖进程变化曲线

4 结论与建议

传统数值方法仅能模拟岩体的单一状态，且受软件本身计算能力的影响，以溶洞的空间二维分布形态为主，精度较差，计算耗时多。本文采用能够模拟隔水岩体渐进破坏的数值分析软件GDEM-DAS，辅助以GPU(显卡)加速，建立了百万网格的三维隧道数值模型，通过计算界面上弹簧的强度，模拟了隔水岩体的渐进破坏过程。通过监测掌子面的位移演化规律，获得隔水岩体的安全厚度，并利用算法的优势提出了定量描述岩体渐进破坏过程的无量纲指标——破裂度，分析了不同水压下隔水岩体的破裂度演化规律，以及与位移之间的对应关系。主要结论和建议如下：

1)突水为隔水岩体从裂纹萌生到贯通直至碎裂的渐进破坏过程，位移突变是破裂到灾变状态的表征现象。

2)本文给定的计算工况下，溶洞水压为1、2、3 MPa时隔水岩体的安全厚度均为7 m，由此可见溶洞水压非影响隔水岩体稳定性的首要因素，需结合围岩强度、地应力特征、溶洞大小、隧道跨度和埋深等其他因素进行分析。

3)由于地质体本身的非均匀性、非连续性，现场监测所得的地表位移和勘察数据大都是表面的、局部的，很难获取岩体内部复杂的破裂状态和真实的材料参数，由此开展的岩体稳定性计算的可靠性大大降低。因此，可考虑将地表位移监测结果、地表破裂状态与数值模拟结合，通过数值计算结果修正关键参数，获得与现场观测最为接近的一组数据作为材料的当前参数，并进一步计算地应力场和隔水岩体破裂状态，得到当前破裂面积，从而开展对隔水岩体危险程度的定量化评价。

本文仅针对均质围岩开展了隧道全断面的开挖模拟，未考虑岩体结构的复杂性和开挖方式的多样性，以及流固耦合作用的影响，具有一定的局限性。后续将分别针对节理裂隙隔水岩体破坏模式和溶洞水入渗诱发的隔水岩体渐进破坏效应开展进一步研究。