花涼亭水庫潰壩洪水演進研究

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

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

0 前 言

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

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

1 工程概況

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

2 壩址潰口流量計算

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

2.1 潰壩方案設定

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

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

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

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

2.2 潰壩模式

(1) 管涌潰壩模式

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

B0=Bry

(1)

(2)

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

(2) 漫頂潰壩模式

漫頂導致潰壩的水流侵蝕，河渠中的水流流量用寬頂堰關系如下：

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個管涌導致潰壩的潰口洪水流量過程線呈現(xiàn)相同的變化趨勢，即在潰壩瞬間，洪水流量迅速達到次高峰，再迅速下降一段后上漲至最高峰，最終逐漸下降并趨于平緩，洪峰流量均出現(xiàn)在潰壩后約5 h處，其中，遭遇1000年一遇設計洪水的潰口洪峰流量為47 091 m3/s，小于10000年一遇洪水的潰口洪峰流量為55 443 m3/s；漫頂導致潰壩的洪水流量過程線則是在庫水位未達到防浪墻頂高程前，潰口無流量；當庫水位超過防浪墻頂高程后，過壩洪水流量迅速達到峰值，之后迅速下降并逐漸趨于平緩，潰口洪峰流量達66 213 m3/s。

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

圖1 遭遇設計洪水工況下管涌導致潰壩的潰口流量過程線

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

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

3 潰壩下泄洪水演進計算

3.1 計算方程

采用丹麥水利科學研究院(DHI)研發(fā)的MIKE21軟件構建二維水動力模型。模型基本方程為二維淺水方程組，方程組包含應力方程見公式(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ù)水系情況在花涼亭水庫壩址至皖河入長江口的計算范圍中設置5個控制斷面，分別為太湖水文站(距壩址3.20 km)、黑河匯合口(距壩址31.86 km)、長河河口(距壩址48.38 km)、石牌水文站(距壩址54.10 km)、江鎮(zhèn)(距壩址66.78 km)，洪水淹沒范圍內控制斷面見圖4。

圖4 潰壩洪水計算淹沒范圍內控制斷面

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

潰壩洪水下泄演進計算需要收集的資料主要包括研究區(qū)域河道數(shù)據(jù)與地形數(shù)據(jù)、研究區(qū)域高程數(shù)據(jù)、上游邊界洪水過程線、下游邊界水位過程線等。在完成資料處理后，水庫下泄洪水演進計算模型采用非結構網(wǎng)格進行網(wǎng)格劃分，不受網(wǎng)格節(jié)點約束，能合理分布網(wǎng)格并契合地形，適合復雜地形區(qū)域的網(wǎng)格劃分。本次對研究區(qū)域進行網(wǎng)格剖分，對地形網(wǎng)格文件進行高程插值，得到研究區(qū)域的最終地形網(wǎng)格文件。網(wǎng)格在庫區(qū)周圍區(qū)域進行加密，對影響較小區(qū)域進行適度稀疏，以保證計算精度的同時提高計算效率。地形網(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ù)研究內容對區(qū)域的邊界進行劃分，設定上游邊界條件的花涼亭水庫入庫洪水過程線。根據(jù)水庫工程情況，當遭遇洪水需要排洪時，花涼亭水庫發(fā)電輸水洞、泄洪洞和溢洪道將按照防汛水庫控制運用原則實行泄洪，此時庫區(qū)淹沒計算模型也將花涼亭水庫調洪演算獲得的壩前對應水位過程線設置為上游邊界條件，具體設置根據(jù)模擬工況確定。

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

3.4 模型耦合分析

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

3.5 潰壩下泄洪水演進計算結果及分析

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

3.5.1洪峰流量

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

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

3.5.2潰壩洪水到達時間

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

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

3.5.3最大淹沒水深

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

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

3.5.4潰壩后果估算

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

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

4 結 論

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

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

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

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

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