梁 建，祁 濤，盧 俊

（安徽省·水利部淮委水利科學研究院，安徽合肥230088）

變步長龍格庫塔方法在安徽中小水庫調(diào)洪中的應用

梁 建，祁 濤，盧 俊

（安徽省·水利部淮委水利科學研究院，安徽合肥230088）

傳統(tǒng)的水庫調(diào)洪方法原理是通過圖解和試算求解相關(guān)水量平衡方程，該方法計算工作量大、效率低，在編程靈活性、調(diào)算精度方面存在一定的局限性。針對傳統(tǒng)水庫調(diào)洪方法存在的問題，將水量平衡方程與泄流曲線方程歸納成為一個一階常微分方程，根據(jù)邊界條件，用變步長龍格庫塔方法，求出該方程的數(shù)值解。該方法易于程序設計，計算精度高。基于C＃，結(jié)合安徽省“84辦法”，用EXCEL做前后處理，將龍格庫塔方法和樣條曲線插值用于安徽中小水庫調(diào)洪中。采用變步長龍格庫塔方法，可以更接近峰值。算例表明，本方法精度高，可以用于安徽中小水庫調(diào)洪計算。

調(diào)洪；變步長龍格庫塔；84辦法

水庫調(diào)洪計算在水庫除險加固、水庫運行管理中起著至關(guān)重要的作用。傳統(tǒng)的調(diào)洪計算具有工作量大、效率低、編程靈活性弱、計算精度低等明顯缺陷。姜樹海［1］通過引入隨機數(shù)學模型，提出運用隨機微分方程進行水庫的調(diào)洪演算；陳元芳［2］采用Monta－Carlo方法對隨機模擬法及傳統(tǒng)方法推求設計防洪庫容優(yōu)劣進行對比研究，結(jié)果表明其計算結(jié)果比傳統(tǒng)的同頻率法和同倍比法精度要高；李傳奇等［3］采用LHS－MC方法計算得到最高調(diào)洪水位隨機分布，以此衡量漫壩風險。以上基于概率統(tǒng)計的研究成果提高調(diào)洪演算的計算精度，但其計算準確的前提在于前期樣本的大小和準確性，而現(xiàn)有中小型水庫，由于管理不善往往缺少水文資料，因此上述計算手段精度無法保證，甚至無法使用。

本文在前人工作的基礎(chǔ)上［4－8］，結(jié)合安徽省水利水電勘測設計院1984年編制的《安徽省暴雨參數(shù)等值線圖、山丘區(qū)產(chǎn)匯流分析成果和山丘區(qū)中、小面積設計洪水計算辦法》（安徽省“84辦法”主要用于淮河以南，集水面積在10 km2～300 km2之間的山丘區(qū)中、小型水庫計算設計洪水之用，以下簡稱“84辦法”）［9］，將變步長龍格庫塔方法用于安徽省中小型水庫調(diào)洪演算，計算結(jié)果表明該方法能很好反應實際情況，可對水庫除險加固、水庫運行起到指導作用。

1 “84辦法”計算洪水過程線

“84辦法”計算洪水過程線的方法屬于綜合瞬時單位線法，是安徽省針對中、小面積設計洪水進行計算的一種辦法?！?4辦法”是根據(jù)安徽省50個站的305次洪水，用納希瞬時單位線為模型，分析并綜合計算瞬時單位線的地區(qū)經(jīng)驗。納希模型主要有兩個參數(shù)，N和k，納希認為流域上瞬時單位凈雨經(jīng)過N個相同的串聯(lián)線性水庫調(diào)蓄后出流過程就是流域出口斷面處的流量過程線。k為每個水庫的調(diào)節(jié)系數(shù)或調(diào)蓄滯時。在地形圖上量算某水庫流域特征值流域面積F、流域平均寬度 B、河道平均坡度 J，在“84辦法”附圖中查讀設計流域中心處的某設計頻率24 h、1 h暴雨點—面折扣系數(shù)和年最大24 h、1 h點雨量的變差系數(shù)，經(jīng)過計算可得該設計頻率下的洪水過程線。

2 水庫調(diào)洪的數(shù)值解析解法

2.1 泄水建筑物匯流公式

水利工程中，常修建水閘或溢流壩等建筑物方便泄洪，以控制河流或水庫的水位及流量，此時的水流狀態(tài)分為閘孔出流或者堰流［10］。

（1）閘孔出流的計算公式

式中：M1為流量常數(shù)；C1為泄流孔口底檻高程；Z為流量常數(shù)；A1為泄流孔口高；B1為泄流孔口寬。

（2）堰流的計算公式

式中：M2為流量常數(shù)；C2為泄流孔口底檻高程；B2為泄流孔口寬。

2.2 數(shù)值解析法調(diào)洪計算原理

數(shù)值解析法用數(shù)值方法對水庫水量平衡微分方程［11－12］進行求解。微分方程可以采用以下兩種方式：

（1）采用水位—庫面面積曲線

式中：Z為水庫水位；Z0為起調(diào)水位；F（Z）為水庫水位對應的水面面積；Q（t）為入庫流量；q（Z）為庫水位對應的出庫流量。

（2）采用水位—庫容曲線

式中：V為水庫庫容；V0為起調(diào)庫容；Q（t）為入庫流量；q（V）為庫水位對應的出庫流量。

根據(jù)“84辦法”或者已知洪水過程線、泄流曲線以及水位 ～庫容關(guān)系曲線，就可以求解水庫各時段的泄流過程線，根據(jù)泄水過程線，通過樣條曲線插值計算，可以得到庫水位的變化過程線。常用的求解偏微分方程的方法很多，四階龍格庫塔方法是精度較高的一種［13］。

2.3 變步長龍格庫塔方法簡述

（1）龍格庫塔方法

（2）變步長龍格庫塔方法

固定步長的龍格庫塔在步長比較大的時候，容易錯過峰值，取步長比較小的時候，又增大了計算量。通過加倍或折半處理步長的方法稱為變步長方法，表面上看，為了選擇步長，對初始步長Δt0進行修正，通過變步長方法，實現(xiàn)快速獲取滿足精度的下一時間步的函數(shù)值。每一步的計算量增加了，但總體考慮往往是合算的［14］。

Δ為檢查步長折半前后兩次計算結(jié)果的偏差：

①假設精度ε，如果Δ＞ε，反復將步長折半進行計算，直至Δ ＜ε為止，這時取最終得到的作為結(jié)果；

②相反當Δ＜ε，反復將步長加倍，直到Δ＞ε為止，這時再將步長折半一次，就得到所要的結(jié)果。

3 程序設計

本程序在設計中可以一次性計算設計、校核和施工期洪水三種工況下的調(diào)洪計算。計算中首先由“84辦法”計算洪水過程線，并將結(jié)果寫入Excel［15］。當調(diào)洪時，由Excel讀入洪水過程線、庫容曲線、泄流曲線等，由變步長龍格庫塔方法進行調(diào)洪，最后將結(jié)果寫入Excel，在Excel中完成后處理（見圖1）。

圖1 程序流程圖

4 算 例

4.1 安徽某水庫主要參數(shù)

以安徽省東部某小一型水庫為例，水庫集水面積2.83 km2，總庫容138.5×104m3。溢洪道為寬頂堰，寬4m。水庫水位庫容關(guān)系見表1。根據(jù)規(guī)范，按30年一遇洪水設計，300年一遇洪水校核。用本程序計算的水庫30年一遇洪水過程線見表2。程序界面設計見圖2。

圖2 程序界面設計圖

4.2 程序計算結(jié)果分析

用本程序和傳統(tǒng)試算方法進行計算，結(jié)果見表3，通過對比可以看出，本程序與傳統(tǒng)程序可以精確到毫米單位，滿足工程需要精度。

表1 水庫水位庫容關(guān)系

表2 水庫30年一遇洪水過程線

入庫洪水過程線，下泄洪水過程線，水位曲線。用本文程序計算結(jié)果見圖3，從表3可以看出，該水庫30年一遇的高水位為54.985 m，對應的下泄流量為5.968m3／s，對應庫容為116.37×104m3，發(fā)生在第13個時段。

5 結(jié) 語

通過以上水庫調(diào)洪實例，可以得出，在數(shù)值解調(diào)洪計算中，當調(diào)整泄洪建筑物參數(shù)時，不需要重新繪制工作曲線。使用這種數(shù)值解析法可以節(jié)省機間，對優(yōu)化水庫溢洪道規(guī)模及參數(shù)等帶來極大的方便。本程序計算結(jié)果精度上滿足工程要求。將變步長龍格庫塔方法用于中小水庫調(diào)洪中，可以解決中小水庫除險加固和運行的計算問題。變步長可以比試算更接近峰值。結(jié)合“84辦法”，可將程序用于安徽省中小水庫調(diào)洪中。本程序尚存缺陷，將在使用過程中，不斷完善。

圖3 調(diào)洪計算結(jié)果圖

表3 調(diào)洪計算結(jié)果對比表

［1］ 姜樹海，范子武.水庫防洪預報調(diào)度的風險分析［J］.水利學報，2004，5（11）：102-107.

［2］ 陳元芳.隨機模擬法與傳統(tǒng)方法推求設計防洪庫容優(yōu)劣的初步研究［J］.水科學進展，2000，3（11）：64-69.

［3］ 李傳奇，王 帥，王 薇，等.LHS－MC方法在漫壩風險分析中的應用［J］.水力發(fā)電學報，2012，31（2）：5-9.

［4］ 陳守煜.水庫調(diào)洪數(shù)值—解析解法［J］.大連理工大學學報，1996，11（6）：721-724.

［5］ Chen Shouyu.Relativemembership function and new frame of fuzzy sets theory for pattern recognition［J］.The Journal of Fuzzy Mathematies，1997，5（2）：401-411.

［6］ Chen Shouyu.Non-structured decision making analysis and fuzzy optimum seeking theory formulti-objective systems［J］.Journal of Fuzzy Mathematies，1996，4（2）：835-842.

［7］ 陳守煜.系統(tǒng)模糊決策理論與應用［M］.大連：連理工大學出版社，1994.

［8］ 陳守煜，周惠成.黃河防洪決策支持系統(tǒng)多目標多層次對策方案的模糊優(yōu)選［J］.水電能源科學，1992，10（2）：94-101.

［9］ 安徽省暴雨參數(shù)等值線圖、山丘區(qū)產(chǎn)匯流分析成果和山丘區(qū)中、小面積設計洪水計算辦法［M］.合肥：安徽省水利水電勘測設計院，1984.

［10］ 吳持恭.水力學（第3版）［M］.高等教育出版社，2003.

［11］ 郭生練.水庫調(diào)度綜合自動化系統(tǒng)［M］.武漢：武漢水利電力大學出版社.

［12］ 水電規(guī)劃設計院主編.水利動能設計手冊（防洪分冊）［M］.水利電力出版社，1988.

［13］ 付 磊，張洪明，姚 激.水庫調(diào)洪演算的數(shù)值解析法［J］.水利與建筑工程學報，2006，4（4）：75-77.

［14］ 李慶揚，王能超，易大義.數(shù)值分析（第5版）［M］.北京：清華大學出版社，2008.

［15］ 張能立，萬 欲，王 睿.ASP.NET技術(shù)及其在企業(yè)辦公信息系統(tǒng)中的應用［J］.計算機與現(xiàn)代化，2005，（8）：121-122，126.

Application of Varying-step Runge-Kutta Method in Small and Medium Reservoir Flood Regulation of Anhui

LIANG Jian，QITao，LU Jun
（Anhui Hydraulic Research Institute，HuaiheWater Resources Commission of MWR，Hefei，Anhui 230088，China）

The traditional principle of flood routing through reservoirmethod is solving thewater balance equation through diagramming and trialmethod，which would have heavy workflows，low efficiency and limitations both on programming flexibility and simulation accuracy.A newmethod is proposed by combining thewater balance equation and flow equation into a first-order ordinary differential equation so as to calculate its numerical solution according to the varying-step Runge-Kuttamethod and boundary conditions.Thismethod is easy to program and with high precision of calculation.Here，based on C＃，combined with the 84 regulation in Anhui Province and using the EXCEL for pre-post-processing，the Runge-Kuttamethod and spline interpolation are applied to a small andmedium reservoir flood routing in Anhui.The results show that the Runge-Kuttamethod can effectively avoid themissing of peak value.The calculation examples show that the proposedmethod has a high precision and it is applicable to the calculations of small andmedium reservoir flood routing.

reservoir flood routing；varying-step Runge-Kuttamethod；84 regulation in Anhui Province

TV68

A

1672—1144（2014）01—0193—04

10.3969／j.issn.1672－1144.2014.01.040

2013-09-10

2013-10-02

梁 建（1984—），男，山東聊城人，碩士，主要從事水工結(jié)構(gòu)及其模擬方面的研究。