降雨作用下的土质边坡变形破坏颗粒流仿真模拟

李龙起，赵瑞志，王 滔，赵皓璆，王梦云

(成都理工大学地质灾害防治与地质环境保护国家重点实验室，成都 610059)

滑坡是指斜坡岩土体在各种内外营力作用下沿着贯通破裂面发生的以水平位移为主的滑动现象，具有突发性和极强的破坏性。同时，由于坡体本身具有不同的结构和复杂的物质组成情况，使滑坡体的变形破坏模式多种多样。特别是中国幅员辽阔，地形复杂，经常发生滑坡灾害，例如四川宣汉天台乡滑坡[1]、重庆武隆县鸡尾山滑坡[2]、四川茂县新磨村滑坡[3]等均造成人员及财产损失。

自然界多数滑坡的发生都与降雨密切相关，因此，众多学者通过不同维度对降雨型滑坡进行了不断的探索与总结[4-6]。对于滑坡的研究，一般是通过物理试验、数值模拟并结合理论分析来进行，其中物理试验通常会受到边界效应、尺寸效应等因素的影响，且效率低、工作量大；而数值模拟凭借高效率、低成本和可重复性好的优点应用越来越广泛。前人对数值模拟在滑坡中的应用进行了不断地探索，但大多都是基于有限元的基础上进行[7-8]，而滑坡体具有不均匀性、非连续性和离散性的特点，通过有限元方法在求解位移不连续及大变形问题上存在不足。Cundall[9]提出离散元分析方法为滑坡数值模拟提供了新途径，并由此产生了颗料流(PFC)、离散单元法程序(UDEC)、三维离散单元法程序(3DEC)等离散元软件。PFC2D相对于其他离散元软件有效率高和可以有效模拟大变形的优点，因此得到了学者广泛关注。毕钰璋等[10]应用PFC2D离散元软件对牛眠沟高速远程滑坡全过程进行模拟，得到不同摩擦系数下的滑坡变形及堆积特征；陈达等[11]运用PFC2D以陕西省洛南县刘涧滑坡为研究对象,得到不同位置特征点的加速度随时间的变化情况；周健等[12]运用颗粒流PFC2D分别对黏性边坡和砂性边坡进行数值模拟，结果表明土的性质对边坡的破坏形式有很大的影响，随着颗粒黏性的减小边坡的破坏形式从脆性向塑性转化；赵洲等[13]运用PFC2D对杨家湾堆积层滑坡的破坏过程进行数值模拟，得到堆积层滑坡的破坏模式是一种典型的牵引式渐进破坏过程。目前，学者们运用颗粒流对滑坡进行的研究主要集中在滑坡的破坏过程及各种因素对堆积形式的影响，并没有对处于时效变形阶段的坡体在降雨或其他工况下的发展趋势进行研究。

以四川省古蔺县的竹林沟滑坡为例，通过现场勘查、理论分析和数值模拟对滑体进行综合研究。运用颗粒流离散元软件PFC2D对滑动体在持续降雨作用下的变形破坏过程进行仿真模拟，并通过在模型上设置监测点，得到坡内不同位置的变形特征，由此分析预测滑坡的发展趋势，为滑坡防治提供科学依据。

1 竹林沟滑坡状况

1.1 竹林沟滑坡地质概况

图1 竹林沟滑坡纵剖面

古蔺县位于亚热带季风气候区，云雨天气较多，月平均降雨量如图2所示。据已有资料[14]显示，在2000年一次降雨中，滑坡后缘一户居民房出现裂缝，裂缝走向大致与主滑方向垂直，并于2007年一次连续两天降雨中(降雨量约200 mm)墙体倒塌，同时另外十几户居民房出现不同程度的开裂(图3)。2008年汶川地震后，当地政府组织人员对滑坡进行过安全评估，结果显示滑体在地震情况下并未出现大的变形。2014年当地政府对滑坡进行治理，在变形较大部位布设抗滑桩(桩埋深7 m)，但据当地村民介绍，坡体在连续降雨或大暴雨情况下仍在变形，持续至今。由此可见降雨是坡体产生变形的主要原因。滑坡已经过多年变形，累积变形较大，存在重大隐患，对其在降雨作用下的稳定性研究显得尤为必要。

图2 古蔺县月平均降雨量

图3 现场房屋受损情况

1.2 滑坡滑动原因分析

竹林沟滑坡的成因与众多因素有关，滑坡全貌如图 4 所示。根据已有资料[14]及现场勘查分析，竹林沟滑坡的成因主要有以下3个方面。

图4 竹林沟滑坡全貌

(1)当地居民在坡体上修建房屋、开垦荒地和开展其他工程活动，造成坡体表面土质疏松裸露，为降雨条件下的雨水下渗创造了有利条件。雨水下渗的过程中使得坡体饱水加载，导致坡体下滑力增加而抗剪强度降低。

(2)滑坡体整体坡度较大(其中中后缘坡度为15°～30°，中部坡度为15°～40°，中前缘坡度相对较缓，为10°～20°)，自重应力沿滑动方向的分力较大，对坡体稳定性不利。坡脚处由于河谷下切，使得前缘临空，缺乏有效的支挡结构。

(3)由于人类活动频繁，生活排水、稻田灌溉及工程活动等均对滑体产生不利影响。而在坡体中后缘由于道路修建进行的大面积的开挖，破坏了坡体的整体性，同样也加剧了坡体的破坏进程。

2 微观参数选取过程

2.1 颗粒流离散元理论简介

(1)

(2)

式中：kn和ks分别为接触点法向和切向刚度；un和us分别为接触点法向和切向的微小位移。

(3)

式(3)中：c为颗粒接触点的黏聚力；φ为颗粒接触点的摩擦角。

颗粒在整个计算过程中，循环运用牛顿运动定律(Newton’s laws of motion)和胡克定律(Hooke’s law)，根据牛顿运动定律不断变换颗粒之间和颗粒与墙体之间的接触关系，根据胡克定律确定颗粒新的接触关系中的接触力(图5)。当颗粒之间或颗粒与墙体之间的作用力大于黏结强度时，颗粒就会发生相对错动。

图5 颗粒计算循环过程

2.2 模型微观参数的选取

离散元 PFC2D采用的是微观的力学参数，通过建立正确的微观和宏观之间的力学关系，把滑坡体宏观现象用微观参数表达出来。微观参数的选取通常是在 PFC2D中进行反复的虚拟双轴压缩试验，得到应力-应变曲线与室内压缩试验曲线进行匹配，最终得出滑坡体的微观力学参数。由于坡体滑动的主要诱发因素是降雨，因此直接对坡体土样在饱和状态下的强度参数进行双轴试验标定以模拟持续降雨工况，坡体的宏观力学参数如表 1 所示。

表1 土的宏观力学参数

通过 PFC2D软件建立双轴压缩试验的模型，随后在模型内部填充颗粒。考虑到实际颗粒大小具有不均匀性，同时保证计算效率，最终选择粒径介于0.36～0.54 m的颗粒进行填充，孔隙率为0.155。为了避免尺寸效应对试验造成影响，模型的短边尺寸还应大于40倍的颗粒平均半径[16]。各粒径颗粒数量从Rmin到Rmax服从正态分布。

模型建立完成对其施加围压，并通过顶、底板施加轴向压力进行试验。经过大量虚拟双轴压缩试验得到多条应力-应变曲线，与被模拟材料饱和状态下的室内试验所得结果匹配，最终得出匹配程度最高的应力-应变曲线如图6所示。对应微观参数的取值如表2所示。

表2 双轴数值模型细观参数

图6 双轴压缩应力应变曲线

图7为模型赋予参数后在不同围压下(50、100、150 kPa)模拟双轴压缩实验所得出的莫尔破坏包络线，最终反算出的宏观参数内摩擦角φ为16.3°，黏聚力c为21.7 kPa，得到的结果与被模拟滑坡土体在持续降雨工况下的c、φ相近，这也论证了所选微观参数数值的可行性。

图7 莫尔破坏包络线

3 数值模拟结果分析

3.1 模型的建立

PFC 提供两种建模方法，ball-ball 法和 ball-wall 法。ball-ball 法是滑坡体和基岩均采用颗粒进行填充，因此所用颗粒较多，计算过程较长，适用于滑动面未知的情况。ball-wall 法是滑坡体采用颗粒填充而基岩采用墙体进行模拟，所用颗粒较少，运算效率高，适用于滑动面已知的情况。根据文献[14]和现场勘查已确定滑坡潜在滑动面的位置，模型采用 ball-wall 法进行建模。具体建模过程如下。

(1)在 AutoCAD 中绘制滑体和滑床的边界线，尺寸与被模拟滑坡真实尺寸一致，将图层存储为 DXF 文件格式。

(2)在 PFC 中使用 geometry import 命令和wall import geometry 命令将 DXF 文件导入后根据滑体和滑床边界线生成边界墙。

(3)在边界墙形成的空间内填入足够多的 Brick，并根据表 2 赋予颗粒接触刚度、摩擦系数等微观参数。

(4)删除超出边界线以外的颗粒，进行多次迭代计算使得模型达到平衡，并在坡体上设置需要的监测区域。

最终所建模型如图 8 所示，模型运用 PFC 软件自带的块体填充法(brick)在墙体内部填充足够多的颗粒[17]，该方法规避了膨胀法由于膨胀系数难以准确选择，而导致颗粒接触处应力过大，颗粒飞出的现象，且操作简单。

图8 滑坡ball-wall模型及测点分布

模型建立完成后，将模型的初始位移和速度置零并去除边界约束墙，对坡体施加重力使坡体在重力作用下开始缓慢变形。

3.2 滑坡破坏过程分析

滑坡体不同时步的位移云图如图9所示。模型坡体边界线只作为坡体变形的参考。

图9 位移云图

由图9(a)可知，4 000 时步时坡体并未出现较大变形,仅在中部一区域内产生小变形，说明滑坡的变形是从中部开始的。坡体上、下两部分覆盖层厚度及滑床坡度均存在差异，导致 12 000 时步时滑体表现出二级滑动特征，一级滑动面高程约 570 m，二级滑动面高程约 680 m。一级滑动发生后，坡体中部产生拉张裂缝，为二级滑动创造了临空条件。20 000 时步时坡脚处有部分土体在上部挤压作用下剪出，坡体中部和后缘拉张裂缝增大，后缘出现“下坐”现象图9(c)，这是因为后缘裂缝增大到一定程度，使得裂缝后土体失去有效支挡从而向前倾倒，伴随着裂缝的缩小和土体的“下坐”现象，这种现象也表明坡体进入累进性破坏阶段[18]，进入该阶段后坡体已不能再承受过多变形，其之后的变化发展也具有不可控性，当潜在剪切破坏面的锁固段被贯通就会发生滑坡。图9(d)为滑坡在 100 000 时步的发展变化情况，此时滑坡剪切面已经被贯通，坡体发生整体滑动，坡脚土体滑出数十米。根据滑坡破坏特征可知，该滑坡的破坏模式是典型的蠕滑-拉裂。

3.3 滑坡运动特征分析

为了能够精确地研究其滑动过程中各部位的发展变化及力学响应情况，在坡体中分别布置 7 个控制测量点及监测区域如图 8 所示。控制点用来监测坡体位移，监测区域用来测量滑动过程的应力变化情况，控制点从右到左序号分别为 1～7，若以坡顶点在水平面上投影为坐标原点，监测点坐标分别为(437.28，61.34)、(362.75，80.83)、(273.95，111.06)、(212.24，137.82)、(141.65，162.89)、(86.67，186.81)、(39.58，198.86)。模型将控制点布置在监测区域中，通过分析同一位置位移及应力的发展趋势对坡体变形破坏特征进行研究。

由于滑坡一旦进入累进性破坏，其稳定性无论在工程措施还是时间上都极难甚至不可能得到控制[19]，因此滑坡防治的最佳时间是在时效变形阶段，所以对时效变形阶段进行研究更具意义。由此，仅对滑坡进入累进性破坏前(20 000 时步前)的位移及应力的变化过程进行研究，以向滑坡治理提供科学依据，对滑坡堆积过程不做详细分析。

图10为各监测点的位移时程曲线。由图10可知，在 5 000 时步之前滑坡各测点位移曲线斜率均较小，且位移的增长基本保持一致，这表明在滑坡浸水饱和后稳定性降低，滑坡发生整体的微小位移。在 5 000 时步后 7 号测点位移曲线的增长速率明显减缓，且之后该点位移量与其余监测点位移差距逐渐增大，这是因为坡体在滑动过程中存在不协调变形[20]，滑坡前部滑动较快，而坡体后缘坡度缓、滑动慢，促使后缘产生拉张裂隙，随着裂隙不断增大，向坡内延伸至基岩层面后，滑坡后缘完全失去了前部的牵引作用，仅在自重作用下滑动，因此 7 号测点区域土体与前部的位移差越来越大。12 000 时步之后，坡体各区域土体位移变化量差异明显增大，表现为从后缘到前缘位移量渐增的特征，且1、2、3号成一组变化，4、5、6号成一组变化，一方面是由于此时坡体出现二级滑动趋势，另一方面是坡体中-后段变形受前部土体抗滑力约束，而中-前段土体发生移动后产生临空面，为坡体后缘土体的失稳提供变形空间，即表现为从下到上的牵引式破坏。在20 000时步时7号测点位移量达到 1.05 m，1～6 号测量圆位移量达到1.65 m左右。

图10 滑坡测量点位移曲线

图11为各测量圆 20 000 时步内的应力曲线图，变形初期，坡体中部位移较前部大，致使 1、2 号测量圆区域应力水平不断积累，在 6 000 时步左右，坡体整体出现较大位移，1、2 号测量圆区域应力得到释放，应力值陡降。其中，1 号测量圆位于坡脚处，坡度较缓，因此X方向应力明显大于Y方向，约为Y方向的1.5倍；2号测量圆所在位置坡度较1号陡，因此X方向应力略大于Y方向，1、2号测量圆的峰值应力均产生在5 500时步。由于滑床坡度整体表现为前缓后陡，因此3、4、5号测量圆X、Y方向的应力差依次缩小，应力大小表现为σ3>σ4>σ5。6号测量圆位于滑坡后缘受拉区附近，受到中部坡体产生的拉应力作用，表现为X方向应力大于Y方向应力。7号测量圆滑床坡度缓，因此测量圆Y方向的应力比X方向大。

图11 各测量圆内应力变化曲线

4 结论

根据现场勘查资料，结合颗粒流 PFC2D技术对竹林沟土质滑坡在持续降雨工况下的发展破坏进行模拟研究，得到结论如下。

(1)坡体在持续降雨工况下的破坏是从坡脚开始的，随后发生从下到上的牵引式渐进破坏。中-前段坡体破坏为中-后段坡体提供了临空条件，因此坡体破坏表现出二级滑动特征，滑坡整体上表现为蠕滑-拉裂破坏模式。

(2)滑坡各区域变形存在明显的不协调，后缘张拉裂缝贯通后，失去前部的牵引作用，仅在自重应力下产生位移。二级滑动也使得中-前部与中-后部呈两组不同发展趋势。

(3)坡体前缘应力和位移均大于后缘，一级滑动发生后坡脚应力得到释放，前缘应力水平有所下降，二级滑动受到前部土体阻挡，应力水平下降较少。