梯級橡膠壩塌壩調度及泄流特性數值研究

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

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

0 引 言

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

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

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

1 計算原理

1.1 控制方程

平面二維水動力學模型將雷諾時均的Navier-Stokes 方程進行深度平均，形成二維淺水運動方程，考慮紊動引入渦黏性系數，同時也考慮風應力和科氏力作用，其方程如下：

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

1.2 數值解法

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

2 模型構建

2.1 模擬范圍

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

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

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

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

2.2 網格劃分

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

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

2.3 初始條件及邊界條件

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

模型上游邊界條件給定蘭村流量變化過程，下游邊界條件給定汾河二壩水位—流量關系。模型應用時邊界條件設置如圖4所示。

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

2.4 參數設置

曼寧系數M：利用觀測資料，通過數值驗證進行參數率定，最終選定M=35 m1/3/s。

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

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

模型計算過程中時間步長取1 s。

3 模型驗證及參數率定

3.1 橡膠壩模擬效果校核計算

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

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

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

圖5 不同壩高流量系數驗證結果Fig.5 Verification results on flow coefficient of different dam heights

3.2 河道恒定流水面線模擬計算

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

圖6 恒定流水面線計算結果Fig.6 Constant flow water surface line calculation results

3.3 洪水演進過程模擬結果驗證

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

圖7 1號橡膠壩模型驗證Fig.7 Validation results at No.1 rubber dam

圖8 6號橡膠壩模型驗證Fig.8 Validation results at No.6 rubber dam

圖9 7號橡膠壩模型驗證Fig.9 Validation results at No.7 rubber dam

4 模型應用

4.1 計算方案

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

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

4.2 方案計算結果分析

4.2.1 洪水演進過程分析

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

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

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

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

為保證洪水能在主河道內安全通過，洪峰到達時，各橡膠壩工程區內水位應低于內堤高程，對各塌壩方案下橡膠壩處最高水位線進行比較分析，繪制圖11 所示各方案最高水位包絡線圖。

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

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

4.2.2 工程區流場特性分析

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

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

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

4.2.3 方案比選分析

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

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

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

5 結 論

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

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

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