某堰塞壩沖刷潰決數值模擬分析

呂 霞，楊興國，2，范 剛，2，林子鈺，劉大瑞

(1.四川大學水利水電學院，四川 成都 610065；2.四川大學水力學與山區河流開發保護國家重點實驗室，四川 成都 610065)

0 引言

堰塞壩(體)是在一定的地形地貌條件下，由滑坡、泥石流、熔巖流等堵塞山谷和河道形成的天然壩體[1-3]。堰塞壩是自然作用的產物，具有壩體形態不規則、物質組成隨機性強、內部結構及空間分布復雜等特點，受上游來流、滲流侵蝕等外載荷的影響，容易在短期內發生失穩破壞，對下游社會和自然環境造成重大危害[4-7]。1974年秘魯Mayunmarca村發生滑坡-碎屑流災害，導致450人死亡，隨后碎屑物質堵塞河流形成堰塞壩，潰決后造成下游巨大災難[8]。2008年汶川8.0級大地震誘發多處堰塞壩，其中唐家山是規模最大、危險性最高的堰塞湖，嚴重威脅綿陽市城區及下游130萬人的生命安全[9]。2018年8月6日，漢源縣富泉鎮西溝發生150萬m3的大規?；?，并在主溝內形成體積約40萬m3的堰塞壩[10]。2018年10月11日及11月3日，金沙江白格滑坡兩次堵塞金沙江，形成的堰塞湖淹沒了上游村鎮，堰塞壩潰決后洪水沖擊下游造成了巨大損失[11]。2020年6月17日，四川省丹巴縣梅龍溝發生特大泥石流災害，泥石流局部堵塞小金川河河道，形成堰塞壩。壩體潰決后，導致電站大壩和大量房屋被淹，直接經濟損失約2000萬元[12]。

美國地質調查局研究發現，漫頂導致的潰壩案例高達90%，而滲流破壞導致的潰壩數量僅占10%[13]。因此有必要深入研究堰塞壩漫頂潰決過程以及潰決機理，為潰壩洪水風險分析和應急預案的編制提供理論和技術支持，以期削弱潰決洪水對下游的破壞程度。目前，已有眾多學者對堰塞壩潰決過程和機理分析做了大量研究。趙天龍等通過不同來流條件下粗粒級配土料的沖刷性能實驗，揭示壩體宏觀潰口的形成機制[14]。鄧明楓等通過壩體漫頂潰決模型實驗，研究了堰塞壩漫頂潰決的動力機制[15]。傅旭東等采用Osman和Thorne模型發展了堰塞壩潰決過程模擬模型，并將其應用于唐家山堰塞湖潰決過程的模擬，再現了潰決洪水流量過程線、潰口展寬和下切過程[16]。Zhang等在天然河道中開展了大型物理模型，提出了堰塞壩縱向以及橫向潰決演化模型，并利用白格滑坡的兩次潰壩事件，驗證了該演化模型的準確性[17]。劉若星等將NWS BREACH模型應用于唐家山堰塞湖泄流過程的模擬，驗證了改模型的合理性[18]。楊興國等構建了考慮水流侵蝕與潰口間歇性崩塌的堰塞壩潰決演化模型，初步揭示了冰磧土滑坡-泥石流運移與堵江機制[19]。目前，針對堰塞壩潰決機理和潰決過程的研究大多未使用堰塞壩的實地模型，而對壩體做了相應的概化，其研究成果可為堰塞壩潰壩機制研究提供重要借鑒，基于概化模型的研究成果需要實際堰塞壩模擬結果的檢驗與驗證。

2020年8月31日四川省甘洛縣黑西洛溝發生滑坡-泥石流災害，溝內沖出大量固體物源堵塞尼日河并形成堰塞壩。本論文通過等高線地形圖，建立了黑西洛堰塞壩實地模型，使用Flow-3D軟件對堰塞壩的潰決過程進行數值模擬，得到壩體泄流時的流速特征和潰口變化過程，以進一步了解潰口的發展規律，為堰塞壩風險分析和應急預案的編制提供技術支撐。

1 黑西洛滑坡-堰塞湖災情概況

2020年8月31日，四川省甘洛縣阿茲覺鄉發生山體滑坡-泥石流-堰塞湖災害鏈，尼日河被大量泥石流物質堵塞，形成堰塞壩，嚴重威脅當地人民生命財產安全。黑西洛溝道兩側松散物源分布廣泛，植被不發育。災害發生時，溝內的松散物質被山洪裹挾帶走，并在搬運過程中不斷鏟刮底層物源，從而演變為滑坡和崩塌災害。在溝內通道持續下切的同時，兩側岸坡持續垮塌，坍塌體進入溝道后，滑坡規模急速擴大，最終演變為泥石流災害。大量的泥石流碎屑以近乎垂直的方式沖入尼日河，阻塞尼日河，形成堰塞壩。經現場實測，堰塞壩壩體縱向方向長度約200m，順河向長度約為400m，高度約30m，堰塞壩體積約100萬m3。堰塞壩形成約20min之后，尼日河上游水位迅速上漲，達到了772.0m，上游河水漫過堰塞壩，并在壩體中偏左側低洼地帶形成天然的泄流槽，由于壩體材料松軟，堰塞壩在10min左右即潰決完畢。堰塞壩潰決后，該段河道截彎取直，近乎順流，如圖1所示。壩體潰決后沖毀下游場鎮、學校、房屋和公路，造成阿茲覺鄉受災人口1730人，3人失聯。溝口的成昆鐵路大橋被沖毀，壩體下游約1.2km的道路以及多處橋梁被沖毀，經濟損失嚴重[20]。

圖1 2020年黑西洛溝滑坡-泥石流-堰塞湖災害情況

通過現場調查和資料收集，本次滑坡-泥石流-堰塞湖災害鏈是典型的“小水大災”災害，災害發生時大量松散堆積物被沖出溝道，堵塞尼日河，形成具有松散結構的堰塞壩，并且堰塞壩在短時間內潰決，壩體上游庫區內的洪水快速下泄，造成下游嚴重的生命財產損失。本文采用現場實測地形，使用數值模擬方法重現和分析黑西洛堰塞壩沖刷潰決過程，以期為今后類似堰塞壩的應急處置和編制避險預案提供技術支持。

2 數值計算模型

本文以黑西洛溝堰塞湖為研究對象，利用Flow-3D數值方法，開展堰塞湖潰決過程數值模擬研究。Flow-3D是基于CFD解算技術的仿真模擬軟件，其計算原理是有限差分法，采用質量連續性方程和Navier-Stokes方程作為控制方程。其FAVOR技術既能進行網格內的流場計算，又能進行固態流場的數值模擬，特別是對自由表面流的模擬具有顯著優勢。另外，該技術還可以通過簡單的矩形網格來構建任何復雜的形狀，并通過VOF對自由面進行追蹤，從而更接近于真實的流體運動，并能夠對流體的多種特性進行模擬。堰塞壩潰決后，洪水迅速沖出，對壩體物質產生強烈的沖蝕作用，考慮到下游彎曲河道處存在紊流現象和劇烈變形等特點，采用RNG湍流模型對潰壩洪水進行模擬，該模型能很好地模擬潰決過程中洪水的復雜流態。此外，將壩體設置成泥沙沖刷模式，以模擬其在水流沖擊作用下逐漸潰決的過程。Flow-3D軟件中的泥沙沖刷模型，通過對泥沙起動、堆積、推移質輸移等計算來描述泥沙運動，與壩體沖刷破壞機理一致。

本文根據獲取的等高線地形圖，建立了黑西洛堰塞壩的三維實體模型，并將其導入Flow-3D軟件中。模型的計算范圍長900m，寬約460m，高約125m，如圖2所示。本文使用正交的矩形網格對模型的計算區域進行劃分，網格大小為2.5m×2.5m×2.5m，計算網格總數為331萬個。進一步對計算區域進行邊界條件的設定，將模型入口邊界設為壓力邊界，出口邊界設為自由出流邊界，底面和兩側邊界設為Wall(墻)邊界，自由液面上方設一個標準大氣壓。并將壩體上游初始水位條件設為772.0m，為上游庫區水位達到的最大值。堰塞體物質來源于溝內松散堆積體，顆粒級配曲線如圖3所示[20]，將6個監測點顆粒中值粒徑的平均值作為本次模擬計算的典型壩體粒徑，即d50=11.75mm。本次設定堰塞壩堆積區域泥沙區域，泥沙密度為2000kg/m3，臨界Shields系數為0.32，挾帶系數設為0.013。根據現場調查和模型試算，壩體在600s后基本穩定，故設定求解時間為600s，最小步長為1×10-7s。

圖2 堰塞壩的數值模擬模型

圖3 堰塞體顆粒級配測點位置及級配曲線

3 模擬結果分析

3.1 流速

圖4是堰塞壩潰決過程中水流流速的分布情況。堰塞壩形成泄流通道后流速較高，流道中后段流速最大為17.5m/s。因流道內的水流流速較大，水流的侵蝕能力較強，潰口隨時間不斷擴展。壩體潰決持續76s后，堰塞體形成貫通的流道，流道內流速高達17.5m/s；潰決持續203s后，流速在流道內分布不均勻，流速相對較高，水流的侵蝕作用強烈，流道不斷下切變寬；潰決持續時間417s后，壩體泄流過程中，上游水位在逐漸降低，導致泄流槽內流速開始緩慢下降，此時流速最大達到15.3m/s。如圖4所示，泄流過程中流速分布有一定相似性，水流進入泄流槽之前，流速基本為0，上游庫區處于靜壓平衡狀態；水流流入泄流槽后，進口處水流流速明顯上升，整個潰決過程中，泄流槽中后段水流流速一直較快。泄流槽末端，槽內斜坡道上的水流流速急速上升，水流侵蝕作用增強，不斷侵蝕壩體，導致泄流槽不斷擴展，且末端斜坡道上的跌坎不斷向上游推進，泄流后期，跌坎消失。

圖4 堰塞壩沖刷潰決過程水面流速動態演進模擬結果

3.2 潰口擴展過程

3.2.1沖淤形態

水流運動時會對泥沙顆粒產生上舉力和拖曳力，當水流強度達到泥沙起動條件時，能使泥沙顆粒躍起，并被水流帶走，壩體產生沖刷現象，導致泄流槽底部逐漸下切。圖5是壩體潰決過程中泥沙的沖淤演進變化過程，其中正值表示淤積，負值表示沖刷。如圖5所示，潰決持續發展過程中，壩體的淤積高度和沖刷深度不斷增加。持續沖刷76s后，泄流槽整個連通并不斷擴展，水流流速高達17.5m/s，水流的挾沙能力強，沖刷速度快，此時堰塞體的最大沖刷深度達22.7m；壩體潰決持續203s后，水流不斷帶動泥沙向下運動，堰塞壩的沖刷深度持續增大，泄流槽底部下切和兩側侵蝕明顯，此時堰塞體的沖刷深度達25.2m；潰決持續417s后，上游水位降低，泄流槽內流量變小，流道中水面寬度和水深逐漸降低，泄流槽底部下切和側壁擴展速度變慢，壩體的沖刷深度增長變緩，此時壩體最大沖刷深度達25.7m。在泄流槽內，各位置的沖蝕狀況也有差別，泄流槽進口處以及壩體中上游區域的沖刷深度變化緩慢，是由于該位置的水流流速較小，而泄流槽下游段流速變高，水流攜帶大量泥沙運動，導致下游壩體的沖刷深度和潰口展寬都較上游面更大，泄流槽明顯下切，并呈向下的“喇叭口”形狀。

圖5 堰塞壩沖刷潰決過程泥沙沖淤動態演進模擬結果

堰塞壩下游，受重力以及床面摩擦等因素的影響，被水流帶走的泥沙在此處沉積，出現淤積現象。泄流初始階段，水流攜帶著流道內的泥沙向下游運輸，由于水流流速不大而堆積在近壩體河段。如圖5所示，堰塞壩持續泄流76s后，河道內出現淤積現象，最大淤積高度達到4.5m；隨著水流持續下泄以及潰口不斷發展，靠近壩體的河段也被水流沖刷，初期沉積在此處的泥沙被沖向下游。泄流后期，流道內流速基本達到穩定，流道內沖刷深度變化不明顯，泥沙不斷淤積在壩體下游，堆積高度最大值為8.5m。泥沙在壩體下游左岸堆積，是因為黑西洛堰塞壩處在彎曲河段，潰決洪水產生橫向環流對岸坡形成頂沖之勢，使得凹岸沖刷而凸岸淤積。水流流經彎曲河道，由于離心作用，導致下游左岸水深變淺，流速降低，水流的搬運能力下降，泥沙在河道左側逐漸堆積。

3.2.2潰口斷面歷時變化

對黑西洛堰塞壩潰決的橫、縱剖面進行截取，進而揭示堰塞壩漫頂潰決的破壞規律。截面具體位置如圖6所示，其中1-1′截面是縱剖面，坐標為x=200m，2-2′截面、3-3′截面、4-4′截面為泄流槽不同位置的橫斷面，坐標分別為y=250m、y=430m、y=600m。

圖6 截面示意圖

堰塞壩縱剖面(1-1′截面)不同時間的沖淤變化如圖7所示，色條表示堰塞壩的沖淤情況，其中正值代表淤積，負值代表沖刷。由圖7可知，在堰塞壩泄流槽斜坡道上，沖刷現象最顯著，泄流槽末端水流流速急速上升，沖蝕能力明顯增強。泄流槽末端不斷被水流侵蝕破壞，最大沖刷深度達到25.0m。漫頂水流侵蝕斜坡段，使其產生跌坎，水流流經跌坎，流速陡然增加，挾沙能力變強，致使跌坎處的潰決速度明顯更快。在跌坎的上游，水流較穩定，流速小，其沖刷強度較低；跌坎的下游，水流湍急，水流流速較大，水流對壩體產生強烈的沖蝕作用，使跌坎附近的泥沙不斷被沖走，泄流過程中跌坎不斷向上游推進，產生明顯的溯源沖刷。溯源侵蝕會在短期內對堰塞壩產生強烈的沖蝕作用，由圖7可知，潰決持續381s后，溯源侵蝕已經發展到壩體上部，跌坎消失。

圖7 1-1′截面沖淤演變過程

圖8—10分別展示了堰塞壩不同橫斷面(2-2′截面、3-3′截面、4-4′截面)在不同時間的沖淤演進過程。如圖8所示，該斷面侵蝕程度較弱，由于潰口上游段流速較小，水流的侵蝕能力不強，因此流道發展比較緩慢，壩體的沖刷深度最大為4.7m；圖9為3-3′橫斷面潰口的演變過程，該斷面流速增大，水流沖刷能力變強，水流攜帶更多的泥沙向下運動，持續沖刷417s之后，此截面最大沖刷深度達到了21.3m；如圖10所示，水流在泄流槽末端侵蝕作用最強烈，此斷面位置的流速明顯增大，沖刷能力強，最大沖刷深度達25.7m，潰口擴展最為迅速，且泄流槽寬度不斷拓寬。

圖8 2-2′截面沖淤演變過程

圖9 3-3′截面沖淤演變過程

圖10 4-4′截面沖淤演變過程

整個泄流過程中潰口的發展過程如下，泄流槽下游端沖刷最強烈，潰口向底部下切和兩側侵蝕，之后逐漸向上游發生侵蝕，當侵蝕擴展到壩體上游后，又開始向壩體兩側緩慢侵蝕。水流挾帶泥沙向下游運動不斷侵蝕坡腳，側壁發生破壞，導致泄流槽的寬度逐漸加寬?？偟恼f來初始潰口發展較為緩慢，然后逐漸加快，最后流道內發展達到平穩狀態。

4 結論

本文采用Flow-3D軟件對黑西洛溝堰塞湖自然泄流潰決過程進行模擬和分析，得出以下結論：

(1)堰塞壩漫頂沖刷是一個復雜過程，壩體潰決過程中流速變化大且分布不均勻，在泄流槽的下游段流速明顯增大，水流的沖刷能力增強。

(2)泄流過程中，泄流槽末端斜坡道上流速最大，壩體的侵蝕破壞最明顯，流道下切深度達25.7m，并在堰塞壩下游出現淤積現象，最大淤積高度達8.5m。

(3)堰塞壩泄流過程中出現了溯源侵蝕現象，水流不斷使潰口加寬、加深，泄流跌坎不斷向壩體上游移動。

(4)黑西洛堰塞壩自然泄流過程中，初期潰口發展緩慢，水流對壩體的沖刷破壞主要發生在泄流槽末端，隨后向其兩側和上游侵蝕。