2010年數(shù)學建模B題(儲油罐問題)

1、摘 要于地基變形而引起的油液面傾斜等因素對罐容表的影響。將理論推導和數(shù)據(jù)擬合情況綜合分析，在理論推導方面，創(chuàng)新性的運用祖規(guī)則變換，最終求得不規(guī)則立體的體積。探討了使用 SURFER 軟件對體積網格化求不規(guī)則立體體積的方法。對兩端平頭的橢圓柱體形小橢圓型儲油罐無變位和傾斜（傾斜角 4.1=）可求出理論推導式，再用 Matlab對實際所測數(shù)據(jù)進行擬合得出近似方程。對近際情況的方程。對實際的儲油罐變位情況（縱向傾斜角度 橫向傾斜角度）建立罐容，一端使其轉化成規(guī)則體元，計算出總的體積，減去虛擬體積。采用 Matlab符號 與體積 ， ，Vh際測量數(shù)據(jù)擬合公式做比較，求出體積微小差異量，進行誤差分析。結

2、果表明，此模型與實際測量數(shù)據(jù)吻合程度較好。關鍵詞：祖暅原理截面轉化;等效變換;虛擬體積;體積網格化1一 問題重述通常加油站都有若干個儲存燃油的地下儲油罐，并且一般都有與之配套的“油位計量管理系統(tǒng)”，采用流量計和油位計來測量進/出油量與罐內油位高度時計算，以得到罐內油位高度和儲油量的變化情況。許多儲油罐在使用一段時間后，由于地基變形等原因，使罐體的位置會發(fā)生縱向傾斜和橫向偏轉等變化（以下稱為變位），從而導致罐容表發(fā)生改變。按照有關規(guī)定，需要定期對罐容表進行重新標定。問題一平放置而無變位和傾斜（傾斜角 4.1 時的情況進行分析得出小橢圓型儲油=罐的罐容與標油計的計量 h 同罐容 的函數(shù)關系，并作出

3、它們的罐容表進行比V較，同時求出間隔為 1cm 罐容表標定體積值。問題二由于地基的變化從而引起油罐傾斜而使原來的“油位計量管理系統(tǒng)”對傾標油計所測的實際數(shù)據(jù)測量儲油罐變位即縱向傾斜角度 和橫向傾斜角度 同 時變化情況同油標記的計量 h 與油罐體積v隔為 10cm 罐容表標定體積值。,h 的函數(shù)關系，并求出求出間二 問題分析問題一（1）由題意可知對于小橢圓型儲油罐無變位情況我們只需要找出油標記的計量 h公式即可求出罐容與小橢圓型儲油罐無變位的油標記的計量 的關系。h（2）對于小橢圓型儲油罐變位情況即傾斜角 縱向變位（如圖3不規(guī)則，這就需要我們對其進行分割，利用積分求解，得出理論方程。同時利用Ma

4、tlab對實測數(shù)據(jù)擬合得出擬合方程，再進行比較和修正。 n n1（3）利用以上所得方程，帶入hnN* n120即可求出罐體變且2位后油位高度間隔為 1cm的罐容表標定值。問題二2割，切割一小塊利用微元，利用 Matlab符號運算工具箱，推導出變位油罐液面高 與體積 hV進行誤差分析，通過合理猜想判斷誤差來源，進而找到最優(yōu)結果。三 模型假設1 由于溫度的變化影響油的體積變化較小，我們將其忽略不計。2 由于油對罐面具有一定的粘度，但是在實際情況中罐壁上的粘度與底座是不同的，并且它還受溫度的影響，在這里我們將其忽略，進行簡化。3 油罐傾斜時，在傾斜腳出有一定的體積用油標記是無法測到的，我們也將其簡化

5、，求出最大限度的體積看作該處的體積。4 因為油浮子，進出口管都是有一定體積的，我們利用積分法求油的體積的時候是沒有考慮的，如果考慮比較復雜，我們將其忽略不計，最后進行修正。5 由于油罐可看作是一剛體，所以其形狀不發(fā)生改變。四 符號說明在沒有標明情況下，長度單位默認為（m m3角度單位默認為弧度（）油面高度測量值hS 油面與圓形截得的面積（圖 2）1油面與橢圓里截得的面積（圖 2）S2V ,a,l, 體積油高函數(shù)H 縱變位后油液面的高端液高gH 縱變位后的低端液高dH 縱變位后油面等效于無變位液高H 縱變位后的低端液高k3五 模型的建立與求解lh面橢圓的長半軸為a米，短半軸為 米。以橢圓的中心為

6、坐標原點，長、短半軸b所在的直線為x 軸、y軸，建立如圖 1 所示的空間直角坐標系zybhlOaxb4ybDAOBxbhb得的面積 S 與在橢圓里截得的面積圖 2 S 之比例也是b2計算油面與圓形截得的面積S1bbb；211 b b 222222222利用比例關系，計算油面與橢圓面截得的面積S2a SR22Sa22綜上，體積便很容易得到： 21.78l m,a 22式，并化簡可得54361 4361 V h 32000 1200 此式即為理論上最簡單，最理想化的無位變臥式橢圓罐部分體積函數(shù)V h從這一函數(shù)表達式出發(fā)，可以進一步討論油罐總體積，即最大儲油量V Vmax對該函數(shù)求一階導數(shù)并描繪出體

7、積函數(shù) V h首先曲線滿足實際的增函數(shù)要求，其次從一階導數(shù)看，通過它反映出來的體積函數(shù)的增長情況正好也符合實際的：先曾得快，到半短軸（即 0.6m）時， 進油結束時的狀態(tài)。但理論公式還是上面推導出來的那個V h 函數(shù)。632V3232b) 近似計算 V h 函示意圖如右邊所示（已在圖上標注出關鍵尺寸和相關假設的長度）由圖 3 可知，OABC 即為傾圖3 7hHg區(qū)域：0 H Ltan；g區(qū)域：Ltan H b ；g區(qū)域：b H 2b0.4tang當H b0.4tan g討論的意義了。由于H 并非直接測量，為了更加方便，運用簡單的直角三角形g的角、邊關系，可以導出不同區(qū)域下高端液高H 與顯示液高

8、h的關系為g區(qū)域：H h htan 1g區(qū)域：H h htan 1g區(qū)域：H h htan 1g容易想到，罐內油的體積不會因為罐體的位變而發(fā)生變化，所以只要有一個高端液高H 值就一定存在一個與之對應gg標，以H為橫坐標描點2，發(fā)現(xiàn)：gg當Ltan H b 時H 與H的關系也是近似一條直線，如上圖中的AB 線段；gggg知這條直線與OA 段斜率一樣。g為H kH （k Hgkg 再把H帶入到1.1中導出的無位變臥式橢圓罐部分體積函數(shù)V h 中。從而即可8得到現(xiàn)在我們需要的縱向位變臥式橢圓罐部分體積函數(shù)V ,k,l, ,a,b 如下2111 1V b1bh hhh kk 0hb 1.78帶縱變位因

9、子 2.45 ,a 0.89m,lm21.22b0.6m帶入上式可得5h tan 1 43614361 tan 1h 13083100001308310000arccos1V k20001200k 6h tan15k20h1.20.4tanhk41在本實例中已知了 4.1 ，再代入上式，即得18005 1 1 hhVk 1200k 6h 15k20h0.4tan 1 hk從上式看出還有一個待定系數(shù)k ，我們擬通過實驗實測的V h數(shù)據(jù)，反解出一系列的k值，再利用統(tǒng)計學的方法，從而能很輕松地確定系數(shù)k。的分析，能求出k。再回帶到這里的表達式，即可表出這里的函數(shù)表達式。ii. 0.14694h1.2

10、段與前一分段OA在拐點A處不相交，只與其延長線相交，交點為A，且不經過圓點。還應強調的是，點A的縱坐標不等于Ltan ，通過驗證，點A的縱坐標H 1.0765Ltan 。通過計算gA直線的方程的可以把H表出為： bkltan H ltanltanb 1k H gk bltan（這里H 1 h）g 同上，再把H帶入到1中導出的無位變臥式橢圓罐部分體積函數(shù)V h 中。9從而即可得到現(xiàn)在我們需要的縱向位變臥式橢圓罐部分體積函數(shù) ,V h kV ,k,l,a,b 。然后把本實例中已知的參數(shù)l。,a,b, 帶入上式可得出我們擬通過實驗實測的 數(shù)據(jù)，來確定系數(shù) 。V hk i ，然后求得 S 對應nn找

11、的算法是先求出 V h k,1,2kSiijjij1的i是多少，這個i對應的k 具體步驟為：i Step1把實測的 數(shù)據(jù)帶入V ,k 函數(shù)內，解得若干個 值；V hk Step2 把上面求得的第i個 值與第 個實測的液面高度h 帶入V ,k 的函kjj數(shù)內，求得第j個測量的油的體積V ；jStep3先對V 取絕對值，再對j求和，得到S ；jiStep4換一個i值，重復上面步驟，直到隸遍及所有的 ；kStep4找求和后的最小值所對應的i值；Step5此時的k 即為我們所求。i這樣最終得到 1.934575156142711.934575kiii. 1.2 1.22867h對稱性此時仍采用關系式H

12、 kH g式變換：低端液高H H ltan ；低端液面空高H bH ，再用H 代替H ，dgkdkgbH l則H 1中導出的無位變臥式橢圓罐部分體積gk 函數(shù)V h 中。從而即可得到現(xiàn)在我們需要的縱向位變臥式橢圓罐部分體積函數(shù), ,V ,k,l, a b l,a,b, ,以及ii中求得的參數(shù) k 即能得到此時的縱位變臥式橢圓罐部分體積函數(shù)V h 。這里需要說明一點：用Matlab 的符號運算器進行求解。詳情，請參考附錄。至此，即可由兩種方法 來表出V h 函數(shù)：其一為多項式曲線擬合出來的；其二為近似計算得到的（分10標高標高31 8113191128252323 233323394494947

13、43259595558 96767668 78767678 986818581 995299 9161619 7292125 3383135 9474142 155352 261626767472737828781891998問題二：油灌變位問題可以簡單的看做為油管先沿水平方向傾斜了一個角度 ，然后沿xoy軸旋轉了一個角度 。yxx y 軸， 為z沿傾斜指向上升方向。先計算球缺最大圓弧的半徑R：R 沿xoy面對油面體積進行切割，得到的圓方程為：x y1.5 1.5222 222圓心的為圓心的兩組平行于 yoz圓系方程2 y22282222 沿xoy方向去掉圓缺后，得到截面圓的方程為：x y1.

14、5 1.5222y z h （雖隨著 的變化，可以看出，油罐液面直線方程會隨著高h0h66 h3232 h3首先，沿 yoz面對油面體積進行切割，得到一組近似扇形的油面體積微圓由于油面體積是沿x軸方向對稱的，所以V 2v，只要計算出v就可以間接得到油罐剩余油量的體積V 。 已知dvsdx，而對于切割后的截面面積 s，有ds f y f y dy （設12f z 為沿 yoz 面方向切割得到左邊同心弧長的函數(shù)，則12yzh h y 2 R x y 228計算 的積分區(qū)間：對于任意微小圓 來說，它的積分區(qū)間下限為：得到 離 軸的距離x， y 帶入沿xoyx y1.5 1.5 222解得到 1.5

15、1.5即為積分下限。yx22同樣對于積分上限：積分上限為以球缺最大圓弧的圓心為圓心的平行于yoz1182 圓系方程 1.5 R x 與直線y ztanhcot的交點，于是y2z22hcot y求解直線z 帶入圓系方程得到：tancot y 11h2 ，所以y1.5 R x2求解得到積分上限： y f h22tan81 h yfh s R x y 2228y 2x2計算體積v：已知 dvsdx討論dv的積分上限： 2x y1.5 1.522對于圓系方程因為 y方向的坐標有y h h替 換 y得到得到 f x,h 0求解得到x f h 的方程，即為的積分上限。所以2 h yfhfhf h 2221

16、v R x y 22800y 2x2最后得到 h yhf hfh 2f221V 2v2 2R x y 22800y 2x20h6tan ；2.求解當6tan h32tan油罐的體積，同樣采用補體積的方法，先求解如圖所示的陰影部分的體積V ，然后減去彌補的體積V ，則體積V V V 。1212146tan h2tan這一部分的體積分為三個區(qū)間，為 1.52tan h1.56tan1.56tan h32tan當 6tan 1.52tan 時V 的體積與第一步計算的體積推導公式相h1同： h yhfhf h 2f221V v2 2R x y 228100y 2x2 f y f y dy計算V 的體積

17、與第一步計算 的推導過程相同 dsV212 f y 的方程發(fā)生變化： f y 的方程變?yōu)橐郧蛉弊畲髨A弧的圓心為圓心的平行11于 yoz圓系方程 2y1.5 z5.375 R x222則求解得后兩函數(shù)想減積分得到：h y 2 R x y 22然后計算的積分區(qū)間：對于任意微小圓來說，它的積分區(qū)間下限為：得到 離 y 軸的距離x， 帶入沿xoyx y1.5 1.5 222解得到 y 1.5 x 1.5即為積分下限。22同樣對于積分上限：積分上限為以球缺最大圓弧的圓心為圓心的平行于yoz圓系方程 2y1.5 z5.375 R x 與直線 y zh 的交點，于是求222hcot y解直線z 帶入圓系方程

18、得到：tan hcot y 5.375 R x 求解得到積分上限：y f h2y1.5 222tan1 hcot y1fh 所以s R x y1.5 5.375dy222tany 2x21.5計算體積v：已知dvsdx討論dv 的積分上限：15 2x y1.5 1.522對于圓系方程因為 y方向的坐標有y h h替 的方程：h換 得到得到y(tǒng), 0求解得到 f x hx f2即為 的積分上限。所以dv cot ytan f2hf2h1 h 2hvsdxR x y1.5 5.375dydx2200y 2x21.5最終得到 cot y1 fhfhfh 2h22V 2v2sdx2R x y1.5 5.

19、375dydx22tan200y 2x21.5最終得到體積的計算公式：V V V12 11 hcot yfhf h 21 2 2R x y1.5 dydx2228tan0y 2x2 cot yh h R x y1.5 5.375hf212dydx22tan0y 2x26tanh1.52tan同理計算當1.52tan h 1.5 6tan油罐的體積：V V V12 11 hcot yf h 11 2R x y1.5 dydx2228tan0y 2x2 hcot yfhfh 21 2R x y1.5 5.375dydx222tan0y 1.52x21.52tanh1.56tan當1.56tan h

20、 3 2tan油罐的體積：V V V12 11 hcot y1f h 1 2R x y1.5 dydx2228tan0y 2x2 hcotyfh 11 2R x y1.5 5.375dydx222tan0y 2x21.56tanh32tan163.推導32tan 3時不規(guī)則體的體積：h推導32tan 3時不規(guī)則體的體積，可以通過轉化先求的油罐的整個h體積（正常情況下時可以通過推導圓柱體體積與兩個圓缺的體積推導出來的）H rH2V R L滿233*1.625*1.5 *623得到油罐體積后通過計算上班部分剩余的體積V，就可以得到油罐中的體積公式?，F(xiàn)在計算剩余部分的體積公式：利用第一步得到的公式，

21、進行適當變化的到剩余體積公式： cot ytanfhfh f h 2h221Vv2 sdx2R2x2 y1.5 5.375dydx f h003 f h x y1.51.5 求解得到；2221 6tanf h x hcot 1.5 求解得到；2222cot ytanh2 f h 由5.375 y1.5 R x 求解得到。2223則體積V V V滿3*1.625 fhcot yh f h R x y 5.3751.5*1.5 *6221dydx22223tan f h0332tanh3下面考慮油罐的扭曲角 ，17 cosh實標 通過觀察上圖，可以發(fā)現(xiàn)實標標所以通過置換體積公式中的h 為h 即可得

22、到油罐變體后的體積公式：標實當0 h 6tan 時：實fhfhfh 2221V 2v2 2R x y 22800當6tan h 32tan時：實V V V12 11 hcot y fhf h 21 2R x y1.5 dydx2228tan0y 2x2 hcot yf2h h 2R x y1.5 5.375dydx222tan0y 2x26tanh1.52tan V V V12 11 hcot y f h 1 2R x y1.5 dydx2228tan0y 2x2 cot yfhfh h21 2R x y1.5 5.375dydx222tan0y 2x21.52tanh1.56tanV V V

23、12 cottan11h yf h 1 2R x y1.5 dydx22280y 2x2 fh hcot y1 2R x y1.5 5.375dydx222tan0y 2x21.56tanh32tan當 32tan 3 h實V V V滿 3*1.625 cot ytanh f h hf*1.5 *6221R x y1.55.375dydx22223 f h0332tanh3利用 MATLAB 符號運算工具箱，對建立起來的數(shù)學模型進行化簡求得, ) 0，帶入數(shù)據(jù)計算 。fV,h ,得到的曲線與實際描點符合程度比較好，擬合得到的數(shù)據(jù)曲線方程為：V h 12.297h h3219 h h h h32

24、六 模型優(yōu)化與改進數(shù)學模型更加符合實際。積，在查閱文獻后發(fā)現(xiàn)這種方法是可行的，遺憾的是，因為時間的急迫，無法再這次建模的論文中建立起科操作的模型。七 模型評價本模型的最大特點是盡量避免了復雜的積分的運算，通過近似計算的方法來得出液高與體積的關系。雖然 Matlab也帶有相當完備的 symbol符號計算功能。機返回的符號運算結果帶來了巨大困難，本題模型一油罐平放式通過使用祖縆公式，放棄使用積分公式，使得體積的計算更加簡化，大大提高了計算效率，不失為一種好方法。k,遺憾的是由于所給的數(shù)據(jù)量而無法達到想要的精度。積。21參考文獻1 田鐵軍，傾斜臥式罐橢直圓筒部分容積的近似計算J，現(xiàn)代測量與實驗室管理

25、，2004 年 01 期：26-27，2004.2 羅亦泳，張立亭，陳竹安，楊偉，基于 Surfer 的數(shù)據(jù)網格化與體積計算精度分析J，測繪科學，2009 年 05 期：97-99，2009.3 張超英，李華兵，譚惠麗，劉慕仁，孔令江，橢圓柱體在牛頓流體中運動的格子 Boltzmann 方法模擬J，物理學報，2005 年 05 期：1982-1987，2005.4 年 03期：30-34，2010.5 年 04期：16-17，2001.6 王若鵬, 吳國民 認識，2007 年 16 期：1-6，2007.7 張偉，祖暅原理的由來及證明J，重慶教育學院學報，2010 年 03 期：113-115

26、，2010.8 同濟大學數(shù)學系，高等數(shù)學M，北京： 高等教育出版社，2008.22附錄1Matlab核心計算程序clc;clear;load V_h_in1;Vtheory=ones(length(h),1);syms yV2=13083/10000*pi-13083/10000*acos(-1+5/3*y)-4361/1200*(3/5-y)*(6/5*y-y2)(1/2);for i=1:length(h)Vtheory(i)=subs(V2,y,h(i);enda=polyfit(h,V,3);%下面代碼用于繪制出差值曲線、擬合曲線與實測數(shù)據(jù)點hi=h(1):0.0001:h(length(h); %設置插值點Vi2=a(1)*hi.3+a(2)*hi.2+a(