水塔水流量估计.docx
《水塔水流量估计.docx》由会员分享，可在线阅读，更多相关《水塔水流量估计.docx（16页珍藏版）》请在冰豆网上搜索。
水塔水流量估计
实验水塔水流量的估计
实验目的
本次实验的主要目的是让学生会用数学软件进行插值计算并解决一些具体的实际问题。
介绍一些经典的插值方法，包括拉格朗日插值法、埃尔米特插值法、分段插值法、三次样条插值法等等。
实验内容
1实验问题
美国某州的各用水管理机构要求各社区提供以每小时多少加仑计的用水率以及每天所用的总水量。
许多社区没有测量流入或流出水塔的水量装置,他们只能代之以每小时测量水塔中的水位,其误差不超过5%。
更重要的是,当水塔中的水位下降到最低水位L时水泵就启动向水塔输水直到最高水位H,期间不能测量水泵的供水量。
因此,当水泵正在输水时不容易建立水塔中水位和用水量之间的关系。
水泵每天输水一次或两次,每次约二小时。
试估计任何时刻（包括水泵正在输水时间）从水塔流出的水流量f（t）,并估计一天的总用水量。
已知该水塔是一个高为40英尺（ft），直径为57英尺（ft）的正圆柱，表12.1给出了某个小镇一天水塔水位的真实数据，水位降至约27.00ft水泵开始工作,水位升到35.50ft停止工作。
（注：
1英尺（ft）=0.3024米（m））
表12-1某小镇某天水塔水位
时间/s
水位/0.01ft
时间/s
水位/0.01ft
0
3175
46636
3350
3316
3110
49953
3260
6635
3054
53936
3167
10619
2994
57254
3087
13937
2947
60574
3012
17921
2892
64554
2927
21240
2850
68535
2842
25223
2795
71854
2767
28543
2752
75021
2697
32284
2697
79254
水泵开动
35932
水泵开动
82649
水泵开动
39332
水泵开动
85968
3475
39435
3550
89953
3397
43318
3445
93270
3340
2问题分析
流量是单位时间内流出水的体积,由于水塔是正圆柱形,横截面积是常数,所以在水泵不工作时段,流量很容易根据水位相对时间的变化率算出。
问题的难点在于如何估计水泵供水时段的流量。
水泵供水时段的流量只能靠供水时段前后的流量经插值或拟合得到。
作为用于插值或拟合的原始数据,我们希望水泵不工作时段的流量越准确越好。
这些流量大体上可由两种方法计算,一是直接对表12.1中的水量用数值微分算出各时段的流量,用它们拟合其它时刻或连续时间的流量;二是先用表中数据拟合水位—时间函数,求导数即可得到连续时间的流量。
有了任何时刻的流量,就不难计算一天的总用水量。
其实,水泵不工作时段的用水量可以由测量记录直接得到,由表12.1中下降水位乘以水塔的截面积就是这一时段的量这个数值可以用来检验数据插值或拟合的结果。
模型建立
为了表示方便,我们将问题中所给表12-1中的数据全部化为国际标准单位（表12-5）,时间用小时（h）,高度用米（m）:
表12-5一天内水塔水位记录
时间（h）
水位（m）
时间（h）
水位（m）
0
9.68
12.95
10.21
0.92
9.45
13.88
9.94
1.84
9.31
14.98
9.65
2.95
9.13
15.90
9.41
3.87
8.98
16.83
9.18
4.98
8.81
17.93
8.92
5.90
8.69
19.04
8.66
7.00
8.52
19.96
8.43
7.93
8.39
20.84
8.22
8.97
8.22
22.02
水泵开动
9.98
水泵开动
22.96
水泵开动
10.93
水泵开动
23.88
10.59
10.95
10.82
24.99
10.35
12.03
10.50
25.91
10.18
模型假设
（1）流量只取决于水位差,与水位本身无关，故由物理学中Torricelli定律：
从小孔流出的液体的流速正比于水面高度的平方根。
题目给出水塔的最低和最高水位分别是8.1648m（27×0.3024）和10.7352m（35.50×0.3024）（设出口的水位为零），因为sqrt（10.7352/8.1648）≈1.1467，约为1,所以可忽略水位对流速的影响。
（2）将流量看作是时间的连续函数，为计算简单,不妨将流量定义成单位时间流出水的高度,即水位对时间变化率的绝对值（水位是下降的）,水塔截面积为S=（57*0.3048）^2*π/4=237.8（m^2）.得到结果后乘以S即可。
流量估计方法
首先依照表12-5所给数据,用MATLAB作出时间—水位散点图12.7。
Fig.12.7
下面我们来计算水箱流量与时间的关系。
根据表12-5中数据散点图12.7，一种简单的处理方法为先将表12-5中的数据分为三段，然后对每一段的数据做如下处理：
设某段数据为
，相邻数据中点的平均流速用下面的公式（流速=（左端点的水位-右端点的水位）/区间长度）：
，
每段数据首尾点的流速用下面的公式计算：
，
，
用以上公式求得时间与流速之间的数据表12-6如下:
表12-6
时间（h）
流速（cm/h）
时间（h）
流速（cm/h）
0
29.89
12.49
31.52
0.46
21.74
13.42
29.03
1.38
18.48
14.43
26.36
2.395
16.22
15.44
26.09
3.41
16.30
16.37
24.73
4.425
15.32
17.38
23.64
5.44
13.04
18.49
23.42
6.45
15.45
19.50
25.00
7.465
13.98
20.40
23.86
8.45
16.35
20.84
22.17
8.97
19.29
22.02
水泵开动
9.98
水泵开动
22.96
水泵开动
10.93
水泵开动
23.88
27.09
10.95
33.50
24.43
21.62
11.49
29.63
25.45
18.48
25.91
13.30
由表12-6作出时间-流速散点图如下:
Fig.12.8
（1）插值法
由表12-6,对水泵不工作时段1,2采取插值方法,可以得到任意时刻的流速,从而可以知道任意时刻的流量.我们分别采取拉格朗日插值法,分段线性插值法及三次样条插值法；对于水泵工作时段1应用前后时期的流速进行插值，由于最后一段水泵不工作时段数据太少，我们将它与水泵工作时段2合并一同进行插值处理（该段简称混合时段）。
我们总共需要对四段数据（第1，2未供水时段，第1供水时段，混合时段）进行插值处理,下面以第1未供水时段数据为例分别用三种方法算出流量函数和用水量（用水高度）。
下面是实现该过程的MATLAB程序。
t=[0,0.46,1.38,2.395,3.41,4.425,5.44,6.45,7.465,8.45,8.97];
v=[29.89,21.74,18.48,16.22,16.30,15.32,13.04,15.45,13.98,16.35,19.27];
t0=0:
0.1:
8.97;
lglr=lglrcz（t,v,t0）;/*注:
lglrcz为一函数,程序同lglrcz.m*/
lglrjf=0.1*trapz（lglr）
fdxx=interp1（t,v,t0）;
fdxxjf=0.1*trapz（fdxx）
scyt=interp1（t,v,t0,'spline'）;
sancytjf=0.1*trapz（scyt）
plot（t,v,'*',t0,lglr,'r',t0,fdxx,'g',t0,scyt,'b'）
gtext（'lglr'）
gtext（'fdxx'）
gtext（'scyt'）
运行结果为lglrjf=145.6231
fdxxjf=147.1430
sancytjf=145.6870
图12.9是对第1未供水段数据用三种不同方法得到的插值函数图，图中曲线lglr、fdxx和scyt分别表示用拉格朗日插值法,分段线性插值法及三次样条插值法得到的曲线。
Fig.12.9
由表12-5知，第1未供水时段的总用水高度为146（=968-822），可见上述三种插值方法计算的结果与实际值（146）相比都比较接近。
考虑到三次样条插值方法具有更加良好的性质，建议采取该方法。
其它三段的处理方法与第1未供水时段的处理方法类似，这里不再详细叙述，只给出数值结果和函数图象，图中曲线标记同Fig.12.9。
表12-7各时段及一天的总用水量（用水高度）
第1
未供水段
第2
未供水段
第1
供水段
混合
时段
全天
拉格朗日
插值法
145.6231
258.8664
54.0689
92.1337
550.6921
分段线性
插值法
147.1430
258.9697
49.6051
76.4688
532.1866
三次样条
插值法
145.6870
258.6547
53.3334
81.7699
539.4450
Fig.12.10（第一供水段时间-流速示意图）
Fig.12.11（第2未供水段时间-流速示意图）
Fig.12.12（混合时段段时间-流速示意图）
下图12.13是用分段线性及三次样条插值方法得到的整个过程的时间-流速函数示意图。
Fig.12.13
下表是对一天中任取的4个时刻分别用三种方法得到的水塔水流量近似值。
表12-8
6.88
10.88
15.88
22.88
拉格朗日
15.98266712348514
33.74260090853463
25.56622418180477
34.70996210551694
分段插值
14.82724137931034
32.99762626262627
25.44655913978495
25.4715789473684
三次样条
15.05278581965820
33.70895535953259
25.54908920557423
29.41733175863551
（2）拟合法
1）拟合水位--时间函数
从表12-5中测量记录看，一天有两次供水时段和三次未供水时段，分别对第1,2未供水时段的测量数据直接作多项式拟合，可得到水位函数（注意，根据多项式拟合的特点，此处拟合多项式的次数不宜过高，一般以3-6次为宜）。
对第3未供水时段来说,数据过少不能得到很好的拟和。
设t,h分别为已输入的时刻和水位测量记录（由表12-5提供，水泵启动的4个时刻不输入），这样第1未供水时段各时刻的水位可由如下MATLAB程序完成：
t=[0,0.92,1.84,2.95,3.87,4.98,5.90,7.00,7.93,8.97,10.95,
12.03,12.95,13.88,14.98,15.90,16.83,17.93,19.04,19.96,
20.84,23.88,24.99,25.66]
h=[9.68,9.48,9.31,9.13,8.98,8.81,8.69,8.52,8.39,8.22,10.82,
10.50,10.21,9.94,9.65,9.41,9.18,8.92,8.66,8.43,8.22,10.59,
10.35,10.18];
c1=polyfit（t（1:
10）,h（1:
10）,3）;
tp1=0:
0.1:
8.9;
x1=polyval（c1,tp1）;
plot（tp1,x1）
下图给出的是第1未供水时段的时间-水位拟合函数图形。
Fig.12.14
变量x1中存放了以0.1为步长算出的各个时刻的水位高度。
同样地，第2未供水时段时间-水位图可由如下MATLAB程序完成，读者可自己上机运行查看。
c2=polyfit（t（11:
21）,h（11:
21）,3）;
tp2=10.9:
0.1:
20.9;
x2=-polyval（c2,tp2）;
plot（tp2,x2）
2）确定流量—时间函数
对于第1,2未供水时段的流量可直接对水位函数求导，程序如下：
c1=polyfit（t（1:
10）,h（1:
10）,3）;
c2=polyfit（t（11:
21）,h（11:
21）,3）;
a1=polyder（c1）;
a2=polyder（c2）;
tp1=0:
0.01:
8.97;
tp2=10.95:
0.01:
20.84;
x13=-polyval（a1,tp1）;
x113=-polyval（a1,[0:
0.01:
8.97]）;
wgsysl1=100*trapz（tp1,x113）;*/该语句计算第1未供水时段的总用水量*/
x14=-polyval（a1,[7.93,8.97]）;*/该语句仅为下面的程序准备数据*/
x23=-polyval（a2,tp2）;
x114=-polyval（a2,[10.95:
0.01:
20.84]）;
wgsysl2=100*trapz（tp2,x114）;*/该语句计算第2未供水时段的总用水量*/
x24=-polyval（a2,[10.95,12.03]）;*/该语句仅为下面的程序准备数据*/
x25=-polyval（a2,[19.96,20.84]）;*/该语句仅为下面的程序准备数据*/
subplot（1,2,1）
plot（tp1,x13*100）*/与Fig.12.13单位保持一致*/
subplot（1,2,2）
plot（tp2,x23*100）*/与Fig.12.13单位保持一致*/
程序运行得到第1,2未供水时段的-时间流量图如下，可以看到与Fig.12.13中用插值给出的曲线比较吻合。
Fig.12.15
如果用5次多项式拟合则得下面的图形，显然较3次拟合的效果好。
Fig.12.16
而第1供水时段的流量则用前后时期的流量进行拟合得到。
为使流量函数在t=9和t=11连续，我们只取4个点，用3次多项式拟合得到第1供水时段的时间-流量图形如下，可以看到与Fig.12.13中的相应部分比较吻合。
dygsdsj=[7.93,8.97,10.95,12.03];
dygsdls=[x14,x24];
nhjg=polyfit（dygsdsj,dygsdls,3）;
nhsj=7.93:
0.1:
12.03;
nhlsjg=polyval（nhjg,nhsj）;
gssj1=8.97:
0.01:
10.95;
gs1=polyval（nhjg,[8.97:
0.01:
10.95]）;
gsysl1=100*trapz（gssj1,gs1）;*/该语句计算第1供水时段的总用水量*/
plot（nhsj,100*nhlsjg）
Fig.12.17
在第2供水时段之前取t=19.96,20.84两点的流量，用第三未供水时段的3个记录做差分得到两个流量数据21.62,18.48，然后用这4个数据做3次多项式拟合得到第2供水时段与第3为供水时段的时间-流量图如下，可以看到与Fig.12.13中的相应部分也比较吻合。
t3=[19.96,20.84,t（22）,t（23）];
ls3=[x25*100,21.62,18.48];
nhhddxsxs=polyfit（t3,ls3,3）;
tp3=19.96:
0.01:
25.91;
xx3=polyval（nhhddxsxs,tp3）;
gssj2=20.84:
0.01:
24;
gs2=polyval（nhhddxsxs,[20.84:
0.01:
24]）;
gsysl2=trapz（gssj2,gs2）;*/该语句计算第2供水时段的总用水量*/
plot（tp3,xx3）
Fig.12.18
3）一天总用水量的估计
分别对供水的两个时段和不供水的两个时段积分（流量对时间）并求和得到一天的总用水量约为526.8935（此数据是总用水高度cm）。
下表列出各段用水量，与插值法算得的表12-7相比，二者较为吻合。
表12-9
时段
第1未
供水段
第2未
供水段
第1
供水段
第2
供水段
全天用水
用水高度
145.67
260.66
46.60
73.9635
526.8935
4）流量及总用水量的检验
计算出各时刻的流量可用水位记录的数值微分来检验，各时段的用水高度可以用实际记录的水位下降高度来检验。
例如，算得第1未供水段的用水量高度是145.67，而实际记录的水位下降高度为968-822=146，两者是吻合的；同样地，算得第2未供水段的用水量高度是260.66，而实际记录的水位下降高度为1082-822=260，两者也是吻合的。
从算法设计和分析可知,计算结果与各时段所用的拟合多项式的次数有关。
下表给出的是对第1,2未供水时段分别用5、6次多项式拟合后得到的用水量结果。
表12-10
时段
第1未
供水段
第2未
供水段
第1
供水段
第2
供水段
全天用水
用水高度
146.5150
257.7605
46.1317
76.3076
526.7148
（3）结果分析
对方法
（1），由表12-7可以看出,使用三次样条插值法得到的结果（145.6870，258.6547）与表12-5中记录的下降高度146,260相差不大，说明插值结果与原始数据比较吻合。
按三次样条插值法估计出全天的用水量约为539.4450×237.8×10=1282800.21升。
对方法
（2），由表12-10可以得全天的用水量约为526.7148×237.8×10=1252527.7944升。
内容小结
本实验主要进行水塔水流量的估算。
第一种估算方法为插值方法,我们用了三种不同的插值法进行估计。
在求解的过程中,可以熟悉数据插值的理论和方法；第二种估算方法为数据拟合法，用多项式进行拟合，得到水塔水流量的估计。
上机作业
水道测量数据
表1给出了在以码为单位的直角坐标为X，Y的水面一点处以英尺计的水深Z。
水深数据是在低潮时测得的。
船的吃水深度为5英尺。
在矩形区域（75，200）×（-50，150）内的哪些地方船要避免进入。
表1
X129.0140.5103.588.0185.5195.0105.5
157.5107.577.081.0162.0162.0117.5
Y7.5141.523.0147.022.5137.585.5
-6.5-81.03.056.5-66.584.0-33.5
Z4868688
9988949
升级会员 