侯 精 明，馬 利 平，陳 祖 煜2，齊 文 超，王 琳

(1.西安理工大學(xué) 省部共建西北旱區(qū)生態(tài)水利國(guó)家重點(diǎn)實(shí)驗(yàn)室，陜西 西安 710048; 2.中國(guó)水利水電科學(xué)研究院，北京 100038)

堰塞湖是由滑坡、地震、泥石流等因素導(dǎo)致河道阻塞而形成的湖泊[1]。2018年10月10日22：06，西藏自治區(qū)昌都市江達(dá)縣和四川省甘孜藏族自治州白玉縣交界處發(fā)生山體滑坡，堵塞金沙江干流河道，形成堰塞湖。10月12日17：20，堰塞湖開始漫頂溢流，10月13日9:00，湖水水位下降20多米，潰壩隱患完全消除，10月13日20：00堰塞湖險(xiǎn)情基本排除。通過(guò)對(duì)無(wú)人機(jī)拍攝的三維點(diǎn)云圖進(jìn)行測(cè)量(見(jiàn)圖1)，得到堰塞體右岸被沖出的梯形水槽頂寬約160 m，底寬約70 m，深約55 m。2018年11月3日17：40，原滑坡地點(diǎn)發(fā)生二次滑坡，由于二次滑坡的土石方堵塞了初次滑坡后自然泄流沖出的水槽，同時(shí)在初次滑坡的土石方上再堆積，致使堰塞壩高度升高，危險(xiǎn)增加。2018年11月12日10：50，人工開挖的泄流槽開始過(guò)流，2018年11月15日08：00，堰塞湖水位基本穩(wěn)定，險(xiǎn)情緩解。

堰塞湖形成后，由于阻塞河道的堆積物顆粒松散，且堰塞體沒(méi)有泄洪設(shè)施，在上游不斷來(lái)流、湖水位持續(xù)升高的條件下，堰塞體潰壩風(fēng)險(xiǎn)逐漸升高[2]。為減小堰塞體潰決后對(duì)下游造成的災(zāi)害損失，準(zhǔn)確計(jì)算潰壩流量及下游洪水演進(jìn)過(guò)程就顯得尤為重要。目前，比較經(jīng)典且使用較多的潰口流量計(jì)算模型有DAMBRK[3]、BEED[4]、BREACH[5]模型等。張大偉[6]等基于物理試驗(yàn)研究了兩組粒徑差別明顯的砂樣形成的堰塞壩潰壩過(guò)程，并建立了具有物理意義的概念性潰口流量計(jì)算模型；在洪水演進(jìn)方面，楊志[7]等采用側(cè)向連接方式，建立了一維河道二維淹沒(méi)區(qū)耦合模型，模擬研究了黑河金盆水庫(kù)潰壩洪水演進(jìn)過(guò)程，但由于潰壩水流具有明顯的二維特性，且流態(tài)復(fù)雜，一維河道模型并不能準(zhǔn)確描述該過(guò)程；肖瀟[8]等基于非結(jié)構(gòu)網(wǎng)格，采用有限體積法建立了蓄滯洪區(qū)洪水演進(jìn)數(shù)學(xué)模型；王曉玲[9]等針對(duì)潰壩洪水在復(fù)雜淹沒(méi)區(qū)域中的演進(jìn)過(guò)程，建立了耦合VOF法的三維紊流數(shù)學(xué)模型，但該模型計(jì)算效率比較低，模擬8h的洪水演進(jìn)過(guò)程CPU耗時(shí)480 h。

為高效高精度地模擬潰壩洪水演進(jìn)過(guò)程，本文采用潰口演變模型DB-IWHR計(jì)算潰口流量過(guò)程，并采用基于GPU加速技術(shù)的二維水動(dòng)力模型模擬金沙江白格堰塞湖潰壩洪水演進(jìn)過(guò)程。通過(guò)與下游葉巴灘及蘇洼龍電站斷面實(shí)測(cè)流量過(guò)程進(jìn)行對(duì)比，證明了該模型可用來(lái)模擬計(jì)算堰塞湖潰壩洪水演進(jìn)過(guò)程。

圖1 初次滑坡后堰塞體三維點(diǎn)云圖Fig.1 Three-dimensional point cloud map after initial landslide

1 模型簡(jiǎn)介

1.1 潰口演變模型

采用中國(guó)水科院陳祖煜等于2014年提出的潰口演變模型DB-IWHR[10]計(jì)算潰口流量。該模型主要包括潰口流量計(jì)算和潰口擴(kuò)展兩部分。

潰口流量采用寬頂堰公式進(jìn)行計(jì)算，公式如下：

Q=C1C2B(H-z)3/2

(1)

C2=1.0-27.8(m-0.67)3

(2)

(3)

式中，Q為潰口流量，m3/s；C1為流量系數(shù)，建議取值范圍為1.43～1.69 m1/2/s[11]；C2為考慮尾水淹沒(méi)的系數(shù)；B為潰口寬度，m；H為水庫(kù)水位，m；Z為潰口底高程；h為潰口水深,m。

潰口擴(kuò)展過(guò)程為巖土工程中已經(jīng)被廣泛接受的簡(jiǎn)化的Bishop法。通過(guò)搜索可能的滑裂面，從而得到下切深度[12]。

1.2 洪水演進(jìn)模型

1.2.1控制方程

本文應(yīng)用具有守恒格式的平面二維淺水方程(簡(jiǎn)稱SWEs)來(lái)模擬潰壩洪水的演進(jìn)過(guò)程，忽略了運(yùn)動(dòng)黏性項(xiàng)、紊流黏性項(xiàng)、風(fēng)應(yīng)力和科氏力。二維非線性淺水方程的守恒格式可用如下的矢量形式[13]來(lái)表示：

(4)

(5)

式中，q為變量矢量，包括水深h，兩個(gè)方向的單寬流量qx和qy；g為重力加速度，u和v分別為x、y方向的流速；F和G分別為x、y方向的通量矢量；S為源項(xiàng)矢量；zb為河床底面高程，Cf=gn2/h1/3為謝才系數(shù)，n為曼寧系數(shù)。

1.2.2數(shù)值方法

本模型采用Godunov格式的有限體積法求解二維淺水方程。單元界面通量應(yīng)用HLLC格式的近似黎曼求解器求解[14]。通過(guò)靜水重構(gòu)來(lái)修正干濕邊界處負(fù)水深。底坡源項(xiàng)使用模型作者提出的底坡通量法處理。摩阻源項(xiàng)的計(jì)算使用半隱式法來(lái)提高穩(wěn)定性[15]。采用二階顯式Runge Kutta[16]方法進(jìn)行時(shí)間步長(zhǎng)的推進(jìn)。構(gòu)造具有二階時(shí)空精度的MUSCL型格式，有效解決復(fù)雜地形干濕界面處的負(fù)水深和偽高流速等現(xiàn)實(shí)中不存在的物理現(xiàn)象所造成的計(jì)算失穩(wěn)和物質(zhì)動(dòng)量的不守恒問(wèn)題。因潰壩洪水演進(jìn)過(guò)程模擬一般尺度較大，為提高計(jì)算效率，模型采用GPU并行計(jì)算技術(shù)實(shí)現(xiàn)高速運(yùn)算[17-18]。

2 模型構(gòu)建

2.1 研究區(qū)域概況

本研究區(qū)域概況已于前節(jié)描述。本文模擬研究的范圍為從白格堰塞體至下游235 km長(zhǎng)的河段(見(jiàn)圖2)。

圖2 研究區(qū)域位置Fig.2 Regional map

2.2 潰口流量計(jì)算

采用潰口演變模型DB-IWHR對(duì)“10·10”與“11·3”兩次堰塞湖潰壩事故進(jìn)行模擬計(jì)算，具體參數(shù)見(jiàn)表1，模擬計(jì)算結(jié)果見(jiàn)圖3。

2.3 水動(dòng)力模型設(shè)置

模型的輸入資料為地形數(shù)據(jù)、入流數(shù)據(jù)等。地形資料為從白格堰塞體到蘇洼龍段的數(shù)字地形高程(DEM)數(shù)據(jù)，由于該區(qū)域精細(xì)地形難以獲取，故基于美國(guó)對(duì)地觀測(cè)衛(wèi)星Terra提供的30 m分辨率DEM地形數(shù)據(jù)(見(jiàn)圖4)，利用遙感影像提取河道寬度，并下挖一定深度概化出水下河道地形。研究區(qū)域高程落差大約500 m，河段長(zhǎng)235 km，平均坡度0.002，平均寬度300～600 m。模擬網(wǎng)格462萬(wàn)個(gè)，其中葉巴灘段87萬(wàn)個(gè)，蘇洼龍段375萬(wàn)個(gè)。入流數(shù)據(jù)為DB-IWHR潰壩模型計(jì)算的兩次潰壩事故的潰口流量過(guò)程。上游邊界為入流邊界，下游邊界為自由出流的開邊界，其余邊界定義為封閉邊界。由于金沙江斷流后下游河道仍有大約200 m3/s的流量，故給予河道200 m3/s的恒定流量作為初始條件。根據(jù)現(xiàn)場(chǎng)實(shí)際情況選取曼寧系數(shù)為0.02[19]，庫(kù)朗數(shù)為0.50，共模擬40 h洪水演進(jìn)過(guò)程。

圖4 研究區(qū)域數(shù)字地形高程Fig.4 DEM of the study area

表1 金沙江白格堰塞湖潰口流量計(jì)算參數(shù)Tab.1 Calculation parameters of dam break discharge at Baige barrier lake in Jinsha River

圖3 所選模擬潰壩洪水過(guò)程線Fig.3 Selected simulated inflow hydrograph of the dam break

3 模擬結(jié)果分析

3.1 “10·10”白格堰塞湖結(jié)果分析

針對(duì)“10·10”白格堰塞湖潰壩事故，應(yīng)用二維水動(dòng)力模型對(duì)堰塞體下游235 km長(zhǎng)的河段進(jìn)行洪水演進(jìn)模擬。由于下游葉巴灘與蘇洼龍?zhí)幎加性诮ɑ蛞呀ǔ傻乃娬?，各水電站的防洪安全是人們最為關(guān)心的，故選取這兩處的流量過(guò)程進(jìn)行模擬對(duì)比。其中葉巴灘距堰塞體54 km，蘇洼龍距堰塞體224 km。模擬耗時(shí)61 min，模擬結(jié)果見(jiàn)圖5。從圖中可看出，葉巴灘處模擬的流量過(guò)程在上升和下降階段都比實(shí)測(cè)快，洪峰流量比實(shí)測(cè)大5 247 m3/s，洪峰出現(xiàn)時(shí)間基本一致，誤差主要在于潰口演變模型DB-IWHR模擬計(jì)算的潰壩流量峰值比實(shí)測(cè)要大；蘇洼龍?zhí)幠M的流量過(guò)程與實(shí)測(cè)流量過(guò)程吻合較好，洪峰流量比實(shí)測(cè)大220 m3/s，在流量過(guò)程下降階段比實(shí)測(cè)要快。

圖5 “10·10”潰壩洪水演進(jìn)模擬與實(shí)測(cè)流量過(guò)程Fig.5 Simulated and measured discharges of the "10·10" dam break flood propagation

3.2 “11·3”白格堰塞湖結(jié)果分析

二維水動(dòng)力模型對(duì)“11·3”白格堰塞湖潰壩洪水演進(jìn)模擬的結(jié)果見(jiàn)圖6，模擬耗時(shí)74 min。

圖6 “11·3”潰壩洪水演進(jìn)模擬與實(shí)測(cè)流量過(guò)程Fig.6 Simulated and measured discharge of the "11·3" dam break flood propagation

模擬研究結(jié)果表明，在葉巴灘處，模擬的洪峰流量比實(shí)測(cè)小1 012 m3/s，流量過(guò)程整體延遲2 h；在蘇洼龍?zhí)幠M流量過(guò)程與實(shí)測(cè)吻合較好，洪峰流量比實(shí)測(cè)小614 m3/s，在下降階段比實(shí)測(cè)稍大。在葉巴灘處產(chǎn)生誤差的主要原因在于地形數(shù)據(jù)的精度不高，難以準(zhǔn)確描述主要河道地形，尤其是水下地形。

4 結(jié) 論

應(yīng)用二維水動(dòng)力模型對(duì)“10·10”與“11·3”金沙江白格堰塞湖潰壩洪水演進(jìn)進(jìn)行模擬，通過(guò)將模擬流量過(guò)程與下游葉巴灘、蘇洼龍電站斷面處的實(shí)測(cè)流量過(guò)程進(jìn)行對(duì)比，可得到如下結(jié)論。

(1) 對(duì)于無(wú)高精度地形資料的山區(qū)，該二維水動(dòng)力模型可以較好地模擬潰壩洪水的演進(jìn)過(guò)程。在對(duì)“10·10”與“11·3”兩次潰壩洪水的模擬中，蘇洼龍?zhí)幠M結(jié)果較好，而葉巴灘處誤差較大，主要原因是該處所用的terra 地形和概化的河道地形數(shù)據(jù)精度不高，較實(shí)際有較大偏差。

(2) 在堰塞體至蘇洼龍共462萬(wàn)網(wǎng)格的地形數(shù)據(jù)上進(jìn)行洪水演進(jìn)模擬，共模擬40 h洪水演進(jìn)過(guò)程，由于該二維水動(dòng)力模型基于GPU并行計(jì)算技術(shù)進(jìn)行加速，采用英偉達(dá)RTX 2080 GPU，兩次潰壩洪水演進(jìn)模擬耗時(shí)分別為61 min和74 min，可見(jiàn)該二維水動(dòng)力模型在模擬洪水演進(jìn)時(shí)是非常高效的，對(duì)于爭(zhēng)分奪秒的洪水應(yīng)急搶險(xiǎn)事件，該模型可實(shí)現(xiàn)快速預(yù)測(cè)，為決策者提供有力數(shù)據(jù)支撐。