王文杰，張曉波，王 攀

（浙江省水利水電勘測設(shè)計(jì)院，浙江 杭州 310002）

1 問題的提出

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

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

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

2 汊點(diǎn)水位預(yù)測 — 校正法

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

圖1 河網(wǎng)基本計(jì)算單元示意圖

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

對于首末斷面：

（2）提取計(jì)算結(jié)果中每一河汊的進(jìn)出流量值，并計(jì)算對應(yīng)汊點(diǎn)水位校正值Δh，計(jì)算公式為：

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

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

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

（5）依次逐時段推求，直至計(jì)算結(jié)束。

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

3 平原含閘河道的處理

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

圖2 含閘河道單元示意圖

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

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

3.1 堰流狀態(tài)

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

3.1.1 自由出流

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

式中：ε為側(cè)收縮系數(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,∞）范圍內(nèi)有且僅有唯一解，因此可采用牛頓迭代法得解。已知Hl2后，分別代入式（6）～（8），并考慮Ql2= Ql3，即可解得閘上水位Zl2和閘下水位Zl3。

3.1.2 淹沒出流

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

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

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

3.2 孔流狀態(tài)

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

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

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

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

4 應(yīng)用實(shí)例

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

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

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

圖3 溫黃平原河網(wǎng)概化示意圖

圖4 溫黃平原“菲特”臺風(fēng)期間驗(yàn)證水位過程圖

5 結(jié) 語

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

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

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

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