基于PFC2D的贵阳西二环滑坡破坏过程模拟

■ 贵州大学资源与环境工程学院 周瑜

贵州同力岩土工程咨询有限公司 郑家金 熊赟

近年来，随着计算机技术的发展，数值模拟方法在滑坡领域得到广泛应用。目前的数值模拟方法中，主要包括确定性分析方法和非确定性分析方法两类，而确定性分析方法又可分为连续介质分析与非连续介质分析方法。其中，连续介质数值分析方法有有限单元法、边界元法、有限差分法等，非连续介质分析方法有块体离散元法、颗粒离散元法、关键块体理论、不连续变形分析（DDA法）等[1]-[2]。颗粒流法（PFC）作为离散单元法的一种，是基于离散元的方法来模拟圆盘或球颗粒介质的运动及其相互作用，可以真实地模拟滑坡的失稳破坏过程。目前已广泛应用于滑坡破坏运动分析之中[3]-[5]。

本文以贵州省贵阳市西二环滑坡为研究对象，采用PFC2D软件模拟了滑坡从变形失稳到运动终止的动态过程，对滑坡不同关键部位颗粒进行位移、速度监测，综合阐明其破坏运动特征，分析其形成演化机制。研究成果可对贵州地区类似滑坡的形成机制提供一定的参考和借鉴，并为工程决策提供依据。

1.工程概况

滑坡部位于贵阳市西二环北段大关立交桥南1000m，路堑工点位于西二环路线左侧，西部苗香对面，此工点处路线走向基本为南北向，地貌类型为低山丘陵。根据现场的地质调查，滑坡区内地层主要有第四系耕植土（Q4pd）、残坡积层（Q4del）滑坡堆积体和三叠系安顺组组（Tla）白云岩构成，地层呈单斜构造，岩层产状140°∠ 25°，为顺向边坡，岩体受节理裂隙切割、破坏，完整性较差。

滑坡底部高程为1270.10m～1272.16m，滑坡顶部高程为1344.79m。滑坡距离路面垂直高度最大约74.5m，滑坡长度约200m、宽约70m、体积约20.0×104m3。滑坡主滑方向为138°，平面上呈现喇叭状形态。滑坡初期失稳后，现状平面地形呈一不规则的中间鼓肚的圆弧状，前翘后躺。滑坡周界明显，滑动面具有明显的滑动擦痕，滑床与滑坡体有明显辨识度。滑体为第四系含碎石耕植土及强风化、中风化白云岩、角砾状白云岩。滑面位于中风化白云岩中的软弱夹层，滑带土厚度为10mm～35mm。滑动带上可以清晰地看到磨光面和擦痕，局部有明显的扰动和拖曳褶皱痕迹，实地调查滑动面擦痕明显。滑床为灰色中风化白云岩、角砾状白云岩，薄层状，节理裂隙发育，层面泥质胶结，受构造切割影响，岩体较破碎、呈块状，节理裂隙较发育，其表层擦痕的方向为138°，擦痕方向与主滑方向一致，滑面整体平直（图1、图2）。

图1 西二环滑坡全貌图

图2 滑坡1-1′剖面图

2.计算模型的建立

本文根据滑坡剖面数据，采用PFC中ball-wall方法建立滑坡滑前模型。该滑坡模型长304m、高111m，其中，滑源区长192m、高68m，颗粒最小半径0.5m、最大半径1m，颗粒总数为1144个。此次模拟在滑坡前缘、滑坡中部、滑坡后缘共设置了3个监测点，分别用于监控滑坡过程中滑体的速度、位移的变化（图3）。

图3 滑坡数值计算模型

3.细观力学参数的确定

PFC滑坡模拟参数主要包括本征参数和接触参数，本文主要选取现场反演的方法确定颗粒单元的接触参数，通过现场参数反演，对所需微观参数进行赋值使颗粒流的运动状态与实际状态基本一致，可认为参数为合理值[6]。具体取值见表1。

表1 滑坡岩土体颗粒细观参数

4.滑坡破坏过程数值模拟

4.1 变形破坏特征

模拟计算到1000时步时（图4a、图5a），坡体中前部颗粒的位移比速度较中后部颗粒大，边坡整体上位移和速度较小，这是由于滑坡区前缘道路建设切脚开挖形成临空面，开挖坡向与岩层倾向一致，在重力荷载作用下，致使坡体向坡前临空方向发生剪切蠕变坡体出现变形；到达10000时步时（图4b、图5b），滑坡前缘产生剪切破坏，坡脚颗粒开始剪出，此时，滑坡滑动带已基本贯通，在坡体压力作用下滑坡开始滑移；到达20000时步时（图4c、图5c），岩体在自重应力的作用下上覆硬岩沿下伏软弱层面向临空面方向蠕动滑移，随着位移量的不断增加，滑移体逐渐被拉裂解体破坏；到达30000时步时（图4d、图5d），中部岩体完全被拉断，中前部和中后部滑体呈现出明显位移差和速度差；到达40000时步时（图4e、图5e），中前部和中后部滑体位移差继续增大，滑体前缘岩体向前滑移至公路，与路面发生碰撞发生解体，后部滑体速度逐渐降为0；50000时步时（图4f、图5f）滑坡基本停止运动并堆积于坡脚平缓处，滑坡整体最大位移为50m。

图4 不同迭代步数位移云图

图5 不同迭代步数速度云图

综上可知，该滑坡破坏模式主要是因为滑坡区前缘由于道路建设切脚开挖形成临空面，在重力荷载作用下，致使坡体向坡前临空方向发生剪切蠕变，边坡的坡体前缘出现了明显的位移变化与应力集中，说明边坡的坡体后缘和上部主要承受拉应力作用，容易使边坡产生张拉裂缝。坡体中存在倾外的弱面，对于层状结构顺向坡，弱面通常为软弱夹层，一般条件下滑体可沿弱面发生蠕滑；在强降雨条件下大量雨水入渗坡内，滑体容增大，重力势能增加，滑带土抗剪强度减小，将沿弱面产生进一步滑移变形，随着位移量的不断增加，滑移体逐渐被拉裂解体破坏，形成滑移－拉裂型滑坡。模拟结果与实际变形情况基本一致。

4.2 运动特征分析

通过观察坡体各部分位移变化曲线（图6）可知，颗粒曲线在20000时步之内变化基本一致，在20000时步曲线斜率开始增长，曲线由下凹过渡到上凸，滑体进入加速阶段，监测颗粒均在25000时步到42000时步之间达到位移最大值，其中，1号颗粒在25000时步附近达到位移峰值15m，2号颗粒在37000时步附近达到位移峰值35m，3号颗粒在42000时步附近达到位移峰值45m。3个监测点位移大小顺序为3号＞2号＞1号，滑坡前部位移大于后部位移。

图6 位移时程曲线图

通过观察坡体各部分速度变化曲线，如图7所示，不同位置监测点达到速度峰值的时刻明显不同，1号颗粒在20000时步附近率先达到速度峰值，2号颗粒在21000时步附近达到速度峰值，3号颗粒在31000时步附近达到速度峰值。各点速度峰值介于6m/s～8m/s之间。滑动初期，20000时步之前，3个颗粒速度曲线整体变化趋势一致，在20000时步后，1号颗粒到达峰值后速度迅速减小，在25000时步附近减小为0。2号颗粒在21000时步到达峰值后，速度呈现逐渐减小的趋势，在38000时步附近减小为0。运行到42000时步3号颗粒速度趋于0，滑坡运动逐渐静止。

图7 速度时程曲线

5.结论

贵阳市西二环滑坡长度约200m、宽约70m、体积约20.0×104m3，滑坡主滑方向为138°，滑面位于白云岩中间软弱夹层，为中型顺层滑坡。通过利用颗粒流方法对西二环滑坡失稳破坏过程的数值模拟分析，得出该滑坡变形破坏的力学机制表现为滑移—拉裂。该滑坡整体位移约40m～50m，滑坡前部位移大于后部位移，颗粒最大位移50m，位于滑体前缘；颗粒最小位移15m，位于滑体后缘；整体速度峰值介于6m/s～8m/s之间，速度和位移变化总趋势为初始高程越低，达到的最大速度和位移越大。

模拟结果表明：降雨为西二环滑坡的直接诱发因素，斜坡变形破坏模式为滑移—拉裂。滑坡滑动最高时速8m/s，最大滑移50m。模拟结果与实际情况基本吻合。由此可见，利用颗粒流法对滑坡的破坏运动过程模拟具有较好的适用性，可为工程决策提供依据。