基于COMSOL的顶板运移规律及渗流特性研究

许大洋，李晓泉，王兄威，李明菊

(广西大学资源与冶金学院， 广西南宁530004)



基于COMSOL的顶板运移规律及渗流特性研究

许大洋，李晓泉，王兄威，李明菊

(广西大学资源与冶金学院， 广西南宁530004)

为探讨水体下煤层开采复杂条件下顶板的运移规律及渗流特性，以流固耦合理论为基础，利用COMSOL Multiphysics多物理场耦合分析软件建立了水体下开采的应力场—渗流场耦合数值模型，针对不同开采条件下顶板岩体的运移规律及渗流特性进行定性研究。数值模拟结果表明：工作面推进长度的增大及水压的增大均会导致顶板岩体变形程度的增加，进而导致顶板岩体的渗流速度相应增大。因此，水体下开采时十分有必要对水压情况进行实时监测和预警，并综合考虑工作面的安全开采长度等因素，以防止导水裂隙带发展到水体，从而避免工作面与水体间形成导水通道而造成透水事故。

采动覆岩；数值模拟；流固耦合；渗流特性

煤炭支持着我国经济的快速增长和社会的全面进步，是我国能源的基石。据不完全统计，仅我国煤矿生产矿井“三下”压煤量就达140×108t左右，其中水体下(包括承压废岩水上)压煤占28%左右[1]。因此，水体下的煤炭资源开采对我国煤炭科技的进步和经济建设具有十分重要的意义。水体下开采一般包括地面水体下和地下水体下的开采。地面水体包括江河湖海、水库池塘、山沟小溪以及地表积水等；地下水体包括表土层的砂层水、砂岩含水层、岩溶水以及老窟水等[2]。对于赋存于水体下的煤层，开采扰动使上覆岩层发生变形和破坏，形成的裂隙将有可能波及至水体，使工作面与水体之间形成导水裂隙，从而对工作面的安全造成极大威胁。

迄今为止，国内外学者围绕采动后岩体的运移规律、渗流特性及流固耦合理论进行了大量研究，取得了许多研究成果和经验[3-13]。其中,早期针对岩体渗流特性的研究大多以理论研究为主，用解析方法往往很难求出其精确解，因而不能很好地解决实际工程问题。随着现代计算机技术的不断提高以及有限元法相关理论的不断发展与成熟，有限元法逐渐成为解决这类问题的最主要也是最有效的途径。目前，针对底板承压含水层及岩溶水流固耦合问题的研究方法相对成熟，而对顶板含水层流固耦合问题的研究成果相对较少，这主要是由于裂隙带范围难以控制、裂隙带发育高度难以预测，以及裂隙带与采场、水源之间的空间关系难以判断等因素造成的。理论研究表明，矿山采动围岩破坏透水过程中，水流经历了在含水层及弯曲带中的Darcy流、破碎带中的非Darcy快速流以及巷道中的Navier-Stokes紊流这3个物理过程[14]。其中，Darcy方程以流体压力驱动为主，适用于低渗透的多孔介质材料，地下水在低渗透多孔介质中的渗流可采用该方法计算[15-16]。因此，本研究主要针对强含水层下压煤的开采，利用COMSOL Multiphysics软件建立了采动覆岩渗流的二维流固耦合数值分析模型，研究强含水层压煤开采条件下顶板岩层的破坏与地下水在低渗透多孔介质中的渗流—应力耦合作用。

1 多孔介质的流固耦合机理

多孔介质由两部分构成，分别为由固体物质组成的骨架以及由骨架分隔开而形成的密集成群的大量微小空隙。流固耦合的本质是通过流体的本构方程来反映材料的体应变对孔隙压力的影响，而孔隙压力的变化又能反过来引起多孔介质产生变形，从而影响其力学性质。由多孔介质骨架和流体所组成的耦合系统在受到外力作用时其骨架将发生变形，流体也要发生流动状态的改变。其中，骨架的变形是在外部荷载与缝隙流体的共同作用下形成的，即控制骨架变形的是流固耦合作用下的有效应力。

有效应力的概念及原理最早由Terzaghi于1923年研究饱和土的固结及水与土壤的相互作用时提出，它反映了孔隙流体压力与骨架结构变形之间的关系[17]。Biot在建立多孔介质本构模型时，在Terzaghi有效应力原理的基础上，将孔隙压力与水容量变化联系起来，同时考虑水的不可压缩，认为水容量的变化与多孔介质孔隙率增量相等，给出了如下一般意义上的等效连续介质耦合关系式：

①平衡方程：

(1)

②几何方程:

εv=εx+εy+εz,

(2)

(3)

③本构方程:

σij=(λεv+αp)δij+2Gεij;

(4)

④渗流方程:

(5)

式中，ρ为多孔介质整体密度；Xj为体力密度；ε、u分别为应变、应力；λ、G分别为拉梅常数和剪切模量；Kij为渗透系数张量；Q为多孔介质总体积保持不变的情况下受到孔隙水压作用而被挤进多孔介质中的水量。

式(1)～(5)的联合便实现了流固耦合。

2 数值模型的建立

2.1 模型基本假设

数值模拟过程中，由于地层组成、地质构造运动等的局限性，不可能完全反映出实际情况。因此，通常需要对实际工程问题进行合理简化，使数值模拟过程既简单又尽可能反映实际情况。本次数值模拟基于以下假设条件：

①不考虑相邻岩层之间的界面影响，且把同一岩层看成是均匀、各向同性介质；

②忽略地质构造运动产生的构造应力；

③忽略开挖前岩体中存在的节理和裂隙，同时考虑采动煤层顶板覆岩的岩移三带分布特征，选取岩层参数时进行适当的弱化处理；

④地下水在多孔介质中的运动符合Darcy定律。

2.2 计算模型及边界条件

图1 计算模型Fig.1 Computational model

本数值研究计算中，模型长度350 m，模拟开挖长度180 m，左右各预留85 m以消除边界影响；近水平煤层厚度4 m，采全厚，模拟采深260 m；模拟工作面回采共分6个开挖步，每步开挖30 m；距煤层顶部110 m处有一承压含水层，水压均匀分布；开采模型顶部取至承压含水层，最终高度为200 m。利用COMSOL Multiphysics中的自由三角形网格生成器剖分得到5 407个单元，其中采空区用材料开关控制，计算模型网格如图1所示。

力学边界条件：模型底部限制垂直移动；模型左右两侧水平方向施加由自重应力产生的侧向应力(σx=0.45σy)，限制水平移动；模型顶部按岩体自重施加等效均布荷载(σy=γh，其中，h为模型顶部距地表的距离，γ为上覆岩体的容重，本研究模拟取平均值为22 kN/m3)。渗流边界条件：模型左右、底部边界及开挖边界给定不透水边界；隔水层初始孔隙水压力为零；承压含水层水压进行梯度变化赋值。本构模型采用D-P模型，匹配摩尔—库仑屈服准则。由于本数值模拟计算中膨胀角不会产生多大影响，在计算时取膨胀角等于0。数值模型岩石力学参数见表1。

表1 岩石力学参数

3 工作面不同推进长度条件下顶板的运移规律及渗流特性

以含水层水压为2 MPa时为例，模拟工作面不同推进长度条件下(30、60、90、120、150、180 m)顶板岩层的运移规律及渗流特性。

3.1 不同推进长度对y方向位移的影响

图2所示为工作面不同推进长度条件下顶板岩体的y方向位移云图。从图2中可以看出，不同推进长度条件下，位于采空区正上方的岩体y方向位移均较大，走向两侧时均逐渐减小；非充分采动时，推进长度越大，局部沉降漏斗范围越大，覆岩最大下沉值也不断增加；达到充分采动时，随着工作面的继续推进，局部沉降漏斗范围继续扩大，而覆岩最大下沉值不再增加。

单位：m

(a) 推进长度为30 m

(b) 推进长度为60 m

(c) 推进长度为90 m

(d) 推进长度为120 m

(e) 推进长度为150 m

(f) 推进长度为180 m

图2 工作面不同推进长度时覆岩y位移云图

Fig.2 Displacement ofydirection at different working face position

选取模型上边界(y=0处)为研究对象，当工作面分别推进30、60、90、120、150、180 m时，模型上边界的最大沉降量分别为0.262、0.329、0.374、0.408、0.456、0.459 m，拟合曲线如图3所示。可以看出，随着工作面不断推进，模型上边界的最大沉降量也逐渐增加，增加的速率呈现逐渐减小的趋势。这表明达到充分采动后，模型上边界的最大下沉值将不再增加。

图3 模型上边界最大沉降量变化

3.2 不同推进长度对应力的影响

图4所示为工作面不同推进长度条件下顶板岩体的Von mises应力云图。从图4中可以看出，开挖后的力学响应使得初始地应力场发生重分布，在采空区两端具有明显的应力集中现象。

单位：Pa

(a) 推进长度为30 m

(b) 推进长度为60 m

(c) 推进长度为90 m

(d) 推进长度为120 m

(e) 推进长度为150 m

(f) 推进长度为180 m

图4 工作面不同推进长度时应力云图

Fig.4 Stress at different working face position

图5 应力最大值变化Fig.5 Change of the maximum stress

推进长度分别为30、60、90、120、150、180 m时，采空区两端的最大Von mises等效应力分别为12.6、19.1、22.8、23.2、25.9、28.0 MPa，拟合曲线如图5所示。可以看出，采空区两端的最大Von mises等效应力呈现随推进长度增加而增大的趋势。

3.3 不同推进长度对塑性区高度的影响

图6所示为工作面不同推进长度条件下的塑性区云图。可以看出，随着工作面的不断推进，塑性区高度不断增加。具体表现为：当工作面分别推进30、60、90、120 m时，顶板塑性区较小，高度分别为12.76、28.68、49.73、73.65 m，尚未发展到承压含水层；当工作面推进到150m时，煤层顶板塑性区高度增至116.77 m，塑性区刚好发展到含水层；当工作面继续推进到180 m时，塑性区在含水层内部向上及两侧发展，高度为126.81 m。

(a) 推进长度为30 m

(b) 推进长度为60 m

(c) 推进长度为90 m

(d) 推进长度为120 m

(e) 推进长度为150 m

(f) 推进长度为180 m

图6 工作面不同推进长度时塑性区云图

Fig.6 Equivalent plastic strain zone at different working face position

图7所示为上覆岩层塑性区高度拟合曲线。通过拟合分析得出塑性区高度与推进长度之间的关系式为：

y=0.001 7x2+0.572 23x-7.964，

(6)

式中，y为塑性区高度(m)；x为推进长度(m)。趋势线拟合度为0.964 25，说明拟合效果良好。

3.4 不同推进长度对隔水层界面渗流速度的影响

由“1节”的分析可知，应力场的改变将会导致渗流场发生变化，而渗流场的改变又能反过来影响多孔介质的变形。因此，有必要对应力场—渗流场耦合作用下的渗流特性进行研究。图8所示为工作面不同推进长度条件下隔水层界面(y=-50 m处)的最大渗流速度变化情况。从图8中可以看出，随着工作面的不断推进，隔水层界面的最大渗流速度也不断增大，且增大的速率也逐渐增加。当工作面分别推进30、60、90、120、150、180 m时，隔水层最大渗流速度分别为1.201×10-10、1.438×10-10、1.907×10-10、2.506×10-10、3.257×10-10、4.070×10-10m/s。

图7 塑性区高度变化

Fig.7 Change of the height of equivalent plastic strain zone

图8 最大渗流速度变化

Fig.8 Change of the maximum seepage velocity

通过拟合分析得出隔水层界面最大渗流速度与工作面推进长度之间的关系式为：

y=0.000 08x2+0.002 73x+1.024 4，

(7)

式中，y为最大渗流速度(10-10m/s)；x为工作面推进长度(m)。趋势线拟合度为0.999 01，说明拟合效果良好。

由分析可知，当工作面推进150 m后，隔水层产生的塑性区将导通工作面与含水层，形成宏观的导水通道。而此时顶板的变形破坏使得隔水层具有较大的渗流速度，因而极可能引发透水事故，造成群死群伤的严重后果。因此，必须采取有效防范措施(如采用预留煤柱或充填法等)，通过降低顶板的变形及破坏程度，防止工作面与水体之间形成直接的导水通道，从而实现水体下煤层的安全开采。

4 含水层不同水压作用下顶板的运移规律及渗流特性

以工作面推进长度为150 m时为例，研究含水层不同水压(0、0.5、1.0、2.0、3.0、4.0 MPa)条件下顶板岩层的运移规律及渗流特性。

4.1 不同水压对y方向位移的影响

由“3.1节”分析可知，位于采空区正上方的岩体具有较大的沉降变形。图9所示为工作面推进150 m时，模型上边界(y=0处)在承压含水层不同水压条件下的最大沉降量。可以看出，随着承压含水层水压的增加，模型上边界的最大沉降量也不断增加，增加的速率呈现逐渐减小的趋势。

4.2 不同水压对应力的影响

由“3.2节”分析可知，开挖后造成的应力集中现象主要发生于采空区的两端。在承压含水层不同水压条件下，采空区两端Von mises等效应力的最大值变化如图10所示。由图10可以看出，应力最大值随着承压含水层水压的增大而增大，二者具有较好的线性关系。

图9 模型上边界最大沉降量变化

Fig.9 Change of the maximum settlement

at the upper boundary

图10 应力最大值变化

Fig.10 Change of the maximum stress

4.3 不同水压对塑性区高度的影响

图11所示为工作面推进150 m时含水层不同水压条件下的覆岩塑性区高度变化。从图11中可以看出，塑性区高度随着水压的不断增加而增加。具体表现为：当水压分别为0、0.5、1.0 MPa时，顶板塑性区高度分别为96.72、103.34、106.69 m，塑性区尚未发展至承压含水层；当水压增至2 MPa时，塑性区发展到含水层中6.77 m，此时含水层下部的岩层塑性区仍未发展至承压含水层；当水压继续增加至3.0、4.0 MPa时，塑性区在含水层内部向上及两侧发展，塑性区高度分别为123.30、126.59 m，此时含水层内部的塑性区与底部隔水层向上发展的塑性区导通，工作面与水体间形成直接的导水通道。

通过拟合分析得出塑性区高度与承压含水层水压之间的关系式为：

y=-1.194x2+12.299x+96.733，

(8)

式中，y为塑性区高度(m)；x为含水层水压(MPa)。趋势线拟合度为0.994 84，说明拟合效果良好。

4.4 不同水压对隔水层界面渗流速度的影响

图12所示为承压含水层不同水压条件下隔水层界面(y=-50 m处)的最大渗流速度变化情况。从图12中可以看出，随着承压含水层水压的增加，隔水层界面的最大渗流速度也不断增大。当水压为0～2 MPa时，塑性区尚未波及至含水层或未形成连续的塑性区，隔水层仍具有一定的隔水性能，此时隔水层界面最大渗流速度随着承压含水层的增加呈现缓慢增加的趋势；当水压增至3 MPa以上时，隔水层界面最大渗流速度随着承压含水层的增加也增加，但由于含水层内部的塑性区与底部隔水层向上发展的塑性区导通，形成了直接的导水裂隙，因而渗流速度增加速率较水压为0～2 MPa时有明显提升。

通过拟合分析得出隔水层界面最大渗流速度与承压含水层水压之间的关系式为：

y=0.053 7x2+0.194 3x+2.724 4，

(9)

式中，y为最大渗流速度(10-10m/s)；x为含水层水压(MPa)。趋势线拟合度为0.982 85，说明拟合效果良好。

由分析可知，含水层水压的增大同样会导致上覆岩层的变形、破坏程度和渗流速度的增大。因此，在水体下煤层的开采中必须对含水层水压及渗流参数进行实时监测和预警，并综合考虑工作面推进长度等因素所造成的顶板变形和破坏，以防止工作面与水体之间形成直接的导水裂隙，从而实现水体下煤层的安全开采。

图11 塑性区高度变化

Fig.11 Change of the height of equivalent

plastic strain zone

图12 最大渗流速度变化

Fig.12 Change of the maximum seepage velocity

5 结 语

本文以流固耦合理论为基础，针对水体下煤层开采的复杂条件，利用COMSOL Multiphysics多物理场耦合分析软件建立了水体下开采的应力场—渗流场耦合数值模型，对采动影响下煤层顶板的运移规律及渗流特性进行了研究，得出如下主要结论：

①随着工作面的不断推进，顶板岩层的变形程度、应力、塑性区高度及隔水层界面渗流速度也逐渐增大；当工作面推进到一定长度后，顶板岩层塑性区可能波及至承压含水层，此时工作面与水体之间将形成连续的塑性区，使得渗流速度大大增加，从而给工作面的安全开采带来威胁。

②随着顶板覆岩中含水层水压的增大，顶板岩层的变形程度、应力、塑性区高度及隔水层界面渗流速度也逐渐增大；当水压增加到一定值后，隔水层中产生的塑性区将导通工作面与水体，使二者之间形成直接的导水裂隙，从而使得渗流速度大大增加，增加工作面透水事故发生的概率。

③水体下开采时有必要对承压含水层的水压情况进行实时监测和预警，并结合实际工程问题，综合考虑工作面的安全开采长度等因素，通过合理预留煤柱或充填采空区等方法减小煤层开采对上覆岩层造成的破坏，保证导水裂隙带不会波及至水体，从而实现水体下煤层的安全开采。

本文的研究结果能为矿区顶板透水灾害的防御和治理提供一些参考。在实际工程问题中，由于岩体的复杂性，需要进一步结合现场实际情况定量研究采动影响下顶板岩层的渗流速度及渗水量，以便更好地为矿区的安全防御提供理论指导。

[1] 邓宁.“三下”采煤技术发展现状及前景展望[J]. 中州煤炭，2011(12)：56-58.

[2] 朱红飞，冉德立.建筑物下、铁路下和水体下采煤技术分析[J]. 内蒙古煤炭经济，2013(9)：126-133.

[3] 冉启全，顾小芸.考虑流变特性的流固耦合地面沉降计算模型[J]. 中国地质灾害与防治学报，1998，9(2)：99-103.

[4] 任春辉，李文平，李忠凯，等.巨厚岩层下煤层顶板水突水机理及防治技术[J]. 煤炭科学技术，2008，36(5)：46-48.

[5] 李利平，李术才，李树忱，等.松散承压含水层下采煤的流固耦合模型试验与数值分析研究[J]. 岩土工程学报，2013，35(4)：679-690.

[6] 陈占清，李顺才，浦海，等.采动岩体蠕变与渗流耦合动力学[M]. 北京：科学出版社，2010：208-225.

[7] 张金才，张玉卓，刘天泉.岩体渗流与煤层底板突水[M]. 北京：地质出版社，1997：56-75.

[8] 李学山，田微微，李芒原.渗流对深基坑工程的影响研究[J]. 广西大学学报(自然科学版)，2008，33(S1)：192-194.

[9] 苏国韶，燕柳斌，李海，等.高地应力与强渗透水压下软弱围岩支护[J]. 广西大学学报(自然科学版)，2009，34(5)：613-617.

[10]彭世龙，荣传新，程桦.潘三煤矿西风井井壁突水机理分析[J]. 广西大学学报(自然科学版)，2015，40(4)：1038-1043.

[11]CAI M，HORRI H.A constitutive model of highly jointed rock masses[J]. Mech Mater，1992，13(3)：217-246.

[12]MARTINS R.Turbulent seepage flow through rockfill structures[J]. International Water Power and Dam Construction，1990，42(3)：41-45.

[13]KIRKHAM D.Explanation of paradoxes in Dupuit-Forchheimer seepage theory[J]. Water Resources Research，1967，3(2)：609-622.

[14]杨天鸿，陈仕阔，朱万成，等.矿井岩体破坏突水机制及非线性渗流模型初探[J]. 岩石力学与工程学报，2008，27(7)：1411-1416.

[15]杨天鸿，唐春安，谭志宏，等.岩体破坏突水模型研究现状及突水预测预报研究发展趋势[J]. 岩石力学与工程学报，2007，26(2)：268-277.

[16]YUAN S C，HARRISON J P.Numerical modeling of progressive damage and associated fluid flow using a hydro-mechanical local degradation approach[J]. International Journal of Rock Mechanics and Mining Sciences，2004，41(Sup.1)：417-418.

[17]缪协兴，刘卫群，陈占清.采动岩体渗流理论[M]. 北京：科学出版社，2004：323-327.

(责任编辑 张晓云 裴润梅)

Study on migration rules and seepage characteristics of roof based on COMSOL

XU Da-yang， LI Xiao-quan， WANG Xiong-wei， LI Ming-ju

(College of Resources and Metallurgy， Guangxi University，Nanning 530004，China)

In order to explore the migration rules and seepage characteristics of coal roof when mining under water body with complex conditions, basing on the fluid-structure interaction theory and using COMSOL Multiphysics more physical field coupling analysis software, it has built the stress-seepage coupling numerical model of mining under water body to conduct on the qualitative research of migration rules and seepage characteristics of coal roof under different mining conditions.The simulation results show that: the increase of the length of working face advancing and the hydraulic pressure will lead to the increase of deformation degree of coal roof, which leads to the seepage velocity of coal roof increased correspondingly.Therefore, it’s very necessary to real-time monitor and warn the hydraulic pressure condition when mining under water body, and consider comprehensively the safety mining length of working face and other factors to prevent the water transmitting fissure zone from developing into the water and the formation of water channel between working face and water body, so as to avoid permeable accident.

mining overburden rock; numerical simulation; fluid-solid interaction; seepage characteristics

2016-03-12；

2016-04-26

国家自然科学基金资助项目(41402306)；广西安全生产科技项目(gxaj201405)

李晓泉(1972-)，男，黑龙江七台河人，广西大学副教授，博士； E-mail：lxq0927@gxu.edu.cn。

许大洋，李晓泉，王兄威，等.基于COMSOL的顶板运移规律及渗流特性研究[J].广西大学学报(自然科学版)，2016，41(5)：1714-1723.

10.13624/j.cnki.issn.1001-7445.2016.1714

TD73；TD745

A

1001-7445(2016)05-1714-10