2024年3月9日发(作者:中国十大麦克风品牌排行榜)
借助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- ZS)1.5 。
即 Zi=(qi/M/B)+ZS (5)
为了与原假定的Zi加以区别,将(5)式改为
Zi=(qi/M/B),,2/32/3+ZS, 式中B=45m,M=1.70,于是:
2/3Zi= ZS+(0.013072×qi) (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列以下各行的公式可以类推。
6、从单元格T6开始输入“增值修正系数”:0.025、0.025、0.05、0.05、0.05、0.1、0.1、0.15、0.15……,直至最后。这些是经验值,不同类型的水库可能有调整。
7、从第6行开始逐行试算。先给L6单元格付初值,因开始时入库流量小,水位微涨,第一次初估水位增值尽量取小一些,例如取L6=0.08。
① 给“预估水位”单元格M6输入公式=R5+L6,其含义是上一行的“采用水位”加本行“水位增值”。
② 给“库容V”单元格N6输入公式:
=IF(M6>=B$25,C$25+(M6-B$25)*(C$26-C$25),IF(M6>=B$24,C$24+(M6-B$24)*(C$25-C$24),IF(M6>=B$23,C$23+(M6-B$23)*(C$24-C$23),C$22+(M6-B$22)*(C$23-C$22))))
该公式使用了“IF逻辑函数”。使用该函数的规则(称为“语法”)是:在IF的后面必须加括号(),括号里的内容分三段,第一段是逻辑判断式,第二段是当逻辑为真时应当执行的指令,第三段是当逻辑为假时应当执行的指令,各段用“,”隔开。在IF函数之内可以再嵌套IF函数,最多可套7层。嵌套的IF函数应当放在外面一层IF函数括号内第三段,属于逻辑假执行的指令。
该公式是从B列与C列红岩水库水位与库容量关系表中,根据单元格M6的水位值内插单元格N6的值(库容量)。它的含义是:假若单元格M6所代表的“预估水位Z1”值 >=B$25,则在C$25库容量的基础上,按(M6-B$25)水位差的比例加一个库容量增值(M6-B$25)*(C$26-C$25),否则执行下一个IF逻辑函数,再作第二次逻辑判断。
③ 单元格O6代表出流量q,它的公式是=11.111*(N5-N6)-O5+K5+K6,这是根据计算出流量的公式qi=11.111*(Vi-1-Vi)-qi-1+Qi-1+Qi套用的。
④ 单元格P6代表计算水位Z2,它的公式是=$B$29+(0.013072*O6)^(2/3),这是根据计算水位的公式Zi=(qi/M/B)2/3+ZS=(0.013072*qi)2/3+ZS套用的。 其中ZS是溢洪道堰顶高程492.36m。
⑤ 单元格Q6代表计算水位Z2与预估水位Z1的水位差值。它的公式是=P6-M6。
⑥ 单元格R6代表采用水位,它的公式是=ROUND(P6,2),称为ROUND函数,括号内的
P6其含义是对单元格P6的值进行四舍五入,括号内的数据2,是取小数后2位。
⑦ 单元格S6代表水位增值提示△Z2,它的公式是:
=IF(O6<0,0.85*L6,IF(ROUND(P6,2)=ROUND(M6,2),L6,L6+T6*Q6))
单元格S6的作用不直接参与调洪演算,但它根据计算水位Z2与预估水位Z1差值的大小来计算出一个提示值,为下一步试算的水位增值Z1的修改起到参考作用。也就是说可以参考S6的值来修改单元格L6的值。S6使用了IF函数,它的含义是:如果代表出库流量q单元格O6出现负值,则S6应当等于L6乘以0.85,否则执行嵌套其中的第二个IF函数。后者也分为三段:第一段(ROUND(P6,2)=ROUND(M6,2)是“逻辑判别式”,第二段L6是逻辑真的执行指令,第三段是执行逻辑假的指令,L6+T6*Q6就是S6的取值。把这三段连起来说就是:“如果P6单元格的值取2位小数以下4舍5入,等于M6单元格同样处理后的值,那么将L6付值给S6。也就是说认可原来输入到L6的值正确,不需要修改了。否则,S6应当等于原来的L6再加一个修正值,此值等于代表水位误差的单元格Q6乘以代表修正系数的单元格T6”。当S6的公式编辑好后,单击编辑栏ƒχ前面的√符号,S6的值0.008313便出现在该单元格内。此值与L6第一次输入的值0.008不一样,说明L6应当参照S6的值修改为0.008313。
⑧ 将单元格L6的值修改为0.008313后,发现单元格S6仍然等于0.008313。再检查M6=
492.3683,P6= 492.3665,二者取2位小数四舍五入后都等于492.37。检查代表“采用水位Z3”的单元格R6也等于492.37,说明这次试算成功。
8、接着试算第7行:
① 先在单元格L7输入一个初值,因为时隔0.5小时后,入库流量继续在增加,水位将继续上升,水位涨幅也会相应增大,上一个时段的水位增幅是0.008313cm,本时段可适当加大一些,例如取L7=0.015cm。
② 用鼠标按住单元格M6不放,向右拖动到S6,此时M6右边各单元格都变成蓝色,然后点击“复制”图标,再将鼠标按住单元格M7不放,向右拖动到S7,此时M7右边各单元格都变成蓝色,然后点击“粘贴”图标。于是第6行从M6到S6各单元格的公式都复制到第7行M7~S7各单元格中,而且各有关单元格的数值也呈现出来。如果S7的值不等于原置入在单元格L7的数值,那就参照S7显的数值修改L7,直到二者相等为止。
9、试算以下各行:与第7行的试算方法相同,不再赘述。
10、这种试算方法在事先设计好的计算公式下运行。在L~S列各行之间,经过“复制”和“粘贴”,巧妙完成试算任务。试算30个时段的出库流量过程,大约需要一个半小时。
四、通过VBA语言编程来实现全自动化的水库调洪演算
(一)准备工作
1、使用与前面Excel电子表格同一工作簿,在当前工作表Sheet1的情况下,点击该工作表左下角的Sheet2标签,于是工作表Sheet2被打开。按照前面第三项“准备工作”在B2~G34区域内输入相同的内容。见附图(二)
2、合并单元格I2~S3区域。输入“P=0.1﹪红岩水库入、出库流量调洪演算过程”作为该部分表格的标题。
3、在I4~S4各单元格依次输入:序号、时间t、入流Q、水位增值△Z1、预估水位Z1、库容V、出流q、计算水位Z2、采用水位Z3、增值修正系数k、试算次数n共11列的标题。
4、从I5到I35输入1、2、3、……31的序号。
5、从J5到J35依次输入-0.4、0.1、0.6、1.1、……14.6,即每两行之间增加0.5小时。
6、从单元格K5到Q5依次输入1.62、0、492.36、479.57、0、492.36、492.36。其中1.62是基流,492.36是红岩水库溢洪道堰顶高程,497.57是其相应的库容量。这些是洪水开始入
库前有关数据的初始状态。
7、在代表“水位增值修正系数k”的R列,从单元格R6开始,依次输入0.025、0.025、0.05、0.05、0.05、0.1、0.1、0.15、0.15、……0.15,一直到R35。这些是针对红岩水库的经验系数,用于其他水库可能有所变化,估计不加修改地直接用这些数也无妨,可能会增加一定的试算次数。但作为当前计算机的高速运算性能来说,多试算几次甚至几十次,在时间上不会发生质的变化。
(二)编写VBA代码
1、VBA编程代码的功能:什么是VBA?VBA是Visual Basic Application的缩写。意即“可视化的Basic应用程序语言”,它是内置于Office软件中的程序设计语言。具体地说,我们当前使用的VBA不是一般范指的VBA,而是依托在Excel电子表格中的VBA。由于它具有循环结构功能,又有逻辑判断能力,它就超越了单纯使用Excel电子表格的公式和函数的局限性,给一些十分复杂的技术难题找到解决办法。
2、编制红岩水库VBA代码:右击当前工作表左下方的“Sheet2”标签出现一个对话框,点击“查看代码”命令,打开Visual Basic编辑器,它是编辑VBA代码的场所。红岩水库调洪演算的代码就在标有“(通用)”条下面的窗口编写。作为范例,“红岩水库调洪演算”的代码已经编成,展现在代码窗口内。见附图(三)之1.及附图(三)之2。
(三)运行代码
1、在Visual Basic编辑器上面的“菜单”栏点击“运行”按钮,接着点击VBA代码编辑窗口左边的工程资源管理器窗口上方“工程—VBA Project”栏目下面3个图标中的第2个“查看对象”图标,于是Sheet2工作表立即展现在眼前。原来该工作表K6~Q35区域各单元格是空的,是等待试算的内容,现在该区域已全部充满数据。
2、上述区域的试算目标项目:入库流量、出库流量、坝前水位(即采用水位栏)、库容量与Sheet1工作表中相同项目各行数据对照,除个别尾数略有差别外,从总体上看是完全一致的。这就说明上述VBA代码的编辑完全成功,而且这7列30行试算数据,大约完成在1~2秒之间,其运算速度相当惊人。见附图(二)。
(四)解析代码(见【附件】红岩水库调洪演算代码表)
1、编写代码首行以“Private Sub 红岩水库调洪演算()”开头,以“End Sub”结束。Sub代表子程序,也称为过程。Private是关键字,它代表局部过程。因该程序是专为红岩水库设计的,使用Private关键字具有“专用”的含义。如果使用Pablic 关键字就成为公共过程,具有通用的含义。“红岩水库调洪演算”是给该过程取一个名字。()内是参数列表,如果没有参数,可让它空着,但括号必须有。这是死规定。
2、中间共有75行,按顺序编列了行号。行号可有可无,此次使用行号是便于解析各行的作用。
3、1~6行是有关数据类型的声明。①其中第1行是定义所涉及的区域变量名myRange的数据为单精度浮点型Single,即带有小数点的32位数据。其语法形式为“Dim某变量As数据类型”。②第2~3行分别对y和x两个变量定义为整型Integer,即16位整数。③4~6行给3个一维数组F(5 To 21)、G(5 To 21)、J(5 To 35)变量定义为单精度型Single。
4、第7行是对Zyd变量付值。此变量由3个英文字母组合而成,代表溢洪道堰顶高程。这是自取的变量名称,我们习惯于用Z表示水位,也含高程的意思,yd是汉语拼音“堰顶”二字的首字母。红岩水库溢洪道堰顶高程为492.36m。此行的意思是将数据492.36付值给变量Zyd。Let代表付值语句。
5、8~11行代表一个循环控制过程。①第8行以For开头,y =5 To 21意思是y的初值为5,终值是21。最初用5付值给变量y,依次去执行第9行、第10行语句,到第11行Next y完成第一次循环,再回到第8行。Next是“下一个”的意思,此时y=y+1,将原来的y加上1
付值给新y变成6。再去执行第9、第10行语句,到第11行完成第2次循环,再回到第8行,此时y的值又增加了1,变成7,……最后y变成21,当执行到第11行已完成了全部循环过程。此时不再返回第8行,程序进入第12行。②第9行F(y) = Cells(y, "F")将F列5~21行各单元格数据付值给一维变量F(y)对应行。其中Cells(y,“F”)表示F列某行单元格中的数据,当y=9时,就是指F列第9行单元格中的数据,即F9=1.46。Cells就是单元格的意思。③第10行G(y) = Cells(y, "G")将G列5~21行各单元格数据(流量Qt)付值给一维变量G(y)对应行。④8~11行的实质是通过变量y从初值5到终值21的多次循环,分别让两个一维数组变量F(y)和G(y)分别在F列及G列5~21行读取设计洪水过程线上各节点的时间t和流量Q的数据。
6、12~14行:又是一个循环功能控制系统,其实质是用一维数组变量J(y)去读取J列5~35行的时间数据t。
7、15~75行是变量y从初值6到终值35的大循环,其中还嵌套了51~71行变量x从初值1到终值500的小循环。
8、16~46行使用了IF函数,根据J列各行的值从F列和G列所代表的设计洪水过程线列表的时间t和对应的流量Q各行的节点值内插与J列对应的M列的入库流量Q。
① 逻辑判断应先用最大的内插值作比较,查看J列35行的时间t=14.6小时,它位于F列19~20行数据12.47~15.33小时之间,因此程序第16行编写为:“If J(y)>= F(19 )Then”,意思是:“假若一维变量J(y)>=一维变量F(19),那么……”就去执行第17行的指令。假若第16行的判断不成立,那么就转到第18行进行第二次逻辑判断,如果该逻辑判断成立,就去执行第19行的指令。以此类推……,将一维变量J(y)在F列各行由大到小依次比较,直到该逻辑判断为真时,就去执行下一行的指令。随后转到第46行“End IF”表示IF函数执行完毕,转到下面执行第47行的指令。如果直到第42行的判断也不成立,就转到第44行“Else”的意思是“那么就……”去执行下一行第45行的指令。
② 内插的方法是用线性比例内插。现以内插K列12行流量Q为例,由于J列12行的时间为3.1小时,它位于F列12~13行数值2.94~3.56之间,在进入VBA代码第30行“ElseIf
J(y)>= F(12)Then”的逻辑判断为真,于是去执行第31行的指令:
Cells(y, "K") = G(12) + (G(13) - G(12)) * (J(y) - F(12)) / (F(13) - F(12)),其中Cells(y, "K")代表单元格K列y行,此时的y=12。等式右端应当在G(12)其值为208的基础上加一个比例增值(G(13) – G(12)) * (J(y) –
F(12)) / (F(13) – F(12)),其中(G(13) – G(12))代表两行之间流量增值△Q=219-208=11,其对应的时间增值是(F(13) – F(12)),即△t=3.56-2.94=0.72小时,(J(y) – F(12))式中的循环变量y=12,该式表示J列12行的时间比F列12行的时间增加了0.16小时,因此这个比例增值(G(13) – G(12)) * (J(y) – F(12)) / (F(13) – F(12))=(219-208)*(3.1-2.94)/(3.56-2.94)=2.84。于是等式右端的值为208+2.84=210.84,取3位有效数,单元格K列12行的流量应当是211m3/s。
9、 47~50行又是一个If函数,它的作用是为单元格Cells(y, "L")即水位增值△Z1设初值。其中47~48行的意思是“假若y=6,就给单元格L6付值0.005”,否则就执行第49行:
Else:
Cells(y, "L") = 2 * Cells(y - 1, "L") - Cells(y - 2, "L"),此行的意思是:“当y>6,把单元格Cells(y, "L")的上一行Cells(y - 1, "L")作基数,再加上Cells(y, "L")的上一行与上两行之差。
10、第51~71行是在第15~75行变量y的大循环中嵌套的一个小循环。其中又包含了两个If函数段,第53~61行,及64~70行。
① 第51行:For x = 1 To 500,给循环变量x定义的初值为1,终值为500,其步长为1,即每循环一次,x的值增加1。为什么要把终值定那么大?是否一定要达到500次才能结束循环?500次是主观估计数,目的是定大一点更保险,由于试算过程比较复杂,如果终值定小了,有可能循环到了终点而试算还没有完成。根据已有实践经验,本人在此前编的程序试算次数有达到150次以上的。但并非一定要达到500次才能结束循环,只要试算成功就通过第67行跳出循环。
② 第52行是给代表“预估水位Z1”的M列y行单元格付值。Cells(y, "M") = Cells(y - 1,
"Q") + Cells(y, "L")意思是将上一行代表“采用水位Z3”的单元格Q列y-1行的值,加代表“水位增值△Z1”的单元格L列y行的和,付值给Cells(y, "M")。
③ 第53~61行是一组If函数。其目的是根据代表“预估水位Z1”M列y行的值,从“水位~库容量关系”表B列和C列对应数据中内插N列y行的库容量。与B列有关的范围在溢洪道堰顶到坝顶高程所在水位区间,即492~496m之间,单元格B22~B26范围内。第53行逻辑判断应当从B列25~26行开始比较,因此第53行为 If Cells(y, "M") >= Cells(25, "B") Then,意即“假若单元格M列y行的值>=单元格B列25行,那么……”,接下来就是去执行第54行:Cells(y, "N")
= Cells(25, "C") + (Cells(y, "M") - Cells(25, "B")) * (Cells(26, C) - Cells(25, "C")), 该行等式右端以单元格C25为基数,加一个比例增值(Cells(y, "M") - Cells(25, "B")) * (Cells(26, C) - Cells(25, "C")),此增值省略了“/(Cells(26, "B")- Cells(25, "B"))这一部分,因为(Cells(26, "B")- Cells(25, "B")=1,以1作除数故可省略。
④ 第62行是给单元格Cells(y, "O")付值。按照计算出库流量的公式qi=11.111*(Vi-1-Vi)-qi-1+Qi-1+Qi编写为:Cells(y, "O") = 11.111 * (Cells(y - 1, "N") - Cells(y, "N")) - Cells(y - 1, "O") +
Cells(y - 1, "K") + Cells(y, "K")。
⑤ 第63行是给代表“计算水位Z2”的Cells(y, "P")付值。按照计算水位的公式为:Zi=(qi/M/B)2/3+ZS=(0.013072qi)2/3+ZS,编写为:Cells(y, "P") = Zyd + (0.013072 * Cells(y, "O")) ^ (2
/ 3)。
⑥ 第64~70行又是一个If函数程序组,它是本程序最核心的试算部分。第64行If Cells(y,
"O") < 0 Then,单元格O列y行代表“出库流量q”,它不应当< 0,如果出现此情况,那是代表“水位增值Z1” Cells(y, "L")单元格的初值假设过大。因此在第65行令Cells(y, "L") = 0.85 * Cells(y, "L"),意即将原来的Cells(y, "L")乘以0.85的系数,再付值给Cells(y, "L")。
⑦ 第66行 ElseIf Round(Cells(y, "P"), 2) = Round(Cells(y, "M"), 2) Then意思是“除非假若代表计算水位Z2的单元格Cells(y, "P")取2位小数四舍五入的值等于代表水位预估值Z1的单元格Cells(y, "M")” 取2位小数四舍五入的值时,那就去执行下一行”
⑧ 第67行GoTo Outforx。此行的意思是:“跳出循环转到以行标识符为Outforx的那一行”。因为第66行的逻辑判断如果为真,那就表明第64~70行的试算过程已经取得成功,就可以结束循环了。行标识符可以随意写几个英文字母,这里用Outforx表示“跳出x变量循环的意思”。
⑨ 如果第66行的逻辑判断不成立,它就不执行67行,转到第68行Else。意即“那么就……”去执行第69行。
Cells(y, "L") = Cells(y, "L") + (Cells(y, "P") - Cells(y, "M")) * Cells(y,
"R"),此行的意思是将原来的代表“水位增值△Z1”的单元格Cells(y, "L")增加一个修正值(Cells(y, "P") - Cells(y, "M")) * Cells(y, "R"),该数值等于“计算水位Z2”减“预估水位Z1”之差乘以一个修正系数Cells(y, "R"),再付值给新的Cells(y, "L")。于是转到70行End If,结束第64~70的If函数组。下到第71行Next x,意即下一个x,让程序再转回到第51行重新再一次循环。
11、 第72行Outforx:,,它没有什么意义,只是简单的行标识符。其作用是引导跳出循环的程序转到该行。
12、 第73行Cells(y, "Q") = Round(Cells(y, "P"), 2)此行的意思是将最终的计算水位Z2经过Round函数处理后,取小数2位以下四舍五入,付值给“采用水位Z3” Cells(y, "Q")。
13、 第74行Cells(y, "S") = x此行是把循环变量x付值给“计算次数n”栏单元格Cells(y,
"S")。x是在第51行设置的循环变量,当第66行的逻辑判断条件为真时跳出循环,此时的x值不再增加,保留了原来的循环次数。因此将x付值给单元格Cells(y, "S")也就代表在Sheet2工作表上y行的调洪演算试算次数。
14、第75行Next y。循环变量y是在第15行设置的,它的初值是6,终值是35。y在第15~75行之间循环,每循环一次y的值增加1,当y=35后,从第15行运行到75行y的值不再增加,脱离循环,进入下面一行。
15、VBA代码的最后一行End Sub表示该子程序的结束。这一行必须有,也是死规定。
(五) 讨论
本代码是专为红岩水库泄洪流量过程调洪演算编写的,它没有涉及到预留防洪库容的问题,它的溢洪道没有安装闸门,不属于有闸调洪问题,因而它不具备多种调洪功能。今后作者准备推出其他功能调洪演算的例子。目前以它作为范例推广应用到其他同类型水库VBA代码编程。现在要讨论的问题是如何将本代码移植到别的水库?一句话:可以照套不可照搬。由于各水库的流域持征值不同,水库大坝及溢洪道各部位结构尺寸不同,水位~库容量关系不同,凡是有差别的部分必须用实际数据来代替。例如VBA代码第7行就必须用具体水库的溢洪道堰顶高程给变量Zyd付值,而不应当照抄492.36。
1、按照前面所做的准备工作,一定要落实到具体水库的实际情况,包括:①水位~库容量关系表,②大坝的坝顶高程,③溢洪道净宽、堰顶高程、堰型、第二流量系数、④设计洪水过程线的水位~流量关系表,⑤在Sheet2工作表B31~G34区域内有关的计算公式等。
2、关于Sheet2工作表J列“时间t栏目”的时段长,红岩水库用0.5小时为一时段。对于一般中小型水库同样可用0.5小时为一时段。红岩水库的时间t采用-0.4、0.1、0.6、1.1、1.6、……的原因:一是与设计洪水过程线入库洪峰所在时间相对应,二是为了便于与毛坝水库1995年工程整治中红岩水库调洪演算成果相比较。本人建议采用0、0.5、1.0、1.5、2.0、……的时间段为好,即或与设计洪水过程线入库洪峰所在时间不相对应,对出库流量的调洪演算成果的影响也不大。
3、关于Sheet2工作表R列“水位增值修正系数k”如何取值的问题:把红岩水库用的k值照套到别的水库是否可以的问题,对于不同水库的k值可能有变化,对此缺乏经验,可进一步探讨。但这不是主要问题,即使用得不对,可能增加一些试算次数,由于当前的计算机运转迅速,可能不会有大的影响。
(六) 注意
编写VBA代码 是一项复杂又细致的工作。有了范例就有了模仿的蓝本 ,给读者尝试编程提高了信心。但在实践过程中总会遇到这样那样的问题,事所难免。在此要提醒读者注意的是:第一要有信心,第二要谨慎。有信心才不怕困难与挫折,谨慎是克服困难的法宝。作者曾多次遇到过这样的问题:当VBA代码编写完毕交付运行时,常常遇到一些莫明其妙的“警告!”,不知从何下手纠正错误。例如“使用了无效的变量或对象域”,后来发现在R列没有事先输入“增值修正系数k”,使得第69行代码执行出错。还有一次警告出现“编辑错误”,“子程序或函数未定义”,“缺少语句结束”,“应用程序定义或对象定义错误”等,最后发现已经编好的代码中两处变成蓝色的小块,属于Cells拼写错误,打了重码,写成“Cellls”了!改正后就成功了。从这些例子来看既不是编程的思路有问题,也不是VBA语法错误,一个英文字母打错了,就全盘皆输。最后送读者一句话:“战略上自信,战术上谨慎”。祝你成功!
附图(一)
利用Excel电子表格的公式和函数实现半自动化的水库调洪演算
附图(二)
通过VBA语言编程来实现红岩水库全自动化的调洪演算
附图(三)之1.
红岩水库调洪演算的VBA代码编写窗口
附图(三)之2.
红岩水库调洪演算的VBA代码编写窗口
【附件】 红岩水库调洪演算VBA代码表
Private Sub 红岩水库调洪演算()
1 Dim myRange As Single
2 Dim y As Integer
3 Dim x As Integer
4 Dim F(5 To 21) As Single
5 Dim G(5 To 21) As Single
6 Dim J(5 To 35) As Single
7 Let Zyd = 492.36
8 For y = 5 To 21
9 F(y) = Cells(y, "F")
10 G(y) = Cells(y, "G")
11 Next y
12 For y = 5 To 35
13 J(y) = Cells(y, "J")
14 Next y
15 For y = 6 To 35
16 If J(y) >= F(19) Then
17 Cells(y, "K") = G(19) + (G(20) - G(19)) * (J(y) - F(19)) / (F(20) - F(19))
18 ElseIf J(y) >= F(18) Then
19 Cells(y, "K") = G(18) + (G(19) - G(18)) * (J(y) - F(18)) / (F(19) - F(18))
20 ElseIf J(y) >= F(17) Then
21 Cells(y, "K") = G(17) + (G(18) - G(17)) * (J(y) - F(17)) / (F(18) - F(17))
22 ElseIf J(y) >= F(16) Then
23 Cells(y, "K") = G(16) + (G(17) - G(16)) * (J(y) - F(16)) / (F(17) - F(16))
24 ElseIf J(y) >= F(15) Then
25 Cells(y, "K") = G(15) + (G(16) - G(15)) * (J(y) - F(15)) / (F(16) - F(15))
26 ElseIf J(y) >= F(14) Then
27 Cells(y, "K") = G(14) + (G(15) - G(14)) * (J(y) - F(14)) / (F(15) - F(14))
28 ElseIf J(y) >= F(13) Then
29 Cells(y, "K") = G(13) + (G(14) - G(13)) * (J(y) - F(13)) / (F(14) - F(13))
30 ElseIf J(y) >= F(12) Then
31 Cells(y, "K") = G(12) + (G(13) - G(12)) * (J(y) - F(12)) / (F(13) - F(12))
32 ElseIf J(y) >= F(11) Then
33 Cells(y, "K") = G(11) + (G(12) - G(11)) * (J(y) - F(11)) / (F(12) - F(11))
34 ElseIf J(y) >= F(10) Then
35 Cells(y, "K") = G(10) + (G(11) - G(10)) * (J(y) - F(10)) / (F(11) - F(10))
36 ElseIf J(y) >= F(9) Then
37 Cells(y, "K") = G(9) + (G(10) - G(9)) * (J(y) - F(9)) / (F(10) - F(9))
38 ElseIf J(y) >= F(8) Then
39 Cells(y, "K") = G(8) + (G(9) - G(8)) * (J(y) - F(8)) / (F(9) - F(8))
40 ElseIf J(y) >= F(7) Then
41 Cells(y, "K") = G(7) + (G(8) - G(7)) * (J(y) - F(7)) / (F(8) - F(7))
42 ElseIf J(y) >= F(6) Then
43 Cells(y, "K") = G(6) + (G(7) - G(6)) * (J(y) - F(6)) / (F(7) - F(6))
44 Else
45 Cells(y, "K") = G(5) + (G(6) - G(5)) * (J(y) - F(5)) / (F(6) - F(5))
46 End If
47 If y = 6 Then
48 Cells(y, "L") = 0.005
49 Else: Cells(y, "L") = 2 * Cells(y - 1, "L") - Cells(y - 2, "L")
50 End If
51 For x = 1 To 500
52 Cells(y, "M") = Cells(y - 1, "Q") + Cells(y, "L")
53 If Cells(y, "M") >= Cells(25, "B") Then
54 Cells(y, "N") = Cells(25, "C") + (Cells(y, "M") - Cells(25, "B")) * (Cells(26, C) - Cells(25, "C"))
55 ElseIf Cells(y, "M") >= Cells(24, "B") Then
56 Cells(y, "N") = Cells(24, "C") + (Cells(y, "M") - Cells(24, "B")) * (Cells(25, "C") - Cells(24, "C"))
57 ElseIf Cells(y, "M") >= Cells(23, "B") Then
58 Cells(y, "N") = Cells(23, "C") + (Cells(y, "M") - Cells(23, "B")) * (Cells(24, "C") - Cells(23, "C"))
59 Else
60 Cells(y, "N") = Cells(22, "C") + (Cells(y, "M") - Cells(22, "B")) * (Cells(23, "C") - Cells(22, "C"))
61 End If
62 Cells(y, "O") = 11.111 * (Cells(y - 1, "N") - Cells(y, "N")) - Cells(y - 1, "O") + Cells(y - 1, "K") + Cells(y, "K")
63 Cells(y, "P") = Zyd + (0.013072 * Cells(y, "O")) ^ (2 / 3)
64 If Cells(y, "O") < 0 Then
65 Cells(y, "L") = 0.85 * Cells(y, "L")
66 ElseIf Round(Cells(y, "P"), 2) = Round(Cells(y, "M"), 2) Then
67 GoTo Outforx
68 Else
69 Cells(y, "L") = Cells(y, "L") + (Cells(y, "P") - Cells(y, "M")) * Cells(y, "R")
70 End If
71 Next x
72 Outforx:
73 Cells(y, "Q") = Round(Cells(y, "P"), 2)
74 Cells(y, "S") = x
75 Next y
End Sub
发布者:admin,转转请注明出处:http://www.yc00.com/num/1709931000a1672818.html
评论列表(0条)