1、儲油罐的變位識別與罐容表標(biāo)定摘要 本文對儲油罐的變位識別和罐容表的標(biāo)定問題進(jìn)行了深入探討，建立了儲油量和油位高度以及變位參數(shù)之間的數(shù)學(xué)模型，并利用matlab軟件對數(shù)據(jù)進(jìn)行處理。 對于問題一，我們利用積分法分別求出小橢圓儲油罐在無變位和發(fā)生縱向傾斜變位之間時的一般公式。并求出儲油罐縱向變位對罐容表的影響。分別給出了小橢圓儲油罐在無變位和縱向傾斜變位的罐容表。如下：表一 罐體變位前油位高度間隔為10cm的罐容表標(biāo)定值表油量高度/cm油量/L油量高度/cm油量/L油量高度/cm油量/L10163.59501621.0903306.620450.27602055.11003659.930803.63

2、702489.11103946.6401199.3802910.81204110.1表二 罐體變位后油位高度間隔為10cm的罐容表標(biāo)定值表油量高度/cm油量/L油量高度/cm油量/L油量高度/cm油量/L1070.1501371.9903072.420281.9601798.61003450.730595.3702232.51103776.740965.7802661.41204110.1針對問題二，我們將儲油罐分為5個區(qū)域進(jìn)行討論，在球罐出的運(yùn)算過于復(fù)雜，我們將球罐部分大致默認(rèn)為球缺，省略一小部分體積進(jìn)行了近似求解，再次得出關(guān)內(nèi)儲油量與油位高度以及變位參數(shù)一般關(guān)系的數(shù)學(xué)模型。 關(guān)鍵詞：積分求

3、解模型，區(qū)域分割法, 最小二乘法一、 問題重述通常加油站都有若干個儲存燃油的地下儲油罐，并且一般都有與之配套的“油位計量管理系統(tǒng)”，采用流量計和油位計來測量進(jìn)/出油量與罐內(nèi)油位高度等數(shù)據(jù)，通過預(yù)先標(biāo)定的罐容表（即罐內(nèi)油位高度與儲油量的對應(yīng)關(guān)系）進(jìn)行實時計算，以得到罐內(nèi)油位高度和儲油量的變化情況。許多儲油罐在使用一段時間后，由于地基變形等原因，使罐體的位置會發(fā)生縱向傾斜和橫向偏轉(zhuǎn)等變化（以下稱為變位），從而導(dǎo)致罐容表發(fā)生改變。按照有關(guān)規(guī)定，需要定期對罐容表進(jìn)行重新標(biāo)定。圖1是一種典型的儲油罐尺寸及形狀示意圖，其主體為圓柱體，兩端為球冠體。圖2是其罐體縱向傾斜變位的示意圖，圖3是罐體橫向偏轉(zhuǎn)變位的

4、截面示意圖。運(yùn)用數(shù)學(xué)建模方法研究解決儲油罐的變位識別與罐容表標(biāo)定的問題。 （1）為了掌握罐體變位后對罐容表的影響，利用如圖4的小橢圓型儲油罐（兩端平頭的橢圓柱體），分別對罐體無變位和傾斜角為a=4.10的縱向變位兩種情況做了實驗，實驗數(shù)據(jù)如附件1所示。請建立數(shù)學(xué)模型研究罐體變位后對罐容表的影響，并給出罐體變位后油位高度間隔為1cm的罐容表標(biāo)定值。（2）對于圖1所示的實際儲油罐，試建立罐體變位后標(biāo)定罐容表的數(shù)學(xué)模型，即罐內(nèi)儲油量與油位高度及變位參數(shù)（縱向傾斜角度a和橫向偏轉(zhuǎn)角度b ）之間的一般關(guān)系。請利用罐體變位后在進(jìn)/出油過程中的實際檢測數(shù)據(jù)（附件2），根據(jù)你們所建立的數(shù)學(xué)模型確定變位參數(shù)，并

5、給出罐體變位后油位高度間隔為10cm的罐容表標(biāo)定值。進(jìn)一步利用附件2中的實際檢測數(shù)據(jù)來分析檢驗?zāi)銈兡Ｐ偷恼_性與方法的可靠性。二 、模型假設(shè)1. 假設(shè)所給數(shù)據(jù)均為儲油罐內(nèi)壁測量值，且測量值均正確；2. 不考慮外界除角度因素以外的因素（如溫度，壓強(qiáng)，人為擠壓）引起儲油罐體積變化；3. 油位探針被固定在儲油罐上，油浮子始終垂直于儲油罐且測量結(jié)果準(zhǔn)確；三 、符號說明為了便于描述問題，我們用一些符號來代替問題中涉及的一些基本變量，如表3所示。 表三：主要符號說明a橢圓半長軸長b橢圓半短軸長L儲油罐總長l油位探針到油罐底部左側(cè)距離VT儲油罐總體積儲油罐縱向傾斜角度儲油罐橫向偏轉(zhuǎn)角度h油位高度R球冠體半徑

6、R0圓柱體的底面半徑注：未注明符號在文中提到時注明四、問題一模型建立利用積分的方法求出小橢圓儲油罐罐體變?yōu)榍昂蟮膬τ土亢陀臀桓叨戎g的關(guān)系。對于縱向傾斜的小橢圓儲油罐，考慮分段求出其儲油量和油位高度之間的關(guān)系，從而得到重新標(biāo)定后的罐容表。4.1 小橢圓儲油罐無變位模型 首先建立如圖所示坐標(biāo)系，橢圓半長軸長為a，半段軸長為b，以罐底為原點，建系，橢圓方程為，圖形如下圖所示：圖一 對橢圓的積分示意圖在y方向上取橢圓面中微元dy積分得到油的側(cè)面積為儲油罐內(nèi)油的體積為故利用matlab進(jìn)一步計算得根據(jù)此模型，我們求出小橢圓儲油罐無變位時罐容表表標(biāo)定值4.2 小橢圓儲油罐縱向傾斜變位時的模型儲油罐縱向變

7、為傾斜后，油位過高或過低都無法起作用，如下圖瑣事v1和v5區(qū)域，講儲油罐按液面高低分為五個部分，來球儲油量和油位高度之間的關(guān)系。如下圖二所示： 圖2 儲油罐分區(qū)示意圖4.2.1 對區(qū)域v1的討論在區(qū)域v1，其油位低于油位探針的油浮子，故油位顯示高度始終為0。將區(qū)域v1放大至圖3：圖3 區(qū)域v1放大圖如圖中，以原點（油罐正視圖左側(cè)頂點）紙面向里為x軸，利用三重積分得到其中l(wèi)為油探針到儲油罐左側(cè)的距離4.2.2 對區(qū)域v2的討論將區(qū)域V2放大至圖4，如下圖所示：圖4 區(qū)域v2示意圖 由于區(qū)域v1和區(qū)域v2的相似性，可以得到區(qū)域v2的儲油量和油位高度的變化關(guān)系，得出4.2.3 對區(qū)域v3的討論將區(qū)域

8、V3放大至圖5，如下圖所示： 圖5 區(qū)域v3示意圖將區(qū)域v3分割為兩部分，一部分為va另一部分為v，其中v又包含了vb，利用積分遠(yuǎn)離我們可一分別計算出va和vb，并可以得出油量為v3=va+v-vb其中4.2.4 對區(qū)域v4的討論將區(qū)域V4放大至圖6，如下圖所示：圖6 區(qū)域v4示意圖 利用總體積減去油罐內(nèi)無油部分的體積vb即可求出v4的體積其中4.2.5 對區(qū)域v5的討論 如下圖所示，由于此時油浮子到達(dá)油位探針的頂點，無法進(jìn)一步測量油位高度。故對v5不再進(jìn)行討論。圖7 區(qū)域v5示意圖4.2.6 綜合各區(qū)域的罐容表標(biāo)定的數(shù)學(xué)模型根據(jù)對各區(qū)域的討論我們可以得到罐體縱向傾斜變位（=4.1。）后油位高

9、度間隔為1cm的罐容表標(biāo)定值（見附錄一）。 當(dāng)=4.1。時，各區(qū)域油位高度及體積變化范圍為如下表所示：表四 各區(qū)域油位高度及體積變化范圍區(qū)域油位高度 h/m儲油量V/m31h=00.002520<h<=0.1470(0.0025,0.151430.1470<h<=1.1713(0.1514,3.958841.1713<h<=1.2(3.9588,4.0127 5h=1.2(4.0127,4.11014.3 罐體變位后對罐容表的影響為了能更準(zhǔn)確刻畫罐體的縱向傾斜變位對罐容表的影響，我們分別對罐體變?yōu)榍昂蟮睦碚撝岛蜏y量值進(jìn)行多方面的比較，所得結(jié)果分別如下圖所示：

10、 圖九4.3.1 罐體變位前理論值與測量值的比較根據(jù)附錄一中所給數(shù)據(jù)，我們計算出在附件所給的油位高度下的理論值和實際值，并劃出其曲線，如圖十所示：圖十同過對比我們發(fā)現(xiàn)對于任意h，儲油量理論值和實際值始終成比例，比例保持為0.9663.（具體代碼見附錄）4.3.2 罐體變位后理論值與實際值比較經(jīng)過MATLAB軟件的作圖，圖片如下所示：圖十一如圖所示測量值始終小于理論值，但保持在極小的誤差下。4.3.3 罐體縱向傾斜變位前后理論值的比較表五：變位前后理論值隨高度變化的差值表高度 h/cm變位后與變位前儲油量差值/L1093.4920168.3730208.3340233.650249.160256

11、.570256.680249.490234.2100209.2110169.91200作圖如下所示：由以上數(shù)據(jù)可得其平均影響率：算得結(jié)果如表六所示： 表六高度 h/cm變位后與變位前儲油量相對誤差100.5714897200.373931197300.259236216400.194780289500.153670574600.124811445700.10308947800.085680912900.0708280411000.0571600321100.0430497141200經(jīng)Matlab畫圖可的如圖十二所示： 圖十二通過以上數(shù)據(jù)得平均影響率為：五、問題二的模型建立如圖十三所示的儲油罐

12、示意圖，其主體為圓柱體，兩端為球冠體，在儲油罐無變位時我們計算其各部分的體積圖十三 圓柱體積計算公式一端球缺體積計算公式為5.1 考慮不發(fā)生變位時儲油量和油位高度的關(guān)系圓柱內(nèi)的油體積隨油位高度變化關(guān)系：一端球缺內(nèi)油體積隨油位高度變化關(guān)系：R為圓柱體地面半徑R0為球缺對應(yīng)半徑r圓為球內(nèi)小圓半徑總儲油量與油位高度變化關(guān)系為5.2 只考慮橫向偏轉(zhuǎn)變位是的儲油量和油位高度關(guān)系 只考慮儲油罐橫向偏轉(zhuǎn)為時：對實際的油位高度沒有影響，但此時的油位探針隨儲油罐發(fā)生偏轉(zhuǎn)圖十四 只考慮橫向偏轉(zhuǎn)示意圖由油浮子測量得到的油位高度與實際油位高度的關(guān)系為所以對罐容表影響轉(zhuǎn)換公式：5.3 只考慮縱向傾斜變位時的儲油量和油位

13、高度的關(guān)系 利用問題一中的思想，將儲油量分成5個區(qū)域，分別求儲油量和油位高度的關(guān)系，如圖10圖十五 實際儲油罐分區(qū)示意圖5.3.1 對區(qū)域1的討論油位探針測得的油位高度h始終未0，將其分成三部分來計算這部分體積，如區(qū)域放大圖11圖十六 區(qū)域1放大圖區(qū)域1的總體積應(yīng)為：對圓柱體三重積分得到：球缺部分計算公式得到由于Va部分體積對于整體來說較小，故可忽略5.3.2 對區(qū)域2討論圖十七 區(qū)域2示意圖 在區(qū)域2種，油位高度h變化范圍為0，（L-l）tan各部分儲油體積和油位高度的變化關(guān)系：5.3.3 對區(qū)域3的討論圖十八 區(qū)域3示意圖 如圖，在區(qū)域3中，油位高度h的變化范圍（L-l）tan，2R-lt

14、an儲油量為利用問題一的類似解法5.3.4 對區(qū)域4的討論圖十七 區(qū)域4示意圖 在區(qū)域4中，油位高度h的范圍2R-ltan,2R.儲油量與油位高度的關(guān)系為其中5.3.5 對區(qū)域v5的討論 如下圖所示，由于此時油浮子到達(dá)油位探針的頂點，無法進(jìn)一步測量油位高度。故對v5不再進(jìn)行討論。圖十八 區(qū)域5示意圖5.4 綜合考慮儲油罐縱向傾斜和橫向偏轉(zhuǎn)根據(jù)如上討論，把5.3中各區(qū)域V-h公式中的h替換為即可得到儲油罐縱向傾斜和橫向偏轉(zhuǎn)的V（h，）一般關(guān)系式。依然考慮將儲油罐分成五個區(qū)域進(jìn)行求解5.4.1 對區(qū)域1的討論當(dāng)發(fā)生橫向偏轉(zhuǎn)后，油位探針測得油位高度為0，h=R-Rcos，故本屬于第二區(qū)域的部分橫向偏

15、轉(zhuǎn)后進(jìn)入了第一區(qū)域，使得l變?yōu)閘+故只需將5.3中的V1中l(wèi)替換為 l+故5.4.2 對區(qū)域2的討論由可以得到h=R-（R-h）cos在第二區(qū)域又要滿足計算得0然后只需做下式變換即可得到V（h，）一般關(guān)系式V2（h）=V2（R-）5.4.3 對v3區(qū)域的討論由可以得到h=R-（R-h）cos在第三區(qū)域又要滿足計算得然后只需做下式變換即可得到V（h，）一般關(guān)系式V3（h）=V3（R-）5.4.4 對v4區(qū)域的討論由可以得到h=R-（R-h）cos在第四區(qū)域又要滿足2R然后只需做下式變換即可得到V（h，）一般關(guān)系式V4（h）=V4（R-）5.4.5 對區(qū)域v5的討論 由于此時油浮子到達(dá)油位探針的頂

16、點，無法進(jìn)一步測量油位高度。故對v5不再進(jìn)行討論。5.5 變位參數(shù)的確定和罐容表的計算5.5.1 用最小二乘參數(shù)估計法確定參數(shù) 最小二乘參數(shù)估計法：根據(jù)V（h，）的關(guān)系表達(dá)式求得幾組油量高度h，計算出相鄰高度油量體積之差，通過與附件的實際儲油量進(jìn)行比較，通過對，進(jìn)行等間距窮舉最終求得理論值與實際值的差值平方和S，當(dāng)S取最小是時，為所求最佳值，求解公式如下用最小二乘參數(shù)估計法得到變位參數(shù)為：=2.51°，=4.85°，角度符合實際情況。5.5.2 實際儲油罐罐容表制定估計出參數(shù)，后，我們可以根據(jù)所建立的儲油量V以及油位高度h以及變位參數(shù)，的一般模型得到罐容表表七 實際儲油罐罐

17、容表（縱向變位=2.51°，橫向偏轉(zhuǎn)=4.85°）油位高度h/m0.10.20.30.40.50.60.70.80.91.0儲油量V/m30.31960.93026.17111.0911.69312.42713.29414.29515.43416.71油位高度h/m1.11.21.31.41.51.61.71.81.92.0儲油量V/m318.12419.67721.3723.1925.1627.2529.4731.8134.2736.83油位高度h/m2.12.22.32.42.52.62.72.82.93.0儲油量V/m339.5142.2745.1248.0551.0

18、354.0757.1359.2160.2761.5625.6數(shù)據(jù)檢驗根據(jù)推導(dǎo)數(shù)據(jù)V（h，）關(guān)系式對第一次補(bǔ)充僅有后的數(shù)據(jù)進(jìn)行檢驗，得到相同油位高度下，理論儲油值與實際儲油值關(guān)系圖： 圖十九 理論儲油值與實際儲油值量差別對比圖圖中顯示理論值與實際值基本重合，表示算的，較為合理。六、模型評價1. 明確給出變位參數(shù)的辨識準(zhǔn)則，由簡單到復(fù)雜建立模型，邏輯性強(qiáng)2. 數(shù)據(jù)參數(shù)化，是模型更具有靈活性，應(yīng)變能力高。3. 運(yùn)用matlab程序加強(qiáng)數(shù)據(jù)的科學(xué)性與準(zhǔn)確性。七、模型改進(jìn)1. 由于公式復(fù)雜，程序運(yùn)行較慢，將模型參數(shù)進(jìn)行近似簡化，會減少模型運(yùn)行時間；2. 將模型公式進(jìn)行更細(xì)致的分割，分開運(yùn)算各部分，來減少

19、運(yùn)行的錯誤以及程序運(yùn)行時間。參考文獻(xiàn) 1陳光亭 裘哲勇 數(shù)學(xué)建模高等教育出版社 2010八、附錄第一部分：1.1、畫出變位前理論值與測量值的曲線圖：h=linspace(0,1.2,78);m= linspace(0,1200,78);>> b=0.6;>> a=1.78/2;>> L=2.45;>> for i=1:78v(i)=1000*a*b*L*(h(i)-b)/(b2)*sqrt(h(i)*(2*b-h(i)+asin(h(i)-b)/b)+pi/2);end>> 

20、plot(m,v,'r*')>> hold on>> xlabel('油位高度 h/mm');>> ylabel('儲油量/L')>> A=xlsread('C:UsershfDesktop問題A附件1：實驗采集數(shù)據(jù)表.xls','G2:G79');>> B=xlsread('C:UsershfDesktop問題A附件1：實驗采集數(shù)據(jù)表.xls','D2:D79&

21、#39;);>> plot(B,A,'b-');legend('理論值','實際值');title('罐體變位前的 V-h曲線對比')1.2、得出上述變位前理論值與測量值的曲線圖中理論與測量的比值：A=xlsread('C:UsershfDesktop問題A附件1：實驗采集數(shù)據(jù)表.xls','D2:D79');C=A./1000;b=0.6;a=1.78/2;L=2.45;for i=1:78V(i)= a*b*L*(C(i)-b)/(b2)*sqrt(C(i)*(2*b-C(

22、i)+asin(C(i)-b)/b)+pi/2)*1000;EndD=xlsread('C:UsershfDesktop問題A附件1：實驗采集數(shù)據(jù)表.xls','G2:G79');for i=1:78T(i)=D(i)/V(i);end·1.3 正常情況下小橢圓罐的罐容表代碼：h=linspace(0.1,1.2,12);>> b=0.6;>> a=1.78/2;>> L=2.45;>> for i=1:12v(i)=1000*a*b*L*(h(i)-b)

23、/(b2)*sqrt(h(i)*(2*b-h(i)+asin(h(i)-b)/b)+pi/2);end結(jié)果：1.4、變位后的相關(guān)公式積分結(jié)果：V1:c =tan(4.1*pi/180); l=0.4; d=c*l;L=2.45;syms h a b x y z;2*int(int(int(1,z,0,y/c),x,0,a*sqrt(1-(y/b-1)2),y,0,d)H=0, V1= 0.0025V2:c =tan(4.1*pi/180); l=0.4; d=c*l;L=2.45;syms h a b x y z;int(int(int(1,z,0,(h+d-y)/c),x,0,a*sqrt(

24、1-(y/b-1)2),y,0,h+d)*2h=0.1470 V=0.1514V3:c =tan(4.1*pi/180); l=0.4; d=c*l;L=2.45;syms h a b x y z;V=a*b*L*(h/(b2)-1/b)*sqrt(2*b*h-h2)+asin(h/b-1)+pi/2);Va=int(int(int(1,z,0,(h+d-y)/c),x,0,a*sqrt(1-(y/b-1)2),y,h,h+d)*2;Vb=2*int(int(int(1,z,(h-y+d)/c,L),x,0,a*sqrt(1-(y/b-1)2),y,h-(L-l)*c,h);V3=V+Va-V

25、bH=1.1713 V= 3.9588V4:c =tan(4.1*pi/180); l=0.4; d=c*l;L=2.45;syms h a b x y z;V4=a*b*L*(h-b)/(b2)*sqrt(h*(2*b-h)+asin(h-b)/b)+pi/2)-2*int(int(int(1,z,(h-y+d)/c,L),x,0,a*sqrt(1-(y/b-1)2),y,h-(L-l)*c,2*b);H=1.2 V4= 4.0127區(qū)域油位高度 h/m儲油量V/m31h=00.002520<h<=0.1470(0.0025,0.151430.1470<h<=1.17

26、13(0.1514,3.958841.1713<h<=1.2(3.9588,4.0127 5h=1.2(4.0127,4.11011.5、畫出變位后測量值與理論值的曲線：得到V2區(qū)間的相對應(yīng)點的體積值：h=linspace(0,0.1470,3);for i=1:3B(i)=vpa(Untitled2(h(i),5);endB = 1.6744, 39.597, 151.37畫圖：h=linspace(0,0.1470,3);B = 1.6744, 39.597, 151.37H=h*1000;>> plot(H,B,'r*')>>

27、 hold on得到V3區(qū)間的相對應(yīng)點的體積值：V = 0.45027, 0.80354, 1.1993, 1.621, 2.0551, 2.4891, 2.9108, 3.3066, 3.6599, 3.9466;Va = 0.007749, 0.0089308, 0.0096788, 0.01009, 0.010205, 0.010035, 0.0095638, 0.0087429, 0.0074578, 0.0053847;Vb = 0.17616, 0.21722, 0.24332, 0.25921, 0.26675, 0.26668, 0.25898, 0.242

28、92, 0.21661, 0.1753V3 = 1.0e+03 * Columns 1 through 7 0.2819 0.5953 0.9657 1.3719 1.7986 2.2325 2.6614 Columns 8 through 10 3.0724 3.4507 3.7767畫圖：s= linspace(0.2,1.1,10);V3=0.2819,0.5953,0.9657,1.3719,1.7986 ,2.2325,2.66143.0724 ,3.4507,3.7767; S=s*1000; V3=V3*1000;>> plot(S,V3,'r*')&

29、gt;>hold on得到V4區(qū)間的相對應(yīng)點的體積值：t=linspace(1.1716,1.2,6);for i=1:6D(i)=Untitled4(h(i);end D=vpa(D,5)畫圖：D= 3.9342, 3.9531, 3.9707, 3.9867, 4.0009, 4.0127t=linspace(1.1716,1.2,6);T=t*1000;D=D*1000;>> plot(T,D,'r*')根據(jù)附錄一畫出測量值：>> xlabel('油位高度 h/mm');>> 

30、;ylabel('儲油量/L')>> Y=xlsread('C:UsershfDesktop問題A附件1：實驗采集數(shù)據(jù)表.xls','傾斜變位進(jìn)油','G2:G79');>> Z=xlsread('C:UsershfDesktop問題A附件1：實驗采集數(shù)據(jù)表.xls','傾斜變位進(jìn)油','D2:D79');>> plot(Z,Y,'b-');legend('理論值','實際值');

31、title('罐體變位后的 V-h曲線對比')1.6、變位前后儲油罐理論值的比較：A=xlsread('C:UsershfDesktop問題A附件1：實驗采集數(shù)據(jù)表.xls','G2:G79');B=xlsread('C:UsershfDesktop問題A附件1：實驗采集數(shù)據(jù)表.xls','D2:D79');plot(B,A,'b-');>> hold on>> h=linspace(0,0.1470,3);B = 1.6744, 39.597, 151.37H=h*1000

32、;plot(H,B,'r*')B = 1.6744 39.5970 151.3700>> hold on>> s= linspace(0.2,1.1,10);V3=0.2819,0.5953,0.9657,1.3719,1.7986 ,2.2325,2.6614 3.0724 ,3.4507,3.7767; S=s*1000; V3=V3*1000;plot(S,V3,'r*')>> hold on>> D= 3.9342, 3.9531, 3.9707, 3.9867, 4.0009, 4.0127t=linsp

33、ace(1.1716,1.2,6);T=t*1000;D=D*1000;plot(T,D,'r*')D = 3.9342 3.9531 3.9707 3.9867 4.0009 4.0127>> hold on>> legend('變位前','變位后');title('罐體變位前后理論值的 V-h 曲線對比')xlabel('油位高度 h/mm');ylabel('儲油量/L')1.7、變位前后理論值隨高度變化的差值由上述罐體變位前油位高度間隔為10cm的罐容表標(biāo)定值表以及罐

34、體變位后油位高度間隔為10cm的罐容表標(biāo)定值表，可得相應(yīng)的差值為：變位前后理論值隨高度變化的差值表 高度 h/cm變位后與變位前儲油量差值/L1093.4920168.3730208.3340233.650249.160256.570256.680249.490234.2100209.2110169.91200高度 h/cm變位后與變位前儲油量相對誤差100.5714897200.373931197300.259236216400.194780289500.153670574600.124811445700.10308947800.085680912900.0708280411000.0571

35、600321100.0430497141200第二部分代碼：2.1 根據(jù)最小二乘法以及窮舉法，算得=2.51。，=4.85。以下為相應(yīng)的計算不同區(qū)域的體積的代碼：syms x y z n l m h;r=sqrt(1.52-(y-1.5)2);h =1.5-(1.5-h)/n;l =(1.5-1.5*n)/m+l;V球缺：syms h y R0 R Hr =sqrt(R02-(h-R)2);int(r2*acos(R0-H)/r) -sqrt(r2-(R0-H)2)*(R0-H),h,0,h)ans =int(H - R0)*(R02 - (H - R0)2 - (R - h)2)(1/2)

36、+ (R02 - (R - h)2)*(pi - acos(H - R0)/(R02 - (R - h)2)(1/2), h, 0, h)>> R0=1.625;H=1;R=1.5;>> int(H - R0)*(R02 - (H - R0)2 - (R - h)2)(1/2) + (R02 - (R - h)2)*(pi - acos(H - R0)/(R02 - (R - h)2)(1/2), h, 0, h)ans =int(- (5*(9/4 - (h - 3/2)2)(1/2)/8 - acos(5/(8*(169/64 - (h - 3/2)2)(1/2)*

37、(h - 3/2)2 - 169/64), h, 0, h)>> h=0:0.004:2.9for i=1:length(h)z(i)= - (5*(9/4 - (h(i) - 3/2)2)(1/2)/8 - acos(5/(8*(169/64 - (h(i) - 3/2)2)(1/2)*(h(i) - 3/2)2 - 169/64);endsum(z(:)V2:syms x y z n l m h;r=sqrt(1.52-(y-1.5)2);Vc=int(int(int(1,z,0,(h+l*m-y)/m),x,0,r),y,0,h+l*m)*2;Vb=int(-acos(5/(

38、8*(169/64 - (y - 3/2)2)(1/2)*(y - 3/2)2 - 169/64) - (5*(9/4 - (y - 3/2)2)(1/2)/8, y, 0, h)*(h+l*m);>> V2=Vb+VcV2 =int(- acos(5/(8*(169/64 - (y - 3/2)2)(1/2)*(y - 3/2)2 - 169/64) - (5*(9/4 - (y - 3/2)2)(1/2)/8, y, 0, h)*(h + l*m) + 2*int(l + (h - y)/m)*(9/4 - (y - 3/2)2)(1/2), y, 0, h + l*m)>

39、;> h =1.5-(1.5-h)/n; l =(1.5-1.5*n)/m+l;>> int(- acos(5/(8*(169/64 - (y - 3/2)2)(1/2)*(y - 3/2)2 - 169/64) - (5*(9/4 - (y - 3/2)2)(1/2)/8, y, 0, h)*(h + l*m) + 2*int(l + (h - y)/m)*(9/4 - (y - 3/2)2)(1/2), y, 0, h + l*m)ans =2*int(9/4 - (y - 3/2)2)(1/2)*(l + (h - 3/2)/n - y + 3/2)/m - (3*n)

40、/2 - 3/2)/m), y, 0, m*(l - (3*n)/2 - 3/2)/m) + (h - 3/2)/n + 3/2) + int(- acos(5/(8*(169/64 - (y - 3/2)2)(1/2)*(y - 3/2)2 - 169/64) - (5*(9/4 - (y - 3/2)2)(1/2)/8, y, 0, (h - 3/2)/n + 3/2)*(m*(l - (3*n)/2 - 3/2)/m) + (h - 3/2)/n + 3/2)syms x y z n l m h;2*int(9/4 - (y - 3/2)2)(1/2)*(l + (h - 3/2)/n

41、- y + 3/2)/m - (3*n)/2 - 3/2)/m), y, 0, m*(l - (3*n)/2 - 3/2)/m) + (h - 3/2)/n + 3/2) + int(- acos(5/(8*(169/64 - (y - 3/2)2)(1/2)*(y - 3/2)2 - 169/64) - (5*(9/4 - (y - 3/2)2)(1/2)/8, y, 0, (h - 3/2)/n + 3/2)*(m*(l - (3*n)/2 - 3/2)/m) + (h - 3/2)/n + 3/2)ans =2*int(9/4 - (y - 3/2)2)(1/2)*(l + (h - 3

42、/2)/n - y + 3/2)/m - (3*n)/2 - 3/2)/m), y, 0, m*(l - (3*n)/2 - 3/2)/m) + (h - 3/2)/n + 3/2) + int(- acos(5/(8*(169/64 - (y - 3/2)2)(1/2)*(y - 3/2)2 - 169/64) - (5*(9/4 - (y - 3/2)2)(1/2)/8, y, 0, (h - 3/2)/n + 3/2)*(m*(l - (3*n)/2 - 3/2)/m) + (h - 3/2)/n + 3/2)>> l=2;>> 2*int(9/4 - (y -

43、3/2)2)(1/2)*(l + (h - 3/2)/n - y + 3/2)/m - (3*n)/2 - 3/2)/m), y, 0, m*(l - (3*n)/2 - 3/2)/m) + (h - 3/2)/n + 3/2) + int(- acos(5/(8*(169/64 - (y - 3/2)2)(1/2)*(y - 3/2)2 - 169/64) - (5*(9/4 - (y - 3/2)2)(1/2)/8, y, 0, (h - 3/2)/n + 3/2)*(m*(l - (3*n)/2 - 3/2)/m) + (h - 3/2)/n + 3/2)ans =2*int(9/4 -

44、 (y - 3/2)2)(1/2)*(h - 3/2)/n - y + 3/2)/m - (3*n)/2 - 3/2)/m + 2), y, 0, (h - 3/2)/n - m*(3*n)/2 - 3/2)/m - 2) + 3/2) + int(- acos(5/(8*(169/64 - (y - 3/2)2)(1/2)*(y - 3/2)2 - 169/64) - (5*(9/4 - (y - 3/2)2)(1/2)/8, y, 0, (h - 3/2)/n + 3/2)*(h - 3/2)/n - m*(3*n)/2 - 3/2)/m - 2) + 3/2)V3syms x y z n

45、 l m h;r=sqrt(1.52-(y-1.5)2);Vc1=int(int(int(1,z,0,(h+l*m-y)/m),x,0,r),y,0,h+l*m)*2;Vc2=2*int(8*r,y,0,h);Vc3=2*int(int(int(1,z,(h-y+l*m)/m,8),x,0,r),y,h-(8-l)*m,h);>> VQ1=int(-acos(5/(8*(169/64 - (y - 3/2)2)(1/2)*(y - 3/2)2 - 169/64) - (5*(9/4 - (y - 3/2)2)(1/2)/8, y, 0, h)*(h+l*m);>> VQ

46、2=int(-acos(5/(8*(169/64 - (y - 3/2)2)(1/2)*(y - 3/2)2 - 169/64) - (5*(9/4 - (y - 3/2)2)(1/2)/8, y, 0, h)*(h-(8-l)*m);>> V3=Vc1+Vc2-Vc3+VQ1+VQ2V3 =9*pi + 18*asin(2*h)/3 - 1) + 8*(9/4 - (h - 3/2)2)(1/2)*(h - 3/2) + int(- acos(5/(8*(169/64 - (y - 3/2)2)(1/2)*(y - 3/2)2 - 169/64) - (5*(9/4 - (y -

47、 3/2)2)(1/2)/8, y, 0, h)*(h + l*m) - 2*int(-(9/4 - (y - 3/2)2)(1/2)*(l + (h - y)/m - 8), y, h + m*(l - 8), h) + int(- acos(5/(8*(169/64 - (y - 3/2)2)(1/2)*(y - 3/2)2 - 169/64) - (5*(9/4 - (y - 3/2)2)(1/2)/8, y, 0, h)*(h + m*(l - 8) + 2*int(l + (h - y)/m)*(9/4 - (y - 3/2)2)(1/2), y, 0, h + l*m)h =1.5

48、-(1.5-h)/n; l =(1.5-1.5*n)/m+l;>> 9*pi + 18*asin(2*h)/3 - 1) + 8*(9/4 - (h - 3/2)2)(1/2)*(h - 3/2) + int(- acos(5/(8*(169/64 - (y - 3/2)2)(1/2)*(y - 3/2)2 - 169/64) - (5*(9/4 - (y - 3/2)2)(1/2)/8, y, 0, h)*(h + l*m) - 2*int(-(9/4 - (y - 3/2)2)(1/2)*(l + (h - y)/m - 8), y, h + m*(l - 8), h) + i

49、nt(- acos(5/(8*(169/64 - (y - 3/2)2)(1/2)*(y - 3/2)2 - 169/64) - (5*(9/4 - (y - 3/2)2)(1/2)/8, y, 0, h)*(h + m*(l - 8) + 2*int(l + (h - y)/m)*(9/4 - (y - 3/2)2)(1/2), y, 0, h + l*m)ans = 9*pi + 18*asin(2*(h - 3/2)/(3*n) + 2*int(9/4 - (y - 3/2)2)(1/2)*(l + (h - 3/2)/n - y + 3/2)/m - (3*n)/2 - 3/2)/m)

50、, y, 0, m*(l - (3*n)/2 - 3/2)/m) + (h - 3/2)/n + 3/2) + int(- acos(5/(8*(169/64 - (y - 3/2)2)(1/2)*(y - 3/2)2 - 169/64) - (5*(9/4 - (y - 3/2)2)(1/2)/8, y, 0, (h - 3/2)/n + 3/2)*(h - 3/2)/n - m*(3*n)/2 - 3/2)/m - l + 8) + 3/2) + int(- acos(5/(8*(169/64 - (y - 3/2)2)(1/2)*(y - 3/2)2 - 169/64) - (5*(9/

51、4 - (y - 3/2)2)(1/2)/8, y, 0, (h - 3/2)/n + 3/2)*(m*(l - (3*n)/2 - 3/2)/m) + (h - 3/2)/n + 3/2) - 2*int(-(9/4 - (y - 3/2)2)(1/2)*(l + (h - 3/2)/n - y + 3/2)/m - (3*n)/2 - 3/2)/m - 8), y, (h - 3/2)/n - m*(3*n)/2 - 3/2)/m - l + 8) + 3/2, (h - 3/2)/n + 3/2) + (8*(9/4 - (h - 3/2)2/n2)(1/2)*(h - 3/2)/nsy

52、ms x y z n l m h;9*pi + 18*asin(2*(h - 3/2)/(3*n) + 2*int(9/4 - (y - 3/2)2)(1/2)*(l + (h - 3/2)/n - y + 3/2)/m - (3*n)/2 - 3/2)/m), y, 0, m*(l - (3*n)/2 - 3/2)/m) + (h - 3/2)/n + 3/2) + int(- acos(5/(8*(169/64 - (y - 3/2)2)(1/2)*(y - 3/2)2 - 169/64) - (5*(9/4 - (y - 3/2)2)(1/2)/8, y, 0, (h - 3/2)/n

53、+ 3/2)*(h - 3/2)/n - m*(3*n)/2 - 3/2)/m - l + 8) + 3/2) + int(- acos(5/(8*(169/64 - (y - 3/2)2)(1/2)*(y - 3/2)2 - 169/64) - (5*(9/4 - (y - 3/2)2)(1/2)/8, y, 0, (h - 3/2)/n + 3/2)*(m*(l - (3*n)/2 - 3/2)/m) + (h - 3/2)/n + 3/2) - 2*int(-(9/4 - (y - 3/2)2)(1/2)*(l + (h - 3/2)/n - y + 3/2)/m - (3*n)/2 -

54、 3/2)/m - 8), y, (h - 3/2)/n - m*(3*n)/2 - 3/2)/m - l + 8) + 3/2, (h - 3/2)/n + 3/2) + (8*(9/4 - (h - 3/2)2/n2)(1/2)*(h - 3/2)/nans =9*pi + 18*asin(2*h - 3)/(3*n) + 2*int(9/4 - (y - 3/2)2)(1/2)*(l + (h - 3/2)/n - y + 3/2)/m - (3*n)/2 - 3/2)/m), y, 0, m*(l - (3*n)/2 - 3/2)/m) + (h - 3/2)/n + 3/2) + int(- acos(5/(8*(169/64 - (y - 3/2)2)(1/2)*(y - 3/2)2 - 169/64) - (5*(9/4 - (y - 3/2)2)(1/2)/8, y, 0, (h - 3/2)/