劉恒光，郝瑞霞，何曉燕，胡孜軍

(1.太原理工大學水利科學與工程學院，山西 太原 030024； 2.中國水利水電科學研究院，北京 100038；3.水利部防洪抗旱減災工程技術研究中心，北京 100038)

橡膠壩具有美化城市景觀等功能。近年來許多城區(qū)河道建設了橡膠壩。水流在河道及河漫灘的演進過程可用水動力模型模擬[1-3]，但橡膠壩對河道行洪具有復雜的阻水影響，給含橡膠壩河道的洪水演進數(shù)值模擬帶來了困難?？紤]橡膠壩影響的水動力學模型可以有效地預測河道水流過程，模擬河網(wǎng)中復雜水利工程的調(diào)度和控制，對刻畫水位、流量等水力要素的時空變化具有重要的實用價值，能夠為相關部門提供科學精準的技術支撐，對保障流域經(jīng)濟社會平穩(wěn)發(fā)展具有重要意義[4]。

目前已有很多學者開展了閘堰影響下的水流數(shù)值模擬研究。例如：Feng等[5]建立了一種洪水演算方法，提出在閘堰調(diào)控的多分支河網(wǎng)中建立虛設單元河段，然后采用雙追趕方法求解；姜恒志等[6]采用泰勒級數(shù)展開法，進行堰出流的模擬，建立了具有內(nèi)部邊界堰存在時的河網(wǎng)非恒定流水動力學模型；施勇等[7]提出了基于水動力學的閘壩調(diào)度計算格式并應用于長江中下游防洪系統(tǒng)；張曉波等[8]利用替代法處理過閘流量并模擬了河網(wǎng)閘堰過流和調(diào)度影響；李大鳴等[9]考慮河網(wǎng)中復雜的水閘控制條件，提出了雙向迭代內(nèi)邊界控制法，對具有閘堰的特殊河段進行數(shù)值模擬；吳松柏等[10]提出了改進雙向迭代內(nèi)邊界法求解含閘堰方程，并應用于多閘聯(lián)合調(diào)度方案的數(shù)值模擬，可以有效模擬河網(wǎng)地區(qū)的閘堰過流和調(diào)度影響；劉芹等[11]全面比較了對閘堰過流能力公式線性化處理的泰勒法和替代法兩種計算方法，得出替代法具有處理方式更便捷、計算公式更簡潔、更易程序化的優(yōu)點，在河網(wǎng)模型中可優(yōu)先采用；陳煉鋼等[12]根據(jù)已知的閘壩運行水位流量資料、調(diào)度規(guī)則和設計參數(shù)等內(nèi)邊界條件給出了閘堰所在河道斷面的追趕系數(shù)，進行了閘壩調(diào)度過程模擬；楊甜甜等[13]為提高洪水預報精度建立了水文水動力耦合模型，但未考慮汛期遇洪水時研究流域內(nèi)橡膠壩等工程對河道洪水演進過程中水位和流量等水力參數(shù)的時空變化影響。可見，關于橡膠壩運行過程對河網(wǎng)水流運動影響的相關研究較少，專門考慮橡膠壩內(nèi)邊界條件的河網(wǎng)水動力學模型并不多見，需要展開深入研究。

汛期橡膠壩與普通閘堰在實際運用過程中有一定的差別，通常在大量蓄水后開始運用橡膠壩，然后橡膠壩逐漸降壩，直至完全塌壩。本文構(gòu)建考慮橡膠壩的不同過流情況的一維水動力學模型，并將該模型應用于淠河流域干流下龍爪到隱賢集河段，對近年來淠河發(fā)生的實測典型洪水過程進行模擬驗證，對橡膠壩不同運行工況進行數(shù)值模擬，并對計算結(jié)果進行對比分析。

1 控制方程與數(shù)值方法

1.1 控制方程

一維河道水動力學計算采用圣維南方程組，其連續(xù)性方程和動力方程形式如下：

(1)

(2)

式中：B為過水斷面寬度；Z為斷面水位；t為時間；Q為過水斷面平均流量；x為空間坐標；q為單寬旁側(cè)入流；α為動能校正系數(shù)；A為過水斷面面積；g為重力加速度；C為謝才系數(shù)；R為水力半徑。

1.2 控制方程的離散與數(shù)值計算方法

采用Preissmann四點隱式差分格式離散方程(1)(2)并作線性化處理，得到任一河段的差分方程為

(3)

(4)

忽略n+1，式(3)(4)可寫為

Qj+1-Qj+CjZj+1+CjZj=Dj

(5)

EjQj+GjQj+1+FjZj+1-FjZj=Φj

(6)

式中：上標n與下標j分別為時間與空間步數(shù)；Cj、Dj、Ej、Fj、Gj、Φj為計算系數(shù)，均由初值計算，具體見參考文獻[14]。

一維河網(wǎng)非恒定流水力計算的關鍵在于汊點水力要素的求解，本文采用河網(wǎng)三級聯(lián)解法[14]對河網(wǎng)模型進行求解。基本思路為：①將河網(wǎng)中兩節(jié)點間的河道劃分為若干小河段，在計算斷面上對圣維南方程組進行有限差分計算，離散后得到河道各斷面水位、流量的差分方程組；②對方程組進行自消元，形成河道首末斷面的流量與水位之間的線性關系；③根據(jù)汊點的水量平衡關系得到汊點水位方程組，求解該方程組可計算出各汊點的水位值；④根據(jù)各汊點水位值回代求解各河道斷面的水位和流量。

2 含有橡膠壩的特殊河段處理方法

在河道上布置閘壩等水工建筑物，會使河道上出現(xiàn)幾何形狀或水力特性的間斷點，河道水流不再具有連續(xù)性[14]。因此，傳統(tǒng)的圣維南方程組不再適用于此類復雜河網(wǎng)，必須根據(jù)其水力特性進行特殊處理。而汛期橡膠壩是一類具有復雜水利工程調(diào)度規(guī)則的擋水建筑物，需要合理、科學地在這類特殊河段上考慮橡膠壩的內(nèi)在調(diào)度運行方式，并保證水流運動計算的連續(xù)性和穩(wěn)定性。

2.1 橡膠壩數(shù)值模擬計算流程

水動力學模型橡膠壩數(shù)值模擬計算分為控制對象、控制條件、控制目標3個部分?？刂茖ο笫且粋€或者多個實際工程控制對象(如單級橡膠壩或梯級橡膠壩等)，汛期橡膠壩調(diào)度控制起始條件是基于實際的水利工程調(diào)度運用規(guī)則確定的，常用的規(guī)則包括壩上水位與壩頂最高高程的水位差、區(qū)間降水量與上游水庫是否超過汛限水位、水利工程出流量、關注點水位及流量、多組條件綜合等，進而達到設定的控制目標(如調(diào)節(jié)橡膠壩的塌壩速度，延時控制時間等)，最終形成完整的橡膠壩數(shù)值模擬計算流程。

當模型需要精細化模擬含有橡膠壩河段的洪水演進過程時，模型進入橡膠壩模塊，根據(jù)計算流程中的判斷條件判定是否開始塌壩，結(jié)合工程實際情況給模型設計一個合理的塌壩速度并均勻降壩，壩高在每個時間步長內(nèi)線性地降低，使橡膠壩在上游洪水到達壩址前完成降壩。確保模型的運行邏輯與實際工程安全運行下的情形相一致。模型計算出當前時刻含有橡膠壩河段的追趕方程，聯(lián)合河網(wǎng)方程組，完成求解。當前時刻計算完成后，進行下一時刻的計算，迭代直至整個計算歷時完成。

圖1為橡膠壩模塊數(shù)值模擬計算流程圖。在實際工程中，由于橡膠壩通過排出壩袋中的水來實現(xiàn)塌壩泄流，隨著壩高的降低，壩頂寬度逐漸增加，壩袋形狀由最初的實用堰型塌至近似于寬頂堰型[15]。所以在設計模型時，橡膠壩在塌壩過程中以實用堰堰型勻速進行控制塌壩，在完全塌壩后以寬頂堰的過流方式進行計算。計算完成后，結(jié)束并退出該模塊。如果當前時刻的雨水工情不滿足橡膠壩模塊的判定條件，則說明實際橡膠壩未觸發(fā)降壩，按照設計壩高的實用堰進行過流計算。

圖1 橡膠壩數(shù)值模擬計算流程

2.2 含有橡膠壩河段的計算方程

由于過壩水流為急變流，且流態(tài)變化十分復雜，模擬的關鍵技術是如何準確描述堰壩的過流能力和確定計算方程的數(shù)值計算格式，使模擬計算穩(wěn)定，不易發(fā)生震蕩[12]。一維水動力學模擬中堰壩調(diào)控的線性化處理方法主要有泰勒級數(shù)展開取一階項[6]、前時段計算結(jié)果近似替代(簡稱替代法)[8]、轉(zhuǎn)化微分形式[12]、能量方程構(gòu)建調(diào)度計算河段方程[7]以及雙向迭代法[9]等。由于替代法的方法處理方式簡潔、易于程序化實現(xiàn)，本文采用替代法對含有橡膠壩的特殊河段進行處理。

當流態(tài)為自由出流和淹沒出流[14]時，過流流量分別為

(7)

式中：Q為過壩流量，m3/s；m為自由出流系數(shù)；B為壩寬，m；Zj、Zj+1分別為當前時刻壩上游水位和下游水位，m；Zd為壩頂高程，m；φ為淹沒出流系數(shù)。

在模型數(shù)值模擬計算中，水流流態(tài)處于自由出流與淹沒出流過渡段時，可能會出現(xiàn)數(shù)值異常，產(chǎn)生數(shù)值振蕩或發(fā)散等問題。采用線性插值連接自由出流和淹沒出流這兩種流態(tài)過渡區(qū)間段的方法，可以保證流態(tài)轉(zhuǎn)換過程中的連續(xù)性。經(jīng)過本文數(shù)值模擬驗證，采用這種方法，既可以有效提高模型的穩(wěn)定性，又可以保證模型的精確性。

根據(jù)替代法思想[8]，在計算時間步長不大的前提下，認為前后兩個時段斷面的水力要素變化較小，可以利用前一時刻的部分變量替代當前時刻的部分變量，以達到分離未知量的目的，這樣的處理方法在文獻[8,11]都有應用和驗證。

以河道首末斷面水位為未知量構(gòu)建含有橡膠壩計算河段的雙追趕方程[17]：

Qj=αj+βjZj+ζjZj+1

(8)

Qj+1=θj+1+ηj+1Zj+1+γjZj

(9)

式中：Qj、Qj+1分別為當前時段含有橡膠壩河段j斷面與j+1斷面的流量，m3/s；αj、βj、ζj、θj、ηj、γj為式(5)(6)相應的追趕系數(shù)。

本文只研究壩的過流在水力計算中的特殊處理，故將過壩流態(tài)分為自由出流和淹沒出流兩大類。推求替代法構(gòu)建的含有橡膠壩計算河段的追趕系數(shù)αj、βj、ζj、θj、ηj、γj，具體計算公式如表1所示。

表1 替代法推導的含有橡膠壩河段的追趕系數(shù)計算公式

3 實例驗證

3.1 流域概況及模型建立

為檢驗受汛期橡膠壩調(diào)度影響所建立的水動力學模型的穩(wěn)定性和有效性，選取位于淮河支流淠河流域佛子嶺、響洪甸兩座水庫至迎河集區(qū)間作為試驗河段，重點研究位于下龍爪站至隱賢集站河段六安城區(qū)的六安橡膠壩和六安城北橡膠壩兩座串聯(lián)橡膠壩，如圖2所示。

圖2 淠河流域研究區(qū)域

淠河是淮河中游南岸的一條較大支流，全長260 km，流域面積6 000 km2，是六安市防汛重點河流。六安城區(qū)兩座橡膠壩對流域防洪安全十分重要。將研究流域整個河道共分為現(xiàn)東淠河河段、西淠河河段、淠河干流上游河段、月亮島東側(cè)河段、月亮島西側(cè)河段和淠河干流下游河段等6個河段，其中橫排頭、下龍爪、隱賢集3站具有實際觀測資料。本文采用斷面資料為2019年實測地形數(shù)據(jù)，其中東淠河上有42個斷面，西淠河上有31個斷面，淠河干流上有86個斷面，共159個大斷面，斷面平均間距約1 km，滿足模型計算要求。采用近年來發(fā)生的典型歷史洪水過程進行演算作為驗證，時間步長Δt取200 s，空間步長采用斷面間距，由于河道區(qū)分灘地和河槽不同地形條件，糙率n的取值范圍為0.035～0.050。以響洪甸水庫和佛子嶺水庫的實測調(diào)度出庫流量作為上邊界條件，迎河集站水位流量關系作為下邊界條件，在此基礎上建立水動力學計算模型進行洪水演算。為有效模擬汛期六安兩座橡膠壩受到的影響，在壩上和壩下分別設置上下兩個特殊節(jié)點作為橡膠壩模塊的特殊計算河段。上游六安橡膠壩充滿狀態(tài)時的壩頂高程P1為36.00 m，下游六安城北橡膠壩充滿狀態(tài)時的壩頂高程P2為34.40 m，兩座橡膠壩概化圖見圖3所示。

圖3 淠河串聯(lián)橡膠壩概化(單位：m)

根據(jù)《六安市淠河城區(qū)段橡膠壩群聯(lián)合運用調(diào)度辦法》調(diào)度運用規(guī)則條件及要求，當壩上水位超壩頂最高高程0.3 m，橫排頭下泄流量達到200 m3/s，橫排頭至橡膠壩區(qū)間降水量達到50 mm，且上游水庫超汛限水位，并有大暴雨時串聯(lián)橡膠壩群進行降壩運用。橡膠壩群自下而上開始逐級泄水，即在上游洪水到來之前，首先泄掉六安城北橡膠壩的蓄水，然后開始泄六安橡膠壩的蓄水。由于上游水庫放水到達六安橡膠壩處時間約為2 h，洪水到達壩前時，橡膠壩均已完成降壩。根據(jù)兩座橡膠壩的實際壩高，計算出模型中兩座串聯(lián)橡膠壩在每個時間步長內(nèi)的降壩高度達到0.23 m，即降壩速度v至少為0.23/Δt，方可達到實際工程要求。所以設計降壩過程中每個時間步長橡膠壩的壩高會隨著時間變化，直至上游洪水到達橡膠壩前，兩座橡膠壩已經(jīng)按照調(diào)度運用規(guī)則完成塌壩，則兩座橡膠壩均按照完全塌壩高度的寬頂堰堰型進行過流計算。綜上所述，本模型中的控制對象是六安城區(qū)的兩座串聯(lián)橡膠壩?？刂茥l件為同時滿足：①橡膠壩群壩上水位大于0.30 m；②橫排頭下泄流量大于200 m3/s；③橫排頭至橡膠壩區(qū)間降水量大于等于50 mm時。開始執(zhí)行淠河串聯(lián)橡膠壩對象組合控制目標，即下游的六安城北橡膠壩開始執(zhí)行以0.23/Δt的啟閉速度勻速降壩，最終控制目標為完全塌落，完全塌落后，按寬頂堰堰型過流計算。上游的六安橡膠壩在達到控制條件時，延遲相應的時間(本文設為3 400 s)再執(zhí)行塌壩控制，以0.23/Δt的啟閉速度勻速降壩，最終控制目標為完全塌落。直至上游橡膠壩降至壩底高程后，也按照寬頂堰堰型過流計算，直至整個計算歷時完成。

3.2 情景模擬及結(jié)果分析

下龍爪和隱賢集水位站分別位于兩座橡膠壩的上游和下游。以2020年淠河流域發(fā)生的典型洪水“20200719”為例進行情景模擬，以下龍爪和隱賢集兩個水位站的實測資料作為水動力學模型的驗證資料，對洪水期間未塌壩運行(情景Ⅰ)和考慮橡膠壩影響(情景Ⅱ)兩種情景分別進行數(shù)值模擬，通過比較典型洪水不同情景下的數(shù)值模擬結(jié)果與兩個水位站實測資料的擬合程度，進行數(shù)據(jù)分析。計算結(jié)果對比過程線如圖4所示，橡膠壩所在斷面的模型計算流量和水位過程如圖5所示。

圖4 下龍爪站和隱賢集站“20200719”洪水水位模擬過程對比

圖5 “20200719”洪水橡膠壩所在斷面的模型計算流量和水位過程

模型在兩種情景下的計算過程中均沒有出現(xiàn)不收斂的情況，說明本文構(gòu)建的水動力學模型算法穩(wěn)定。情景Ⅱ的計算結(jié)果與下龍爪、隱賢集水位站的實測數(shù)據(jù)擬合程度較高，漲落趨勢一致，洪峰水位及峰現(xiàn)時間基本吻合。相較于情景Ⅰ，情景Ⅱ的計算精度得到明顯改善。由圖4(a)可知，下龍爪站在2020年7月19日22：36水位達到峰值6.94 m，情景Ⅰ計算水位在2020年7月20日0：54出現(xiàn)洪峰水位7.89 m，絕對水位誤差為0.95 m，情景Ⅱ計算水位在2020年7月20日0：43出現(xiàn)洪峰水位6.71 m，絕對水位誤差為0.23 m；由圖4(b)可知，隱賢集站實測水位在2020年7月20日7：18水位達到峰值10.65 m，情景Ⅰ計算水位于2020年7月20日11：33出現(xiàn)洪峰水位9.60 m，絕對水位誤差為1.05 m，情景Ⅱ計算水位于2020年7月20日9：47出現(xiàn)洪峰水位10.40 m，絕對水位誤差為0.25 m。影響水位精度的主要因素可能是：①淠河下游控制斷面的邊界條件迎河集水位流量關系不穩(wěn)定且呈現(xiàn)復雜的繩套關系；②淮河干流的頂托現(xiàn)象對淮河支流淠河下游河道可能造成實測水位比模擬水位偏高，導致計算結(jié)果與實測結(jié)果有一定誤差，但模擬結(jié)果整體滿足精度要求，具有較高的穩(wěn)定性和可靠性。

另外選取近年來發(fā)生的“20160702”“20180818”和“20200623”3場典型歷史洪水進行分情景演算，并用實測水位資料、峰現(xiàn)時間誤差進行驗證，情景Ⅱ的計算擬合效果均優(yōu)于情景Ⅰ的擬合效果，將4場典型洪水計算結(jié)果整理至表2。由表2可知，情景Ⅱ?qū)崪y水位與水位計算值相差極小，峰現(xiàn)時間也吻合較好，表明根據(jù)本文提出的思路所建立的一維水動力學模型可以較好地模擬含有橡膠壩河段的洪水演進過程。而在實際模擬中，如果對洪水期間橡膠壩塌壩過程考慮較少，導致計算誤差較大，對洪水預報精度產(chǎn)生影響。

表2 歷史典型洪水計算誤差

4 結(jié) 論

a.在Preissmann四點隱式差分格式離散方法基礎上，提出了以橡膠壩作為內(nèi)邊界條件的河網(wǎng)水力計算模型，主要包括含有橡膠壩的特殊河段計算流程和考慮橡膠壩工程過流計算中采用的處理方法等。

b.針對淠河流域下游橡膠壩的運行規(guī)則特點，構(gòu)建淠河流域橡膠壩數(shù)值模擬計算流程，通過多年實測水文資料驗證，對比洪水期間兩種工況下的情景模擬結(jié)果可以看出，過于簡化考慮橡膠壩運行過程會使計算結(jié)果產(chǎn)生較大誤差。而考慮橡膠壩影響的水動力學模型可以提高洪水期間的模擬精度，有效模擬含有橡膠壩河段的洪水演進過程。

c.本文采用的模型可以較好地擬合淠河流域的洪水演進過程，可以穩(wěn)定、精確地模擬六安城區(qū)含橡膠壩河段的洪水演進過程，并在淠河流域洪水預報預警調(diào)度系統(tǒng)中取得較好的應用效果?？蓪⒃撃Ｐ屯茝V應用于其他含有橡膠壩河段的洪水預報。