港里铁矿水下开采流固耦合数值模拟与风险分析

马海涛 李迎佳 王云海

(1.北京科技大学土木与环境工程学院，北京100083;2.中国安全生产科学研究院，北京100012;3.鲁中矿业有限公司，山东莱芜，271113)

据不完全统计，我国目前受河流影响而滞留大量矿石储量的冶金矿山65座，有色金属矿25座，铁矿山最为突出，10亿t以上的矿石储量被压埋在河流之下未能合理开发，仅邯郸、莱芜地区就有5亿多t铁矿石赋存于水体之下，难以开发。因而，解决水体下安全开采问题成为冶金矿山所面临的重大生产安全问题之一［1］。水下开采地质灾害本质上是由水、岩等固体、准流体及流体构成的复杂地质系统对开挖扰动作出的响应或反馈。响应的方式和程度不同，灾害的类型和规模也就不同。从物理形态上硬岩、软岩、水等是可以分开的，但是，无论是天然条件下，还是开挖环境下，它们的变形、运动往往都互相联系、彼此影响。即应力场和渗流场之间存在耦合效应。因此，需要对水下矿山开采风险进行分析，并考虑流固耦合作用的影响。

有关流固耦合的计算方法大致可以分为两类，一类是解析解，一类是数值解。数值解法通常采用的是有限单元法或有限差分法［2-4］。有限差分法FLAC3D分析三维流－固耦合问题，应用较为广泛。

1 三维流固耦合计算方程

相互关联的变形－扩散过程的增量公式是通过线性准静态耦合理论的框架结构完成的。描述多孔介质中的渗流需要的变量包括孔隙压力和3个特殊的释放矢量。变量之间通过3个方程相互联系，3个方程分别是:①流体质量平衡方程;②渗流力学本构定律，体积应变对孔隙压力的影响是通过流体本构定律反映出来的，孔隙压力的变化引起了力学变形的发生;③渗流边界和初始条件，初始条件符合一个特定的压力场，边界条件通常是根据孔隙压力或者根据垂直于边界的比流量的分量来确定。

本研究数值计算考虑了4种情况:①给定孔隙压力;②给定垂直边界的比流量分量;③给定渗流边界;④给定不具有渗透性区域的边界。

2 工程概况

港里河在鲁中矿业有限公司港里1号矿体上部流经，为大汶河支流，属源短流急的间歇性河流，全长46 km，流域面积144 km2，历年最大洪峰流量为180.6 m3/s。

矿区内均有第四系地层覆盖。矿床主要含水层为中奥陶系大理岩、灰岩，含水性和透水性均较弱。矿体顶板为石炭系板岩、角岩，节理裂隙发育，性脆易碎，稳定性差;底板为中奥陶系大理岩，裂隙、溶洞一般不发育，较为稳定。当其含泥质较多时，易破碎坍塌。

根据矿山采用无底柱分段崩落法多年生产实践，本区地应力较大，接近矿体附近时巷道多变形或坍塌，不经支护几乎无法生产，回采难度较大，采准巷道多采用喷锚网支护。因此，矿床开采时第四系地下水和河水可能通过第三系黏土层中形成的导水裂隙带下渗，出现坍塌、突水危害。另一方面，矿体上部第三系为隔水层，如果采用充填法进行开采，限制围岩产生较大的变形，则不会破坏覆岩隔水层，可以保证矿山生产安全。

因此，本研究在矿山地质条件分析、上覆岩层隔水性分析基础上，进行三维流固耦合的数值模拟分析，评估河下矿体开采的安全性。

3 地表隔水层安全性分析

3.1 地质条件分析

港里矿地质条件简单，区内无较大断层、节理切割隔水层，地表水和地下水联系较少。故顶板裂隙带的发育高度，决定了隔水层的安全性。根据国家煤炭工业局制定颁发的《建筑物、水体、铁路及主要井巷煤柱留设与压煤开采规程》，正常地质构造和充分采动条件下，港里矿(中硬和软弱覆岩)裂隙带最大高度

式中，HL为冒落带加裂隙带的高度，m;M为开采矿体真厚度，m;h为开采垂高，m。

计算得港里矿开采后顶板覆岩裂隙带高度为105 m。由于港里矿采用充填法开采，在保证充填体强度和完全接顶的条件下，冒落带高度可以忽略，实际发生的裂隙带高度也远远小于该值。港里矿顶板上覆岩层为厚150～700 m的第三系红板岩，具有较好隔水性。因此，从地质条件上看，井下导水裂隙不会发育到地表，不会与地表水系产生直接联系，满足井下安全采矿的隔水要求。

3.2 岩石物理力学性质试验

通过现场岩石取样室内试验获得的岩石物理力学参数见表1。

表1 矿岩石力学参数Table1 Rock mechanics parameters

港里矿矿体上覆岩层变形量较大，且在采动影响下，其蠕变性明显，即岩体变形随着时间缓慢发展，且持续时间较长。研究结果表明，软岩中一般存在黏塑性区，且随时间而增大，但最终趋于稳定。通过实验室三轴试验结果表明，上覆软弱岩层——红板岩残余强度为2.0～3.0 MPa。当应力达到6 MPa时为稳定蠕变;当大于6 MPa时转化为加速蠕变。

采用最小二乘法、遗传算法和遗传规划法分别对蠕变本构模型进行参数辨识和模式识别。与红板岩蠕变试验结果进行比对，如图1。

图1 蠕变模型预测结果Fig.1 Forecasting result of themodel of creep behavior

如图1可见遗传规划预测曲线与实验曲线吻合程度较高，参数取值和蠕变模型较为合理。

3.3 地面河床安全性分析

采用三维概率积分法对地表移动变形进行预测，揭示开采沉陷的相应规律，从而判定开采对河床带来的影响。

选取大地坐标(55 500，17 000)为原点(0，0)建立模型，其中X为正轴指向东、Y为正轴指向北，预测分析在最大限度利用矿产资源的情况下地表移动变形对河流的影响，即在矿体开采结束后，开采沉陷影响范围。预测结果如图2。

图2 开采后地表下沉影响范围预测(单位:mm)Fig.2 Region of influence of ground depression after m ining(mm)

由图2经计算，地面河床最大下沉量在8 cm范围内。地表最大倾斜值小于3 mm/m，最大水平变形值小于2 mm/m。满足《采矿设计手册(矿床开采卷下)》和《建筑物、水体、铁路开采规程》规定的允许变形极值，即地表倾斜3 mm/m，地表曲率0.2×10－31/ m，地表水平变形2 mm/m。

4 三维流固耦合数值模拟

4.1 计算模型及边界条件

根据矿区地质特征，地表第四系取8 m，第三系隔水层取186 m，铁矿厚度取平均厚度14.35 m，倾角16°。1号矿体划分为2个阶段开采，阶段高度50 m，划分为－146 m和－196 m 2个中段开采，上向进路充填采矿法方案实际开采范围为－90～－190 m。

建立矿区三维模型，共计79 396个三维渗流计算单元，88 501个网格节点，如图3所示。

图3 数值模拟三维模型Fig.3 3D model for numerical simulation

根据矿区岩性属性分析，采用Mohr－Coulomb屈服准则判断材料的破坏［5］:

式中，σ1、σ3分别是最大和最小主应力;c、分别是黏结力和内摩擦角。当fs＞0时，材料将发生剪切破坏。在通常应力状态下，岩体的抗拉强度很低，因此可根据抗拉强度准则判断岩体是否产生拉破坏。计算中参数选取见表1。

4.2 计算结果及分析

从整体计算过程的不平衡力变化情况来看，模型在自重应力作用下逐渐趋于收敛，最大不平衡力变化率达到10－5以下。说明模型最终稳定，计算数据可信［6］。

4.2.1 无充填开采模拟结果

受到矿体开采扰动的影响，地表以及顶板岩层在不同程度位移作用下，会进入塑性状态，破坏形式主要以拉伸破坏为主。破坏程度较大。

在不采用充填采矿法的情况下，对上向进路采矿进行数值模拟，底板最大下沉17 m，地表下沉7 m。

如果不采用充填法开采，地表会发生很大的塌陷，裂隙带直接与地表水沟通，造成严重的安全事故。因此，港里矿开采必须采用充填法。

4.2.2 充填开采模拟结果

(1)水平地应力场分析。水平地应力的分布服从由上到下逐渐增加的变化趋势，模型底部水平地应力达到极值。随着开挖的进行，受到开采扰动的影响，磁铁矿内部应力卸载，形成应力降低区，围岩有应力集中现象。

(2)破坏场分析。受到矿体开采扰动的影响，地表在不同程度位移作用下，会进入塑性状态，破坏形式主要以拉伸破坏为主。地表稍微有破坏，如图4所示。尽管合理使用进路充填采矿法可以有效控制地应力，在应力集中和应力升高区，可以通过适当提高充填体强度或开采时间调整等来缓解地压的影响，但仍需采取河道治理的措施。

图4 围岩破坏区分布Fig.4 Distribution map of wall rock failure zone

图5 垂直位移分布Fig.5 Distribution map of vertical displacement

(3)位移场分析。围岩的位移方向主要指向采场，由于采用上向水平充填采矿法，铁矿开采后立刻回填，围岩产生的位移量较小。矿区沉降位移不明显，受开采扰动影响，仍会形成下沉盆地(地表移动盆地)。地表整体沉降量不大，矿区稍微有些沉降，在4～8 cm范围。矿区水平位移最大值也主要发生在顶板岩层，最大下沉量15 cm，如图5所示。三维数值计算结果与开采地表沉陷理论分析相符，说明预测结果数据可信。

(4)渗流场分析。孔隙水压力主要集中在第四系底层，极值约为0.1 MPa，矿体围岩孔隙水压力很小。渗流矢量最大值为1.962×10－2m/s。但随着充填开采的进行，孔隙水压力和渗流矢量在不断增加。随着开挖的进行，受到开采扰动的影响，第四系的水会向下流。同时，顶板的水会流向开采区域，成层状分布。

5 结论

(1)由于受开采扰动的影响，地表仍然会发生移动变形，有一定下沉量，但矿区整体沉降量不大，在4～8 cm范围。满足采矿扰动引起的地表最终下沉情况和《建筑物、水体、铁路及主要井巷煤柱留设与压煤开采规程》等相应的安全规定。

(2)矿井不受地表水害威胁，可以满足河下安全开采的要求。围岩产生的位移量较小，矿体围岩孔隙水压力很小，随着充填开采的进行，孔隙水压力和渗流矢量在不断增加。

(3)根据三维渗流数值模拟和概率积分法计算分析，港里铁矿1号矿体的开采，造成地表平缓、均匀的下沉和水平移动，总体下沉量不大，在港里河下开采，不会造成地表河流堤坝变形破坏，不会影响河流储水能力。

(4)采用三维流固耦合(渗流场与力场相互作用)力学数值模拟分析研究港里铁矿1号矿体采用上向进路充填采矿法开采过程中的渗流场变化和力学稳定性，为矿山设计和生产提供了参考依据。

［1］ 何国清，杨 伦，贾凤彩.矿山开采沉降学［M］.北京:中国矿业大学出版社.1991.

He Guoqing，Yang Lun，Jia Fengcai.Mining Subsidence Science［M］.Beijing:China University of Mining and Technology Press，1991.

［2］ 朱 军，刘光廷，蓝青松，等.饱和模型非定常渗流的数值求解方法［J］.岩土力学，2002，23(3):355-357.

Zhu Jun，Liu Guangting，Lan Qingsong，et al.Numericalmethod for saturated unsteady seepage［J］.Rock and Soil Mechanics，2002，23 (3):355-357.

［3］ 王金安，马海涛.河下金属矿合理开采顺序的固流耦合力学分析［J］.金属矿山，2005(S):133-134.

Wang Jinan，Ma Haitao.Numerical analysis of the scheme formining metal ore under the body of water［J］.Metal Mine，2005(S):133-134.

［4］ 马海涛.“11.6”特别重大坍塌事故矿区采场稳定性三维数值模拟分析［J］.中国安全生产科学技术，2007，3(6):68-72.

Ma Haitao.3D numerical simulation of stope stability in 11.6 accident of gypsum mine collapse［J］.Journal of Safety Science and Technology，2007，3(6):68-72.

［5］ 蔡美峰，何满潮，刘东燕.岩石力学与工程［M］.北京:科学出版社，2002.

CaiMeifeng，He Manchao，Liu Dongyan.Rock Mechanics and Engineering［M］.Beijing:Science Press，2002.

［6］ 王金安.岩土工程数值计算方法实用教程［M］.北京:科学出版社，2010.

Wang Jinan.Practical Tutorial of Numerical Calculation Method on Geotechnical Engineering［M］.Beijing:Science Press，2010.