水资源的保护与利用是我国煤矿绿色开采亟待解决的难题之一，特别是在我国西部地区，煤炭资源约占全国总储量的80%，但是气候干旱，蒸发量是降雨量的6倍左右，再加上西部矿区大规模、高强度开采对其生态水环境造成严重的影响^[1-2]。因此针对煤炭资源开采与水资源保护相协调的难题，顾大钊等^[3-5]提出了煤矿地下水库的概念及矿井水地下储存利用的理论与技术，有利于西部矿区煤炭开采与水资源保护。

针对煤矿地下水库的稳定性，诸多学者进行了大量的研究，取得了丰富的成果。姚强岭等^[6-7]研发了无损浸水实验装置，分析了不同含水状态的煤柱坝体损伤特征，并结合理论分析研究了煤矿地下水库煤柱坝体的合理宽度计算方法及稳定性；张村等^[8]、HAN等^[9]开发了采动水浸作用下煤矿地下水库流固耦合数值模型，并通过定义煤柱损伤度研究了工程尺度下煤柱坝体的失稳特征；师维刚等^[10]通过对防水隔离煤柱进行合理分区提出了计算煤柱宽度理论方法；白东尧等^[11]通过数值模拟及理论分析得出了地下水库人工坝体的极限水头，并对地下水库人工坝体的稳定性展开研究；刘少伟等^[12]通过理论计算及数值模拟研究了影响缓倾斜煤层沿含水采空区掘巷时影响护巷煤柱稳定性的因素敏感性；王文才等^[16]通过实验室测试及数值模拟研究分析了作业条件下矿井地下水库坝体衔接处受力，变形情况及其失效机制；李刚等^[14]通过理论分析及数值模拟相结合研究了采动过程中底板应力场分布特征及岩层破坏损伤演化特征；刘祥龙等^[15]针对圆柱形试件进行不同侧向渗透压力下的煤样压缩试验，研究煤样的破坏特征与渗透特性，将煤样在侧向渗透-轴向承载下的变形分为5个阶段；刘晓飞等^[16]通过理论分析及实验室测试等分析了工作面末采阶段煤柱坝体所受应力环境以及研究了周期应力及水浸作用对媒体的稳定性影响。

综上，对于地下水库煤柱坝体的稳定性研究多集中在实验室尺度下对浸水煤岩样强度及损伤特性研究以及通过理论分析对影响煤柱坝体稳定的临界水头研究，而在工程尺度下，部分学者通过数值模拟提出了力学计算与渗流计算交替进行的间接流固耦合方法，研究采动水浸作用下煤柱坝体的稳定性^[8-9]。然而，工程尺度下煤柱坝体在完全流固耦合作用下对其承载能力的动态劣化特征及损伤破坏规律的影响仍值得深入研究。同时基于煤柱的多孔介质特性，FLAC^3D数值模拟软件在模拟地下水库煤柱坝体在复杂载荷作用下处于弹性未破坏阶段时的稳定性存在一定的局限性。为此，通过数值模拟及二次开发，构建了基于煤矿地下水库煤柱坝体在流固耦合中多场耦合作用下力学参数动态弱化的Fish语言并嵌入FLAC^3D中，实现煤柱坝体在流固耦合作用下力学参数的动态调整以减弱数值模拟进行流固耦合分析时的不足；从而研究了地下水库煤柱坝体在采动及水浸叠加作用下的稳定性，得到了地下水库煤柱坝体在复杂条件下应力分布响应特征及塑性区扩展变化规律，为地下水库煤柱坝体稳定性研究提供借鉴。

1.   工程背景和流固耦合数值模型的构建

1.1   工程背景

以陕西北部某矿主要开采煤层2-2中的2个相邻工作面开采为例，煤层倾角为0°～5°，煤层厚度为8.89～11.73 m，平均厚度为10.41 m。双面开采及采空区蓄水过程如图1所示。

图  1  双面开采及采空区蓄水过程

Figure  1.  Dual-sided mining and goaf water storage process

下载: 全尺寸图片 幻灯片

P122106工作面为首采工作面，走向长度为6 200 m，倾斜长度为350 m，开采高度为6 m；P122106长壁工作面南部为P122108长壁工作面，在2个工作面之间留设安全煤柱，煤柱拟留设宽度为40 m，并在工作面开挖稳定后，开始在工作面采空区进行水库建设并进行蓄水。因此，对地下水库煤柱坝体在采动作用及上覆载荷-水浸流固耦合依次作用下的稳定性研究对地下水库长期安全运行非常重要。

1.2   数值模型的构建

按照前述工程概况，建立FLAC^3D流固耦合数值仿真计算模型。模型以x方向为倾向，y方向为工作面走向，z方向为垂直于xy平面；模型尺寸为x=1 000 m，y=400 m，z=145 m；在模型的四周及底部边界采用位移约束，上部为自由边界，共模拟2个相邻工作面的开采，为了使煤柱的监测区域更为准确，将其网格细化为0.5 m。数值模拟中材料参数取值见表1，长壁开采工作面数值模型如图2所示。

表  1  数值模拟中材料参数取值

Table  1.  Values of material parameters in numerical simulation

岩性	密度/ (kg·m−3)	体积模量/ GPa	剪切模量/ GPa	内摩擦角/ (°)	黏聚力/ MPa	抗拉强度/ MPa	渗透系数/ (nm·s−1)	孔隙率
中粒砂岩	2 500	11.61	9.09	28	1.00	2.30	0.060	0.45
砂质泥岩	2 400	2.57	1.45	26	0.55	0.80	0.006	0.35
粉砂岩	2 500	9.63	8.09	30	2.50	0.84	50.000	0.40
泥岩	2 141	2.21	1.24	24	0.52	0.83	65.000	0.40
细粒砂岩	2 400	2.13	0.93	26	0.45	0.35	0.006	0.30
煤层	1 400	2.13	0.93	25	0.40	0.25	80.000	0.50

下载: 导出CSV 

| 显示表格

图  2  长壁开采工作面数值模型

Figure  2.  Numerical model of longwall mining face

下载: 全尺寸图片 幻灯片

模型的中上覆各岩层及底板均采用Mohr-coulomb模型，煤层采用基于Mohr-coulomb的应变软化模型；同时，由于采空区逐渐压实过程影响工作面周围应力分布，开发了基于salamon应力应变曲线的模拟采空区垮落破碎矸石的承载特性的Fish语言嵌入FLAC^3D中。

1.3   基于煤柱坝体的水浸软化模型

对地下水库煤柱坝体进行流固耦合分析时，应力和水压是影响煤柱稳定的2个关键物理场，两者的耦合作用是导致煤柱坝体失稳破坏的重要原因；而FLAC^3D流固耦合中通过监测模型内部存在孔隙水压力和处于压缩状态时对其力学参数进行软化。实际上，煤柱坝体在弹性状态下，内部处于压密状态下的孔隙闭合会形成弹性压密隔水区，其渗透率显著降低，为煤柱的主要承载区域，此时，这部分煤体应表现为应变硬化特性。

为了更加准确的表现地下水库煤柱坝体在流固耦合作用下的破坏失效形式，对基于煤柱坝体的应变软化模型进行了二次开发，得出煤柱坝体内部应力场、塑性区、能量场多场耦合作用下力学参数动态更新的Fish语言，并嵌入FLAC^3D中。基于煤柱的水浸软化模型计算流程图如图3所示。

图  3  基于煤柱的水浸软化模型计算流程图

Figure  3.  Calculation flowchart of water immersion softening model based on coal pillar

下载: 全尺寸图片 幻灯片

Fish语言能够捕捉煤柱坝体内部各单元的弹塑性状态：①当目标单元的处于塑性状态时，根据现有研究^[17]确定煤柱在饱和状态下的弹性模量弱化系数k_E=0.429 2, 黏聚力弱化系数k_c=0.572 5, 内摩擦角弱化系数k_f=0.960 1，从而模拟煤柱内部塑性破坏区承载能力的劣化及失效特征；②当目标单元处于弹性状态且体积应变小于0，确定煤柱各力学参数的弱化系数为1，从而减弱FLAC^3D中流固耦合效应对其力学参数弱化的影响，最终实现对煤柱内部弹性压密隔水区的模拟。

2.   煤柱坝体稳定性

2.1   多工作面采动影响

由于矿方拟留设煤柱宽度为40 m左右，为了研究采动及水浸叠加作用下煤柱坝体的承载特征及渐进失效规律，以煤柱坝体宽40 m作为稳定性的研究对象，同时，为了对比研究煤柱坝体在水浸及采动叠加作用下的承载特征及失稳规律，还研究了煤柱留设宽度小于40 m（20 m和30 m）时的稳定性。煤柱坝体塑性区和垂直应力分布如图4和图5所示，不同宽度煤柱坝体垂直应力分布特征与弹性核区占比如图6所示。

图  4  煤柱坝体的塑性区分布

Figure  4.  Plastic zone distribution of coal pillar dam

下载: 全尺寸图片 幻灯片

图  5  煤柱坝体垂直应力分布云图

Figure  5.  Vertical stress distribution of coal pillar dam

下载: 全尺寸图片 幻灯片

图  6  不同宽度煤柱坝体垂直应力分布特征与弹性核区占比

Figure  6.  Characteristics of vertical stress distribution and percentage of elastic core in the coal pillar dam with different widths

下载: 全尺寸图片 幻灯片

由图6所可以看出，3种宽度的煤柱坝体垂直应力分布曲线均呈现双峰形分布特征，且煤柱右侧的垂直应力峰值均大于左侧，煤柱宽度为20、30、40 m时其垂直应力峰值分别为44.00、35.15、29.13 MPa，呈现出随着煤柱宽度的增加应力峰值随之减小。煤柱坝体2个垂直应力峰值之间为弹性核区，煤柱宽度为20、30、40 m时，弹性核区占比分别为39.73%、73.68%、83.54%。为了保证煤柱坝体的稳定性，煤柱内部需至少维持一定的弹性核区，现场实践表明隔水煤柱内部弹性核区占比应不小于40%，以保证在复杂载荷作用下保持稳定^[18]，并有效防止动力灾害的发生。因此，在开采过后，煤柱坝体留设宽度为20 m时，承载能力劣化并会发生失稳。表明随着煤柱宽度的增加其承载能力随之增强，煤柱的稳定性越好，且煤柱坝体留设宽度为30 m和40 m时能够保持稳定。

2.2   水浸作用下煤柱坝体稳定性

在工作面开挖稳定后，开始向采空区进行蓄水。此时，地下水库煤柱坝体会在水浸及上覆岩层载荷的流固耦合作用下强度进一步降低，且不同于传统的承载式煤柱，地下水库煤柱坝体在长期水浸作用下更容易失效；因此，在工作面开挖稳定后在采空区通过“watertable”命令建立水平面及各渗流参数，并根据上述煤柱流固耦合数值模拟方法，对不同宽度煤柱坝体在水浸作用下破坏失稳特征展开研究。渗流作用下不同宽度煤柱坝体塑性区如图7所示，渗流作用下不同宽度煤柱坝体垂直应力分布云图如图8所示，渗流作用下不同宽度煤柱坝体垂直应力及塑性区分布特征如图9所示。

图  7  渗流作用下不同宽度煤柱坝体的塑性区分布

Figure  7.  Plastic zone distribution of coal pillar dam with different widths under seepage effect

下载: 全尺寸图片 幻灯片

图  8  渗流作用下不同宽度煤柱坝体垂直应力分布云图

Figure  8.  Vertical stress distribution of coal pillar dam with different widths under seepage effect

下载: 全尺寸图片 幻灯片

图  9  渗流作用下不同宽度煤柱坝体垂直应力及塑性区分布特征

Figure  9.  Characteristics of vertical stress and plastic zone distribution of coal pillar dam with different widths under seepage effect

下载: 全尺寸图片 幻灯片

在水浸作用下，地下水库煤柱坝体宽度为30 m时，垂直应力曲线呈类似“拱形”分布，“拱顶”部分为煤柱坝体的主要承载区域，也是主要隔水区域。应力峰值为26.75 MPa，位于煤柱左侧，相较于采空区蓄水之前，降幅为23.9%；随着煤柱宽度的增加，垂直应力分布曲线中部主要承载区域的宽度也相应增加，应力峰值为24.11 MPa，位于煤柱左侧，相较于煤柱宽度为30 m时有所下降，采空区蓄水之前降幅为17.2%，结果表明，煤柱坝体宽度越小，流固耦合作用对其承载能力的影响越大，增加煤柱宽度有利于保持其稳定性。

在水浸作用下，不同宽度煤柱坝体的塑性区向煤柱内部进一步扩展，塑性区边界煤柱坝体内部为高垂直应力区域构成煤柱坝体的有效阻水区域。当煤柱宽度为30 m时，弹性核区占比为38.33%，相较于采空区蓄水之前降幅为48%；随着煤柱坝体宽度的增加其弹性核区占比显著增加；当煤柱宽度为40 m时，弹性核区占比为59.53%，相较于采空区蓄水之前降幅为28.7%。结果表明，煤柱宽度越小，流固耦合作用对其强度的影响越大，增加煤柱宽度有利于保持其稳定性。然而，当煤柱留设宽度过大时，一方面会浪费大量的煤炭资源，另一方面会产生一些煤柱型动力灾害。因此，为了保证地下水库蓄水后的系统安全，当煤柱宽度为40 m时，内部弹性核区占比将大于40%，能够在采动及水浸依次作用下保持稳定。

采空区积水水浸作用下煤柱内部力学参数大幅弱化，导致塑性区进一步扩展及承载能力显著降低。运行过程中不同宽度煤柱塑性区扩展演化规律如图10所示，运行过程中煤柱宽度为40 m时黏聚力弱化过程如图11所示。

图  10  运行过程中不同宽度煤柱塑性区扩展演化规律

Figure  10.  Evolution law of plastic zone expansion of coal pillar with different widths in the process of operation

下载: 全尺寸图片 幻灯片

图  11  运行过程中煤柱宽度为40 m时黏聚力弱化过程

Figure  11.  Weakening process of cohesion when the width of coal pillar is 40 m during operation

下载: 全尺寸图片 幻灯片

由图11可知：煤柱坝体右侧的损伤范围逐渐由边缘向内部扩展，并呈现出煤柱下方区域先于上方区域，边界煤体先于内部煤体的渐进失效特征，同时煤柱边界区域受开采扰动作用的影响较大，初始及最终损伤破坏程度均大于煤柱内部区域，且在采动及水浸作用下煤柱坝体中部区域始终未受到影响，进一步表明煤柱坝体宽度为40 m时能够保持稳定。

综上所述，采空区地下水库蓄水后，积水将沿着煤柱坝体在采动作用下形成的塑性区高导水通道侵入煤体，并向煤柱内部不断渗透，使煤体内部原有结构的黏结强度减弱以及裂隙面间摩擦系数减弱，从而增加了煤体破坏的可能性，导致煤体承载能力降低，从而进一步导致煤柱塑性区扩展并引起应力重分布向煤柱内部转移，并促使积水向煤柱内部继续渗透，直至煤柱内部高应力区形成稳定阻水区域或煤柱内部积水区域贯通造成煤柱坝体破坏。为此，构建的数值模型能够很好地体现地下水库煤柱坝体在采动及水浸依次作用下的渐进破坏过程。在采空区积水水浸作用下，煤柱坝体塑性区占比增量呈现“稳定-增大-减小-稳定”的变化规律，能够很好地模拟水浸作用下煤柱坝体的“渗流-弱化-损伤-渗流”的渐进破坏特征。

3.   结　语

1）基于FLAC^3D自带的流固耦合分析在模拟煤柱弹性孔隙压密阶段存在的不足，开发了基于煤柱多场耦合作用下参数动态调整的Fish语言并嵌入FLAC^3D中，实现煤柱坝体在流固耦合作用下力学参数的动态调整以实现煤柱在弹性压密阶段的模拟。

2）研究了多工作面开采对地下水库煤柱坝体稳定性的影响规律，结果表明：在相邻工作面开采后，不同宽度的隔水煤柱应力分布曲线均呈现双峰形分布，并随着煤柱宽度的增加其峰值应力逐渐降低，同时煤柱内部弹性核区占比也随煤柱留设宽度的增加而增大，当煤柱宽度为20 m时其弹性核区占比为39.73%，将发生破坏失稳，而煤柱宽度为30 m和40 m时煤柱能够保持稳定。

3）运用流固耦合方法对煤柱蓄水水浸作用下的稳定性进行了分析。结果表明：采空区蓄水后，水浸作用下煤柱坝体塑性区占比增量呈现“稳定-增大-减小-稳定”的变化规律，能够很好地模拟水浸作用下煤柱坝体的“渗流-弱化-损伤-渗流”的渐进破坏特征，水浸作用对煤柱坝体稳定性的影响更为显著，是导致其破坏失效的主要因素；同时煤柱内部应力峰值降低，煤柱宽度越小应力峰值降低幅度越大，且煤柱宽度为30 m时其垂直应力分布曲线从双峰形演化为拱形分布，同时其弹性核区占比降为38.33%，降幅为48%，将发生失稳；而煤柱宽度为40 m时，能够在采动及水浸依次作用下保持稳定。