梯級橡膠壩塌壩調(diào)度及泄流特性數(shù)值研究

張倩楠，周 敏，郝瑞霞，胡孜軍，田 淳

（1.太原理工大學(xué)水利科學(xué)與工程學(xué)院，太原 030024；2.重慶開物工程咨詢有限公司，重慶 401122）

0 引 言

近年來，許多城市為改善城區(qū)生態(tài)環(huán)境和人居環(huán)境，進(jìn)行了河道美化工程建設(shè)，在河道上修建多級橡膠壩，實(shí)現(xiàn)層級蓄水，形成良好的城市景觀［1］。但汛期不合理調(diào)度會造成洪水疊加，影響下游安全，因此確定合理的閘壩調(diào)度策略顯得尤為重要［2］。

對于橡膠壩調(diào)度運(yùn)行，諸多學(xué)者根據(jù)工程運(yùn)行管理經(jīng)驗(yàn)及存在問題［3，4］提出運(yùn)行方式，或是對橡膠壩塌壩泄流量的計(jì)算［5-7］進(jìn)行研究分析。部分學(xué)者基于常規(guī)調(diào)度理論，利用程序搭建模型進(jìn)行調(diào)度尋優(yōu)，如鄧浩等［8］采用改進(jìn)遺傳算法，從基本試算轉(zhuǎn)向智能優(yōu)化，尋求全局最佳調(diào)度方案；張慶華等［9］編制梯級橡膠壩塌壩泄流仿真系統(tǒng)，對不同調(diào)度方案塌壩過程及塌壩泄量進(jìn)行預(yù)測；趙瑜琪等［10］綜合考慮塌壩時(shí)間及速度，對突發(fā)汛情的應(yīng)急調(diào)度方案展開研究等。上述調(diào)度方案的研究大都采用一維模型，對完全塌壩泄空的情況進(jìn)行探索。對于大頻率洪水，需要提前塌壩泄空以保證安全行洪，而對于較小頻率洪水，塌壩泄空會造成水資源浪費(fèi)，不利于洪后蓄水。目前對不完全塌壩的研究較少［11］。

本文以汾河太原城區(qū)段7 座橡膠壩為研究對象，本著安全行洪和蓄水兼顧的原則，進(jìn)行小頻率洪水不完全塌壩方面的研究。由于該段河道建筑物布置比較復(fù)雜，平面上有渾水渠和蓄水渠，加之橡膠壩、中隔堤和閘墩等，一維模型不能很好地反映這些建筑物的特點(diǎn)，而平面二維模型理論上可以有效地描述建筑物附近河道水流形態(tài)以及局部建筑物對水流的影響，且這方面研究成果有限。故利用MIKE21軟件建立能考慮塌壩過程的二維水動力學(xué)模型，分析20%頻率洪水來臨時(shí)橡膠壩不同塌壩方案泄水過程，以期為工程調(diào)度提供一定參考。

1 計(jì)算原理

1.1 控制方程

平面二維水動力學(xué)模型將雷諾時(shí)均的Navier-Stokes 方程進(jìn)行深度平均，形成二維淺水運(yùn)動方程，考慮紊動引入渦黏性系數(shù)，同時(shí)也考慮風(fēng)應(yīng)力和科氏力作用，其方程如下：

式中：h=η+d，為總水深，m；η表示水面到基準(zhǔn)面的高度，m；d為靜水深，m；S為點(diǎn)源流量，m3/s；u˙、v˙為基于水深平均的流速，m/s；us、vs為點(diǎn)源在x、y方向上的水流流速，m/s；g為重力加速度，m/s2；A為渦黏性系數(shù)，m2/s；ρ為水的密度，kg/m3；ρ0為水的相對密度，kg/m3；pa為當(dāng)?shù)卮髿鈮簭?qiáng)，Pa；f=2ωsinφ為科氏力系數(shù)；φ為地理緯度；ω為地球自轉(zhuǎn)角速度，rad/s；Sij為輻射應(yīng)力分量，N/m2；τsx、τsy為表面風(fēng)應(yīng)力，N/m2；τbx、τby為河床底部應(yīng)力分量，N/m2；M為曼寧系數(shù)，m1/3/s。

1.2 數(shù)值解法

計(jì)算區(qū)域劃分為不重疊的三角形及四邊形單元。采用有限體積法對控制方程進(jìn)行離散，變量位于單元中心，跨邊界通量垂直于單元邊，用近似Riemann 解法求得單元界面的對流通量。時(shí)間積分采用了低階顯式Euler方法。

2 模型構(gòu)建

2.1 模擬范圍

模擬范圍選擇為汾河蘭村至二壩段，全長約52 km，平均河寬約350 m，計(jì)算域面積18.371 km2。其中太原城區(qū)段為重點(diǎn)研究區(qū)域，河道平均底坡0.082%，有十條支溝匯入，經(jīng)一二期綜合治理后，共修建7 座橡膠壩，形成6 個(gè)蓄水池。橡膠壩工程區(qū)域建筑物和河道地形資料來源于具體設(shè)計(jì)報(bào)告及相關(guān)圖紙，其余河段基于分辨率為30 m 的DEM 數(shù)據(jù)進(jìn)行插值，模擬區(qū)域見圖1。

圖1 模擬區(qū)域和橡膠壩布置Fig.1 Simulated area and rubber dam arrangement

一期工程總長近6 km，每隔2 km 設(shè)一橡膠壩，共有1～4 號4 座橡膠壩，1 號橡膠壩前斷面尺寸見圖2 所示；二期工程在一期工程基礎(chǔ)上向南延伸7.6 km，有5～7號3 座橡膠壩。1～5號橡膠壩之間采用分槽蓄水，由中隔堤將河槽分為渾水渠和蓄水渠，當(dāng)泄水流量較小時(shí)，由渾水渠下泄，流量較大時(shí)，蓄水渠塌壩過流；5～7號橡膠壩之間為全槽蓄水。

圖2 1號橡膠壩處河道橫剖面Fig.2 Cross-section of the river at No.1 rubber dam

2.2 網(wǎng)格劃分

研究區(qū)域主河道和中隔堤處采用三角形非結(jié)構(gòu)化網(wǎng)格，其中河道網(wǎng)格面積300～1 000 m2之間，中隔堤網(wǎng)格面積5.5～13 m2；橡膠壩附近區(qū)域采用四邊形網(wǎng)格，網(wǎng)格尺寸為8 m×8 m。模型共生成網(wǎng)格88 130 個(gè)，圖3 為1 號橡膠壩處局部網(wǎng)格示意圖。

圖3 1號橡膠壩處網(wǎng)格劃分Fig.3 Grid division at No.1 rubber dam

2.3 初始條件及邊界條件

初始條件采用初始時(shí)刻的恒定流水位和流量。

模型上游邊界條件給定蘭村流量變化過程，下游邊界條件給定汾河二壩水位—流量關(guān)系。模型應(yīng)用時(shí)邊界條件設(shè)置如圖4所示。

圖4 模型邊界條件設(shè)置Fig.4 Model boundary condition setting

2.4 參數(shù)設(shè)置

曼寧系數(shù)M：利用觀測資料，通過數(shù)值驗(yàn)證進(jìn)行參數(shù)率定，最終選定M=35 m1/3/s。

渦黏性系數(shù)：采用Smagorinsky公式，取0.28 m2/s。

科氏力系數(shù)：根據(jù)所處地理緯度推求，為8.95×10-5。

模型計(jì)算過程中時(shí)間步長取1 s。

3 模型驗(yàn)證及參數(shù)率定

3.1 橡膠壩模擬效果校核計(jì)算

模型中橡膠壩處流量采用Villemonte 堰流公式計(jì)算，其描述如下：

式中：q為通過堰的流量，m3/s；W為寬度，m；C為堰流系數(shù)；k為堰指數(shù)；Hus為上游水位，m；Hds為下游水位，m；Hw為堰頂高程，m，堰頂高程等于堰底高程加堰高。

為了分析模型計(jì)算所用公式（5）的合理性，本文選擇工程區(qū)3種典型壩高：2.5、3、4 m 的橡膠壩進(jìn)行矩形水槽橡膠壩控水效果的驗(yàn)證分析。計(jì)算不同水位不同塌壩高度泄流量所得流量系數(shù)，與文獻(xiàn)［12］中模型試驗(yàn)數(shù)據(jù)進(jìn)行對比分析，率定公式中參數(shù)C取值的合理性。模型計(jì)算結(jié)果見圖5，模擬所得流量系數(shù)與文獻(xiàn)平均相對誤差為0.51%，認(rèn)為模擬結(jié)果可行，公式中堰流系數(shù)C選取1.838。

圖5 不同壩高流量系數(shù)驗(yàn)證結(jié)果Fig.5 Verification results on flow coefficient of different dam heights

3.2 河道恒定流水面線模擬計(jì)算

為率定模型曼寧系數(shù)，考慮橡膠壩全塌，泄流量為2 000 m3/s 的情況，進(jìn)行河道恒定流水面線計(jì)算。將模型水面線計(jì)算結(jié)果與傳統(tǒng)算法進(jìn)行比較分析，城區(qū)段河道水面線計(jì)算結(jié)果見圖6所示，兩者水深最大相對偏差4.56%，平均相對偏差2.23%，結(jié)果良好，模型曼寧系數(shù)選定35 m1/3/s。

圖6 恒定流水面線計(jì)算結(jié)果Fig.6 Constant flow water surface line calculation results

3.3 洪水演進(jìn)過程模擬結(jié)果驗(yàn)證

選取模擬河段2016年7月19日0∶00至21日0∶00的水文資料及橡膠壩塌壩運(yùn)行數(shù)據(jù)進(jìn)行模擬結(jié)果驗(yàn)證計(jì)算。河道水文資料基于上下游測站實(shí)測數(shù)據(jù)，各支流由控制流域面積采用類比的方法插補(bǔ)延展得到；橡膠壩塌壩運(yùn)行實(shí)際過程見圖7-9（a）。1 號、6 號及7 號壩前水位和流量變化過程計(jì)算結(jié)果與實(shí)測過程對比見圖7-9（b）、（c）所示。計(jì)算得各壩前水位平均相對偏差2.38%，流量平均相對偏差9.2%，驗(yàn)證結(jié)果良好。

圖7 1號橡膠壩模型驗(yàn)證Fig.7 Validation results at No.1 rubber dam

圖8 6號橡膠壩模型驗(yàn)證Fig.8 Validation results at No.6 rubber dam

圖9 7號橡膠壩模型驗(yàn)證Fig.9 Validation results at No.7 rubber dam

4 模型應(yīng)用

4.1 計(jì)算方案

選擇研究區(qū)域20%頻率洪水（洪峰流量為1 394 m3/s）過程線作為來水條件，分析7 座橡膠壩在不同塌壩組合情況下的過流特性。以洪水能在主河槽內(nèi)安全行洪和節(jié)約水資源，同時(shí)利于洪后蓄水為前提，根據(jù)河道過流能力及洪峰流量進(jìn)行初步計(jì)算，設(shè)計(jì)4種塌壩組合方案，進(jìn)行洪水演進(jìn)模擬。當(dāng)上游來流小于500 m3/s（渾水渠設(shè)計(jì)泄量）時(shí)，洪水直接通過渾水渠下泄；上游來流量大于500 m3/s時(shí)蓄水渠側(cè)開始塌壩。各壩塌壩高度如表1所示，塌壩速度按照設(shè)計(jì)標(biāo)準(zhǔn)為0.03 m/min。

表1 橡膠壩塌壩計(jì)算方案 mTab.1 Calculation scheme of rubber dam collapse

4.2 方案計(jì)算結(jié)果分析

4.2.1 洪水演進(jìn)過程分析

河道中洪水以波的形式傳播，當(dāng)發(fā)生20%頻率洪水時(shí)，4種塌壩方案洪水演進(jìn)過程類似。洪峰由蘭村到達(dá)1號橡膠壩處約3 h，到達(dá)7 號壩約4 h，到達(dá)二壩處約8 h。圖10 為典型方案特征斷面流量隨時(shí)間變化曲線，可知，洪水在傳播過程中，洪峰逐漸坦化，蘭村處洪峰流量為1 394 m3/s，到1 號橡膠壩處為1 381 m3/s，到達(dá)二壩處，約1 256 m3/s。

圖10 典型方案特征斷面流量過程線（方案B）Fig.10 Flow process of characteristic section in typical scheme（scheme B）

但不同塌壩方案下，洪水通過時(shí)1 號橡膠壩處流量分配情況不同，表2為洪峰通過1號橡膠壩時(shí)流量分配情況。隨著1號橡膠壩的塌壩高度的增大，阻水能力減弱，蓄水渠分配流量增大，渾水渠分配流量隨之減少。

表2 1號壩渾水渠及蓄水渠流量分配Tab.2 Flow distribution of muddy channels and storage channels of No.1 dam

為保證洪水能在主河道內(nèi)安全通過，洪峰到達(dá)時(shí)，各橡膠壩工程區(qū)內(nèi)水位應(yīng)低于內(nèi)堤高程，對各塌壩方案下橡膠壩處最高水位線進(jìn)行比較分析，繪制圖11 所示各方案最高水位包絡(luò)線圖。

圖11 不同方案工程區(qū)最高水位包絡(luò)線圖Fig.11 The highest water level envelope diagram in the project area under different schemes

方案A 中，洪峰通過時(shí)，渾水渠側(cè)水位均低于內(nèi)堤高程，蓄水渠側(cè)水位略高于內(nèi)堤高程，1號壩前水位超出內(nèi)堤高度較大，為0.06 m；方案B 與方案C 中，洪峰通過時(shí)渾水渠側(cè)及蓄水渠側(cè)水位均低于內(nèi)堤高程，且方案C 中渾水渠側(cè)水位低于中隔堤高程；方案D 中，各橡膠壩塌壩高度相同，下游橡膠壩對上游頂托作用較弱，除7 號橡膠壩外，各處水位均低于內(nèi)堤高程。綜上，方案B與方案C均滿足安全行洪要求。

4.2.2 工程區(qū)流場特性分析

工程區(qū)1 號壩和5 號壩處為中隔堤始端和末端，均有繞流現(xiàn)象出現(xiàn)，圖12 和圖13 為典型方案不同時(shí)刻1 號和5 號壩處局部流速矢量圖。T=17 h 和18 h 分別為1 號壩和5 號壩相應(yīng)洪峰出現(xiàn)時(shí)刻，相較而言此時(shí)繞流現(xiàn)象不大明顯，如圖12 和圖13（b）所示。在流量較小時(shí)刻，洪水多由渾水渠下泄，1 號壩和5號壩處出現(xiàn)明顯的繞流現(xiàn)象，如圖12和圖13（a）（c）所示。

圖12 典型方案不同時(shí)刻1號橡膠壩處流速矢量（方案B）Fig.12 Velocity vector at No.1 rubber dam of typical scheme at different time（scheme B）

圖13 典型方案不同時(shí)刻5號橡膠壩處流速矢量（方案B）Fig.13 Velocity vector at No.5 rubber dam of typical scheme at different time（scheme B）

4.2.3 方案比選分析

經(jīng)上述分析可知，方案B與方案C均能滿足安全行洪要求。綜合考慮洪后蓄水，方案B塌壩高度較小，雖局部稍有漫過中隔堤現(xiàn)象，但未漫上內(nèi)堤，既能滿足行洪要求，也更利于洪后蓄水，避免水體浪費(fèi)，為較優(yōu)方案。

方案B 洪峰通過時(shí)，各壩前斷面特征參數(shù)見表3所示，表中ΔZ1為壩前水位超出內(nèi)堤高度，ΔZ2為壩前水位超出中隔堤高度。1～5 號壩分槽蓄水，其中1～4 號壩的設(shè)計(jì)壩高、塌壩高度及河道底坡均相同，故在蓄水渠側(cè)水深和流速大致相當(dāng)，水深3.12～3.25 m 之間，流速1.14～1.16 m/s 之間；在渾水渠側(cè)，由于4號～5 號橡膠壩間底坡最陡，且5 號壩在中隔堤末端，洪峰通過時(shí)此處水深最小，流速較大，可達(dá)2.68 m/s。6 號橡膠壩立壩高度最大，為全槽蓄水，此處水深較大，最大4.11 m，流速1.07 m/s。

表3 方案B洪峰通過時(shí)壩前斷面特征參數(shù)匯總表Tab.3 Summary of the characteristic parameters of the pre-dam cross-section during the passage of the flood peak of Scheme B

5 結(jié) 論

本文以汾河上蘭村至二壩段為分析對象，構(gòu)建了考慮橡膠壩塌壩過程的二維洪水演進(jìn)模型。在模型驗(yàn)證的基礎(chǔ)上，針對小頻率洪水，對不完全塌壩組合方案進(jìn)行對比分析，主要結(jié)論如下：

（1）通過對橡膠壩模擬效果進(jìn)行校核，發(fā)現(xiàn)堰流系數(shù)C取1.838 合理，計(jì)算所得流量系數(shù)與試驗(yàn)比較，平均相對偏差0.51%。通過對河道進(jìn)行恒定流水面線模擬計(jì)算，曼寧系數(shù)選定35 m1/3/s，城區(qū)段河道水深平均相對偏差2.23%。在此基礎(chǔ)上進(jìn)行了考慮橡膠壩調(diào)度過程的實(shí)際洪水演進(jìn)過程模擬驗(yàn)證工作，計(jì)算所得各壩前水位和流量過程與實(shí)測過程比較，水位平均相對偏差2.38%，流量平均相對偏差9.2%，吻合良好。

（2）選擇研究區(qū)域20%頻率洪水過程作為來水條件，分析7座橡膠壩在4 種不完全塌壩組合情況下的過流特性，得出以下結(jié)論：洪水以波的形式傳播，且傳播過程中洪峰逐漸坦化，20%頻率洪水到達(dá)工程區(qū)1 號壩前約3 h；1 號橡膠壩的塌壩高度與蓄水渠內(nèi)的流量分配呈正相關(guān)關(guān)系；1 號壩和5 號壩分別為中隔堤始端與末端，會發(fā)生繞流現(xiàn)象，特別是流量較小時(shí)，大部分由渾水渠下泄，繞流現(xiàn)象愈加明顯；基于安全行洪及節(jié)約水資源的原則，對4 組塌壩方案下洪水通過工程區(qū)的過流特性進(jìn)行分析，方案B 為最優(yōu)方案，其1～7 號壩的塌壩高度分別為：1.1、1.1、1.1、1.1、1.8、1.8、2 m。上述結(jié)論有助于太原城區(qū)發(fā)生低頻率洪水時(shí)橡膠壩運(yùn)行規(guī)則制定提供參考，同時(shí)可為類似工程的模擬分析提供參考。