滑坡是一种常见的地质灾害，其危害严重，会对人们的生命、生产、生活造成极大的威胁和影响。对滑坡运动过程进行数值分析可以模拟滑坡的运动和发展过程，计算出滑坡的致灾范围等重要的灾害特征参数，为建筑场地的合理选择和规避滑坡提供建议，进而研究如何采用安全经济且有效合理的措施防治滑坡灾害，并为定量的危险性分析提供依据，因而具有重要的经济、社会效益及广阔的应用前景.
国内外学者对滑坡的发展过程进行深入的研究，国外学者 Sassa对滑坡运动进行了较早的研究，把滑坡近似看作流体，推导出滑坡运动的微分方程，实现了滑距的预测。Hungr开发了模拟滑坡运动的分析软件 DAN,并对一些典型的滑坡进行了全过程的模拟。潘家铮在《建筑物的抗滑稳定和滑坡分析》一书中，采用基于条分的方法，提出了滑坡速度的计算公式，但该方法并未考虑滑坡运动过程中的变形。廖小平等探讨了高速远程滑坡的形成条件和运动规律，认为滑坡运动是一种连续可变的块体运动，提出了冲击性碰撞作用机理和连续可变的块体运动理论。刘忠玉等基于其对高速滑坡发生机理的认识及对运动特征和堆积特征的分析，建立了预测高速远程滑坡的块体运动模型，应用该模型可以预测出滑坡的最大滑速和最大滑距。刘涌江等应用力学原理从理论上分析了大型高速滑坡岩体与阻挡山体的碰撞过程，得出了碰撞后滑坡岩体的平均运动速度和运动方向的关系式，为进一步研究滑坡的运动奠定了基础。黄润秋根据在中国西部地区多年的研究经验，对几类典型的大型岩质斜坡变形及失稳机理进行了较为深入的分析，阐述了这些大型滑坡形成的地质条件、斜坡变形过程及发生机理的概念模型。近年来，针对岩质边坡和高速远程碎屑滑坡运动的离散单元法的运动模拟以及滑坡的治理研究也取得了很大的发展.以上学者在对滑坡运动的速度计算、滑距预测等方面所做的工作为滑坡发展运动过程的研究做出了巨大的贡献。
本文尝试采用滑坡运动模型方程和 PIC（ Particle In Cell） 数值离散方法,通过编写计算程序，对卡拉韦拉斯大坝滑坡破坏后的运动过程进行计算模拟，进而得到其滑动过程、致灾范围和堆积形态等结果，再现滑坡的运动过程，并对过程结果进行分析，旨在为今后具有类似特征滑坡的防治提供参考和依据。
1 模型理论和 PIC 数值方法
描述滑坡运动的模型方程可以表示为：【1】
方程（ 1） 中，h 为滑坡深度，u 代表滑坡速度，z0为原始地形高程，g 表示重力加速度，C 为黏聚力，γm为滑坡体容重，φ 表示摩擦角，sgn（ u） 为表示速度方向的标记函数。
PIC 数值方法为计算流动问题的欧拉 - 拉格朗日混合方法，基本要点是把流动所通过的区域用欧拉法分成有限个网格，每个网格中流体再用一组特定的离散化拉格朗日质点表示。每个质点具有一定的质量，每个网格单元内的质点数目和质点分布都以流体流动的初始状态为依据，而且这些质点具有一定的速度和能量。计算开始后，质点在欧拉网格之间迁移，表示流体在运动。滑坡运动的空间由滑坡体流团连续地、无空隙地充满着，流团的尺寸大小基本符合微观足够大、宏观足够小的条件。每个流团受到体积力及附近流团或边界对其的力作用，根据牛顿力学基本运动定律，流团的加速度、速度与空间坐标可以确定。在计算时间域上对这些参量不断更新，并在计算空间域上运用统计平均原理将这些运动参量表达在网格节点上，即可得到滑坡堆积过程的总体运动规律.
通过对上述模型方程进行有效的简化，并应用 PIC 数值离散方法对其离散，采用的本构模型关系为摩尔库伦准则，得到如下计算方程:【2】
其中，up,t +1与 up,t分别为 t +1 及 t 时刻流团 p 的速度，Δt 为时间步长，Sp,t为 t 时刻流团 p 处于网格单元上任一位置对应的滑坡表面重力坡降，而用于计算网格节点表面高程的流深则由节点单位面积控制流团数量计算得出.令方程（ 2） 右端项为流团合加速度 ap,t,则在计算时间步长内的速度关系可表达为：【3】
流团在计算时间步长内的位置关系为：【4】
方程（ 4） 中，Xp,t +1与 Xp,t分别为流团 p 在两个不同时刻的空间坐标。
2 工程实例与数值模拟
以美国卡拉韦拉斯大坝滑坡为例进行数值模拟。卡拉韦拉斯大坝始建于 1914 年，是当时世界上最大的土坝，其坝高最大达到 73 m.1918 年 3 月 24 日早上，坝体上游面的东侧坝体有 60 000 m3的材料破坏并滑进水库。破坏时，坝的最大高度约有 61 m,水库最深深度约有 23 m.原有 385 m 宽的大坝中有 215 m 的坝体部分发生了破坏。
坝底附近材料饱和重度大约是18. 9 kN/m3,坝顶附近材料的饱和重度大约是 17. 3 kN/m3.组成大坝壳体的填砂中细粒含量在10% ~50%.构成粘土心墙材料粒径在0. 002 ~0. 2 m,重度为15. 7 kN/m3.
相关研究表明，卡拉韦拉斯大坝发生破坏主要是由液化引起的，液化使土的抗剪强度大大地降低，并进一步导致孔隙水压力增加，最终发生破坏.
据 Simpson 等的现场调查和稳定性分析，得到了卡拉韦拉斯大坝滑坡破坏前的剖面形态和破坏后的滑动面及最终的堆积形态，如图 1 所示。【图1】
由于现有的参数资料并不充分，只能参考已有的相关文献来实现滑坡整个运动过程的计算模拟。结合滑坡的实际情况并参考相关文献,得到计算参数如表 1 所示。模拟过程中，取水平距离方向上 0 ~600 m 为计算区域，每 10 m 为一个网格，共划分成 60 份网格。
流团尺寸为 25 cm ×25 cm,时间间隔 Δt 为0. 1 s,共运行11 000 步，即运行时间共为1 100 s,并采用表1计算参数得到计算结果如图 2 所示。从图 2 所示的计算模拟结果与实际破坏后地形线作对比可以发现模拟结果与实际最终堆积曲线轮廓吻合良好，说明了模拟是正确和有效的。图 3 为滑坡运动过程中主要时刻的堆积形态。【图2】
如图 3 所示，滑坡运动到第100 s 时，滑坡最前缘水平距离为288 m,100 s 内运动了28 m,平均速度为0. 28 m/s.而滑坡后部水平位移基本不变，即滑动现象不明显，而开始发生滑动时在100 ~300 s.滑坡最前缘运动的水平距离及平均速度如表 2 和表 3.最终运行到第 1 100 s 时，最前缘水平距离 470 m与实际水平距离 471 m 基本吻合，最后缘水平距离 29 m 也与实际水平距离 26. 5 m 基本吻合。【图3】
卡拉维拉斯坝的破坏是由静态液化造成的，在计算初期滑坡以近似整体的形态向前运动，随后滑坡前部由于所受阻力较小而快速运动，而滑坡后部由于受到前部的阻滞约束作用运动相对较缓，在最终的堆积形态中滑坡沿运动方向产生变形。
从计算结果中可以看出，模拟结果和实际情况还存在一些差异。这一方面是由于模型的简化假设造成了与实际复杂条件不符的情况，另一方面在近似平面应变条件下采用典型剖面分析立体的滑坡运动与实际情况也必然存在差异，主要表现为由于侧向运动造成的滑坡体积的差异。同时，模拟过程中数值离散方法也仍有待进一步的完善。
3 结 语
目前，国内外对滑坡的研究大多为破坏机理和稳定性分析等方面,而对滑坡破坏后的运动过程研究相对较少，而现有对滑坡运动的计算研究大多采用非变形计算方法或针对岩质边坡的离散元方法,对于在液化破坏情况下产生连续大变形运动的研究则相对较少。本文应用滑坡模型理论和PIC 数值离散方法，对在液化破坏情况下的典型连续大变形滑坡美国卡拉韦拉斯大坝滑坡破坏后的运动过程进行计算模拟，并将模拟结果与实际情况进行对比分析，结果表明，计算结果与实际情况基本吻合，从而验证了本文所用方法的有效性和稳定性。
模拟过程中，滑坡最前缘前 100 s 运行了 28 m,平均速度为 0. 28 m/s.最后缘的运动发生在 100 ~300 s.运动过程中运动速度比较低，整体运动比较缓慢，滑坡整体变形主要以慢速蠕动为主，显示了较低的流动性。通过滑坡运动过程模拟可以得到滑坡的运动发展过程、滑坡的致灾范围和堆积形态等，从而为滑坡防治提供科学的参考和依据。【表2.3】
然而，计算结果与实际情况不可避免地存在一定的差异，这些差异一方面是由于模型简化和数值离散方法选用造成的，另一方面平面应变计算假定和关键计算参数的选取也将对滑坡运动堆积的计算模拟结果产生影响。因此，需要进一步研究滑坡运动过程中的复杂机理，发展更加能够表现复杂滑坡发展行为的滑坡模型和准确计算运动堆积的数值方法，同时细化参数研究。
参考文献：
[1] SASSA K. Analysis on slope stability: I. Mainly on the basis of the indoor experiments using the standard sand produced inToyoura,Japan[J]. Journal of the Japan Society of Erosion Control Engineering,1972,25（ 2） : 5-17.
[2] SASSA K. Analysis on slope stability: II. Mainly on the basis of the indoor experiments using the standard sand produced inToyoura,Japan[J]. Journal of the Japan Society of Erosion Control Engineering,1974,26（ 3） : 8-19.
[3] HUNGR O. A model for the runout analysis of rapid flow slides,debris flows,and avalanches[J]. Canadian GeotechenicalJournal,1995,32: 610-625.
[4] HUNGR O,MCDOUGALL S. Two numerical models for landslide dynamic analysis[J]. Computers & Geosciences,2009,35（ 5） : 978-992.
