变步长龙格库塔方法在安徽中小水库调洪中的应用

梁 建，祁 涛，卢 俊

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

变步长龙格库塔方法在安徽中小水库调洪中的应用

梁 建，祁 涛，卢 俊

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

传统的水库调洪方法原理是通过图解和试算求解相关水量平衡方程，该方法计算工作量大、效率低，在编程灵活性、调算精度方面存在一定的局限性。针对传统水库调洪方法存在的问题，将水量平衡方程与泄流曲线方程归纳成为一个一阶常微分方程，根据边界条件，用变步长龙格库塔方法，求出该方程的数值解。该方法易于程序设计，计算精度高。基于C＃，结合安徽省“84办法”，用EXCEL做前后处理，将龙格库塔方法和样条曲线插值用于安徽中小水库调洪中。采用变步长龙格库塔方法，可以更接近峰值。算例表明，本方法精度高，可以用于安徽中小水库调洪计算。

调洪；变步长龙格库塔；84办法

水库调洪计算在水库除险加固、水库运行管理中起着至关重要的作用。传统的调洪计算具有工作量大、效率低、编程灵活性弱、计算精度低等明显缺陷。姜树海［1］通过引入随机数学模型，提出运用随机微分方程进行水库的调洪演算；陈元芳［2］采用Monta－Carlo方法对随机模拟法及传统方法推求设计防洪库容优劣进行对比研究，结果表明其计算结果比传统的同频率法和同倍比法精度要高；李传奇等［3］采用LHS－MC方法计算得到最高调洪水位随机分布，以此衡量漫坝风险。以上基于概率统计的研究成果提高调洪演算的计算精度，但其计算准确的前提在于前期样本的大小和准确性，而现有中小型水库，由于管理不善往往缺少水文资料，因此上述计算手段精度无法保证，甚至无法使用。

本文在前人工作的基础上［4－8］，结合安徽省水利水电勘测设计院1984年编制的《安徽省暴雨参数等值线图、山丘区产汇流分析成果和山丘区中、小面积设计洪水计算办法》（安徽省“84办法”主要用于淮河以南，集水面积在10 km2～300 km2之间的山丘区中、小型水库计算设计洪水之用，以下简称“84办法”）［9］，将变步长龙格库塔方法用于安徽省中小型水库调洪演算，计算结果表明该方法能很好反应实际情况，可对水库除险加固、水库运行起到指导作用。

1 “84办法”计算洪水过程线

“84办法”计算洪水过程线的方法属于综合瞬时单位线法，是安徽省针对中、小面积设计洪水进行计算的一种办法。“84办法”是根据安徽省50个站的305次洪水，用纳希瞬时单位线为模型，分析并综合计算瞬时单位线的地区经验。纳希模型主要有两个参数，N和k，纳希认为流域上瞬时单位净雨经过N个相同的串联线性水库调蓄后出流过程就是流域出口断面处的流量过程线。k为每个水库的调节系数或调蓄滞时。在地形图上量算某水库流域特征值流域面积F、流域平均宽度 B、河道平均坡度 J，在“84办法”附图中查读设计流域中心处的某设计频率24 h、1 h暴雨点—面折扣系数和年最大24 h、1 h点雨量的变差系数，经过计算可得该设计频率下的洪水过程线。

2 水库调洪的数值解析解法

2.1 泄水建筑物汇流公式

水利工程中，常修建水闸或溢流坝等建筑物方便泄洪，以控制河流或水库的水位及流量，此时的水流状态分为闸孔出流或者堰流［10］。

（1）闸孔出流的计算公式

式中：M1为流量常数；C1为泄流孔口底槛高程；Z为流量常数；A1为泄流孔口高；B1为泄流孔口宽。

（2）堰流的计算公式

式中：M2为流量常数；C2为泄流孔口底槛高程；B2为泄流孔口宽。

2.2 数值解析法调洪计算原理

数值解析法用数值方法对水库水量平衡微分方程［11－12］进行求解。微分方程可以采用以下两种方式：

（1）采用水位—库面面积曲线

式中：Z为水库水位；Z0为起调水位；F（Z）为水库水位对应的水面面积；Q（t）为入库流量；q（Z）为库水位对应的出库流量。

（2）采用水位—库容曲线

式中：V为水库库容；V0为起调库容；Q（t）为入库流量；q（V）为库水位对应的出库流量。

根据“84办法”或者已知洪水过程线、泄流曲线以及水位 ～库容关系曲线，就可以求解水库各时段的泄流过程线，根据泄水过程线，通过样条曲线插值计算，可以得到库水位的变化过程线。常用的求解偏微分方程的方法很多，四阶龙格库塔方法是精度较高的一种［13］。

2.3 变步长龙格库塔方法简述

（1）龙格库塔方法

（2）变步长龙格库塔方法

固定步长的龙格库塔在步长比较大的时候，容易错过峰值，取步长比较小的时候，又增大了计算量。通过加倍或折半处理步长的方法称为变步长方法，表面上看，为了选择步长，对初始步长Δt0进行修正，通过变步长方法，实现快速获取满足精度的下一时间步的函数值。每一步的计算量增加了，但总体考虑往往是合算的［14］。

Δ为检查步长折半前后两次计算结果的偏差：

①假设精度ε，如果Δ＞ε，反复将步长折半进行计算，直至Δ ＜ε为止，这时取最终得到的作为结果；

②相反当Δ＜ε，反复将步长加倍，直到Δ＞ε为止，这时再将步长折半一次，就得到所要的结果。

3 程序设计

本程序在设计中可以一次性计算设计、校核和施工期洪水三种工况下的调洪计算。计算中首先由“84办法”计算洪水过程线，并将结果写入Excel［15］。当调洪时，由Excel读入洪水过程线、库容曲线、泄流曲线等，由变步长龙格库塔方法进行调洪，最后将结果写入Excel，在Excel中完成后处理（见图1）。

图1 程序流程图

4 算 例

4.1 安徽某水库主要参数

以安徽省东部某小一型水库为例，水库集水面积2.83 km2，总库容138.5×104m3。溢洪道为宽顶堰，宽4m。水库水位库容关系见表1。根据规范，按30年一遇洪水设计，300年一遇洪水校核。用本程序计算的水库30年一遇洪水过程线见表2。程序界面设计见图2。

图2 程序界面设计图

4.2 程序计算结果分析

用本程序和传统试算方法进行计算，结果见表3，通过对比可以看出，本程序与传统程序可以精确到毫米单位，满足工程需要精度。

表1 水库水位库容关系

表2 水库30年一遇洪水过程线

入库洪水过程线，下泄洪水过程线，水位曲线。用本文程序计算结果见图3，从表3可以看出，该水库30年一遇的高水位为54.985 m，对应的下泄流量为5.968m3／s，对应库容为116.37×104m3，发生在第13个时段。

5 结 语

通过以上水库调洪实例，可以得出，在数值解调洪计算中，当调整泄洪建筑物参数时，不需要重新绘制工作曲线。使用这种数值解析法可以节省机间，对优化水库溢洪道规模及参数等带来极大的方便。本程序计算结果精度上满足工程要求。将变步长龙格库塔方法用于中小水库调洪中，可以解决中小水库除险加固和运行的计算问题。变步长可以比试算更接近峰值。结合“84办法”，可将程序用于安徽省中小水库调洪中。本程序尚存缺陷，将在使用过程中，不断完善。

图3 调洪计算结果图

表3 调洪计算结果对比表

［1］ 姜树海，范子武.水库防洪预报调度的风险分析［J］.水利学报，2004，5（11）：102-107.

［2］ 陈元芳.随机模拟法与传统方法推求设计防洪库容优劣的初步研究［J］.水科学进展，2000，3（11）：64-69.

［3］ 李传奇，王 帅，王 薇，等.LHS－MC方法在漫坝风险分析中的应用［J］.水力发电学报，2012，31（2）：5-9.

［4］ 陈守煜.水库调洪数值—解析解法［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］ 陈守煜.系统模糊决策理论与应用［M］.大连：连理工大学出版社，1994.

［8］ 陈守煜，周惠成.黄河防洪决策支持系统多目标多层次对策方案的模糊优选［J］.水电能源科学，1992，10（2）：94-101.

［9］ 安徽省暴雨参数等值线图、山丘区产汇流分析成果和山丘区中、小面积设计洪水计算办法［M］.合肥：安徽省水利水电勘测设计院，1984.

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

［11］ 郭生练.水库调度综合自动化系统［M］.武汉：武汉水利电力大学出版社.

［12］ 水电规划设计院主编.水利动能设计手册（防洪分册）［M］.水利电力出版社，1988.

［13］ 付 磊，张洪明，姚 激.水库调洪演算的数值解析法［J］.水利与建筑工程学报，2006，4（4）：75-77.

［14］ 李庆扬，王能超，易大义.数值分析（第5版）［M］.北京：清华大学出版社，2008.

［15］ 张能立，万 欲，王 睿.ASP.NET技术及其在企业办公信息系统中的应用［J］.计算机与现代化，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—），男，山东聊城人，硕士，主要从事水工结构及其模拟方面的研究。