徐云乾，袁明道，史永勝，張旭輝，潘展釗，馬妍博

(1.廣東省水利水電科學(xué)研究院，廣東 廣州 510635；2.廣東省大壩安全技術(shù)管理中心，廣東 廣州 510635；3.廣東省水動(dòng)力學(xué)應(yīng)用研究重點(diǎn)實(shí)驗(yàn)室，廣東 廣州 510635)

水庫(kù)大壩的防洪、發(fā)電、灌溉等功能為社會(huì)經(jīng)濟(jì)的發(fā)展提供了重要的基礎(chǔ)保障。水庫(kù)大壩失事將對(duì)下游洪泛區(qū)的人民生命財(cái)產(chǎn)安全及經(jīng)濟(jì)社會(huì)穩(wěn)定造成不可估量的影響。為了對(duì)水庫(kù)大壩潰壩風(fēng)險(xiǎn)進(jìn)行有效的防范，需要對(duì)潰壩洪水演進(jìn)進(jìn)行模擬，確定洪水淹沒(méi)范圍、淹沒(méi)水深、洪水流速、前鋒到達(dá)時(shí)間等洪水風(fēng)險(xiǎn)要素，并進(jìn)行避洪轉(zhuǎn)移分析，最大限度的減輕洪水風(fēng)險(xiǎn)影響。

目前主流通用的洪水?dāng)?shù)值模擬軟件主要有：DHI(丹麥水力研究所)的MIKE11/21，西班牙Next Limit 的RealFlow、荷蘭的Delft3D、美國(guó)陸軍工程兵團(tuán)開發(fā)的HEC-RAS(Hydrologic Engineering Center-River Analysis System)以及我國(guó)陳祖煜院士團(tuán)隊(duì)在2015年基于Excel平臺(tái)利用VB語(yǔ)言編寫了壩體潰口擴(kuò)展程序DB-IWHR。常配合使用的GIS軟件(Geographic Information System)主要有ArcGIS、QGIS等，兩者聯(lián)合用以預(yù)測(cè)洪水演進(jìn)風(fēng)險(xiǎn)和制定撤離路線。孫銳嬌[1]、陳建峰[2]、田景環(huán)[3]、賀娟[4]、周毅[5]、吳博[6]等利用HEC-RAS軟件模擬了大壩潰決后洪水在下游的洪水演進(jìn)過(guò)程，并結(jié)合了ArcGIS或HEC-RAS自帶的地理信息模塊分析研究下游區(qū)域的淹沒(méi)范圍。本文運(yùn)用HEC-RAS的二維模型聯(lián)合QGIS地理信息平臺(tái)模擬漫頂潰壩后洪水演進(jìn)，工程適用性較好，其結(jié)果對(duì)洪水災(zāi)害防御工作具有重要意義。

QGIS(原稱Quantum GIS)是跨平臺(tái)的開源版桌面地理信息系統(tǒng)，QGIS提供與其他開源GIS軟件包的集成，包括PostGIS、GRASS和MapServer，同時(shí)能集成當(dāng)前流行的Python插件擴(kuò)展，其提供的豐富人口、經(jīng)濟(jì)、行政區(qū)劃、交通設(shè)施的基礎(chǔ)地理信息系統(tǒng)，極大的方便了洪水風(fēng)險(xiǎn)圖的編制。

1 洪水分析方法

HEC-RAS的二維水動(dòng)力學(xué)模型，模型精度高，模型構(gòu)建要求簡(jiǎn)單，適用于平原地區(qū)泛濫以及山區(qū)河道洪水等多種場(chǎng)景。本研究主要利用其中的潰壩計(jì)算模塊和洪水演進(jìn)模塊。

在進(jìn)行水庫(kù)潰決模擬時(shí)，需要建立上游庫(kù)區(qū)、壩體以及下游淹沒(méi)區(qū)的聯(lián)合模型，并輸入水位庫(kù)容曲線、大壩潰決參數(shù)以及下游河道邊界條件等。

HEC-RAS二維水動(dòng)力學(xué)模型采用的是淺水方程，即Navier-Stokes方程的二維簡(jiǎn)化形式。

連續(xù)方程

(1)

動(dòng)量方程

(2)

式中，H為水面高程，m；h為水深，m；為流速，m/s；R為水力半徑，m；q為旁側(cè)入流單寬流量，m2/s；g為重力加速度，m/s2；υt為水平方向運(yùn)動(dòng)黏度，m2/s；cf為河床底部糙率；f為科里奧利系數(shù)；k為垂直方向單位矢量；n為糙率。

潰口出流假設(shè)為寬頂堰進(jìn)行計(jì)算，計(jì)算方程為

(3)

式中，Q為流量，m3/s；ε為堰流側(cè)收系數(shù)；m為流量系數(shù)；B為寬頂堰寬度，m；s為淹沒(méi)系數(shù)；Zi，Zi+1為堰前和堰后水位，m。

HEC-RAS二維水動(dòng)力學(xué)模型混合了有限體積法和有限差分法，邊界網(wǎng)格為不規(guī)則多邊形，非邊界網(wǎng)格為正方形。因其每個(gè)網(wǎng)格的邊類似下游河道斷面，所以在較低的網(wǎng)格密度或初始地形精度不高時(shí)，仍能提取地形細(xì)節(jié)，保證模型精度。

HEC-RAS中大壩為線性結(jié)構(gòu)，潰口預(yù)測(cè)模塊可以手動(dòng)輸入，也可以采取其集成的以往學(xué)者根據(jù)歷史潰壩數(shù)據(jù)推導(dǎo)出的潰口尺寸回歸方程輸入相關(guān)參數(shù)預(yù)測(cè)潰口。本次模擬采用能夠預(yù)測(cè)潰口的Von Thun and Gillette方程[7]估計(jì)潰口最大尺寸及其成形時(shí)間。Von Thun和Gillette收集了57座大壩潰壩數(shù)據(jù)，總結(jié)出潰口尺寸和壩前水深、庫(kù)容系數(shù)的關(guān)系，并建議使用1H∶1V的潰口坡比(粘土心墻壩則建議使用1H∶2V至1H∶3V的潰口坡比)，潰口形成時(shí)間根據(jù)大壩材料由兩組方程進(jìn)行計(jì)算。

潰口尺寸計(jì)算方程為

Bavg=2.5hw+Cb

(4)

式中，Bavg為潰口平均寬度，m；hw為潰口底起算水深，m；Cb為庫(kù)容系數(shù)，m。潰口幾何參數(shù)見(jiàn)圖1，庫(kù)容系數(shù)見(jiàn)表1。

圖1 潰口幾何參數(shù)

表1 庫(kù)容系數(shù)

以hw為參數(shù)的潰口形成時(shí)間計(jì)算方程為

(5)

以hw和Bavg為參數(shù)的潰口形成時(shí)間計(jì)算方程為

(6)

2 洪水分析模型構(gòu)建

2.1 工程概況及計(jì)算范圍確定

西山河發(fā)源于陽(yáng)春、信宜兩縣交界的西山山脈的三甲頂，海拔1 211 m，自西北向東南轉(zhuǎn)東北，流經(jīng)永寧、那離、大河、陂面至合水墟匯入漠陽(yáng)江干流。流域面積989 km2，河長(zhǎng)l08 km，河床比降2.03‰。大河水庫(kù)位于廣東省陽(yáng)江市境內(nèi)漠陽(yáng)江的支流西山河上，距陽(yáng)春市23 km，距陽(yáng)江市80 km。壩址以上控制的流域面積438 km2。水庫(kù)總庫(kù)容3.322億m3，水庫(kù)正常蓄水位110.00 m，死水位96.00 m，為Ⅱ等大(2)型水庫(kù)。

為保證模型區(qū)域能覆蓋可能淹沒(méi)區(qū)域，結(jié)合西山河以及漠陽(yáng)江流域的地形特點(diǎn)，經(jīng)過(guò)多次初步模擬，確定被模擬范圍如圖2、圖3所示。

圖2 大河水庫(kù)dem圖與行政區(qū)劃圖疊加

圖3 大河水庫(kù)水位庫(kù)容關(guān)系曲線

2.2模型參數(shù)率定與模型驗(yàn)證

河床及邊灘糙率是影響潰壩洪水向下游演進(jìn)傳播的主要參數(shù)，因此對(duì)其進(jìn)行驗(yàn)證十分必要，本文根據(jù)西山河典型洪水過(guò)程進(jìn)行率定，將典型洪水過(guò)程的流量數(shù)據(jù)(見(jiàn)圖4)結(jié)合泄流曲線(見(jiàn)圖5)輸入建立的計(jì)算模型，通過(guò)調(diào)整河床及邊灘各斷面糙率，通過(guò)河道水流傳播時(shí)間來(lái)率定下游糙率，水庫(kù)在相應(yīng)流量下的泄流到陽(yáng)春站的實(shí)際傳播時(shí)間為8.5 h，本次模擬相對(duì)泄流量傳播時(shí)間為8.2 h，兩者相差為3.5%，驗(yàn)證結(jié)果表明，河床糙率取值范圍為0.02～0.05，邊灘糙率取值范圍為0.03～0.055。

圖4 典型和設(shè)計(jì)洪水過(guò)程線

圖5 泄流能力曲線

2.3 地理模型

本次采取的地理模型為30 m×30 m的數(shù)字高程地形圖(DEM)，利用實(shí)際地理信息數(shù)據(jù)，在數(shù)字高程圖上勾勒出庫(kù)區(qū)，下游河道及框出可能淹沒(méi)范圍等二維網(wǎng)格。

圖6 水庫(kù)地理模型及與地理信息系統(tǒng)疊加

2.4 潰壩參數(shù)

潰壩參數(shù)主要是潰口尺寸、潰壩形式以及潰口成形時(shí)間，根據(jù)式(4)和式(6)，以1號(hào)副壩(土石壩)潰決為例，Bbottom=20 m，hb=15 m，z=0.5，潰壩方式為漫頂破壞，潰口形成時(shí)間為14 min，起潰高程取壩頂高程。

主壩作為面板堆石壩，潰口的發(fā)展與副壩土石壩不一樣[8]，面板堆石壩是逐級(jí)潰決的，洪水淘蝕造成上游面板階段性折斷，面板分段逐級(jí)潰決，隨著潰決速度加快最后發(fā)展為線性模式，主壩Bbottom=50 m，hb=15 m，z=0.5，潰壩方式為漫頂破壞，潰口形成時(shí)間為47 min，起潰高程取壩頂高程。各壩的主要尺寸見(jiàn)表2。

表2 各壩的主要尺寸 m

2.5 邊界條件

入流條件采用P=0.1%的設(shè)計(jì)洪水過(guò)程線進(jìn)行模擬，如圖4所示。

下邊界條件取決于下游水位，模擬時(shí)采用河床比降，本次邊界條件采用陽(yáng)春水文站下游的河床比降。一般來(lái)說(shuō)，當(dāng)發(fā)生大壩潰決事件時(shí)，潰壩洪水下泄流量非常大，下游的流量相對(duì)來(lái)說(shuō)非常小，故在初始條件中假設(shè)下游為干河床。即模型的初始表面高度設(shè)置為模擬區(qū)域的最小高程。

水庫(kù)初始水位設(shè)為校核洪水位115.98 m，模擬入庫(kù)洪水到來(lái)后水庫(kù)水位壅高過(guò)壩頂造成漫頂潰壩的場(chǎng)景。鑒于西山河河床平均比降較大，本次計(jì)算選用擴(kuò)散波格式的淺水方程加快計(jì)算速度。

2.6 洪水計(jì)算結(jié)果分析

潰壩發(fā)生后，壩址下游流量會(huì)在短時(shí)間內(nèi)發(fā)生激增，其峰頂流量向河道下游逐漸衰減，但仍遠(yuǎn)大于校核洪水流量，且控制斷面越靠近下游，其流量過(guò)程線上部越寬，形式呈“矮胖形”，說(shuō)明其洪水持續(xù)時(shí)間越長(zhǎng)，但峰值有所下降。結(jié)果顯示，潰壩發(fā)生后，距離大壩越近的地方，水位下降速率越快(0.29 m/min)，距離大壩越遠(yuǎn)的地方，水位下降速率越慢(0.16 m/min)，與事實(shí)相符。

主壩最大潰壩流速在壩址處為13.207 m/s，副壩最大潰壩流速在壩址處為12.714 m/s，部分受潰壩影響區(qū)域情況見(jiàn)表3。

表3 下游部分村莊受潰壩影響

大河水庫(kù)主壩潰壩洪水(主壩為堆石壩，不考慮全部潰決)發(fā)生后，洪水沿西山河主河道依次淹沒(méi)壩后電站、沿岸村莊(大河村、長(zhǎng)村、禾里村)以及工況企業(yè)、西山陂樞紐，并在西山河支流交匯處(禾里村)往上游回溯到圭崗鎮(zhèn)，往下游匯入漠陽(yáng)江主河道，淹沒(méi)合水鎮(zhèn)、陂面鎮(zhèn)以及部分陽(yáng)春市市區(qū)，并先后影響X600縣道、汕湛高速(羅陽(yáng)高速)以及三茂鐵路。

大河水庫(kù)副壩(副壩為土石壩，模擬全潰)下游為圭崗鎮(zhèn)，副壩潰壩將直接導(dǎo)致鎮(zhèn)區(qū)受淹。此后洪水將沿西山河支流經(jīng)陽(yáng)春合水匯入漠陽(yáng)江，出合水后淹沒(méi)范圍與主壩一致。

圖7為潰壩洪水風(fēng)險(xiǎn)分布，通過(guò)與QGIS聯(lián)合，展示了演進(jìn)過(guò)程中淹沒(méi)的下游城鎮(zhèn)、公路等設(shè)施，以及影響范圍，同時(shí)可以繪制出應(yīng)急轉(zhuǎn)移路線，為洪水防御提供決策參考。

圖7 潰壩洪水風(fēng)險(xiǎn)分布

3 結(jié) 論

本文利用 HEC-RAS 軟件建立了基于大河水庫(kù)DEM地形圖的二維潰壩模型進(jìn)行潰壩模擬，通過(guò)對(duì)潰口、下游及庫(kù)區(qū)斷面流量過(guò)程線進(jìn)行分析，證明了計(jì)算結(jié)果的合理性；通過(guò)QGIS軟件對(duì)模擬演算的潰壩淹沒(méi)范圍與社會(huì)行政、經(jīng)濟(jì)區(qū)劃圖層疊加，繪制出信息豐富直觀的洪水風(fēng)險(xiǎn)圖，為防洪決策提供有效支持。

HEC-RAS軟件可操作性強(qiáng)，二維洪水模擬兼顧了一維和三維模擬的優(yōu)點(diǎn)，其模擬速度快，模擬精度能滿足實(shí)際運(yùn)用要求，應(yīng)用該模型進(jìn)行洪水演進(jìn)計(jì)算可減輕工作量，同時(shí)與QGIS地理信息系統(tǒng)相結(jié)合，能快速直觀地反映洪水模擬結(jié)果。

對(duì)于模擬存在的誤差可能由于DEM精度不夠及潰壩流量過(guò)程的概化沒(méi)有完全符合實(shí)際情況，需要進(jìn)一步獲得下游精確斷面信息及精細(xì)模擬潰壩流量過(guò)程，以取得更好的模擬效果，為防汛減災(zāi)工作提供參考。