滑坡是一種常見的地質災害，其危害嚴重，會對人們的生命、生產、生活造成極大的威脅和影響。對滑坡運動過程進行數值分析可以模擬滑坡的運動和發展過程，計算出滑坡的致災范圍等重要的災害特征參數，為建筑場地的合理選擇和規避滑坡提供建議，進而研究如何采用安全經濟且有效合理的措施防治滑坡災害，并為定量的危險性分析提供依據，因而具有重要的經濟、社會效益及廣闊的應用前景.
國內外學者對滑坡的發展過程進行深入的研究，國外學者 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 工程實例與數值模擬
以美國卡拉韋拉斯大壩滑坡為例進行數值模擬?？ɡf拉斯大壩始建于 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 所示?！緢D1】
由于現有的參數資料并不充分，只能參考已有的相關文獻來實現滑坡整個運動過程的計算模擬。結合滑坡的實際情況并參考相關文獻,得到計算參數如表 1 所示。模擬過程中，取水平距離方向上 0 ~600 m 為計算區域，每 10 m 為一個網格，共劃分成 60 份網格。
流團尺寸為 25 cm ×25 cm,時間間隔 Δt 為0. 1 s,共運行11 000 步，即運行時間共為1 100 s,并采用表1計算參數得到計算結果如圖 2 所示。從圖 2 所示的計算模擬結果與實際破壞后地形線作對比可以發現模擬結果與實際最終堆積曲線輪廓吻合良好，說明了模擬是正確和有效的。圖 3 為滑坡運動過程中主要時刻的堆積形態?！緢D2】
如圖 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 基本吻合?！緢D3】
卡拉維拉斯壩的破壞是由靜態液化造成的，在計算初期滑坡以近似整體的形態向前運動，隨后滑坡前部由于所受阻力較小而快速運動，而滑坡后部由于受到前部的阻滯約束作用運動相對較緩，在最終的堆積形態中滑坡沿運動方向產生變形。
從計算結果中可以看出，模擬結果和實際情況還存在一些差異。這一方面是由于模型的簡化假設造成了與實際復雜條件不符的情況，另一方面在近似平面應變條件下采用典型剖面分析立體的滑坡運動與實際情況也必然存在差異，主要表現為由于側向運動造成的滑坡體積的差異。同時，模擬過程中數值離散方法也仍有待進一步的完善。
3 結 語
目前，國內外對滑坡的研究大多為破壞機理和穩定性分析等方面,而對滑坡破壞后的運動過程研究相對較少，而現有對滑坡運動的計算研究大多采用非變形計算方法或針對巖質邊坡的離散元方法,對于在液化破壞情況下產生連續大變形運動的研究則相對較少。本文應用滑坡模型理論和PIC 數值離散方法，對在液化破壞情況下的典型連續大變形滑坡美國卡拉韋拉斯大壩滑坡破壞后的運動過程進行計算模擬，并將模擬結果與實際情況進行對比分析，結果表明，計算結果與實際情況基本吻合，從而驗證了本文所用方法的有效性和穩定性。
模擬過程中，滑坡最前緣前 100 s 運行了 28 m,平均速度為 0. 28 m/s.最后緣的運動發生在 100 ~300 s.運動過程中運動速度比較低，整體運動比較緩慢，滑坡整體變形主要以慢速蠕動為主，顯示了較低的流動性。通過滑坡運動過程模擬可以得到滑坡的運動發展過程、滑坡的致災范圍和堆積形態等，從而為滑坡防治提供科學的參考和依據?！颈?.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.
