花涼亭水庫潰壩洪水演進(jìn)研究

徐子滿，周克發(fā)，楊德瑋，張文東

(1.安慶市花涼亭水庫管理處，安徽 太湖 246400；2.南京水利科學(xué)研究院大壩安全與管理研究所，南京 210029；3.福州大學(xué)，福州 350000)

0 前 言

隨著社會經(jīng)濟(jì)發(fā)展、人口數(shù)量增長以及城市化水平不斷提高，水庫大壩一旦潰決，后果極其嚴(yán)重。2015年，水利部發(fā)布了SL/Z 720 -2015《水庫大壩安全管理應(yīng)急預(yù)案編制導(dǎo)則》。該導(dǎo)則規(guī)定為提升水庫大壩安全管理水平，確保大壩安全運行，大中型水庫需針對大壩突發(fā)事件，制定科學(xué)合理的大壩安全管理應(yīng)急預(yù)案[1]。該預(yù)案需要對潰壩洪水事件及其后果進(jìn)行分析，包括壩址潰口洪水計算、下泄洪水演進(jìn)、淹沒范圍(主要包括下游淹沒區(qū)的洪峰流量、洪水達(dá)到時間及淹沒水深等參數(shù))、潰壩后果計算(主要包括生命損失估算、經(jīng)濟(jì)損失估算及社會與環(huán)境影響分析)。為提高花涼亭水庫大壩安全管理應(yīng)急預(yù)案的可行性、有效性，有必要對該水庫潰壩事件開展壩址潰口洪水計算、下泄洪水演進(jìn)及淹沒范圍分析研究。

現(xiàn)階段針對潰壩計算研發(fā)的數(shù)值計算軟件和商業(yè)程序都較為成熟，其中賀娟[2]等采用Hec-Ras建立長河壩水庫漫頂一維潰壩模型，研究其對下游黃金坪、瀘定水電站連續(xù)潰壩的下游風(fēng)險分析；楊忠勇等[3]通過MIKE21FM模型對4座梯級水庫下游城市的淹沒過程進(jìn)行了計算和討論分析；周清勇等[4-5]將BREACH模型成功應(yīng)用于七一水庫風(fēng)險分析；程坤等[6]以西藏藏木水電站為例，采用Mike21軟件進(jìn)行潰壩洪水模擬；楊德瑋等[7]通過BREACH-MIKE21耦合模型分析在保障水庫大壩安全運行及風(fēng)險管理方面發(fā)揮了積極作用。考慮到花涼亭水庫大壩為粘土心墻砂殼壩，本文采用BREACH模型計算出潰口流量與潰口尺寸，利用MIKE21中Flow Model模塊生成潰壩洪水模擬文件，以此分析花涼亭水庫潰壩對下游造成的影響，為水庫大壩安全管理應(yīng)急預(yù)案提供技術(shù)支持。

1 工程概況

花涼亭水庫位于安徽省安慶市太湖縣城西北約5 km處長江流域皖河支流長河上游，壩址以上控制流域面積1 870 km2，總庫容23.66億m3，是一座以防洪、灌溉為主，兼有發(fā)電、供水、養(yǎng)殖、航運、旅游等綜合利用的 Ⅰ 等大(1)型水庫，工程主要建筑物級別為1級；正常蓄水位88.00 m，汛限水位85.50 m，死水位74.00 m；1000年一遇設(shè)計洪水位95.21 m，10000年一遇校核洪水位97.30 m。主要建筑物包括大壩、溢洪道、泄洪隧洞、引水隧洞及電站。大壩為粘土心墻砂殼壩，壩頂高程99.25 m，最大壩高58.0 m，壩頂長566 m，壩頂寬7.57 m，防浪墻頂高程100.45 m。溢洪道控制段為有閘控制WES-Ⅰ 型實用堰，堰頂高程82.80 m，共8孔，單孔凈寬12.0 m，每孔設(shè)1扇12 m×9 m(寬×高)弧形工作鋼閘門。泄洪隧洞進(jìn)口段閘槽2孔，每孔尺寸3.5 m×8.0 m(寬×高)，設(shè)2扇潛孔式事故檢修鋼閘門；出口閘室段底板高程53.60 m，設(shè)1扇5.0 m×5.8 m(寬×高)弧形工作鋼閘門。引水隧洞為灌溉、發(fā)電及放空水庫共用洞，進(jìn)口段渠首設(shè)頂高程51.50 m鋼筋混凝土攔沙坎；進(jìn)口閘槽2孔，底板高程43.00 m，設(shè)2扇事故檢修鋼閘門；隧洞出口閘室段底板高程42.15 m，設(shè)1扇4.5 m×4.5 m(寬×高)弧形工作閘門。水庫設(shè)計灌溉面積7萬hm2，防洪保護(hù)下游安慶市太湖縣、潛山市、懷寧縣、望江縣、大觀區(qū)，包括26個鄉(xiāng)鎮(zhèn)和大觀開發(fā)區(qū)，保護(hù)面積956.44 km2，涉及人口69.45萬人、耕地4萬hm2，長河兩岸圩區(qū)、潛水和皖水下游圩區(qū)、皖河干流右岸望江縣新東隔堤以東的同馬大堤保護(hù)圩區(qū)、皖河干流左岸沿河圩區(qū)，以及合安九高鐵、G50滬渝高速、G35濟(jì)廣高速等重要基礎(chǔ)設(shè)施。

2 壩址潰口流量計算

根據(jù)花涼亭水庫入庫洪水過程線、大壩縱橫斷面、壩體填筑材料與壩基巖土層的物理性質(zhì)指標(biāo)與力學(xué)參數(shù)等，采用BREACH軟件計算潰口流量過程[9]。該軟件是一套專門用于土石壩滲透破壞和漫頂破壞潰口洪水(包括潰口流量過程線及潰口尺寸)模擬計算的專業(yè)軟件，由D.L.Fread于1988年提出，是目前應(yīng)用較為廣泛的大壩潰決模型，在實際工程應(yīng)用中可預(yù)測最終潰口發(fā)展過程和最終尺寸，并計算得到潰口洪水流量過程線和水庫水位過程線。

2.1 潰壩方案設(shè)定

根據(jù)花涼亭水庫現(xiàn)場查看情況和已失事的土石壩失事經(jīng)驗，花涼亭水庫可能發(fā)生的洪水突發(fā)事件包括超標(biāo)準(zhǔn)泄洪和潰壩洪水。本文主要研究潰壩洪水演進(jìn)，故考慮花涼亭水庫主要可能發(fā)生的潰壩洪水，設(shè)定3種計算方案：

(1) 庫水位82.80 m時，遭遇1000年一遇設(shè)計洪水，泄洪設(shè)施正常泄洪，大壩發(fā)生管涌，導(dǎo)致潰壩；

(2) 庫水位82.80 m時，遭遇10000年一遇校核洪水，泄洪設(shè)施正常泄洪，大壩發(fā)生管涌，導(dǎo)致潰壩；

(3) 庫水位82.80 m時，遭遇10000年一遇校核洪水，泄洪設(shè)施不能正常泄洪，逼高庫水位，大壩發(fā)生漫頂，導(dǎo)致潰壩。

2.2 潰壩模式

(1) 管涌潰壩模式

模擬管涌潰決時，須假設(shè)初始矩形河渠遭受侵蝕管涌后尺寸的中心線高程，并保證水庫水位大于此中心線高程。同時，管涌通道底部受到向下的垂直向侵蝕，其頂部也存在向上的同樣大小的垂直向侵蝕。潰口寬度計算見公式(1)，管涌通道的流量計算見公式(2)：

B0=Bry

(1)

(2)

公式(1)～(2)式中:Br為基于最合適河渠水力有效作用的一個因子,對于漫頂破壞，參數(shù)Br值為2，對于管涌破壞，參數(shù)Br值設(shè)置為1.0;Qb為通過管涌通道的流量,m3/s；A為潰口橫斷面面積,m2；g為重力加速度,取9.8 m/s2；Hp為中心線高程,m；(H-Hp)為潰口靜態(tài)水頭,m；L為管涌通道長度,m；D為管涌通道直徑或?qū)挾?m；f為摩擦因數(shù)。

(2) 漫頂潰壩模式

漫頂導(dǎo)致潰壩的水流侵蝕，河渠中的水流流量用寬頂堰關(guān)系如下：

Qb=3B0(H-Hc)

(3)

公式(3)中：Qb為潰口河渠中的流量，m3/s；B0為初始矩形形狀河渠的瞬時寬度，m；H為壩前水位高程,m；Hc為潰口底部高程,m。

2.3 潰口流量分析

經(jīng)計算得到3種計算方案的潰口尺寸、最大下泄流量(見表1)及洪水流量過程線(見圖1～3)，由圖1～3可見2個管涌導(dǎo)致潰壩的潰口洪水流量過程線呈現(xiàn)相同的變化趨勢，即在潰壩瞬間，洪水流量迅速達(dá)到次高峰，再迅速下降一段后上漲至最高峰，最終逐漸下降并趨于平緩，洪峰流量均出現(xiàn)在潰壩后約5 h處，其中，遭遇1000年一遇設(shè)計洪水的潰口洪峰流量為47 091 m3/s，小于10000年一遇洪水的潰口洪峰流量為55 443 m3/s；漫頂導(dǎo)致潰壩的洪水流量過程線則是在庫水位未達(dá)到防浪墻頂高程前，潰口無流量；當(dāng)庫水位超過防浪墻頂高程后，過壩洪水流量迅速達(dá)到峰值，之后迅速下降并逐漸趨于平緩，潰口洪峰流量達(dá)66 213 m3/s。

表1 潰口尺寸及潰口最大下泄流量計算結(jié)果

圖1 遭遇設(shè)計洪水工況下管涌導(dǎo)致潰壩的潰口流量過程線

圖2 遭遇校核洪水工況下管涌導(dǎo)致潰壩的潰口流量過程線

圖3 遭遇校核洪水工況下漫頂導(dǎo)致潰壩的潰口流量過程線

3 潰壩下泄洪水演進(jìn)計算

3.1 計算方程

采用丹麥水利科學(xué)研究院(DHI)研發(fā)的MIKE21軟件構(gòu)建二維水動力模型。模型基本方程為二維淺水方程組，方程組包含應(yīng)力方程見公式(4)、連續(xù)方程見公式(5)和動量方程見公式(6)，具體如下：

h=η+d

(4)

(5)

(6)

3.2 計算控制斷面

長河發(fā)源于岳西縣多枝尖，流經(jīng)冶溪、桃陽兩鄉(xiāng)后，在桃陽鄉(xiāng)的梅子林入太湖縣境，在牛鎮(zhèn)龍灣注入花涼亭水庫；長河與潛水在陶灣匯合形成皖河干流，于橫壩頭納皖水，過石牌鎮(zhèn)，在皖河閘(距皖河口約20 km)納武昌湖來水，最終在安慶市上游西郊注入長江。根據(jù)水系情況在花涼亭水庫壩址至皖河入長江口的計算范圍中設(shè)置5個控制斷面，分別為太湖水文站(距壩址3.20 km)、黑河匯合口(距壩址31.86 km)、長河河口(距壩址48.38 km)、石牌水文站(距壩址54.10 km)、江鎮(zhèn)(距壩址66.78 km)，洪水淹沒范圍內(nèi)控制斷面見圖4。

圖4 潰壩洪水計算淹沒范圍內(nèi)控制斷面

3.3 地形網(wǎng)格剖分處理及邊界條件和參數(shù)設(shè)定

潰壩洪水下泄演進(jìn)計算需要收集的資料主要包括研究區(qū)域河道數(shù)據(jù)與地形數(shù)據(jù)、研究區(qū)域高程數(shù)據(jù)、上游邊界洪水過程線、下游邊界水位過程線等。在完成資料處理后，水庫下泄洪水演進(jìn)計算模型采用非結(jié)構(gòu)網(wǎng)格進(jìn)行網(wǎng)格劃分，不受網(wǎng)格節(jié)點約束，能合理分布網(wǎng)格并契合地形，適合復(fù)雜地形區(qū)域的網(wǎng)格劃分。本次對研究區(qū)域進(jìn)行網(wǎng)格剖分，對地形網(wǎng)格文件進(jìn)行高程插值，得到研究區(qū)域的最終地形網(wǎng)格文件。網(wǎng)格在庫區(qū)周圍區(qū)域進(jìn)行加密，對影響較小區(qū)域進(jìn)行適度稀疏，以保證計算精度的同時提高計算效率。地形網(wǎng)格文件共剖分得到網(wǎng)格73 495個。計算面積302.70 km2，網(wǎng)格平均面4 119 m2，最大網(wǎng)格面積48 346 m2，最小網(wǎng)格面積621 m2，網(wǎng)格劃分如圖5所示。

圖5 網(wǎng)格劃分

根據(jù)研究內(nèi)容對區(qū)域的邊界進(jìn)行劃分，設(shè)定上游邊界條件的花涼亭水庫入庫洪水過程線。根據(jù)水庫工程情況，當(dāng)遭遇洪水需要排洪時，花涼亭水庫發(fā)電輸水洞、泄洪洞和溢洪道將按照防汛水庫控制運用原則實行泄洪，此時庫區(qū)淹沒計算模型也將花涼亭水庫調(diào)洪演算獲得的壩前對應(yīng)水位過程線設(shè)置為上游邊界條件，具體設(shè)置根據(jù)模擬工況確定。

水庫下泄洪水演進(jìn)計算模型在計算前需要分別確定多類參數(shù)。參數(shù)設(shè)定在MIKE21軟件界面進(jìn)行設(shè)置，依次設(shè)定時間參數(shù)包括時間步數(shù)(10 800)、時間步長(1 s)、總時長(30 h)；選定水動力和內(nèi)陸洪水模塊。內(nèi)陸洪水的水動力求解器不包括那些對海洋和海岸非常重要的功能，如科氏力的影響、潮汐的潛在影響、二維區(qū)域風(fēng)應(yīng)力的影響、波輻射的影響、湍流的影響等。水動力模塊參數(shù)則包括干濕邊界(干水深0.005 m；淹沒水深0.050 m；濕水深0.100 m)、底床阻力(曼寧系數(shù)取25 m1/3/s)、初始條件(隨空間變化的水位)、邊界條件的設(shè)定。

3.4 模型耦合分析

基于水庫大壩的相關(guān)參數(shù)，通過BREACH模型計算獲得潰口流量過程與潰口形態(tài)變化，計算所得結(jié)果導(dǎo)入MIKE21模型中作為洪源要素。耦合分析的計算次序是迭代的，進(jìn)入潰口的水流取決于潰口的底部高程及其寬度，而潰口屬性特征取決于潰口的尺寸和流量。BREACH初次計算出潰口流量與潰口尺寸，作為洪源數(shù)據(jù)導(dǎo)入MIKE21中模擬運行，提取mesh文件中潰口所在網(wǎng)格淹沒水深并得出計算歷時下的潰口流量，若此流量與BREACH計算所得潰口流量差值小于0.5 m3/s，即滿足精度要求，所得結(jié)果可輸出后處理；若差值大于0.5 m3/s，須將潰口流量與初次計算潰口尺寸迭代至所得結(jié)果滿足精度要求。

3.5 潰壩下泄洪水演進(jìn)計算結(jié)果及分析

針對花涼亭水庫潰壩下泄洪水演進(jìn)過程，本文采用不同工況下淹沒面積和生命財產(chǎn)損失及5個控制斷面處的洪峰流量、洪水到達(dá)時間、最大淹沒水深反映下泄洪水的破壞性。

3.5.1洪峰流量

從表2可以看出，對比同為管涌潰壩模式的計算方案1和方案2，隨著洪水量級增加，洪峰流量增加；對比洪水量級相等的計算方案2和方案3，漫頂潰壩模式的洪峰流量多于管涌潰壩模式的洪峰流量，與實際相符。對比同一計算方案下不同控制斷面的洪峰流量變化發(fā)現(xiàn)，從太湖水文站至長河河口，潰壩下泄洪水洪峰流量從47 074 m3/s降至27 144 m3/s，表明隨著潰壩洪水向下游推進(jìn)，河道調(diào)蓄作用致使洪峰流量呈現(xiàn)遞減趨勢；從長河河口至石牌水文站，潰壩下泄洪水洪峰流量從27 114 m3/s增至32 571 m3/s，位于下游的石牌水文站的洪峰流量反而更高，這是由于汛期下游長江水位的頂托作用造成江水倒灌，導(dǎo)致江鎮(zhèn)至石牌水文站的河段洪峰流量遞增，出現(xiàn)石牌水文站處的洪峰流量高于長河河口的現(xiàn)象。

表2 不同潰壩洪水計算方案下控制斷面洪峰流量

3.5.2潰壩洪水到達(dá)時間

從表3可以看出，對比同為管涌潰壩模式的計算方案1和方案2，隨著洪水量級增加，潰壩洪水到達(dá)控制斷面的時間提前；對比洪水量級相等的計算方案2和方案3，漫頂潰壩模式下潰壩洪水到達(dá)控制斷面的時間早于管涌潰壩模式，與實際相符。對比同一計算方案下潰壩洪水到達(dá)不同控制斷面的時間變化發(fā)現(xiàn)，從太湖水文站至江鎮(zhèn)，洪水到達(dá)各控制斷面的時間從0.15 h增至7.75 h，與控制斷面至花涼亭水庫壩址的距離成正比，距壩址越遠(yuǎn)，洪水到達(dá)時間越長。

表3 不同潰壩洪水計算方案下控制斷面洪水達(dá)到時間

3.5.3最大淹沒水深

從表4可以看出，對比同為管涌潰壩模式的計算方案1和方案2，隨著洪水量級增加，同一控制斷面處最大淹沒水深增加；對比洪水量級相等的計算方案2和方案3，漫頂潰壩模式下同一控制斷面處的最大淹沒水深高于管涌潰壩模式，與實際相符。對比同一計算方案下不同控制斷面處的最大淹沒水深變化發(fā)現(xiàn)，由于不同控制斷面底部地形不一致，最大淹沒水深無明顯規(guī)律。

表4 不同潰壩洪水計算方案下控制斷面最大淹沒水深

3.5.4潰壩后果估算

基于模型計算淹沒范圍、水深和淹沒歷時等信息，結(jié)合淹沒區(qū)域內(nèi)各市、縣、區(qū)的社會經(jīng)濟(jì)情況，綜合估算洪水影響程度并評估洪水淹沒損失如表5所示。由表5可以看出，從工況1至工況3，潰壩洪水淹沒面積從916.45 km2遞增至956.44 km2；預(yù)估受災(zāi)人口從63.87萬人遞增至69.45萬人；預(yù)估損失GDP從265.14億元遞增至287.54億元。可見在同一潰壩模式下，洪水量級越大，淹沒面積、受災(zāi)人口、損失GDP越大；同一洪水量級下，漫頂潰壩下的淹沒面積、受災(zāi)人口、損失GDP多于管涌潰壩，符合自然規(guī)律。從整體評估淹沒損失，3種潰壩工況下洪水淹沒面積均大于900 km2，預(yù)估受災(zāi)人口均多于60萬人，預(yù)估損失GDP均高于250億元，可見花涼亭水庫潰壩將對下游造成及其嚴(yán)重的后果。

表5 不同潰壩洪水計算方案下?lián)p失估算

4 結(jié) 論

采用BREACH-MIKE21耦合模型，對花涼亭水庫3種潰壩洪水工況的壩址潰口流量及潰壩下泄洪水演進(jìn)進(jìn)行分析計算，得到結(jié)論如下：

(1) 花涼亭水庫遭遇10000年一遇校核洪水，泄洪設(shè)施不能正常泄洪，逼高庫水位，大壩發(fā)生漫頂導(dǎo)致潰壩為最不利潰壩工況，其潰口最大下泄流量達(dá)到66 213 m3/s。

(2) 在最不利潰壩工況下，由于水庫下游長江水位頂托作用導(dǎo)致江水倒灌，導(dǎo)致距離壩址較遠(yuǎn)處石牌水文站的洪峰流量大于距離壩址較近處長河河口的洪峰流量。

(3) 在最不利潰壩工況下，水庫下游潰壩洪水淹沒面積為956.44 km2，預(yù)估受災(zāi)人口接近69.45萬人，預(yù)估損失GDP達(dá)到287.54億元。

計算結(jié)果符合花涼亭水庫工程實際情況，并已在花涼亭水庫大壩安全管理應(yīng)急預(yù)案編制中得到了采納,這對于潰壩洪水災(zāi)害預(yù)防以及提高應(yīng)急預(yù)案的可行性、有效性，具有重要的技術(shù)支撐作用。