大型地震滑坡高速滑动堵江机制的离散元数值模拟

董金玉, 赵志强, 郑珠光, 杨继红

(华北水利水电大学 资源与环境学院，河南 郑州 450045)



大型地震滑坡高速滑动堵江机制的离散元数值模拟

董金玉, 赵志强, 郑珠光, 杨继红

(华北水利水电大学 资源与环境学院，河南 郑州 450045)

摘要：汶川大地震触发了大量的大型滑坡，这些滑坡体在峡谷河流地带堵塞河道形成了堰塞湖.本文以唐家山滑坡形成的堰塞湖为例，通过离散元数值方法对地震作用下唐家山滑坡的滑动堵江机制进行了模拟,结论如下：唐家山滑坡确实为一高速滑坡，滑坡从启动到停止，其速度变化曲线具有显著的非线性特征，滑坡的最大滑速达27 m/s；直观地再现了唐家山滑坡的滑动以及形成堰塞体的全过程运动特征和滑坡堵江机制，并把其划分为5个阶段,即滑动启动阶段、加速滑动阶段、减速滑动阶段、遇阻堆积阶段和自稳成坝阶段.

关键词：唐家山滑坡;离散元数值模拟;滑坡堵江机制;高速滑坡;全过程运动特征

地震是最为严重的自然灾害之一,不仅直接造成人员伤亡和财产损失，还能引发其他类型的自然灾害，尤其是山体滑坡、崩塌.2008年5月12日在四川省汶川县发生的8.0级地震诱发了大量的滑坡.这些滑坡滑动速度快、能量大、冲击破坏力强，在高速滑动过程中常堵塞河道形成堰塞湖，其中备受人们关注的就是北川县的唐家山滑坡形成的堰塞湖，一度威胁到下游30万人民的生命财产安全.

针对高速滑坡的成因机制、运动过程等，一些学者进行了一些研究，如：殷跃平[1]对汶川地震中高速远程滑坡的特征进行了研究，认为高速远程滑坡具有明显的抛掷效应和气垫效应特征，并导致滑动距离为滑体长度的数倍，堆积成坝后堵江易形成堰塞湖；魏欣等[2]分析了汶川地震灾区高速远程滑坡的空间分布特征，认为随着滑体剪出口相对高差的增大，更多的势能转化为动能，导致滑体滑动速度快，滑动距离远；苏生瑞等[3]以谢家店子滑坡为例，对汶川地震引发高速远程滑坡的运动机理进行了数值模拟研究，得到了地震震级越大，滑体启动的加速度和速度也越大，越易形成高速远程滑坡的结论；陈禄俊等[4]以头寨沟滑坡为例，使用数值模拟方法研究高速滑体凌空飞行的运动规律，认为考虑空气动力学效应分析高速滑坡的滑动时，滑体凌空飞行的时间和距离会增加，更加符合实际滑动情况；邬爱清等[5]采用DDA对唐家山堰塞坝的形成机制进行了模拟分析，计算得到了滑体的滑动速度和块体的主应力分布规律；笔者亦曾结合案例对汶川地震灾区滑坡的成因进行了讨论[6].在这些研究的基础上，笔者通过对唐家山滑坡堰塞湖的现场调查并结合收集到的地质资料，建立离散元数值模型，输入四川卧龙台记录到的实际地震波，对汶川地震作用下唐家山滑坡的动力破坏过程进行模拟，直观地再现唐家山滑坡的高速滑动特征以及形成堰塞坝的全过程运动特征.

1唐家山滑坡概况

唐家山滑坡位于北川上游5 km的湔江右岸，距主中央断裂带的直线距离为2.3 km，位于断裂带的上盘.滑坡地段斜坡陡峻，下部坡陡，坡度约50°，基岩裸露，上部较缓，坡度30°左右，上、下游各分布1条冲沟，上游为大水沟，下游为小水沟，滑坡前、后缘高差650 m左右.唐家山滑坡堰塞体长803 m，宽611 m，厚82.65～124.40 m，方量约2 037×104m3，堰塞体坝顶最低部位高程为752 m，最高部位高程790 m，唐家山滑坡的地貌形态特征如图1所示.

图1　唐家山滑坡地貌形态

唐家山滑坡部位出露的地层岩性为寒武下统清平组的砂岩、泥灰岩、泥岩、灰岩，岩层软硬互层.唐家山堰塞坝的坝体由原山坡上部残坡积的碎石土和寒武系下统清平组上部基岩经下滑、挤压、破碎形成的碎裂岩组成，碎裂岩呈现出似层状结构.唐家山堰塞坝地质横剖面如图2所示.从图2中可以看出，唐家山滑坡滑动前岩层产状为330°∠40°，坡面产状330°∠30°～50°，属于顺层结构岩质边坡.节理裂隙主要发育2组高倾角节理，第1组产状为220°∠75°，第2组产状为335°∠80°，分别对岩体起到侧向和横向切割作用.

图2　唐家山堰塞坝地质横剖面图

2数值模型建立及参数选取

2.1　数值模型

UDEC(Universal Distinct Element Code)是一个处理不连续介质的二维离散元程序，可用于模拟非连续介质(如岩体中的节理裂隙等)承受静载或动载作用下的响应.本文通过UDEC计算程序，首先建立数值模型，然后将实际地震作用力施加到模型中，计算分析唐家山滑坡在地震作用下边坡岩体的运动破坏过程.

根据现场调查和地质资料，建立唐家山边坡离散元计算模型,如图3所示.河流右侧高程500～1 300 m，高差为800 m；河流左侧高程300～700 m，高差为400 m；计算模型断面沿水平方向延伸1 600 m.滑体部分视为刚体，不考虑自身的变形量，滑床为变形体，采用弹塑性本构模型，屈服准则采用Mohr-Coulomb屈服准则[7].结构面主要有两组，一组为岩层层面，另一组为岩体中发育的高倾角节理，层面倾角40°，节理倾角80°.

图3　唐家山边坡的几何模型及网格划分

2.2　岩体物理力学参数选取

根据室内试验和工程类比，将计算所需的各岩土介质及节理的岩体物理力学参数列于表1，阻尼比取岩土体常用阻尼比2%.

表1　计算参数

2.3　地震波输入

地震荷载采用汶川地震卧龙台记录到的地震波，其水平加速度峰值为9.58 m/s2，竖向加速度峰值为9.48 m/s2，动力持续时间为45 s，计算输入的加速度时程曲线如图4所示.

图4　输入的加速度时程曲线

3计算结果及分析

3.1　最大不平衡力

地震作用下边坡的最大不平衡力监测曲线如图5所示.从图中可以看出：在地震波作用下，模型中的最大不平衡力不断震荡变化，随着地震加速度值的增大而增大，并在10 s时达到最大值，之后减小，和加速度值相对应，在37 s左右达到第二次峰值，之后逐渐减小，并在48 s之后迅速减小，趋近于零.这也反映了坡体从滑动启动、高速滑动以及逐渐堆积稳定的演化过程.

图5　计算过程中的不平衡力监测曲线

3.2　高速滑动堵江机制

为了研究唐家山滑坡在地震作用下的详细滑动过程和堵江机制，在滑体的前、中和后部设置了速度监测点.图6为滑体前、中、后部位单元体滑动速度随时间的变化过程曲线.从图中可以看出：滑坡体前缘监测块体在7 s时滑到河床，水平速度为-17 m/s，垂直速度为-9 m/s，其中负号表示x轴和z轴的负向；之后从7 s到12 s块体沿着河床滑移，垂直速度逐渐减小到0 m/s左右，水平速度增加到绝对值最大值，为-22 m/s;再后由于滑体前缘受坡体阻挡，但滑体在后缘的推动作用下向对岸爬坡，垂直速度开始变为正值，水平速度绝对值开始减小；在30 s时，滑动速度减小至较小值；50 s时趋近于0 m/s，滑动基本停止.而在滑体后缘，在15 s时滑动速度的绝对值达到最大值，水平速度为-22 m/s，垂直速度为-16 m/s，之后滑动速度的绝对值逐渐减小，在30 s时，滑动速度的绝对值减小至较小值，50 s时趋近于0 m/s，滑动基本停止.

唐家山滑坡确实为一高速滑坡，滑坡从启动到停止，其速度变化曲线具有显著的非线性特征.滑坡滑动总的持续时间约30 s，其中，在0~15 s内滑坡处于滑动加速阶段，在15~30 s处于滑动减速阶段.滑坡的最大总滑动速度达到27 m/s，最大滑动距离约为750 m.

图6　滑坡体前、中、后部位速度监测曲线

本文对唐家山滑坡的整个滑动堵江过程进行了全过程的离散元模拟，整个模拟过程用了70 s，模拟结果如图7所示，直观地再现了唐家山滑坡滑动堵江的全过程及其运动特征.

图7　地震滑坡堵江全过程模拟结果

根据滑动的破坏状态，结合前面的滑动速度监测曲线，将该滑坡的变形破坏过程分为5个阶段.

1)滑动启动阶段：0.0～0.1 s.由于地震力及滑体重力的作用，滑体沿着岩体层面(层间软弱带)产生应力集中，使层面产生塑性变形并逐渐贯通.

2)加速滑动阶段：0.1～15 s.在地震力和滑体势能转化动能的作用下，滑体滑动速度急剧增加，在此阶段末滑动速度达到最大值27 m/s.

3)减速滑动阶段：15～30 s.由于河谷宽仅约100 m，滑体在滑动过程中，受到对岸的阻挡作用力，同时由于滑床与滑体之间的摩擦阻力，滑动速度绝对值逐渐减小，在此阶段末，滑动速度绝对值减小到0 m/s左右.

4)遇阻堆积阶段：30～50 s.滑体速度绝对值较小，滑动基本停止，在地震力的作用下，局部地带发生蠕滑或者震荡变形.

5)自稳成坝阶段：50 s至计算完成.在50 s以后，地震力作用已经结束，滑体在自重作用下逐渐变得密实、稳定，堵塞河流，形成堰塞坝.

4结语

基于现场调查和地质资料，本文建立了唐家山边坡离散元数值计算模型，通过在滑床和滑体不同部位设置监测点，对唐家山滑坡在实际地震波作用下滑动破坏的堵江过程进行了模拟分析，得到了以下结论：

1)唐家山滑坡确实为一高速滑坡，滑坡从启动到停止，其速度变化曲线具有显著的非线性特征，滑坡的最大滑速达27 m/s.

2)对唐家山滑坡的整个滑动堆积过程进行了全过程的离散元模拟重现，直观地再现了唐家山滑坡的全过程运动特征，并把其划分为5个阶段：滑动启动阶段、加速滑动阶段、减速滑动阶段、遇阻堆积阶段、自稳成坝阶段.

参考文献

[1]殷跃平.汶川八级地震滑坡高速远程特征分析[J].工程地质学报,2009,17(2):153-166.

[2]魏欣,胡瑞林,李丽慧,等.强震条件下高速滑坡的空间分布特征研究[J].工程地质学报,2010,18(4):490-496.

[3]苏生瑞,张永双,李松,等.汶川地震引发高速远程滑坡运动机理数值模拟研究—以谢家店子滑坡为例[J].地球科学与环境学报,2010,32(3):277-287.

[4]陈禄俊,邢爱国,陈龙珠,等.高速远程滑坡飞行数值分析[J].水文地质工程地质,2008,35(5):1-6.

[5]邬爱清,林绍忠,马贵生,等.唐家山堰塞坝形成机制DDA模拟研究[J].水文地质工程地质,2008,39(22):91-95.

[6]董金玉,杨国香,杨继红,等.汶川地震灾区滑坡的成因及典型实例分析[J].华北水利水电学院学报,2011,32(5):10-13.

[7]Itasca Consulting Group,Inc..Universal distinct element code:theory and background[R].Minneapolis:Itasca Consulting Group Inc.,2005.

(责任编辑：乔翠平)

Discrete Element Numerical Simulation of the Mechanism of the Large-scale

Earthquake High-speed Landslide′s Blocking the River

DONG Jinyu, ZHAO Zhiqiang, ZHENG Zhuguang, YANG Jihong

(School of Water Resources and Environment, North China University of Water Resources and Electric Power, Zhengzhou 450045, China)

Abstract:A large number of landslides were triggered by Wenchuan earthquake, and many barrier lakes were formed for the slumping bodies blocking the rivers. Taking the Tangjiashan landslide as an example, the mechanism of its blocking the river is simulated based on the discrete element numerical simulation method. The results show that the Tangjiashan landslide is indeed a high-speed landslide, the variation curve of its velocity from starting to stopping is obviously nonlinear, and the maximum sliding velocity reaches 27 m/s. The movement and formation mechanism of the barrier lake is reappeared intuitively, and the whole process can be divided into five phases: the starting sliding stage, the speed-up sliding stage, the decelerated sliding stage, the accumulation stage and the self-stabilization stage.

Keywords:Tangjiashan landslide; discrete element numerical simulation; mechanism of landslide′s blocking the river; high-speed landslide; the whole process of motion features

文献标识码：A

文章编号：1002-5634(2015)06-0047-04

中图分类号：P642

DOI:10.3969/j.issn.1002-5634.2015.06.012

作者简介：董金玉(1977—)，男，河南济源人，副教授，博士，主要从事地质工程方面的研究.

基金项目:国家自然科学基金青年基金资助项目(41102203); 水利部公益性行业科研专项经费项目(201301034)；河南省科技创新人才计划；河南省基础与前沿技术研究项目(122300410146)；华北水利水电大学青年科技创新人才支持计划项目.

收稿日期:2015-07-16