估计水塔的水流量
摘要:建立模型提供水流量及一天的总用水量的较为准确的估计值。本文采用多项式拟合的方法，利用数学软件MATLAB对水塔流量进行计算，两个供水时段水泵功率相当，可见计算结果与实际相吻合。
关键词：建模流量拟合
模型假设
(1)影响水从水塔中流出的流量的唯一因素是公众对水的传统要求。因为表给出的数据没有提及任何其他的影响因素，我们假定所给数据反映了有代表性的一天，而不包括任何特殊情况，如自然灾害、火灾、水塔溢水、水塔漏水等对水的特殊要求。
(2)流量只取决于水位差。水塔中的水位不影响水流量的大小，气候条件、温度变化等也不影响水流量．因为物理学的Torricelli定律指出：水塔的最大水流量正比于水位高度
的平方根，由表中数据有，说明最高水位和最低水位的两个流量几乎相等，
可忽略水位对流速的影响。
(3)水泵工作起止时间由它的水位决定，每次充水时间大约为两个小时．水泵工作性能效率总是一定的，没有工作时需维修、使用次数多影响使用效率问题，水泵充水量远大于水塔水流量。
(4)假设流量是对时间的连续函数。
(5)水塔的水流量与水泵是否工作的状态无关，并不因水泵工作而增加或减少水流量的大小。
(6)定义流量为单位时间流出的水的高度，即水位下降对时间变化率的绝对值，由于水塔截面积是常数s，最后给出结果时乘以s即可。任何时刻的流量：v(t)=-h’(t)*s
流量(水位变化率)估计方法
(1)拟合水位-时间函数，从表中的测量记录看，一天有两次供水时段和三次未供水时段，分别对第1，2未供水时段的测量数据直接作多项式函数拟合，可得到水位函数，根据多项式拟合的特点，为使拟合曲线比较光滑，此处拟合多项式的次数不宜过高，一般为3～6次。对第3未供水时段来说，3个测量记录的数据过少无法得到很好的拟合，故采用外推的方法解决。
T1=[0.000.921.842.953.874.985.907.017.938.97];
H1=[968948931913901881869849839822];
p1=polyfit(T1,H1,5)%用三次多项式拟合第1时段水位，p1输出多项式的系数
a1=polyder(p1)%输出多项式(系数为p1)导数的系数
tp1=0:0.01:8.97;
h1=-polyval(a1,tp1);%输出多项式(系数为a1)在tp1点的函数值(取负后变为正值)，即tp1时刻的流量
hh1=-polyval(a1,T1);
plot(t1,hh1,'s',tp1,h1,'k')
执行命令后得到p1=-0.00540.1380-1.29455.5908-26.7446968.1889
a1=-0.02710.5518-3.883611.1817-26.7446
第一未供水时间段流量图如下所示，
同时得到如下多项式：
H1=-0.0054*t^5+0.1380*t^4-1.2945*t^3+5.5908*t^2-26.7446*t+968.1889
h1=-0.0271*t^4+0.5518*t^3-3.8836*t^2+11.1817*t-26.7446
第二未供水时间段流量，类似可得，
(2)确定流量-时间函数，对于第1，2未供水时段只需将水位-时间函数求导数即可，而对于两个供水时段的流量，则用供水时段前后期(水泵不工作时段)的流量拟合得到，为使流量函数在t=9和t=11连续，只取4个点，用三次多项式拟合得到第一供水时段的时间-流量图。在第二供水时间段之前取t=19.96,20.84两点的流量，用第3未供水时段的3个记录做差分得到两个流量数据，然后用这4个数据点做三次多项式拟合得到第二供水时段与第3未供水时段的时间-流量图。
jt1=-polyval(a1,[7.938.97]);%取第1时段在t=7.93,8.97的流量
jt2=-polyval(a2,[10.9512.03]);%取第2时段在t=10.95,12.03的流量
djt=[7.938.9710.9512.03];
djh=[jt1jt2];
pjt=polyfit(djt,djh,5)%4个记录拟合三次多项式
tpj=7.93:0.01:12.03;
hj=polyval(pjt,tpj);TT=[7.938.9710.9512.03];
plot(TT,djh,'o',tpj,hj,'b')
执行命令后得到：pjt=-0.00500.1530-1.53245.265700
即得到一个五次多项式，djh=-0.0050*t^5+0.1530*t^4-1.5324*t^3+5.2657*t^2
第一供水时段的时间-流量图第二供水时段的时间-流量图
(3)一天总用水量的估计，总用水量等于两个供水时段和两个未供水时段的用水量之和，它们可由流量对时间的积分得到。
第一用水时段的用水量，其中h1通过梯
形公式计算，第二用水时段的用水量及第一，二供水时段的用水量同理可得。
HH1=0.01*trapz(h1)%第一用水时段的用水量(仍按高度计算)，0.01为步长
HH2=0.01*trapz(h2)%第二用水时段的用水量(仍按高度计算)，0.01为步长
HH3=0.01*trapz(hj)%第一供水时段的用水量
HH4=0.01*trapz(h3)%第二供水时段的用水量
s=pi*[((57/3.281)/2)^2]*10000;
V=(HH1+HH2+HH3+HH4)*s%一天的总用水量(立方厘米)
计算结果
五次多项拟合的时间-流量(以高度计)曲线图
从算法设计和分析可知,计算结果与各时段所用的拟合多项式的次数有关。下表给出的是对第1,2未供水时段分别用5次多项式拟合后得到的用水量结果。
时段
第1未
供水段
第2未
供水段
第1
供水段
第2
供水段
全天用水
用水高度
145.93
260.05
41.67
78.726
526.376
全天的用水量为V=1.2478e+009立方厘米
流量及总用水量的检验
计算出各时刻的流量可用水位记录的数值微分来检验，各时段的用水高度可以用实际记录的水位下降高度来检验。例如，算得第1未供水段的用水量高度是145.93cm，而实际记录的水位下降高度为968-822=146cm，两者是吻合的；同样地，算得第2未供水段的用水量高度是260.05cm，而实际记录的水位下降高度为1082-822=260cm，两者也是吻合的。供水时段的用水量加上水位上升值260是该时段泵入的水量，除以时段长度得到水泵的功率(单位时间泵入的水量)，而两个供水时段水泵的功率应大致相等。第1、2时段水泵的功率可以计算如下：
P1=[V1+V11]/2%第1供水时段水泵的功率
P2=(s*HH4+V22)/2.12%第2供水时段水泵的功率
P1=3.4600e+008，P2=3.7874e+008
模型分析
优点：1.模型可以用于一个标准水箱的小镇使用，容易推广；
2.模型用到的只是简单易懂，模型容易完成；
3.模型不仅提供了水流量及一天用水量的较为准确的估计，还可以估计任何时刻的水流量，包括水泵工作的水流量。
缺点：
无法准确估计结果的误差。
参考文献：
(l)张德丰，《MATLAB数值计算方法》，北京：机械工业出版社，，p.162-168
附录：
表：某小镇某天的水塔水位(1m=3.281英尺)
时间(秒)
水位(英尺)
时间(秒)
水位(英尺)
时间(秒)
水位(英尺)
0
31.75
35932
水泵工作
68535
28.42
3316
31.10
39332
水泵工作
71854
27.67
6635
30.54
39435
35.50
75021
26.97
10619
29.94
43318
34.45
79154
水泵工作
13937
29.55
46636
33.50
82649
水泵工作
17921
28.92
49953
32.67
85968
34.75
21240
28.50
53936
31.56
89953
33.89
25223
27.87
57254
30.81
93270
33.40
28543
27.52
60574
30.12
32284
26.97
64554
29.27
时间(小时)
水位(厘米)
时间(小时)
水位(厘米)
时间(小时)
水位(厘米)
0.00
968
10.95
1082
20.84
822
0.92
948
12.03
1050
23.88
1059
1.84
931
12.95
1021
24.99
1033
2.95
913
13.88
996
25.91
1018
3.87
901
14.98
962
4.98
881
15.90
939
5.90
869
16.83
918
7.01
849
17.93
892
7.93
839
19.04
866
8.97
822
19.96
843
详细程序：
T1=[0.000.921.842.953.874.985.907.017.938.97];
H1=[968948931913901881869849839822];
p1=polyfit(T1,H1,5);
a1=polyder(p1);tp1=0:0.01:8.97;
h1=-polyval(a1,tp1); hh1=-polyval(a1,T1);
T2=[10.9512.0312.9513.8814.9815.9016.8317.9319.0419.9620.84];
H2=[108210501021996962939918892866843822];
p2=polyfit(T2,H2,5);
a2=polyder(p2);
tp2=10.95:0.01:20.84;
h2=-polyval(a2,tp2);
hh2=-polyval(a2,T2);
jt1=-polyval(a1,[7.938.97]);
jt2=-polyval(a2,[10.9512.03]);
djt=[7.938.9710.9512.03];
djh=[jt1jt2];
pjt=polyfit(djt,djh,5);
tpj=8.97:0.01:10.95;
hj=polyval(pjt,tpj);TT=[7.938.9710.9512.03];
TT3=[19.9620.8423.8824.99];
T3=[23.8824.9925.91];
H3=[105910331018];
dt3=diff(T3);%最后三个时刻的两两之差
dh3=diff(H3);%最后三个水位的两两之差
dth3=-dh3/dt3;%T(22)和T(23)的流量
ht3=[-polyval(a2,TT3,dth3)];%TT3各时刻的流量
p3=polyfit(TT3,ht3,5);%拟合三次多项式
tp3=19.96:0.01:24;
h3=polyval(p3,tp3);%输出第二供水时段外推至T=24各时刻的流量
hh3=polyval(p3,TT3);plot(T1,hh1,'s',tp1,h1,'k',TT,djh,'o',tpj,hj,'b',T2,hh2,'*',tp2,h2,'m',TT3,hh3,'p',tp3,h3,'k')
xlabel('时间-h')
ylabel('流速cm/h')
title('5次拟合的流量曲线')
legend('第一未供水时刻水流速','第一未供水时段','第一供水时刻水流速','第一供水时段','第二未供水时刻水流速','第二未供水时段','第二供水时刻水流速','第二供水时段')
HH1=0.01*trapz(h1)
HH2=0.01*trapz(h2)
HH3=0.01*trapz(hj)
HH4=0.01*trapz(h3)
s=pi*[((57/3.281)/2)^2]*10000;
>>V=(HH1+HH2+HH3+HH4)*s
V1=HH1*s
V2=HH2*s
V11=(968-822)*s
V22=(1082-822)*s
P1=[V1+V11]/2%第1供水时段水泵的功率
P2=(s*HH4+V22)/2.12%第2供水时段水泵的功率
