基于汊點水位預測校正的閘控河網水動力模擬

王文杰，張曉波，王 攀

（浙江省水利水電勘測設計院，浙江 杭州 310002）

1 問題的提出

平原河網水動力數(shù)值模擬廣泛應用于平原地區(qū)防洪排澇、生態(tài)調度、水環(huán)境整治等諸多水利工作，是水利規(guī)劃設計的重要研究手段。河網水動力數(shù)值模擬的原理是求解離散化的圣維南方程組，但由于系數(shù)矩陣的稀疏非帶狀特征，直接求解河網整體離散方程組往往耗時較長，已無法滿足智慧水利等相關工作對模擬速度、精度的高要求。因此，目前主流思路是對河網進行分級求解，如通過特殊的河道（汊點）編碼，把河網系數(shù)矩陣轉換為類似單一河道的帶狀矩陣進行求解[1-3]，或先將河道首末斷面變量與汊點連接條件及邊界條件聯(lián)合建立整體矩陣，得解后再反推河道內各斷面水力要素值[4-5]。但無論是哪種方法，在斷面數(shù)較多、初始誤差較大時都存在誤差累計效應，若首末斷面變量關系存在失真，也極易導致求解計算崩潰。為此，近年來有學者提出汊點水位預測 — 校正法[6]，主要思路是利用河網汊點回水效應，建立汊點處水位校正量與流量的關系，先對未知時層的汊點水位進行迭代收斂，求得汊點水位后再反推汊點間內河道各斷面水力要素值，由于無需建立和求解整理連接矩陣，因此能有效避免分級解法中出現(xiàn)的各種不穩(wěn)定現(xiàn)象。圍繞該方法，不少學者又相繼開展水流 — 水質耦合模型、一二維河湖系統(tǒng)耦合模型等研究，豐富了該方法的應用場景[7-8]。

但是，當河網內建有水閘、泵站等控制性泄水建筑物時，由于含閘河道水流連續(xù)狀態(tài)被破壞，通過求解非線性離散方程回推進汊流量的方法失效，導致有關汊點水位無法像普通河道一樣直接迭代校正。針對該問題，朱德軍等提出2點改進[9]：一是建筑物局部為緩流時，可近似處理成局部水頭損失；二是將建筑物視為連接上下游河道的特殊汊點，與普通汊點一起參與計算，并在迭代步中不斷修正建筑物下游水位直至滿足流量守恒條件。但在實際應用中，上述改進均存在一定局限，如局部水頭損失的量化十分復雜，且計算精度無法保證；作為特殊汊點參與計算則僅適用于局部流態(tài)為急流的情況，且需要事先明確建筑物上游側的邊界條件。由于平原河網水位相差有限，泄水建筑物局部多以堰流為主，其過流能力通常與上、下游水位相關，有的甚至沒有固定的水流流向，以上改進均難以應用于存在較多閘站工程的平原河網水動力模擬中，限制了推廣范圍。

針對以上缺陷，本文以平原河網最為常見的平底式閘門為例，引入閘門不同流態(tài)過流流量公式作為建筑物與上下游河道的連接條件，通過建立上下游河道與水閘過流條件聯(lián)合方程并求解，可有效解決原方法中含泄水建筑物河道進汊流量無法推求的問題，且對泄水建筑物的控制形式可根據(jù)實際調度需要靈活設置，使汊點水位預測校正法能夠有效地應用于含較多控制性閘站的平原緩流河網水動力模擬中，拓展該方法的應用范圍。本文以浙江省溫黃平原河網為例，驗證該方法的可行性和穩(wěn)定性，分析模型的計算效率。

2 汊點水位預測 — 校正法

樹狀或環(huán)狀河網的基本單元可視為由汊點和河段組成（見圖1）的結構，其中汊點表現(xiàn)為不同河段的交接點或邊界條件。當河網水動力過程由一維圣維南方程組描述時，對于包含n個斷面的河段，每個未知時層均可聯(lián)立2（n -1）個方程式，但未知水力要素有2n個（Z1，Q1，…Zn，Qn），因此若給定兩端汊點處未知時層水位，則齊次方程組成立，河段未知時層全部水力要素得解。汊點水位預測校正法的思路便是通過不斷的預測 — 校正 — 迭代過程，優(yōu)先求解汊點處未知時層水位，再反推至河段內各斷面。其求解流程如下（假設T時刻全部水力要素已知，欲求T+1時刻要素值）。

圖1 河網基本計算單元示意圖

（1）以T時刻水位作為各汊點的初始迭代水位，代入河段方程組中求得各斷面該預測水位下的水力要素值。朱德軍等提出可采用Newton - Raphson 法求解形成的非線性方程組[5]。本文建議在河段不包含控制閘站時，采用追趕系數(shù)法進行求解，即：

對于首末斷面：

（2）提取計算結果中每一河汊的進出流量值，并計算對應汊點水位校正值Δh，計算公式為：

式中：Qi、Ai、Bi分別為流入汊點斷面的流量（m3/s）、過流面積（m2）和面寬（m）；Qo、Ao、Bo分別為流出汊點斷面的流量、過流面積和面寬；α為近似系數(shù)，一般取1.0 ～2.0[4]。可以證明Δh具有穩(wěn)定的收斂性。

（3）若Δh小于預設容差（一般取0.005 ～ 0.010 m），則認為當前汊點滿足汊點連接條件，水位無需修正，否則繼續(xù)校正汊點水位：+Δh，并以作為下一迭代步的汊點水位代入河段方程組進行計算。

（4）若某迭代步中所有汊點均能滿足汊點連接條件，則認為該狀態(tài)下的汊點水位組合符合實際，并以該狀態(tài)下各斷面的水力要素值作為T+1時刻的要素值。

（5）依次逐時段推求，直至計算結束。

可以看出，該方法的求解思路類似于分級解法，即通過分步求解河汊和河道2部分的水力要素值，避免面對大型稀疏系數(shù)矩陣，從而降低方程組求解的復雜性，提高模型計算速度。但區(qū)別于傳統(tǒng)分級解法，該方法在構建模型時不用區(qū)分干支流、無需優(yōu)化汊點編碼結構，因此具有較強的靈活性和易于模塊化編程的特點，且在每個迭代步推求河汊進出流量時，有條件實現(xiàn)多單元并行計算，進一步提升模型響應速度，在大規(guī)模河網模擬中具有較好的應用前景。根據(jù)對溫黃平原實例的測試驗證，相比傳統(tǒng)三級解法，該方法在同等計算條件下完成時長可減少34%，效率提升較為顯著。

3 平原含閘河道的處理

當河段包含閘（泵）等建筑物時，雖然閘室局部水流流態(tài)復雜，但建筑物前后河道仍滿足緩流條件。若設想一個包含閘（泵）建筑物的河道單元，由建筑物、上游河段、下游河段及兩端汊點組成，并記2河段首末斷面分別為l1、l2、l3、l4（見圖2），則根據(jù)式（1）及式（2），可建立關系：

圖2 含閘河道單元示意圖

當采用汊點水位預測校正法時，上述方程僅汊點處水位Zl1、Zl4已知，因此無法直接求解河段各斷面水位、流量值。但是，建筑物處通常滿足流量連續(xù)條件，即Ql1= Ql3，且根據(jù)閘門過流經驗公式，該流量一般是閘門上下游水位的單值或多值函數(shù)，即滿足Ql1= Ql3= f閘( Zl2)或Ql2= Ql3=f閘（Zl2，Zl3）。若將該連接條件與式（6）、式（7）聯(lián)立，則可構建河段未知量的齊次方程組。

本文以平原河網最為常見的平底閘門為例，分別引入堰流和孔流2種不同流態(tài)過流經驗公式作為建筑物處的連接條件，基本涵蓋平原水閘的常見調度場景，使得在已知汊點水位條件下，含閘河道單元內各項水力要素可解。

3.1 堰流狀態(tài)

對于平原平底閘門，當閘門啟出水面，不影響閘壩泄流量時的水流狀態(tài)為堰流，一般以閘門開啟高度與堰上水頭的比值大于0.65作為判別條件。在堰流狀態(tài)下，根據(jù)閘上下游水位關系又可分為自由出流和淹沒出流2種流態(tài)，其中自由出流時閘下水流對泄流能力影響較小，淹沒出流需要考慮閘下水流的頂托影響，因此2種流態(tài)的求解方法存在顯著區(qū)別。

3.1.1 自由出流

根據(jù)SL 265 — 2016《水閘設計規(guī)范》[11]，當堰流流態(tài)為自由出流時，過閘流量僅與閘上水頭有關，對于平原緩流河網，若忽略閘上行近流速，其過流能力公式為（以圖2為例，下同）：

式中：ε為側收縮系數(shù)；m為流量系數(shù)；b為水閘總凈寬（m）；d為閘底板高程（m）；g為重力加速度，可采用9.81 m/s2。

令Hl2= Zl2-d，并將式（8）與式（6）聯(lián)立，可以得到Hl2的一元函數(shù)。整理后函數(shù)形式為：

由于追趕系數(shù)ηl2＜0恒成立且Hl2＞0，因此可證函數(shù)系數(shù)b＞0、c＞0，則式（9）在Hl2∈（0,∞）范圍內有且僅有唯一解，因此可采用牛頓迭代法得解。已知Hl2后，分別代入式（6）～（8），并考慮Ql2= Ql3，即可解得閘上水位Zl2和閘下水位Zl3。

3.1.2 淹沒出流

根據(jù)SL 265 — 2016《水閘設計規(guī)范》，一般當上下游水頭比（Zl3-d）/（Zl2-d）≥0.9時認為水閘流態(tài)為淹沒出流，此時過流流量與閘上下游水位均有關系。平原河網水級相差有限，因此淹沒出流在閘站運行工況中較為常見?！端l設計規(guī)范》中推薦平底閘門淹沒出流流量采用下式計算：

式中：μ′為淹沒出流綜合流量系數(shù)，可采用下式計算，當計算時間步長較小時，前后時刻水位變化有限，Zl2、Zl3可取前一時刻計算成果代替；其余符號意義同前：

已知追趕系數(shù)ηl2＜0、βl3＞0，因此有m-m′＜0，又根據(jù)式（6）和式（7）可以證得n-n′＜0，則式（15）在（Zl2-Zl3）∈（0,∞）范圍內有且僅有唯一解。采用牛頓迭代法求得（Zl2-Zl3）后，代入式（13）即可分別求得閘上下游水位值，回代河段追趕方程則單元內所有未知水力要素得解。

3.2 孔流狀態(tài)

當閘門門板對過閘水流有不可忽視的阻擋影響時的水流狀態(tài)為孔流，可采用閘門開啟高度與堰上水頭的比值小于0.65作為判別條件。根據(jù)《水閘設計規(guī)范》，當忽略閘上行近流速時，孔流狀態(tài)下自由出流和淹沒出流流態(tài)時的過閘流量均可采用下式計算：

式中：σ為孔流淹沒系數(shù)，μ為孔流流量系數(shù)，均可查表或根據(jù)經驗公式確定；he為孔高（m）。各參數(shù)單位與堰流自由出流時類似，令Hl2= Zl2-d，并將式（16）與式（6）聯(lián)立，可得到Hl2的一元函數(shù)。整理后函數(shù)形式為：

同理，可證式（17）在Hl2∈（0,∞）范圍內有且僅有唯一解，因此可采用牛頓迭代法得解。已知Hl2后，分別代入式（6）、式（7）和式（16），并考慮Ql2= Ql3，即可解得閘上水位Zl2和閘下水位Zl3。

此外，閘門關閉或啟用泵站是已知建筑物處流量的特殊情況。當閘門關閉時，有Ql2= Ql3= 0；當啟用泵站時，有Ql2= Ql3= q，其中q為泵站當前時刻提水流量。將上述關系代入式（6）和式（7）則方程得解。采用以上方法基本可處理平原平底式水閘運行中的各類工況，在實際模擬中，可先根據(jù)上一時刻閘門上下游水位或調度指令判斷當前時刻過流狀態(tài)，再根據(jù)堰流或孔流選取對應的推求方程進行求解。此外，由于堰流淹沒出流和孔流的流量系數(shù)公式中均采用上一時刻水位作為本時刻水位的替代，因此需注意替代成立條件，時間步長的選取不宜過大，避免在流量系數(shù)中引入較大誤差。

4 應用實例

為了進一步分析論證以上方法的可行性，將上述方法用于浙江省溫黃平原河網水動力數(shù)值模擬。溫黃平原位于浙江省臺州市，北臨椒江、東部南部面海，域內人口經濟集聚，是浙江省重點治澇區(qū)域。為分隔水級和抵御外潮，在其河網內部和排海出口處建有眾多節(jié)制閘站，部分閘站處規(guī)劃新建排澇泵站，是較為理想的分析對象。2013年“菲特”臺風期間，平原發(fā)生嚴重內澇，是近年來發(fā)生的較為典型的洪澇災害，且站點實測資料完整，因此選取該場次洪水過程對模擬結果進行驗證。

模型共概化有313個汊點、1 370個斷面、42個邊界，模擬100座節(jié)制閘站，由模型根據(jù)其上下游水位關系和設定的閘站控制水位自主調度。河道糙率和計算時間步長均在正常范圍內取值。模型河網網絡概化見圖3。計算所得各代表斷面水位過程與實測結果對比見圖4。

其中，巖頭閘為外海擋潮閘，閘下即為潮位邊界，需根據(jù)潮位漲落控制開關調度，金清閘為內河節(jié)制閘，洪水期間閘門全開泄洪。由計算結果可以看出，2種閘門類型計算水位過程與實測過程均擬合良好，最高水位差值在0.02 m以下，說明模型對閘門過流過程的模擬符合實際。桐嶼站、溫嶺站均為平原內部控制站點，由結果可見，模型在不同洪水階段、不同來水條件下均能獲得良好的模擬效果，體現(xiàn)了上述方法的穩(wěn)定性。

圖3 溫黃平原河網概化示意圖

圖4 溫黃平原“菲特”臺風期間驗證水位過程圖

5 結 語

（1）汊點水位預測校正法在確保模擬精度的前提下，具有通用性好，計算效率高，易于編程實現(xiàn)等突出優(yōu)勢，在大型復雜河網水動力數(shù)值模擬中有較強的應用價值，但原方法及現(xiàn)有改進對含泄水建筑物河網的模擬仍存在不足，難以實際應用。

（2）引入閘門過流能力公式作為建筑物與前后河道的連接條件，通過聯(lián)立閘門過流能力方程和前后河道的追趕方程，并將方程組轉化為關于閘前后水力要素的一元函數(shù)，可解決原方法中含閘河段進汊流量無法直接內推的問題，使原方法在大型復雜閘控河網中的應用成為可能，拓展了該方法的應用場景。

（3）以溫黃平原“菲特”臺風期間洪水過程為例進行實例驗證。結果表明，該改進可良好模擬平原不同調度規(guī)則下的泄水建筑物，具備可行性和穩(wěn)定性，且模擬誤差滿足水動力模擬精度要求。

（4）后續(xù)可結合該方法汊點結構和計算特點進行并行計算開發(fā)，進一步提升模型的響應速度，同時優(yōu)化流量系數(shù)的確定方式，減少替代誤差對計算結果的影響，為開展平原地區(qū)智慧水利建設、閘站聯(lián)合優(yōu)化調度等復雜調度計算奠定基礎。