防水煤柱安全评价的流-固耦合分析

战玉宝，刘焕芝，尤春安

( 1.山东科技大学资源与土木工程系，山东 泰安 271019；2.山东科技大学泰安校区，山东 泰安 271019 3. 山东科技大学土建学院，山东 青岛 266510)

突水是煤矿防水煤柱失效的重要形式，引起的损失是巨大的。突水事故轻者会冲垮工作面、淹没设备、增加矿井的排水负担，严重者会造成人员伤亡、淹井等。全国已经发生许多起近松散层开采上部含水层造成的水砂突涌和淹井事故，造成重大经济损失。防水煤柱失稳突水，不仅涉及到煤柱在水压力作用下的力学性质与长期稳定性问题，还涉及到渗流场与应力场的耦合问题。[1]

含水松散层下采煤，是水体下采煤的重要方面。近年来，华东、华北、西北等地的许多煤矿，开采上限不断提高，己经在不同矿区、不同富水程度的松散含水层下开展近松散层开采。对防水煤柱的合理留设进行了深入研究，[2-4]并取得一批重要理论成果。于进广等[5]对中等含水层下留设防砂煤柱的安全开采机理进行了研究，廖学东、疏开生[6]对松散含水层下采煤合理煤岩柱高度的确定进行了探讨，蔡荣[7]对水体下采煤岩层的移动与保护层理论及应用进行了研究，尹先尚[8]对防水煤柱留设对围岩变形的影响进行了数值模拟方面的研究。中国矿业大学的狄乾生[9]和隋旺华[10]，结合“三下”采煤和开采覆岩移动的有限元分析，对工程地质模型进行了深入研究，建立了厚松散含水层下提高开采上限的工程地质模型。

岩土体地下渗流的流-固耦合问题的研究已经开展多年，但流-固耦合理论应用在防水煤柱安全评价方面的研究非常少，这方面的文献能检索到的几乎没有。本文就是应用有限差分软件FLAC3D，对某矿煤矿四采二层煤防水安全煤柱进行渗流场与应力场的耦合分析。

1 计算分析方法与原理[11-13]

FLAC3D模拟多孔介质(如煤岩体)中流体流动时，流体的模拟独立于结构计算。其主要通过孔隙水压力的消散，引起岩体中位移的变化，这一过程包含两种力学效果。第一、孔隙水压力的变化，引起结构体中有效应力的变化；第二、孔隙水压力的变化，又引起流体区域的变化。流体在孔隙介质中的流动依据Darcy定律，流-固耦合过程满足Biot方程。

FLAC3D中，实现流-固耦合分析的几个关键方程为[13]：

(1)平衡方程

对于小变形，流体质点平衡方程为：

(1)

式中：qi是渗流速度(m/s)；qv是被测体积的流体源强度(1/s)；ζ是单位体积孔隙介质的流体体积变化量。而

(2)

式中：M是Biot模量(N/m2)；p是孔隙压力；α是Biot系数；ε是体积应变；T是温度；β是考虑流体和颗粒热膨胀系数(1/°C)。

(2)运动方程

流体的运动用Darcy定律来描述。对于均质、各向同性固体和流体密度是常数的情况，这个方程具有如下形式：

(3)

式中：k是介质的渗透系数(m2/pa·s)；ρf是流体密度(kg/m3)；gi(i=1，3)是重力加速度的三个分量(m/s2)。

(3)本构方程

流体的流动导致孔隙介质中，孔隙压力p、饱和度s、体积应变e和温度T的改变，孔隙流体方程为

(4)

体积应变的改变，引起流体孔隙压力的变化，反过来，孔隙压力的变化也会导致体积应变的发生。

孔隙介质本构方程的增量形式为

(5)

(4)相容方程

应变率和速度梯度之间的关系为

(6)

式中：υi是介质中某点的速度。

(5)边界条件

在计算中，有四种类型的边界条件，他们分别是：①给定孔隙水压力；②给定边界外法线方向流速分量；③透水边界；④不透水边界。

不透水边界程序中默认，透水边界采用下式给出：

(7)

式中：qn是边界外法线方向流速分量；h是渗漏系数(m3/N·s)；p是边界面处的孔隙水压力；pe是渗流出口处的孔隙水压力。

2 数值计算方案

2.1 工程概况

山东某矿矿井田西翼四采二层煤煤层露头区，位于该矿井四采风井以西，煤层露头东西走向长450m，走向92°，平均倾角25°。露头区地面位于一条河的北岸，东部为该河所覆盖。由于附近砖厂取土，村民挖砂，致使煤层露头上覆的第四系表土覆盖层受到不同程度的破坏。加上煤层浅部被非法开采、乱采和盗采，近年来四采区二层煤露头区域地表明显地下沉，形成连续的斑裂区，与地表形成直接的水力联系。斑裂区位于该河北岸，是雨季洪水淹没范围，因此存在地表水泄漏威胁。

2.2 计算模型与参数

数值模拟时，模型取为340m×160m×1m，水平面内向右为X轴，垂直向上为Y轴，X、Y、Z轴成右手系。计算范围取为-50m≤X≤290m，0m≤Y≤160m，-1m≤Z≤0m。对模型进行划分，总共得到2106个单元，4388个节点。

采用摩尔-库仑材料模拟岩石与煤层，计算参数根据该矿提供的岩石的物理力学性质选定。具体的参数见表1。

2.3 初始地应力场

根据煤矿提供的资料，初始地应力在垂直方向按岩体自重考虑，水平地应力根据勘测单位提供的侧压力系数为0.6考虑。

3 数值计算结果分析

3.1 不考虑流-固耦合时计算结果分析

此种情况下，只考虑采空区顶板以及防水煤柱受自重和上覆岩层的作用，分别对20m、24m、28m、32m不同长度的防水煤柱进行数值模拟，得到了开挖后采空区顶板的塑性区分布以及防水煤柱的塑性区分布(图1)。

图1 未考虑流-固耦合时不同长度的安全防水煤柱的塑性区

图2 流-固耦合时不同长度的安全防水煤柱的塑性区分布

(1)采空区顶板上端20m范围内塑性区较大，有向上发展的趋势；在采空区下端30m范围出现较大塑性区；中间塑性区较小，比较平缓；整个塑性区呈现为两端高、中间低的“马鞍”形状。

(2)当煤柱大于24m时，煤柱上方采空区下端的塑性区比留设20m煤柱时的塑性区小，随着煤柱长度的增加，塑性区分布开始趋于稳定。

(3)安全煤柱长度从20m一直到32m的情况下，整个煤柱均未全部进入塑性状态。

3.2 考虑流-固耦合的计算结果分析

当上部采空区与地表形成水力连接时，考虑采空区顶板以及安全煤柱受自重和上覆岩层的作用，同时考虑水压力的作用。分别对20m、24m、28m、32m不同长度的防水煤柱进行数值模拟，得到了开挖后采空区顶板的塑性区分布以及防水煤柱的塑性分布(图2)

分析图2可以得出：

(1)上部采空区顶板上端20m范围内，塑性区有继续向上发展的趋势，在采空区下端30m范围，也出现更大的塑性区，中间塑性区较小，比较平缓；整个塑性区呈现为两端高中间低的“马鞍”形状。

(2)当煤柱留设20m和24m时，煤柱全部进入塑性状态，并且水能够透过煤柱及其上方的裂隙带，形成水力通道，会出现水溃入下部采空区的现象，对安全生产构成威胁。

(3)当煤柱留设28m、32m时，煤柱未全部进入塑性状态，虽然非塑性区范围不大，但是可以通过采取治理措施加以控制，保证生产安全。

(4)考虑流-固耦合后，底板出现了不同程度的塑性区。

4 结论

对比不考虑流-固耦合与考虑流-固耦合作用两种情况下防水煤柱的数值分析结果，可以得出如下结论：

(1)考虑流-固耦合时与不考虑耦合时，塑性区都呈现为两端高、中间低的“马鞍”形状，但考虑耦合时塑性区的范围比较大。

(2)考虑流-固耦合作用时，能很好的模拟底板的受力状态与变形破坏规律，与不考虑流-固耦合相比，煤柱留设需加长。

(3)采用FLAC3D，可以很好地对防水煤柱进行流-固耦合数值模拟，模拟结果与实际情况更加接近。

[1] 张金才, 张玉卓, 刘天泉. 岩体渗流与煤层底板突水[M]. 北京: 地质出版社, 2001.

[2] 许传峰, 王思栋, 纪宪江，等. 柴里煤矿留12米煤柱的开采上限研究[J] .河北建筑科技学院学报,2001, 18(2):68-70.

[3] 隋旺华. 煤层开采防水煤岩柱工程地质预计[J]. 江苏煤炭, 1994,(2):21-25.

[4] 马福祥, 隋旺华. 山东横河煤矿分层开采防水煤柱留设[J].中国地质灾害与防治学报, 1994, 5(增):15-17.

[5] 于进广, 陈风祥,等. 中等含水层下留设防砂煤柱的安全开采机理[J].淮南工业学院学报, 2001, (2): 15-17.

[6] 廖学东, 疏开生. 松散含水层下采煤合理煤岩柱高度的确定[J] .淮南矿业学院学报, 1998.8(2):29-33.

[7] 蔡 荣. 水体下采煤岩层的移动与保护层理论及应用[D].徐州: 中国矿业大学, 2003.

[8] 尹先尚. 陷落柱防水煤柱留设对围岩变形影响的数值模拟[J].煤炭学报, 2006, 4(2):179-182.

[9] 狄乾生，隋旺华，黄山民. 开采岩层移动工程地质研究[M]. 北京：中国建筑工业出版社, 1992.

[10] 隋旺华, 董青红,等. 工程地质模型在防水煤岩柱研究中的应用[J].中国矿业大学学报, 1999，28(4):417-420.

[11] FLAC-3D,Fast Lagrangian Analysis of Continua in 3Dimensions, Version 3.0,User’s Manual, Itasca Consulting Group, Inc. USA.

[12] 刘 波，韩彦辉. FLAC原理、实例与应用指南[M] . 北京：人民交通出版社, 2005.

[13] 李廷春，李术才，陈卫忠，等. 厦门海底隧道的流固耦合分析[J] .岩土工程学报, 2004, 26(3):397-401.