《水塔流量估计》word版
水塔流量的估计
一.问题的提出
某居民区有一供居民用水的圆柱形水塔,一般可以通过测量其水位来估计其
流量。但面临的困难是,当水塔水位下降到设定的最低水位时,水泵自动启动向
水塔供水,到设定的最高水位时停止供水,这段时间无法测量水塔的水位和水泵
的供水量。通常水泵每天供水一两次,每次约Array 2h(小时)。
水塔是一个高为12.2m,直径为17.4m是正圆
10.8m 柱。按照设计,水塔水位降至约8.2m时,水泵自动
启动,水位升到约为10.8m时水泵停止工作。
表1是某一天的水位测量记录(符号“//”表示水
泵启动),试估计任何时刻(包括水泵正供水时)从水塔流出的水流量及一天的
总用水量。
表1:水位测量记录(时刻:h,水位:cm)
二、问题分析
流量是单位时间流出的水的体积,由于水塔是正圆柱形,横截面积是常数,
在水泵不工作的时段,流量很容易从水位对时间的变化率算出,问题是如何估计
水泵供水时段的流量。
水泵供水时段的流量只能靠供水时段前后的流量拟合得到,作为用于拟合的
原始数据,我们希望水泵不工作的时段流量越准确越好。
这些流量大体可由两种方法计算: 一是直接对表1中的水位用数值微分算
出各时段的流量,用它们拟合其它时刻或连续时间的流量。
二是先用表中数据拟合水位~时间函数,求导数即可得到连续时间的流量。 一般说来数值微分的精度不高,何况测量记录还是不等距的,数值微分的计算尤其麻烦。下面我们用第二种方法处理。
有了任何时刻的流量,就不难计算一天的总用水量。
其实,水泵不工作时段的用水量可以由测量记录直接得到,如表1可知从
t=0到t=8.97(h)水位下降了968 –822=146(cm),乘以水塔的截面积就是这一时段的用水量。这个数值可以用来检查拟合的结果。
三、模型假设
1. 流量只取决于水位差,与水位本身无关。按照Torricelli (托里切利, 1608-1647, 意大利数学家、物理学家、气压计原理发现者)定律从小孔流出的流体的流速正比于水面高度的平方根,题目给出水塔的最低和最高水位分别是
8.2m 和10.8m(设出口的水位为零), 1.151=≈,所以可忽略水位对速度的影响。
2. 根据最低和最高水位分别是8.2m 和10.8m 及表1的水位测量记录,假设水泵第1次供水时段为9t =到11t =,第2次供水时段为20.8t =到23t =。
其中前3个时刻取自实测数据(精确到0.1h ),最后1个时刻来自每次供水约两小时的已知条件(从记录看,每2次供水时段应在有记录的22.96h 之后不久结束)。
3. 水泵工作时单位时间的供水量大致是常数,此常数大于单位时间的平均流量。
4. 流量是对时间的连续函数。
5. 流量与水泵是否工作无关。
6. 由于水塔截面积是常数, ()2
2217.4/2237.787S r m ππ==?≈,为简单起见,计算中将流量定义为单位时间流出的水的高度,即水位对时间变化率的绝对值(水位是下降的),最后给出结果时再乘以S 即可。
即:水位是时间的连续函数 ()h h t =
水位对时间的变化率(流量) ()
dh t h dt
'= 任何时刻的流量: ()()V t h t S
'=
四、模型建立
1.拟合水位~时间函数
从表1 测量记录看,一天有两个供水时段(以下称第1供水时段和第2供水时段)和3个水泵不工作时段(以下称第1用水时段0t =到8.97t =,第2用水时段10.95t =到20.48t =第3用水时段23t =以后)。
对第1、2用水时段的测量数据分别作多项式拟合,得到水位函数 11()h h t = 和22()h h t =。为使拟合曲线比较光滑,多项式次数不要太高,一般用3~6次。由于第3时段只有3个测量记录,无法对这一时段的水位作出比较好的拟合,可采用外推的办法解决。
2. 确定流量~时间函数
对于第1、2用水时段,只需将水位函数(),1,2i i h h t i == 求导数即可,对于两个供水时段的流量,则用供水时段前后(水泵不工作时段)的流量拟合得到,并且将拟合得到的第2供水时段流量外推,将第3用水时段流量包含在第2供水时段内,需要拟合四个流量函数。
3. 一天总用水量的估计
总用水量等于两个水泵不工作时段和两个供水时段用水量之和,它们都可以由流量对时间的积分得到。
t
t
t t V V dt S h dt ''==-?
?
五、模型求解
根据表一,可以对各时段的数据进行拟合。建立时间和水位向量
1212[,,],[,,]t t t h h h ==的函数关系。
为使拟合曲线比较光滑,多项式次数不要太高,一般用3~6次。先用3次函数进行拟合第一时段[0,9]的水位、流量。
设t 、h 为已输入的时刻和水位测量记录(水泵启动的4个时刻不输入):
>> h=[968 948 931 913 898 881 869 852 839 822 0 0 1082 1050 1021 994 965 941 918 892 866 843 822 0 0 1059 1035 1018];
>> t=[0 0.92 1.84 2.95 3.87 4.98 5.90 7.01 7.93 8.97 9.98 10.92 10.95 12.03 12.95 13.88 14.98 15.90 16.83 17.93 19.04 19.96 20.84 22.01 22.96 23.88 24.99 25.91];
>> x1=t(1:10);x2=x1.^2;x3=x1.^3;
>> n=10;m=3;
>> x=[ones(n,1),x1',x2',x3'];
>> [b,bint,r,rint,stats]=regress(h(1:10)',x,0.05) b =
967.7356
-22.1079
1.3586
-0.0785
bint =
966.5367 968.9345
-23.3536 -20.8622
1.0241 1.6932
-0.1030 -0.0541
r =
0.2644
-0.4851
-0.1676
0.6754
0.0260
-0.6328
0.5374
-0.4676
0.3089
-0.0590
rint =
-0.3216 0.8505
-1.5934 0.6232
-1.3634 1.0282
-0.2883 1.6392
-1.2427 1.2947
-1.7032 0.4376
-0.5339 1.6088
-1.5560 0.6207
-0.8696 1.4874
-0.6534 0.5353
stats =
1.0e+004 *
0.0001 2.3222 0.0000 0.0000
>> rcoplot(r,rint)
残差分析可以看到,拟合比较符合,拟合效果较好。经过实验,选取3次较合适。故以下选用三次函数进行拟合。
1、拟合第1用水时段各时刻的流量,可由如下程序代码得到:
c1=polyfit(t(1:10),h(1:10),3); %用3次多项式拟合第1用水时段水位h1=h1(t),c1输出
3次多项式的系数
a1=polyder(c1); %a1输出多项式(系数为c1)导数的系数,
h1=-polyval(a1,t(1:10)); %给出水位变化率h1=h1(t)在 t(1)-t(10)上的离散值,即流量
t1 = 0:0.1:9; %将第一用水时段[0 , 9]细分
h11= -polyval(a1,tp1); %h11输出多项式a1在t11点的函数值(取负后边为正值),即t11
时刻的流量(水位下降的速率)。
用程序进行计算,求得在[0,9]内各时刻的流量值(水位变化率)如下表: i t 0 0.92 1.84 2.95 3.87 4.98 5.90 7.01 7.93 8.97 i h
22.11
19.81
17.91
16.14
15.12
14.42
14.28
14.64
15.38
16.69
第一用水段流量(水位变化率)曲线图
2. 拟合第2用水时段[11,20.8]的水位、流量
c2=polyfit(t(13:23),h(13:23),3); %用3次多项式拟合第2用水时段水位水位h2=h2(t) ,
c2输出3次多项式的系数 a2 = polyder(c2); %a2输出多项式(系数为c2)导数的系数,
h2=-polyval(a2,t(13:23)); %给出水位变化率h2’=h ’2(t) t2 = 11:0.1:20.8; %将第二用水时段[11 , 20.8]细分
h22 = -polyval(a2,t2); %x2输出多项式(a2)在tp2点的函数值(取负后边为正 值),即tp2
时刻的流量(水位下降的速率)。
进行计算,求得在[11,20.8]内各时刻的流量值(水位变化率)如下表:
i t
10.95 12.03 12.95 13.88 14.98 15.90 16.83 17.93 19.04 19.96 20.84 i h
22.11
19.81 17.91 16.14 15.12 14.42 14.28 14.64 15.38 16.69 67.81
第二用水段流量(水位变化率)曲线图
3. 拟合第1供水时段[9,11]的流量
在第1供水时段(t = 9~11)之前(即第1用水时段)和之后(第2用水
时段)各取几点,其流量已经得到,用它们拟合第1供水时段的流量。为使流量
函数在t =9 和t =11连续我们简单地只取4个点,拟合3次多项式(即曲线必过这4个点),实现如下:
x1 = -polyval(a1,[8 9]) %取第1时段在t=8,9的流量
x2=-polyval(a2,[11 12]) %取第2时段在t=11,12的流量
x12 = [x1 x2]; %将四个点合并
c12 = polyfit([8 9 11 12], x12, 3); %拟合3次多项式
t12 =9:0.1:11; %将第一供水时段[11 , 20.8]细分 xx12 = polyval(c12,t12); %xx12输出第一供水时段各时刻的流量
进行计算,求得在[9,11]内各时刻的流量值(水位变化率)如下表:
t9.9810.92
i
h 25.3932.45
i
第一供水段流量(水位变化率)曲线图
4.拟合第2供水时段[20.8,24]的流量
在第2供水时段之前取t =20,20.8两点的水流量,在该时刻之后(第3用水段)仅有3个水位记录,我们用差分得到流量,然后用这4个数值拟合第2供水时段的流量如下:
dt3 = diff(t(26:28)); %最后3个时刻之后的两两之差
dh3 = diff(h (26:28)); %最后3个水位的两两之差
dht3 = -dh3./dt3; % t(22)和t(23)的流量(差商代替微商)
t3 = [20 20.8 t(26) t(27)]; %取第二供水时段前后各两点
x3=[-polyval(a2,t3(1:2)),dht3]; %求得t3各时刻的流量
c3 = polyfit(t3,x3,3); % 拟合3次多项式
t3 = 20.8:0.1:24; %将第二供水时段和第三用水时段细分
xx3=polyval(c3,t3); % xx3输出第2供水时段及第三用水时刻的流量
用程序进行计算,求得在[20.8,24]内各时刻的流量值(水位变化率)如下表:
t22.0122.9623.8824.9925.91
i
h 25.54723.93721.62218.47816.08
i
第二供水段流量(水位变化率)曲线图