基于Rhino3D-FLAC3D 耦合技术的库区下伏综放采煤沉陷区稳定性分析

黄健丰，吴 璋，张振振，2，李永强，张 浩

（1.中煤科工集团西安研究院有限公司，陕西 西安710077；2.煤炭科学研究总院，北京100013）

随着我国国民经济的持续发展，区域性交通工程、能源工程、水利水电工程等选址规划常与采煤沉陷区产生叠压现象，尤其部分采煤沉陷区工程地质条件复杂、地表移动监测数据缺失，使外部荷载作用下采煤沉陷区稳定性评价存在巨大困难，对建（构）筑物的建设及运营构成重大安全隐患[1-3]。国内外相关学者对采煤沉陷区稳定性开展了大量研究工作，取得了系列研究成果。李东阳等[4-5]采用砌筑法建立了浅部采煤沉陷区三维相似实验模型，研究了地面荷载作用下浅部采煤沉陷区覆岩破坏演变过程；杨利民等[6]提出了矩形截面采空区上方路基瞬时沉降预计公式，计算结果与现场监测数据一致性良好；邓喀中等[7-8]基于采空区覆岩移动规律建立空洞空隙预计模型，可为老采空区稳定性评价提供参考；张志祥等[9-10]通过相似模拟实验研究了交通荷载作用下采煤沉陷区覆岩“活化”规律；张向东等[11]基于弹性板理论建立了动静荷载耦合作用下采煤沉陷区活化力学模型，研究了地震与采动耦合作用下采煤沉陷区覆岩变形特征；Seryakov[12]研究了采空区覆岩应力状态，提出了一种围岩破坏状态估计方法；Rezaei[13]基于最小势能原理研究了导水断裂带稳定条件，推导了长壁开采地表最大下沉量预计公式。上述研究推动了采煤沉陷区在外加荷载作用下稳定性研究的进展，但关于我国黄土高塬综放采煤沉陷区在外部荷载作用下稳定性研究还较少，不能完全满足工程实践需求。

针对上述问题，以红岩河水库下伏综放采煤沉陷区为研究对象，联合地表等高线数据及勘探地质信息通过Rhino3D快速构建三维数值模型，将其导入FLAC3D数值模拟软件，结合钻孔冲洗液漏失量实测法，开展了综放采煤沉陷区覆岩裂隙发育特征及库区水体荷载作用下采煤沉陷区稳定性模拟研究，以期为类似采煤沉陷区稳定性评价及防治技术提供参考。

1 工程概况

红岩河水库位于彬州市西北部，是以城镇居民生活、工业生产供水和拦沙为主要用途的Ⅲ等中型水利工程，能够有效缓解居民生活供水及工业用水保证率低的现状，为彬州市打造关中能源基地发挥巨大作用。勘查表明，库区淹没线与火石咀煤矿8712 综放采煤沉陷区存在叠压现象，采煤沉陷区覆岩变形对水库蓄水及运营构成重大安全隐患，亟需研究水体荷载作用下采煤沉陷区稳定性，研究区域地理位置图如图1。

图1 研究区域地理位置图Fig.1 Location map of the study area

火石咀煤矿采用倾斜长壁综放开采技术，主采煤层为侏罗系延安组4-2煤层，4-2煤为近水平煤层，平均厚度为厚7.1 m，结构简单。8712 工作面采用自然垮落法管理顶板，倾向长度为210 m，走向长度为2 250 m。地表起伏多变，沟壑纵横，大部分被第四系、新近系松散层覆盖，红岩河沟内侧有华池组、洛河组岩层出露，井田内地层由老至新依次为：侏罗系富县组（J1f）、延安组（J2y）、直罗组（J2z）及安定组（J2a）；白垩系宜君组（K1y）、洛河组（K1l）、华池组（K2h）；新近系（N）；第四系（Q）。

2 三维数值模型可视化建模方法

2.1 模型构建流程

FLAC3D是Itasca 公司基于连续介质理论及显式有限差分法开发的三维数值模拟软件，广泛应用于岩土工程、矿业工程、地质工程等领域，尤其适用于复杂应力条件下岩土体变形破坏机理、地质灾害孕灾机制及安全风险评估等研究，然而该软件存在前期建模能力较弱、灵活性差等问题，难以满足构建大型复杂地质体三维数值模型的需求。针对上述问题，众多学者[14-18]对三维数值模型构建方法展开研究，取得系列研究成果。但多依托数个建模软件，操作复杂，且通过计算机语言编写数据转化接口，对使用者有较高编程水平要求。将矢量化地形图、钻孔地质信息及采掘工作面轮廓导入Rhino3D工作空间，通过布尔运算对已有信息进行处理、融合，构造相应的地质体和结构面，三维地质体建模流程如图2。

图2 三维地质体建模流程Fig.2 3D geological body modeling process

2.2 三维数值模型构建

从不连续的原始数据点源重新建立连续的地质面，通常采用插值法或近似拟合法。插值法通过给定数据点刚性地形成数据曲面，不改变数据点的位置坐标，在原始点源数据较少时外推精确度较差；近似拟合法是通过某种数学方法构建最大近似平滑表面，这种方法会改变部分数据点的坐标，不依赖于点源数量，但难以达到高精度标准。因此，采用2 种方法相结合构建复杂地质体数值模型。

以研究区域地表模型构建为例，首先通过三维GIS 软件收集研究区域数字化地形图，提取地表矢量化等高线数据，将等高线数据导入Rhino3D工作空间离散提取点源三维坐标，地表等高线点源如图3。

图3 地表等高线点源Fig.3 Point sources of contour lines

提取点源初始数据集合为T=｛Rs，（s=0，1，2，…，m-1）｝，分布均匀子集为T1=｛Ki，（i=0，1，2，…，r）｝，分布不均匀子集为T2=｛Kj，（j=0，1，2，…，m-1-r）｝。集合T 采用插值蒙皮法构造曲面，不改变点源的坐标，其中子集T1可形成均匀曲面，由于子集T2点源不均造成曲面存在尖锐区域，不利于后续数值计算分析，插值法构建曲面如图4。

图4 差值法构建曲面Fig.4 Difference method to construct a surface

Rhino3DNURBS 技术通过曲线起点、终点、控制点及编辑点完成曲线的编辑，曲线的形状主要由曲线的控制点和编辑点决定，通过Rhino3DNURBS 技术利用给定的精度以及控制点源数据可实现地表插值-逼近拟合完整的NURBS 拟合地表曲面如图5。

图5 NURBS 拟合地表曲面Fig.5 NURBS fitting surface surface

从地表某端点沿着z 轴负方向建立体模型，结合工程地质资料及现场钻孔勘查数据对地质体进行切割形成地层界面，通过切割、组合构建的研究区域三维数值模型如图6。x 轴方向1 640 m，y 轴方向910 m，z 轴方向430～662 m，8712 开采工作面走向长度为600 m，倾向长度为210 m，采厚为7.1 m，六面体网格共686 504 个，节点共713 460 个。

图6 研究区域三维数值模型Fig.6 3D numerical model of study area

3 综放采煤沉陷区覆岩裂隙发育特征

基于开采损害理论，采煤沉陷区覆岩裂隙主要分为边缘裂隙和动态裂隙。边缘裂隙一般位于沉陷盆地的边缘区域，以“带状”、“椭圆状”的形态分布在沉陷盆地的边缘；动态裂隙则位于工作面上覆岩层中，并随着工作面的推进不断扩展，直至覆岩形成新的平衡结构，其中导水断裂带覆岩存在较大自由空间，为矿区地下水资源流失提供通道。库区水体长期作用下沿沉陷盆地边缘拉张裂隙带补给含水层，通过原生裂隙、与导水断裂带间接沟通形成的水力通道，对水库运营构成重大安全威胁。因此，必须查明采煤沉陷区导水断裂带分布特征。

采用钻孔冲洗液漏失量观测法获取导水断裂带发育高度，判别标准为：①钻孔冲洗液漏失量急剧增加；②吸风现象明显。在采煤沉陷区地表布设探测孔K-14，地面高程为910 m，孔底高程为468.7 m，通过统计分析钻孔冲洗液消耗曲线图（图7），孔深为306 m 时漏失量由0.14 L/（m·s） 急剧上升为18.49 L/（m·s），吸风现象明显，认为已经到达导水断裂带顶端。确定8712 工作面导水断裂带高度为135.3 m，裂采比为19.05。

图7 钻孔冲洗液消耗曲线图Fig.7 Consumption curve of drilling fluid

4 库区下伏综放采煤沉陷区稳定性

4.1 边界条件及力学参数选取

采用Mohr-Coulomb 弹塑性本构破坏准则，利用FLAC3D软件新功能对数值模型表面进行分组，将模型East 和West 面作为x 轴方向的边界条件，模型South 和North 作为y 轴方向的边界条件，Bottom 面约束z 轴边界条件，TOP 面为自由边界。参考工程地质勘查资料及相邻井田岩石物理力学实验数据，确定模拟的岩体物理力学参数见表1。

4.2 原岩应力场反演情况

初始应力场计算采用弹性本构准则，计算数值模型仅在自重作用下的应力状态，得到研究区域的初始应力分布如图8。沿工作面倾向中线y=505 m布置剖面，最大垂直应力为13.97 MPa，应力场呈“波浪形”层状分布，是由于黄土塬凸起埋深增加造成的。

4.3 8712 工作面导水断裂带发育特征

为深入分析随工作面推进覆岩变形破坏全过程，采用分步开挖的方式对8712 工作面开采过程进行模拟，赋予开挖部分为空本构模型，每次工作面推进40 m，共开挖10 步，累计开采400 m。沿工作面倾向中线y=505 m 布置剖面，8712 工作面塑性区演化分布如图9。

表1 岩体物理力学参数Table 1 Rock physical and mechanical parameters

图8 初始应力分布Fig.8 Initial stress field distribution

图9 8712 工作面塑性区演化分布Fig.9 Evolution distribution of 8712 working face plastic zone

通过分析8712 工作面推进过程的塑性区演化分布，采动覆岩先以剪切破坏为主，导致顶板裂隙初步发育，逐步发展成为拉伸破坏产生垮落带。随着工作面的推进，拉伸破坏区的范围逐渐增大，上部剪切破坏区域也逐渐扩大，工作面前方煤壁和开切眼处的围岩主要发生了剪切破坏，覆岩塑性区分布整体成“马鞍”型。塑性区发育高度为113.92～138.24 m，最大裂采比为19.47，与钻孔冲洗液漏失量观测法实测导水断裂带高度（135.3 m）基本一致，反映出三维数值模型、本构条件及岩体力学参数选取的合理性。

4.4 库区水体荷载作用下采煤沉陷区稳定性

在模型上表面按照正常蓄水位（903 m）施加静水压力，库区水体荷载分布如图10。模型采用Mohr-Coulomb 破坏准则，将模型四周及底部作为边界条件，TOP 面为自由边界。

图10 库区水体荷载分布Fig.10 Water load distribution in the reservoir area

库区水体荷载对采煤沉陷区的覆岩位移及塑性区影响如图11、图12。

图11 蓄水前后覆岩位移对比Fig.11 Comparison of overburden displacement before and after water storage

图12 蓄水前后覆岩塑性区对比Fig.12 Comparison of plastic area of overburden before and after water storage

库区水体荷载作用下，采煤沉陷区覆岩破坏范围增大，破坏逐渐向上发展，使原本处于稳定状态的采煤沉陷区产生“活化”现象，沿工作面倾向中线y=505 m 布置剖面，水库蓄水后，地表垂向位移最大值为686.4 mm，比蓄水前垂直位移增加334.7 mm。导水断裂带最大发育高度为142 m，比蓄水前增加3.76 m。

5 结 语

1）采用Rhino3D-FLAC3D耦合技术结合8712 综放采煤沉陷区地形特征及地层岩性快速构建了研究区域三维数值模型。

2）8712 综放采煤沉陷区覆岩导水断裂带整体呈“马鞍”状，最大发育高度为138.24 m，裂采比为19.47，与钻孔冲洗液漏失量观测法实测值（135.3 m）基本一致。

3）库区水体荷载作用下，综放采煤沉陷区地表垂向位移增加334.7 mm，导水断裂带发育高度增加3.76 m。