水塔水流量估计

《水塔水流量估计》由会员分享，可在线阅读，更多相关《水塔水流量估计（12页珍藏版）》请在装配图网上搜索。

1、实验 水塔水流量的估计实验目的本次实验的主要目的是让学生会用数学软件进行插值计算并解决一些具体的实际 问题。介绍一些经典的插值方法，包括拉格朗日插值法、埃尔米特插值法、分段插值法 三次样条插值法等等。实验内容1 实验问题 美国某州的各用水管理机构要求各社区提供以每小时多少加仑计的用水率以及每天所用的总水量。许多社区没有测量流入或流出水塔的水量装置,他们只能代之以每小 时测量水塔中的水位，其误差不超过5%。更重要的是，当水塔中的水位下降到最低水位L 时水泵就启动向水塔输水直到最高水位H,期间不能测量水泵的供水量。因此，当水泵正 在输水时不容易建立水塔中水位和用水量之间的关系。水泵每天输水一次或两

2、次，每次 约二小时。试估计任何时刻（包括水泵正在输水时间）从水塔流出的水流量f（t）,并估计一天 的总用水量。已知该水塔是一个高为40英尺（ft）,直径为57英尺（ft）的正圆柱，表 12.1给出了某个小镇一天水塔水位的真实数据，水位降至约27.00ft水泵开始工作，水 位升到35.50ft停止工作。（注：1英尺（ft） =0.3024米（m）表 12-1 某小镇某天水塔水位时间/s水位/0.01f t时间/s水位/0.01f t03175466363350331631104995332606635305453936316710619299457254308713937294760574301

3、217921289264554292721240285068535284225223279571854276728543275275021269732284269779254水泵开动35932水泵开动82649水泵开动39332水泵开动8596834753943535508995333974331834459327033402 问题分析流量是单位时间内流出水的体积，由于水塔是正圆柱形，横截面积是常数，所以在水 泵不工作时段，流量很容易根据水位相对时间的变化率算出。问题的难点在于如何估计 水泵供水时段的流量。水泵供水时段的流量只能靠供水时段前后的流量经插值或拟合得到。作为用于插值 或拟合的原始数

4、据,我们希望水泵不工作时段的流量越准确越好。这些流量大体上可由 两种方法计算,一是直接对表 12.1 中的水量用数值微分算出各时段的流量,用它们拟合 其它时刻或连续时间的流量;二是先用表中数据拟合水位时间函数 ,求导数即可得到 连续时间的流量。有了任何时刻的流量,就不难计算一天的总用水量。其实,水泵不工作时段的用水量 可以由测量记录直接得到,由表12.1中下降水位乘以水塔的截面积就是这一时段的量这 个数值可以用来检验数据插值或拟合的结果。模型建立为了表示方便,我们将问题中所给表12-1中的数据全部化为国际标准单位(表 12-5), 时间用小时(h)，高度用米(m):表12-5 天内水塔水位记录

5、时间(h)水位(m)时间(h)水位(m)09.6812.9510.210.929.4513.889.941.849.3114.989.652.959.1315.909.413.878.9816.839.184.988.8117.938.925.908.6919.048.667.008.5219.968.437.938.3920.848.228.978.2222.02水泵开动9.98水泵开动22.96水泵开动10.93水泵开动23.8810.5910.9510.8224.9910.3512.0310.5025.9110.18模型假设(1) 流量只取决于水位差，与水位本身无关，故由物理学中 Tor

6、ricelli 定律：从小孔 流出的液体的流速正比于水面高度的平方根。题目给出水塔的最低和最高水位分别是 8.1648m (27 X 0.3024 )和 10.7352m ( 35.50 X 0.3024 )(设出 口的水位为零)，因为 sqr t(10.7352/8.1648p1.1467，约为1,所以可忽略水位对流速的影响。(2) 将流量看作是时间的连续函数，为计算简单，不妨将流量定义成单位时间流出水 的高度,即水位对时间变化率的绝对值(水位是下降的), 水塔截面积为 S=(57*0.3048厂2*n/4=237.8(m2).得到结果后乘以 S 即可。流量估计方法首先依照表12-5所给数据

7、，用MATLAB作出时间一水位散点图12.7。051015202530Fig.12.7 下面我们来计算水箱流量与时间的关系。根据表12-5中数据散点图 12.7，一种简 单的处理方法为先将表12-5 中的数据分为三段，然后对每一段的数据做如下处理：设某段数据为(X y ),(x y ) ,.,(x y )，相邻数据中点的平均流速用下面的公式(流速=0 0 1 1 n n(左端点的水位-右端点的水位)/区间长度)：v(x + x、i+4i-)=2(y厂 +i)/( x+i每段数据首尾点的流速用下面的公式计算：v(x )二(3y -4y + y )/(x -x ),0 0 1 2 2 0v(x )

8、二(-3y + 4y - y)/(x - x ),n nn -1n - 2nn - 2用以上公式求得时间与流速之间的数据表12-6 如下:表 12-6时间(h)流速(cm/h)时间(h)流速(cm/h)029.8912.4931.520.4621.7413.4229.031.3818.4814.4326.362.39516.2215.4426.093.4116.3016.3724.734.42515.3217.3823.645.4413.0418.4923.426.4515.4519.5025.007.46513.9820.4023.868.4516.3520.8422.178.9719.29

9、22.02水泵开动9.98水泵开动22.96水泵开动10.93水泵开动23.8827.0910.9533.5024.4321.6211.4929.6325.4518.4825.9113.30由表 12-6 作出时间-流速散点图如下:10 051015202530Fig.12.8(1) 插值法由表 12-6,对水泵不工作时段1,2 采取插值方法,可以得到任意时刻的流速,从而可以 知道任意时刻的流量.我们分别采取拉格朗日插值法 ,分段线性插值法及三次样条插值 法；对于水泵工作时段1 应用前后时期的流速进行插值，由于最后一段水泵不工作时段 数据太少，我们将它与水泵工作时段2 合并一同进行插值处理(该

10、段简称混合时段)。我们总共需要对四段数据(第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

11、*/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和scy

12、t分别表示用拉格朗日插值法，分段线性插值法及三次样条插值法得到的曲线。由表 12-5 知，第 1 未供水时段的总用水高度为 146（=968-822），可见上述三种插 值方法计算的结果与实际值（146）相比都比较接近。考虑到三次样条插值方法具有更 加良好的性质，建议采取该方法。其它三段的处理方法与第1 未供水时段的处理方法类似，这里不再详细叙述，只给 出数值结果和函数图象，图中曲线标记同Fig.12.9o表 12-7 各时段及一天的总用水量（ 用水高度）第1未供水段第2未供水段第1供水段混合 时段全天拉格朗日 插值法145.6231258.866454.068992.1337550.6921分

13、段线性 插值法147.1430258.969749.605176.4688532.1866三次样条 插值法145.6870258.654753.333481.7699539.4450Fig.12.10 （第一供水段时间-流速示意图）101920212223242526Fig.l2.12（混合时段段时间-流速示意图）下图12.13是用分段线性及三次样条插值方法得到的整个过程的时间-流速函数示意 图。10051015202530Fig.12.13下表是对一天中任取的4 个时刻分别用三种方法得到的水塔水流量近似值。表 1286.8810.8815.8822.88拉 格 朗 日15.982667123

14、4851433.7426009085346325.5662241818047734.70996210551694分 段 插 值14.8272413793103432.9976262626262725.4465591397849525.4715789473684三次样条15.0527858196582033.7089553595325925.5490892055742329.41733175863551(2) 拟合法1) 拟合水位-时间函数从表 12-5 中测量记录看，一天有两次供水时段和三次未供水时段，分别对第 1,2 未供水时段的测量数据直接作多项式拟合，可得到水位函数(注意，根据多项式拟合

15、的 特点，此处拟合多项式的次数不宜过高，一般以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.66h=9.68,9.48,9.31,9.13,8.98,

16、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未供水时段的时间-水位拟合函数图形。变量x1中存放了以0.1为步长算出的各个时刻的水位高度。同样地，第2未供 水时段时间冰位图可由如下MATLAB程序完成，读者可自己上机运行查看。c2=polyfit(t(11:21),h(11:

17、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*tra

18、pz(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

19、,1)plo t(t p1,x13*100)*/与 Fig.12.13 单位保持一致*/subplot(1,2,2)plo t(t p2,x23*100)*/与 Fig.12.13 单位保持一致*/程序运行得到第1,2未供水时段的-时间流量图如下，可以看到与Fig.12.13中用插 值给出的曲线比较吻合。Fig.12.15如果用 5次多项式拟合则得下面的图形，显然较3次拟合的效果好。Fig.12.16而第1供水时段的流量则用前后时期的流量进行拟合得到。为使流量函数在t=9和 t=ll连续，我们只取4个点，用3次多项式拟合得到第1供水时段的时间-流量图形如 下，可以看到与Fig.12.13中的相

20、应部分比较吻合。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);gssj18=.97:0.01:10.95;gs1=polyval(nhjg,8.97:0.01:10.95);gsysl1=100*trapz(gssj1,gs1); */该语句计算第1供水时段的总用水量*/ plot(nhsj, 100*nhlsjg)3230282624222018168910111213147Fi

21、g.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=polyv

22、al(nhhddxsxs, 20.84:0.01:24);gsysl2=trapz(gssj2,gs2); */该语句计算第2供水时段的总用水量*/plot(tp3,xx3)3) 一天总用水量的估计 分别对供水的两个时段和不供水的两个时段积分 (流量对时间 )并求和得到一天的总用水量约为526.8935(此数据是总用水高度cm)。下表列出各段用水量，与插值法 算得的表12-7相比，二者较为吻合。表12-9时段第1未 供水段第2未 供水段第1 供水段第2 供水段全天用水用水咼度145.67260.664 6.6073.9635526.89354) 流量及总用水量的检验计算出各时刻的流量可用水位记

23、录的数值微分来检验，各时段的用水高度可以用实 际记录的水位下降高度来检验。例如，算得第1未供水段的用水量高度是145.67，而 实际记录的水位下降高度为 968-822=146，两者是吻合的；同样地，算得第2 未供水段 的用水量高度是260.66，而实际记录的水位下降高度为1082-822=260，两者也是吻合 的。从算法设计和分析可知,计算结果与各时段所用的拟合多项式的次数有关。下表给 出的是对第1,2未供水时段分别用5、6次多项式拟合后得到的用水量结果。表 12-10时段第1未 供水段第2未 供水段第1 供水段第2 供水段全天用水用水咼度146.5150257.760546.131776.

24、3076526.71483)结果分析对方法 (1)，由表 12-7 可以看出 ,使用三次样条插值法得到的结果( 145.6870 ， 258.6547)与表12-5中记录的下降高度146,260相差不大，说明插值结果与原始数据比 较吻合。按三次样条插值法估计出全天的用水量约为539.4450 X 237.8 X 10= 1282800.21升。对方法,由表12-10可以得全天的用水量约为526.7148 X 237.8 X 10=1252527.7944 升。内容小结 本实验主要进行水塔水流量的估算。第一种估算方法为插值方法,我们用了三种不 同的插值法进行估计。在求解的过程中,可以熟悉数据插值的理论和方法；第二种估算 方法为数据拟合法,用多项式进行拟合,得到水塔水流量的估计。上机作业水道测量数据表1给出了在以码为单位的直角坐标为X, Y的水面一点处以英尺计的水深Z。水深 数据是在低潮时测得的。船的吃水深度为5英尺。在矩形区域（75, 200）X（-50, 150）内的哪些地方船要避 免进入。表1X129.0140.5103.588.0185.5195.0105.5157.5107.577.081.0162.0162.0117.5Y7.5141.523.0147.022.5137.585.5-6.5-81.03.056.5-66.584.0-33.5Z48686889988949