地震作用下黄土滑坡离散元数值模拟研究

关晓迪

(西安理工大学土木建筑工程学院，陕西 西安 710048)

黄土因其形成环境的特殊性，表现出大孔隙、水敏性、弱胶结、力学强度低、湿陷性强等特点，使得黄土地区滑坡、崩塌等地质灾害频发，形成机理复杂，严重制约了黄土地区的社会经济发展。据调查，自20世纪60年代以来，在陕西省范围内发生的滑坡就多达1 131处[1]，甘肃东部新老滑坡有4万余处[2]。作为黄土滑坡代表，甘肃黑方台滑坡一直备受关注，对于黑方台黄土滑坡的形成机理、滑坡类型、分布规律、监测预警及时间预测预报等方面，学者们做了大量的研究[3-6]，为该滑坡周围的人们生产生活及区域规划等提供理论依据和指导。

数值技术克服了传统方法不能对几何形态复杂的边坡进行计算，不能体现材料的各向异性等缺陷，可更加全面地考虑影响滑坡运动破坏的各种因素。近些年来由于计算机技术的成熟及软件的不断完善，数值模拟方法在滑坡领域得到了广泛应用，模拟效果逐步趋于实际，已成为滑坡运动过程分析中的重要手段。胡炜等[7]以黑方台焦家崖头黄土滑坡为研究对象，采用有限差分数值模型模拟了滑坡的运动过程，证明了此类灌溉渗透型黄土滑坡具有高速远程滑动的特征。贾俊等[8]采用离散元数值方法，对黑方台地区饱水黄土软弱层诱发滑坡的形成过程和滑体的运动学特征进行了研究，分析了滑体高速运动及稳定堆积过程中的相关特征。王念秦等[9]采用PFC2D软件模拟了泾阳南塬滑坡滑体的冲击过程，分析了滑体的刮铲堆积及其速度场的变化情况。通过扫描岩土体颗粒轮廓进行建模以及FLAC3D耦合PFC3D的真实三维地形建模已经成为了多数学者常用的研究方法[10-12]，滑坡模拟的研究方向也扩展到了滑体能量分析、滑带摩擦热研究等更多领域[13-16]。对于边坡地震动力分析的离散元模拟工作，石崇等[17-18]、李祥龙等[19-20]、胡训健等[21]、卞康等[22]学者做了较多研究。施凤根[23]运用PFC3D对文家沟滑坡的高速远程运动特征进行了研究，通过数值分析与现场调查结果进行对比，揭示了文家沟滑坡的运动机理。

当下，尽管国内外专家学者对黄土滑坡的破坏过程、运动机理等方面做了大量的工作，取得了丰富的研究成果，但仍有一些问题需要进一步研究，如地震荷载作用下黄土滑坡的运动滑移特征等。本文主要以甘肃黑方台滑坡为例，对不同地震荷载下以及滑源区不同位置颗粒的运动特征进行分析，找出地震荷载作用下滑源区颗粒的运动规律，为黄土地区地震滑坡的防治和规划提供一定的参考。

1 模型建立

黑方台位于中国甘肃省中部的临夏永靖县，是祁吕贺山字型构造与陇西旋卷构造的复合部位，区内新构造运动强烈，以差异性上升为主要特征，由于受黄河侵蚀，形成了多级阶地地貌。黑方台属黄河Ⅳ级基座阶地。黑方台发育的滑坡多达77处，结合黑方台滑坡所处的地质环境条件，该地区的滑坡大致分为6个区段，分别为新源端、党川段、黄茨段、焦家崖段、焦家段和陈家段。本文以2019年2月发生的甘肃黑方台焦家6号滑坡为例展开研究，图1为焦家6号滑坡的无人机影像，根据AA′剖面做出滑坡的纵剖面见图2。然后利用PFC2D程序，建立滑坡的二维模型(图3)，模型底部宽540 m，高140 m，模型中共包含颗粒8 005个，分为滑源区颗粒(黄色)和非滑源区颗粒(蓝色)。由于滑坡处于黄土地区，黄土颗粒之间的接触与PFC中自带的接触黏结接触模型较为接近，故而选取接触黏结模型作为本次模拟的接触模型。通过试错法以及大量的数值试验来对细观参数进行修改，最终得到模型的模拟参数，见表1。

图2 焦家6号滑坡纵剖面(m)

图3 焦家6号滑坡PFC2D模型

表1 模拟参数

2 模型边界以及动荷载设置

自然界中的边坡均为半无限介质，但受限于计算负荷，数值模型的计算区域只能是有限的，因此对滑坡模型的底部、左侧和右侧的边界颗粒进行设定，对边界颗粒固定位移，图4所示，在滑源区对12个点的颗粒进行监测，从左向右、从上向下依次编号为1—12号，监测颗粒的位移、速度以及运行轨迹，其中1、2、5、8、11、12号点为近坡面点，其中1号点位于滑坡后缘，11、12号点位于滑坡前缘，12号点相比于11号点位置更靠前，滑源区中部有三列点，每一列又有3个监测点，每一列点的x坐标是一致的，主要目的是监测滑源区颗粒随深度增加，对地震动荷载的响应规律。绿色颗粒为边界颗粒，可传递地震荷载，地震荷载为汶川地震时的加速度时程曲线，峰值加速度1.0g，振动持时50 s，见图5。

图4 模型边界以及监测颗粒设置

图5 地震波加速度时程曲线

3 初始加速度下模拟结果

3.1 滑坡运动过程力链分布

滑坡滑动过程主要是滑源区颗粒的运动，也体现在滑源区颗粒间接触的破坏。图6所示，不同时刻滑坡的运动状态，0 s时，滑坡处于稳定状态；30 s时，滑坡已经开始滑动，在地震波的影响下，滑源区颗粒经过流通区，向下滑移并开始堆积，可以看到，在滑坡滑动过程中，滑源区、流通区和堆积区都有颗粒，60 s时，滑坡运动已经结束，可以看到，此时流通区基本上没有颗粒，颗粒主要在堆积区，滑源区也有一些残积土，但是残留在滑源区的土体从力链结构分布上看，较为松散，堆积区的土体也重新形成新的结构连接，但是力链结构也较为松散。

a)0 s时滑坡处于稳定状态

3.2 滑源区颗粒位移分析

通过对滑源区不同位置颗粒位移的分析可以获得整体滑源区颗粒的位移情况，见图7。由图7a可以看出，滑源区最外层颗粒的位移也是有很大差异性的，12号点颗粒，即滑坡体最前缘颗粒的位移量是最小的，位移量基本为0；而1号点颗粒作为滑坡体后缘颗粒，其位移量也是非常小的。可以看出，在滑坡滑动过程中并非是整个滑坡体整体下滑，滑坡体前缘和后缘都会有部分土体存在残留的情况。而11号点颗粒，作为滑坡体前缘颗粒，在整个滑坡体的推动下，位移量最大为392 m，其次是位于滑坡体中部位置的5号点颗粒位移量为345 m，两者相差13.6%，而且在最开始的滑坡过程中，5号点颗粒的位移量要高于11号点颗粒的位移量，说明在地震荷载加载过程中，滑坡体中部位置的颗粒最先滑动，随即带动滑坡前缘的颗粒进行滑动，进而整个滑坡体开始滑动。由图7b可以看出，位于滑坡体中部位置的5、6、7号点的初始位移是不一样的，说明滑源区不同深度的颗粒对地震荷载的响应程度也不一样，近滑坡面的颗粒对于地震荷载的响应程度更大，而且随着地震持时的增加，7号点的位移量基本不变，6号点颗粒的位移量是缓慢增加，然后逐渐平稳，而5号点颗粒的位移量是急剧增加，然后在30 s后开始逐渐趋于平稳。

a)滑源区近坡面颗粒

3.3 滑源区颗粒速度分析

对滑源区不同位置颗粒速度进行分析可以看到整个滑坡体在滑动过程中各个位置颗粒的运动状态，见图8。由图8a可以看出，12号点颗粒(滑坡体最前缘的颗粒)以及1号点颗粒(滑坡体最后缘的颗粒)的速度都是较小的，在滑坡运动过程中，速度也变化不大，而其他近坡面位置颗粒的速度较大，而且速度变化也比较大，和位移变化规律类似，11号点颗粒的速度变化幅度是最大的，速度峰值达到28 m/s，其次是5、8号点。由图8b可以看出，滑源区不同深度的颗粒速度时程曲线差异性较大，7号点的速度时程曲线在0附近波动，而6号点的速度变化不大，其速度峰值为6.5 m/s，但是要大于7号点颗粒的速度，而位于近坡面位置的5号点颗粒速度最大为25 m/s，其速度峰值是6号点速度峰值的3.8倍，表明随着深度的增加，颗粒所受到的限制和约束越多，进而导致其速度变化越小，而近坡面的颗粒由于所受的限制和约束较少，速度波动较大。

a)滑源区近坡面颗粒

4 结论

采用离散单元法对地震动荷载下黄土滑坡进行数值模拟，分析地震荷载作用下滑源区黄土运动特征，主要结论如下：①地震荷载作用下的黄土滑坡存在滑源区、流通区以及堆积区3个部分，滑源区仍然会有一定黄土滑坡的残积土存在，但是结构相较于未滑动前更加松散；②滑源区颗粒的位移差异性较大，处于滑源区前端的颗粒的位移量最大，随着深度的增加，颗粒的位移量逐渐减小，且差异性较大。