山區(qū)河流潰壩洪水演進(jìn)分析

王 雯， 董嘉銳， 楊 杰， 李 鵬， 李占斌， 劉基興， 朱戰(zhàn)齊

(1.西安理工大學(xué) 省部共建西北旱區(qū)生態(tài)水利國家重點(diǎn)實(shí)驗(yàn)室， 陜西 西安 710048; 2.大唐陜西發(fā)電有限公司石泉水力發(fā)電廠， 陜西 石泉 725200)

1 研究背景

山區(qū)洪水具有流速急，水量集中，突發(fā)性強(qiáng)的特點(diǎn)。一直以來，山洪災(zāi)害對(duì)世界上諸多山地國家的公共安全和社經(jīng)濟(jì)發(fā)展構(gòu)成巨大威脅[1]；山區(qū)又是建造大壩的天然基地，我國各類水庫數(shù)量從新中國成立前的1 200多座，增加到98 000多座，已形成初步的防洪減災(zāi)工程體系，當(dāng)壩體存在潰壩風(fēng)險(xiǎn)時(shí)，突發(fā)的潰壩洪水造成的安全威脅不可忽視。

針對(duì)潰壩后的洪水演進(jìn)過程，周興波等[2]基于潰口擴(kuò)展模型和一維潛水波模型對(duì)白格堰塞湖潰決洪水進(jìn)行了分析和對(duì)比，為堰塞湖應(yīng)急處置提供了參考。馬利平等[3]通過集成 HLLC 近似黎曼求解器的二維水動(dòng)力模型，對(duì)支溝潰壩洪水進(jìn)行了模擬，說明支溝潰壩洪水對(duì)主河道行洪具有一定影響。楊志等[4]通過耦合模型對(duì)黑河金盆水庫潰壩過程進(jìn)行了模擬，結(jié)果表明，耦合模型可提高模型計(jì)算的準(zhǔn)確性。王曉玲等[5]通過三維水動(dòng)力模型，對(duì)東武仕水庫潰壩洪水演進(jìn)進(jìn)行了數(shù)值模擬，并將傳統(tǒng)的壩區(qū)附近模擬范圍擴(kuò)大到下游城市建筑群整體范圍。這些研究均側(cè)重于潰決下泄流量過程[6]，數(shù)值計(jì)算方法[7-8]以及洪水波演進(jìn)[9-10]等方面，由于山區(qū)河流蜿蜒曲折，有明顯的區(qū)間匯流現(xiàn)象，河道寬窄變化，邊界條件復(fù)雜，目前對(duì)于長距離山區(qū)河流潰壩洪水演進(jìn)的過程及相關(guān)問題研究尚顯不足。鑒于此，本文以石泉大壩至喜河大壩段40 km長度河道為模擬對(duì)象，通過擬定大壩潰口形狀及潰壩時(shí)刻，獲得潰壩下泄流量過程線，建立區(qū)段平面二維計(jì)算模型，對(duì)典型山區(qū)型河道潰壩洪水演進(jìn)過程進(jìn)行了模擬，為防洪規(guī)劃和調(diào)洪決策提供依據(jù)。

2 數(shù)學(xué)模型及數(shù)值求解方法

2.1 平面二維洪水演進(jìn)數(shù)學(xué)模型

本文基于不可壓縮Reynolds值均布的Navier-Stokes方程，滿足Boussinesq假定和靜水壓力假定，建立洪水演進(jìn)模型。

二維非恒定淺水方程組[11]為：

(1)

(2)

(3)

式中：t為時(shí)間，s；h=η+d為總水深，m，其中d為靜止水深； m、η為水面高度，m；u,v分別為x,y方向速度分量，m/s；f為柯氏力系數(shù)，f=2ωsinφ，其中ω為自轉(zhuǎn)角速度、φ為緯度；g為重力加速度，m/s2；ρ為水的密度，kg/m3；sxx、sxy、syy分別為輻射應(yīng)力通量；Txx、Txy、Tyx、Tyy為黏滯應(yīng)力項(xiàng);S為源項(xiàng)矢量;us,vs為源項(xiàng)流速。

(4)

水平黏滯應(yīng)力Tij的表達(dá)式為：

(5)

在空間上對(duì)控制方程采用非結(jié)構(gòu)化網(wǎng)格有限體積法離散，可保證連續(xù)性方程和動(dòng)量方程的守恒；采用顯性歐拉法對(duì)時(shí)間進(jìn)行離散。

模型基本方程可表示為如下形式：

(6)

式中：U為守恒型物理向量；F為通量向量；S為源項(xiàng)；I為無黏性通量;V為黏性通量。

2.2 模型參數(shù)

計(jì)算河段劃分采用非結(jié)構(gòu)三角形網(wǎng)格，采用動(dòng)邊界處理技術(shù)[12]處理干濕邊界。河道糙率反映了計(jì)算河段的形態(tài)變化、邊界條件等因素的綜合影響，計(jì)算所采用的河道糙率主要由實(shí)測(cè)水流資料率定確定，分段、分高程對(duì)河道糙率進(jìn)行了調(diào)整，使得5年一遇及20年一遇洪水模擬條件下的水位與實(shí)測(cè)水位偏差小于1%，調(diào)整后河道糙率范圍為0.020～0.035。水平渦黏系數(shù)采用Smagorinsky公式計(jì)算，模型中采用值為0.28。垂向渦黏系數(shù)采用對(duì)數(shù)定律公式計(jì)算。

3 潰壩洪水演進(jìn)過程模擬

3.1 潰壩方案及流量確定

石泉大壩主壩壩型為混凝土空腹重力壩，最大壩高65 m，壩頂長度353 m，壩基巖石為石英片巖。壩體從左至右編號(hào)：1～3壩段為左岸非溢流壩段，4～6壩段為廠房壩段，7～23壩段為泄洪、排沙、溢流壩段，24～29壩段為右岸非溢流壩段。

本次計(jì)算采用瞬時(shí)潰壩模式，根據(jù)大壩安全監(jiān)測(cè)及運(yùn)行管理人員建議，潰壩范圍擬定為8～23壩段，潰壩缺口底部高程擬定為387 m，潰口寬度為170.5 m，潰口深度控制為29 m；另一方案為非溢流壩段潰決，潰口寬度為52 m，潰口控制深度為39 m，深度到達(dá)壩體建基面。潰決時(shí)刻為庫區(qū)水位達(dá)到擬定潰壩水位時(shí)刻，此時(shí)庫區(qū)來流為入庫洪水(設(shè)計(jì)洪水位410.29 m，校核洪水位413.67 m)。最大下泄流量采用寧利中等[13]提出的解析解求得，疊加入庫洪水及庫容計(jì)算潰壩下泄流量過程線，計(jì)算結(jié)果見圖1。計(jì)算工況見表1。

3.2 研究區(qū)域及網(wǎng)格劃分

二維數(shù)學(xué)模型計(jì)算河段總長約40 km，起于石泉水庫大壩，上至城關(guān)鎮(zhèn)，途經(jīng)石泉縣城、池河鎮(zhèn)、后柳鎮(zhèn)、上至城關(guān)鎮(zhèn)，喜河鎮(zhèn)，老喜河鎮(zhèn)、下至喜河水電站，見圖2。計(jì)算區(qū)域地形數(shù)據(jù)根據(jù)區(qū)段河道帶狀1∶1000地形圖進(jìn)行構(gòu)建，局部地區(qū)采用實(shí)測(cè)斷面數(shù)據(jù)進(jìn)行修正，計(jì)算網(wǎng)格節(jié)點(diǎn)數(shù)125 006個(gè)，網(wǎng)格單元數(shù)235 905個(gè)。為了盡可能準(zhǔn)確的反映區(qū)域流場(chǎng)，對(duì)河岸及模型邊界處網(wǎng)格進(jìn)行加密處理，計(jì)算區(qū)域局部河段網(wǎng)格分布見圖3。

圖1 各工況潰壩下泄流量過程線

表1 各工況計(jì)算參數(shù)

圖2 計(jì)算范圍水系圖及重點(diǎn)防護(hù)對(duì)象分布圖

圖3 研究區(qū)域局部網(wǎng)格結(jié)構(gòu)圖

3.3 結(jié)果分析

3.3.1 潰壩洪水演進(jìn)過程 潰壩洪水演進(jìn)的關(guān)鍵研究?jī)?nèi)容之一是洪峰的傳播到達(dá)時(shí)間，精確掌握潰壩洪水到達(dá)時(shí)間，可有效進(jìn)行防洪預(yù)警，減少洪水損失。工況2和6的特征斷面平均流速隨時(shí)間變化圖見圖4和5。由圖4、5可看出，溢流壩段潰壩時(shí)，由于潰壩增加流量相對(duì)正常泄洪流量較小，下游特征斷面流速增大不顯著；當(dāng)非溢流壩段潰壩，斷面越靠近大壩，流速峰值越明顯，潰壩發(fā)生于洪水過程第45 h，石泉大壩下游G210線石泉漢江大橋處斷面平均流速可達(dá)5.28 m/s，但由于石泉縣城段河道順直且河寬較寬，壩體下游5 km處河道自然束窄且形成彎道，為一天然卡口河道，壩體至彎道間形成滯洪區(qū)，彎道后十天高速以及安陽鐵路漢江6號(hào)橋河段流速增大較為平緩，斷面平均流速峰值相比潰壩下游1 km范圍內(nèi)削減明顯。由各斷面峰值出現(xiàn)的時(shí)間可獲得潰壩傳播過程，潰壩水流在研究區(qū)域整體傳播時(shí)長為43 min，由于河道支流匯入較多，且大洪水條件下主流河水倒灌入支流，遲滯并削減了洪峰向下游的傳播。

3.3.2 水位特征分析 選取具有代表性的石泉縣城城區(qū)段河道進(jìn)行說明，圖6、7分別為工況1、工況2潰壩發(fā)生40 min時(shí)洪水水位分布圖。工況1最高水位為378.0 m，工況2最高水位為382.0 m，高水位均位于饒峰河與漢江交匯口左岸區(qū)域。主要原因?yàn)椋菏髩涡购樗黜敍_點(diǎn)位于該部位，部分水流動(dòng)能轉(zhuǎn)化為勢(shì)能，并逐漸向兩側(cè)擴(kuò)散；下泄洪水在饒峰河河口形成逆時(shí)針環(huán)流，環(huán)流頂托饒峰河支流水位，并在支流形成倒灌；石泉縣城段下游河道急劇束窄，并形成彎道，彎道凹岸頂沖點(diǎn)水位明顯抬升，并在彎段下游逐漸降低，石泉縣城段潰壩產(chǎn)生的洪水波波峰在經(jīng)過下游彎道卡口段后被削減。

圖4工況2特征斷面潰壩洪水流速隨時(shí)間變化過程 圖5工況6特征斷面潰壩洪水流速隨時(shí)間變化過程

圖6工況1潰壩后40min石泉縣城洪水水位分布圖 圖7工況2潰壩后40min石泉縣城洪水水位分布圖

圖8、9分別為工況2潰壩發(fā)生40、43 min計(jì)算區(qū)域水位分布圖。圖8與9表明，由于石泉大壩下游11 km處池河匯入，池河匯入方向與漢江主河道正交，且交匯口位于漢江彎段左岸凹岸，一方面使得潰壩洪水流量倒灌池河，另一方面頂沖點(diǎn)處山體嚴(yán)重頂托下泄洪水，壅高上游段漢江河道水位及池河水位。計(jì)算模型整體分析，由于山區(qū)河道的彎曲和不規(guī)則形態(tài)，洪水演進(jìn)沿程阻力大，受邊界條件影響顯著，使得潰壩洪峰在傳播過程中迅速坦化，對(duì)下游喜河電站影響作用逐漸減小，表明河道的蜿蜒能夠很好地降低遠(yuǎn)處的受災(zāi)。

3.3.3 流速分布 以工況2為例進(jìn)行整體分析，工況2潰壩發(fā)生40 min計(jì)算區(qū)域流速分布見圖10。

由圖10可看出，潰壩發(fā)生后大壩下游河道流速迅速增加，但由于河流較長、彎道眾多且主河道兩側(cè)存在有支流及支毛溝，潰壩洪峰流量及能量在傳播過程中由于支流倒灌以及交匯口處的環(huán)流消耗水量及能量，使洪峰流量過程逐漸坦化，水流流速減小。除潰口段水流流速較大外，其余河段并無出現(xiàn)6 m/s以上斷面平均流速，各工況河段最大流速介于4～6 m/s之間。

以老喜河鎮(zhèn)段河道為例，工況5潰壩發(fā)生40 min老喜河鎮(zhèn)段河道洪水流速矢量分布如圖11。此河段洪水由于支流與主流間產(chǎn)生眾多的擴(kuò)散回流區(qū)，支溝水流與主流進(jìn)行動(dòng)量交換，削減主流能量，但水流橫向切割沖刷岸坡，可能導(dǎo)致邊灘、河岸崩塌，河道拓寬，應(yīng)注意重點(diǎn)位置處的防護(hù)。局部河段兩岸束窄處岸坡存在有突出的陡峭巖體，造成岸坡段形成較高的流速值，由圖11可見，局部流速可達(dá)5 m/s。

圖8工況2潰壩發(fā)生40min圖9工況2潰壩發(fā)生43min圖10工況2潰壩發(fā)生40min

計(jì)算區(qū)域水位分布圖 計(jì)算區(qū)域水位分布圖 計(jì)算區(qū)域流速分布圖

圖11 工況5潰壩發(fā)生40 min老喜河鎮(zhèn)段河道洪水流速矢量分布圖

4 結(jié) 論

本文通過建立山區(qū)河道平面二維水動(dòng)力學(xué)模型，對(duì)模型參數(shù)進(jìn)行了率定，模擬了6種工況下的潰壩洪水演進(jìn)過程，對(duì)復(fù)雜的地形進(jìn)行了精細(xì)化建模并考慮區(qū)間匯流，得出以下結(jié)論：

(1)山區(qū)河流形態(tài)對(duì)潰壩水流具有顯著的影響作用，卡口段上游河段易形成局部的滯洪區(qū)，通過調(diào)蓄洪水，影響下泄流量過程。河流彎段對(duì)調(diào)整水流結(jié)構(gòu)起著至關(guān)重要的作用，彎道的存在，極大削減了潰壩洪水波波峰；

(2)山區(qū)河流的支流以及存在的眾多支毛溝，在潰壩洪水發(fā)生的條件下，能有效吸納洪峰流量，在交匯處形成環(huán)流，削減主流能量，但應(yīng)注意回流區(qū)可能導(dǎo)致的岸坡沖刷失穩(wěn)；

(3)各潰壩工況斷面平均流速最大值未超過6m/s，但局部范圍內(nèi)由于河道邊界形狀突變，可能產(chǎn)生局部高流速區(qū)，工程防護(hù)時(shí)應(yīng)加以注意。