基于离散单元法的
三维滑坡过程数值模拟分析
赵
川1，2，
付成华1，2，邹海明3
，何欢2，
钟学梅2
（1．西华大学流体及动力机械教育部重点实验室，四川
成都610039；
2．西华大学能源与动力工程学院，四川
成都610039；3．巴中市水务局，四川
巴中636000）
摘
要：为进一步分析滑坡体的运动过程，提出了一种采用离散单元法模拟滑坡运动过程的方法，通过与有限元模
型对比验证其合理性，并用于三维滑坡数值模拟。计算结果表明：滑坡体在t =3s 时平均速度最大，为6.4m /s ，此
时滑坡体动能达9．79ˑ108
J ，滑坡体在运动过程中产生大量飞石，最大速度可达21．3m /s ，
在t =9s 后运动基本停止；滑坡运动停止后，滑坡体前缘和后部土体结构发生破坏，中部基本完整，确定了滑坡影响范围距离坡脚20m 左右。计算结果对滑坡灾害评估和边坡防护工程设计有一定参考价值。关键词：滑坡体；离散单元法；数值模拟；动能；飞石中图分类号：P642．22文献标识码：A 文章编号：1001-9235（2015）02-0012-04基金项目：西华大学流体及动力机械省部共建教育部重点实验室学术成果培育项目（SBZDPY －11－9）；西华大学研究生创新基金项目（ycjj2014067）
收稿日期：2015-02-03作者简介：赵川，男，四川内江人，主要从事水利水电工程、岩土工程以及边坡工程加固设计等工作。
前言
滑坡是我国一种常见的地质灾害，特别是西南地区，每年都造成了数以亿计的财产损失，严重地制约着区域经济的可持续发展。受到自然环境和人为因素的影响，近年来我国大型滑坡灾害发生的频率和规模呈现明显的上升趋势
［1－2］
。
目前，
边坡稳定性分析方法主要以定量分析法为主，定量分析法又可分为极限平衡法和数值分析法，极限平衡法主要适用于土质边坡，而针对结构复杂的岩质边坡，需采用数值分析法，其中又以有限元法和离散元法为主
［3］
。近年来，
有限元法在边坡工程中得到了广泛应用，有限元法基于等效连续介质基本理论，对于一般均质边坡和含有少数结构面的边坡分析已经取得了较好的分析效果
［4－5］
。但是，对于含有
大量不同构造、产状和特性的不连续结构面的节理化岩质边坡，该方法不适用。因此，基于非连续介质理论的离散单元法应运而生。本文基于离散单元法，建立三维离散元滑坡模型，模拟整个滑坡过程，分析得到滑坡体的运动状态和力链结构特征，对滑坡防灾减灾工作有积极意义［6－7］
。
1
离散单元法简介
离散单元法（DEM ）是由美国学者Cundall P. A.［8］
教授
在1971年提出的一种颗粒离散体分析方法，此法主要应用于岩土领域。离散单元法基本思想是把不连续体分离为刚性元素的集合，各刚性元素满足运动方程，采用时步迭代的方法求解各个刚性元素的运动方程，每一步结束后更新各元素运动状态，重复以上过程，继而求得整个非连续体的运动
形态。离散单元法基于非连续介质基本理论，与基于连续介质理论的有限单元法（FEM ）相比，更适合求解岩体内大量的节理与不连续面问题，为求解复杂岩体大变形、大位移和非线性问题提供了一种有效的技术手段［9－10］
。
2
离散元模型合理性验证
目前，边坡工程中常采用有限元强度折减法对边坡进行稳定性进行分析，当计算不收敛时的折减系数即为边坡的稳定安全系数，并可以得到边坡滑裂面。强度折减法同样可以应用于边坡离散元分析中。为验证离散单元法的准确性，分别建立了二维边坡有限元模型和离散元模型，模型长150m ，高100m ，坡高50m ，坡角60ʎ。同时采用强度折减法对有限元模型和离散元模型进行折减计算，直到边坡发
生破坏，得到边坡滑裂面和滑坡体，见图1、2。分析可以发现，
离散单元
图1
有限元法计算的滑坡体
人民珠江2015年第2期·PEAＲL ＲIVEＲdoi ：10.3969/j.issn.1001-9235.2015.02.004
图2离散元法计算的滑坡体
法计算得到的滑坡体与有限单元法计算的结果很接近，由此认为采用离散单元法对边坡进行稳定性分
析计算是可行的。3
三维离散元模型滑坡分析
由于离散元计算效率受颗粒数的影响很大，当颗粒数量太多时计算速度很慢，甚至无法计算，本文旨在研究滑坡体的运动过程。因此，只对滑坡体部分建立离散元颗粒模型，离散元模型由4456个颗粒组成，颗粒半径1m 。而对滑坡体下部岩土体采用刚性材料模拟，与离散元滑坡体之间的接触作用采用输入相关摩擦系数来模拟，本模型中摩擦系数取0.45。基于图2的二维离散元模型，建立了三维边坡模型，见图3a 。离散元模型采用粘结接触模型，设置粘结参数见表1。
表1
离散元计算初始力学参数
颗粒模型法向刚度/（N ·m －3）剪切刚度/（N ·m －3）临界法向强度/Pa 临界剪切强度/Pa 粘结半径/
m 粘结接触模型
2ˑ109
1ˑ109
5ˑ107
1.6ˑ107
1.3
3.1滑坡体运动过程速度分析
为从宏观上分析三维滑坡体的运动过程，计算得到了不
同时刻，滑坡体的运动状态，见图3。t =0s 时，边坡处于临界状态，有下滑的趋势，之后滑坡体开始滑动；t =1.5s 时，坡体上部土体首先破坏，上部颗粒速度较大，
部分土体颗粒达
a ）t =0s
图3
不同时刻三维离散元模型滑坡速度分布
b ）t =1.5
s
c ）t =3.5
s
d ）t =6
s
e ）t =9s
f ）t =15s
续图3不同时刻三维离散元模型滑坡速度分布
到7.52m /s ，而下部保持完整；当t =3.5s 时，滑坡体前缘部分抵达边坡底部平面，造成前缘土体结构破坏，前缘土体颗粒速度可达10m /s ；直到t =9s 后，滑坡体运动速度逐渐减小，只有滑坡体表面上还有部分散落的土体颗粒在运动，整个滑坡过程持续了大约15s 。
以上从宏观上分析了滑坡体的运动过程，为进一步从微观上分析土体颗粒运动过程，图4给出了滑坡体所有土体颗粒的平均速度和颗粒的最大速度随时间的变化曲线。分析可知：滑坡体所有土体颗粒平均速度在t =3s 的时候达到最大值，为6.4m /s ，与上述滑坡体宏观速度分析结果一致，说明滑坡体前缘刚抵达底部平面时运动速度最大，之后与底面逐渐接触碰撞，滑坡体整体速度随之下降，并在t =10s 过后，除了少数颗粒仍在运动，滑坡体整体运动基本停止。从图4中分析可知，土体颗粒最大速度在t =5.6s 的时候取得最大值，速度达到21.3m /s ，而整个过程波动较大，说明在滑坡体运动过程中，上部会产生部分破坏的土体颗粒，在坡体上部发生滚动，由于滚动摩擦系数远远小于滑动摩擦系数，所以脱离母体的颗粒可以以很大的速度迅速滚下，这些颗粒可以理解为真实滑坡过程中产生的飞石，十分危险
。
图4滑坡体速度
时间曲线
图5滑坡体动能 时间曲线
3.2
滑坡体动能分析
滑坡体的动能反映了坡体运动过程中所携带的能量，动
能分析可以评估其造成破坏的程度。图5给出了滑坡体动能随时间变化曲线，可以看出，在t =3s 的时候，滑坡体动能出现峰值，这与上述分析t =3s 土体颗粒平均速度取得最大
值反映的结果相同，此时滑坡体动能可达9.79ˑ108
J ，滑坡
体动能巨大，若事前不做好防灾减灾措施，滑坡体将在10s 之内携带巨大能量冲击坡脚附近建筑，将对人民生命财产安全造成巨大损失。3.3
滑坡体微观结构分析
为进一步分析滑坡体运动过程中，岩土体内部结构特征变化，单独将滑坡体运动前后的力链图导出，见图6。滑坡体中蓝色线条代表完好的土体力链结构，部分红色线条代表破坏的力链结构，可以理解为岩土体中天然存在的局部细小节理裂隙。滑坡发生后，滑坡体前后部分土体均发生破坏，滑坡体尾部土体主要受拉破坏，前缘部分则主要是挤压破坏，而滑坡体中部土体结构几乎保持完整，前缘部分由于力链结构破坏，部分颗粒将脱离母体，向前滑出，以此可确定滑坡灾害影响范围为距离坡脚20m 左右，部分颗粒可能滚动得更远。模拟结果可为滑坡灾害的评估提供参考
。
a ）滑坡前
图6滑坡前后滑坡体力链结构
b）滑坡后
续图6滑坡前后滑坡体力链结构
4结语
a）基于离散单元法，建立了二维边坡离散元模型，采用强度折减法计算得到了滑坡体，通过与有限元法计算结果进行比较，验证了离散单元法模拟滑坡的可靠性。
b）建立三维离散元滑坡模型，计算得到了不同时刻滑坡体的速度分布，在t=3s时平均速度达到最大值6.4m/s，此时坡体动能最大，达到9.79ˑ108J，滑坡过程不时有颗粒飞出，速度较大。
c）滑坡停止后，滑坡体前后部分土体均发生破坏，中部保持完整。最后得到滑坡体堆积范围大约20m。
d）以上分析表明本文采用的离散单元法模型滑坡过程是可行的，计算结果可为防灾减灾提供有力支撑。参考文献：
［1］郑颖人，陈祖煜，王恭先，等.边坡与滑坡工程治理［M］.北京：人民交通出版社，2007.
［2］皱丽春，王国进，汤献良，等.复杂高边坡整治理论与工程实践［M］.北京：中国水利水电出版社，2006.
［3］丁参军，张林洪，于国荣，等.边坡稳定性分析方法研究现状与趋势［J］.水电能源科学，2011，29（8）：112－114，212.［4］陈浩.极限分析上限法在边坡稳定性问题中的应用［J］.人民珠江，2014，35（3）：67－70.
［5］郑颖人，赵尚毅，张鲁渝.用有限元强度折减法进行边坡稳定分析［J］.中国工程科学，2002，4（10）：57－62．
［6］雷远见，王水林.基于离散元的强度折减法分析岩质边坡稳定性［J］.岩土力学，2006，27（10）：1693－1698.
［7］崔溦，刘学昆，戚蓝.大型复杂堆积边坡稳定性的离散元分析［J］.工程地质学报，2012，20（2）：222－227.
［8］Cundall P A，Strack ODL.A discrete numerical model for granular asmbles［J］.Geotechnique，1979，29（1）：47－65.
［9］周健，王家全，曾远，等.土坡稳定分析的颗粒流模拟［J］.岩土力学，2009，30（1）：86－90.
［10］朱颖彦，崔鹏，陈晓晴.泥石流堆积体边坡失稳机理的试验与稳定性分析［J］.岩石力学与工程学报，2005，24（21）：3927－3931.
（责任编辑：李泽华）
Analysis of Nuerical Simulation of Three－dimensional
Landslide Process with Discrete Element Method
ZHAO Chuan1，2，FU Chenghua1，2，ZOU Haiming3，HE Huan2，ZHONG Xuemei2
(1.Key Laboratory of Fluid and Power Machinery under Ministry
of Education，Xihua University，Chengdu610039，China;
2.School of Energy and Power Engineering of Xihua University，Chengdu610039，China;
3.Ba Zhong Water Authority，Bazhong636000，China)
Abstract:For further analysis of landslide motion process，discrete element method was ud to simulate the movement of landslide，verified its rationality by contrast with the finite element model，and DEM was ud for the numerical simulation of3D landslide.The results showed that:the average maximum speed of the landslide reached6.4m/s at t=3s，and the landslide kinetic energy was about 9.79ˑ108J，landslide body produced a lot flyrock in the process of movement，the maximum speed was21.3m/s，the landslide stopped at t=9s;when the landslide stopped，the front and rear soil structure damaged，the middle part kept integrity，determined the influence scope of the landslide was about20m of the toe.The results could provide some reference for the landslide hazard asssment and slope protection engineering design.
Keywords:landslide;discrete element method;numerical simulation;kinetic energy;flyrock
