基于PFC3D的滑坡堰塞壩堆積過程與形態模擬

吳 建 川,張 世 殊，吳 爽，冉 從 彥

(1.中國電建集團 成都勘測設計研究院有限公司，四川 成都 610072; 2.中國地質大學(武漢)工程學院，湖北 武漢 430074)

滑坡堰塞壩是指由于降雨、地震、火山噴發等原因引起山體滑坡而截堵山谷、河道形成的壩體[1]。依據成因的不同，堰塞壩可以劃分為滑坡型堰塞壩、崩塌型堰塞壩以及滑坡、崩塌運動過程中形成的泥石流型堰塞壩，其中滑坡型堰塞壩約占全部堰塞壩的70%以上[2]。在西部高山峽谷區域的河道，滑坡失穩堵江極易形成較高的堰塞壩，堰塞壩的物質組成通常由松散堆積物組成，在河水漫頂、雨水沖刷的條件下極易發生潰壩并引發洪水，對水電工程的選址、建設、安全運行和下游居民的生命財產安全造成巨大威脅。如2018年10月10日、11月3日，西藏自治區境內的金沙江白格村在同一地點先后兩次發生大規模的山體滑坡，堵塞金沙江形成堰塞湖，嚴重威脅著堰塞壩上下游人民群眾的生命財產和下游已建及在建水電站的安全[3]。

滑坡型堰塞壩的形成過程和幾何形態特征與河道縱剖面形狀、滑坡體積以及滑坡體的運動形態有關[4]。為了探尋滑坡型堰塞壩的形成過程和幾何形態特征，成本低、效果明確的數值模擬是一種有效的手段。連續力學數值模擬方法方面：原俊紅[5]、羅剛[6]運用連續力學數值模擬方法，模擬邊坡的破壞特征并預測堰塞壩高度；徐文杰[7]對肖家河滑坡穩定性做了分析，并用ABAQUS/Explicit動力有限元分析滑坡在地震工況下的堵江過程。非連續力學數值模擬方法方面：Chang和 Taboada[8]利用二維顆粒流離散元程序模擬了臺灣九份二山滑坡力學過程及堆積行為；孟云偉和柴賀軍[9]認為顆粒流離散元能充分考慮實際巖體的大位移、大變形，能夠更有效、更真實地模擬滑坡運動的優點；張龍[10]認為PFC3D軟件能較好地模擬高速遠程滑坡的運動過程及堆積情況，對于災害堆積形態與影響范圍可以直接模擬得到；Poisel和Preh[11]利用PFC3D對滑坡運動導致的水庫涌浪進行了模擬，并取得非常好的效果；梁承洋[12]用PFC2D模擬了川藏交通廊道冰磧物滑坡堵江堆積物的形態，同時研究了單因素變量對堵江的影響；李子隆[13-14]用PFC2D開展了覆蓋層庫岸邊坡在靜水位條件下的失穩堰塞模擬和滑坡失穩形成堰塞壩的高度預測模擬；Dal Sasso[15]采用地貌指數法對滑坡堰塞壩幾何特征進行了研究；Zhao[16]等采用CFD-DEM耦合的方法模擬滑坡堵江形成堰塞壩的過程。針對滑坡堰塞壩幾何特征的研究，Dong[17]等指出準確的數字地形數據在絕大多數情況下無法短時間內獲得。

綜上，連續力學數值模擬方法雖有模擬速度快等突出優勢，但由于模擬的塊體不能破裂，模擬的效果并不明顯；而于非連續力學方法的顆粒流離散元方法能較好地模擬大變形、大位移，對模擬滑坡失穩堰塞過程及堰塞壩堆積形態具有較好的適用性。

針對堰塞壩的潰決，蔡耀軍[18]研究了白格堰塞湖形成過程及成因，并分析了堰塞體結構及形態特征、潰決發展階段、潰決特征值；陳祖煜[19]運用DB-IWHR[20]潰壩洪水分析程序反演分析了金沙江“10.10”白格堰塞湖漫頂自然泄流過程；王敏[21]采用MIKEl l模型和一維潰壩洪水模型(DBFM)，對金沙江“11·3”白格堰塞湖潰決洪水演進過程進行了對比計算。

本文以雅礱江上游水電規劃、開發的重點關注對象唐古棟滑坡為研究對象，運用三維顆粒流離散元PFC3D建立三維滑坡模型，對唐古棟滑坡1967年失穩形成的堰塞壩堆積形態進行反演模擬，校核標定的細觀參數。在此基礎上，開展唐古棟滑坡變形強烈部位失穩形成堰塞壩的預測模擬，并采用DAMBRK潰壩洪水計算數學模型對堰塞壩潰壩造成的影響進行了探討。

1 唐古棟滑坡

1.1 簡 介

唐古棟滑坡位于雅礱江中游雅江縣孜河鄉雨日村，位于雅礱江兩河口至卡拉中游規劃河段力邱河口-蒙古山段的右岸，為一巨型高速巖質滑坡，主滑方向145°，與雅礱江流向近垂直[22]。滑坡體橫河長約1.4 km，順河寬約1.2 km，展布面積超過1.4 km2，前緣高程2 450 ～2 480 m，后緣高程3 500 m，高差約1 000 m。滑坡區地層巖性為三疊系上統侏倭組(T3zh)變質砂巖與印支-燕山期侵入的花崗偉晶巖脈(γρ)，巖層產狀為N50°～75°W/NE∠40°～60°，傾向坡內[23]。

1.2 滑坡分區

唐古棟滑坡主要包括上游側變形體和下游側已滑坡區，根據岸坡變形破壞程度及殘留物規模，由上游至下游可分為A、B、C、D四個區(見圖1)，A區為未發生滑動的強變形區，B、C、D區為已發生滑動區，各區特征及規模見表1[24]。

1.3 滑坡體物理力學參數

唐古棟滑坡物質組成由外至內可分為滑坡堆積物、崩坡堆積體、強風化巖體、弱風化巖體、微-新巖體。基于滑坡巖土體物理力學試驗成果，并通過抗剪強度參數反演，綜合確定各層巖土體物理力學參數取值，見表2[25]。

表1 唐古棟滑坡分區特征及殘余方量匯總

表2 唐古棟滑坡各層巖土體抗剪強度參數取值

1.4 滑坡失穩-堰塞-潰決

1967年6月8日，唐古棟滑坡發生大規模整體滑移失穩破壞，且在5 min之內堵斷雅礱江并沖向對岸岸坡；滑坡體在高速滑動過程中充分解體，堆積形成順河向呈梯形的堰塞壩，堰塞壩方量約6.41×107m3；堰塞壩底面沿江長1 500 m,頂面沿江長860 m，高225～270 m；對岸壩高335 m，本岸壩高175 m[26]。

堰塞壩阻斷雅礱江達10 d，黃潤秋[27]等調查發現，堰塞壩回水53 km，壅高水位高程達2 575 m，庫容達6.8億 m3；堰塞壩潰決后，洪峰流量達5 3000 m3/s，壩下游10 km處水位上漲達48 m，洪水沿雅礱江而下，沖毀房屋、田地、公路、橋梁等，給下游造成巨大損失。

2 PFC3D顆粒流參數標定

2.1 PFC3D軟件

PFC3D是根據Cundall提出的細觀離散元理論(粒子流理論)開發的一種數值計算軟件。該軟件模擬單元為三維球體顆粒(granular)，從微觀結構的角度出發，研究顆粒集合體的破裂、破裂發展和大位移的顆粒流問題[28]。

模型的建立基于以下假設：① 顆粒單元被視為剛性體；② 顆粒間的接觸發生在可忽略不計的點上；③ 顆粒間的接觸視作柔性接觸，顆粒允許在接觸點上相互重疊；④ 重疊量與接觸力大小有關，接觸力通過力-位移定律計算，需要指出的是重疊量相比于顆粒尺寸非常小；⑤ 顆粒間的接觸點可以存在連接(bonds)；⑥ 所有顆粒的形狀都是球形，也可以用到clump生成任意的形狀，在每個clump中顆粒相互重疊但沒有相互作用力。

基本計算原理是運用顯示中心差分法進行動態松弛，計算過程是在每個時步內交替利用力-位移定律與牛頓第二定律，時步迭代并遍歷整個顆粒集合體(見圖2)。

圖2 PFC3D計算過程

材料的本構行為通過接觸點上的本構模型來模擬，顆粒間的接觸本構模型有3種：接觸剛度模型、滑動模型和接觸連接模型。

巖土體并非剛性體，系統的能量并不僅僅通過單元之間的摩擦方式進行耗散，顆粒與顆粒及顆粒與墻體之間的碰撞能耗也是系統能量損失的最重要途徑之一。為了更真實地還原巖土體介質在運動過程中的碰撞問題，PFC3D提供了局部阻尼、黏性阻尼兩類基本阻尼模型。局部阻尼的存在不適用于滑坡在重力作用下的動力學分析，模擬過程中將局部阻尼系數統一設為零；黏性阻尼通過模擬顆粒自由下落的“碰撞-反彈”，分析顆粒的速度恢復系數，確定出法向阻尼系數為0.4，切向阻尼系數為0.1。

2.2 參數標定

PFC3D模擬采用的是顆粒之間的微觀參數，主要包括：球最小半徑Rmin、球半徑比Rmax/Rmin、球-球接觸模量Ec、球剛度比kn/ks、球摩擦因數f、平行黏結半徑乘子λ、平行黏結模量Ec、平行黏結剛度比kn/ks、平行黏結法向強度與切向強度σ和τ。試驗獲取的巖土體參數是宏觀參數，巖土體的宏觀參數與微觀參數之間存在一定的聯系，但目前尚沒有成熟的理論和方法來界定這種聯系。針對PFC3D微觀參數的獲取，筆者運用軟件自帶的虛擬三軸測試程序進行參數標定與反演，通過將反演參數與唐古棟滑坡綜合物理力學參數進行對比，確定出微觀物理力學參數。

參數標定按照以下過程進行：① 首先根據實際方量與計算機計算能力確定顆粒最小半徑為3 m，最大顆粒與最小顆粒半徑比為1.67；② 根據崩坡積物的變形參數,如彈性模量與泊松比反演出對應的法向剛度與切向剛度為1×106N/m；③ 采用平行連接模型模擬崩坡積物力學特性，根據強度參數凝聚力與內摩擦角反演出連接強度，分別為法向強度20 MPa，切向強度15 MPa；④ 顆粒摩擦系數會影響最后的堆積體形態，試算不同摩擦系數最終導致的堆積體形態，并與1967年唐古棟形成的堆積形態進行對比，摩擦系數的增加會導致最終堆積壩堆積坡度的增高，通過多次調整，最終得到相似堆積壩堆積坡度對應的摩擦系數為0.324，并將其用于本次數值計算中。模型參數見表3。

2.3 堰塞壩堆積反演分析

根據1967年唐古棟滑坡堆積堰塞壩形態反演，以驗證標定細觀參數的可靠性。

2.3.1構建模型

根據A區地形恢復唐古棟滑坡1967年以前的地形，1967年主要失穩范圍B、C、D區的失穩方量約為9 260萬m3。通過提取一定密度的高程點，運用三維地質建模，將每3個點組合成一個三角形墻面，最后將所有墻面有序連接。地形由18 409個三角形組成，滑體由半徑3～5 m的171 343個球形單元組成，球與球之間設置平行連接，球與地表之間設置為接觸剛度模型中的彈性模型(見圖3)。

表3 唐古棟滑坡PFC3D數值模擬細觀參數

圖3 1967年失穩反演三維模型

2.3.2堰塞壩形態反演分析

在模擬進行到100 000時步時，滑坡運動達到收斂平衡，堆積形成堰塞壩(圖4)，分別沿河流流向截取一條橫剖面線，沿主滑方向截取一條縱剖面線，得到堰塞壩的橫剖面形態圖(圖5)和縱剖面形態圖(圖6)。

圖4 1967年失穩反演堰塞壩堆積形態

從堰塞壩橫剖面形態圖(圖5)中可以看出，堰塞壩在河道內呈中部高兩邊低的梯形形態分布，壩底寬為1 483 m，壩頂寬為516 m，最大高度178 m，下游坡度為18°，上游坡度為22°。

圖5 1967年失穩反演堰塞壩橫剖面形態(單位：m)

從堰塞壩縱剖面形態圖(圖6)中可以看出，堰塞壩在滑坡對岸一側高度明顯較高，在靠近滑坡一側高度較低，對岸壩高245m，本岸壩高178m。

圖6 1967年失穩反演堰塞壩縱剖面形態(單位：m)

對比分析反演堰塞壩與1967年堆積堰塞壩的形態，結果見表4。沿河流流向的橫剖面，堰塞壩底寬偏差1.13%，頂寬偏差40%，最大壩高偏差20.89%～34.07%；沿滑坡主滑方向的縱剖面，對岸壩高偏差26.87%，本岸壩高偏差1.71%。

通過上述對比，驗證出標定的細觀參數是合理的。1967年堆積堰塞壩的形態是通過調查殘余痕跡獲得，且壩頂寬殘余形態保留較差，模擬堰塞壩與堆積堰塞壩上述3個數據的對比偏差達到20%～40%。

表4 堰塞壩堆積形態對比

3 堰塞壩堆積過程與形態分析

根據對唐古棟滑坡的地質調查及穩定性分析，滑坡A區(見圖7)蠕滑體在地震工況下可能發生整體失穩，失穩方量達1 700萬m3，規模巨大，故需對A區失穩后堆積形成的堰塞壩開展預測分析，為潰決分析、應急治理等提供依據。

3.1 構建模型

構建三維模型，地形由19 572個三角形組成，滑體由半徑為3～5 m的35 109個球形單元組成，球與球之間設置為平行連接，球與地表之間設置為接觸剛度模型中的彈性模型(見圖8)。參數選取反演標定的細觀參數。

圖7 滑坡A區平面示意

圖8 A區三維模型

3.2 失穩堆積過程分析

根據失穩過程中不同時刻滑坡形態的改變及運動的特點，分別選取6 000,10 000,20 000時步截取位移云圖和縱剖面圖進行堆積過程分析。

3.2.16000時步堆積分析

6 000時步滑坡整體下滑約500～600 m，部分后緣滑體繞過基巖沿著滑坡邊界沖溝滑動，前緣部分巖土體有部分已滑入河谷處但是堆積方量不大。由于地形原因前緣滑體滑動速率較后緣滑體慢，大部分后緣及中部滑體集中于整個滑坡中部。滑動前期由于動能較大，此時巖土體移動速率較大(見圖9,10)。

圖9 6000時步位移云圖(單位：m)

圖10 6000時步縱剖面圖(單位：m)

3.2.210000時步堆積分析

10 000時步滑坡整體下滑約1 000 m，此時滑體整體下滑速率依舊很快。中后部滑體在前緣地形較陡處發生飛躍的現象。下部河谷處堆積體方量仍不大，河谷處堆積體明顯沿著河谷方向向兩邊擴展，堆積體主要由前緣及中部巖土體堆積而成，堆積高度約20～30 m左右(見圖11,12)。

圖11 10000時步位移云圖(單位：m)

圖12 10000時步縱剖面示意(單位：m)

3.2.320000時步堆積分析

20 000時步滑體運動速度開始明顯減小，逐步趨于穩定狀態；滑坡整體下滑約1 500～1 600 m左右。滑體主要集中在前緣及河谷處，在河谷處形成較為穩定的堆積形態，堆積體高度約為60 m。堆積體主要由滑坡前緣及中部的滑體形成，有少量后緣巖土體。由于受到深沖溝的限制，堆積體主要集中于A區、B區正下方，沿著河谷的擴展寬度并不是很大(見圖13,14)。

圖14 20000時步縱剖面示意(單位：m)

3.3 堆積形態分析

滑坡失穩經過80 000時步后達到穩定狀態，在河谷處形成穩定的堰塞壩。巖土體整體滑移1 600 m左右。分別沿河流流向截取一條橫剖面線，沿主滑方向截取一條縱剖面線(圖15)，得到堰塞壩的橫剖面形態圖(圖16)和縱剖面形態圖(圖17)。

圖15 堰塞壩堆積形態

從堰塞壩橫剖面形態圖(圖16)中看出:堰塞壩形態呈中部高兩邊低的梯形形態分布，壩頂寬437.91 m，壩底寬994.39 m，壩高92.65 m，下游坡度為18.6°，上游坡度為18.7°。

圖16 堰塞壩橫剖面形態(單位：m)

從堰塞壩縱剖面形態圖(見圖17)中看出:堰塞壩在滑坡對岸一側堆積高度較高，在靠近滑坡一側堆積高度較低，對岸壩高100.80 m，本岸壩高92.65 m。

圖17 堰塞壩縱剖面形態(單位:m)

4 潰壩影響分析

基于模擬分析得到堰塞壩壩高為92.65 m，堰塞壩壩頂最低高程約為2 483 m，當上游持續來流導致水位漫頂時，堰塞壩前最大水深約為93 m，回水長度約為33 km。

采用DAMBRK潰壩洪水計算數學模型[29]，進行楞古水電站施工期和運行期的唐古棟滑坡A區失穩堆積堰塞壩潰壩影響預測分析。潰決起始水位取為堆積體壩頂高程2 483 m，潰口尺寸參照1967年唐古棟垮山滑坡潰壩，考慮為局部潰決，潰口深度按相同比例推算為55.8 m，潰口頂寬155 m，潰口底寬50 m，潰口邊坡坡度取1967年潰壩自然發展的左岸邊坡坡度，坡比為1∶0.94。計算結果見表5。

(1)施工期潰壩洪水計算結果表明：在各種流量方案下，強變形A區失穩堆積堰塞壩潰壩洪峰流量均遠大于楞古水電站擬選壩址校核洪水流量；上游不同來流方案對潰口最大洪峰的影響較小，各方案差值均在3 700 m3/s以內；當遭遇5a一遇洪水流量時，潰口洪峰流量為25 200 m3/s，約為楞古水電站擬選壩址校核洪水流量8 690 m3/s的2.9倍。

表5 A區失穩堆積堰塞壩潰壩計算方案及參數

(2)運行期潰壩洪水計算結果表明：在各種流量方案下，強變形A區失穩堆積堰塞壩潰壩洪峰流量均遠大于楞古水電站擬選壩址校核洪水流量；上游不同來流方案對潰口最大洪峰的影響較小，方案間差值均在3 700 m3/s以內；當遭遇5a一遇洪水流量時潰口洪峰流量為13 000 m3/s，約為楞古水電站擬選壩址校核洪水流量8 690 m3/s的1.5倍。

5 結 論

本文運用三維離散元軟件PFC3D，以唐古棟滑坡為研究對象，通過反演模擬1967年滑坡失穩堰塞壩堆積形態，驗證標定了離散元細觀參數。并基于此細觀參數進行強變形A區失穩堰塞壩的預測模擬。在此基礎上，開展堰塞壩潰壩影響預測分析，主要研究結論如下。

(1)基于滑坡室內物理力學試驗和抗剪強度反演得到的唐古棟滑坡體綜合物理力學參數，運用PFC3D軟件進行虛擬三軸測試，標定了一組細觀力學參數。

(2)通過對比分析1967年失穩堆積堰塞壩與PFC3D反演堰塞壩的形態，校核標定的細觀參數是合理的。堰塞壩底寬和本岸壩高偏差較小，頂寬、最大壩高、對岸壩高的偏差20%～40%，是由于1967年堆積堰塞壩的形態是通過調查殘余痕跡獲得，且壩頂寬殘余形態保留較差。

(3)強變形A區失穩堆積過程如下：A區整體下滑約1 600 m；滑動前期，巖土體移動速率較大，前緣滑體滑動速率較后緣滑體慢；滑動后期，堰塞壩主要由滑坡前緣及中部的滑體形成，有少量后緣巖土體；受沖溝限制，堰塞壩沿河谷的擴展寬度并不大。

(4)強變形A區失穩堆積堰塞壩的形態，順河向呈中部高兩邊低的梯形形態分布，壩頂寬437.91 m，壩底寬994.39 m，壩高92.65 m，下游坡度為18.6°，上游坡度為18.7°；橫河向對岸一側堆積高度較高，對岸壩高100.8 m，本岸壩高92.65 m。

(5)采用DAMBRK潰壩洪水計算數學模型，選取年平均流量、汛期平均流量和5 a一遇洪水流量，進行楞古水電站施工期和運行期的潰壩影響分析。在各種流量方案下，堰塞壩潰壩洪峰流量均遠大于楞古水電站擬選壩址校核洪水流量，上游不同來流方案對潰口最大洪峰的影響較小；強變形A區失穩導致的堰塞是影響楞古水電站壩址選擇的重要因素。

(6)通過模擬發現，PFC3D軟件用于模擬滑坡堰塞壩堆積過程及堆積形態有較好的適用性，尤其可以獲取堰塞壩堆積的三維形態，且堰塞壩高度較一般經驗公式更為準確，可有效運用于潰壩影響分析。