借助Excel电子表格实现水库泄流全过程调洪演算
【内容提要】 借助Excel电子表格实现水库泄流全过程调洪演算，以井研县毛坝水库枢纽病害整治工程初步设计中，其上游水库红岩水库p=0.1%校核洪水过程的调洪演算为算例，探索出两种不同的计算方法：一种是利用Excel电子表格自身的函数和公式运算，实现半自动化的水库调洪演算；二是以Microsoft Excel2010以上版本为依托，通过VBA语言编程来实现全自动化的水库调洪演算。两种方法的计算结果，与本人当年使用日本夏普便携式PC-1500计算机通过Basic语言编程的水库多功能调洪演算程序的计算结果完全一致。作者推荐第二种方法，它能在大约1至2秒之内实现多行多列上百数据数百次试算的入库和出库洪水流量演算全过程。由于它的快速和智能化设计，给水库洪水预报盈得时间，是编制已成水库防洪预案和水利水电工程设计及病害整治工作的有力助手。
【正 文】现以井研县毛坝水库1995年工程整治中红岩水库调洪演算为算例，展现更先进的计算方法。本算例在计算入库洪水过程中，所涉及暴雨洪水参数来源有二：一是四川省水利电力局水文总站1979年10月编制的《四川省水文手册》；二是四川省水利电力厅1984年6月编制的《四川省中小流域暴雨洪水计算手册》。虽然后者的暴雨等值线图已经被近20多年来新增实测
暴雨资料的补充分析成果所代替。四川省水文水资源勘测局已增补到2000年的暴雨资料，于2010年12月出版了新的《四川省暴雨统计参数等值线图》，按理应以新的暴雨资料为依据，但为了保持该工程整治设计的原貌，为了与当时本人使用PC-1500计算机和原来的调洪演算程序的运算结果进行比较，仍然使用《四川省中小流域暴雨洪水计算手册》所附暴雨等值线图计算该水库设计洪水过程线。
一、水库的基本情况
毛坝水库位于岷江左支泥溪河上游，井研县天云乡境内。水库大坝的地理坐标为E103°57′24″，N29°45′37″。大坝以上总集水面积34.3km2，其上游还有一座小（一）型水库名红岩水库，集水面积10.8km2，，主河道长6.40km，河道平均比降12.2‰ 。红岩水库的洪水计算标准采用与毛坝水库相同的防洪标准，即五十年一遇（p=2%）设计，千年一遇（p=0.1%）洪水校核。为了节省篇幅，本文只以红岩水库p=0.1%校核洪水的调洪演算过程为算例，不涉及其他的计算内容。
红岩水库坝址在毛坝水库大坝上游沿主河道长11.3km。大坝在1970年建成后，因坝坡较陡，在1979年6月天云乡发生地震的同时，产生大滑坡（主要在内坡），后根据乐山地区水电局（81）120号批文，决定将大坝高程由498m降低到495m，并削坡整治，其溢洪道堰顶高程相应降低到492.4m，该水库在1981年11月动工整治，到1982年2月经费用光，劳动力下马，但工程并未达到原批文要求。根据1994年11月27日实测大坝最低点高程为495.13m。溢洪道为9孔浆砌条石弧形薄壁堰，每孔净宽（弦长）5m，其弧长为6.04m，各孔平均堰顶高程为492.36m，总净宽B=5×9=45m，按此堰型取第二流量系数M=1.7，其泄流量计算公式为Q=MBH1.5=76.5H1.5，式中H为堰顶以上水头，以米计。
根据《四川省中小流域暴雨洪水计算手册》计算得红岩水库p=0.1%的入库洪峰流量为217m3/s，单峰洪水总量为463.32万m3，其矩形概化历时为5.94h，基流q0=1.62m3/s 。
二、水库调洪演算原理及计算步骤
水库调洪演算的目的，是从水库的入库流量过程，运用水量平衡原理，逐时段演算出库流量过程及其相应的库内水位。
设时间ti-1及ti的入、出库流量各为Qi-1、Qi及qi-1、qi；相应的库容为Vi-1、Vi，根据水量平衡原理，可以建立以下方程式：
0.5（Qi-1+Qi）△t-0.5（qi-1+qi）△t=Vi-Vi-1                                 （1）
即qi=2（Vi-1-Vi）/△t+Qi-1+Qi-qi-1                                            （2）
式中△t——时段历时，△t=ti-ti-1  若将（2）式中的各项按常用的单位计算，即流量单位为m3/s，库容单位为万m3 ，时间单为小时，则（2）式改变为：
qi=5.556（Vi-1-Vi）/△t-qi-1+Qi-1+Qi                                          （3）
若取△t=0.5小时，则（3）式可改写为
qi=11.111（Vi-1-Vi）-qi-1+Qi-1+Qi                                            （4）
（3）～（4）式中因Vi为未知项，只能先假定时段末水位Zi，在事先置入Excel电子表格中的水位～库容量关系序列中插补计算Vi，再由（3）或（4）式计算qi。又由于出库流量qi与水位Zi的关系同时受到溢洪道堰顶高程ZS、净宽B及与堰型有关的第二流量系数M的制约，qi=MB（Zi- ZSvim命令列间编辑）1.5 。
即  Zi=（qi/M/B）2/3+ZS                                                      （5）
为了与原假定的Zi加以区别，将（5）式改为
Zi，=（qi/M/B）2/3+ZS， 式中B=45m，M=1.70，于是：
Zi，= ZS+（0.013072×qi）2/3                                              (6)
坝前水位的单位为m，若计算水位Zi，与原假定水位Zi各取小数后2位（第3位四舍五入）之值相等，则Zi，即为试算终值。其对应的Zi、Vi及qi即为时间ti的库内水位、库容量、出库流量。
上述试算方法实际操作起来很费时间，最熟练人员也要两个工作日，使用“蓄率中线”图解法需要一个工作日，但水位的试算精度达不到厘米级。本人经过长期实践，提出以下两种方法，可以准确快速得到试算成果。
三、利用Excel电子表格的公式和函数实现半自动化的水库调洪演算
（一）准备工作：
1、 打开Excel电子工作簿，其“sheet1”工作表立即展现出来，见附图（一）。合并B2～G2单元格，输入第一行表头“p=0.1%红岩水库校核洪水过程线”
2、 分别合并B3与C3、D3与E3、F3与G3，输入第二行表头“水位～库容关系”、“洪水概化过程线”、“设计洪水过程”。
3、 分别在B4～G4单元格输入第三行表头“水位”、“库容”、“y”、“x”、“时间t”、“流量Qt”。
4、 将红岩水库的水位～库容关系依次输入到B5～C28区域。
5、 将《手册》94页附表3-2“ 四川省小流域典型洪水概化过程线综合成果表”中的“东部地区”
单峰（2）型洪水的相对坐标的y与x对应值输入到D5～E21区域的单元格内。
6、 将“集水面积”及“基流”等数据输入到D22～E26单元格区间。
7、 分别将“溢洪道堰顶”高程及其库容量、“坝顶”高程及其库容量输入到B29～D30单元格区间。
8、 将有关计算公式输入到B31～H34单元格区间。
9、 用公式t=x*TP计算入流过程线各节点的时间。具体作法是先点击F5单元格，在有ƒχ的“公式编辑栏”内输入=E5*$E$25,接着点击“公式编辑栏”ƒχ前面的“√”，于是F5单元格的值出现为“0”，而该公式仍保留在“公式编辑栏”内。此时单元格F5仍被“黑框”包围，表示F5仍处于被“点击”状态，接着点击“复制”图标，再单击单元格F6并将鼠标按住不放，一直拖到F21，于是从F6至F21的单干元格全部变成蓝，然后点击“粘贴”图标，从此F5单元格的计算公式全部被“相对引用”到这些单元格，它们的计算值全部被展现出来。现以F8单元格为例，其值为1.19，点击F8，它的公式就出现在“公式编辑栏”内，=E8*$E$25,其中E8是F8所在第8行的相对引用，而$E$25则是对E25单元格数值5.94的绝对引用。
10、 用公式Qt=У*QP+Q0计算入流过程线各节点的入库流量。点击G5单元格，在有ƒχ的“公式编辑栏”内输入=D5*$E$23+$E$26,其中$E$23是对E23单元格“入库洪峰流量QP”的绝对引用，其值为217。$E$26是对E26单元格“基流Q0”的绝对引用，其值为1.62。接着点击“公式编辑栏”ƒχ前面的“√”，于是G5单元格的值出现为1.62。按照第⑨项的办法将单元格G5的公式“相对复制”到单元格G6～G21，则G列6～21行各单元格的值全部出现到位。
（二）调洪演算
1、在sheet1工作表中合并I2～U3区域各单元格，输入“P=0.1﹪红岩水库入、出库流量调洪演算过程”，作为该部分第一行标题。另在I4～U4区域各单元格按照附图（一）的内容输入各列的单项标题。
2、从单元格I5～I35输入0～30的序号。
3、在单元格J4“时间t”栏目下，单元格J5～J35依次输入-0.4、0.1、0.6、1.1、……14.6，即以0.5小时为1时段的时间值。
4、在单元格K5～R5内依次输入1.62、0、492.36、479.57、0、492.36、0、492.36，它们
分别为洪水入库前坝前水位蓄到溢洪道堰顶时的开始状态。其中1.62是洪水入库前的基流量，492.36是蓄满到溢洪道堰顶的水位，479.57是溢洪道堰顶的库容量。
5、根据J列6～35行的时间t值在“设计洪水过程”栏目的F列和G列各行t和Qt节点值内插出相应的入库流量Q，单元格K6的内插公式为=G$5+(J6-F$5)*(G$6-G$5)/(F$6-F$5)，K6的值为3.28。内插的方法是采用线性内插，由于J6=0.1，它位于F列时间节点0.00与0.65之间，J6的值＞F5，又＜F6，因此在G$5数值的基础上，再增加一个比例增值，这个增值就是（G$6-G$5）乘以（J6-F$5）/(F$6-F$5)的比值。在行号的前面加$符号，是对该行的“绝对引用”。同样K7的公式要看J7的值位于F列哪两个单元格之间，因J7=0.6，它仍然位于F5和F6之间，则K7的公式同样可从K6的公式直接复制下来，成为=G$5+(J7-F$5)*(G$6-G$5)/(F$6-F$5)。此式与K6的公式比较，只是将J6改换成J7，这是对J6的“相对引用”，即从6行变成7行。K8的公式是=G$7+(J8-F$7)*(G$8-G$7)/(F$8-F$7)，因为J8=1.1位于F7与F8之间，故以G7的流量为基础，再加一个比例增值(J8-F$7)*(G$8-G$7)/(F$8-F$7)。K8的公式同样可从K6的公式复制下来，只修改相关的行号，可以节省工作量。K列以下各行的公式可以类推。

版权声明：本站内容均来自互联网，仅供演示用，请勿用于商业和其他非法用途。如果侵犯了您的权益请与我们联系QQ:729038198，我们将在24小时内删除。