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

(1.西安理工大學(xué) 省部共建西北旱區(qū)生態(tài)水利國家重點(diǎn)實(shí)驗(yàn)室，陜西 西安 710048; 2.中國水利水電科學(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堰塞湖險情基本排除。通過對無人機(jī)拍攝的三維點(diǎn)云圖進(jìn)行測量(見圖1)，得到堰塞體右岸被沖出的梯形水槽頂寬約160 m，底寬約70 m，深約55 m。2018年11月3日17：40，原滑坡地點(diǎn)發(fā)生二次滑坡，由于二次滑坡的土石方堵塞了初次滑坡后自然泄流沖出的水槽，同時在初次滑坡的土石方上再堆積，致使堰塞壩高度升高，危險增加。2018年11月12日10：50，人工開挖的泄流槽開始過流，2018年11月15日08：00，堰塞湖水位基本穩(wěn)定，險情緩解。

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

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

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

1 模型簡介

1.1 潰口演變模型

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

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

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為考慮尾水淹沒的系數(shù)；B為潰口寬度，m；H為水庫水位，m；Z為潰口底高程；h為潰口水深,m。

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

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

1.2.1控制方程

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

(4)

(5)

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

1.2.2數(shù)值方法

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

2 模型構(gòu)建

2.1 研究區(qū)域概況

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

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

2.2 潰口流量計算

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

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

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

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

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

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

3 模擬結(jié)果分析

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

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

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

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

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

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

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

4 結(jié) 論

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

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

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