基于Excel的半图解法在水库调洪计算中的应用

2021-12-08 01:24王新华张慧颖黄辉曦郭美华
人民黄河 2021年11期

王新华 张慧颖 黄辉曦 郭美华

摘 要:相对于列表试算法而言,半图解法避免了大量的试算工作量,在人工调洪计算中广泛应用。使用该方法时要事先绘制水位—库容曲线、泄流量—库容曲线、单辅助曲线,并反复查算曲线,工作量仍然较大,效率不高。在分析半图解法基本原理和调洪计算步骤基础上,提出实现曲线查算的关键步骤是将平面坐标系下各种曲线进行分段,根据x值所处的横坐标区间段,利用线性内插法查出分段曲线上的纵坐标y值。利用Excel中内置的函数trend、offset、match嵌套使用,将曲线查算转化为求分段线性内插值,从而在Excel中实现半图解法调洪计算,大大提高了计算效率。

关键词:半图解法;调洪计算;trend函数;offset函数;match函数; Excel

中图分类号:TV214; TV697   文献标志码:A

doi:10.3969/j.issn.1000-1379.2021.11.014

引用格式:王新华,张慧颖,黄辉曦,等.基于Excel的半图解法在水库调洪计算中的应用[J].人民黄河,2021,43(11):76-80,158.

Application of Semi-Graphic Method Based on Excel in Reservoir Flood Regulation Calculation

WANG Xinhua1, ZHANG Huiying1, HUANG Huixi2, GUO Meihua1

(1. College of Water Conservancy, Yunnan Agricultural University,Kunming 650201, China;

2.Kunming Branch of Yunnan Provincial Hydrology and Water Resources Bureau, Kunming 650051, China)

Abstract: Compared to the list-try algorithm, the semi-graphic method, avoiding a lot of trial calculating work, is widely used in artificial flood regulation calculation. The water level storage curve, discharge and storage curve, auxiliary curve should be drawn in advance when using this method. Curve drawing and repeated search, the workload is still large and the efficiency is not high. Based on the analysis of the basic principle of semi-graphic method and the calculation steps of flood regulation, this paper put forward that the key step to realize the curve calculation and search was to segment all kinds of curves in the plane coordinate system and find out the y-coordinate on the segmented curve according to the x-coordinate interval and x value by linear interpolation method. By using the embedded three functions commands of trend, offset and match, the semi-graphic flood regulation calculation could be realized in Excel by transforming the curves search into finding the linear interpolation of segmented line segments, which greatly improved the calculation efficiency of the method.

Key words: semi-graphic method;flood regulation calculating; trend function; offset function; match function; Excel

水庫的调洪库容、设计洪水位、校核洪水位、泄洪建筑物型式及尺寸拟定和已建水库防洪安全复核等都要进行调洪计算。调洪计算的基本原理是联合运用水量平衡方程和动力平衡方程(水库蓄泄方程),求出逐时段末的水库水位、库容及下泄流量。常用的调洪计算方法有列表试算法、图解法或半图解法[1]。随着计算机的普及以及多种算法的实现,水库调洪出现了一些新的解法,如龙格-库塔数值解法、人工神经网络法、迭代法、双斜率法等[2-4]。数值解法、人工神经网络法、迭代法、双斜率法等适合于计算机编程,但要求开发者有较高的计算机基础和编程能力,且程序的编制和调试十分复杂,对大多数使用者来说,需要花费不少资金购买程序,处理程序运行中出现的错误比较困难。列表试算法概念清晰,容易掌握,但试算工作量较大,效率不高,现在大多是通过软件编程计算代替试算迭代过程。半图解法(又称单辅助线法)通过将水量平衡方程适当变形,结合水位—库容曲线、水位—下泄流量曲线,绘制一条下泄流量q与(VΔt+q2)辅助关系曲线,通过查算该辅助曲线和泄流量与库容关系曲线,求出逐时段末的下泄流量,避免了大量的试算工作量,手工计算时效率较高,在计算机普及之前得到了广泛的应用,但半图解法需要手工绘制q=f(VΔt+q2)辅助线及由V=f(Z)和q=f(Z)转化而来的q=f(V)关系曲线,绘图比例尺选用不当时查图值误差可能较大,如果入库洪水时段改变,需要重新绘制辅助线,工作量依然较大。Excel在日常工作中应用非常普遍,不少学者尝试利用matlab软件、VB程序或Excel自带的VBA编制调洪计算程序[5-7],这仍然需要开发者有一定的编程能力。

本文通过利用Excel中常用的3个函数组合实现曲线查算功能,使得半图解法可以在Excel 中轻松实现,从而快速高效地解决水库调洪计算问题。使用者只要拥有一台安装有Office 软件的电脑,在打开的Excel表格中完成相应公式函数输入,将调洪计算需要的资料数据录入表格中相应位置就可完成调洪计算任务。编制好的调洪计算表格只需再更换一下水位库容资料、溢洪道宽度及流量系数数据、入库洪水过程数据,就可以用于其他水库的调洪计算。

1 半图解法调洪计算基本原理和步骤

水库的蓄泄方程由泄水建筑物型式決定,一般为堰流公式或闸口出流公式,出库流量为堰顶水头H或闸上下游水头差H的函数,而水头H取决于水库水位Z,水位Z又与库容V成单值函数关系,因此下泄流量q可以转化为库容V的单值函数 q=f(V)。

在给定的时段内,水库的水量平衡方程为

Q1+Q22Δt-q1+q22Δt=V2-V1(1)

将水量平衡方程整理移项后可写为

(V2Δt+q22)=Q1+Q22-q1+(V1Δt+q12)(2)

式中:Q1、Q2分别为计算时段初、末的入库流量,m3/s;q1、q2分别为计算时段初、末的出库流量,m3/s;V1、V2分别为计算时段初、末的水库蓄水量,m3;Δt为计算时段长,s。

式(2)中等号右端项是已知项,左端项是未知项,两个括号中都包含V/Δt和q/2,它们都是q的函数,进而可以建立q=f(VΔt+q2) 函数关系,并绘制出q—(VΔt+q2)辅助曲线。

调洪计算时入库洪水过程Q—t、水位库容曲线Z—V、泄洪建筑物型式及尺寸、水库初始蓄水量V1和初始下泄流量q1、计算的时段Δt都是已知的,可以计算出式(2)等号右端各项的值。查辅助曲线就可得到时段末的下泄流量q2,进而得到时段末的蓄水量V2。将计算的时段末q2、V2值作为相邻的下一时段的初值q1、V1,重复以上计算方法,逐时段进行计算、查图就可得到整个下泄流量过程及水库蓄水量的变化过程,再通过查水位库容曲线就可得到水库的水位变化过程。

根据以上半图解法原理可知,调洪计算的核心步骤有4个:

(1)几种曲线的绘制。根据水位库容资料绘制水位库容曲线Z—V;根据泄流公式得到水位与下泄流量关系曲线q—Z;结合水位库容曲线将水位与下泄流量曲线转化为下泄流量与库容关系曲线q—V;结合给定的固定计算时段Δt和曲线q—V,计算出(VΔt+q2),并绘制辅助曲线q—(VΔt+q2)。

(2)时段末下泄流量q2查算。计算式(2)等号右端已知项的值Q1+Q22-q1+(V1Δt+q12),并将该值赋给(V2Δt+q22);查辅助曲线,由(V2Δt+q22)的值查出对应的下泄流量q2。

(3)时段末库容V2推求。 由(V2Δt+q22)的值、q2的值、固定时段Δt的值计算出V2。

(4)时段末水位Z2的查算。由V2的值查水位库容曲线得到相应的Z2值。

2 Excel函数套用实现曲线查算功能

半图解法步骤中数值计算都较为简单,繁琐的任务在于二维坐标系下曲线查算,需要根据水位库容曲线、水位下泄流量曲线转化而成的泄流量库容曲线、单辅助曲线,逐时段查算相应的q2、V2、Z2。曲线查找不外乎是根据已知的x值,查函数值y=f(x)。Excel表格中有大量的内置函数,通过一些函数的套用,可以实现这些查找、计算功能,使繁琐的人工查算由计算机替代,从而高效快捷地完成调洪计算。

为了实现平面直角坐标系下曲线的查找功能,需要在Excel表格中将辅助曲线、下泄流量与库容曲线细分成若干区间段(离散化),每一区间段内的曲线可以用直线近似代替,这样已知x值,先判断x值所处的横坐标处于哪个区间段,然后在该区间段内利用线性内插函数就可以求出y值。

在Excel软件中线性内插函数命令为trend, 其功能是根据已知的y系列数组和x系列数组,得到以y为因变量、x为自变量的线性回归方程,并由新的x值代入回归方程中求取新的y值,其语法为:trend(known_ys,known_xs,new_xs)。其中:known_ys 为已知的y数组系列;known_xs是已知的x数组系列;new_xs是新的x值,由x值代入线性方程可求出待插值y。该函数返回值为y=y1+y2-y1x2-x1(x-x1)。例如在Excel某单元格中输入公式=trend({1,6},{ 9},7),返回值为4.33,相当于已知y1、y2分别为1和6,x1、x2分别为3和9。换句话说,已知直线的两个点的坐标(xi,yi)分别为( 1)、(9,6),则当x等于7 时,该直线上的y坐标值y=y1+y2-y1x2-x1(x-x1)=1+6-19-3(7-3)=4.33。

offset函数功能是以指定单元格为参照系,通过给定偏移量得到新的引用数组区域。该函数返回的引用可以是一个单元格或单元格区域。其语法为:offset(reference,rows,cols,[height],[width])。reference指的是作为偏移引用的参照系;rows为相对于参照单元格向下偏移的行数;cols为相对于参照单元格向右偏移的列数;height为所要引用的区域的行数;width为所要引用区域的列数,省略时视为与参照系的列数相同。例如在某单元格中输入公式=offset(C 2, 2,1),返回的是位于单元格C3下方2行、右侧3列、高度(行数)为2、宽度(列数)为1的单元区域(即F5:F6单元区域)。

match函数功能是在指定单元格区域内搜索指定项,返回该项在单元格区域中的相对排列位置。其语法为match(lookup_value,lookup_array,[match_type])。其中:lookup_value为要在lookup_array区域中查找的值;lookup_array为要搜索的单元格区域;match_type为搜索的类型,取值为1,0,或-1。当match_type取1或省略时,要求lookup_array区域中的数值按升序排列,match函数会查找小于或等于lookup_value值的最大值并给出该最大值在lookup_array数组中的相对位置。例如在某单元格中输入“=match(2.2,{ 2.1, 4})”,则其返回值为3,表示在数组( 2.1, 4)中小于等于2.2的最大值(为2.1)在该数组中的相对排序位置为3。

通过这3个函数的套用就可解决已知曲线的查算问题。例如某水库的水位库容曲线数据存放在Excel表格A、B两列(见表1),在表2的H列中输入新的水位系列值,想要通过A、B两列中的水位库容数据线性内插出H列水位对应的库容值,并将结果存放在表2的J列对应的位置上(比如H4单元格的水位Z为571.9 m,其对应的库容值放在J4单元格中),则只需在J4单元格输入套用函数“=trend(offset(B$1,match(H4,A$2:A$100,1),0,2,1),offset(A$1,match(H4,A$2:A$100,1),0,2,1),H4)”,即可显示库容值为51.27,然后把该单元格中的公式粘贴到该列其他单元格中就可得到水位系列相应的库容值。

3 Excel函数套用在半图解法水库调洪中的应用

已知某小(二)型水库的水位库容曲线、20 a一遇入库洪水过程线,分别存放在同一张表格内A、B、C、D列中。水库汛限水位为正常蓄水位570.90 m,与溢洪道堰顶高程齐平;泄洪建筑物为无闸门控制的宽顶堰,堰顶高程570.90 m,堰宽B=2 m,流量系数m=0.36,调洪计算采用等时段间隔长为1 h,这些数据分别存放在F列单元格$F$1:$F$5区域中。原始数据见表1。

(1)各种曲线绘制数据的生成 。在H列中产生水库起调水位以上等高差的水位序列,在I列生成与H列对应的下泄流量,J列生成与H列水位对应的库容值,K列生成(VΔt+q2)值。具体操作如下:在H2单元格输入起调水位值570.90 m,H3单元格输入公式“=H2+0.50”,然后选中H3单元格,鼠标放到该单元格右下角并下拉若干行,得到以0.50 m为间隔的递增的水位序列;在I2单元格输入堰流计算公式“=$F$4*$F$2*(2*9.8)^0.5*(H3-$F$3)^1.5”, 然后选中该单元格,鼠标放到该单元格右下角并下拉若干行,得到水位序列对应的下泄流量q值;在J2单元格输入套用函数“=trend(offset(B$1,match(H2,A$2:A$100,1),0,2,1),offset(A$1,match(H2,A$2:A$100,1),0,2,1),H2)”, 然后选中该单元格,鼠标放到该单元格右下角并下拉若干行,可得到不同水位Z(或不同下泄流量q)对应的库容V。在K2单元格输入公式“=J2/(0.36*$F$5)+I2/2”, 然后选中该单元格,鼠标放到该单元格右下角并下拉若干行,可得到不同下泄流量q对应的(VΔt+q2)。生成的各种曲线数据成果见表2。

(2)调洪计算表格编制。半图解法调洪计算表格中包含等时段离散化的时段值、逐时段入库洪水流量值Q、时段平均入库流量Q-、逐时段下泄流量q、逐时段末(V2Δt+q22)值、式(2)等号右端项Q1+Q22-q1+(V1Δt+q12)值、逐时段末的水库蓄水量值V2和对应水位值Z2。本文中依次将其放入Excel表格工作簿M、N、O、P、Q、R、S、T列中相应位置。

当水库洪水过程是不等时段间隔时,难以绘制辅助曲线,不能应用半图解法;或者即使是等时段间隔洪水过程,但时段长Δt较大,调洪时可能会错过最大下泄流量与入库洪水过程线的交点,不能给出准确的坝前最高水位。为了解决这两种情况下遇到的问题,需要将入库洪水转化为等时段间距的洪水过程线,时段长度在单元格$F$5中给出,并可根据需要由使用者随时调整,在M、N两列中完成该功能。在M3单元格输入0,表示从0时刻(第1时段初)开始,M4单元格输入公式“=M3+$F$5”,然后选中该M4单元格,鼠标放到该单元格右下角并下拉若干行,完成按等时段长递增的时间序列;在N3单元格输入函数套用公式“=trend(offset($D$1,match(M C$2:C$100),0,2),offset($C$1,match(M C$2:C$100),0,2),M3)”,然后选中该N3单元格,鼠标放到该单元格右下角并下拉若干行,即可求出给定洪水过程的等时段间隔内插出来的流量值。

O列中为计算平均入库流量,即N列中上下相邻两个单元格的平均值,从O4单元格开始录入公式“=(N3+N4)/2”,在该列中其余各行粘贴O4单元格中的公式,就可得到各时段平均入库流量。

P列、Q列分别为下泄流量q及其对应的辅助线值(VΔt+q2),其中P列中的值由同行的Q列值查算得到(即由q—(VΔt+q2)辅助曲线的已知横坐标(VΔt+q2)值查求曲线的纵坐标q值)。在单元格P3中为调洪计算开始时的第一个时段初的下泄流量(为已知值),本例中P3单元格“=$I$2”,即汛限水位对应的下泄流量值,在P4单元格中输入套用函数“=trend(OFFSET(I$1,match(Q4,K$2:K$100),0,2,1),offset(K$1,match(Q4,K$2:K$100),0,2,1),Q4)”,相當于在该单元格完成辅助曲线查找功能,将该单元格公式复制粘贴到该列P4以下其他单元格中就完成了由Q列的(V2Δt+q22)值查辅助曲线得到q2值的任务。Q列存放(VΔt+q2)的值,R列存放Q1+Q22-q1+(V1Δt+q12)的值,根据式(2)得知这两列数值是相等的,但前者要比后者错后一个单元格(一个时段)。在Q3单元格输入已知值(V1Δt+q12)=$K$2,即汛限水位对应的V1Δt+q12值,在R3单元格中输入“=O4-P3+Q3”,相当于计算Q1+Q22-q1+(V1Δt+q12)的值,并将该值赋给Q4单元格,即在Q4单元格中输入“=R3”,表示水量平衡方程式(V2Δt+q22)=Q1+Q22-q1+(V1Δt+q12)。将Q4单元格公式复制粘贴到Q列中Q4单元格以下其余行中,R3单元格中公式复制粘贴到R列中R3单元格以下其余行中,就完成了已知入库洪水过程对应的下泄流量计算。

S列存放时段末库容值、T列存放时段末水位值。根据已经求得的(V2Δt+q22)值(位于Q列)和q2值(位于P列),很容易求得V2=[(V2Δt+q22)-q22]×Δt;水位由水位库容曲线查算得到。在Excel中具体实现方法是:在S3单元格中输入“=(Q3-P3/2)*(0.36*$F$5)”,并下拉到该列其他单元格中完成公式复制,得出各时段末的库容值;在T3单元格中输入套用函数“=trend(offset($A$1,match(S4,B$2:B$100),0,2,1),offset($B$1,match(S4,B$2:B$100),0,2,1),S4)”,并下拉到该列其他单元格中完成公式复制,得到各时段末的水位值。

至此就完成了半图解法水库的调洪计算表格编制和计算。按以上步骤完成单元格数据及公式输入后,几乎在同时就出现调洪计算结果,程序运行时间不超过1 s。表3为本例中时段长为1 h的调洪计算成果。从表1中看出,在给定的入库洪水下,在第5时段末水库达到最高库水位573.09 m,最大下泄流量为10.35 m3/s。但最大下泄流量与该时段末的入库流量8.00 m3/s并不相等,相对差22.7%。可见以固定1 h间隔时段长进行调洪计算时,最高洪水位出现时刻入库流量与下泄流量不相等,成果不尽合理,主要原因是计算时段间隔较大,造成一定的误差。根据这个调洪结果可以判断出水库最高水位理论上应该在第4 h和第5 h之间出现。

为了得到较为精确的最高洪水位出现时刻及其对应的最高洪水位,将单元格F5中的计算时段间隔由1 h改为0.45 h,程序会自动重新计算(VΔt+q2)值,得到新的调洪计算结果,见表4。从调洪结果中可以看出,计算时段缩短为0.45 h后,在4.5 h时入库流量和出库流量分别为10.50、10.43 m3/s,已经非常接近,二者相差不足1%,可以认为二者近似相等,水库最高水位为573.10 m。该时段的调洪得到坝前最高水位和时段为1 h的调洪最高坝前水位相差1 cm,最大下泄流量相差0.08 m3/s,虽然和时段为1 h的调洪计算成果相差不大,但其结果解释起来更合理。

4 结 语

利用Excel函数套用编制的半图解法水库调洪计算程序,不需要掌握复杂的编程语言和程序语句,只需要trend、oeeset、match三个函数简单套用,按本文中的步骤建立好各单元格计算的公式,就可很快得到调洪计算结果。编制好的调洪计算表格程序可以用于其他水库的调洪,只需在表格A、B两列中更换水位、库容数据,C、D两列更换入库洪水过程数据,F列F1~F5中依次更换汛限水位、溢洪道宽度、堰顶高程、流量系数、拟采用的计算时段间隔长度,在不到1 s的时间内就可得到水库调洪计算结果。该方法计算效率很高,程序的通用性强,适合于溢洪道上无闸门控制且以堰顶高程为水库起调水位的水库调洪计算。当溢洪道上有闸门控制出流时,同一水位下可以有不同的闸门开启度,因而可以有多个下泄流量对应某一水位,水位泄流曲线不是单值函数关系,要想使用半图解法,必须事先拟定好调洪调度方式,确定各水位下对应的唯一的下泄流量,得到水位与下泄流量的关系曲线,罗列于H列、I列对应位置,就可得到调洪结果。

参考文献:

[1] 雒文生,宋星原.工程水文及水利计算[M].2版.北京:中国水利水电出版社,2010:265-271.

[2] 郭世兴,王光社.水庫调洪计算方法的应用研究[J].人民黄河,2016,38(7):24-26.

[3] 颜婷莉,钟平安,刘伟莉.水库调洪演算方法比较与改进[J].水力发电,2007,33(3):26-28.

[4] 廖小龙,陈俊贤,朱毅峰.一种简单实用的水库调洪数值-解析法[J].水力发电学报,2014,33(5):44-47.

[5] 黄启有.单辅助线法在水库调洪计算中的应用研究[J].人民珠江,2018,39(3):54-56.

[6] 雷静思,郭纯青.VB程序与Excel结合在水库调洪演算中的应用[J].广西水利水电,2015(4):48-51.

[7] 安航永,赵文龙,张雅萍.基于半图解法的MATLAB程序和Excel在水库调洪演算中的联合运用[J].陕西水利,2019(4):69-71,76.

【责任编辑 许立新】