徐 寅 ，陳勝宏

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

1 引 言

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

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

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

2 離散元軟件開發(fā)

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

3 離散單元法的基本原理

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

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

3.1 接觸判斷

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

3.2 塊體運動

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

圖2 計算循環(huán)[5]Fig.2 Calculation cycle[5]

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

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

用中心差分法可得

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

4 滑速的計算

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

4.1 能量守恒法

如圖3所示，根據(jù)動能定理有

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

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

又有

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

求得

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

4.2 潘家錚法

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

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

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

4.3 Scheidegger法

由于內(nèi)摩擦系數(shù)值在實際應用中不易選取，Scheidegger[12]在搜集了大量滑坡數(shù)據(jù)的基礎(chǔ)上，進行了統(tǒng)計分析，發(fā)現(xiàn)滑坡體的等效內(nèi)摩擦系數(shù)與滑坡體體積間有一定相關(guān)關(guān)系，給出的經(jīng)驗公式為[12]

式中：f為等效內(nèi)摩擦系數(shù)；V為滑坡體體積。

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

4.4 離散元法

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

4.5 算例考核

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

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

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

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

5 涌浪的計算

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

5.1 Noda法

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

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

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

5.2 潘家錚法

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

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

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

波速c可按下式計算：

式中：ζ為波高。

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

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

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

5.3 算例考核

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

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

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

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

表1 涌浪計算結(jié)果比較Table 1 Comparative results of surge

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

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

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

6 滑坡堆積計算

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

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

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

6.1 基本原理

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

6.2 算例考核

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

7 結(jié) 論

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

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

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

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

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

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

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

[1]黃潤秋. 20世紀以來中國的大型滑坡及其發(fā)生機制[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]鄭穎人, 趙尚毅, 鄧衛(wèi)東. 巖質(zhì)邊坡破壞機制有限元數(shù)值模擬分析[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]潘家錚. 建筑物的抗滑穩(wěn)定與滑坡分析[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]龐明勇. 基于離散模型的二維水波實時動態(tài)模擬方法[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.