基于FLAC3D的贵州马达岭煤炭采动型滑坡机理研究

2019-12-24潘网生傅良同张晓周韩国栋

煤炭工程 2019年12期

潘网生，傅良同，李鑫，张晓周，韩国栋

(1.黔南民族师范学院旅游与资源环境学院，贵州都匀 558000；2.长安大学环境科学与工程学院，陕西西安 710054)

贵州省喀斯特出露面积约11×104km2，约占全省国土面积61.9%[1-3]，其喀斯特地貌形态以锥状为主，类型多样、形成发育过程相对复杂，其形态结构独特、锥形明显，单体高度为60～250m，坡度为45°～ 47°[4-6]。近几年，贵州喀斯特区域滑坡地质灾害频发，生命与财产损失巨大，如关岭县山体滑坡(2010年6月28日)，死亡99人；福泉市道坪镇英坪村大型山体滑坡(2014年8月27日)，死亡23人；大方县理化乡偏坡村山体滑坡(2016年7月1日)，死亡21人；玉屏县田坪镇岩屋口村山体滑坡(2017年7月7日)，死亡2人；纳雍县张家湾镇大型山体滑坡(2017年8月28日)，死亡15人；三荔高速荔波段滑坡(2017年9月23日)，死亡3人，六盘水六枝特区关寨镇大型山体滑坡(2018年7月9日)，因预警及时以致无人员伤亡。上述滑坡地质灾害已经成为制约贵州社会经济可持续发展的重要因素之一。由于喀斯特滑坡灾害机理研究可以揭示滑坡灾害成因、发展及演变规律，也可以为早期预警和地质灾害防治提供科学依据和理论支撑。因此，研究煤炭采动型滑坡机理具有重要的实践意义。

国内学术界关于贵州喀斯特滑坡机理研究，总体概括为岩性和岩土结构等(内因)，降雨和不科学的人类工程活动等(外因)。分以下两种类型：其一是由自然地质沉积所形成的软弱结构层在地表渗流和地下水溶蚀作用下，其粘聚力和内摩擦角减小，成为天然滑动面，导致软弱结构层上方岩土体静力平衡被打破，进而失稳致滑。如贵州关岭“6.28” 喀斯特山体滑坡地质灾害等[7-13]；其二是人为斜坡开挖、矿产资源开采等人类工程活动导致边坡应力平衡被打破，在降雨或震动条件下诱发山体滑坡。如贵州印江岩口大型喀斯特山体滑坡地质灾害等[14-18]。

贵州都匀马达岭滑坡属于第二种类型—人为开挖煤层导致采空区上方岩体滑坡，也称煤炭采动型滑坡[19]。赵建军[20]基于物理模型研究其地质力学模式；李禹霏[21]通过自动化监测研究其稳定性；亓星[22]研究其滑坡碎屑流形成泥石流的特征；王玉川[23]对其开展岩石力学试验。本文在地质野外调查基础上，拟采用有限差分方法研究采动滑坡的变形过程及特征，以揭示其滑坡机理和驱动模式，相关研究成果可以应用于类似喀斯特矿山开采的安全管理和地质灾害防治。

1 数据来源及研究方法

研究区域地质资料来源于前人研究成果及贵州省区域地质志(1987版、2017版)[24，25]，地形地貌数据来源于国家地理空间数据云服务平台及野外现场勘察、测绘。基于CAD-ANSYS构建原始三维边坡模型，导入FLAC3D2600软件后模拟煤层开采对边坡稳定性的影响。

FLAC3D由美国ITASCA公司研究开发，是一款基于拉格朗日算法有限差分的显式数值分析模拟软件。该软件包含11种弹塑性材料本构模型，非常适合模拟岩土体的屈服、塑性流动、软化及大变形问题，可以快速求解场的力学方程，准确分析岩土边坡的稳定性[26]。拉格朗日差分计算原理是基于显式差分求解偏微分方程[27，28]，即通过运动方程计算节点运动速率，根据运动时间Δt，计算节点(ni)时步位移，由节点位移推算单元(ei)之间的相对位移(单元ei由边长lj，lj+1，…lj+k构成，边长lj由节点nk，nk+1和构成)，进而求解单元应变增量(△εi，j)。由单元材料本构方程获取应变条件下的应力增量(△σi，j)，在此基础上求解各节点的不平衡力(即合力)，根据中心差分原理[29]，计算节点的速度与加速度，进而再一次获得各节点的时步位移。按照上述步骤往复迭代，扩展到整个区域，直至各节点的不平衡力足够小或者各节点的位移趋于平衡[30]。

2 研究区域概况及地质模型

马达岭滑坡位于贵州省都匀市江洲镇富溪村干坝马达岭沟，沟口地理坐标为26°10′18.9″N，107°17′23.4″E[22](图1)。马达岭沟流域面积0.86km2，流域全长1.62km，最大高程1570m(北部山脊)，最小高程1150m(南部沟口)。沟道走向大致由北向南，中游呈“V”型，上、下游呈“U”型，沟道平均纵向比降约28.2%，上陡下缓。灾害点地处亚热带季风湿润气候，雨热同期，4—9月为雨季，年平均降雨量约1400mm，且主要集中在6—8月[22]。灾害点地层属滨海沉积相，地质构造相对稳定，沟口上方有一条东西走向的正断层。地层地质时代划归早石炭世大塘早期，为祥摆组(C1x)。地层结构由灰黄色的石英砂岩、黄褐色砂质页岩、灰黑色炭质页岩及黑色煤层组成，其中夹少量泥灰岩。地层主要特征为砂页岩夹煤层，且被认为是贵州石炭系出露最广、唯一含煤的岩组。由Ⅰ-Ⅰ′地质剖面(图2)可知，岩层产状280°∠18°，缓倾坡内。煤层主要有3层，自下而上分别为M1、M2和M3，其中M1和M2煤层厚度1.9～2.4m，M3煤层厚度0.2～0.9m不等(图3)。因煤炭开采，坡顶产生纵横交错、宽度不等的裂隙，发育形成两组优势结构面[20]，大致产状分别为34°∠87°、118°∠86°，受其影响滑坡体与山体拉裂，形成“圈椅状”滑坡后壁。滑坡位于马达岭坡顶，滑坡体后缘肩部高程约1573.6m，滑坡前缘坡脚高程约1481.5m，相对高差约92.1m(图3)。滑坡滑动主方向约204.6 °，滑坡体沿该主方向长度约88m，平均宽度102m，平均厚度92m，滑坡体体积达8.26×105m3。此外，滑坡物堆积于马达岭沟中部，堆积体平均厚度约6m[23]。

图1 滑坡区域平面图(滑坡前)

图2 Ⅰ-Ⅰ′剖面地层结构

图3 滑坡全貌图

由于FLAC3D构建不规则三维模型具有一定的难度，因此本文以Ⅰ-Ⅰ′剖面为基础，在ANSYS软件平台上恢复原始三维边坡概化模型，之后再导入FLAC3D(图4)。该边坡模型x轴方向长0～426.49m，y轴方向宽0～100m，z轴方向高1400.00～1577.09m，模型共产生381786个节点和365950个单元，分19个地层，分别为砂岩、炭质页岩和煤岩互层。岩性参考文献[23][31]，本文数值模拟计算时假设各层砂岩的岩性一致，各层炭质页岩的岩性一致，各层煤岩的岩性也一致(表1)。由于FLAC3D关于岩体变形模量仅接受体积模量K和剪切模量G，因此本文设定砂岩、炭质页岩及煤岩均弹塑性材料，服从莫尔-库伦屈服准则，并依据文献[32]确定体积模量K、剪切模量G与弹性模量E及泊松系数μ的关系。

图4 初始概化三维边坡模型

3 煤层开挖过程的边坡变形分析

为深入分析煤层开挖对边坡稳定性的影响，本文设计5种工况，分别探讨开挖进深、开挖作业面宽度及开挖煤层的先后次序等对马达岭边坡变形的影响。

表1 岩石物理力学参数

3.1 工况设计

1)工况一：开挖M1、M2煤层至x=300m，作业面宽度8m。通过FLAC3D模拟M1、M2煤层开挖过程，并监测边坡应力应变状态。M1、M2煤层开挖最终各形成5个工作面，每个工作面宽8m，工作面之间以煤柱作为支撑，其中，M1煤层开挖至x=300m，开挖进深为245m，M2煤层也开挖至x=300m，开挖进深为187m，工作面y轴方向的坐标分别为：4～12m；24～32m；44～52m；64～72m；84～92m，开挖位置及尺寸如图5(a)所示。模拟过程中，在边坡体中选取有代表性的点位A(256，0，95)、B(280，0，130)和C(360，0，140)(见图2)作为监测点，以记录节点竖向位移、节点合速率等参数随运行步的变化，其在格网中的节点编号分别为A(GP4104)、B(GP1815)、C(GP1414)。

2)工况二：开挖M1、M2煤层至x=400m，作业面宽度8m。M1煤层开挖至x=400m，开挖进深345m，M2煤层也开挖至x=400m，开挖进深287m，其他条件同开挖进深至x=300m一致，开挖位置及尺寸如图5(a)所示。

3)工况三：开挖M1、M2煤层至x=400m，作业面宽度增加至12～16m。M1煤层开挖至x=400m，开挖进深为345m，M2煤层也开挖至x=400m，开挖进深为287m，工作面12～16m不等，y轴方向坐标分别为：4～20，24～40，44～60，64～80，84～92，开挖位置及尺寸见图5(b)各工作面之间以煤柱作为支撑。其他条件同工况二。

4)工况四：在工况二基础上开挖M3煤层，开挖至x=300m。在工况一基础上开挖M3煤层，开挖至x=300m，M3煤层工作面在y轴方向分别为：4～12m；24～32m；44～52m；64～72m；84～92m，开挖位置及尺寸如图5(c)所示。

5)工况五：同时开挖M1、M2、M3煤层，其他开挖条件同工况四，开挖位置及尺寸如图5(c)所示。

图5 各工况煤层开挖位置及工作面宽度(m)

3.2 边坡变形结果分析

3.2.1 开挖进深对边坡变形的影响

根据工况一和工况二的边坡变形计算结果(图6)，当M1、M2煤层开挖进深由300m增加至400m，边坡整体变形增加了约2倍。开挖之初，坡脚前缘部分最先产生塑性变形，但由于坡脚前缘煤层上方的岩体厚度较薄，较小的势能和载荷使得坡脚变形相对较小。随着开挖进深增加，边坡变形越来越大，且变形最大区域逐步由坡脚前缘向坡顶和边坡内部转移。

图6 不同开挖进深下的边坡合位移

根据边坡内部A、B、C节点监测数据，节点产生位移的先后次序分别为A、B、C，且由于节点在y轴方向的水平位移最小，在x轴方向的水平位移次之，而在z轴竖向位移最大，因此本文选取竖向位移作为节点位移监测指标，另选取合速率作为节点速率监测指标。由工况一的边坡变形和节点竖向位移、合速率计算结果(图6、图7)可知，当M1、M2煤层开挖进深由300m增加至400m，节点合速率增加了约2倍。节点竖向位移与节点合速率与开挖进深均呈正比例关系(图7)，即随开挖进深的增加而增大。坡顶部位合位移将沿Z轴负向推进，其位移大小在任一Y剖面上沿X轴正向呈递减分布，即形成不同的合位移区间。因此，边坡岩体存在合位移差异，导致边坡垂向节理和裂隙形成并发育，且随着开挖进深的增加而增大，逐步向边坡内部推进，直至下切到M1、M2和M3煤层，形成潜在的优势滑动结构面。

图7 不同开挖进深下的节点竖向位移与合速率监测

3.2.2 开挖作业面宽度对边坡变形的影响

根据工况三的边坡合位移计算可知，当M1、M2煤层的作业面开挖宽度由8m增加至12～16m时(2倍)，边坡合位移增加了约10倍(图8)。

图8 增加作业面开挖宽度下的边坡合位移

根据边坡内部A、B、C节点竖向位移监测数据记录，增加作业面宽度后，A、B节点的竖向位移增大了约11倍，C节点的竖向位移增加了约12.3倍。此外，增加作业面宽度后(2倍)，峰值合速率增大约4倍(图9)。上述数据说明，边坡内部节点竖向位移与作业面宽度呈正比关系，即作业面宽度越大，节点竖向位移就越大，即作业面宽度对边坡变形的影响要远大于开挖进深。

图9 增加开挖宽度下的节点竖向位移与合速率监测

由此表明，随着M1、M2煤层开挖作业面宽度的增加，边坡合位移、节点竖向位移和节点合速率均呈不同程度增加，且节点竖向位移增幅远大于节点合速率增幅。

3.2.3 开挖煤层先后次序对边坡变形的影响

由工况四的边坡变形计算(图10a，b)可知，边坡最终合位移由两次开挖叠加导致，且开挖M3煤层所产生的合位移主要集中于M3煤层上方岩体，合位移约为单独开挖M1、M2煤层合位移的10 %～20%，说明深部煤层开挖对边坡变形的影响要远大于浅层煤层开采。此外，根据工况五的边坡合位移计算结果(图10c)，同时开挖煤层M1、M2和M3(工况五)产生的最大变形区域主要分布于边坡后缘中下部，其整体变形较分步骤开挖煤层产生的变形(工况四)稍大，节点竖向位移和节点合速率记录也验证此现象(图11)。总体而言，开挖煤层先后次序对边坡变形存在一定的影响，且同时开挖各煤层对边坡变形的影响程度较分步骤开挖各煤层要大，但相对于开挖进深和开挖作业面宽度而言，其对边坡变形的影响程度又相对较小。

图10 不同开挖次序情况下的边坡合位移

图11 不同开挖次序下的节点竖向位移与合速率监测

4 基于数值分析的滑坡机理推定

本文在现场调查和定期的边坡变形监测(GPS)基础上，开展基于FLAC3D的数值模拟分析，得到与现场调查、监测相一致的结论。说明本文选取边坡合位移、节点竖向位移和节点合速率等指标监测边坡变形过程具有可操作性和实际意义，为揭示边坡滑动机理和本质特征提供了重要的观测和量化工具。由此，借助该工具推定马达岭滑坡机理和过程如下：

单独开挖M1、M2煤层时，随煤层开挖进深的推进，边坡内部因变形位移差异产生沿x轴方向的若干条垂向优势结构面(e、f、g、h)(图12)。该类优势结构面受煤层开挖作业面宽度控制，随作业面宽度的增加而拉裂扩张。其向上延伸至坡顶，形成若干条地裂缝，成为坡面雨水入渗边坡的优势通道，向下则切割边坡，直至M1、M2煤层。

图12 滑坡机理示意图

在此基础上若进一步开挖M3煤层，边坡上方前缘再次遭受采掘扰动。此时位于M3开挖进深范围内、前期所形成的垂向优势结构面，将随M3煤层作业面宽度的增加而进一步拉裂扩张。受M3煤层采空影响，M3煤层上方岩体沿z轴负向产生位移，同时M3下方岩体因采空卸荷而反弹，产生若干卸荷拉张裂隙面，即横向优势结构面(d的下部)。在满足一定的降雨条件时，M3煤层采空区汇集大量裂隙水，一部分裂隙水沿M3煤层开挖工作面入口溢出，另一部分则沿横向优势结构面继续下渗。当M3煤层上方岩体发生失稳垮塌，在重力作用下冲击M3煤层下方岩体时，横向优势结构面锁固段被打通，优势结构面(d)在水、岩等综合作用下实现连续和贯通。此时，上方岩体沿优势结构面在剖面D处临空剪出，剪出岩体受重力作用抛向马达岭沟谷，遭受撞击后形成大小不等块体散落于沟谷中段，滑坡地质灾害由此发生。

马达岭滑坡地质灾害发生后，该边坡灾害隐患依然存在，甚至更为严峻。一方面，抛落在沟谷中的大小岩块堆积物在极端天气条件下可能会再次发生滑动，威胁沟谷下游安全。据2018年6月22日强降雨后的现场观测记录数据，该原滑坡堆积区(D下方)边缘产生5～10cm宽度不等的蠕滑裂缝。另一方面，开挖M1、M2过程中所产生的一系列垂向优势结构面e、f、g、h等，已经将边坡切割为若干部分，这些被切割岩体之所以能处于相对稳定状态，是因为边坡岩层反倾向于坡内，并缺少足够水文地质条件。然而，该边坡及后缘山体岩层产状为NW282°∠18°，在水、岩长期相互作用下，山体垂向优势结构面和岩层水平节理、裂隙会持续扩张。满足一定降雨条件后，地表水将沿后缘山体垂向优势结构面渗入岩体层状节理、裂缝，引起岩层顺层错动，直接挤压和推动前述边坡的被切割岩体沿优势结构面(如：g-M1)剪出，到时可能会形成规模更大的滑坡地质灾害。