深埋隧道前方承压溶洞隔水岩体最小安全厚度研究

2021-09-26房忠栋杨为民王旌石锦江巴兴之王浩

中南大学学报（自然科学版） 2021年8期

房忠栋，杨为民，王旌，石锦江，巴兴之，王浩

(1.山东大学岩土工程中心，山东济南，250061；2.山东大学齐鲁交通学院，山东济南，250002；3.济南轨道交通集团有限公司，山东济南，250101；4.山东高速集团有限公司，山东济南，250098)

我国山区面积占国土总面积的2/3，并且是世界上岩溶分布最广的国家，西南地区最发育。近年来，交通事业的快速发展促进了隧道工程的蓬勃发展。在岩溶山区将不可避免地修建大量深埋隧道，岩溶对于隧道的危害主要是诱发突水涌泥，突水涌泥诱发地面次生灾害，如造成地面塌陷、井泉干涸等。岩溶突水已经成为隧道、矿山巷道等地下工程建设中常见的地质灾害之一。深埋隧道导致围岩自重应力处于较高水平，围岩所处环境的流−固强耦合效应不可忽略。一般认为，岩溶区隧道与岩溶构造间岩层防突安全厚度不够，导致隔水岩体在开挖扰动作用下垮塌或者在高水压作用下被击穿是引发岩溶隧道突水涌泥等工程事故的主要因素[1−5]。因此，必须对深埋岩溶隧道最小防突岩层安全厚度进行研究。

岩溶隧道施工中经常遇到不同发育程度的溶洞，隐伏溶洞的隐蔽性导致灾害往往具有突发性，严重威胁工程和人员安全[6]。在隐伏溶洞安全性研究方面，国内外学者针对溶洞的尺寸、分布位置、充填物性质、水压等因素对隧道稳定性的影响进行了研究，取得了一些实质成果。赖金星等[7]研究了不同位置的溶洞对隧道结构的影响，发现溶洞位置对于收敛值的影响由大至小依次为隧底溶洞、隧顶溶洞和侧壁溶洞，对于沉降值的影响由大至小依次为侧壁溶洞、隧顶溶洞和隧底溶洞，该研究未考虑溶洞含水状态下对隧道稳定性的影响。一些研究者通过构建耦合计算模型，借助数值模拟和模型试验等方法，针对富水隧道围岩稳定性与防突岩层安全厚度开展了大量研究，如：李术才等[8]开展了多种工况下隐伏溶洞突水灾变演化过程模拟试验，将承压隐伏溶洞突水过程划分为群裂隙萌生扩展、优势导水通道形成、隔水岩体破裂失稳3 个阶段；潘东东等[9]综合模型试验和数值模拟手段，针对强充填滞后型溶洞突水模式，基于隧道围岩稳定性受不同充水条件隐伏溶洞的影响规律，揭示了溶洞突水致灾机理；张桥[10]将围岩的位移变化情况和塑性区分布情况作为评价隧道围岩突水的判据，得出了隧道围岩防突岩层的最小安全厚度；李涛等[11]通过分析掌子面位移演化规律获得不同水压下隔水岩体的最小安全厚度，认为位移突变是破裂到灾变状态的表征现象；SHAN等[12]分析了成—桂铁路上隧道与周围隐伏洞之间的安全厚度和岩石破碎机理，提出了基于玉景山隧道工程的综合数值分析方法；王志杰等[13]通过研究隧道前方存在正交高压溶洞时中间岩墙的承压能力和破坏模式，建立了掌子面失稳破坏的圆锥台模型，并通过势能判据的尖点突变理论得到掌子面失稳时的溶洞临界压力。

现有研究主要集中在隧道洞身范围外存在隐伏溶洞时围岩结构的稳定性、防突岩层安全厚度及突水机理等方面。但当溶洞位于隧道前方时，随着掌子面的推进，溶洞必然会被揭露，更易导致突水灾害的发生。深埋岩溶隧道面临高应力和强水压的耦合作用，对隔水岩体最小安全厚度研究尚不够深入，考虑隧道开挖扰动与流−固耦合综合作用的分析较少。针对上述问题，本文以重庆歇马隧道为工程背景，利用COMSOL Multiphysics软件对流−固耦合计算和隧道开挖过程进行模拟分析，确定深埋隧道前方隐伏承压溶洞隔水岩体最小安全厚度，并分析溶洞尺寸、水压以及隧道埋深、跨度等因素对最小安全厚度的影响规律，最后通过多元线性回归方法推导出岩溶隧道前方承压水溶洞掌子面预留安全厚度公式，以期对类似工程施工安全提供技术支撑。

1 工程背景

歇马隧道是重庆快速路工程中“一横线”工程的控制性工程，设计为分离式双洞双向六车道行驶隧道，隧道开挖净宽为14.5 m，净高为9.6 m。隧道全长4 200 m，最大隧道埋深为392 m，属于深埋特长隧道。隧道穿越以观音峡背斜轴部隆起为主的中梁山，山顶部受碳酸盐岩溶蚀作用的影响，东西方向形成了刘家槽和窑湾共2 条岩溶槽谷，呈“一山两槽三岭”形态。在2 条岩溶槽谷段，地表及地下水丰富，隧道存在岩溶突水涌水等地质灾害，地质剖面图如图1所示。

图1 歇马隧道地质剖面图Fig.1 Geological profile of Xiema Tunnel

当隧道右线开挖至里程桩号YK9+268 时，围岩体较破碎，节理裂隙较发育，经超前钻探得掌子面前方有1 个体积约(长×宽×高)为5 m×10 m×15 m的充水溶洞。根据地质勘探报告显示，YK9+266—YK9+296 段隧道平均埋深235 m，饱和抗压强度为34.36 MPa，完整性系数为0.68，综合判定为Ⅲ级围岩。静水压力为1.8～2.0 MPa，地下水活动状态多为大股涌水或喷射状涌水。

2 突水计算模型的构建

深埋隧道前方存在承压水溶洞时，随着开挖扰动作用的影响，围岩内部渗流场和应力场变化明显，且相互作用，表现为2 个物理场的强耦合[14]：一方面，高水压作用引起溶洞附近围岩的应力场变化；另一方面，隧道的开挖导致围岩的应力释放，围岩孔隙率和渗透率发生改变，从而改变渗流场的分布。在开挖扰动和渗流作用的双重作用下，隧道围岩强度逐渐降低，当开挖扰动损伤带H1和水压致裂带H2连通时，即可认为隧道发生突水，此时的岩层厚度为最小安全厚度HS，最小安全厚度示意图如图2所示。

图2 最小安全厚度示意图Fig.2 Schematic diagram of minimum safe thickness

2.1 模型构建与参数选取

本文使用COMSOL Multiphysics软件中的固体力学模块和达西定律模块对流−固耦合作用时隧道开挖过程进行模拟计算分析。引入以下基本假设：1)围岩材料服从Mohr−Coulomb屈服准则；2)地下水渗流发生时服从Darcy定律；3)溶洞经过长期溶蚀作用已达到稳定状态，溶腔内充满水体。

以歇马隧道突水段为工程背景，建立数值计算模型。为了消除边界效应，整个模型长为120 m，宽为60 m，高为80 m；隧道开挖跨度为14.5 m，高度为9.6 m，高跨比为0.66。模型采用自由四面体进行网格剖分，总单元数10 218 个，具体计算模型如图3所示。根据地质勘察资料和岩石物理力学参数测定试验，确定计算模型中围岩物理力学参数如表1所示。

表1 围岩的物理力学参数Table 1 Physical and mechanical parameters of surrounding rock

图3 计算模型Fig.3 Calculation model

应力场边界条件为：在岩体应力场变形计算中，模型的四周边界设为辊支撑约束，底部边界为固定约束，模型顶部和隧道开挖临空面均设定为自由面；模型内不考虑构造应力，初始地应力按照自重应力计算，侧压力系数取1；隧道埋深为235 m，构建的模型上表面至隧道拱顶为40 m，剩余195 m 岩体自重作为外部荷载作用于模型上表面，计算后为5 MPa。

渗流场边界条件为：在多孔介质渗流场计算中，模型的外边界均设定为不透水边界，质量通量为0 kg/(m2·s)，溶洞内部水压设置为2 MPa，隧道开挖轮廓边界设定为透水边界，压力水头为0 m。

2.2 开挖过程实现

将隧道轮廓内开挖岩体沿隧道轴线分割为多个柱体，设置多个计算步骤，每个计算步骤包含一组物理场(固体力学和达西定律)。每步计算开始时选择计算范围并设定参数。将待开挖的计算域钝化，使其不参与计算，以实现隧道的开挖过程模拟。计算域分布见图4。具体操作如下：

图4 计算域分布图Fig.4 Computational domain distribution map

1)求解固体力学和达西定律物理场信息，计算范围为整个模型体，模拟隧道开挖前的平衡状态，即计算模型体的初始应力状态。

2)以步骤1)的计算结果作为该步骤的初始值，计算范围为去掉计算域2的模型体，即计算域2不参与计算，模拟第一步开挖。

3)以步骤2)的计算结果作为该步骤的初始值，计算范围为去掉计算域2和计算域3的模型体，即计算域2和计算域3 不参与计算，模拟第二步开挖。

4)依此类推，直至隧道开挖范围内的所有计算域全部删除，完成隧道的开挖过程模拟。

3 突水计算结果分析

通过COMSOL Multiphysics计算模拟歇马隧道突水段开挖过程，分析隧道开挖完成后围岩位移变形分布规律，同时进一步分析考虑耦合作用下开挖后围岩渗透压力、渗流速度和塑性区的分布特征，总结隔水岩层最小安全厚度确定方法，并通过模型试验结果对模拟结果的有效性进行验证。

3.1 数值计算结果

综合分析隧道开挖过程中围岩渗流场与应力场的变化规律，结合开挖扰动和水压作用下主应力差和塑性区分布特征，建立隔水岩层最小安全厚度确定方法。

3.1.1 渗流场变化规律

取隧道轴线所在竖直截面作为研究对象，围岩渗流场分布如图5所示(图中箭头的大小代表渗流速度的大小)。由图5(a)可见：隧道开挖后导致掌子面临空，此时隧道轮廓为透水边界，故在掌子面附近出现较大的渗流速度集中区，但由于掌子面距离溶洞较远，掌子面与溶洞间的岩体内渗流速度很小，故在掌子面部位不会发生突水灾害。由图5(b)～(f)可见：随着隧道逐步开挖，隔水岩体厚度逐渐降低，溶洞水流动所受阻力减小，围岩内的渗流速度逐渐增加，溶洞内承压水引起的渗流场逐渐扩展至掌子面，同时，溶洞除了靠近掌子面一侧的边界外，其他方向的渗流速度均逐渐减少，说明溶洞内承压水逐渐通过掌子面前方岩体向隧道内流入。

图5 开挖过程中围岩渗流场分布图Fig.5 Distribution of seepage field in surrounding rock during excavation

3.1.2 塑性区分布规律

溶洞水压和隧道开挖扰动均导致溶洞周围或隧道边界附近形成塑性区，图6所示为隧道轴线截面围岩塑性区云图。由图6(a)～(c)可见：随着掌子面逐渐靠近溶洞，原本开挖扰动控制的塑性区逐渐受溶洞的影响，该部分塑性区有增大的趋势；溶洞承压水控制的塑性区首先在溶洞边界附近形成，在隧道的开挖过程中逐渐向掌子面方向发展。由图6(d)可见：当隧道开挖至溶洞前方4 m 位置时，掌子面附近塑性区和溶洞附近塑性区发生贯通，认为此时突水通道已经形成，即该施工工况对应的隔水岩体最小安全厚度在4～5 m之间。

图6 开挖过程中围岩塑性区分布图Fig.6 Distribution of plastic zone in surrounding rock during excavation

3.1.3 隔水岩体最小安全厚度确定

围岩主应力差与围岩破坏的发生密切相关，相比于最大(小)主应力更能反映围岩状态[15]。通过模拟计算得到的隧道开挖过程中掌子面主应力差(最大主应力和最小主应力之差)与掌子面最大水平位移，其量纲和取值范围存在较大差异。为了更方便分析这2个物理量变化产生的影响，将它们归一化处理，并绘制于同一坐标系中，如图7所示。由图7可知：当掌子面与溶洞距离d>12 m 时，掌子面上主应力差和最大水平位移增长缓慢，即当掌子面距离溶洞较远时，掌子面受溶洞影响较小，围岩内的流−固耦合作用不明显，隧道开挖过程的扰动作用是引起应力和位移变化的主要因素；当5

图7 开挖过程中掌子面信息曲线图Fig.7 Curves of tunnel face information during excavation

模拟过程中靠近溶洞时的每循环开挖进尺为1 m，综合以上模拟结果，认为隔水岩体最小安全厚度在4～5 m 之间。因此，改变最后一步开挖进尺，以0.1 m 为最小分度值，设计多组计算工况，重新模拟以上过程，以开挖扰动引起的塑性区和溶洞周围塑性区贯通时掌子面与溶洞边界距离作为隔水岩体的最小安全厚度，最终确定该工况最小安全厚度为4.7 m，此时的塑性区分布如图8所示。

图8 塑性区分布图Fig.8 Distribution of plastic zone

3.2 溶洞突水模型试验对比

文献[16]依托歇马隧道工程突水实例，采用1:50的模型相似比尺，开展了隧道前方充水溶洞突水灾变演化模型试验，总结了溶洞突水的3个阶段即导水裂隙萌生、导水裂隙扩展和导水裂隙贯通，隧道开挖完成后，掌子面先后经历了滴状渗水—股状流水—突水等过程，见图9。

图9 试验突水过程图Fig.9 Diagrams of water inrush process during test

综合试验模型体位移场、应力场和渗流场信息变化规律，最终确定该工况下隔水岩体最小安全厚度为4.5 m。隔水岩体最小安全厚度数值模拟结果与该试验结果相对误差为4.4%，因此，本文所选数值模拟方法与计算模型可用于类似工程隔水岩体最小安全厚度的研究。

4 不同因素对隔水岩体最小安全厚度的影响

隔水岩体最小安全厚度的影响因素较多，并且需考虑具体的施工工况。为研究不同因素的影响效果，本文采用正交试验进行研究。正交试验设计是根据正交性从全面试验中挑选部分具有代表性的点进行研究，是研究多因素多水平的试验设计方法[17]，具有均匀分散、齐整可比的优点，即通过正交试验设计的各因素水平组合在全部组合中的分布是均匀的，并且每个因素的各水平之间具有可比性。

4.1 正交模拟试验设计

基于岩溶隧道隔水岩层最小安全厚度的研究成果[18−20]，研究溶洞水压P、溶洞直径D、隧道跨度L和隧道埋深h这4 个因素对隔水岩层最小安全厚度的影响规律，各因素分别选择4个水平，采用如表2所示的L16(44)正交试验表设计模拟工况，具体参数设置见表3。

表2 L16(44)正交试验表Table 2 L16(44)orthogonal test table

表3 因素水平设置Table 3 Factor level setting

4.2 模拟结果分析

4.2.1 不同工况最小安全厚度计算结果

基于4.1节中正交试验表(表2)和各因素水平设置(表3)规则，确定各模拟工况的参数取值，开展16 组模拟试验。根据隔水岩层最小安全厚度确定方法，确定各模拟工况隔水岩体最小安全厚度，如表4所示。

表4 最小安全厚度模拟结果Table 4 Simulation results of minimum safe thickness

4.2.2 各因素对最小安全厚度的影响规律

各个因素对最小安全厚度的影响程度各异，通过对上述模拟结果进行极差分析，确定其主次顺序和重要程度，并分析各因素单独作用效应。

1)极差分析。最小安全厚度的模拟结果极差分析如表5所示，其中，Ki为某个影响因素的第i个水平所对应模拟结果的总和，ki为某个影响因素的第i个水平模拟结果的平均值，Rj为第j列因素的各个水平所对应模拟结果的平均值的极差。Rj反映的是第j列因素的水平变动引起的观测指标的变化幅度，Rj越大，则说明该因素对试验指标的影响越大。通过对深埋隧道前方隐伏承压溶洞的隔水岩体最小安全厚度模拟结果的极差分析，可得出如

表5 模拟结果极差分析Table 5 Range analysis of simulation results

下关系：R溶洞水压>R隧道埋深>R溶洞直径>R隧道跨度，即认为各因素中，对隔水岩体最小安全厚度影响最大的是溶洞水压，其次是隧道埋深，影响最小的是隧道跨度。

2)各因素影响趋势。根据隔水岩体最小安全厚度极差分析结果分别绘制各因素对最小安全厚度的影响规律曲线，如图10所示。由图10(a)可见：当其他因素确定时，最小安全厚度随溶洞水压增大而增大，并且水压越大，最小安全厚度增大幅度越大。这是因为隧道突水通道是溶洞周围岩体受水压作用产生的破碎区和隧道开挖扰动产生的破碎区贯通而形成的，溶洞内的水压越大，则水压致裂产生的破坏区域越大，从而需要更大的隔水岩体厚度来抵抗洞内水压的作用。

由图10(b)可见：当其他因素确定时，最小安全厚度随溶洞直径增大而增大；隧道开挖前的地层各物理场均处于平衡状态，每一步隧道开挖步都对隧道周围岩体造成一定的扰动，围岩通过一定变形形成新的平衡状态，而溶洞的存在放大了该效应，溶洞越大，隧道开挖引起的溶洞周边物理场变化范围越大，对隔水岩体的稳定性越不利，从而导致最小安全厚度越大。

图10 各因素对最小安全厚度的影响Fig.10 Influence of various factors on the minimum safe thickness

由图10(c)可见：当其他因素确定时，在隧道高跨比不变时，随着隧道跨度增大，最小安全厚度也有增大趋势，但相较于溶洞水压和溶洞直径而言，隧道跨度对最小安全厚度的影响效应并不突出。在通常情况下，由于岩体结构面的存在，围岩压力随着隧道跨度增大而增大。本文采用的有限元模拟方法将岩体假定为连续体，未考虑结构面的影响，因此，最小安全厚度并未随隧道跨度增大而大幅度增加。

由图10(d)可见：当其他因素确定时，隔水岩体最小安全厚度随着隧道埋深增大而减小。这是因为隧道埋深越大，开挖前隧道原平衡状态中围岩的地应力越大，对溶洞水压致裂效应起到抑制作用。因此，模拟结果显示，当溶洞内水压和溶洞直径一定时，隧道埋深越大，溶洞周围围岩因水压致裂导致的破坏区越小，从而导致最小安全厚度越小。

4.3 隔水岩体最小安全厚度预测模型

多元线性回归分析方法是指对2个或2个以上的事件因素设置成自变量，而将事件结果设置为因变量，因变量与每一个自变量呈线性变化关系，通过相关性分析，建立因变量与若干自变量的线性关系式，由自变量变化进行回归分析。通过采用多元回归分析方法对正交数值模拟结果进行分析。

4.3.1 一元回归分析

以影响因素溶洞水压P、溶洞直径D、隧道跨度L、隧道埋深h分别作为自变量，最小安全厚度HS作为因变量，选择基本的回归方程(线性函数、指数函数、对数函数、幂函数等)分别进行相关性分析，并取拟合优度最大的函数作为其对应的回归方程。

1)溶洞水压P与最小安全厚度HS的关系：

拟合优度R²=0.916 4。

2)溶洞直径D与最小安全厚度HS的关系：

拟合优度R²=0.930 6。

3)隧道跨度L与最小安全厚度HS的关系：

拟合优度R²=0.934 9。

4)隧道埋深h与最小安全厚度HS的关系：

拟合优度R²=0.981 6。

以上公式的拟合优度均大于0.900 0，拟合效果较好，可认为所得经验公式有意义。

4.3.2 多元回归分析

根据上述的单因素回归关系，建立岩溶隧道掌子面前方隔水岩层最小安全厚度的预测模型。首先假设最小安全厚度HS与e0.0937P，ln(L)，D和h之间存在如下线性关系：

式中：ai(i=1，2，…，5)为待定系数。

通过多元线性回归分析方法求出各待定系数，得最小安全厚度预测公式：

式中：1.5≤P≤6.0 MPa；9≤L≤18 m；4≤D≤13 m；300≤h≤900 m。

对式(6)进行显著性检验，验证回归方程的线性相关性。式(6)对应的方差分析结果表明：F服从自由度为(4,11)的分布，由F分布检验表查得F0.05(4,11)=3.357，F0.01(4,11)=5.668；F>F0.01(4,11)，并且P<0.001，说明在显著性水平α=0.01 下，HS与e0.0937P，lnL，D和h之间存在显著的线性关系，即回归方程显著，F检验通过。