云南黄坪库区滑坡运动及其失稳模式的离散元模拟

蒋 权，陈希良，肖江剑,邱向东,李攀峰

(1.中国科学院宁波材料技术与工程研究所，浙江 宁波 315201；2.中国科学院大学，北京 100049;3.中国电建集团成都勘测设计研究院有限公司，四川 成都 610072)

0 引言

库区滑坡失稳破坏后，滑坡体入水激起巨浪而产生滑坡涌浪，若激起的涌浪足够强还可能会造成坝顶过水甚至冲毁水工建构筑物，堵塞河道、威胁船只航行及沿岸居民生命财产安全[1]。特别是正在蓄水或蓄水不久的库区更容易发生此类地质灾害。一方面，水库的修建会不可避免地改变库区岸坡原有的自然地质条件，如库岸岩土边坡的物理性质，进而影响边坡的稳定。另一方面，库水会对库岸已存在的不稳定地质体和滑坡体、崩塌体产生浸润和托浮作用，其变化往往对库岸坡体的稳定性进一步造成影响，在受到暴雨、库水位突然变化等外界影响因素的作用下，一些处于不太稳定状态的库岸滑坡体就可能失稳[2]。

极端气候异常，工程扰动增强，地质条件复杂等，是地质灾害防治工作需要面对的严峻挑战[3]。针对滑坡体失稳后滑坡体运动的问题，相关研究在许多国家的不同学科和领域都在开展。2002年1月成立的国际滑坡协会及每3年举办一次的世界滑坡论坛更是大力推进了对滑坡运动的研究[4-5]。在分析滑坡运动时，对不同滑坡体物理形态有不一样的方法。就小块且松散堆积的滑坡体而言，许多学者用离散单元法、非连续变形法、颗粒流等非连续介质数值方法对坡体滑动过程进行了模拟分析，并取得不少成果[6-8]。但目前的一些工作通常有很多简化，如将滑坡体简化成简单的块体，或者只研究二维情况下的离散坡体滑动情况，这难以真实刻画三维复杂的滑坡体力学特性，距实际工程应用较远。

本文以溪洛渡水电站黄坪滑坡为例，结合现场地质调查，对滑坡特征进行了详细的介绍，并分析了滑动成因和机制；在不考虑空气与江水作用的情况下，利用离散元对三维离散滑坡体的滑坡失稳过程进行了模拟，并详细的分析了其各阶段的运动特征。

1 滑坡概况及成因机制分析

1.1 滑坡概况

2013年7月27日下午，云南省昭通市永善县黄华镇黄坪村大石包附近发生滑坡。滑坡位于金沙江右岸永善县黄华镇老黄坪集镇上游约1 km大石包附近，即黄码公路桩号K42+000～K42+200 m段。滑塌范围：顺江方向长约200 m，江水位高程554 m以上约110 m，后缘至黄码公路附近高程约660 m，部分区域高于路面50 m范围，滑体厚度3～5 m，滑坡体积数十万立方米。滑坡区位于预测滑塌区范围内的上游部分，现场及断面示意分别见图1、图2。

图1 黄坪滑坡全貌Fig.1 The panorama of the Huangping landslide

图2 黄坪滑坡剖面示意图Fig.2 Cross section of the Huangping Landslide

该岸坡地形较陡，天然坡度在34°～38°。堆积体下部地层岩性为奥陶系上统与志留系O3+S(S)灰色、深灰色、灰黄色砂岩、粉砂岩、泥岩、页岩、砂质页岩夹泥灰岩、泥质灰岩，岩性较软弱，产状N30°～ 45°E/SE∠17°～20°；上部陡坡地形为阳新灰岩；堆积体为崩坡积与冲洪积形成的块碎石土层，块碎石粒径较小，一般为5～20 cm，表部含量较高约40%～50%，下部含量较低约25%～35%，细粒部分主要由细粉砂与黏粒组成，在堆积体表层可见钙质胶结的硬壳。

1.2 成因机制分析

据现场调查和已有的资料分析，黄坪滑坡形成的原因主要有以下两点：(1)坡体物质组成、结构、坡度等自然条件造成该处岸坡天然状态下稳定性较差。岸坡下部地层物质组成决定了其较软弱的岩性；该部位岸坡地形较陡，中间被一条较小的无水干沟切割成上、下游两部分，岸坡岩体结构完整性差；岸坡由于受长期的水流冲刷作用，为表层破碎岩体各类重力地质现象的产生提供了有利条件。(2)水库蓄水是黄坪滑坡的重要诱发因素，加速了该处库岸失稳。滑坡发生前水库已经开始蓄水，岸坡前缘坡脚处于蓄水位以下，由于水的软化弱化作用，前缘坡脚岩土体的强度被进一步降低，同时也导致其有效应力的减小。

根据滑坡概况及成因分析其滑动机制为：(1)库水位的快速上升，导致了中下部坡体的软化，抗剪强度降低；水位线抬升产生的浮托力还减小了岸坡前缘岩土体的有效应力导致抗滑能力降低。(2)在上部坡体的重压下，水下坡体发生突发式破坏，从而导致库水位上下坡体快速失稳，坠入库内。

2 离散元模拟及模型参数

2.1 离散元数值模拟方法

离散元法的基本原理是基于离散化的思想，对研究对象进行单元划分，离散化的各个单元之间相互独立。根据单元之间的相互作用和牛顿运动定律，采用动态松弛法和静态松弛法等迭代方法，对各单元进行循环迭代计算，得出每一个时间步长内各单元的受力及位移。其控制方程主要包括平动与转动两方面，具体形式如下:

(1)

(2)

2.2 数值模型

基于黄坪滑坡的实际数据，以黄坪滑坡为中心选取了1 280 m×200 m×416 m 模拟计算区域，即模型底部边界长1 280 m，两侧边界高为416 m，沿江长200 m。基于滑坡剖面建立了黄坪滑坡所在区域的离散元计算模型，模型分为滑坡体、基岩滑床、滑坡边界三部分。为监测滑坡体随时间的变化，在滑坡体中心取截面A，并在截面A上设置11个监测点：后缘和前缘各分布1个点，滑坡后部、中部以及前部按照不同深度在表层、中层、底层分别设置了监测点。由于模型较大，取右边关键局部示意图见图3。

图3 黄坪滑坡离散元计算模型Fig.3 Discrete element model of the Huangping Landslide

在离散元模拟中，滑体采用离散颗粒，滑床及滑坡边界采用同滑体一样参数的不考虑变形的刚体。对于土质滑坡和破碎岩石滑坡来说，要想得到和实际较为一致的模拟结果，可将滑坡运动看成颗粒流运动，也可将滑坡体等效为拟流体来进行数值模拟[9]。本文用正态分布的不同大小圆球颗粒表征破碎岩土，并组成仿真中所需要的滑坡体，滑坡体体积约为3.0×105m3，组成滑坡体的颗粒总数近4.0×105颗。

离散元法的接触模型有多种，接触力的计算方法也各不相同，但是整体计算的原理都是相同的。本文中，颗粒间的接触力采用Herz-Mindlin无滑移接触模型求解[10]。其中，模型中法向接触力的求解基于Hertzian 接触理论，而切向接触力的求解基于Middlin-Deresiewicz 接触理论[11]。

2.3 计算参数

计算中用到的滑坡体与滑床或滑坡边界的材料属性、相关接触属性及运行参数见表1。

表1 黄坪滑坡涌浪计算参数

3 滑坡运动模拟结果与分析

3.1 滑坡运动过程的整体特征

从滑坡体外部出发：首先，取滑坡体失稳后运动过程中典型时刻的滑坡体表面形态见图4，可知滑坡体表面在运动过程中由凸形态逐渐转为凹形态。

图4 不同时刻黄坪滑坡运动形态Fig.4 Movement states of landslide at different time

然后，从滑坡体上方俯视向下看可见滑坡失稳后的运动过程中，滑坡体后缘边界线趋于凹形态(图5)。

图5 滑坡后缘边界运动Fig.5 Trailing edge movement of landslide

因此滑坡运动过程中，其表面形态会发生转变。而滑坡体后缘边界线趋于凹形态是由于失稳滑坡体两侧受到来自周围相邻未失稳滑体的支撑作用，滑坡体呈现出中间失稳体积较大，两侧失稳体积较小的特征。这与一般滑坡现场所表现出的形式一致。

在滑坡体失稳后的滑动过程中，整个滑坡体的平均速度、最大速度及最小速度如图6所示。

图6 滑坡运动速度Fig.6 The velocity of landslide movement

从图6中平均速度曲线的趋势可知滑坡体运动的大致情况：0～5 s为滑坡启动后滑坡体整体运动的加速阶段，其中加速较快的是后缘及后部受压力与阻力较小的滑坡体，这可从压力分布图7中看出；5～15 s，在中部前部凸形滑坡体颗粒阻力的作用下，滑坡体整体开始减速；15～25 s，在后缘及后部滑体冲击与不断堆积作用下，滑坡体整体再次缓慢加速，并开始滑向坡底；25 s以后，滑坡体在前方坡底堆积的滑坡体颗粒阻力作用下开始缓慢减速运动，直到70 s，滑坡体趋于稳定并停止运动。

从图6中最大速度曲线可以看到：在启动时刻，有局部滑坡体颗粒达到12 m/s的速度；在高速滑动阶段，有局部滑坡体颗粒达到32 m/s的最大速度；高速滑动阶段过后，滑坡体以似整体颗粒流形式流动，滑坡体颗粒的最大速度在缓慢减小。

从图6中最小速度曲线可看到，在滑坡启动、高速滑动及自稳堆积阶段，存在速度较小甚至没有运动的滑坡体颗粒；而在整体流动阶段，滑坡体颗粒几乎都在运动。

图7 A截面上的压力分布Fig.7 Compressive force distribution on cross section A

3.2 滑坡体内部运动特征

为深入的分析滑坡运动特征，取滑坡体内截面A的8个典型运动时刻来对滑坡体内部颗粒分布、速度分布及压力分布做出分析。

3.2.1滑坡体内颗粒分布特征

如图8不同时刻截面A上不同大小滑坡体颗粒分布情况可知：在滑坡运动启动前或初始阶段，滑坡体1大颗粒(灰色)分布在滑坡体底面，滑坡体2小颗粒(橙色)分布在滑坡体表面。随着滑坡运动的进行，颗粒之间碰撞与振动加剧，部分小颗粒慢慢运动到大颗粒滑坡体1内部及底部，大颗粒运动至滑坡体表面，而小颗粒会运动到滑坡体下层。这说明，在滑坡体运动过程中产生了“巴西果效应”，小颗粒通过大颗粒间的缝隙逐渐往下层渗流并填充到大颗粒下部空穴中，而大的滑坡体颗粒则移动到上部并在一定程度上被阻止向下运动[13]。

图8 A截面上颗粒分布Fig.8 Particles distribution on cross section A

3.2.2滑坡体速度分布特征

与将滑坡体看成刚体不同的是，离散滑坡体颗粒的速度在不同位置是不同的，且速度的变化是分阶段进行的。典型时刻截面A上滑坡体颗粒的速度分布见图9。

图9 A截面上的速度分布Fig.9 Velocity distribution on cross section A

从图9中可看到滑坡运动启动-高速滑动-碰撞-滑至坡底-堆积自稳的典型阶段：0～20 s，滑坡运动启动后，以中后部为主的高速滑动区滑坡向中前部运动，并发生碰撞堆积；20～40 s，滑坡体在前面阶段的冲击与堆积重力作用下，开始以似整体颗粒流的形式滑向坡底；40～60 s，进一步滑向坡底的滑坡体颗粒发生碰撞堆积，并进行自稳过程；在70 s时完成自稳过程，滑坡体整体停止运动。值得注意的是，在20 s之后，滑体似整体颗粒流的形式滑向坡底，其整体性逐渐变好，速度分布也较为均匀。

3.2.3滑坡体内压力分布特征

滑坡体内压力的分布是随着运动的变化而变化的。从图7 A截面上不同时刻压力分布图中可看到：0～20 s，滑坡运动启动时中后部滑坡体所受压力较小，在初始动能和势能作用下开始向中前部运动，并发生碰撞堆积，滑坡体内颗粒所受压力累积增大；20～40 s，中前部滑坡体在中后部滑坡体的冲击与挤压作用下滑向坡底，此时滑坡体内部压力开始释放减小；40～60 s，进一步滑向坡底的滑坡体颗粒发生碰撞堆积，随着颗粒逐步的聚集，滑坡体内颗粒所受压力再一次累积增大；到70 s时，滑坡体内部压力分布不再变化。另外，与连续性刚体不同的是，由于滑坡体颗粒具有离散性与不均匀性，滑坡体内压力的分布是不均匀的，且内部颗粒所受压力也不是在最厚重的地方就一定最大。

3.3 滑坡内部监测颗粒运动的速度与距离

为更进一步分析滑坡体内部颗粒运动特征，对截面A上的11个典型监测点颗粒的运动进行了研究。值得注意的是，虽然在滑坡运动过程中监测颗粒的相对位置会发生改变，但还是能反映出一些滑坡体的运动特点。不同时刻典型监测颗粒的运动速度及运动距离分别见图10和图11。

从图10(a)、(b)、(c)、(d)中的曲线1～11可知，在0～5 s滑坡运动启动后的加速阶段，典型监测颗粒的加速度从后缘、后部、中部、前部到前缘，表层、中层至底层大致由大到小分布。这可从图7 A截面上的压力分布得到解释，初始阶段滑坡体后缘表层所受压力及阻力相对较小，其加速度也较大。因此滑坡体内同一横剖面不同深度的滑坡体颗粒，其运动速度大致表现出离表层越近速度越大的特点；在5 s之后，滑坡体内同一横剖面不同深度的滑坡体颗粒速度关系发生了变化，这是因为颗粒位置及受力情况都发生了改变，这也可从图8 A截面上的颗粒分布图和图7 A截面上的压力分布图得到解释。

图10 滑坡监测点速度Fig.10 Velocity of monitor particles of landslide

图11 滑坡监测点中心运动距离Fig.11 Distance of monitor particles of landslide

从图10(a)中曲线1看到前缘颗粒速度在25～38 s的时候有较大不连续的波动，这说明它正经过不同曲率的滑坡边界运动至坡底。同样地，由图11曲线2～10可知：前部、中部滑坡体颗粒大约分别在25 s、30 s经过不同曲率的边界开始向坡底运动；后部、后缘滑坡体颗粒则在35 s、40 s经过滑床底部拐点运动至滑坡边界上。

在图10(b)中，原来位于前部表层的颗粒4的速度在55 s后处在了原来位于前部中层的颗粒3下面；图10(c)中，原来位于中部表层的颗粒5的速度在55 s后处在了原来位于中部底层的颗粒7下面。这说明滑坡体内部颗粒位置的确发生了很大的变化，滑坡体运动无法以完全的整体式颗粒流运动，这从图8滑坡体截面上颗粒的分布变化也可以说明。

整个滑坡运动过程中监测颗粒1～11的运动距离分别为：150.643 m、146.576 m、125.662 m、171.64 m、217.484 m、203.217 m、186.862 m、229.388 m、227.097 m、240.755 m、277.274 m。其中，后缘监测颗粒11的运动距离最大，即277.274 m。从图11中的曲线1～11可知，在0～5 s滑坡启动后的加速阶段，监测颗粒的运动距离从滑坡后缘至前缘由大到小分布，从滑坡体同一横剖面的表层至底层由大到小分布，即随着滑坡体深度的增加，滑动距离逐渐减小；5～15 s，依旧符合运动距离从后缘到前缘，表层至底层由大到小分布的规律。但滑坡体内部颗粒位置已经发生了很大变化，从图中曲线2、3就可以看到，前部中层监测颗粒3相对于前部底层监测颗粒2运动到了滑坡体内更深的位置。

4 结论

本文总结了黄坪滑坡的基本特征，利用离散元方法对滑坡过程进行了模拟。通过由表及里，从整体到局部，深入细致的分析滑坡启动-高速滑动-碰撞-滑至坡底-堆积自稳运动过程中的主要特点后，得到了如下结论：

(1)滑坡体表面在运动过程中会发生较大变化，且滑坡体会呈现出中间失稳体积较大，两侧失稳体积较小的特征。

(2)在滑坡启动、高速滑动及自稳堆积阶段，滑坡体内存在速度较小甚至没有运动的滑坡体颗粒，当滑坡体处于整体流动阶段后，颗粒几乎都在运动。

(3)与将滑坡体看成刚体不同的是，离散滑坡体运动中：滑坡内部颗粒分布会发生改变，存在大颗粒运动至滑坡体表面，小颗粒运动到滑坡体底面的现象；颗粒的速度和压力在不同位置也会发生改变，且存在不同颗粒间速度相差很大，滑坡体内压力分布不均匀的情况。

(4)滑坡体颗粒的运动速度、运动距离表现出从后缘到前缘，同一横剖面中表层至底层由大到小分布的规律。