基于Comsol 的矿井顶板水突出模拟研究

孙 凯，王 哲，吕 飞

（1.山东省煤炭泰山疗养院，山东 泰安271000；2.兖煤万福能源有限公司，山东 菏泽 274900）

0 引 言

矿井顶板突水问题一直是煤矿安全生产的重大隐患。采煤工作面顶板状况对于煤矿安全生产具有重大意义。传统的矿压理论有：预成裂隙梁、铰接岩梁、悬臂梁以及砌体梁模型等，运用这些理论研究采场顶板破坏也取得一定成果[1-3]。然而煤矿顶板破坏是一个极为复杂的过程，简单的理论分析不能精确的分析这一过程。而这些理论无法深入研究顶板水突出时，流体的分布情况。研究表明在出现顶板水突出时，工作面会有大量水涌[4-5]。运用数值模拟的方法对这一过程进行模拟研究对于研究顶板突出问题具有一定的意义。前人在这方面做过大量的研究[6-9]。

Comsol Multiphysics 是一款耦合模拟仿真软件，由于去具有开放界面，可结合MATLAB 对求解复杂多物理场模型，模拟出整个物理过程[10-11]。本文运用Comsol Multiphysic 中的PDE 模块进行模型镶嵌，结合MATLAB 编程，对复杂变量进行设定。模拟出矿井顶板突水过程中水分布情况已经整体变形。并将模拟结果和已做的实验做对比，二者的位移及渗流规律基本一致。

1 耦合模型控制方程

矿井水突出是一个复杂的动力学过程。突水过程中，岩层断裂。本文模型基于渗流—损伤力学。

1.1 应力场方程

在受力初始阶段，弹性多孔介质孔隙流体压力p及位移的本构方程为:

式中：G 为介质的剪切模量；v 为介质的泊松比；Fi和ui（i=x,y,z）分别为体力在位移i 方向的分量。

图1 单轴应力状态下的岩石的本构方程

设，，当介质的应力状态满足最大拉应力准则和摩尔库伦准则时，其分别发生拉伸损伤和剪切损伤：

1.2 渗流场方程

由流体质量守恒方程和Darcy 定律得：

式中：εv为体积应变，I 为介质的孔隙率；βl为孔隙流体的体积模量（Pa）；t 为时间（s）；k 为连续介质固有的渗透率（m2）；μl为流体的动力粘性系数（Pa.s），ρl为流体密度（kg/m2）；g 为重力加速度（m/s2）。

2 顶板突水过程数值模拟

以上数学模型在时间和空间域都是非线性的，因此运用有限元软件进行数值求解。本文应用运用MATLAB 强大的编程功能，联合Comsol Multiphysics对以上模型进行求解。

2.1 几何模型及参数

设计的几何模型尺寸为：160m（长）×190m（高），如图2 所示：

图2 几何模型

表1 模型岩石力学参数

表2 模型耦合参数

实验过程中整个岩层左边、右边及底面均受固定约束。对顶面进行加压，分析在压力变化过程中煤层顶板突水的情况。

2.2 数值模拟结果

模拟得出整个岩层加压至8MPa 时，开始出现顶板断裂，有大量水突出，整个模型流体速度场分布如图3、图4 所示。

图3 水速度场分布

图4 顶板水速度场曲线图

由图3 得，采空区出现顶板跨落时，工作面的水速度场分布集中，表明在此位置有大量水涌出。顶板水流速度分布如图4 所示。其模拟结果和前人的研究基本吻合。

图5、图6 分别为2MPa 和8MPa 外力作用下整个模型的应变。

图5 2MPa 下应变

图6 8MPa 下应变

由图5、图6 可以看出，在工作面推进到一定程度，直接暴露于采空区的上覆岩层不再受到下部煤层的支撑作用，在重力作用下发生弯曲、沉降、离层。煤层顶板覆岩变形破坏存在一个渐进的演变过程，是伴随工作面向前推进，采空区范围的不断扩大，应力的不断调整而发展变化的。当外界压力达到一定值时，出现顶板断裂和局部弯曲下沉。

图7 垂直位移分布

图8 顶板最终垮塌情况

图7 为顶板跨落时，整个模型竖直方向位移。由于受工作面回采的影响，工作面竖直方向位移变化较大，位移变化和理论研究基本一致。工作面前方岩层下沉梯度明显大于工作面后方岩层下沉梯度，岩层下沉呈非对称分布。

实际实验结果如图8 所示，模拟研究同实验研究结果基本吻合。

3 结 论

本文根据流固耦合相关理论，分析得出煤层顶板突水流固耦合模型。运用Comsol Multiphysics 结合MATLAB 进行数值模拟，得出顶板突水时，煤层、岩层模型体积应变、竖直方向上的位移及水速度场分布。模型的模拟变形同实际实验结果基本吻合，验证了数学模型的正确性。