基于离散元循环荷载作用下的边坡稳定性分析

长期以来，边坡稳定性是土力学中的重要课题之一，由地震引发的边坡失稳在全国乃至全球范围内并不少见。针对这一工程问题，众多学者从试验、理论及数值计算等多个角度开展了相关研究。20 世纪90 年代中后期，随着计算机技术的飞速发展，数值计算方法在边坡稳定性研究中发挥了重要作用，很大程度上降低了地震模拟的难度，众多学者先后应用ANSYS，FLAC，Geostudio，PFC等数值软件进行边坡稳定性研究，取得了丰硕的成果。郑颖人等利用ANSYS 建立了边坡模型，模拟了边坡开挖和加固工程。王伟兴等结合工程实例，研究了应力场与渗流场耦合作用下的边坡失稳，并基于Geostudio 分析了多种水力梯度下的边坡安全系数，运用FLAC 分析了抬升塑性区的变化。王培涛等结合工程案例，基于颗粒流离散元法研究了边坡开挖扰动前后的变形行为，并计算了相应的安全系数。丁辰等以层状土边坡为研究对象，通过模型试验与Geostudio 模拟结果的相互验证，提出了强降雨条件下层状土边坡的失稳机制。周健等引入了强度折减法和重力增加法，运用颗粒流离散元软件分析了边坡稳定性，并对比了有限元法和条分法的计算结果，验证了颗粒流计算边坡稳定性的有效性。J.B.Wei 等通过无人机航拍技术获取了三维地形数据，并借助PFC 模拟了滑坡启动，监测出的颗粒位移及速度与边坡实体基本一致，该项研究证明PFC 可作为模拟边坡的有效手段。C.Shi 等运用PFC 建立了软化微力学接触模型，研究了阻尼、应变速率对地震滑坡的影响。

上述研究主要运用颗粒流软件分析了滑坡灾害。目前，大多数边坡稳定性研究往往只针对特定工程，受到了区域地质条件、降水等因素的限制，形成的理论具有一定的局限性，加上不同区域的物源组成差异较大，相应的细观参数设计也存在诸多不同之处，因此很难对具有特殊气候、地理地貌等特征的区域滑坡灾害形成客观认识。鉴于此，本文以贺兰山东麓沿线分布的碎石土竖直边坡为研究对象，运用颗粒流模拟软件建立边坡物理模型，研究动力条件下的边坡稳定性，从细观上揭示滑坡灾害的演化特征，为研究区及相似地质条件的区域滑坡、泥石流灾害防治提供参考。

贺兰山位于宁夏回族自治区与内蒙古自治区交界处，绵延200 余公里，主峰高3 506 m。贺兰山既是宁夏回族自治区的行政边界，也是银川市发展成为“塞上江南”的天然屏障。受地质构造、气候条件的影响，贺兰山东麓沿线分布了数十条泥石流冲沟与历史泥石流堆积扇。堆积扇受山洪冲刷和人类活动的影响，形成了典型的垂直边坡。泥石流发育时有流体侵蚀，可形成补给物源，同时在地震荷载作用下，堆积物原有的沉积结构被破坏，形成松散的泥石流物源，极易引发泥石流灾害，为贺兰山东麓沿线道路、古文化遗址、种植园区埋下安全隐患，也给疏浚工程、山洪治理带来极大的不便。典型的碎石土竖直边坡如图1 所示，从揭露的剖面来看，坡体内部粗、细颗粒混杂，表层土体发生了崩塌，堆积在坡脚处。

研究区地处鄂尔多斯板块边缘位置，属于强震带。贺兰山地震带如图2 所示。据资料记载，宁夏区内先后发生过多次地震，引发的滑坡、泥石流等次生灾害曾造成了重大的经济损失。

本文运用PFC 2D 软件生成边坡物理模型，并使用wall 命令生成模型边界。模型宽为25 m，高为10 m，坡面坡度为90°，坡高为5 m，共由40 706 个颗粒组成。本文参照现场取样的级配曲线，对颗粒尺寸进行了相似设计；并根据前期研究成果，设定研究区内物源粒径大于60 mm 的颗粒占7.61%，粒径大于2 mm 的颗粒占65.31%。根据《建筑地基基础设计规范》（GB 50007—2011），可将该类型岩土定名为碎石土，属级配不良土。

由于模型中颗粒的尺寸相差较大，颗粒间会出现应力集中现象和边界效应，导致计算机无法运算，这也是运用颗粒流软件进行模拟时出现的难题之一。由相关研究可知，当岩土材料的内尺寸比小于0.01 时，边界效应会显著降低，同时本文考虑研究区内碎石土中粗颗粒的质量分数较高，其不但是土体的主要组成部分，也是灾变演化中主要的受力体，因此模型颗粒采用了2 个尺寸范围（主体部分由3～5 cm 的颗粒组成，选用2～3 cm 的颗粒进行填充），并且在坡脚位置设置了宽为1.5 m，深为0.5 m的粗化层（用于模拟沟床上覆粗化层的压实作用），其粒径为4～5 cm，模型颗粒分布图如图3 所示。为了体现边坡的脆性破坏，同时考虑到细颗粒的黏性，本文采用了接触粘结模型；并且项目组参考了颗粒组成相似的碎石土案例，将法向刚度和切向刚度分别设定为2.5×107，1.0×107 N/m。

本文将抗震烈度设为8 级，相应的水平峰值加速度为0.2～0.3 g，竖向峰值加速度为0.1～0.17 g。考虑到竖向加速度作用下的边坡滑动效果并不明显，因此本文将模型的地震峰值加速度均设置为0.2 g，荷载作用时间设定为20 s，并且对边坡底部墙体分别施加了水平向和竖直向速度。峰值加速度变化曲线如图4 所示。

本文运用ball distribute 命令分区生成颗粒，颗粒间有一定程度的重叠，因此需要通过循环求解命令，求解至颗粒间的相互作用力达到平衡，使颗粒均匀且无重叠地分布在模型内，平衡后的竖直边坡颗粒分布状态如图5 所示。天然边坡在自重作用下可以进行压实固结，为了模拟这一行为，该模型为平衡后的颗粒设置了重力加速度（9.8 m/s²），待土体固结完成后施加动力荷载。

本文设置了2 种方向的地震波来监测坡面颗粒位移。坡体加载后，到达稳定状态时的位移分布如图6 所示。由图6 可知：颗粒位移较为接近，最大颗粒位移达到5 m；颗粒滑动后，堆积在坡脚处。加载形式不同时，颗粒的移动方式差别较大。竖向加载时，坡面颗粒在临空状态产生的横向位移较小，颗粒滑落后仍有位移产生；横向加载时，坡面颗粒在临空状态产生的横向位移较大，颗粒滑落后产生的位移较小。本文经分析后发现：纵向地震波加载时，颗粒受竖向力作用，产生的水平速度较小；横向地震波加载时，颗粒的初速度、加速度均大于竖向加载时，且受力产生位移的颗粒其范围更大（图6）。

制样完成后，通过平衡求解，可知边坡处于不平衡状态。由荷载作用下坡内不平衡力分布状态图（图7）可知，坡肩、坡面及坡脚在震后的不平衡力最大并向坡内逐步减小，说明震后的坡面颗粒接触力小，颗粒间相互作用有所减弱。横向地震波作用时，不平衡区域相对较大，边界线呈弧状。在横向地震波作用下，粗化层的“增重”效果并不明显，颗粒间接触力整体减小；而在纵向地震波作用下，坡内存在1 个小的平衡区域，颗粒间接触力较大，滑坡范围相对减小。

图8a 给出了1，2，3 号监测点的竖向应力值。随着地震荷载持续作用，监测点应力值快速降低，不同监测点应力变化速度基本一致，但最终的目标值略有差异，表现出埋深越小，应力目标值越小，应力降低幅度越大的特点。

图8b 给出了2，4，5 号监测点的竖向应力变化规律。3 个监测点的应力减小至某一值后即保持恒定，应力降低的速度和监测点到坡面的距离成正比，说明图8a 中3 个监测点的应力几乎同步降低至目标值。同时，表层颗粒在地震荷载作用下会发生位移，导致坡面形态发生变化。结合图10，可知沿着垂直坡面方向，应力由外向内逐步增大，4，5 号监测点在坡体重新稳定后，两者相对坡面的距离差距不大，因此应力降低的幅度无明显差异（图8b）。1，2，3 号监测点在坡体重新稳定后，其中3 号监测点的埋深最大，所以应力变化幅度最小。

图8c 给出了5 个监测点的水平应力变化规律，该图与图8a、图8b 表现出相同的变化规律，距离坡面较近的监测点其应力变化快于距离坡面较远的监测点。4，5 号监测点的应力变化幅度大于3 号监测点，则是由于坡体重新达到稳定状态后，4，5 号监测点到坡面的距离小于3 号监测点。

图9 给出了横向地震荷载作用下5 个监测点的竖向、水平向应力变化特征，其与纵向地震波作用下的应力演化特征基本一致。在2 种地震波作用下，3 号监测点均表现出应力先减小后增大，最后再减小至某稳定范围内的特点。本文认为：在加载初期，由于坡面崩塌，3 号监测点的颗粒迁移，离开坡面，应力降至最小；在加载后期，上部崩塌的颗粒在坡脚堆积，3 号监测点附近的上覆压力增大，进而引起内部应力增大；随着持续加载，颗粒不断滑动，趋于新的平衡状态，应力值也逐渐达到恒定值。从应力降低的幅度来看，在2 种动力条件下，坡体大部分区域的应力降低至初始应力的1/10，显然不利于边坡稳定，特别是在震后降雨、山洪暴发期易引发次生灾害。

对比图8 和图9，可知对于不同的动力荷载，碎石土竖直边坡的响应机制明显不同，其中1，2，3 号监测点的竖向应力和水平应力在2 种地震波作用下无明显差异，这是因为3 个监测点距离临空的竖直坡面较近，颗粒间的相互作用力小于4，5 号监测点，所以在地震波作用下其应力衰减较快。而4，5 号监测点的竖向应力值在横向地震波作用下的衰减速度要小于在纵向地震波作用下的衰减速度，5 号监测点水平应力的变化规律也是如此。上述现象说明纵向地震波作用下的碎石土内部应力的释放要快于横向地震波作用下，在地震波作用初期该现象易引发边坡失稳。

图10 给出了地震荷载作用结束时坡体内力链接触分布特征，力链的形状与颗粒间接触应力的大小有关，力链越细，表明接触力越小。图10 与图8、图9 给出的监测点应力分布特征较为一致，在坡面两侧的坡肩、粗化层区域，接触力很小，说明颗粒间的作用减小，颗粒处于相对松散的状态，在汛期颗粒易随水流迁移，形成较大的冲击力，从而影响沿途边坡稳定性和基础设施安全。

本文利用PFC 2D 软件分别模拟了纵向和横向地震荷载作用下的边坡失稳过程，并从颗粒位移、监测点应力等多个因素的动态变化过程入手，揭示了滑坡体颗粒的运动形式和滑坡机制，得出以下结论。

（1）在横向与纵向地震波动力条件下，碎石土边坡均发生了滑坡，坡面颗粒最大位移达到5 m，坡顶高度随之减小。在2 种地震波作用下，颗粒的运动方式略有差别，其中横向地震波作用时颗粒在临空状态下产生的横向位移较纵向地震波大。

（2）距离临空面越近，应力消散越快；埋深越大，应力降低的幅度越小，且埋深大的颗粒在持续加载过程中因颗粒堆积，会出现短暂的应力增大。

（3）从应力变化幅度来看，在横向和纵向地震波作用下，坡体上部、坡面等区域的应力在震后会下降至初始应力的1/10，显然不利于灾害防治。纵向地震波作用时，其滑动范围小于横向地震波，但应力释放速度要高于横向地震波。