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

1、储油罐的变位识别与罐容表标定摘要 本文对储油罐的变位识别和罐容表的标定问题进行了深入探讨，建立了储油量和油位高度以及变位参数之间的数学模型，并利用matlab软件对数据进行处理。 对于问题一，我们利用积分法分别求出小椭圆储油罐在无变位和发生纵向倾斜变位之间时的一般公式。并求出储油罐纵向变位对罐容表的影响。分别给出了小椭圆储油罐在无变位和纵向倾斜变位的罐容表。如下：表一 罐体变位前油位高度间隔为10cm的罐容表标定值表油量高度/cm油量/L油量高度/cm油量/L油量高度/cm油量/L10163.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的小椭圆型储油罐（两端平头的椭圆柱体），分别对罐体无变位和倾斜角为a=4.10的纵向变位两种情况做了实验，实验数据如附件1所示。请建立数学模型研究罐体变位后对罐容表的影响，并给出罐体变位后油位高度间隔为1cm的罐容表标定值。（2）对于图1所示的实际储油罐，试建立罐体变位后标定罐容表的数学模型，即罐内储油量与油位高度及变位参数（纵向倾斜角度a和横向偏转角度b ）之间的一般关系。请利用罐体变位后在进/出油过程中的实际检测数据（附件2），根据你们所建立的数学模型确定变位参数，并

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

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

7、为倾斜后，油位过高或过低都无法起作用，如下图琐事v1和v5区域，讲储油罐按液面高低分为五个部分，来球储油量和油位高度之间的关系。如下图二所示： 图2 储油罐分区示意图4.2.1 对区域v1的讨论在区域v1，其油位低于油位探针的油浮子，故油位显示高度始终为0。将区域v1放大至图3：图3 区域v1放大图如图中，以原点（油罐正视图左侧顶点）纸面向里为x轴，利用三重积分得到其中l为油探针到储油罐左侧的距离4.2.2 对区域v2的讨论将区域V2放大至图4，如下图所示：图4 区域v2示意图 由于区域v1和区域v2的相似性，可以得到区域v2的储油量和油位高度的变化关系，得出4.2.3 对区域v3的讨论将区域

8、V3放大至图5，如下图所示： 图5 区域v3示意图将区域v3分割为两部分，一部分为va另一部分为v，其中v又包含了vb，利用积分远离我们可一分别计算出va和vb，并可以得出油量为v3=va+v-vb其中4.2.4 对区域v4的讨论将区域V4放大至图6，如下图所示：图6 区域v4示意图 利用总体积减去油罐内无油部分的体积vb即可求出v4的体积其中4.2.5 对区域v5的讨论 如下图所示，由于此时油浮子到达油位探针的顶点，无法进一步测量油位高度。故对v5不再进行讨论。图7 区域v5示意图4.2.6 综合各区域的罐容表标定的数学模型根据对各区域的讨论我们可以得到罐体纵向倾斜变位（=4.1。）后油位高

9、度间隔为1cm的罐容表标定值（见附录一）。 当=4.1。时，各区域油位高度及体积变化范围为如下表所示：表四 各区域油位高度及体积变化范围区域油位高度 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 罐体变位后对罐容表的影响为了能更准确刻画罐体的纵向倾斜变位对罐容表的影响，我们分别对罐体变为前后的理论值和测量值进行多方面的比较，所得结果分别如下图所示：

10、 图九4.3.1 罐体变位前理论值与测量值的比较根据附录一中所给数据，我们计算出在附件所给的油位高度下的理论值和实际值，并划出其曲线，如图十所示：图十同过对比我们发现对于任意h，储油量理论值和实际值始终成比例，比例保持为0.9663.（具体代码见附录）4.3.2 罐体变位后理论值与实际值比较经过MATLAB软件的作图，图片如下所示：图十一如图所示测量值始终小于理论值，但保持在极小的误差下。4.3.3 罐体纵向倾斜变位前后理论值的比较表五：变位前后理论值随高度变化的差值表高度 h/cm变位后与变位前储油量差值/L1093.4920168.3730208.3340233.650249.160256

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

12、示意图，其主体为圆柱体，两端为球冠体，在储油罐无变位时我们计算其各部分的体积图十三 圆柱体积计算公式一端球缺体积计算公式为5.1 考虑不发生变位时储油量和油位高度的关系圆柱内的油体积随油位高度变化关系：一端球缺内油体积随油位高度变化关系：R为圆柱体地面半径R0为球缺对应半径r圆为球内小圆半径总储油量与油位高度变化关系为5.2 只考虑横向偏转变位是的储油量和油位高度关系 只考虑储油罐横向偏转为时：对实际的油位高度没有影响，但此时的油位探针随储油罐发生偏转图十四 只考虑横向偏转示意图由油浮子测量得到的油位高度与实际油位高度的关系为所以对罐容表影响转换公式：5.3 只考虑纵向倾斜变位时的储油量和油位

13、高度的关系 利用问题一中的思想，将储油量分成5个区域，分别求储油量和油位高度的关系，如图10图十五 实际储油罐分区示意图5.3.1 对区域1的讨论油位探针测得的油位高度h始终未0，将其分成三部分来计算这部分体积，如区域放大图11图十六 区域1放大图区域1的总体积应为：对圆柱体三重积分得到：球缺部分计算公式得到由于Va部分体积对于整体来说较小，故可忽略5.3.2 对区域2讨论图十七 区域2示意图 在区域2种，油位高度h变化范围为0，（L-l）tan各部分储油体积和油位高度的变化关系：5.3.3 对区域3的讨论图十八 区域3示意图 如图，在区域3中，油位高度h的变化范围（L-l）tan，2R-lt

14、an储油量为利用问题一的类似解法5.3.4 对区域4的讨论图十七 区域4示意图 在区域4中，油位高度h的范围2R-ltan,2R.储油量与油位高度的关系为其中5.3.5 对区域v5的讨论 如下图所示，由于此时油浮子到达油位探针的顶点，无法进一步测量油位高度。故对v5不再进行讨论。图十八 区域5示意图5.4 综合考虑储油罐纵向倾斜和横向偏转根据如上讨论，把5.3中各区域V-h公式中的h替换为即可得到储油罐纵向倾斜和横向偏转的V（h，）一般关系式。依然考虑将储油罐分成五个区域进行求解5.4.1 对区域1的讨论当发生横向偏转后，油位探针测得油位高度为0，h=R-Rcos，故本属于第二区域的部分横向偏

15、转后进入了第一区域，使得l变为l+故只需将5.3中的V1中l替换为 l+故5.4.2 对区域2的讨论由可以得到h=R-（R-h）cos在第二区域又要满足计算得0然后只需做下式变换即可得到V（h，）一般关系式V2（h）=V2（R-）5.4.3 对v3区域的讨论由可以得到h=R-（R-h）cos在第三区域又要满足计算得然后只需做下式变换即可得到V（h，）一般关系式V3（h）=V3（R-）5.4.4 对v4区域的讨论由可以得到h=R-（R-h）cos在第四区域又要满足2R然后只需做下式变换即可得到V（h，）一般关系式V4（h）=V4（R-）5.4.5 对区域v5的讨论 由于此时油浮子到达油位探针的顶

16、点，无法进一步测量油位高度。故对v5不再进行讨论。5.5 变位参数的确定和罐容表的计算5.5.1 用最小二乘参数估计法确定参数 最小二乘参数估计法：根据V（h，）的关系表达式求得几组油量高度h，计算出相邻高度油量体积之差，通过与附件的实际储油量进行比较，通过对，进行等间距穷举最终求得理论值与实际值的差值平方和S，当S取最小是时，为所求最佳值，求解公式如下用最小二乘参数估计法得到变位参数为：=2.51°，=4.85°，角度符合实际情况。5.5.2 实际储油罐罐容表制定估计出参数，后，我们可以根据所建立的储油量V以及油位高度h以及变位参数，的一般模型得到罐容表表七 实际储油罐罐

17、容表（纵向变位=2.51°，横向偏转=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数据检验根据推导数据V（h，）关系式对第一次补充仅有后的数据进行检验，得到相同油位高度下，理论储油值与实际储油值关系图： 图十九 理论储油值与实际储油值量差别对比图图中显示理论值与实际值基本重合，表示算的，较为合理。六、模型评价1. 明确给出变位参数的辨识准则，由简单到复杂建立模型，逻辑性强2. 数据参数化，是模型更具有灵活性，应变能力高。3. 运用matlab程序加强数据的科学性与准确性。七、模型改进1. 由于公式复杂，程序运行较慢，将模型参数进行近似简化，会减少模型运行时间；2. 将模型公式进行更细致的分割，分开运算各部分，来减少

19、运行的错误以及程序运行时间。参考文献 1陈光亭 裘哲勇 数学建模高等教育出版社 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：实验采集数据表.xls','G2:G79');>> B=xlsread('C:UsershfDesktop问题A附件1：实验采集数据表.xls','D2:D79&

21、#39;);>> plot(B,A,'b-');legend('理论值','实际值');title('罐体变位前的 V-h曲线对比')1.2、得出上述变位前理论值与测量值的曲线图中理论与测量的比值：A=xlsread('C:UsershfDesktop问题A附件1：实验采集数据表.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：实验采集数据表.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结果：1.4、变位后的相关公式积分结果：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区域油位高度 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区间的相对应点的体积值：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区间的相对应点的体积值：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区间的相对应点的体积值：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*')根据附录一画出测量值：>> xlabel('油位高度 h/mm');>> 

30、;ylabel('储油量/L')>> Y=xlsread('C:UsershfDesktop问题A附件1：实验采集数据表.xls','倾斜变位进油','G2:G79');>> Z=xlsread('C:UsershfDesktop问题A附件1：实验采集数据表.xls','倾斜变位进油','D2:D79');>> plot(Z,Y,'b-');legend('理论值','实际值');

31、title('罐体变位后的 V-h曲线对比')1.6、变位前后储油罐理论值的比较：A=xlsread('C:UsershfDesktop问题A附件1：实验采集数据表.xls','G2:G79');B=xlsread('C:UsershfDesktop问题A附件1：实验采集数据表.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的罐容表标定值表以及罐

34、体变位后油位高度间隔为10cm的罐容表标定值表，可得相应的差值为：变位前后理论值随高度变化的差值表 高度 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 根据最小二乘法以及穷举法，算得=2.51。，=4.85。以下为相应的计算不同区域的体积的代码：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)/