基于PFC3D的鱼鳅坡滑坡运动过程分析

2021-08-26张家勇邹银先杨大山

中国地质灾害与防治学报 2021年4期

张家勇，邹银先，杨大山

（贵州省地质环境监测院，贵州贵阳 550001）

0 引言

贵州省地处西南腹地，是无平原支撑的省份，地形切割剧烈，地质环境脆弱，地质灾害频发[1−3]。2016年5月23日，贵州省开阳县龙岗镇大石板村鱼鳅坡组突发滑坡灾害，造成下方1 户居民2 人死亡，5 间房屋毁坏，10 余亩耕地损毁，直接经济损失90 万元（图1）。鱼鳅坡滑坡为堆积层失稳，是贵州省内典型突发性灾害类型之一。贵州及我国西部山区常见堆积体，在自然、人为因素影响下，容易诱发滑坡地质灾害，给人民的生命财产安全造成了严重的威胁。黄润秋[4]通过大量的研究表明：查清滑坡变形破坏的地质力学模式是滑坡地质灾害防治的基础所在。

图1 滑坡发生前后影像图Fig.1 Image before and after landslide

当前，针对滑坡机理分析主要是借助地质力学方法开展定性分析，判断破坏模式。近年来，随着计算机技术的发展，数值模拟方法在滑坡领域得到了广泛应用。目前的数值模拟方法中，主要包括确定性分析和非确定性分析两类方法，而确定性分析又可分为连续介质分析、非连续介质分析方法以及近年来李世海等[5]提出的基于连续介质力学的离散元方法(CDEM)。其中连续介质数值分析方法有有限单元法、边界元法、有限差分法等，非连续介质分析方法有块体离散元法、颗粒离散元法、关键块体理论、不连续变形分析(DDA 法) 等。CDEM 方法将有限元与离散元进行耦合，在块体内部进行有限元计算，在块体边界进行离散元计算，通过块体内部及块体边界的断裂，不仅可以模拟材料在连续状态下及非连续状态下的变形、运动特性，更可以实现材料由连续体到非连续体的渐进破坏过程。

离散单元法，它的基本思想是把滑体介质看作由一系列离散的、独立运动的单元所组成。20 世纪90年代以后，离散单元法发展十分迅速，基于不规则形状块体单元，如UDEC、3DEC 等软件，以及基于球形和圆盘形散体单元的软件，如PFC、MatDEM 等，得到了不断的完善，使其在许多的领域都得到了广泛的应用。颗粒流法（PFC）作为离散单元法的一种，通过圆形离散单元（Ball）来模拟颗粒介质的运动及其相互作用，克服了传统连续介质力学的宏观连续性假设，能真实直观地模拟滑坡的失稳破坏过程[6]。目前一些学者逐步将其应用于滑坡破坏运动分析之中。吴剑[7]采用PFC 对滑带剪切过程进行了模拟研究，张龙等[8]、施凤[9]、李冬冬[10]、王畯才[11]采用PFC3D对多个高速远程滑坡的运动过程进行了模拟，深入讨论了位移、速度、堆积形态等各方面的因素。杜永彬[12]、王宇等[13]、胡江春等[14]通过采用PFC2D对滑坡形成演化过程进行模拟，得到了良好的效果。赵洲等[15]、陈达等[16]、闫晓娟等[17]、汪华安等[18]采用PFC2D对堆积层滑坡形成演化过程进行模拟，并对堆积层滑坡的机理进行了总结[19−21]。

因此，本文利用颗粒流程序PFC3D模拟研究鱼鳅坡滑坡的破坏运动过程，选用Ball-Wall 建模方法，引入颗粒流(PFC3D)程序中平行黏结模型，进行滑坡模型的建立，对滑坡不同关键部位颗粒进行位移、速度监测，模拟分析鱼鳅坡滑坡破坏运动过程。研究成果可为对该类滑坡影响范围预测，以及工程措施的制定具有一定的参考意义。

1 滑坡概况

鱼鳅坡滑坡位于开阳县龙岗镇大石板村鱼鳅坡组。研究区属亚热带季风气候，冬季干冷少雨，夏季受季风影响，暖湿多雨。四季降雨分布不均，且降雨多集中在汛期4—9月，占全年降雨量的75.8%，滑坡发生前10 天的累计降雨量达到了70 mm。滑坡所处区域属斜坡地貌，总体地势为东北高，西南低，相对高差约210 m，平均地形坡度25°，地形起伏变化较大。滑坡区地层岩性主要为第四系残坡积层（Qel+dl）含碎石黏土和寒武系下统金顶山组(∈1j)粉砂岩、泥岩。鱼鳅坡滑坡地处大翁林背斜南东翼，靠近背斜核部，滑坡前缘靠近一条北西-南东向平移断层，左侧缘靠近一条北东-南西向的正断层。地下水分为碎屑岩裂隙水与松散层孔隙水，根据资料显示该区地下水主要富存在第四系松散堆积层中，其透水性及含水性较好。地下水以大气降水补给为主，地表水补给次之。

鱼鳅坡滑坡在平面上的投影大致呈上小下大的塔式形状，滑坡在剖面上呈现出近似直线状，主滑方向244°，滑坡斜长约140 m，平均横宽约60 m，平均厚度约4 m，总体积约40 000 m3（图2）。滑坡土体为寒武系牛蹄塘-金顶山组地层风化形成的残坡积土，碎石含量较少，碎石粒径较小（图3）。

图2 鱼鳅坡滑坡平面图Fig.2 Plan of Yuqiupo Landslide

图3 滑坡剖面图（I−I′）Fig.3 Section of landslide (I−I′)

鱼鳅坡滑坡的变形破坏是其独特的内外因共同作用的结果，坡体临空条件，上部疏松下部致密的岩土体组合，不良地表排水条件等为内因，短期强降雨入渗为外因。降雨为鱼鳅坡滑坡的直接诱发因素，斜坡变形破坏模式为上部的土体沿基岩面先出现变形，对下部滑体产生推挤作用，直至坡脚产生剪切破坏，并向上牵引发展，导致整体滑动。

2 滑坡破坏过程数值模拟

2.1 计算模型的建立

本文根据滑前1∶1 万地形图、滑后无人机航测DEM，生成滑体和滑床的几何边界模型，并导出为PFC3D能识别的STL 格式，通过“Geometry import”命令导入PFC 中作为两个“Geometry”。选择Ball-Wall 模型建模，其中滑床的“Geometry”作为Wall 模型。然后通过有限步循环计算消除颗粒间的内力，得到滑坡最终计算模型（图4）。

图4 滑坡数值计算模型Fig.4 Numerical calculation model of landslide

该滑坡模型长230 m，宽170 m，高115 m，颗粒最小半径0.25 m，最大半径0.5 m，颗粒总数为2 447 个。

为了从整体和局部两方面共同把握滑体的运动特征，此次模拟在滑坡前缘、滑坡中部、滑坡后缘共设置了6 个监测点，分别用于监控滑坡过程中滑体的速度、位移的变化（图4）。

2.2 细观力学参数的确定

PFC 滑坡模拟参数主要包括本征参数和接触参数，其中本征参数主要包括剪切模量、泊松比和密度，可通过室内岩土体物理力学试验和工程类比确定，颗粒单元间的接触参数主要通过以下两个途径获取[12]：（1）通过PFC 进行参数标定确定物理力学参数。（2）通过现场参数反演，即对所需微观参数进行赋值使颗粒流的运动状态与实际状态基本一致，可认为参数为合理值。本文主要选取第二种现场反演的方法确定颗粒单元的接触参数，具体取值见表1。

表1 滑坡岩土体颗粒细观参数Table 1 Mesoscopic parameters of rock and soil particles in landslide

2.3 破坏过程模拟及分析

模拟计算到1 s 时，坡体出现变形，边坡整体上位移较小，滑坡处于蠕滑变形阶段，见图5（a）。3 s 时，坡体变形进一步增大，且坡体中后部颗粒的位移较中前部颗粒大，见图5（b）。到达5 s 时，滑坡前缘产生剪切破坏，坡脚颗粒开始剪出，见图5（c）。7 s 时，滑坡滑动带已基本贯通，在坡体压力作用下滑坡开始整体滑移，见图5（d）。到达8 s，滑坡整体的速度差进一步减小，滑坡处于整体加速阶段，见图5（e）。到达10 s 时，速度峰值12.4 m/s，此时，近一半的滑体从滑源区滑出，见图5（f）。到达15 s时，由于地形原因，前缘颗粒的速度逐渐减小，见图5（g）。30 s 时，滑坡速度进一步减小，前缘的滑坡已基本停止运动，见图5（h）。到50 s 时,滑坡停止运动并堆积于坡脚平缓处，见图5（i）。

图5 不同时步速度云图（单位：m）Fig.5 Velocity cloud map of different steps (unit: m)

综上可知，该滑坡在破坏初始阶段以蠕滑变形为主，随着变形量的增加，滑坡体不断挤压坡脚，滑坡岩土体到达应力平衡极限，坡脚产生剪切破坏，并向上牵引发展，滑坡发生整体滑动，整个滑坡历时50 s，整体运动距离约为80 m，滑坡运动距离较短。整个滑坡破坏模式表现为牵引式破坏。模拟结果与实际变形情况基本一致。

2.4 运动特征分析

通过观察速度−时间曲线（图6），整体上坡体各部分到达速度峰值的时刻基本一致，均在10 s 到达峰值，各点速度峰值介于5～11.6 m/s，速度峰值最高颗粒是位于滑体前缘的3 号颗粒，峰值11.6 m/s。速度峰值最低颗粒是位于滑体后缘的6 号颗粒，峰值约7 m/s。6 个颗粒速度时程曲线整体上变化趋势基本一致。