梁 建，祁 濤，盧 俊

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

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

梁 建，祁 濤，盧 俊

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

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

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

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

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

1 “84辦法”計(jì)算洪水過程線

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

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

2.1 泄水建筑物匯流公式

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

（1）閘孔出流的計(jì)算公式

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

（2）堰流的計(jì)算公式

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

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

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

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

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

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

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

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

2.3 變步長(zhǎng)龍格庫塔方法簡(jiǎn)述

（1）龍格庫塔方法

（2）變步長(zhǎng)龍格庫塔方法

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

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

①假設(shè)精度ε，如果Δ＞ε，反復(fù)將步長(zhǎng)折半進(jìn)行計(jì)算，直至Δ ＜ε為止，這時(shí)取最終得到的作為結(jié)果；

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

3 程序設(shè)計(jì)

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

圖1 程序流程圖

4 算 例

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

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

圖2 程序界面設(shè)計(jì)圖

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

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

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

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

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

5 結(jié) 語

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

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

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

［1］ 姜樹海，范子武.水庫防洪預(yù)報(bào)調(diào)度的風(fēng)險(xiǎn)分析［J］.水利學(xué)報(bào)，2004，5（11）：102-107.

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

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

［4］ 陳守煜.水庫調(diào)洪數(shù)值—解析解法［J］.大連理工大學(xué)學(xué)報(bào)，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)模糊決策理論與應(yīng)用［M］.大連：連理工大學(xué)出版社，1994.

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

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

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

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

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

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

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

［15］ 張能立，萬 欲，王 睿.ASP.NET技術(shù)及其在企業(yè)辦公信息系統(tǒng)中的應(yīng)用［J］.計(jì)算機(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)及其模擬方面的研究。