储油罐地变位识别与罐容表标定

1、储油罐的变位识别与罐容表标定摘要本文对储油罐的变位识别和罐容表的标定问题进行了深入探讨，建立了储油量和 油 位高度以及变位参数之间的数学模型，并利用matlab软件对数据进行处理。对于问题一，我们利用积分法分别求出小椭圆储油罐在无变位和发生纵向倾斜变 位 之间时的一般公式。并求出储油罐纵向变位对罐容表的影响。分别给出了小椭圆 储油罐 在无变位和纵向倾斜变位的罐容表。如下：表一罐体变位前油位高度间隔为10cm的罐容表标定值表油量高度/cm油量/L油量高度/cm油量/L油量高度/cm油量/L10.59501621.0903306.620450.27602055.11003659.930803.63

2、702489.11103946.6401199.3802910.81204110.1表二 罐体变位后油位高度间隔为10cm的罐容表标定值表油量高度/cm油量/L油量高度/cm油量/L油量高度/cm油量/L1070.1501371.9903072.420281.9601798.61003450.730595.3702232.51103776.740965.7802661.41204110.1针对问题二，我们将储油罐分为 5个区域进行讨论，在球罐出的运算过于复杂, 我们 将球罐局部大致默认为球缺，省略一小局部体积进行了近似求解，再次得出关 储油量与 油位高度以及变位参数一般关系的数学模型。关键词：

3、积分求解模型，区域分割法,最小二乘法问题重述通常加油站都有假设干个储存燃油的地下储油罐，并且一般都有与之配套的“油位计量 管理系统，采用流量计和油位计来测量进 /出油量与罐油位高度等数据，通过 预先标定 的罐容表即罐油位高度与储油量的对应关系进行实时计算，以得到罐 油位高度和储 油量的变化情况。许多储油罐在使用一段时间后，由于地基变形等原因，使罐体的位置会发生纵向倾斜和横向偏转等变化以下称为变位，从而导致罐容表发生改变。按照有关规定， 需要 定期对罐容表进行重新标定。图 1是一种典型的储油罐尺寸及形状示意图，其 主体为圆 柱体，两端为球冠体。图2是其罐体纵向倾斜变位的示意图，图 3是罐体 横向

4、偏转变位 的截面示意图。运用数学建模方法研究解决储油罐的变位识别与罐容表标定的问题。1 为了掌握罐体变位后对罐容表的影响，利用如图4的小椭圆型储油罐两端平 头的椭圆柱体，分别对罐体无变位和倾斜角为=4.10的纵向变位两种情况做了实验，实验数据如附件1所示。请建立数学模型研究罐体变位后对罐容表的影响， 并给出罐体 变位后油位高度间隔为1cm的罐容表标定值。2对于图1所示的实际储油罐，试建立罐体变位后标定罐容表的数学模型，即 罐 储油量与油位高度及变位参数纵向倾斜角度和横向偏转角度之间的一般关系。请利用罐体变位后在进/出油过程中的实际检测数据附件 2，根据你们所 建立的数学 模型确定变位参数，并给

5、出罐体变位后油位高度间隔为10cm的罐容表标定值。进一步 利用附件2中的实际检测数据来分析检验你们模型的正确性与方法的可靠性。模型假设1. 假设所给数据均为储油罐壁测量值，且测量值均正确；2. 不考虑外界除角度因素以外的因素如温度，压强，人为挤压引起储油罐体 积变化；3. 油位探针被固定在储油罐上，油浮子始终垂直于储油罐且测量结果准确；、符号说明 为了便于描述问题，我们用一些符号来代替问题中涉及的一些根本变量，如表所示表三：主要符号说明a椭圆半长轴长b椭圆半短轴长L储油罐总长l油位探针到油罐底部左侧距离Vt储油罐总体积a储油罐纵向倾斜角度储油罐横向偏转角度h油位高度R球冠体半径Ro圆柱体的底面

6、半径注：未注明符号在文中提到时注明四、问题一模型建立利用积分的方法求出小椭圆储油罐罐体变为前后的储油量和油位高度之间的关系对于纵向倾斜的小椭圆储油罐，考虑分段求出其储油量和油位高度之间的关系，从而得到重新标定后的罐容表。4.1小椭圆储油罐无变位模型首先建立如下图坐标系，椭圆半长轴长为a，半段轴长为b，以罐底为图一对椭圆的积分示意dy原点，建系，椭圆方程为，图形如下列图所示:储油罐油的体积为利用matlab进一步计算得根据此模型，我们求出小椭圆储油罐无变位时罐容表表标定值4.2小椭圆储油罐纵向倾斜变位时的模型储油罐纵向变为倾斜后，油位过高或过低都无法起作用，如下列图琐事V1和V5区域，讲 储油罐

7、按液面上下分为五个局部，来球储油量和油位高度之间的关系。如下 图二所示：图2储汕罐分区示意图4.2.1 对区域vi的讨论0。将区域在区域vi，其油位低于油位探针的油浮子，故油位显示高度始终为v1 放大至图3:图3区域vi放大图如图中，以原点油罐正视图左侧顶点纸面向里为x轴，利用三重积分得到其中I为油探针到储油罐左侧的距离422对区域V2的讨论将区域V2放大至图4，如下列图所示:图4区域也示意图和区域V2的相似性，可以得到区域V2的储油量和油位高度的变化关系 得出423对区域V3的讨论将区域V3放大至图5，如下列图所示:械探针图5区域w示意I冬I将区域V3分割为两局部，一局部为 v另一局部为v，

8、其中v又包含了 Vb，利用 积分 远离我们可一分别计算出 Va和Vb，并可以得出油量为V3=Va+V -V b其中424对区域V4的讨论将区域V4放大至图6,如下列图所示:利用总体积减去油罐无油局部的体积Vb即可求出V4的体积其中425对区域V5的讨论如下列图所示，由于此时油浮子到达油位探针的顶点，无法进 一步测量油位高度。故 对V5不再进行讨论。图7区域vs示意图426综合各区域的罐容表标定的数学模型根据对各区域的讨论我们可以得到罐体纵向倾斜变位=4.1 ° 后油位高度间隔为1cm的罐容表标定值见附录一。当=4.1 °时，各区域油位高度及体积变化围为如下表所示：表四 各区

9、域油位高度及体积变化围区域油位高度h/m储油量V/m31h=00.002520<h<=0.1470(0.0025,0.1514130.1470<h<=1.1713(0.1514,3.958841.1713<h<=1.2(3.9588,4.0127 5h=1.2(4.0127,4.11014.3罐体变位后对罐容表的影响为了能更准确刻画罐体的纵向倾斜变位对罐容表的影响，我们分别对罐体变为前后的理论值和测量值进行多方面的比拟，所得结果分别如下列图所示：图九罐体变位前理论值与测量值的比拟 根据附录一中所给数据，我们计算出在附件 所给的油位高度下的理论值和实际值，并戈

10、U出其曲线，如图十所示：同过比照我们发现对于任意h，储油量理论值和实际值始终成比例，比例保持为0.9663.具体代码见附录432罐体变位后理论值与实际值比拟经过MATLAB软件的作图，图片如下所示:n J+fcL/ </ *p K1 M>-1 200400f 兀303：* 2C 'Jmm4500图十一如下图测量值始终小于理论值，但保持在极小的误差下罐体纵向倾斜变位前后理论值的比拟表五：变位前后理论值随高度变化的差值表高度 h/cm变位后与变位前储油量差值/L1093.4920168.3730208.3340233.650249.160256.570256.680249.49

11、0234.2100209.2110169.91200作图如下所示:由以上数据可得其平均影响率：n V 变位前V 变位后iOnV 总算得结果如表六所示:表六高度h/cm变位后与变位前储油量相对误差1100.5714897200.373931197300.259236216400.194780289500.153670574600.124811445700.10308947800.085680912900.0708280411000.0571600321100.0430497141200经 Matlab画图可的如图十二所示:图十二通过以上数据得平均影响率为：4.7215%五、问题二的模型建立如图十

12、三所示的储油罐示意图，其主体为圆柱体，两端为球冠体，在储油罐 无变位时我们计算其各局部的体积V圆柱L.R圆柱体积计算公式一端球缺体积计算公式为V球缺.H2 .(Ro H3)35.1考虑不发生变位时储油量和油位高度的关系圆柱的油体积随油位高度变化关系：hV 圆柱(h) 2 L. R2 (h R)2dhoV 球缺(h) S.dharccos(Ro-Hfclr圆(-Ro -H)( . Ro H)dh一端球缺油体积随油位高度变化关系:r 圆Ro2 ( - h-R) 2R为圆柱体地面半径 Ro为球缺对应半径 r圆为球小圆半径总储油量与油位高度变化关系为V(h) 2V 球缺(h) V 圆柱(h)5.2只考

13、虑横向偏转变位是的储油量和油位高度关系只考虑储油罐横向偏转为时：对实际的油位高度没有影响，但此时的油位探由忙捋针针随储油罐发生偏转图十四只考虑横向偏转示意图由油浮子测量得到的油位高度与实际油位高度的关系为RhCOS所以B对RhV(h ) V(R )cos5.3只考虑纵向倾斜变位时的储油量和油位高度的关系利用问题一中的思想，将储油量分成5个区域，分别求储油量和油位高度的关系，如图10图十五 实际储油罐分区示意图对区域1的讨论 油位探针测得的油位高度 h始终未0，将其分成三部 分来计算这局部体积，如区域放大图11图十六区域1放大图区域1的总体积应为:V1 Va Vb Vc对圆柱体三重积分得到:2

14、2 yI tanR2 ( y R)2 tanVc 2 dy dx dz 0 0 0球缺局部计算公式得到Vb V 球缺(I tan )由于Va局部体积对于整体来说较小，故可忽略Va 0532在区域2种，油位高度h变化围为0 ,( L-I ) tan a|各局部储油体积和 油位高度的变化关系：h I tan yh I tanR2 ( y R ) 2tanVc 2 dy dx dz0 0 0Vb V 球缺(h Itan )Va 0对区域3的讨论如图, 储油量为在区域3匚图十八区域3示意图，油位高度h的变化围Ltan a , 2R-ltanV3 Vc1 V C2 Vc3V球缺1 V球缺2利用问题一的类

15、似解I tanR (y R)2dydxh y I tantandzhR2 (h R)hVC2dxdydz 2 dhdx dy0 02 L R2 (h R)2dh0VcioR2 ( y R)2LV ,dyC3dxdz(L l)tan0hltantanV球缺i V球缺(h I tanV球缺2 V球缺 L ltan对区域/的讨论az图十七区域4示意图在区域4中，油位高度h的围2R-Itan a2R.储油量与油位高度的关系为V4 Vt VbV 球缺其中VtR2LV 球缺 (2R)2R R2 ( y R)2LVb2 dy dx dz h (L I )ta n 0 h y I tan tanV球缺V球缺h

16、 (L l)tan 对区域V5的讨论如下列图所示，由于此时油浮子到达油位探针的顶点，无法进一步测 量油位高度。故 对V5不再进行讨论。图十八区域5示意I冬5.4综合考虑储油罐纵向倾斜和横向偏转根据如上讨论，把5.3中各区域V-h公式中的h替换为R R h即可得到储油罐纵cos向倾斜和横向偏转的V( h, a, ® 般关系式。依然考虑将储油罐分成五个区域进行求解对区域1的讨论当发生横向偏转B后，油位探针测得油位高度为0, h=R-Rcos 3,故本属于第二区rR在第二区域又要满h=R-costanR Rcosl tanVi 2dyh l tan y 2 tan(R y)dz 0dx 0

17、R-h b) cos542论由对区域2的讨可以得3(L l)tan计算l)ta n R(1 cos cos )然后只需做下式变换即可 v ( h, a B 一般关系 域的局部横向偏转后进入了第一区域式使得 l变为l+ R Rcos tan故只需将5.3中的V1中l替换为l+ R RcosRhV2( h B =V2( R- R h )cos543对vs区域的讨论由Rh cos可以得到h=R-( R-h b) cos B在第三区域又要满足(L l) tan h 2R l tan计算得(L l) tan R (1 cos R (1 cos ) l tancos cos然后只需做下式变换即可得到V (

18、 h , a B 般关系式RhV3 ( h B =V3 (R- R h)cos544对V4区域的讨论由Rh cos可以得到h=R- ( R-h b) cos B在第四区域又要满足R (1 cos ) I tanh 2Rcos然后只需做下式变换即可得到V ( h, a, B) 一般关系式RhV4 ( h B) =V4 ( R- R h ) cos对区域V5的讨论 由于此时油浮子到达油位探针的顶点，无法进一步测量油位 高度。故对V5不再进行讨论。5.5变位参数确实定和罐容表的计算用最小二乘参数估计法确定参数最小二乘参数估计法：根据 V (h, a B)的关系表达式求得几组油量高度h,计算出相邻高度

19、油量体积之差，通过与附件的实际储油量进行比拟，通过对a B进行 等间距穷举最终求得理论值与实际值的差值平方和S，当S取最小是时，a B为所求最正确值，求解公式如下min S( , ) ( V i Vi ')i1用最小二乘参数估计法得到变位参数为：a =2.51 ° 3=4.85。，角度符合实际情况。实际储油罐罐容表制定估计出参数a 3后，我们可以根据所建立的储油量V以及油位高度h以及变位 参数a? 3的一般模型得到罐容表表七 实际储油罐罐容表(纵向变位0=2.51。，横向偏转3=4.85 ° )油位高度h/m0.10.20.30.40.50.60.70.80.91.

20、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.0354.0757.1359.2160.2761.5625.6数据检验根据推导数据V(h , a, 3)关系式对第一次补充仅有后的数据进行检验，

21、得到相 油位高度下，理论储油值与实际储油值关系图：39姗M与雪辄猛油量的宴别时比5 DWM頰71图十九理论储油值与实际储油值量差异对 比图图中显示理论值与实际值根本重合，表示算的a,B较为合理六、模型评价1. 明确给出变位参数的辨识准那么，由简单到复杂建立模型，逻辑性强2. 数据参数化，是模型更具有灵活性，应变能力高。3. 运用 matlab 程序加强数据的科学性与准确性。七、模型改良1. 由于公式复杂，程序运行较慢，将模型参数进行近似简化，会减少模型运行时 间；2. 将模型公式进行更细致的分割，分开运算各局部，来减少运行的错误以及程序 运行时间。参考文献1 光亭裘哲勇数学建模高等教育 202

22、1八、附录第一局部：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:78 v(i)=1000*a*b*L*(h(i)-b)/(bA2)*sqrt(h(i)*(2*b-h(i)+asi n( (h(i)-b)/b)+pi/2 ); end>> plot(m,v,'r*')>> hold on>> xlabel(' 油位高

23、度 h/mm');>> ylabel('储油量 /L')>>问题A附件1 :实验采集数据表.xls','G2:G79');>>问题A附件1 :实验采集数据表.xls','D2:D79');>> plot(B A'b-')；legend('理论值','实际值');title('罐体变位前的V-h曲线比照')1.2、得出上述变位前理论值与测量值的曲线图中理论与测量的比值：问题A附件1 :实验采集数据表.xls'

24、;,'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)/(bA2)*sqrt(C(i)*(2*b-C(i)+asi n( (C(i)-b)/b)+pi/2)*1000;End问题A附件1：实验采集数据表.xls','G2:G79,); fori=1:78T(i)=D(i)/V(i);end1.3正常情况下小椭圆罐的罐容表代码：h=linspace(0.1,1.2,12);>> b=0.6;>> a=1.78/2;>> L=2.45;&

25、gt;> for i=1:12 v(i)=1000*a*b*L*(h(i)-b)/(bA2)*sqrt(h(i)*(2*b-h(i)+asin(h(i)-b)/b)+pi/2 ); end结果:1.4、变位后的相关公式积分结果：V1:c =ta n(4.1*pi/180);1=0.4;d=c*l;L=2.45;syms h a b x y 乙 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 =ta n(4.1*pi/180);1=0.4;d=c*l;L=2.45;syms h a b

26、x y 乙 in t(i nt(i nt(1,z,0,(h+d-y)/c),x,0,a*sqrt(1-(y/b-1)八2),y,0,h+d)*2 h=0.1470V=0.1514V3:c =ta n(4.1*pi/180);l=0.4;d=c*l;L=2.45;syms h a b x y 乙 V=a*b*L*(h/(b2)-1/b)*sqrt(2*b*h-hA2)+asi n( h/b-1)+pi/2);Va=i nt(i nt(in t(1,z,0,(h+d-y)/c),x,0,a*sqrt(1-(y/b-1)A2),y,h,h+d)*2; Vb=2*i nt (i nt (int(1,z

27、,(h-y+d)/c,L ) ,x,0,a*sqrt(1-(y/b-1)A2 ) ,y,h-(L-l)*c,h ); V3=V+Va-Vb H=1.1713 V= 3.9588V4:c =ta n(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)/(bA2)*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*s qrt(1-(y/b-1)A2),y,h-(L-l)*c,2*b);H=1.2 V4= 4.0127区域油位高度h/

28、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.1101得到V2区间的相对应点的体积值:1.5、画出变位后测量值与理论值的曲线:h=li nspace(0,0.1470,3);for i=1:3B(i)=vpa(U ntitled2(h(i),5);endB = 1.6744, 39.597, 151.37画图：h=li nspace(0,0.1470,3);B = 1

29、.6744, 39.597, 151.37H=h*1000;>> plot(H,B,'r*')>> hold on得到V3区间的相对应点的体积值：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.

30、21722, 0.24332, 0.25921, 0.26675, 0.26668, 0.25898, 0.24292, 0.21661, 0.1753V3 =1.0e+03Columns 1 through 70.2819 0.5953 0.9657 1.3719 1.7986 2.2325 2.6614Columns 8 through 103.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*

31、1000;V3=V3*1000;>> plot(S,V3,'r*')>>hold on得到 V4 区间的相对应点的体积值： t=linspace(1.1716,1.2,6);for i=1:6 D(i)=Untitled4(h(i);endD=vpa(D,5)画图：D= 3.9342, 3.9531, 3.9707, 3.9867, 4.0009, 4.0127 t=linspace(1.1716,1.2,6); T=t*1000;D=D*1000;>> plot(T,D,'r*')根据附录一画出测量值：>> xl

32、abel ' 油位高度 h/mm' ;>> ylabel ' 储油量 /L'>> 问题 A 附件 1：实验采集数据表.xls',' 倾斜位进油 ','G2:G79' ;变问题A附件1:实验采集数据表.xls','位进油 ','D2:D79' ;>> plotZ,Y,'b-'legend ' 理论值',' 实际值 ' title' 罐体变位后的 V-h 曲线比照 '1.6 、变位前后储油

33、罐理论值的比拟:问题 A 附件 1:实验采集数据表 .xls','G2:G79' 问 题 A 附件 1 :实验采集数据表 .xls','D2:D79'plotB,A,'b-'>> hold on>> h=linspace0,0.1470,3;B = 1.6744, 39.597, 151.37H=h*1000;plot(H,B,'r*')B =1.6744 39.5970 151.3700 >> hold on>> s= lin space(0.2,1.1,10);V

34、3=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.0127 t=li nspace(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

35、4.0127 >> hold on>> legend('变位前,'变位后');title('罐体变位前后理论值的V-h曲线比照')xlabel('油位高度h/mm');ylabel('储油量 /L')1.7、变位前后理论值随高度变化的差值由上述罐体变位前油位高度间隔为10cm的罐容表标定值表以及罐体变位后油位高度间隔为10cm的罐容表标定值表，可得相应的差值为：变位前后理论值随高度变化的差值表高度h/cm变位后与变位前储油量差值/L1093.4920168.3730208.3340233.65024

36、9.160256.570256.680249.490234.2100209.2110169.91200高度h/cm变位后与变位前储油量相对误差100.5714897200.373931197300.259236216400.194780289500.153670574600.124811445700.10308947800.085680912900.0708280411000.0571600321100.0430497141200第二局部代码：2.1根据最小二乘法以及穷举法，算得a=2.51。，3=4.85以下为相应的计算不同区域的体积的代码： syms x y z n I m h;r=sqr

37、t(1.5A2-(y-1.5)A2);h =1.5-(1.5-h)/n;I =(1.5-1.5* n)/m+l;V 球缺：syms h y RO R H r =sqrt(R0A2-(h-R)A2);in t(rA2*acos(R0-H)/r) -sqrt(rA2-(R0-H)A2)*(R0-H),h,0,h) ans = in t(H - R0)*(R0A2 - (H - R0)A2 -(R - h)A2)A(1/2) + (R0A2 - (R - h)A2)*(pi - acos(H - R0)/(R0A2 - (R -小八2)八(1/2), h, 0, h)>> R0=1.62

38、5;H=1;R=1.5;>> in t(H - R0)*(R0A2 - (H - R0)A2 - (R -小八2)八(1/2) + 侶0八2 - (R - h)A2)*(pi-acos(H - R0)/(R0A2 - (R -小八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)* - 3/2)八2 - 169/64),h, 0, h)>> h=0:0.004:2.9for i=1:length(h)z(i)= - (5

39、*(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.5A2-(y-1.5)A2); Vc=i nt(in t(i nt(1,z,0,(h+l*m-y)/m),x,0,r),y,0,h+l*m)*2;Vb=i nt(-acos(5/(8*(169/64 - (y - 3/2)八2)八(1/2)*© - 3/2)八2 - 169/64)-(5*(9/4 -

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

41、- acos(5/(8*(169/64 - (y - 3/2)A2)A(1/2)*(y - 3/2)A2 - 169/64) -(5*(9/4 - (y - 3/2)A2)A(1/2)/8, y, 0, h)*(h + l*m) + 2*int(l + (h - y)/m)*(9/4 - (y - 3/2)A2)A(1/2), y, 0, h + l*m)ans =2*int(9/4 - (y - 3/2)A2)A(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

42、 - 3/2)/n + 3/2) + int(- acos(5/(8*(169/64- (y - 3/2)A2)A(1/2)*(y - 3/2)A2 - 169/64) - (5*(9/4 - (y - 3/2)A2)A(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)A2)A(1/2)*(l + (h - 3/2)/n - y + 3/2)/m - (3*n)/2 - 3/2)/m), y, 0, m*(

43、l - (3*n)/2- 3/2)/m) + (h - 3/2)/n + 3/2) + int(- acos(5/(8*(169/64- (y - 3/2)A2)A(1/2)*(y - 3/2)A2 - 169/64) - (5*(9/4 - (y - 3/2)A2)A(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)A2)A(1/2)*(l + (h - 3/2)/n - y + 3/2)/m - (3*n)/2 - 3/2)/m)

44、, 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)A2)A(1/2)*(y - 3/2)A2 - 169/64) - (5*(9/4 - (y - 3/2)A2)A(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 - 3/2)A2)A(1/2)*(l + (h - 3/2)/n - y

45、+ 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)A2)A(1/2)*(y - 3/2)A2 - 169/64) - (5*(9/4 - (y - 3/2)A2)A(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)A2)A(1/2)*(h - 3/

46、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)A2)A(1/2)*(y - 3/2)A2 - 169/64) - (5*(9/4 - (y - 3/2)A2)A(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 I m h; r=sqrt(1.5八2-(

47、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+I*m)/m,8),x,0,r),y,h-(8-I)*m,h);>> VQ1=i nt(-acos(5/(8*(169/64 - (y - 3/2)A2)A(1/2)*(y - 3/2)八2 - 169/64) - (5*(9/4 - (y - 3/2)A2)A(1/2)/8, y, 0, h)*(h+l*m);>> VQ2=int(-acos(

48、5/(8*(169/64 - (y - 3/2)A2)A(1/2)*(y - 3/2)A2 - 169/64) - (5*(9/4 - (y - 3/2)A2)A(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)A2)A(1/2)*(h - 3/2) + int(- acos(5/(8*(169/64 - (y - 3/2)A2)A(1/2)*(y - 3/2)A2 - 169/64) - (5*(9/4 - (y - 3

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

50、m)h =1.5-(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)A2)A(1/2)*(h - 3/2) + int(- acos(5/(8*(169/64 - (y- 3/2)A2)A(1/2)*(y - 3/2)A2 - 169/64) - (5*(9/4 - (y - 3/2)A2)A(1/2)/8, y, 0, h)*(h + l*m) - 2*int(-(9/4 - (y - 3/2)A2)A(1/2)*(l + (h - y)/m - 8), y, h + m*

51、(l - 8), h) + int(- acos(5/(8*(169/64 - (y- 3/2)A2)A(1/2)*(y - 3/2)A2 - 169/64) - (5*(9/4 - (y - 3/2)A2)A(1/2)/8, y, 0, h)*(h + m*(l - 8) + 2*int(l + (h - y)/m)*(9/4 - (y - 3/2)A2)A(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)A2)A(1/2)*(l + (h- 3/2)/n - y + 3/2)

52、/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)A2)A(1/2)*(y - 3/2)A2 - 169/64) - (5*(9/4 -(y - 3/2)A2)A(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)A2)A(1/2)*(y -

53、 3/2)A2 - 169/64) - (5*(9/4 - (y - 3/2)A2)A(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)A2)A(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

54、/2)A2/nA2)A(1/2)*(h - 3/2)/nsyms x y z n l m h;9*pi + 18*asin(2*(h - 3/2)/(3*n) + 2*int(9/4 - (y - 3/2)A2)A(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)A2)A(1/2)*(y - 3/2)A2- 169/64) - (5*(9/4 - (y -

55、3/2)A2)A(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)A2)A(1/2)*(y - 3/2)A2 - 169/64) - (5*(9/4 - (y- 3/2)A2)A(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*in t(-(9/4 - (y - 3/2)八2)八(1/2)*(1 + (h - 3/2)/n - y + 3/2)/m- (3*n)/2 - 3/2)/m - 8), y, (h