灌溉渠系非恒定流數值模擬及閘門運行設計

古 玉，方天鈺，陳 建，劉曉東，韓 宇

(1.中國農業大學水利與土木工程學院，北京 100083；2.中國農業大學工學院，北京 100083)

0 引 言

目前我國水資源短缺問題十分嚴重，已成為制約我國國民經濟和社會發展的重要因素之一。我國現有灌溉面積0.664 億hm2，農業用水更是占到用水總量的62.4%[1]，因此為應對水資源危機，采用農業用水節水措施是必要的。灌溉渠系進行輸配水時，需要根據灌區用水量要求通過閘門進行調節，此時渠道中的水流受到干擾必定會出現非恒定流，使渠道中的水位、流量既沿程發生變化，又隨時間發生變化。由于無法預測該過渡過程，常發生重復調度，導致渠道配水的可靠性及準確性很差，甚至跑水和棄水，造成水量浪費，降低灌溉質量。為使整個渠系引水、輸水、配水過程科學合理，減少渠道輸水損失，提高水的利用效率，達到節水的目的，數值模擬是最為有效可行的方法。

明渠非恒定流的理論研究起源于18世紀后期Laplace和Lagrange對淺水波的研究，之后Lagrange給出了其波速公式。到1871年，Saint-Venant在做了大量研究的基礎上，提出了明渠非恒定流的理論和通用方程，即Saint-Venant方程，使明渠非恒定流現象得到了數學描述，推動了通過數學方法研究這一問題的發展。之后國內外學者不斷深入研究，提出了有限差分法、特征線法和有限單元法等多種求解方法，其中有限差分法又包括直接解法、分級解法、單元劃分法、混合解法以及松弛迭代法等。相關理論的不斷發展完善的基礎上，隨著計算機的應用，明渠非恒定流的數值模擬也得到了充分發展。法國的CARIMA[2]模型以及荷蘭的MODIS[3]、DUFLOW[4]模型都在田間試驗中得到了較好驗證。Gary P Merkley[5]教授開發出了一個可視化的模擬軟件CanalMan，它允許用戶自己建立所需渠道并定義相關參數，進一步修改后形成RootCanal水力模擬軟件，但在后來的使用中發現此軟件也存在著一些問題。到現在為止，應用較為廣泛且效果較好的一個軟件是丹麥水力研究所(DHI)開發的Mike11。Mike11雖然模擬效果較好，但其無法模擬閘門開度調控過程，且費用很高。在國內，20世紀80年代起，我國的非恒定流數值模擬開始有了進展，針對灌溉渠網，從1982年到2000年，王長德[6-9]建立了上、下游常水位自動控制渠道以閘門運動為動態邊界條件的渠道非恒定流計算模型，并將P+PR算法(比例+比例微分算法)與比威爾算法相結合。2001年，趙競成等[10]建立了較完整的渠系運行模型，編制了具有一定通用性和可擴充性的計算機模擬軟件。到2010年，韓宇等[11]模擬了上游和下游常水深兩種運行方式下梯形斷面渠段由于不同閘門調控方式所引起的非恒定流過渡過程并進行了比較，相對明確地模擬了閘門調節過程。2013年，章少輝等[12]把渠系輸配水控制方程變形處理，在具備高精度和分辨率的標量有限體積法的基礎上建立了末級渠系非恒定流數值模擬模型并驗證。

發展至今，非恒定流模擬的相關理論已經較為完善，Preissmann四點隱格式和特征線最為常用。在數值模擬方面，河網非恒定流數值模擬技術已經較好，Mike軟件得到普遍使用且模擬效果較好，但灌溉渠網中非恒定流數值模擬仍缺乏較好的相關模型或程序。且國內外相較來說，我國針對非恒定流的研究開始的較晚，雖然理論方面發展較完善，但多基于南水北調工程，忽略某一段面的改變對其上下游穩態水流的干擾，并且由于灌溉渠道中閘門的調控方式及其形成的邊界條件不容易合理概化，針對河網的數值模擬模型無法直接用于灌溉渠道，所以在灌溉渠網非恒定流數值模擬技術方面仍有所欠缺。另外，對于數值模擬可視化的研究較少，且多針對河網，2000年，徐小明等[13]用Fortran語言，對計算出的水位變化過程作了可視化顯示研究；2002年，馬洪明等[14]，運用計算機信息平臺，建立了一套圖形顯示水力模型的支持系統；2008年，陳棟[15]基于松弛迭代法，運用混合編程技術，開發了河網非恒定流數值模擬的可視化系統；但這些成果都不夠成熟故沒有得到推廣使用。隨著計算機技術的不斷發展成熟，將計算機技術與非恒定流模擬相結合，建立起一個較為簡便且適用性強的灌溉渠系非恒定流模擬軟件，對渠道各處的水力要素進行模擬，將能夠使工作人員及時了解渠道水流的動態變化過程，給出適當的閘門調控方法，實現對水資源的科學合理調配，從而提高渠系輸配水的利用率，大大減少農業用水的浪費。

本文采用特征線法求解Saint-Venant方程，并建立簡單渠系模型模擬引湯灌區北干渠及九、十支渠系閘后水流過渡過程，并與Mike11模擬結果進行比對并驗證該模型準確性。其后基于該模擬結果利用閘下出流公式建立閘門開度模擬模型，同時對多級閘門同步與逐級調控兩種情況進行模擬，討論其優缺點，對渠系閘門調控進行了簡單設計。

1 灌溉渠系非恒定流數學模型

1.1 控制方程

非恒定流是指過水斷面上的水力要素(流量、水位及流速等)不僅隨空間位置的改變而變化，而且隨時間不斷變化的水流。灌溉渠道中的非恒定流為明渠非恒定流，圣維南方程就是針對這一問題的函數關系式，即明渠非恒定流的控制方程，它由連續性方程和運動方程組成。

從有旁側入流的明渠非恒定流中取一長為ds的微小流段，由質量守恒原理可推出連續性方程如下：

(1)

假定水面波動是漸變的，忽略局部水頭損失，由動量定理或牛頓第二定律可推導出其運動方程如下：

(2)

式中：B為水面寬度，m；Z為水位，m；Q為流量，m3/s；s為流程，m；ql為單位長度上的側向匯流流量[流入為正，流出為負，m3/(s·m)]；A為過水斷面面積，m2；g為重力加速度，m/s2；C為謝才系數；R為水力半徑；h為水深，m。

1.2 求解方法

本文采用特征線法，用矩形網格劃分求解域從而把微分方程的連續問題離散化，以有限個網格點代替原連續區域，用數值積分逼近微分方程，建立網格函數的代數方程組，在離散點上求解網格函數的值，總的來說就是把圣維南偏微分方程組化為常微分方程即特征方程組進行求解[16, 17]：

(3)

(4)

采用矩形網格特征差分法進行離散：首先將求解域依據時間t和距離s進行矩形網格劃分，如圖1，網格的交點稱為結點。其中，i=0,1,2,…,n表示距離s的段號；j=0,1,2,…,n表示時間t的段號。Δsi=si-si-1表示距離步長，即渠段分段長度，可以等步長也可變步長；Δtj=tj-tj-1表示時間步長，需要根據穩定格式的限制和實際要求確定。圖1中，順、逆特征線(ω±)在時層j+1上的交點為點P，與時層j的交點分別為L和R，其間的結點為M，左右兩邊點各為點D和點E，求解方程即求解交點P。當給定初始條件后，則根據時間步長進行循環計算求解。

圖1 矩形網格圖示意圖Fig.1 The Schematic of Rectangle grid diagram

要求點P的水位和流量，則需先求出點L和R，這兩點不在網格結點上，因此采用柯朗格式進行線性內插，之后利用積分和積分中值定理可得到離散結果如下：

(5)

(6)

(7)

同時，由于柯朗格式為一階顯格式，為使計算穩定，需要選擇合理的距離步長和時間步長使其滿足穩定性條件式(8)：

(8)

1.3 邊界條件

灌溉渠道邊界條件包括內、外兩種邊界條件，其中內邊界條件是指閘門等控制建筑物在渠段中的調控所產生的邊界處的變化，外邊界條件是指入流及出流邊界條件。

上述的計算格式只適用于內部點的計算，如圖2所示，對于邊界點來說：僅有一條特征線，即只有兩個方程，但其未知數有3個，故需補充邊界條件。

上、下游邊界條件一般情況下有如下3種情況：①Z=Z(t)，即水位過程已知；②Q=Q(t)，即流量過程已知；③Q=f(Z)，即水位流量關系已知。

圖2 上下游邊界點示意圖Fig.2 The schematic diagram of the upstream and downstream boundary points

1.4 汊點處條件

灌溉渠系多呈樹狀，因此在其布置形式中必然存在眾多分岔。灌溉渠道分為干、支、斗、農四級，它們之間必然由一系列分支、匯交點連接，這些連接點即為汊點。在汊點處，水流需滿足相應的連接條件，即流量相容和能量相容條件[18, 19]。

圖3 分岔渠道示意圖Fig.3 The Schematic diagram of the distribution channel

流量相容條件為：

(9)

式中：n為汊點處的渠道數目；Qi為第i個渠道在汊點處的流量；當渠道水流流入汊點時sli取“+1”，流出時取“-1”。對圖3所示的分岔渠道，其流量相容條件可寫為：

Q1=Q2+Q3

(10)

能量相容條件為：

首先將汊點概化為一個幾何點，如圖3，斷面A-A、B-B、C-C圍出了一個小流域，由于只須計算其邊界斷面上的流量或水位，故可不研究此內的非恒定流狀態，假設流域的面積為零。在理想情況下，即不考慮局部水頭損失，同時一般情況下不存在水位突變，故有：

Z1=Z2=Z3

(11)

式中：Zi(i=1,2,3)為第i個渠道在汊點處的水位。對于灌溉渠道，汊點處經常安有閘門等水工建筑物，故考慮局部水頭損失，可列出能量方程如下：

(12)

(13)

式中：αi(i=1,2,3)為能量修正系數且α1=0。

結合上述汊點處條件并參考文獻[17],本文對分岔渠系的模擬采用逐級追趕的方法。顧名思義，就是對渠系中的各級渠道逐次求解，以前一級渠道模擬得到結果作為下一級渠道的初始條件，分岔處則根據汊點處的流量相容條件，用干渠所得結果減去支渠的結果作為下一級渠道的初始條件，這樣逐級追趕，得出所需結果。這種方法思路清晰，易于理解與掌握，在編程過程中只需注重汊點處的條件變化設置，可以大大減少編程工作量。在實際應用中，若汊點較少且流量變化率不大，則可直接按照此法逐級求解。

2 渠系數據及模擬條件

引湯灌區位于黑龍江省湯原縣，灌區的水源來自湯旺河，設計灌溉區域2.68 萬hm2，工程總的設計引流量為160 m3/s。本文選取部分引湯北干渠及其九、十支渠作為分岔渠系模擬模型，相關模型數據見表1、表2及表3。

表1 北干渠數據表Tab.1 North trunk channel data table

表2 九支渠數據表Tab.2 Data table for nine canals

表3 十支渠數據表Tab.3 Data table for ten canals

模擬條件：模擬總時長200 min，時間步長取值根據式(8)計算得到：干渠1，Δt≤62 s；干渠2，Δt≤61 s；干渠3，Δt≤62 s；九支渠，Δt≤108 s；十支渠，Δt≤112 s，故各渠道模擬時間步長均取為60 s。各渠道均采用下游常水深運行方式，其中干渠1上游流量在20 min內從25.28 m3/s線性增加到30.11 m3/s，九支渠上游流量在20 min內從1.02 m3/s線性增加到1.32 m3/s，干渠2的上游斷面的流量變化則為干渠1模擬得到的下游流量減去九支渠上游斷面的流量。十支渠上游流量在20 min內從0.79 m3/s線性增加到1.02 m3/s，干渠3上游斷面的流量變化則為干渠2模擬得到的下游流量減去十支渠上游斷面的流量。

3 模擬結果分析

在上述工況下，可模擬得到各渠道上游斷面水位變化和下游斷面流量變化，將所得過渡過程數值結果與用Mike11軟件模擬得到的結果進行對比，來驗證其可靠性，見圖4～圖8。

圖4 干渠1上下游斷面水位流量變化過程Fig.4 Process of water level flow change in the upper and lower sections of the main canal 1

由以上模擬結果可知，本研究中所建立模型與Mike11模型模擬結果大致吻合，說明該模型模擬結果較為準確，并由此得出具體結果如下。

(1)在同時調節干渠一、九支渠以及十支渠上游閘門使其按照模擬條件產生相應流量時，各渠道水位、流量在50 min之后都基本達到了穩定。

(2)干渠一的下游流量在一次波動后達到穩定，在流量線性增大的20 min內，上游水位也表現為線性增長，經過小幅波動后過渡到穩定狀態。

(3)在經過之前的傳播及九支渠的分流后，干渠二的下游流量波動更加明顯，同時干渠二水位變化表現為非線性，在經歷一定波動后達到穩定；干渠三的流量水位變化規律與干渠二類似，只是在經過進一步傳播后，波動較干渠二更為明顯。

(4)在給定模擬條件下，九支渠的流量變化幅度很小，故其下游斷面流量波動很小，隨后逐漸達到穩定。同時在流量線性增大的20 min內，其上游水位也表現為線性增長，之后趨于穩定；十支渠的流量水位變化與九支渠類似，故不再贅述。

圖5 干渠2上下游斷面水位流量變化過程Fig.5 Process of water level flow change in the upper and lower sections of the main canal 2

圖6 干渠3上下游斷面水位流量變化過程Fig.6 Process of water level flow change in the upper and lower sections of the main canal 3

圖7 九支渠上下游斷面水位流量變化過程Fig.7 Process of water level flow change in the upper and lower sections of the nine branch canal

圖8 十支渠上下游斷面水位流量變化過程Fig.8 Process of water level flow change in the upper and lower sections of the ten branch canal

4 閘門運行設計

4.1 閘門開度調節設計

實現灌溉渠道自動化是提高渠系輸配水效率，建設節水農業的重要內容。而其自動化設計中的一個重要環節就是合理的閘門調控設計。灌溉渠道運行時，往往需要根據灌溉流量需求來進行水量輸配，所需水流從水源到田間的過程中受各種因素干擾，引起取水口水深變化，就需要通過相應閘門來進行控制，要得到合理的閘門調控措施，就需要知道相應波動下閘門的開度隨時間的變化。根據本文建立的渠道模型可以求出相應于閘門斷面的水位或流量變化，從而可以根據閘孔出流的公式添加閘門開度計算模型，得出閘門的調控過程，為相關工作人員提供決策依據。

閘孔出流可分為自由出流和淹沒出流，如圖9和圖10所示，它們的流量計算都有對應的經驗公式，可統一為以下格式：

(14)

式中：Qz表示通過閘門的流量；bz表示閘孔寬度；Go表示閘門開度；hu表示閘門上游水深；cd表示流量系數，對于平板閘門，其取值可按下式計算：

圖9 淹沒出流Fig.9 Flooded outflow

圖10 自由出流Fig.10 Free flow

自由出流：

(15)

淹沒出流：

(16)

式中：hd表示閘門下游水深。

在上述引湯灌區分岔渠系模型中，渠道采用下游常水深運行方式，其給定的初始條件中的上游流量的線性變化應該是在下游水位保持不變的情況下由上游閘門調控所產生的，因此根據模擬得到的上游斷面的水位變化過程(即閘后水深變化)，反推出上游閘門在流量線性變化情況下對應的開度隨時間的調節過程。

對于九支渠，首先通過試算得出其初始閘門開度：已知初始閘門上、下游水深分別為1.39 m和0.66 m，閘孔寬度為0.8 m，先假定一個閘門開度，判斷其屬于自由出流還是淹沒出流，根據相應公式計算出其流量系數，之后代入式(14)求出流量并與已知初始流量對比，若不一致則繼續試算至一致為止，通過計算得到九支渠初始閘門開度為0.47 m。得到初始閘門開度后，以此為起點，利用上述試算過程針對上游斷面每一時刻都進行求解，則可得到九支渠上游閘門開度調節過程，如圖11(a)所示。由圖可看出：在前5 min閘門開度不變(水深較小，閘門對水流沒有控制)之后閘門開度隨時間線性增加，到20min后保持不變，其線性調控過程可擬合為以下函數：y=0.007 6x+0.419 8；隨后在本研究中又模擬該渠道在流量從1.32減少到1.02時的閘門開度變化過程，如圖11(b)所示，從圖中可以看到，在前7 min開度保持不變(水深較大，不需調節閘門)，隨后線性減小，到20 min時達到穩定，其線性調控過程可擬合為y=-0.007 9x+0.677 1。

圖11 流量增加或減少時九支渠上游閘門調節過程Fig.11 Regulating process of the upstream gate of the nine branch canal when the flow increases or decreases

同理，對于十支渠也可求出其閘門調控過程。十支渠初始閘門上、下游水深分別為2.49 m和0.63 m，閘孔寬度為0.47 m，試算得其初始閘門開度為0.44 m。十支渠閘門調控過程如圖12(a)所示，在前7 min閘門開度不變，之后閘門開度隨時間線性增加，到20 min后保持不變，其線性調控過程可擬合為以下函數：y=0.006 7x+0.381 9；隨后又模擬了十支渠流量從1.02減少到0.79時的閘門調控過程，如圖12(b)所示，從圖12中可以看到，在前7 min閘門開度保持不變，隨后線性減小，到20 min時達到穩定，其線性調控過程可擬合為y=-0.007 9x+0.677 1，十支渠總體規律與九支渠類似。

圖12 流量增加或減少時十支渠上游閘門調節過程Fig.12 Adjustment process of the upstream gate of the ten branch channel when the flow rate increases or decreases

4.2 多級閘門運行設計

在分岔渠系中，往往存在多級閘門，當要通過閘門調控來滿足下游流量需求時，可以同時調節各級閘門，也可以一級一級分別調節，這樣同步與逐級的兩種調節方式下，渠道中的水流動態也會有所差別。在上述引湯灌區渠道模型中，模擬情況為同時調節干渠一、九支渠和十支渠的上游閘門，在此基礎上，再模擬逐級調節情況，即干渠一達到穩定后調節九支渠，干渠二穩定后調節十支渠，由于干渠一和九、十支渠的模擬工況未變化，因此將模擬得到的干渠二與干渠三的流量變化過程進行對比，如圖13、圖14所示。

圖13 干渠2下游流量變化過程對比Fig.13 Comparison of downstream flow changes in the main canal 2

圖14 干渠3下游流量變化過程對比Fig.14 Comparison of downstream flow changes in the main canal 3

由以上兩圖可以看出：當各級閘門同步調節時，渠道中的水流波動不斷傳播，各級渠道中的波動幅度越來越明顯，而各渠段達到穩定的時間雖有所延長，但卻變化不大，均在50～60 min之后穩定，最后一級渠道即北干渠三在60 min之后達到穩定；當各級閘門逐級調節時，下一級閘門是在上一級渠道中水流達到穩定后開始調節，故渠道中的水流波動無法逐級傳播，對于各級渠道來說，其波動幅度都比較小且穩定時間也比較短，在35～40 min之后達到穩定。但在這種情況下，最后一級渠道穩定到目標量的時間不能僅僅看這一渠道，而是應該從上游第一級閘門的調控開始考慮，于是北干渠三中水流波動達到穩定的時間為自身的40 min加上各上級渠道的穩定時間，則其最終的總穩定時間比同步調控下更長。

目前在我國的大部分灌區還無法實現同步調控，在傳統的運行方式中采用逐級調控，即首先調節渠首的閘門，達到穩定后進行下級閘門調節，按此規律依次推進。在這種情況下，渠道中的水流量波動小，對渠道的沖刷也更小，運行也更安全，同時可降低渠道的超高標準，減少工程造價，但是各級閘門的啟閉互相影響，不利于簡化操作，并且渠系達到最終穩定所需的總時長也會比較長。若采用同步調控，則渠系到達最終穩定的總時長較短，且各級閘門可同時調控，不用以上游渠道的過渡過程為依據，相對獨立，啟閉操作簡單方便，但這種情況下渠道中水流波動會比較大，對渠道產生沖刷，也不利于渠道的安全運行。在實際灌溉中，若時間允許，則推薦逐級調控的方式；但如果需要在盡可能短的時間內達到需求量，則推薦采用同步調控的方式。

5 結 語

本文總結了采用特征線法對圣維南方程進行離散的具體過程、邊界條件、汊點處條件以及穩定性條件，建立了渠系非恒定流數值模擬模型以及閘控模型，且所建模型可以通過修改參數來實現不同渠道的模擬。

(1)本文基于特征線法及汊點處條件，建立了渠系模擬模型，并對引湯灌區北干渠及九、十支渠在給定模擬條件下渠道中所產生非恒定流過渡過程進行模擬。將該模擬結果與Mike11所獲結果比較后發現較為吻合，證明本研究所建模型具有一定可靠性。

(2)依據該模型模擬結果，建立了相應閘控模型，結果表明在流量線性變化的條件下閘門開度變化也是線性的。在本文中通過計算得到九支、十支渠的閘門開度調節過程，得到擬合的閘控方程，作為相關工作人員決策依據，同時實現其智能控制提供依據。

(3)基于本文所建立的渠系模型，模擬各級閘門在同步調控與逐級調控兩種不同方式下的非恒定流過渡過程，結果表明在同步調控時，總調控時間短且此時閘門的啟閉操作方便，但渠道中的水流波動較大，不利于渠道安全；逐級調控時，其水流波動則會比較小且針對各部分渠道來說穩定時間較短，但此種情況下閘門啟閉操作則需滿足上級渠段穩定時長，同時渠系調控的總時長會變長。在實際中，若時間允許則可采用逐級調控，若需在盡量短的時間內得到目標流量則可采用同步調控，但目前我國還有很多灌區無法實現同步調控。

在實際當中，一個灌區的灌溉渠道系統多為由干渠、支渠、斗渠和農渠等多級渠道組成的樹狀網絡，其連接非常復雜，因此，在對簡單的分岔渠系模擬的基礎上，進一步擴展到對更加復雜的大型灌溉渠網的模擬仍是未來工作的重要方向。同時，由于實測資料的欠缺，本文未能將模擬結果與實測資料進行對比驗證。另外，本文所建模型未考慮渠道滲漏損失，對于渠段中有較大底坡跌落的模擬也會失真，這些問題都需要在今后的工作中進行不斷修正與改進。