某堰塞壩沖刷潰決數(shù)值模擬分析

呂 霞，楊興國(guó)，2，范 剛，2，林子鈺，劉大瑞

(1.四川大學(xué)水利水電學(xué)院，四川 成都 610065；2.四川大學(xué)水力學(xué)與山區(qū)河流開發(fā)保護(hù)國(guó)家重點(diǎn)實(shí)驗(yàn)室，四川 成都 610065)

0 引言

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

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

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

1 黑西洛滑坡-堰塞湖災(zāi)情概況

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

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

通過現(xiàn)場(chǎng)調(diào)查和資料收集，本次滑坡-泥石流-堰塞湖災(zāi)害鏈?zhǔn)堑湫偷摹靶∷鬄?zāi)”災(zāi)害，災(zāi)害發(fā)生時(shí)大量松散堆積物被沖出溝道，堵塞尼日河，形成具有松散結(jié)構(gòu)的堰塞壩，并且堰塞壩在短時(shí)間內(nèi)潰決，壩體上游庫(kù)區(qū)內(nèi)的洪水快速下泄，造成下游嚴(yán)重的生命財(cái)產(chǎn)損失。本文采用現(xiàn)場(chǎng)實(shí)測(cè)地形，使用數(shù)值模擬方法重現(xiàn)和分析黑西洛堰塞壩沖刷潰決過程，以期為今后類似堰塞壩的應(yīng)急處置和編制避險(xiǎn)預(yù)案提供技術(shù)支持。

2 數(shù)值計(jì)算模型

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

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

圖2 堰塞壩的數(shù)值模擬模型

圖3 堰塞體顆粒級(jí)配測(cè)點(diǎn)位置及級(jí)配曲線

3 模擬結(jié)果分析

3.1 流速

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

圖4 堰塞壩沖刷潰決過程水面流速動(dòng)態(tài)演進(jìn)模擬結(jié)果

3.2 潰口擴(kuò)展過程

3.2.1沖淤形態(tài)

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

圖5 堰塞壩沖刷潰決過程泥沙沖淤動(dòng)態(tài)演進(jìn)模擬結(jié)果

堰塞壩下游，受重力以及床面摩擦等因素的影響，被水流帶走的泥沙在此處沉積，出現(xiàn)淤積現(xiàn)象。泄流初始階段，水流攜帶著流道內(nèi)的泥沙向下游運(yùn)輸，由于水流流速不大而堆積在近壩體河段。如圖5所示，堰塞壩持續(xù)泄流76s后，河道內(nèi)出現(xiàn)淤積現(xiàn)象，最大淤積高度達(dá)到4.5m；隨著水流持續(xù)下泄以及潰口不斷發(fā)展，靠近壩體的河段也被水流沖刷，初期沉積在此處的泥沙被沖向下游。泄流后期，流道內(nèi)流速基本達(dá)到穩(wěn)定，流道內(nèi)沖刷深度變化不明顯，泥沙不斷淤積在壩體下游，堆積高度最大值為8.5m。泥沙在壩體下游左岸堆積，是因?yàn)楹谖髀逖呷麎翁幵趶澢佣?，潰決洪水產(chǎn)生橫向環(huán)流對(duì)岸坡形成頂沖之勢(shì)，使得凹岸沖刷而凸岸淤積。水流流經(jīng)彎曲河道，由于離心作用，導(dǎo)致下游左岸水深變淺，流速降低，水流的搬運(yùn)能力下降，泥沙在河道左側(cè)逐漸堆積。

3.2.2潰口斷面歷時(shí)變化

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

圖6 截面示意圖

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

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

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

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

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

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

整個(gè)泄流過程中潰口的發(fā)展過程如下，泄流槽下游端沖刷最強(qiáng)烈，潰口向底部下切和兩側(cè)侵蝕，之后逐漸向上游發(fā)生侵蝕，當(dāng)侵蝕擴(kuò)展到壩體上游后，又開始向壩體兩側(cè)緩慢侵蝕。水流挾帶泥沙向下游運(yùn)動(dòng)不斷侵蝕坡腳，側(cè)壁發(fā)生破壞，導(dǎo)致泄流槽的寬度逐漸加寬。總的說來初始潰口發(fā)展較為緩慢，然后逐漸加快，最后流道內(nèi)發(fā)展達(dá)到平穩(wěn)狀態(tài)。

4 結(jié)論

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

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

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

(3)堰塞壩泄流過程中出現(xiàn)了溯源侵蝕現(xiàn)象，水流不斷使?jié)⒖诩訉?、加深，泄流跌坎不斷向壩體上游移動(dòng)。

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