云南黃坪庫區(qū)滑坡運動及其失穩(wěn)模式的離散元模擬

蔣 權(quán)，陳希良，肖江劍,邱向東,李攀峰

(1.中國科學(xué)院寧波材料技術(shù)與工程研究所，浙江 寧波 315201；2.中國科學(xué)院大學(xué)，北京 100049;3.中國電建集團成都勘測設(shè)計研究院有限公司，四川 成都 610072)

0 引言

庫區(qū)滑坡失穩(wěn)破壞后，滑坡體入水激起巨浪而產(chǎn)生滑坡涌浪，若激起的涌浪足夠強還可能會造成壩頂過水甚至沖毀水工建構(gòu)筑物，堵塞河道、威脅船只航行及沿岸居民生命財產(chǎn)安全[1]。特別是正在蓄水或蓄水不久的庫區(qū)更容易發(fā)生此類地質(zhì)災(zāi)害。一方面，水庫的修建會不可避免地改變庫區(qū)岸坡原有的自然地質(zhì)條件，如庫岸巖土邊坡的物理性質(zhì)，進而影響邊坡的穩(wěn)定。另一方面，庫水會對庫岸已存在的不穩(wěn)定地質(zhì)體和滑坡體、崩塌體產(chǎn)生浸潤和托浮作用，其變化往往對庫岸坡體的穩(wěn)定性進一步造成影響，在受到暴雨、庫水位突然變化等外界影響因素的作用下，一些處于不太穩(wěn)定狀態(tài)的庫岸滑坡體就可能失穩(wěn)[2]。

極端氣候異常，工程擾動增強，地質(zhì)條件復(fù)雜等，是地質(zhì)災(zāi)害防治工作需要面對的嚴峻挑戰(zhàn)[3]。針對滑坡體失穩(wěn)后滑坡體運動的問題，相關(guān)研究在許多國家的不同學(xué)科和領(lǐng)域都在開展。2002年1月成立的國際滑坡協(xié)會及每3年舉辦一次的世界滑坡論壇更是大力推進了對滑坡運動的研究[4-5]。在分析滑坡運動時，對不同滑坡體物理形態(tài)有不一樣的方法。就小塊且松散堆積的滑坡體而言，許多學(xué)者用離散單元法、非連續(xù)變形法、顆粒流等非連續(xù)介質(zhì)數(shù)值方法對坡體滑動過程進行了模擬分析，并取得不少成果[6-8]。但目前的一些工作通常有很多簡化，如將滑坡體簡化成簡單的塊體，或者只研究二維情況下的離散坡體滑動情況，這難以真實刻畫三維復(fù)雜的滑坡體力學(xué)特性，距實際工程應(yīng)用較遠。

本文以溪洛渡水電站黃坪滑坡為例，結(jié)合現(xiàn)場地質(zhì)調(diào)查，對滑坡特征進行了詳細的介紹，并分析了滑動成因和機制；在不考慮空氣與江水作用的情況下，利用離散元對三維離散滑坡體的滑坡失穩(wěn)過程進行了模擬，并詳細的分析了其各階段的運動特征。

1 滑坡概況及成因機制分析

1.1 滑坡概況

2013年7月27日下午，云南省昭通市永善縣黃華鎮(zhèn)黃坪村大石包附近發(fā)生滑坡?；挛挥诮鹕辰野队郎瓶h黃華鎮(zhèn)老黃坪集鎮(zhèn)上游約1 km大石包附近，即黃碼公路樁號K42+000～K42+200 m段。滑塌范圍：順江方向長約200 m，江水位高程554 m以上約110 m，后緣至黃碼公路附近高程約660 m，部分區(qū)域高于路面50 m范圍，滑體厚度3～5 m，滑坡體積數(shù)十萬立方米?；聟^(qū)位于預(yù)測滑塌區(qū)范圍內(nèi)的上游部分，現(xiàn)場及斷面示意分別見圖1、圖2。

圖1 黃坪滑坡全貌Fig.1 The panorama of the Huangping landslide

圖2 黃坪滑坡剖面示意圖Fig.2 Cross section of the Huangping Landslide

該岸坡地形較陡，天然坡度在34°～38°。堆積體下部地層巖性為奧陶系上統(tǒng)與志留系O3+S(S)灰色、深灰色、灰黃色砂巖、粉砂巖、泥巖、頁巖、砂質(zhì)頁巖夾泥灰?guī)r、泥質(zhì)灰?guī)r，巖性較軟弱，產(chǎn)狀N30°～ 45°E/SE∠17°～20°；上部陡坡地形為陽新灰?guī)r；堆積體為崩坡積與沖洪積形成的塊碎石土層，塊碎石粒徑較小，一般為5～20 cm，表部含量較高約40%～50%，下部含量較低約25%～35%，細粒部分主要由細粉砂與黏粒組成，在堆積體表層可見鈣質(zhì)膠結(jié)的硬殼。

1.2 成因機制分析

據(jù)現(xiàn)場調(diào)查和已有的資料分析，黃坪滑坡形成的原因主要有以下兩點：(1)坡體物質(zhì)組成、結(jié)構(gòu)、坡度等自然條件造成該處岸坡天然狀態(tài)下穩(wěn)定性較差。岸坡下部地層物質(zhì)組成決定了其較軟弱的巖性；該部位岸坡地形較陡，中間被一條較小的無水干溝切割成上、下游兩部分，岸坡巖體結(jié)構(gòu)完整性差；岸坡由于受長期的水流沖刷作用，為表層破碎巖體各類重力地質(zhì)現(xiàn)象的產(chǎn)生提供了有利條件。(2)水庫蓄水是黃坪滑坡的重要誘發(fā)因素，加速了該處庫岸失穩(wěn)。滑坡發(fā)生前水庫已經(jīng)開始蓄水，岸坡前緣坡腳處于蓄水位以下，由于水的軟化弱化作用，前緣坡腳巖土體的強度被進一步降低，同時也導(dǎo)致其有效應(yīng)力的減小。

根據(jù)滑坡概況及成因分析其滑動機制為：(1)庫水位的快速上升，導(dǎo)致了中下部坡體的軟化，抗剪強度降低；水位線抬升產(chǎn)生的浮托力還減小了岸坡前緣巖土體的有效應(yīng)力導(dǎo)致抗滑能力降低。(2)在上部坡體的重壓下，水下坡體發(fā)生突發(fā)式破壞，從而導(dǎo)致庫水位上下坡體快速失穩(wěn)，墜入庫內(nèi)。

2 離散元模擬及模型參數(shù)

2.1 離散元數(shù)值模擬方法

離散元法的基本原理是基于離散化的思想，對研究對象進行單元劃分，離散化的各個單元之間相互獨立。根據(jù)單元之間的相互作用和牛頓運動定律，采用動態(tài)松弛法和靜態(tài)松弛法等迭代方法，對各單元進行循環(huán)迭代計算，得出每一個時間步長內(nèi)各單元的受力及位移。其控制方程主要包括平動與轉(zhuǎn)動兩方面，具體形式如下:

(1)

(2)

2.2 數(shù)值模型

基于黃坪滑坡的實際數(shù)據(jù)，以黃坪滑坡為中心選取了1 280 m×200 m×416 m 模擬計算區(qū)域，即模型底部邊界長1 280 m，兩側(cè)邊界高為416 m，沿江長200 m?；诨缕拭娼⒘它S坪滑坡所在區(qū)域的離散元計算模型，模型分為滑坡體、基巖滑床、滑坡邊界三部分。為監(jiān)測滑坡體隨時間的變化，在滑坡體中心取截面A，并在截面A上設(shè)置11個監(jiān)測點：后緣和前緣各分布1個點，滑坡后部、中部以及前部按照不同深度在表層、中層、底層分別設(shè)置了監(jiān)測點。由于模型較大，取右邊關(guān)鍵局部示意圖見圖3。

圖3 黃坪滑坡離散元計算模型Fig.3 Discrete element model of the Huangping Landslide

在離散元模擬中，滑體采用離散顆粒，滑床及滑坡邊界采用同滑體一樣參數(shù)的不考慮變形的剛體。對于土質(zhì)滑坡和破碎巖石滑坡來說，要想得到和實際較為一致的模擬結(jié)果，可將滑坡運動看成顆粒流運動，也可將滑坡體等效為擬流體來進行數(shù)值模擬[9]。本文用正態(tài)分布的不同大小圓球顆粒表征破碎巖土，并組成仿真中所需要的滑坡體，滑坡體體積約為3.0×105m3，組成滑坡體的顆?？倲?shù)近4.0×105顆。

離散元法的接觸模型有多種，接觸力的計算方法也各不相同，但是整體計算的原理都是相同的。本文中，顆粒間的接觸力采用Herz-Mindlin無滑移接觸模型求解[10]。其中，模型中法向接觸力的求解基于Hertzian 接觸理論，而切向接觸力的求解基于Middlin-Deresiewicz 接觸理論[11]。

2.3 計算參數(shù)

計算中用到的滑坡體與滑床或滑坡邊界的材料屬性、相關(guān)接觸屬性及運行參數(shù)見表1。

表1 黃坪滑坡涌浪計算參數(shù)

3 滑坡運動模擬結(jié)果與分析

3.1 滑坡運動過程的整體特征

從滑坡體外部出發(fā)：首先，取滑坡體失穩(wěn)后運動過程中典型時刻的滑坡體表面形態(tài)見圖4，可知滑坡體表面在運動過程中由凸形態(tài)逐漸轉(zhuǎn)為凹形態(tài)。

圖4 不同時刻黃坪滑坡運動形態(tài)Fig.4 Movement states of landslide at different time

然后，從滑坡體上方俯視向下看可見滑坡失穩(wěn)后的運動過程中，滑坡體后緣邊界線趨于凹形態(tài)(圖5)。

圖5 滑坡后緣邊界運動Fig.5 Trailing edge movement of landslide

因此滑坡運動過程中，其表面形態(tài)會發(fā)生轉(zhuǎn)變。而滑坡體后緣邊界線趨于凹形態(tài)是由于失穩(wěn)滑坡體兩側(cè)受到來自周圍相鄰未失穩(wěn)滑體的支撐作用，滑坡體呈現(xiàn)出中間失穩(wěn)體積較大，兩側(cè)失穩(wěn)體積較小的特征。這與一般滑坡現(xiàn)場所表現(xiàn)出的形式一致。

在滑坡體失穩(wěn)后的滑動過程中，整個滑坡體的平均速度、最大速度及最小速度如圖6所示。

圖6 滑坡運動速度Fig.6 The velocity of landslide movement

從圖6中平均速度曲線的趨勢可知滑坡體運動的大致情況：0～5 s為滑坡啟動后滑坡體整體運動的加速階段，其中加速較快的是后緣及后部受壓力與阻力較小的滑坡體，這可從壓力分布圖7中看出；5～15 s，在中部前部凸形滑坡體顆粒阻力的作用下，滑坡體整體開始減速；15～25 s，在后緣及后部滑體沖擊與不斷堆積作用下，滑坡體整體再次緩慢加速，并開始滑向坡底；25 s以后，滑坡體在前方坡底堆積的滑坡體顆粒阻力作用下開始緩慢減速運動，直到70 s，滑坡體趨于穩(wěn)定并停止運動。

從圖6中最大速度曲線可以看到：在啟動時刻，有局部滑坡體顆粒達到12 m/s的速度；在高速滑動階段，有局部滑坡體顆粒達到32 m/s的最大速度；高速滑動階段過后，滑坡體以似整體顆粒流形式流動，滑坡體顆粒的最大速度在緩慢減小。

從圖6中最小速度曲線可看到，在滑坡啟動、高速滑動及自穩(wěn)堆積階段，存在速度較小甚至沒有運動的滑坡體顆粒；而在整體流動階段，滑坡體顆粒幾乎都在運動。

圖7 A截面上的壓力分布Fig.7 Compressive force distribution on cross section A

3.2 滑坡體內(nèi)部運動特征

為深入的分析滑坡運動特征，取滑坡體內(nèi)截面A的8個典型運動時刻來對滑坡體內(nèi)部顆粒分布、速度分布及壓力分布做出分析。

3.2.1滑坡體內(nèi)顆粒分布特征

如圖8不同時刻截面A上不同大小滑坡體顆粒分布情況可知：在滑坡運動啟動前或初始階段，滑坡體1大顆粒(灰色)分布在滑坡體底面，滑坡體2小顆粒(橙色)分布在滑坡體表面。隨著滑坡運動的進行，顆粒之間碰撞與振動加劇，部分小顆粒慢慢運動到大顆?；麦w1內(nèi)部及底部，大顆粒運動至滑坡體表面，而小顆粒會運動到滑坡體下層。這說明，在滑坡體運動過程中產(chǎn)生了“巴西果效應(yīng)”，小顆粒通過大顆粒間的縫隙逐漸往下層滲流并填充到大顆粒下部空穴中，而大的滑坡體顆粒則移動到上部并在一定程度上被阻止向下運動[13]。

圖8 A截面上顆粒分布Fig.8 Particles distribution on cross section A

3.2.2滑坡體速度分布特征

與將滑坡體看成剛體不同的是，離散滑坡體顆粒的速度在不同位置是不同的，且速度的變化是分階段進行的。典型時刻截面A上滑坡體顆粒的速度分布見圖9。

圖9 A截面上的速度分布Fig.9 Velocity distribution on cross section A

從圖9中可看到滑坡運動啟動-高速滑動-碰撞-滑至坡底-堆積自穩(wěn)的典型階段：0～20 s，滑坡運動啟動后，以中后部為主的高速滑動區(qū)滑坡向中前部運動，并發(fā)生碰撞堆積；20～40 s，滑坡體在前面階段的沖擊與堆積重力作用下，開始以似整體顆粒流的形式滑向坡底；40～60 s，進一步滑向坡底的滑坡體顆粒發(fā)生碰撞堆積，并進行自穩(wěn)過程；在70 s時完成自穩(wěn)過程，滑坡體整體停止運動。值得注意的是，在20 s之后，滑體似整體顆粒流的形式滑向坡底，其整體性逐漸變好，速度分布也較為均勻。

3.2.3滑坡體內(nèi)壓力分布特征

滑坡體內(nèi)壓力的分布是隨著運動的變化而變化的。從圖7 A截面上不同時刻壓力分布圖中可看到：0～20 s，滑坡運動啟動時中后部滑坡體所受壓力較小，在初始動能和勢能作用下開始向中前部運動，并發(fā)生碰撞堆積，滑坡體內(nèi)顆粒所受壓力累積增大；20～40 s，中前部滑坡體在中后部滑坡體的沖擊與擠壓作用下滑向坡底，此時滑坡體內(nèi)部壓力開始釋放減?。?0～60 s，進一步滑向坡底的滑坡體顆粒發(fā)生碰撞堆積，隨著顆粒逐步的聚集，滑坡體內(nèi)顆粒所受壓力再一次累積增大；到70 s時，滑坡體內(nèi)部壓力分布不再變化。另外，與連續(xù)性剛體不同的是，由于滑坡體顆粒具有離散性與不均勻性，滑坡體內(nèi)壓力的分布是不均勻的，且內(nèi)部顆粒所受壓力也不是在最厚重的地方就一定最大。

3.3 滑坡內(nèi)部監(jiān)測顆粒運動的速度與距離

為更進一步分析滑坡體內(nèi)部顆粒運動特征，對截面A上的11個典型監(jiān)測點顆粒的運動進行了研究。值得注意的是，雖然在滑坡運動過程中監(jiān)測顆粒的相對位置會發(fā)生改變，但還是能反映出一些滑坡體的運動特點。不同時刻典型監(jiān)測顆粒的運動速度及運動距離分別見圖10和圖11。

從圖10(a)、(b)、(c)、(d)中的曲線1～11可知，在0～5 s滑坡運動啟動后的加速階段，典型監(jiān)測顆粒的加速度從后緣、后部、中部、前部到前緣，表層、中層至底層大致由大到小分布。這可從圖7 A截面上的壓力分布得到解釋，初始階段滑坡體后緣表層所受壓力及阻力相對較小，其加速度也較大。因此滑坡體內(nèi)同一橫剖面不同深度的滑坡體顆粒，其運動速度大致表現(xiàn)出離表層越近速度越大的特點；在5 s之后，滑坡體內(nèi)同一橫剖面不同深度的滑坡體顆粒速度關(guān)系發(fā)生了變化，這是因為顆粒位置及受力情況都發(fā)生了改變，這也可從圖8 A截面上的顆粒分布圖和圖7 A截面上的壓力分布圖得到解釋。

圖10 滑坡監(jiān)測點速度Fig.10 Velocity of monitor particles of landslide

圖11 滑坡監(jiān)測點中心運動距離Fig.11 Distance of monitor particles of landslide

從圖10(a)中曲線1看到前緣顆粒速度在25～38 s的時候有較大不連續(xù)的波動，這說明它正經(jīng)過不同曲率的滑坡邊界運動至坡底。同樣地，由圖11曲線2～10可知：前部、中部滑坡體顆粒大約分別在25 s、30 s經(jīng)過不同曲率的邊界開始向坡底運動；后部、后緣滑坡體顆粒則在35 s、40 s經(jīng)過滑床底部拐點運動至滑坡邊界上。

在圖10(b)中，原來位于前部表層的顆粒4的速度在55 s后處在了原來位于前部中層的顆粒3下面；圖10(c)中，原來位于中部表層的顆粒5的速度在55 s后處在了原來位于中部底層的顆粒7下面。這說明滑坡體內(nèi)部顆粒位置的確發(fā)生了很大的變化，滑坡體運動無法以完全的整體式顆粒流運動，這從圖8滑坡體截面上顆粒的分布變化也可以說明。

整個滑坡運動過程中監(jiān)測顆粒1～11的運動距離分別為：150.643 m、146.576 m、125.662 m、171.64 m、217.484 m、203.217 m、186.862 m、229.388 m、227.097 m、240.755 m、277.274 m。其中，后緣監(jiān)測顆粒11的運動距離最大，即277.274 m。從圖11中的曲線1～11可知，在0～5 s滑坡啟動后的加速階段，監(jiān)測顆粒的運動距離從滑坡后緣至前緣由大到小分布，從滑坡體同一橫剖面的表層至底層由大到小分布，即隨著滑坡體深度的增加，滑動距離逐漸減小；5～15 s，依舊符合運動距離從后緣到前緣，表層至底層由大到小分布的規(guī)律。但滑坡體內(nèi)部顆粒位置已經(jīng)發(fā)生了很大變化，從圖中曲線2、3就可以看到，前部中層監(jiān)測顆粒3相對于前部底層監(jiān)測顆粒2運動到了滑坡體內(nèi)更深的位置。

4 結(jié)論

本文總結(jié)了黃坪滑坡的基本特征，利用離散元方法對滑坡過程進行了模擬。通過由表及里，從整體到局部，深入細致的分析滑坡啟動-高速滑動-碰撞-滑至坡底-堆積自穩(wěn)運動過程中的主要特點后，得到了如下結(jié)論：

(1)滑坡體表面在運動過程中會發(fā)生較大變化，且滑坡體會呈現(xiàn)出中間失穩(wěn)體積較大，兩側(cè)失穩(wěn)體積較小的特征。

(2)在滑坡啟動、高速滑動及自穩(wěn)堆積階段，滑坡體內(nèi)存在速度較小甚至沒有運動的滑坡體顆粒，當(dāng)滑坡體處于整體流動階段后，顆粒幾乎都在運動。

(3)與將滑坡體看成剛體不同的是，離散滑坡體運動中：滑坡內(nèi)部顆粒分布會發(fā)生改變，存在大顆粒運動至滑坡體表面，小顆粒運動到滑坡體底面的現(xiàn)象；顆粒的速度和壓力在不同位置也會發(fā)生改變，且存在不同顆粒間速度相差很大，滑坡體內(nèi)壓力分布不均勻的情況。

(4)滑坡體顆粒的運動速度、運動距離表現(xiàn)出從后緣到前緣，同一橫剖面中表層至底層由大到小分布的規(guī)律。