水塔水流量估計(jì)教學(xué)文案

精品文檔實(shí)驗(yàn) 水塔水流量的估計(jì)實(shí)驗(yàn)?zāi)康谋敬螌?shí)驗(yàn)的主要目的是讓學(xué)生會(huì)用數(shù)學(xué)軟件進(jìn)行插值計(jì)算并解決一些具體的實(shí)際問題。介紹一些經(jīng)典的插值方法， 包括拉格朗日插值法、 埃爾米特插值法、 分段插值法、三次樣條插值法等等。實(shí)驗(yàn)內(nèi)容實(shí)驗(yàn)問題美國(guó)某州的各用水管理機(jī)構(gòu)要求各社區(qū)提供以每小時(shí)多少加侖計(jì)的用水率以及每天所用的總水量。許多社區(qū)沒有測(cè)量流入或流出水塔的水量裝置,他們只能代之以每小時(shí)測(cè)量水塔中的水位,其誤差不超過5%。更重要的是,當(dāng)水塔中的水位下降到最低水位L時(shí)水泵就啟動(dòng)向水塔輸水直到最高水位H,期間不能測(cè)量水泵的供水量。因此,當(dāng)水泵正在輸水時(shí)不容易建立水塔中水位和用水量之間的關(guān)系。水泵每天輸水一次或兩次,每次約二小時(shí)。試估計(jì)任何時(shí)刻(包括水泵正在輸水時(shí)間)從水塔流出的水流量f(t),并估計(jì)一天的總用水量。已知該水塔是一個(gè)高為40英尺（ft），直徑為57英尺(ft)的正圓柱，表12.1給出了某個(gè)小鎮(zhèn)一天水塔水位的真實(shí)數(shù)據(jù)，水位降至約27.00ft水泵開始工作,水位升到35.50ft停止工作。(注：1英尺（ft）=0.3024米（m）)表12-1某小鎮(zhèn)某天水塔水位時(shí)間/s水位/0.01ft時(shí)間/s水位/0.01ft03175466363350331631104995332606635305453936316710619299457254308713937294760574301217921289264554292721240285068535284225223279571854276728543275275021269732284269779254水泵開動(dòng)35932水泵開動(dòng)82649水泵開動(dòng)39332水泵開動(dòng)859683475394353550899533397433183445932703340問題分析流量是單位時(shí)間內(nèi)流出水的體積,由于水塔是正圓柱形,橫截面積是常數(shù),所以在水泵不工作時(shí)段,流量很容易根據(jù)水位相對(duì)時(shí)間的變化率算出。問題的難點(diǎn)在于如何估計(jì)水泵供水時(shí)段的流量。水泵供水時(shí)段的流量只能靠供水時(shí)段前后的流量經(jīng)插值或擬合得到。 作為用于插值或擬合的原始數(shù)據(jù) ,我們希望水泵不工作時(shí)段的流量越準(zhǔn)確越好。這些流量大體上可由兩種方法計(jì)算,一是直接對(duì)表 12.1 中的水量用數(shù)值微分算出各時(shí)段的流量 ,用它們擬合其它時(shí)刻或連續(xù)時(shí)間的流量 ;二是先用表中數(shù)據(jù)擬合水位—時(shí)間函數(shù) ,求導(dǎo)數(shù)即可得到連續(xù)時(shí)間的流量。有了任何時(shí)刻的流量 ,就不難計(jì)算一天的總用水量。 其實(shí),水泵不工作時(shí)段的用水量可以由測(cè)量記錄直接得到 ,由表12.1中下降水位乘以水塔的截面積就是這一時(shí)段的量這精品文檔精品文檔個(gè)數(shù)值可以用來檢驗(yàn)數(shù)據(jù)插值或擬合的結(jié)果。模型建立為了表示方便,我們將問題中所給表12-1中的數(shù)據(jù)全部化為國(guó)際標(biāo)準(zhǔn)單位（表12-5）,時(shí)間用小時(shí)(h),高度用米(m):表12-5一天內(nèi)水塔水位記錄時(shí)間(h)水位(m)時(shí)間(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水泵開動(dòng)9.98水泵開動(dòng)22.96水泵開動(dòng)10.93水泵開動(dòng)23.8810.5910.9510.8224.9910.3512.0310.5025.9110.18模型假設(shè)流量只取決于水位差,與水位本身無關(guān)，故由物理學(xué)中Torricelli定律：從小孔流出的液體的流速正比于水面高度的平方根。題目給出水塔的最低和最高水位分別是8.1648m（27×0.3024）和10.7352m（35.50×0.3024）(設(shè)出口的水位為零)，因?yàn)閟qrt(10.7352/8.1648)≈1.4671，約為1,所以可忽略水位對(duì)流速的影響。2）將流量看作是時(shí)間的連續(xù)函數(shù)，為計(jì)算簡(jiǎn)單,不妨將流量定義成單位時(shí)間流出水的高度,即水位對(duì)時(shí)間變化率的絕對(duì)值(水位是下降的),水塔截面積為S=(57*0.3048)^2* π/4=237.8(m^2). 得到結(jié)果后乘以 S即可。流量估計(jì)方法首先依照表 12-5所給數(shù)據(jù),用MATLAB作出時(shí)間—水位散點(diǎn)圖 12.7。1110.5109.598.58510152025300Fig.12.7下面我們來計(jì)算水箱流量與時(shí)間的關(guān)系。根據(jù)表 12-5中數(shù)據(jù)散點(diǎn)圖 12.7，一種簡(jiǎn)單的處理方法為先將表 12-5中的數(shù)據(jù)分為三段， 然后對(duì)每一段的數(shù)據(jù)做如下處理： 設(shè)某段數(shù)據(jù)為{(x0,y0),(x1,y1),...,(xn,yn)}，相鄰數(shù)據(jù)中點(diǎn)的平均流速用下面的公式 (流速=（左端點(diǎn)的水位 -右端點(diǎn)的水位） /區(qū)間長(zhǎng)度)：精品文檔精品文檔v(xi1xi)(yiyi1)/(xi1xi)，2每段數(shù)據(jù)首尾點(diǎn)的流速用下面的公式計(jì)算：v(x0)(3y04y1y2)/(x2x0)，v(xn)(3yn4yn1yn2)/(xnxn2)，用以上公式求得時(shí)間與流速之間的數(shù)據(jù)表12-6如下:表12-6時(shí)間(h)流速(cm/h)時(shí)間(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.2922.02水泵開動(dòng)9.98水泵開動(dòng)22.96水泵開動(dòng)10.93水泵開動(dòng)23.8827.0910.9533.5024.4321.6211.4929.6325.4518.4825.9113.30由表12-6作出時(shí)間-流速散點(diǎn)圖如下:353025201510510152025300Fig.12.8插值法由表12-6,對(duì)水泵不工作時(shí)段 1,2采取插值方法 ,可以得到任意時(shí)刻的流速 ,從而可以精品文檔精品文檔知道任意時(shí)刻的流量 .我們分別采取拉格朗日插值法 ,分段線性插值法及三次樣條插值法；對(duì)于水泵工作時(shí)段 1應(yīng)用前后時(shí)期的流速進(jìn)行插值，由于最后一段水泵不工作時(shí)段數(shù)據(jù)太少，我們將它與水泵工作時(shí)段 2合并一同進(jìn)行插值處理 (該段簡(jiǎn)稱混合時(shí)段 )。我們總共需要對(duì)四段數(shù)據(jù)（第 1，2未供水時(shí)段，第 1供水時(shí)段，混合時(shí)段）進(jìn)行插值處理,下面以第 1未供水時(shí)段數(shù)據(jù)為例分別用三種方法算出流量函數(shù)和用水量 （用水高度）。下面是實(shí)現(xiàn)該過程的 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為一函數(shù),程序同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')運(yùn)行結(jié)果為 lglrjf=145.6231fdxxjf=147.1430sancytjf=145.6870圖12.9是對(duì)第1未供水段數(shù)據(jù)用三種不同方法得到的插值函數(shù)圖，圖中曲線lglr、fdxx和scyt分別表示用拉格朗日插值法,分段線性插值法及三次樣條插值法得到的曲線。3028262422201816fdxx14scytlglr121234567890Fig.12.9由表12-5知，第1未供水時(shí)段的總用水高度為146（=968-822），可見上述三種插值方法計(jì)算的結(jié)果與實(shí)際值（146）相比都比較接近。考慮到三次樣條插值方法具有更加良好的性質(zhì)，建議采取該方法。其它三段的處理方法與第 1未供水時(shí)段的處理方法類似， 這里不再詳細(xì)敘述，只給出數(shù)值結(jié)果和函數(shù)圖象，圖中曲線標(biāo)記同 Fig.12.9。表12-7 各時(shí)段及一天的總用水量 (用水高度)精品文檔精品文檔第1第2第1混合全天未供水段未供水段供水段時(shí)段拉格朗日145.6231258.866454.068992.1337550.6921插值法分段線性147.1430258.969749.605176.4688532.1866插值法三次樣條145.6870258.654753.333481.7699539.4450插值法35lglr30scytfdxx2520151089101112137Fig.12.10 (第一供水段時(shí)間 -流速示意圖)3432302826 lglrscytfdxx24222010 12 14 16 18 20 22Fig.12.11 (第2未供水段時(shí)間 -流速示意圖)35lglr30scyt25fdxx2015102021222324252619Fig.12.12(混合時(shí)段段時(shí)間-流速示意圖)下圖12.13是用分段線性及三次樣條插值方法得到的整個(gè)過程的時(shí)間-流速函數(shù)示意精品文檔精品文檔圖。3530scyt25201510

fdxx0 5 10 15 20 25 30Fig.12.13下表是對(duì)一天中任取的 4個(gè)時(shí)刻分別用三種方法得到的水塔水流量近似值。表12-86.8810.8815.8822.88拉15.9826671234851433.7426009085346325.5662241818047734.70996210551694格朗日分14.8272413793103432.9976262626262725.4465591397849525.4715789473684段插值三15.0527858196582033.7089553595325925.5490892055742329.41733175863551次樣條(2)擬合法1)擬合水位--時(shí)間函數(shù)從表12-5中測(cè)量記錄看，一天有兩次供水時(shí)段和三次未供水時(shí)段，分別對(duì)第1,2未供水時(shí)段的測(cè)量數(shù)據(jù)直接作多項(xiàng)式擬合，可得到水位函數(shù)（注意，根據(jù)多項(xiàng)式擬合的特點(diǎn)，此處擬合多項(xiàng)式的次數(shù)不宜過高，一般以3-6次為宜）。對(duì)第3未供水時(shí)段來說,數(shù)據(jù)過少不能得到很好的擬和。設(shè)t,h分別為已輸入的時(shí)刻和水位測(cè)量記錄(由表12-5提供，水泵啟動(dòng)的4個(gè)時(shí)刻不輸入)，這樣第1未供水時(shí)段各時(shí)刻的水位可由如下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未供水時(shí)段的時(shí)間 -水位擬合函數(shù)圖形。109.598.580246810Fig.12.14變量x1中存放了以0.1為步長(zhǎng)算出的各個(gè)時(shí)刻的水位高度。同樣地，第2未供水時(shí)段時(shí)間-水位圖可由如下MATLAB程序完成，讀者可自己上機(jī)運(yùn)行查看。c2=polyfit(t(11:21),h(11:21),3);tp2=10.9:0.1:20.9;x2=-polyval(c2,tp2);plot(tp2,x2)確定流量—時(shí)間函數(shù)對(duì)于第1,2未供水時(shí)段的流量可直接對(duì)水位函數(shù)求導(dǎo)，程序如下：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); */該語句計(jì)算第 1未供水時(shí)段的總用水量 */x14=-polyval(a1,[7.93,8.97]); */該語句僅為下面的程序準(zhǔn)備數(shù)據(jù) */x23=-polyval(a2,tp2);x114=-polyval(a2,[10.95:0.01:20.84]);wgsysl2=100*trapz(tp2,x114); */ 該語句計(jì)算第 2未供水時(shí)段的總用水量 */x24=-polyval(a2,[10.95,12.03]); */該語句僅為下面的程序準(zhǔn)備數(shù)據(jù) */x25=-polyval(a2,[19.96,20.84]); */該語句僅為下面的程序準(zhǔn)備數(shù)據(jù) */subplot(1,2,1)plot(tp1,x13*100) */與Fig.12.13 單位保持一致*/subplot(1,2,2)plot(tp2,x23*100) */與Fig.12.13 單位保持一致 */程序運(yùn)行得到第 1,2未供水時(shí)段的-時(shí)間流量圖如下，可以看到與 Fig.12.13中用插值給出的曲線比較吻合。精品文檔精品文檔2334222132203019182817261624151422051010152025Fig.12.15如果用5次多項(xiàng)式擬合則得下面的圖形，顯然較3次擬合的效果好。2431233022292120281927182617162515241423051010152025Fig.12.16而第1供水時(shí)段的流量則用前后時(shí)期的流量進(jìn)行擬合得到。 為使流量函數(shù)在 t=9和t=11連續(xù)，我們只取4個(gè)點(diǎn)，用3次多項(xiàng)式擬合得到第1供水時(shí)段的時(shí)間-流量圖形如下，可以看到與Fig.12.13中的相應(yīng)部分比較吻合。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); */該語句計(jì)算第 1供水時(shí)段的總用水量 */plot(nhsj,100*nhlsjg)精品文檔精品文檔3230282624222018161478910111213Fig.12.17在第2供水時(shí)段之前取 t=19.96,20.84兩點(diǎn)的流量，用第三未供水時(shí)段的 3個(gè)記錄做差分得到兩個(gè)流量數(shù)據(jù) 21.62,18.48 ，然后用這 4個(gè)數(shù)據(jù)做 3次多項(xiàng)式擬合得到第 2供水時(shí)段與第 3為供水時(shí)段的時(shí)間 -流量圖如下，可以看到與 Fig.12.13中的相應(yīng)部分也比較吻合。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); */該語句計(jì)算第 2供水時(shí)段的總用水量 */plot(tp3,xx3)262422201816142021222324252619Fig.12.18一天總用水量的估計(jì)分別對(duì)供水的兩個(gè)時(shí)段和不供水的兩個(gè)時(shí)段積分總用水量約為526.8935(此數(shù)據(jù)是總用水高度算得的表12-7相比，二者較為吻合。

cm)

(流量對(duì)時(shí)間)并求和得到一天的。下表列出各段用水量，與插值法表

12-9時(shí)段

第1未供水段

第2未供水段

第1供水段

第2供水段

全天用水用水高度

145.67

260.66

46.60 73.9635

526.89354）流量及總用水量的檢驗(yàn)精品文檔精品文檔計(jì)算出各時(shí)刻的流量可用水位記錄的數(shù)值微分來檢驗(yàn)， 各時(shí)段的用水高度可以用實(shí)際記錄的水位下降高度來檢驗(yàn)。例如，算得第 1未供水段的用水量高度是 145.67 ，而實(shí)際記錄的水位下降高度為 968-822=146，兩者是吻合的；同樣地，算得第