巷道前伏承压溶洞突水灾变流固耦合分析

赵延林 ，曹平，万文 ，王卫军 ，张盛国 ，邹声华

(1. 湖南科技大学 煤矿安全开采技术湖南省重点实验室，湖南 湘潭，411201；2. 湖南科技大学 能源与安全工程学院，湖南 湘潭，411201；3. 中南大学 资源与安全工程学院，湖南 长沙，410083)

突水事故造成的直接经济损失一直排在各类煤矿灾害前列。南方煤矿突水事故的水源主要有岩溶水、老窑积水和地表水，而且三者之间连通关系复杂，往往其中一类突水导致其他突水事故并发[1-3]。在各类水源中, 地下岩溶突水压力大, 突发性特别强, 造成的破坏损失最严重。矿井开挖过程中，当击穿大型充水岩溶体(如暗河、填充溶洞及岩溶破碎带等)时，会引发矿井涌水突水地质灾害[4-5]。如2003 年湖南涟源市七一煤矿石坝井的“4·16”特大突水灾害(17 人遇难)、2010 年湖南张家界市鑫茂煤矿“9.28”溶洞透水事故、2004 年贵州思南县天池煤矿“12.12”特大透水事故(21人死亡，15 人失踪)都是晚二叠系茅口灰岩内高承压溶洞水体突破岩层引起的突水事故。对于隐伏溶洞突水机理的研究主要集中在西南岩溶隧道突水机理的研究上。刘招伟等[6-8]分别对圆梁山隧道、沪蓉西高风险岩溶隧道、季家坡隧道的溶洞突水致灾机制和围岩、衬砌破坏机理进行了研究，总结了岩溶突水的不同地质模式、基本力学机制、突水流态演化模型和衬砌压裂致灾机制。目前，岩溶突水的力学机理研究主要针对岩溶隧道，而对于岩溶矿井在采掘过程中由于溶洞承压水击穿岩柱而产生工作面岩溶突水的力学机理和突水过程分析研究较少。在岩溶矿井中，一般通过布置探放水钻孔来探放岩溶水，通过防突岩柱的留设来防治岩溶突水。而《煤矿防治水规定》[9]没有具体给出针对隐伏承压溶洞防突岩柱留设的安全厚度计算方法。进行隐伏承压溶洞防突岩柱的稳定性及安全厚度研究，对于防治矿井岩溶突水具有重要意义[10-11]。本文作者从岩石力学和渗流力学出发采用流固耦合分析方法探讨巷道前伏承压溶洞防突岩柱在采动应力和高水力梯度作用下的流固耦合效应，研究高水压下防突岩柱的失稳演化，同时引入强度折减法理论提出防突岩柱安全系数的概念，分析巷道前伏承压溶洞防突岩柱的安全厚度，为合理设计和预留防突岩柱提供理论依据。

1 承压溶洞突水的流固耦合机理

当巷道前伏承压溶洞时，采动应力诱发围岩应力场与渗流场的急剧变化，在溶洞与工作面之间的围岩部位将会显现应力场和渗流场的强耦合效应。这种耦合作用表现在：高水力梯度(高动水压力)作用于围岩导致应力场的改变；在高静水压力和应力场共同作用下围岩裂隙扩展，其渗透特性大大加强。从岩体流体力学角度出发，承压溶洞突水可认为是采动应力和高动、静水压力作用下，溶洞与工作面之间岩柱的突变失稳并伴随岩溶水大规模涌出。防突岩柱在采动应力与渗流的双重作用下强度不断降低，当系统状态濒临临界状态时，微小的扰动诱发防突层整体破断突水。

1.1 非线性渗流分析

防突岩柱渗流分析微分控制方程为[12]：

式中：p 为渗透压；k 为渗透系数；M 为比奥模量；ε为体积应变；α 为有效应力的比奥系数。

渗透压和应力对岩体渗透性的影响如下[13]

式中：k0为岩体渗透系数初值； k(Θ , p)为耦合分析下岩体的渗透系数；Θ = σ1+ σ2+ σ3，为体积应力；σ1，σ2， σ3分别为第一、第二、第三有效主应力；ξ 为渗透突跳系数。

在耦合分析中，对于处于弹性状态的单元，认为岩体渗透系数是有效体积应力的负指数函数，取ξ=1.0，β=0.5。

对于塑性屈服单元，渗透系数大大增加，渗透突跳系数显著提高，取ξ=5 000，β =1.0。

式(2)反映了岩体应力对渗流的耦合作用，尤其反映了塑性单元对渗流的强耦合效应。

1.2 弹塑性力学分析

应力场与渗流场耦合作用下岩体的弹塑性力学本构关系[13-14]：

式中：(α p),j为渗透水力梯度作为等效体积力作用于岩体骨架，反映了渗流场对应力场的耦合作用。

屈服准则采用Mohr-Coulomb 准则。该准则的剪切破坏判据为：

拉伸破坏判据为：

式中： σt为岩体抗拉强度；φ 为岩体内摩擦角；c 为黏聚力。

1.3 流固耦合程序设计

巷道承压溶洞突水的流固耦合分析，采用间接耦合手段，首先计算出ti时刻岩柱的弹塑性应力场，由式(2)计算出ti时刻岩柱渗透系数，将该渗透系数传递给渗流场计算模块，实现应力场对渗流的耦合；再将渗流场计算得到的渗透体积力施加给弹塑性应力计算单元，完成渗流场对应力场的耦合。如此对时间循环，直至达到预定的渗流计算时刻。为研究防突岩柱的流固耦合效应，在FLAC3D2.1 的平台下，通过FISH 语言成功开发了巷道前伏承压溶洞防突岩柱流固耦合分析程序。该程序分为弹塑性应力计算、渗流计算和耦合分析3 个模块。

基于FLAC3D开发的流固耦合间接分析方法，其优点体现在各耦合参数(如渗透系数、渗透体积力等)随计算时步时空演化和各耦合子系统(渗流子系统和弹塑性力学子系统)的耦合信息随计算时步准确传递。

1.4 流固耦合-强度折减法联动分析

采用强度折减法思想研究防突岩柱的稳定性。

强度折减法的基本原理是将岩体的抗剪参数c,tanφ 同时除以某个折减系数Fi, 得到一组新的参数c′，φ′, 将新的材料参数作为计算参数进行试算。

将承压溶洞突水的流固耦合理论和强度折减法结合起来，形成流固耦合-强度折减法联动分析方法。研究在防突岩柱强度参数不断折减的情况下，岩柱失稳突出和渗流突变的演化过程。

流固耦合-强度折减法联动分析中，首先进行弹塑性应力-渗流的流固耦合分析，判断计算是否收敛，若收敛则对防突岩柱的强度参数进行折减，并记录折减系数，用折减后强度参数作为防突岩柱的强度，进行流固耦合分析，再判断计算是否收敛，如此重复，直至计算不收敛，此时的折减系数即为防突岩柱的安全系数。流固耦合-强度折减法联动分析流程图如图1所示。

图1 流固耦合-强度折减法联动分析流程图Fig.1 Flow diagram of Linkage analysis between fluid-solid coupling and strength reduction method

2 突水案例

以湖南省娄底涟源市七一煤矿石坝井-160 m 水平老水仓扩容掘进工作面发生岩溶突水事故为案例，进行防突岩柱的稳定性分析。该矿井水文地质条件复杂，煤层下伏茅口灰岩层，在茅口灰岩中溶洞往往成串珠状存在，溶洞呈半充填状态，充填物含泥砂。23采区-160 m 水平老水仓扩容掘进工作面位于茅口灰岩之中，2003 年4 月15 日在掘进过程中发生了承压溶洞突水的重大事故，16 人遇难，28 h 共计涌水约12 500 m3，涌泥5 440 m3。清碴后，在工作面前方3～4 m 处露出1 个10 m×6 m×10 m 左右的溶洞。

2.1 计算模型

图2 所示为七一煤矿石坝井-160 m水平老水仓扩容掘进防突岩柱稳定性分析的计算模型，模型尺寸x×y×z=50 m×15 m×30 m，掘进巷道位于模型中部，半圆拱形巷道，模型右边界设有承压溶洞，为简化假定溶洞为近椭球形，老水仓扩容巷道埋深450 m。

力学边界及网格划分见图2，其中，y=0 平面为固定约束，σ1=12.5 MPa，σ3=10 MPa。渗流边界条件：茅口灰岩内溶洞水压p 分别取2，3 和4 MPa，工作面水压为0.1 MPa。

图2 中，1～5 为模拟分步掘进。在数值计算中，假定每步进尺5 m。当工作面距溶洞只有10 m 时，为研究不同厚度防突岩柱的稳定性，假定每步进尺1 m，以防突岩柱厚度分别为10，9，8，7，6，5 和4 m 为例研究不同厚度防突岩柱的稳定性。

对石坝井-160 m 水平的茅口灰岩层进行取样，并进行岩石力学和渗透性实验，得到该类岩石的强度参数和渗流力学参数见表1。

采用流固耦合-强度折减法联动分析法对承压溶洞突水过程展开研究。以溶洞内压p=4 MPa，防突岩柱h=8 m 为例来研究不同强度下防突岩柱失稳演化的突水突变过程。

图2 计算模型Fig.2 Numerical calculation model

2.2 联动分析结果

随着折减系数的增大，防突岩柱强度参数削减，在采动应力和高水压力作用下，防突岩柱的塑性区扩展演化，导致岩柱的渗透系数显著增加。图3 所示为折减系数Fi=1.10 和1.54 时防突岩柱的渗透系数分布图。Fi=1.10 时，岩柱强度参数为c=1.36 MPa，tan φ=0.52，岩柱渗透系数主要受体积应力和孔隙压的影响，岩柱渗透系数自工作面向溶洞形成了工作面附近的渗透增强区、岩柱中部的未扰动渗透区及溶洞附近的渗透增强区。距工作面1.5 m 左右的渗透增强区，渗透系数为3.25×10-10～15.21×10-10m/s；岩柱中部的未扰动渗透区，其渗透系数为0.70×10-10～1.79×10-10m/s；溶洞附近的渗透增强区，渗透系数为2.97×10-10～5.68×10-10m/s。Fi=1.54 时，岩柱强度参数为c=0.97 MPa，tan φ=0.39，此强度下岩柱内85%左右的单元处于塑性状态，防突岩柱渗透系数成数量级增加，为1.47×10-6～2.80×10-6m/s。在此强度下的岩柱不具备防突水的能力。

图3 防突岩柱的渗透系数分布Fig.3 Distribution of hydaulic conductivity of rock pillar for preventing water inrush

表1 岩体力学计算参数Table 1 Mechanical parameters of rock mass

图4 所示为折减系数 Fi=1.10 时，计算模型的渗流场分布。巷道隐伏承压溶洞，掘进扰动了围岩渗流场，岩溶水在高水力梯度作用下向临空面排出，渗流场不断调整以消除承压溶洞高势能不稳定状态。在岩柱内形成了自溶洞向工作面逐渐降低的渗透压，聚集于岩柱的平均渗透压梯度高达0.5 MPa/m，高渗透压梯度产生指向工作面的水平体积力是导致防突岩柱突出失稳的动力。

图4 折减系数Fi =1.10 时渗透场分布Fig.4 Distribution of seepage field when Fi =1.10

图5 不同折减系数下的工作面的渗流量Fig.5 Discharge of working face under different reduction factors

在流固耦合分析中，通过记录工作面的流量，得到不同折减系数下工作面的渗流量见图5。由图5 可知：伴随着防突岩柱强度的减弱，岩柱工作面渗流量也由小流量演化成大规模的涌水量。Fi=1.10 时，工作面流量为6.158×10-4m3/s，此时该强度下的岩柱稳定，岩柱内塑性区范围不大，工作涌水流量少；Fi=1.54～1.58 时，岩柱内塑性单元占85%～87%，与塑性单元强耦合的岩柱渗透系数突跳增长，导致工作涌水量急剧增加，工作面涌水量高达0.753～ 0.811 m3/s，当Fi=1.59 时，防突岩柱全部处于塑性状态，计算不收敛，工程上则表现为承压溶洞突水。

在采动应力和渗透体积力共同作用下岩柱向工作面产生位移，随着折减系数的增加，其位移量也随之增大。图6 所示为当Fi=1.54 时，防突岩柱水平位移等值线图，耦合计算得出：当Fi=1.54 时，防突岩柱向工作面位移高达14.58 cm。而当Fi=1.59 时，力学计算不收敛，此时防突岩柱向工作面突出失稳。

图6 折减系数Fi =1.54 时防突岩柱中剖面上水平位移云图Fig.6 Horizontal displacement diagrams of rock pillar when Fi =1.54

2.3 防突岩柱安全系数

采用流固耦合-强度折减法联动分析，改变溶洞内压p 和岩柱厚度h，得到水压p=2，3 和4 MPa 时，防突岩柱计算厚度h 分别为3，4，5，6，7，8，9 和10 m 时，防突岩柱的安全系数和溶洞水压、防突岩柱厚度的关系曲线(图7)。

由图7 可知：在一定溶洞水压下，随着防突岩柱厚度的增大，岩柱安全系数也逐渐增大，采用指数函数对其进行拟合：

式中：Fs为防突岩柱安全系数；h 为防突岩柱计算厚度；a，b，d 为待定系数。

由式(9)可拟合得到不同溶洞水压下防突岩柱的安全系数与岩柱计算厚度之间关系式，见表2。将Fs=1.5 作为评价防突岩柱稳定性的判断标准，可得出如下判断标准[15]：

(1) Fs≥1.5 时防突岩柱稳定；

(2) Fs＜1.5 时防突岩柱不稳定。

图7 防突岩柱安全系数和岩柱计算厚度的关系曲线Fig.7 Variation relation curves between safe factor of rock pillar for preventing water inrush and rock pillar thickness

表2 不同溶洞水压下防突岩柱安全系数与岩柱计算厚度的关系Table 2 Relationship between safety factor of rock pillar for preventing water inrush and rock pillar thickness under different karst cave water pressures

由表2 可见：随溶洞水压力增大，安全防突岩柱计算厚度要相应增加，如溶洞水压为2 MPa 和4 MPa时，安全防突岩柱计算厚度分别为5.8 m 和8.0 m，后者的安全防突岩柱厚度要比前者多2.2 m。这是由于溶洞水压的增加，聚集在防突岩柱内的渗透体积力增大，岩柱受到向工作面方向的侧向推力增加，因此安全防突岩柱要相应增加留设厚度。

在七一煤矿石坝井-160 m 水平老水仓扩容掘进工作面发生的溶洞高承压水体突破岩柱涌出的重大突水事故中发现：10 m×6 m×10 m 的溶洞距工作面3～4 m，即工作面与溶洞之间只有3～4 m 的防突岩柱，且岩溶水与地表水沟通良好，溶洞内水压为3.5～4.0 MPa。防突岩柱工程留设厚度至少应该达到炮眼长度(约2.2 m)和安全防突岩柱计算厚度(7～8 m)之和，为9.2～10.2 m。3～4 m 厚的防突岩柱完全不具备有防突功能。

3 结论

(1) 巷道隐伏承压溶洞突水的流固耦合模型揭示突水机理在于采动应力和高渗流体积力共同作用于防突岩柱，导致岩柱塑性区扩展而突出失稳，伴随岩柱的力学失稳，防突岩柱发生渗流突变，工作面发生大规模突水。

(2) 采用承压溶洞突水的流固耦合-强度折减法联动分析方法，研究高水力梯度作用下防突岩柱的安全储备，提出防突岩柱安全系数的概念，为安全防突岩柱留设厚度提供了一种新的定量方法。

(3) 将具有安全系数为1.5 的防突岩柱厚度作为防突岩柱的计算安全厚度。随溶洞水压的增加，防突岩柱的计算安全厚度相应增大。

(4) 七一煤矿石坝井-160 m 水平老水仓扩容掘进工作面发生的溶洞突水是由于防突岩柱留设过小，不具有安全储备，导致溶洞高承压水体突破岩柱涌出。

[1] 李敬敬, 常聚才, 陈阳. 厚松散含水层下特厚煤层综放开采相似模拟[J]. 矿业工程研究, 2013, 28(4): 47-50.LI Jingjing, CHANG Jucai, CHEN Yang. Similar simulation research on caving mining of extra thick coal seam under thick loose aquifer[J]. Mineral Engineering Research, 2013, 28(4):47-50.

[2] 李文敏. 承压水上采煤的流固耦合数值模拟[J]. 矿业工程研究, 2013, 28(2): 37-41.LI Wenmin. The fluid-solid coupling numerical simulation of confined water mining[J]. Mineral Engineering Research, 2013,28(2): 37-41.

[3] 李建明, 朱斌, 武强. 赵各庄矿大倾角煤层综放开采突水条件的力学机制分析[J]. 煤炭学报, 2007, 32(5): 453-457.LI Jianming, ZHU Bin, WU Qiang. Mechanical mechanism analysis of water inrushcondition on mechanized top-coal caving in steep seam in Zhaogezhuangcoal mine[J]. Journal of China Coal Society, 2007, 32(5): 453-457.

[4] 杨天鸿, 唐春安, 谭志宏, 等. 岩体破坏突水模型研究现状及突水预测预报研究发展趋势[J]. 岩石力学与工程学报, 2007,26(2): 268-277.YANG Tianhong, TANG Chunan, TAN Zhihong, et al. State of the art of inrush models in rock mass failure and developing trend for prediction and forecast of groundwater inrush[J].Chinese Journal of Rock Mechanics and Engineering, 2007,26(2): 268-277.

[5] Zhang J C. Investigations of water inrushes from aquifers under coal seam[J]. International Journal of Rock Mechanics and Mining Sciences, 2005, 42(2): 350-360.

[6] 刘招伟, 何满潮, 王树仁. 圆梁山隧道岩溶突水机理及防治对策研究[J]. 岩土力学, 2006, 27(2): 228-232.LIU Zhaowei, HE Manchao, WANG Shuren. Study on karst water burst mechanism and prevention countermeasures in Yuanliangshan tunnel[J]. Rock and Soil Mechanics, 2006, 27(2):228-232.

[7] 李利平, 李术才, 张庆松. 岩溶地区隧道裂隙水突出力学机制研究[J]. 岩土力学, 2010, 31(2): 523-528.LI Liping, LI Shucai, ZHANG Qingsong. Study of mechanism of water inrush induced by hydraulic fracturing in karst tunnels[J]. Rock and Soil Mechanics, 2010, 31(2): 523-528.

[8] 许振浩, 李术才, 李利平, 等. 一种典型的岩溶隧道衬砌压裂突水灾害成因与防治[J]. 岩石力学与工程学报, 2011, 30(7):1396-1404.XU Zhenhao, LI Shucai, LI Liping, et al. Cause, disaster prevention and controlling of a typical kind of water inrush and lining fracturing in karst tunnels[J]. Chinese Journal of Rock Mechanics and Engineering, 2011, 30(7): 1396-1404.

[9] 国家安全生产监督管理总局, 国家煤矿安全监检局. 煤矿防治水规定[M]. 北京: 煤炭工业出版社, 2009.Sate Administration of Work Safety, Sate Administration of Coal Safety. Rule of mine prevent ion and cure water disaster[M]. Beijing: China Coal Industry Publishing House,2009.

[10] 高保彬, 高佳佳, 袁东升. 基于UDEC 的大采高覆岩破裂的模拟与分析[J]. 湖南科技大学学报(自然科学版), 2013, 28(2):83-87.GAO Banbin, GAO Jiajia, YUAN Dongshen. Simulation and analysis of overlying strata fracture in lame mining height top-coal caving face of coal seams based on UDEC[J]. Journal of Hunan University of Science & Technology (Natural Science Edition), 2013, 28(2): 83-87.

[11] 刘超, 张义平, 齐飞, 等. 地面直流电法在煤矿防治水中的应用[J]. 矿业工程研究, 2012, 27(3): 64-67.LIU Chao, ZHANG Yiping, QI Fei, et al. The application of ground DC method in mine water prevention and control[J].Mineral Engineering Research, 2012, 27(1): 54-58.

[12] Itasca Consulting Group Inc. FLAC3D users manual[R].Minneapolis: Itasca Consulting Group Inc., 2004.

[13] 王卫军, 赵延林, 李青锋, 等. 矿井岩溶突水灾变机理[J]. 煤炭学报, 2010, 35(3): 443-448.WANG Weijun, ZHAO Yanlin, LI Qingfeng, et al. Disaster mechanism of karst water bursting in mine[J]. Journal of China Coal Society, 2010, 35(3): 443-448.

[14] 赵延林, 彭青阳, 万文, 等. 高水压下岩体裂纹扩展的渗流-断裂耦合机制与数值实现[J]. 岩土力学, 2014, 35(2): 556-564.ZHAO Yanlin, PENG Qingyang, WAN Wen, et al. Seepagefracture coupling mechanism of rock masses cracking propagation under high hydraulic pressure and numerical verification[J]. Rock and Soil Mechanics, 2014, 35(2): 556-564.

[15] 赵延林, 吴启红, 王卫军, 等. 基于突变理论的采空区重叠顶板稳定性强度折减法及应用[J]. 岩石力学与工程学报, 2010,29(7): 1412-1513.ZHAO Yanlin, WU Qihong, WANG Wenjun, et al. Strength reduction method to study stability of goaf overlapping roof based on catastrophe theory[J]. Chinese Journal of Geotechnical Engineering, 2010, 29(7): 1412-1513.