基于離散單元法的滑坡堆積及其涌浪計算

徐 寅 ，陳勝宏

（1.武漢大學 水資源與水電工程科學國家重點實驗室，武漢 430072；2.南昌工程學院 水利與生態工程學院，南昌330029）

1 引 言

邊坡分布廣泛，對邊坡穩定性的研究一直是工程地質和巖土工程界研究的一個難點和熱點問題。我國是一個滑坡災害極為頻繁的國家，尤其是隨著西部大開發戰略的進一步實施，一批大型或特大型的水電工程在我國西南高山峽谷地區興建，滑坡更是表現出規模大、危害大和機制復雜等特點，在全世界范圍內具有典型性和代表性[1]。由于巖體中往往含有斷層、節理、裂隙、軟弱夾層等大量不同構造、產狀和特性的不連續結構面，給巖質邊坡的穩定性分析帶來了巨大的困難[2]。

目前，研究邊坡穩定的方法主要有剛體極限平衡法及有限單元法。剛體極限平衡法以原理簡單、有大量的工程經驗且有與之相配套的規程規范，在巖土工程中應用較廣。有限單元法能考慮真實的應力-應變關系，為準確地模擬邊坡的變形破壞提供了強有力的手段。盡管如此，上述兩種方法對分析和預測邊坡的整個破壞過程還無能為力。然而，給人類帶來巨大災害的往往是大型高速滑坡的失穩。它的失穩會導致巖土體高速入水，不僅會沖毀水工建筑物、航行船只、堵塞河道，而且還會激起巨大的涌浪，威脅庫區壩體及沿岸居民生命財產安全[3]。要想分析和預測這類高速滑坡的失事影響范圍及處理措施，需要研究滑坡失穩的整個過程。

20世紀70年代初由Cundall[4]提出的離散單元法是一種將弱面所切割的巖體視為復雜的塊體的集合體，允許各個塊體可以平移或轉動，甚至相互分離。離散單元法最初用于模擬巖石邊坡的漸進破壞過程，之后經歷從二維到三維，從剛性塊體到變形塊體的發展過程[5]。目前，離散單元法已廣泛應用于邊坡[6]、地基[7]、地下洞室[8]等的破壞分析。離散單元法的最大特點是能實現邊坡從穩定狀態到極限平衡狀態再到邊坡分解、脫離、滑動的過程。因此，采用離散單元法對邊坡失穩的整個過程進行研究是可行的且非常有意義的。本文采用離散單元法研究庫岸邊坡失穩后的滑坡速度并對滑坡體滑入水庫中所產生的涌浪進行了評估，并通過離散單元法研究了滑坡的堆積形狀，以便對滑坡的災害進行評估。

2 離散元軟件開發

前后處理相對于數值計算來說，雖然是服務性質的，但它們的重要性不言而喻，是數值程序是否完善的一個重要標志，由于離散單元法允許塊體與塊體可以脫離開來，塊體的位置變化可能較大，因此，前后處理器對離散單元法來說顯得尤為重要[9]。本文的計算流程如圖1所示，圖中的“讀入數據”和“生成DEM模型”模塊相當于前處理部分，“輸出結果”模塊相當于后處理部分。為了直觀地描繪離散元模型及其計算后所得到的塊體運動，本文根據圖1所給出的流程圖并采用VC++與OpenGL相結合的方法，開發了具有前后處理功能的離散單元法軟件CORE_DEM。

3 離散單元法的基本原理

離散單元法的實質就是將求解區域離散成有限個完全獨立的塊體，使求解區域成為一個塊體的系統，再對塊體之間的接觸關系進行判斷，以確定塊體的受力，最后采用差分方法求解基于牛頓第二定理的塊體運動方程。

圖1 計算流程圖Fig.1 Flow chart of calculation

3.1 接觸判斷

模擬不連續系統破壞過程的復雜力學行為的關鍵之一是快速有效檢測和識別塊體的接觸關系[10]。一般地，接觸判斷分兩個階段，即粗判與細判。粗判的目的是排除不可能接觸的塊體，避免處理沒有相互接觸的塊體的檢索，以節省計算時間；細判的目的就是對粗判得到的相關塊體進行進一步的接觸判斷，以提供明確的接觸位置和接觸法線方向，以便塊體之間接觸力的確定。

3.2 塊體運動

離散單元法采用顯式差分方法求解塊體運動系統。在每個時間步上，用運動方程和本構方程來描述塊體的運動。當把塊體視為剛體時，本構方程用接觸力與侵入深度間的關系來表示。首先通過對運動方程的積分來更新塊體的位置、速度及侵入深度增量；然后通過接觸力與侵入深度的關系來確定新的接觸力，在新接觸力的基礎上開始下一步的循環。計算循環見圖2。

圖2 計算循環[5]Fig.2 Calculation cycle[5]

整體坐標系下的平動方程為

式中：m為塊體質量； Fi為施加在塊體上的合力分量(除重力外)，gi為重力加速度分量；x˙i為塊體形心的平動速度；˙x˙i為塊體形心的平動加速度；α為黏性阻尼常數。

用中心差分法可得

式中： D1=1 -(α Δt /2)；D2= 1/[1+ (α Δ t/2)]；t為計算時刻；Δt為時間步長。另外，塊體的轉動方程也可采用類似的方法求解，只要將平動速度換成轉動速度，力換成力矩，質量換成轉動慣量即可。

4 滑速的計算

滑坡體滑速的計算一直是個比較復雜的問題。在計算過程中，不僅要考慮滑坡體所在邊坡的地形地貌、巖性特征和材料強度參數等條件，還要考慮滑坡體在運動過程中滑塊受力條件和滑動面上的強度參數的變化等。因此，在實際工程應用中，常采用簡化方法或經驗公式法求滑坡體滑速，如能量守恒法、潘家錚法[11]、Scheidegger法[12]。

4.1 能量守恒法

如圖3所示，根據動能定理有

式中：g為重力加速度；H為塊體的質心下落的豎直高度；v為塊體速度；Fr為塊體阻力；S為塊體滑行長度。

圖3 塊體滑動示意圖Fig.3 Diagram of block sliding

又有

式中：c為黏聚力；f為內摩擦系數；A為塊體與坡面的接觸面積；N為塊體受到坡面的反力。

求得

計算滑坡的速度采用能量守恒法是合適的，但得到的是某一平均高度滑坡巖土體的勢能轉化為動能的平均速度，而對滑坡體滑入水庫時所激起的最大涌浪高度的估計往往要依賴于獲得滑坡的最大速度，即需要求出滑坡局部巖土體運動到入水口前緣的速度，此時仍采用能量守恒法就不太合適。

4.2 潘家錚法

潘家錚法[11]把滑坡體作為質點系來研究其運動，考慮了加速度的變化。取滑坡開始下滑的瞬間為t0(t0= 0)，則t=ti時的滑坡水平速度為[11]

式中：vx為塊體水平速度；ax為塊體水平加速度；ΔL 為塊體水平滑距，i為滑坡開始運動后的計算時刻標識（如共有n個時段，則i的取值可為0，1，…，n）。滑坡豎直速度的求法與水平速度類似，此時有ay=axtanαs，αs為坡面傾角。

潘家錚法在其推導公式的過程中假定各條塊間完全無摩擦及其他阻力，這一假定可能導致偏于安全的結果。另外，要將該方法推廣至三維情況則較繁瑣。

4.3 Scheidegger法

由于內摩擦系數值在實際應用中不易選取，Scheidegger[12]在搜集了大量滑坡數據的基礎上，進行了統計分析，發現滑坡體的等效內摩擦系數與滑坡體體積間有一定相關關系，給出的經驗公式為[12]

式中：f為等效內摩擦系數；V為滑坡體體積。

Scheidegger法避開了內摩擦系數值難以確定的問題，但其所求得的等效內摩擦系數只反映了滑坡規模的影響，而沒有考慮到滑面材料參數的作用，如黏聚力和內摩擦角等，因此，與實際情況有出入。

4.4 離散元法

盡管滑坡體的下滑速度可通過能量守恒法或經驗公式得到。但同時可以看出：滑坡體失穩后滑速發展過程極為復雜，各種計算方法都含有很多不確定因素，所得成果僅供定性研究之用。然而，離散單元法可以根據所需求解的基本公式（2）直接得到塊體的速度，然后對時步進行積分，就可得到塊體的位置，這樣用離散單元法模擬滑坡運動時可直接獲得任意時刻、任意塊體的位置及滑速，為后續的涌浪計算提供了基礎數據，與上述文獻中給出的經驗公式的方法相比，用離散單元法計算滑坡體的滑速具有更方便、更高效的特點。

4.5 算例考核

算例在xz剖面的幾何尺寸如圖4所示，y方向寬為20 m，塊體尺寸為m× 1 m× 1 m。計算參數為時間步長 Δt = 10-5s ，總時間步長為 106步，邊界面的法向剛度 Kn= 1010N/m3，切向剛度Ks=4.0× 109N/m3，塊體重度 γ= 20.0 kN/m3，接觸面上的黏聚力為2 kPa，內摩擦角為10°。

圖4 滑坡算例 (單位: m)Fig.4 A landslide example (unit: m)

圖5給出了能量法、潘家錚法和離散單元法計算塊體在坡面上滑行前3 s的速度，離散單元法計算中每1 000時步（即0.1 s）輸出一次結果，輸出的塊體位置結果當作能量法、潘家錚法的計算條件，能量法和潘家錚法所需要的其他計算參數見圖4及計算參數說明，速度方向為沿著坡面方向。從圖中可以看出，離散單元法的結果與能量法和潘家錚法一致性較好，這說明程序中塊體與邊界面的接觸判斷、接觸力的計算是正確的。同時可以在計算過程中看出，能量守恒法通常只能得到塊體在某一具體位置時的滑速，而要得到任意時刻的塊體速度則計算較繁瑣；潘家錚法則需要將整個計算時段劃分成很多小的時段以精確地模擬加速度的變化，采用手工計算時計算量較大。

圖5 各方法的滑坡體速度比較(c = 2 kPa, φ = 10°)Fig.5 Comparative results of block velocity (c = 2 kPa, φ = 10°)

5 涌浪的計算

當失穩后的滑坡體進入水庫中時，將產生涌浪。涌浪以滑坡體入水處為源點，先向對岸傳播，然后通過波的擴散和反射，引起一個向上游和下游傳播的過程[11]。要評價庫岸滑坡涌浪對附近建筑物的危害性，首先需要估算出涌浪到達各建筑物處的浪高。然而，這是一個相當復雜的過程，不容易求得精確解。在實際工程應用中，估算法用得較多，廣泛采用的有 Noda[13]提出的推算法和潘家錚[11]提出的近似估算法。

5.1 Noda法

Noda法[13]假定滑坡體滑落于半無限的水體中，且下滑高程大于水深。根據重力表面波線性理論，推導出滑坡引起的涌浪高度計算公式。其最大涌浪高度的表達式為[13]

式中：hmax為最大涌浪高度；h為水深；vs為滑坡體進入水庫時的速度。

從式（8）中可以看出，Noda法估算最大涌浪高度時主要考慮的因素是滑坡滑速，兩者呈線性關系。此外，最大涌浪高度還與水深有一定關系，水越深，滑坡體所激起的涌浪高度越高。

5.2 潘家錚法

潘家錚[11]認為，對于各種類型的滑坡只要確定其岸坡變形速率和過程以及反射特性后，均可用數值法求出涌浪過程。其基本假定為：①涌浪首先在滑坡入水處發生，產生初始波，然后向周圍傳播，假定傳播過程中涌浪能量損耗為已知；②不考慮邊界條件的非線性影響，整個涌浪過程可以視為在源點處產生的一系列小波的線性疊加；③每個小波都是孤立波，且以常波速在水面上傳播；④庫岸對涌浪的反射系數為已知值；⑤水庫庫岸為兩條平行陡壁，滑坡范圍內的庫岸斷面一致，滑坡變形率dA/dt或滑速v(t)為常數。

根據以上假定，可求得涌浪正對岸某點A的最大涌浪高度為

式中：ζmax為最大涌浪高度；ζ0為初始波高；k為波的反射系數；l為滑坡半長；B為滑坡體入水口至對岸點的距離。

波速c可按下式計算：

式中：ζ為波高。

初始波高ζ0可以通過下面的方法得到，對于水平岸坡變形情況：對于垂直岸坡變形情況：當時，當時，ζ0/λ呈曲線變化；當時， ζ0/λ =1；λ為滑坡的厚度。

在國內，潘家錚法用得較多，但潘家錚法只能給出某點處某一時刻所能達到的涌浪最大高度，而不能給出任意時刻涌浪的高度及形態。鑒于塊體的規模比水庫的庫容小得多，為了更好地模擬涌浪的形態，本文采用潘家錚法計算滑坡入水口處的初始涌浪，然后將其作為初始條件直接用二維波動方程求解涌浪在各個時刻的形態，邊界條件為如圖4所示的水庫邊界，其涌浪值一開始為 0，后續的值通過差分計算得到，這樣可以模擬水波的反射。含黏性阻尼的小振幅二維波動方程為[14]

式中：z為波高；c1為波速；μ為流體的黏度。采用顯式差分格式進行求解。

5.3 算例考核

計算參數同4.5節的算例，取圖4中的A、B、C、D為涌浪計算的特征點。

圖6中給出了一個塊體滑入水庫時所產生的涌浪在特征點處的涌浪高度的歷程曲線。從圖可以看出，特征點A處最早達到第1次涌浪峰值為0.098 m（在6.09 s時），因為A點的x坐標為0 m，距離塊體入水口的x坐標4 m最近，但同時可以看出，特征點A處在前10 s的最大涌浪高度是在7.29 s時的0.187 m，這說明涌浪在水庫某處取得最大值不一定在滑坡體剛失穩時，而有可能在涌浪反射后疊加時取得。

圖6 一個塊體失穩時特征點處涌浪高度-時程曲線Fig.6 Relationships of surge height and time of characteristic points after one block failure

表1中給出了一個塊體失穩時潘家錚法與本文的涌浪計算結果。從表中可以看出，本文的涌浪計算結果與潘家錚法計算結果較接近，能滿足實際需求。另外，本文中的涌浪計算可獲得任意時刻的涌浪計算結果，比潘家錚法更方便。

表1 涌浪計算結果比較Table 1 Comparative results of surge

圖7中給出了3個塊體滑入水庫時所產生的涌浪在特征點處的涌浪高度的歷程曲線。與圖6相比，由于涌浪的疊加，3個塊體所激起的涌浪高度比 1個塊體所激起的涌浪高很多，這說明在考慮滑坡所激起的涌浪時要考慮到多塊體下落的疊加效應，涌浪的最大值很可能出現在幾個涌浪值不大的塊體失穩后的疊加時刻。

圖7 3個塊體失穩時特征點處涌浪高度-時程曲線Fig.7 Relationships of surge height and time of characteristic points after three blocks failure

圖8、9分別給出了1個塊體和3個塊體失穩后7 s時的涌浪位置及形態（由于涌浪值較水庫的尺寸較小，圖中顯示涌浪高度值時采用的放大比例為5倍）。從圖中可以看出，結合離散單元法所提供的滑坡體的運動速度等基礎數據，采用潘家錚法和波動方程相結合對涌浪進行估算能給出滑坡體所激起的較合理的涌浪形態，為滑坡失穩所造成的災害評估提供了理論依據。

6 滑坡堆積計算

大型滑坡失穩后會堵塞河道，從而影響航運的正常運行，如果滑坡離電站進水口較近時，還會影響電站的正常工作。為了更好地對滑坡災害進行評估，對滑坡的運動軌跡及堆積形狀的模擬是至關重要的。

圖8 1個塊體失穩后產生的涌浪 (t=7.0 s)Fig.8 Surge shape after one block failure (t=7.0 s)

圖9 3個塊體失穩后產生的涌浪 (t=7.0 s)Fig.9 Surge shape after three blocks failure (t=7.0 s)

6.1 基本原理

離散單元法基于牛頓第二定律進行求解，每次求解前需要根據前一時刻給出的塊體位置進行接觸判斷，然后求得本時步各個塊體的位置信息。根據離散元所求解的基本公式（2）可以直接得到塊體的速度，然后對時步積分就可得到塊體的位置，這樣就得到了塊體的堆積形狀。

6.2 算例考核

算例在xz平面的幾何尺寸同4.5節中的圖4，y方向寬20 m；塊體尺寸為：x方向長4 m，y方向寬為4 m，z方向高為5 m。計算參數為：時間步長Δt = 10-5s，總時間步長為 1.5×106步，邊界面的法向剛度 Kn= 1010N/m3，切向剛度 Ks=4.0 ×109N/m3，塊體重度 γ= 20.0 kN/m3，接觸面上的黏聚力為2 kPa，內摩擦角為10°，黏性阻尼常數α=1.0。不同時刻塊體的堆積形狀如圖10所示。從圖中可以看出，離散單元法能給出滑坡運動后較合理的堆積形狀。

7 結 論

（1）在水庫庫岸滑坡涌浪的分析中，失穩后滑坡體的滑速及規模是分析的基礎，其值準確與否的關鍵在于如何正確獲得滑坡體的運動速度、運動軌跡及堆積形狀等數據，離散單元法可以直接得到上述數據。

圖10 不同時刻塊體的堆積形狀Fig.10 Heap shape of blocks at different times

（2）離散單元法不僅能達到與能量法和潘家錚法同樣的精度，而且計算滑坡體的滑速具有更方便、更高效的特點，且能獲得任意時刻塊體的位置與滑速。

（3）當滑坡體滑入水庫中時，結合離散單元法所提供的滑坡體的運動速度等基礎數據，并采用涌浪估算的潘家錚法和求解波動方程能給出滑坡體所激起的較合理的涌浪形態。

（4）從涌浪的計算結果來看，在涌浪分析中要特別注意涌浪的疊加現象，涌浪的最大高度很可能出現在幾個規模不大的滑坡體失穩后的疊加時刻。

（5）由邊坡的動態破壞過程可知，邊坡在失穩前一般都會表現出一定的預兆，如沿坡面向外的位移或速度的持續增大。因此，對邊坡加強安全監測和動態跟蹤是非常必要的。

（6）由于初步研究，涌浪的計算主要是根據經驗公式得到的，若要更好地模擬滑坡入水后涌浪的形狀、衰減規律和反射情況應采用離散單元法與計算流體力學耦合的方法，這有待于進一步的研究。

[1]黃潤秋. 20世紀以來中國的大型滑坡及其發生機制[J].巖石力學與工程學報, 2007, 26(3): 433－454.HUANG Run-qiu. Large-scale landslides and their sliding mechanisms in China since the 20th century[J]. Chinese Journal of Rock Mechanics and Engineering, 2007,26(3): 433－454.

[2]鄭穎人, 趙尚毅, 鄧衛東. 巖質邊坡破壞機制有限元數值模擬分析[J]. 巖石力學與工程學報, 2003, 22(12):1943－1952.ZHENG Ying-ren, ZHAO Shang-yi, DENG Wei-dong.Numerical simulation on failure mechanism of rock slope by strength reduction FEM[J]. Chinese Journal of Rock Mechanics and Engineering, 2003, 22(12): 1943－1952.

[3]汪洋, 殷坤龍. 水庫庫岸滑坡初始涌浪疊加的攝動方法[J]. 巖石力學與工程學報, 2004, 23(5): 717－720.WANG Yang, YIN Kun-long. Perturbation method of superposing initial surge height of landslide along reservoir shoreline[J]. Chinese Journal of Rock Mechanics and Engineering, 2004, 23(5): 717－720.

[4]CUNDALL P A. A computer model for simulating progressive, large-scale movements in blocky rock systems[C]//Proceedings of the International Symposium on Rock Mechanics. Nancy: [s. n.], 1971.

[5]ITASCA Consulting Group, Inc. 3DEC－3 Dimensional distinct element code (version3.0), user’s manual[R].Minneapolis: ITASCA Consulting Group, Inc., 2007.

[6]焦玉勇, 葛修潤, 劉泉聲, 等. 三維離散單元法及其在滑坡分析中的應用[J]. 巖土工程學報, 2000, 22(1): 101－104.JIAO Yu-yong, GE Xiu-run, LIU Quan-sheng, et al.Three-dimensional discrete element method and its application in landslide analysis[J]. Chinese Journal of Geotechnical Engineering, 2000, 22(1): 101－104.

[7]LEMOS J A. A distinct element model for dynamic analysis of jointed rock with application to dam foundation and fault motion[D]. St. Paul: University of Minnesota, 1987.

[8]周曉青, 王元漢, 哈秋舲. 離散單元法與邊界單元法的外部耦合計算[J]. 巖石力學與工程學報, 1996, 15(3):231－235.ZHOU Xiao-qing, WANG Yuan-han, HA Qiu-ling.External coupling computation between distinct element method and boundary element method[J]. Chinese Journal of Rock Mechanics and Engineering, 1996,15(3): 231－235.

[9]王泳嘉, 邢紀波. 離散單元法及其在巖土力學中的應用[M]. 沈陽: 東北工學院出版社, 1991.

[10]張楚漢. 巖石和混凝土離散-接觸-斷裂分析[M]. 北京:中國水利水電出版社, 2008.

[11]潘家錚. 建筑物的抗滑穩定與滑坡分析[M]. 北京: 水利出版社, 1980.

[12]SCHEIDEGGER A E. On the predication of the reach and velocity of catastrophic landslides[J]. Rock Mechanics,1973, 5(4): 231－236.

[13]NODA E. Water waves generated by landslides[J].Journal of Waterways, Harbors and Coastal Engineering Division, ASCE, 1970, 96(4): 835－855.

[14]龐明勇. 基于離散模型的二維水波實時動態模擬方法[J].水利學報, 2007, (11): 1358－1363.PANG Ming-yong. Realtime dynamic simulation of 2D water waves based on discrete model[J]. Journal of Hydraulic Engineering, 2007, (11): 1358－1363.