参考格式-全国二等奖论文A题

更新时间：2024-05-30 09:28:01 阅读量：综合文库文档下载

说明：文章内容仅供预览，部分内容可能不全。下载后的文档，内容与下面显示的完全一致。下载之前请确认下面内容是否您想要的，是否完整无缺。

论文参考格式参考文献格式推荐度：
相关推荐

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

张三李四王五

摘要：本文针对油罐的变位识别和罐容表标定的具体问题，通过建立模型

将储油分成无数个小微元，研究小微元的规律从而利用微积分的方法推导出，当纵向倾斜?角度和横向偏移?角度一定时，通过建立数学模型将储油罐中剩余油的体积与油位探针显示高度的函数关系表示出来，并制定罐容表标定值。

针对问题一：研究罐体变位对灌容表的影响，首先利用了逼近法计算油罐未变位前的储油体积，在MATLAB中对实际与理论数据进行拟合利用误差分析公式

E?理论值?真实值真实值，求得其误差E范围为0.014866?E?0.034917，求储油罐

中储油量时将油罐分成多部分考虑，利用微元思想和积分方法求得其储油量的体积与油位探针的读数h及变位角?，?的函数关系V?f?h,?,??在此问中

????0,??4.1时得出V?f?h,4.1,0?，求出其一定h时的V。用模型求出的理论

值与题目附表1中的实际值相比较，得出其误差1.42%?E?3.37%。并标出变位后间距为1cm罐容表。

针对问题二：对实际储油罐建立罐体变位罐容表的标定模型。在问题一的理论和方法的基础上加上球冠中的储油体积即可得到实际储油罐储油体积，并采用最小二乘法推导出所求变位参数?及?。并得出当?及?一定时油位探针的高度与剩余储油体积的关系V?f?h?，进而制定间距为10cm罐容表标定值。本文充分运用了数学分析、高等数学等知识对储油体积积分，并通过MATLAB软件模拟的方法对理论数据进行了误差分析，以及运用最小二乘法估计其变位参数值?，?。最后对模型的优缺点进行了评价，并给出了改进方向。

关键词：MATLAB数据拟合、微元法、最小二乘法、罐容表标定

1 问题重述

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

通常加油站都有若干个储存燃油的地下储油罐，并且一般都有与之配套的“油位计量管理系统”，采用流量计和油位计来测量进/出油量与罐内油位高度等数据，通过预先标定的罐容表（即罐内油位高度与储油量的对应关系）进行实时计算，以得到罐内油位高度和储油量的变化情况。

许多储油罐在使用一段时间后，由于地基变形等原因，使罐体的位置会发生纵向倾斜和横向偏转等变化（以下称为变位），从而导致罐容表发生改变。按照有关规定，需要定期对罐容表进行重新标定。图1是一种典型的储油罐尺寸及形状示图3是罐体横向偏转变位的截面示意图。

请你们用数学建模方法研究解决储油罐的变位识别与罐容表标定的问题。

（1）为了掌握罐体变位后对罐容表的影响，利用如图4的小椭圆型储油罐（两端平头的椭圆柱体），分别对罐体无变位和倾斜角为??4.1°的纵向变位两种情况做了实验，实验数据如附件1所示。请建立数学模型研究罐体变位后对罐容表的影响，并给出罐体变位后油位高度间隔为1cm的罐容表标定值。

（2）对于图1所示的实际储油罐，试建立罐体变位后标定罐容表的数学模型，即罐内储油量与油位高度及变位参数（纵向倾斜角度?和横向偏转角度? ）之间的一般关系。请利用罐体变位后在进/出油过程中的实际检测数据（附件2），根据你们所建立的数学模型确定变位参数，并给出罐体变位后油位高度间隔为10cm的罐容表标定值。进一步利用附件2中的实际检测数据来分析检验你们模型的正确性与方法的可靠性。

2 模型假设

[1] 假设油浮子不会超过油位探针颈口，在要到油位探针劲口时就停止注油； [2] 假设温度和压强对油的体积不影响；

[3] 假设油浮子，油位探针和出入油导管的体积不计；

[4] 假设问题一中的储油罐为标准的椭圆柱体，没有凹凸地方；

[5]假设储油罐只沿着题目中所示的纵向倾斜，当其沿反方向倾斜时其处理方法

类似，在此题中我们不予考虑； [6] 假设纵向倾斜变位的角度??9.5?。

3 符号说明

? 油罐纵向倾斜的角度

? 油罐横向偏转的角度

油位探针的显示高度 L 油罐的长度 V 储油的体积 a 椭圆的长半轴 b 椭圆的短半轴

h4 问题的分析

本文主要研究油罐变位对罐容和罐容表的影响，首先我们利用微元思想与体积积分的算法建立油面高度与体积关系的模型。然后对问题一利用附件1的油位高度所对应的实际油量的体积与模型求出的理论油量的体积进行数据拟合。比较拟合曲线得到相应的误差范围。

在问题一里面由于椭圆柱的纵向偏离4.1°,罐容和罐容表都将受到影响。我们将倾斜的椭圆柱体分成三部分，并对每一部分通过微元和积分的思想和方法对其进行计算。从而推导纵向变位（倾斜）对罐容的影响，并标定出以1cm为间隔的罐容表。

在问题二里面我们利用第一题的思想方法，纵向倾斜为例，将球冠和圆柱都沿其轴线切截，从而油面形成弓形，利用微积分，分别对两边的球冠和倾斜的圆柱积分，从而得出储油的体积。

5模型的建立与求解

本节主要研究变位对椭圆柱油罐罐容表影响，分别对它在变位前后油罐中油体积的计算，从而推导倾斜对罐容表的影响，标定间隔为1cm的罐容表。

由于油罐是不规则容器，因此油罐中储油量的体积可用微元法的思想来求解。

微元法求几何体体积的基本原理

如图1所示的一个空间立体物体，设它在x 处截面面积为S(x)，利用极坐标下推导面积公式的思想求出它的体积？

如果纵向切把它切成薄片，则每个薄片可近似看作直柱体，其体积等于底面积乘高，所有薄片体积加在一起就近似等于该立体的体积。

在[a, b]上任取一点x，并且任给x的一个增量?x，这样就得到一个非常薄的薄片，这个小薄片我们可以近似地把它看成柱体，于是这个微小的柱体体积为：

S(x) a ⊿x b dv?s?x??x?A?x?dx 图1

把这些小体积加起来，就是我们要求的体积即V?

?S?x?dx

ab5.0油罐中理论油量与实际油量的数据拟合

如图2所示建立坐标系，设截面椭圆的方程为

Xa22?Yb22?1,椭圆弓形的高

为h，图2（储油截面图）中带阴影部分为储油截面，用定积分求椭圆的体积。设弓形的面积为S(h)则

S?h??2ab?h?b?b??h?b22b?ydy???b?2??h?b?1???b??2?arcsinh?bb??ab?? (1)

设x?h?bb?由0?h?2b,可知?1?x?1?油罐的长为L，储油的体积为V?x?，可得

???2V?x????x1?x?arcsinx?abL?2? (2)

图2 储油截面图

利用公式（2）所求的结果和附件1的油面高度得到对应的储油量体积，将这些点拟合得到理论的储油量变化曲线，并与题目给出的实际的储油量的拟合曲线比较，如图3所示。

450040003500300025002000150010005000-500 0200400600800100012001400理论数据理论拟合曲线实验数据实验拟合曲线

图3 无变位时的数据拟合图

由图像可知理论值大于实际数据，利用相对误差E?理论值?实际值实际值表示

它们的误差，计算出误差E的范围为：0.014866?E?0.034917。产生误差的原因可能是：

(1)油浮子，探针，以及油的出入管道的体积

(2)由于油的出入导致罐内压强改变，从而导致理论数据与实际数据的误差 (3)进出油环境温度不一致

5.1模型的建立与求解

如图以椭圆柱底线为Y轴，过椭圆心，垂直于Y轴的直线为Z轴，过椭圆左端面的最低点垂直于YOZ面为X轴，其空间坐标的圆心为椭圆左端面的最低点。当椭圆柱形油罐纵向倾斜时，我们将油罐的体积分为三个阶段。

图4 倾斜的椭圆柱

1.求每一截面的面积建立的坐标系如图4，设h1为椭圆柱端截面与油面的交线与Z轴的交点为A，则该点的坐标为A（0，h1）。油面与椭圆柱面的交线与Y轴的交点为B，则该点的坐标为B（

Z?h?y?tan?1h1tan?,0）。由A，B两点可得油面的方程为：

我们沿平行于XOZ面切截椭圆柱，则所切截面均为椭圆。

沿着平行于XOZ面切截椭圆柱，其横截面必为椭圆，其椭圆的表达式为：

xa22??z?b?2b2?1

油面与z轴的交点即油面的高为h1，由图5和题意可知切截椭圆面中的油面所对应的Z轴坐标即为椭圆中储油的高，我们设其对应的高为h11则储油部分的

面积为：S?2?h110?z?b?1???dz?b?2，由于Z?h1?ytan?，所以 h11?h1?ytan?

图5 变位椭圆柱的切截图

2.求每一对应液面高度的体积

2.1当油面与斜放的椭圆柱的下柱面相交时即：h1?Ltan?

此时我们对每一个切截面的油面积从y1?0到y2?果即为储油量的体积。

V1?h1tan?进行积分，其积分结

?y2y12ab?h110b??z?b?dzdy?22?y2y12ab?h?ytan?10b??z?b?dzdy22

图 6油面上升立体图

2．2当油面交于椭圆柱两端面时

即Ltan??h1?1.2此时油面方程不变仍为：

z?h?ytan?1,h11?h1?ytan?

同理对其切截面及积分可得其体积，但对油面积分的上限变为了L，即：

V2??L2ab0?h?ytan?10b??z?b?dzdy22

图7油面淹没左端面最高点示意图

2.3当油面对斜椭圆柱的上柱面和端面相交时如图7

即2b?h1?2b?Ltan?，由油面的方程可得：设其右端面油面的高为h1则

h1?h?Ltan?1

设h2为油面与右端面交线到椭圆柱顶部的高：

h2?2b?h1?2b?Ltan??h1

由图7可知上下两阴影部分的体积相等。

所以油的体积V3=椭圆柱体总体积-对右上角没油部分对称的左下部分的体积。

对其右上没油部分的体积由对称性可将h2代入式子V1中计算。又因为椭圆

h2柱体的体积为?abL，所以V3??abL?代入h2?2b?Ltan??h1得

2b?Ltan??h1?tan?2ab0?h2?tan?0b??z?b?dzdy22

V3??abL??tan?2ab0?2b?Ltan??h?ytan?10b??z?b?dzdy

22设其油位探针的显示高度为h，则由几何关系可得h1?h?0.4tan?。将此关系带入V?f(h1)可得变位后探针高度与储油体积的关系如下 (1)当h?Ltan??0.4tan?时

V1?y2y1?2ab?h??y?0.4?tan?0b??z?b?dzdy22

（2）当?L?0.4?tan??h1?1.2?0.4tan?时

V2??L2ab0?h??y?0.4?tan?0b??z?b?dzdy22

（3）当2b?0.4tan??h1?2b??L?0.4?tan?时

2b??0.4?L?tan??hV3??abL??tan?2ab0?2b??0.4?L?tan??h?ytan?0b??z?b?dzdy

225.2椭圆柱体的储油罐的罐容表

在MATLAB的环境下计算油浮子的所显示的高度每改变1cm时的储油体积（即变位后的储油罐的罐容表）。其罐容表见表1。纵向变位时的数据拟合结果如下。

5000实实理理际际论论数数数数据据拟合曲线据据拟合曲线 40003000200010000-1000 0200400600800100012001400

图8纵向变位时的数据拟合曲线表1 椭圆柱纵向变位后的罐容表油高(mm) 油量(L) 油高(mm) 油量(L) 油高(mm) 油量(L) 油高(mm) 油量(L) 0 10 20 30 40 50 60 70 80 90 100 110 120 130 140 150 160 170 180

<=1.67 3.53 5.94 9.97 14.76 21.04 27.85 36.32 47.56 63.16 78.09 84.4 100.25 117.4 136.9 157.8 180.3 204 228.9 310 320 330 340 350 360 370 380 390 400 410 420 430 440 450 460 470 480 490 630.1 665.6 701.5 738 774.9 812.2 850 888.2 926.7 965.7 1005 1044.6 1084.5 1124.8 1165.3 1206.2 1247.2 1288.6 1330.1 8

610 620 630 640 650 660 670 680 690 700 710 720 730 740 750 760 770 780 790 1841.8 1885.1 1928.5 1971.9 2015.4 2058.8 2102.3 2145.7 2189.1 2232.5 2275.8 2319.1 2362.3 2405.4 2448.4 2491.3 2534 2576.6 2619.1 910 920 930 940 950 960 970 980 990 1000 1010 1020 1030 1040 1050 1060 1070 1080 1090 3112 3151.2 3190.1 3228.6 3266.7 3304.4 3341.7 3378.5 3414.9 3450.7 3486.1 3520.9 3555.1 3588.8 3621.8 3654.2 3685.9 3716.9 3747.2 190 200 210 220 230 240 250 260 270 280 290 300 254.9 281.9 309.8 338.5 368.1 398.5 429.7 461.5 494 527.1 560.9 595.2 500 510 520 530 540 550 560 570 580 590 600 1371.9 1413.9 1456 1498.4 1540.9 1583.5 1626.3 1669.2 1712.2 1755.3 1798.5 800 810 820 830 840 850 860 870 880 890 900 2661.4 2703.6 2745.5 2787.2 2828.7 2870 2911.1 2951.8 2992.3 3032.5 3072.4 1100 1110 1120 1130 1140 1150 1160 1170 1180 1190 1200 3776.6 3805.3 3833 3859.8 3885.6 3910.3 3933.9 3956.1 3975.3 3995.5 4016.7 二、问题二的模型的建立与求解：

由题目可得，罐体变位后对标定罐容表的影响主要是罐体纵向倾斜和横向偏转。而要想根据罐内储油量与油位高度来建立数学模型，我们分为两种情况来讨论。

1只考虑纵向倾斜

由于在实际储油罐变位问题中我们利用了问题一的理论基础将储油罐分为5个部分，分别以储油罐左端面的最高点，圆柱右端面的最低点，两端球冠的顶点为边界，共4个分界点，如图9所示。由于储油罐在安放时是水平卧放的，经过土质的不均匀的压缩等原因导致出现纵向倾斜角度?，并且其?不会过大。所以我们假设??9.5?，根据此角度，通过计算可以定出油面上升时淹没4个分界点的顺序，其分段区域如图9中虚线所示：

图9分段区域示意图

在求储油罐的体积时，首先要求出两端球冠体的圆心坐标，为此，我们建立了球的俯视图，如图10所示，利用勾股定理求出该球冠的半径。我们仍沿着平行于XOZ面切截球冠，从而得到该球冠的每个切面的表达式与该球半径的关系

R1?R??0.625?y?

22假设球冠所在球的半径为R，则由图10可得：R2??R?1??1.52。代入数据，

2计算可得球冠所在球的半径R=1.625m。

图10俯视图

因为油的高度可以在图9的分段区域的任何位置，所以求储油罐内油的体积时就要分为5中情况考虑，分别如下：

图11切割后截面油面积示意图

(1)当油面即没有超过左端球冠的左顶点又没有超过右端面的下底点时，所截的切面如图11所示，此时h1?8tan?。我们在求此部分的体积时，利用了问题一的原理，把图从坐标z轴分为两部分，一部分是由球冠体的一部分组成，另一部分是由所截的圆柱体的一部分组成。设球冠体中油量的体积V11，圆柱体中油量的体积为V1。

在求图11的体积时，我们采用垂直于y轴并平行于z轴的切面去切图，用三重积分的方法来求其体积。由于截面截球冠部分投影到z轴的阴影部分为圆，

R1为所以我们用圆的面积作为变量来求球冠的体积。假设所截的圆的半径为R1，

2变量。由勾股定理，得R2??R?1?y??R1，其中R=1.625m.即

22R1?1R??0.625?y?，由于截面是水平的，且倾斜角度为?，所以截面方程

22z?h?ytan?。每个切面的油面高度h11?h1?ytan?，所以切面的面积：

S?2?h111.5?R12R1??z?1.5?dz。

2假设油面中轴线与左球冠相交点为A，其相对于Ｙ轴的坐标为y0，则由三重积分的原理得到球冠的体积为：

V110h?ytan?1??y0.2?1.5?R1R??0.625?y???z?1.5?dzdy

222由于实际储油罐的主体为圆柱体，而问题一中的主体的椭球体，所以我们采用和问题一同样的方法求圆柱体的体积，表达式为：

h1V1?2?tan?02?h?ytan?011.5??z?1.5?dzdy22

即储油罐内油的总体积为：V1?V1?V112

图12第二种情况的示意图

(2)当油面在低于倾斜球冠的左顶点且高于圆柱右端面的最低点时所截的图如图12阴影部分所示。此时8tan??h1?1.5?tan?，储油罐内油量的体积可以分为三个部分讨论，分别为：左边球冠储油量体积V21,圆柱体部分的储油量体积

V22，右边球冠的储油量体积V23。

因为截面没有过左端球冠的左顶点，投影到z轴的切面也是圆的一部分，所

11以我们采用和第(1)种求球冠体积相同的方法，其表达式和V11相同，即： V2=V1。

在求圆柱体的体积时，其分析和问题一一样，只是这里被积函数不一样而已，因为这里的实际球罐是圆柱体，而第一问是椭球体，利用三重积分建立的表达式为：

V2?2?L02?h?tan?011.5??z?1.5?dzdy22

在求右端球冠的储油量体积时，其方法和左端一样，也是采用垂直于y轴并平行于z轴的切面去切球体，同理，利用三重积分可以建立体积的表达式，如下所示：

V2?3?y182?h?ytan?11.5?R2R??y?7.375???z?1.5?dzdy

222即储油罐内的油量的总体积为：V2?V21?V22?V23

图13第三种情况的示意图

(3)当油面在高于左球冠顶点低于右球冠顶点时所截的图如图13所示，此时：

1.5?tan??h?1.5?9tan?1。计算储油量的体积时和第二种情况一样，也是分成

23V3，V3。三个部分讨论，即球冠体的左右边和圆柱体中储油量的体积分别为V31，

在建立体积的表达式时，此时由于左端油面超出了球冠的左端顶点，在计算

这部分的面积时要分为两部分讨论，因为左端顶点到截点这段距离投影到z轴的阴影部分为圆的整体，所以这部分的面积的表达式和前面的不一样，同样用三重积分的方法，建立表达式为：

V3?1?y0?1?R1dy?2?0y02?h?ytan?11.5?R1R1??z?1.5?dzdy22

求圆柱体的体积和前面的方法一样，建立的表达式为：

V3?2??L02?h?ytan?011.5??z?1.5?dzdy22

同理，由于右端油面没有过右端最高顶点，所以其表达式和前面一样，即：

V3?3y182?h?ytan?11.5?R2R??y?7.35???z?1.5?dzdy222

所以储油罐内的油量的总体积为：V3?V31?V32?V33.

图14第5种情况的示意图

(4)当油面高于右球冠的顶点且低于圆柱左端面最高点时所截的图如图14所示，此时：1.5?9tan??h1?3.求储油量的体积仍采用分为三部分讨论。求左端球冠的体积时，其方法和第(3)问的一样，其表达式为：

V4?1?y0?1?R1dy?2??y00h?ytan?11.5?R1R1??z?1.5?dzdy?V3221

同理，圆柱体的体积也和前面的一样，即：V42＝V32。因为右端油面过了右端顶点，求其体积的时候要和左端求体积的方法一样，也是分为两部分讨论，只是这时右端和左端的球冠的方程不一样，被积函数也不一样，利用三重积分，建立的表达式如下：

V4?3?y182?h?ytan?11.5?R2R??y?7.35???z?1.5?dzdy?222?9y1?R2dy2

所以储油罐内的油量的总体积为：V4?V41?V42?V43

(5)当油面高于圆柱左端面低于油位探针的颈口时，此时：3?h1?3?2tan? 储油量的体积我们分为四部分讨论。首先，由于此时左端油面已经高过左端口，其体积即为球冠的体积，利用三重积分建立的表达式为：

V5?1?0?1?R1dy?2?0?122??R??0.625?y??dy

??由于油面过了左端球冠最高点，还没有到油位探针的位置，所以在求圆柱体的体积时，分为两部分讨论，一部分为左端最高点到油面的最高点，这部分的体积就为圆柱体的体积了，即：

V5???1.5?y1

22另一部分为圆柱体从油面最高点到右端球冠部分的圆柱体的体积，这部分和前面的一样，建立的表达式为：

V5?3?Ly12?h?ytan?011.5??z?1.5?dzdy22

最后右端球冠的体积和前面的方法一样，建立的表达式为：

V5?4?y182?h?ytan?11.5?R2R??y?7.35???z?1.5?dzdy?222?9y122??R??y?7.35??dy??所以球罐内的油量的总体积为：V5?V51?V52?V53?V54

这里我们对建立的模型取了不同的高度和倾斜角度来检验模型的准确。其中纵向倾斜角度??3?时高度间隔10cm时的罐容表标定值见如表2所示。

表2.储油罐为主体圆柱两端为球冠其只纵向偏移3?时的罐容表

油位高度 0 100 200 300 400 500 600

油量(纵向和横向) 74.31491 383.669051 1002.000078 1988.832141 3287.159609 4883.629707 6715.762293

油位高度 1500 1600 1700 1800 1900 2000 2100 油量(纵向和横向) 29284.88497 31108.07876 34943.96173 37761.79985 40550.41694 43293.44186 45961.35349

700 800 900 1000 1100 1200 1300 1400

2横向偏转倾斜时

8745.41834 10943.80179 13286.75517 15752.90982 18322.72985 20977.93662 23701.13933 26475.57014

2200 2300 2400 2500 2600 2700 2800 2900 3000 48574.88207 510775.6497 53462.58567 55708.65462 57791.88932 59684.30884 61386.90111 65062.23466 64956.99592

当储油罐有横向偏转倾斜时，设h0为没有横向偏转倾斜时，油位探针应该显示的高度；设h为有横向偏转时，油位探针应该显示的高度。因为储油罐的主体为圆柱体，两端为球冠体，所以该储油罐具有对称结构。由此可得当储油罐横向转倾斜后，其油面的水平高度没有变化，只是油位探针的位置发生了变化。设R为圆柱的半径，如图所示为横向转倾斜后正截面图：

图15横向倾斜与无横向倾斜时的关系图

由平面几何关系可得：（1）当h0?R?1?cos??时，h=0；（2）当R?1?cos???h0?2R?Rcos?时，

h?cos??h0?R?1?cos??即：h?h0?R?1?cos??cos?；

（3）当2R?Rcos??h0时，横向转倾斜后的油位探针位置达到最大值2R，即

h=2R。

综上所述，没有横向偏转倾斜的油位探针应该显示的高度h0与有横向偏转倾斜的油位探针应该显示的高度h的关系如下：

?0??h?R?1?cos?h??0cos????2R?h0?R?1?cos??

R?1?cos???h0?2R?Rcos?2R?Rcos??h0在此题中，由题意可知该圆柱体的半径R=1.5（米）。所以此题中h和h0的关系为：

?0??h?1.5?1?cos?h??0cos????3?h0?1.5?1?cos??

1.5?1?cos???h0?3?1.5cos?3?1.5cos??h0当储油罐既有纵向倾斜又有横向偏转时，又上表达式可得：当1.5?1?cos???h0?3?1.5cos?时

h0?h?1.5?cos??1.5

设该储油罐体只纵向倾斜时的油位探针的显示高度为h1则：

?h0??h?1.5?cos??1.5?1 h1?h?0.4tan???h0?h1?解此方程组可得：h1??h?1.5?cos??1.5?0.4tan?

h?h?0.4tan??1.5cos?1?1.5

因为h?0，所以

h?0.4tan??1.5cos?1?1.5?0

解此方程可得：h1?0.4tan??1.5?1?cos?? 3当具有横向和纵向变位时

将只有纵向变位时油位探针显示高度h1与相应纵向变位时再横向变位?角的

?0.4ta?n带入变位时的油位探针显示高度h的对应关系h1??h?1.5?cos??1.5只有纵向倾斜时V=（h1）中即得到既有纵向变位也有横向变位时的油位探针的显示高度与储油罐所剩的体积的表达式如下： 1.当0.4tan??1.5?1?cos???h?7.6tan??1.5?cos??1?cos?时，

V?V1?211?0y0.2??h?1.5?cos??1.5?0.4?y?tan?1.5?R12R??0.625?y???z?1.5?dzdy

22?h?1.5?cos??1.5?0.4tan??tan?022??h?1.5?cos??1.5??0.4?y?tan?01.5??z?1.5?dzdy22

V1?V1?V11；

?h?1.5cos??1.4tan?cos2.当

17.6tan??1.5cos??1.5cos?时，

V2=V11

V2?V2?32?L02??h?1.5?cos??1.5?0.6tan?01.5??z?1.5?dzdy

22?y182??h?1.5?cos??1.5??0.4?y?tan?1.5?R223R??y?7.375???z?1.5?dzdy

222V2?V2?V2?V21；

?h?1.5cos??8.6tan?cos?23.当

V3?11.5cos??1.4tan?cos?时，

2?y0?1?R1dy?2?2?0y02??h?1.5?cos??1.5??0.4?y?tan?1.5?R1R1??z?1.5?dzdy

2V3?V3?32??L0y1?h?1.5?cos??1.5??0.4?y?tan?01.5??z?1.5?dzdy

2812??h?1.5?cos??1.5??0.4?y?tan?1.5?R223R??y?7.35???z?1.5?dzdy

222V3?V3?V3?V3；

?h?14.当

V4?V428.6tan??1.5cos?cos?1.5?1?cos???0.4tan?cos?2时，

211?y0?1?R1dy?2??y00?h?1.5?cos??1.5??0.4?y?tan?1.5?R1R1??z?1.5?dzdy?V3

＝V32

V4?3?y182??h?1.5?cos??1.5??0.4?y?tan?1.5?R2R??y?7.35???z?1.5?dzdy?222?9y12?R2dy

V4?V4?V4?V4123；

5.当1.5??h?1.5?cos??0.4tan??h1?1.5??h?1.5?cos??1.6tan?时，

V5?21?0?1?R1dy?22?0?122??R??0.625?y??dy

??V5???1.5?y1

V5?V5?43??Ly1y12??h?1.5?cos??1.5??0.4?y?tan?01.5??z?1.5?dzdy

22812??h?1.5?cos??1.5??0.4?y?tan?1.5?R2234R??y?7.35???z?1.5?dzdy?222?9y122??R??y?7.35??dy??V5?V5?V5?V5?V5。

在此模型中，我们取双变位??3?，??4.3?时，由编辑程序（见附录）制的油位高度间隔为10cm的罐容表如表3所示。

表3.纵向横向同时变位时，即??3，??4.3

油位高度/mm 0 100 200 300 400 500 600 700 800 900 1000 1100 1200 1300 1400 1500

4变位参数的确定

根据题目给出的罐体变位后在进/出油过程中的实际检测数据，我们通过matlab描出各对数据的对应点，如图17所示。从图上可以看出，这些点的连线大致接近于一条直线，如果要想确定纵向倾斜角度?和横向偏转角度?，必须采用一些辨识方法来对参数进行估计，这里我们就采用了线性最小二乘法原理。

我们假设罐体的油量的体积和油位高度的关系函数是V?f(h),

f(h)?ah?b，其中

油量/L 69.89519 377.7966 996.0623 1939.966 3241.934 4835.118 6656.727 8697.624 10898.89 13245.99 15717.34 18293.2 20955.14 23685.61 26467.68 29284.88 油位高度/mm 1600 1700 1800 1900 2000 2100 2200 2300 2400 2500 2600 2700 2800 2900 3000 油量/L 31198.0788 34959.9221 37745.2732 40581.6748 43331.7512 46018.8368 48625.3466 51132.7997 53521.4626 55769.814 57853.6804 59744.6812 61406.9011 65057.9031 64952.7523 a和b是待定常数。常数a和b确定了之后我们就可以根据

前面所建立的表达式来求出?和?。根据最小二乘法原理，要想使我们建立的函

7x 1046543210050010001500200025003000

图17

数尽可能过这些点，必须使函数值和实际数据的偏差最小，即：vi?f(hi)很小。因为每个偏差都很小，其中有正有负，所以我们对偏差取绝对值求和来保证每个

n偏差都很小，即：?|vi?f(hi)|很小。

i?1 由于取绝对值不便于进一步分析讨论，而任何实数的平方都是正数或零，所以我们建立了根据偏差的平方和为最小的条件来选择常数a、b。

nniM??(vi?1?f(hi))2??(vi?1i?(ahi?b))2

利用多元函数的极值及最大值、最小值原理，我们把求偏差的平方和最小归结为求函数

M?M(a,b)在哪点取得最小值。方程组为：

??Ma?a,b??0 ???Mb?a,b??0的解来解决，即令：

n??M??2???Vi??ahi?b???hi?0??a?i?1 ?n??M??2???Vi??ahi?b????0??bi?1?亦即

?n?Vi??ahi?b????0??hi??i?1 ?n????Vi??ahi?b????0??i?1将括号内各项进行整理合并，并把未知数a和b分离出来，便得

nn?n2?a?hi?b?hi??hiVi?i?1i?1i?1?nn?a?hi?8b??Vi?i?1i?1?

b?16.77679。然后通过matlab编程解上试方程组，得到a?231488.45176，

在对模型二中的表达式求积分,得到一个函数表达式,根据上面求出的a和b来对

应的求出?和?，其中??2.1?,??4.2?.并且当???4.2?时?都不变，因为横向左右偏转角度一样结果不变。由于当油位高度和体积不同时，所求的?和?也不同，我们做了多次计算，然后对结果进行了分析，发现当纵向倾斜角度增大时，偏转角度逐渐下降，它们是成反向的关系。同时通过前面对偏转倾斜所建立的模型可以看出来，它们之间的关系也是存在此种关系的。当然由于最小二乘法估计参数的实时性较差，所以会出现一些误差，导致计算出来的数据不是很精确，可能用另外的方法估计参数会使误差更小，这也是我们模型的不足之处。

7 模型的评价

模型的优点：

（1）本文模型充分利用了微元思想和积分的方法，将储油分割成许多微小薄面单元，从中选取一个小单元进行研究，再对其积分找出整体的规律。

（2）模型在MATLAB环境下计算，并对问题一中的实际数据与理论数据的拟合曲线比较，得到误差范围0.014866?E1?0.034917，充分提高了模型结果的准确性与精确性。问题一中当纵向倾斜??4.1?时，由模型一中的函数关系算出的理论数值与其记录的实际数值一一对应的相对误差1.42%?E2?3.37%。其误差值相对较小。

（3）在问题二中我们采用了分段考虑、合理切截和微积分方法进行了模型的建立，用递推的最小二乘法进行参数估计。模型的不足：

本文模型忽略环境的改变对储油体积的影响以及油位探针，注油口管，出油管对储油罐中剩余油体积的影响，从而使得模型的误差变大。模型的模拟实验不可能和实际完全一致，因为很多未量化数据收集误差等因素因此不能完全准确的进行分析。模型仍然需要修正和完善。

8 模型的改进与推广

本文假设实验环境的温度、压强、密度保持不变，但实际应用中环境的温度以及压强随时间不断改变的，因此储油罐中油的体积在不同时间会随着温度等其它未经量化的因素的影响而不同。 8.1 模型的改进

将环境的温度变化对储油量的体积的影响关系引入，从而减少理论数据与实际数据的误差。比如汽油的胀系数B1=12×T(1/℃) 其中，即气温每上升1度，体积增大为原来的1.0012。 8.2模型的推广

本文通过设计储油罐的变位识别与罐容表标定模型有效的模拟出储油量与油位高度以及变位参数之间关系的数学模型，具有很好的预测和拟合效果，因此有较强的现实意义，可以在现实生活中广泛应用。

9 参考文献

[1]刘卫国.MATLAB程序设计与应用.北京：高等教育出版社，2006年第二版. [2]同济大学应用数学系编.高等数学.北京：高等教育出版社，2006年第三版. [3]华东师范大学数学系编.数学分析.北京：高等教育出版社，2008年第三版. [4]韩中庚.数学建模竞赛获奖论文精选与点评.北京：科学出版社，2007.

[5]用逼近法计算横截面为椭圆形圆形储油罐的储油体积，百度文库，[2010-9-11].http://wenku.http://www.wodefanwen.com//view/3ecedebfc77da26925c5b00f.html

程序：

附录一：

用来求无变位进/出油的体积

%进油油位高度

hi=[159.02,176.14,192.59,208.5,223.93,238.97,253.66,268.04,282.16,296.03,309.69,323.15,...

336.44,349.57,362.56,375.42,388.16,400.79,413.32,425.76,438.12,450.4,462.62,474.78,...

486.89,498.95,510.97,522.95,534.9,546.82,558.72,570.61,582.48,594.35,606.22,618.09,...

629.96,641.85,653.75,665.67,677.63,678.54,690.53,690.82,702.85,714.91,727.03,739.19,...

751.42,763.7,764.16,776.53,788.99,801.54,814.19,826.95,839.83,852.84,866,879.32,...

892.82,892.84,906.53,920.45,934.61,949.05,963.8,978.91,994.43,1010.43,1026.99,1044.25,...

1062.37,1081.59,1102.33,1125.32,1152.36,1193.49]; %累加进油量

xi=[50,100,150,200,250,300,350,400,450,500,550,600,650,700,750,800,850,900,950,1000,1050,...

1100,1150,1200,1250,1300,1350,1400,1450,1500,1550,1600,1650,1700,1750,1800,1850,1900,...

1950,2000,2050,2053.83,2103.83,2105.06,2155.06,2205.06,2255.06,2305.06,2355.06,2404.98,...

2406.83,2456.83,2506.83,2556.83,2606.83,2656.83,2706.83,2756.83,2806.83,2856.83,2906.83,...

2906.91,2956.91,3006.91,3056.91,3106.91,3156.91,3206.91,3256.91,3306.91,3356.91,3406.91,...

3456.91,3506.91,3556.91,3606.91,3656.91,3706.91]; % 累加出油量

xo=[52.72,102.72,152.72,202.72,252.72,302.72,352.72,402.72,452.72,502.72,552.72,602.72,652.72,...

702.72,752.72,802.72,852.72,902.72,952.72,1002.72,1052.72,1102.72,1152.72,1202.72,1252.72,...

1302.72,1352.72,1402.72,1452.72,1502.72,1552.72,1602.72,1652.72,1702.72,1752.72,1802.72,...

1852.72,1902.72,1952.72,2002.72,2052.72,2102.72,2152.72,2202.72,2252.72,2302.72,2352.72,...

2402.72,2452.72,2502.72,2552.72,2602.72,2652.72,2702.72,2752.72,2802.72,2852.72,2902.72,...

2952.72,3002.72,3052.72,3102.72,3152.72,3202.72,3252.72,3302.72,3352.72,3402.72,3452.72,...

3502.72,3552.72,3602.72,3652.72,3702.72]; %出油油位高度

ho=[1150.72,1123.99,1101.15,1080.51,1061.36,1043.29,1026.08,1009.54,993.57,978.08,962.99,948.26,...

933.84,919.69,905.78,892.1,878.61,865.3,852.15,839.14,826.27,813.52,800.87,788.33,775.88,...

763.51,751.21,738.98,726.81,714.7,702.64,690.61,678.63,666.68,654.75,642.84,630.96,619.08,...

607.21,595.35,583.48,571.61,559.72,547.82,535.9,523.95,511.97,499.96,487.9,475.8,463.65,...

451.43,439.15,426.8,414.36,401.84,389.22,376.49,363.64,350.67,337.55,324.27,310.82,297.18,...

283.33,269.24,254.88,240.21,225.21,209.81,193.94,177.54,160.48,142.62];

hi=hi/1000; ho=ho/1000; xi=xi+262;

xo=xi(end)-xo;

a=0.89;b=0.6;l=2.45; s=[];

for i=1:length(hi)

s=[s 2*a/b*quad('sqrt(0.6^2-y.^2)',-b,hi(i)-b)]; end

v=s*l*1000;

p1=polyfit([hi*1000 ho*1000],[xi xo],5); p2=polyfit(hi*1000,v,5); x1=0:1300;

y1=polyval(p1,x1);

y2=polyval(p2,x1);

plot(hi*1000,v,'b.',x1,y2,'b',[hi*1000 ho*1000],[xi xo],'r*',x1,y1,'r');

legend('理论数据','理论拟合曲线','实验数据','实验拟合曲线');

%求误差

v2=abs(v-xi)./xi;

ave=sum(v2)/length(v2); m=max(v2); n=min(v2);

disp(['max=',num2str(m),'min=',num2str(n),'average=',num2str(ave)]);

附录二：

求纵向倾斜时的储油罐内油量的体积和高度间隔为1cm的罐容表标定值：

syms y;

a=0.89;b=0.6;l=2.45;

xi=[747.86,797.86,847.86,897.86,947.86,997.86,1047.86,1097.79,1147.79,1197.73,1247.73,1297.73,1347.73,1397.73,1447.73,1497.73,1547.73,1597.73,1647.73,1697.73,1747.73,1797.73,1847.73,1897.73,1947.73,1997.73,2047.73,2097.73,2147.73,2197.73,2247.73,2297.73,2347.73,2397.73,2447.73,2497.73,2547.73,2597.73,2647.73,2697.73,2747.73,2797.73,2847.73,2897.73,2947.73,2997.73,3047.73,3097.73,3147.73,3197.73,3247.73,3297.73,3299.74];

hi=[411.29,423.45,438.33,450.54,463.9,477.74,489.37,502.56,514.69,526.84,538.88,551.96,564.4,576.56,588.74,599.56,611.62,623.44,635.58,646.28,658.59,670.22,680.63,693.03,704.67,716.45,727.66,739.39,750.9,761.55,773.43,785.39,796.04,808.27,820.8,832.8,844.47,856.29,867.6,880.06,892.92,904.34,917.34,929.9,941.42,954.6,968.09,980.14,992.41,1006.34,1019.07,1034.24,1035.36];

xo=[50,100,150,200,250,300,350,400,450,500,550,600,650,700,750,800,850,900,950,1000,1050,1100,1150,1200,1250,1300,1350,1400,1450,1500,1550,1600,1650,1700,1750,1800,1850,1900,1950,2000,2050,2100,2150,2200,2250,2300,2350,2400,2450,2500,2550];

ho=[1020.65,1007.73,994.32,980.96,967.1,956.01,941.54,929.69,916.44,904.14,891.9,879.23,868.99,855.13,844.02,831.64,820.47,808.16,796,785.04,773.07,762.09,750.81,739.42,727.09,715.32,705.43,693.52,682.5,671.02,658.68,647.74,635.76,624.61,612.53,600.69,589.4,577,564.58,554.33,540.76,528.65,517.19,504.87,490.78,478.06,465.97,452.4,439.98,425.83,411.73]; xi=xi+215; hi=hi/1000;

h1=hi+0.4*tan(4.1*pi/180); vi=[];

%进油量的理论值 for i=1:length(h1)

f=(b^2/2*asin((-y*tan(4.1*pi/180)+h1(i)-b)/b)+((-y*tan(4.1*pi/180)+h1(i)-b)/2)*sqrt(b^2-(-y*tan(4.1*pi/180)+h1(i)-b)^2)+b^2*pi/4)*2*a/b; d=int(f,0,2.45); vi=[vi eval(d)]; end

vi=vi*1000; vi'

ho=ho/1000;

h2=ho+0.4*tan(4.1*pi/180); vo=[];

%出油量的理论值 for i=1:length(h2)

f=(b^2/2*asin((-y*tan(4.1*pi/180)+h2(i)-b)/b)+((-y*tan(4.1*pi/180)+h2(i)-b)/2)*sqrt(b^2-(-y*tan(4.1*pi/180)+h2(i)-b)^2)+b^2*pi/4)*2*a/b; d=int(f,0,2.45); vo=[vo eval(d)]; end

vo=vi(end)-vo*1000; vo'

p1=polyfit([hi*1000 ho*1000],[xi xi(end)-xo],5); p2=polyfit([hi*1000 ho*1000],[vi vi(end)-vo],5); x1=0:1300;

y1=polyval(p1,x1); y2=polyval(p2,x1); plot([hi*1000 ho*1000],[xi xi(end)-xo],'b.',x1,y1,'b',[hi*1000 ho*1000],[vi vi(end)-vo],'r*',x1,y2,'r');

legend('实际数据','实际数据拟合曲线','理论数据','理论数据拟合曲线');

ti=abs(vi-xi)./xi; to=abs(vo-xo)./xo; figure; ti' to'

p=polyfit(hi,ti,5); x=1:1.5;

y1=polyval(p,x);

plot(hi,ti,x,y1); average=mean(ti); s=std(ti);

h2=0.1+0.4*tan(4.1*pi/180);

f=(b^2/2*asin((-y*tan(4.1*pi/180)+h2-b)/b)+((-y*tan(4.1*pi/180)+h2-b)/2)*sqrt(b^2-(-y*tan(4.1*pi/180)+h2-b)^2)+b^2*pi/4)*2*a/b; d=int(f,0,h2/tan(4.1*pi/180)); eval(d);

vv=[];

hii=0.15:0.01:1.17;

h11=hii+0.4*tan(4.1*pi/180); for i=1:length(h11)

f=(b^2/2*asin((-y*tan(4.1*pi/180)+h11(i)-b)/b)+((-y*tan(4.1*pi/180)+h11(i)-b)/2)*sqrt(b^2-(-y*tan(4.1*pi/180)+h11(i)-b)^2)+b^2*pi/4)*2*a/b;

d=int(f,0,2.45); vv=[vv eval(d)]; end

vv=vv*1000; vv'

vvv=[];

hiii=0:0.01:0.15;

h111=hiii+0.4*tan(4.1*pi/180); for i=1:length(h111)

f=(b^2/2*asin((-y*tan(4.1*pi/180)+h111(i)-b)/b)+((-y*tan(4.1*pi/180)+h111(i)-b)/2)*sqrt(b^2-(-y*tan(4.1*pi/180)+h111(i)-b)^2)+b^2*pi/4)*2*a/b;

d=int(f,0,h111(i)/tan(4.1*pi/180)); vvv=[vvv eval(d)]; end

vvv=vvv*1000; format long; vvv'

vvvv=[];

hiiii=1.18:0.01:1.2;

h1111=1.2-hiiii+2.05*tan(4.1*pi/180); for i=1:length(h1111)

f=(b^2/2*asin((-y*tan(4.1*pi/180)+h1111(i)-b)/b)+((-y*tan(4.1*pi/180)+h1111(i)-b)/2)*sqrt(b^2-(-y*tan(4.1*pi/180)+h1111(i)-b)^2)+b^2*pi/4)*2*a/b;

d=int(f,0,h1111(i)/tan(4.1*pi/180)); vvvv=[vvvv eval(d)]; end

vvvv=pi*a*b*l*1000-vvvv*1000; format long; vvvv'

附录三：

用来求纵向倾斜时的体积：

syms y;

a=3;R=1.625;

%第一部分体积 h=0:0.1:0.3;

h1=h+2*tan(a*pi/180); y0=[];

for i=1:length(h1)

f1=(y-0.625)^2+(-y*tan(a*pi/180)+h1(i)-1.5)^2-1.625^2; x0=solve(f1);

y0=[y0 min(eval(x0))]; end y0'

v11=[]; v12=[];

for i=1:length(h1) %求左球冠体积

f2=((R^2-(0.625-y)^2)/2*asin((h1(i)-y*tan(a*pi/180)-1.5)/(sqrt(R^2-(0.625-y)^2)))+((h1(i)-y*tan(a*pi/180)-1.5)/2)*sqrt(R^2-(0.625-y)^2-(h1(i)-y*tan(a*pi/180)-1.5)^2)+pi*(R^2-(0.625-y)^2)/4)*2; t=y0(i);

d2=int(f2,t,0);

v11=[v11 eval(d2)];

%求圆柱体积

f3=((1.5^2)/2*asin((-y*tan(a*pi/180)+h1(i)-1.5)/1.5)+((-y*tan(a*pi/180)+h1(i)-1.5)/2)*sqrt(1.5^2-(-y*tan(a*pi/180)+h1(i)-1.5)^2)+pi*1.5^2/4)*2;

d3=int(f3,0,(h1(i))/(tan(a*pi/180))); v12=[v12 eval(d3)]; end

v1=v11+v12;

disp('第一部分的体积'); v1'

%第二部分体积 h=0.4:0.1:1.3;

h1=h+2*tan(a*pi/180); y0=[];

%求水平面与左球冠的交点 for i=1:length(h1)

f1=(y-0.625)^2+(-y*tan(a*pi/180)+h1(i)-1.5)^2-1.625^2; x0=solve(f1);

y0=[y0 min(eval(x0))]; end %y0

v21=[]; v22=[]; v23=[];

for i=1:length(h1) %求左球冠体积

d2=int(f2,t,0);

v21=[v21 eval(d2)];

%求圆柱体积

f3=((1.5^2)/2*asin((-y*tan(a*pi/180)+h1(i)-1.5)/1.5)+((-y*tan(a*pi/180)+h1(i)-1.5)/2)*sqrt(1.5^2-(-y*tan(a*pi/180)+h1(i)-1.5)^2)+pi*1.5^2/4)*2;

d3=int(f3,0,8);

v22=[v22 eval(d3)];

%求水平面与右球冠的交点

f4=(y-7.375)^2+(-y*tan(a*pi/180)+h1(i)-1.5)^2-1.625^2; x00=solve(f4);

y00=max(eval(x00));

%求右球冠体积

f5=((R^2-(y-7.375)^2)/2*asin((h1(i)-y*tan(a*pi/180)-1.5)/(sqrt(R^2-(y-7.375)^2)))+((h1(i)-y*tan(a*pi/180)-1.5)/2)*sqrt(R^2-(y-7.375)^2/2-(h1(i)-y*tan(a*pi/180)-1.5)^2)+pi*(R^2-(y-7.375)^2)/4)*2; y00

d4=int(f5,8,y00); v23=[v23 eval(d4)]; end

v2=v21+v22+v23;

disp('第二部分的体积'); v2'

%第三部分体积 h=1.4:0.1:1.8;

h1=h+2*tan(a*pi/180); y0=[];

%求水平面与左球冠的交点 for i=1:length(h1)

f1=(y-0.625)^2+(-y*tan(a*pi/180)+h1(i)-1.5)^2-1.625^2; x0=solve(f1);

y0=[y0 min(eval(x0))]; end y0'

v311=[]; v312=[]; v31=[]; v32=[]; v33=[]; v3=[];

for i=1:length(h1) %求左球冠体积

d2=int(f2,t,0);

v311=[v311 eval(d2)];

f6=pi*(R^2-(0.625-y)^2); t=y0(i);

d5=int(f6,-1,t);

v312=[v312 eval(d5)];

%求圆柱体积

f3=((1.5^2)/2*asin((-y*tan(a*pi/180)+h1(i)-1.5)/1.5)+((-y*tan(a*pi/180)+h1(i)-1.5)/2)*sqrt(1.5^2-(-y*tan(a*pi/180)+h1(i)-1.5)^2)+pi*1.5^2/4)*2;

d3=int(f3,0,8);

v32=[v32 eval(d3)];

%求水平面与右球冠的交点

f4=(y-7.375)^2+(-y*tan(a*pi/180)+h1(i)-1.5)^2-1.625^2; x00=solve(f4);

y00=max(eval(x00));

%求右球冠体积

f5=((R^2-(y-7.375)^2)/2*asin((h1(i)-y*tan(a*pi/180)-1.5)/(sqrt(R^2-(y-7.375)^2)))+((h1(i)-y*tan(a*pi/180)-1.5)/2)*sqrt(R^2-(y-7.375)^2/2-(h1(i)-y*tan(a*pi/180)-1.5)^2)+pi*(R^2-(y-7.375)^2)/4)*2; y00

d4=int(f5,8,y00); v33=[v33 eval(d4)]; end

v31=v311+v312; v3=v31+v32+v33;

disp('第三部分的体积'); v3'

%第四部分体积

h=[1.9:0.1:2.8 2.63]; h1=h+2*tan(a*pi/180); y0=[];

%求水平面与左球冠的交点 for i=1:length(h1)

f1=(y-0.625)^2+(-y*tan(a*pi/180)+h1(i)-1.5)^2-1.625^2; x0=solve(f1);

y0=[y0 min(eval(x0))]; end %y0

v411=[]; v412=[];

v41=[]; v42=[]; v431=[]; v432=[]; v43=[]; v4=[];

for i=1:length(h1) %求左球冠体积

d2=int(f2,t,0);

v411=[v311 eval(d2)];

f6=pi*(R^2-(0.625-y)^2); t=y0(i);

d5=int(f6,-1,t);

v412=[v312 eval(d5)]; v41=v411+v412;

%求圆柱体积

f3=((1.5^2)/2*asin((-y*tan(a*pi/180)+h1(i)-1.5)/1.5)+((-y*tan(a*pi/180)+h1(i)-1.5)/2)*sqrt(1.5^2-(-y*tan(a*pi/180)+h1(i)-1.5)^2)+pi*1.5^2/4)*2;

d3=int(f3,0,8);

v42=[v42 eval(d3)];

%求水平面与右球冠的交点

f4=(y-7.375)^2+(-y*tan(a*pi/180)+h1(i)-1.5)^2-1.625^2; x00=solve(f4);

y00=max(eval(x00));

%求右球冠体积

f5=((R^2-(y-7.375)^2)/2*asin((h1(i)-y*tan(a*pi/180)-1.5)/(sqrt(R^2-(y-7.375)^2)))+((h1(i)-y*tan(a*pi/180)-1.5)/2)*sqrt(R^2-(y-7.375)^2/2-(h1(i)-y*tan(a*pi/180)-1.5)^2)+pi*(R^2-(y-7.375)^2)/4)*2; d4=int(f5,8,y00);

v431=[v431 eval(d4)];

f9=pi*(R^2-(y-7.375)^2);

d9=int(f9,y00,9);

v432=[v432 eval(d9)]; end

v43=v431+v432; v4=v41+v42+v43;

disp('第四部分的体积'); v4'

%第五部分体积 h=2.9:0.1:3;

h1=h+2*tan(a*pi/180); y0=[];

%求水平面与左球冠的交点 for i=1:length(h1)

f1=(y-0.625)^2+(-y*tan(a*pi/180)+h1(i)-1.5)^2-1.625^2; x0=solve(f1);

y0=[y0 min(eval(x0))]; end %y0

v51=[]; v52=[]; v53=[]; v541=[]; v542=[]; v54=[]; v5=[];

for i=1:length(h1) %求左球冠体积

f6=pi*(R^2-(0.625-y)^2); d5=int(f6,-1,0); v51=[v51 eval(d5)];

%求圆柱体积

f4=(y-7.375)^2+(-y*tan(a*pi/180)+h1(i)-1.5)^2-1.625^2;%求水平面与右球冠的交点

x00=solve(f4);

y00=max(eval(x00));

v52=[v52 pi*1.5^2*y00];

f3=((1.5^2)/2*asin((-y*tan(a*pi/180)+h1(i)-1.5)/1.5)+((-y*tan(a*pi/180)+h1(i)-1.5)/2)*sqrt(1.5^2-(-y*tan(a*pi/180)+h1(i)-1.5)^2)+pi*1.5^2/

4)*2;

d3=int(f3,y00,8); v53=[v53 eval(d3)];

%求右球冠体积

v541=[v541 eval(d4)];

f9=pi*(R^2-(y-7.375)^2); d9=int(f9,y00,9);

v542=[v542 eval(d9)]; end

v54=v541+v542;

v5=v51+v52+v53+v54;

disp('第五部分的体积'); v5'

附录四：

用来求纵向和横向同时变位的第一部分的体积：

syms y;

a=3;R=1.625;b=4.3;

%第一部分体积 h=0:0.1:0.3;

h1=h*cos(b*pi/180)+1.5*(1-cos(b*pi/180))+2*tan(a*pi/180); y0=[];

for i=1:length(h1)

f1=(y-0.625)^2+(-y*tan(a*pi/180)+h1(i)-1.5)^2-1.625^2; x0=solve(f1);

y0=[y0 min(eval(x0))]; end y0'

v11=[]; v12=[];

for i=1:length(h1) %求左球冠体积

d2=int(f2,t,0);

v11=[v11 eval(d2)];

%求圆柱体积

f3=((1.5^2)/2*asin((-y*tan(a*pi/180)+h1(i)-1.5)/1.5)+((-y*tan(a*pi/180)+h1(i)-1.5)/2)*sqrt(1.5^2-(-y*tan(a*pi/180)+h1(i)-1.5)^2)+pi*1.5^2/4)*2;

d3=int(f3,0,(h1(i))/(tan(a*pi/180))); v12=[v12 eval(d3)]; end

v1=v11+v12;

disp('第一部分的体积'); format long; v1'

附录五：

用来求纵向和横向同时变位的第二部分的体积：

syms y;

a=3;R=1.625;b=4.3;

%第二部分体积 h=0.4:0.1:1.3;

h1=h*cos(b*pi/180)+1.5*(1-cos(b*pi/180))+2*tan(a*pi/180); y0=[];

%求水平面与左球冠的交点 for i=1:length(h1)

f1=(y-0.625)^2+(-y*tan(a*pi/180)+h1(i)-1.5)^2-1.625^2; x0=solve(f1);

y0=[y0 min(eval(x0))]; end y0

v21=[]; v22=[]; v23=[];

for i=1:length(h1) %求左球冠体积

d2=int(f2,t,0);

v21=[v21 eval(d2)];

%求圆柱体积

f3=((1.5^2)/2*asin((-y*tan(a*pi/180)+h1(i)-1.5)/1.5)+((-y*tan(a*pi/180)+h1(i)-1.5)/2)*sqrt(1.5^2-(-y*tan(a*pi/180)+h1(i)-1.5)^2)+pi*1.5^2/4)*2;

d3=int(f3,0,8);

v22=[v22 eval(d3)];

%求水平面与右球冠的交点

f4=(y-7.375)^2+(-y*tan(a*pi/180)+h1(i)-1.5)^2-1.625^2; x00=solve(f4);

y00=max(eval(x00));

%求右球冠体积

f5=((R^2-(y-7.375)^2)/2*asin((h1(i)-y*tan(a*pi/180)-1.5)/(sqrt(R^2-(y-7.375)^2)))+((h1(i)-y*tan(a*pi/180)-1.5)/2)*sqrt(R^2-(y-7.375)^2/2-(h1(i)-y*tan(a*pi/180)-1.5)^2)+pi*(R^2-(y-7.375)^2)/4)*2; y00

d4=int(f5,8,y00); v23=[v23 eval(d4)]; end

v2=v21+v22+v23;

disp('第二部分的体积'); format long; v2'

附录六：

用来求纵向和横向同时变位的第三部分的体积：

syms y;

a=3;R=1.625;b=4.3;

%第三部分体积 h=1.4:0.1:1.8;

h1=h*cos(b*pi/180)+1.5*(1-cos(b*pi/180))+2*tan(a*pi/180);

y0=[];

%求水平面与左球冠的交点 for i=1:length(h1)

f1=(y-0.625)^2+(-y*tan(a*pi/180)+h1(i)-1.5)^2-1.625^2; x0=solve(f1);

y0=[y0 min(eval(x0))]; end y0'

v311=[]; v312=[]; v31=[]; v32=[]; v33=[]; v3=[];

for i=1:length(h1) %求左球冠体积

d2=int(f2,t,0);

v311=[v311 eval(d2)];

f6=pi*(R^2-(0.625-y)^2); t=y0(i);

d5=int(f6,-1,t);

v312=[v312 eval(d5)];

%求圆柱体积

f3=((1.5^2)/2*asin((-y*tan(a*pi/180)+h1(i)-1.5)/1.5)+((-y*tan(a*pi/180)+h1(i)-1.5)/2)*sqrt(1.5^2-(-y*tan(a*pi/180)+h1(i)-1.5)^2)+pi*1.5^2/4)*2;

d3=int(f3,0,8);

v32=[v32 eval(d3)];

%求水平面与右球冠的交点

f4=(y-7.375)^2+(-y*tan(a*pi/180)+h1(i)-1.5)^2-1.625^2; x00=solve(f4);

y00=max(eval(x00));

%求右球冠体积

f5=((R^2-(y-7.375)^2)/2*asin((h1(i)-y*tan(a*pi/180)-1.5)/(sqrt(R^2-(y-7.375)^2)))+((h1(i)-y*tan(a*pi/180)-1.5)/2)*sqrt(R^2-(y-7.375)^2/2-(h1(i)-y*tan(a*pi/180)-1.5)^2)+pi*(R^2-(y-7.375)^2)/4)*2; y00

d4=int(f5,8,y00); v33=[v33 eval(d4)]; end

v31=v311+v312; v3=v31+v32+v33;

disp('第三部分的体积'); format long; v3'

附录七：

用来求纵向和横向同时变位的第四部分的体积：

syms y;

a=3;R=1.625;b=4.3;

%第四部分体积 h=[1.9:0.1:2.8];

h1=h*cos(b*pi/180)+1.5*(1-cos(b*pi/180))+2*tan(a*pi/180); y0=[];

%求水平面与左球冠的交点 for i=1:length(h1)

f1=(y-0.625)^2+(-y*tan(a*pi/180)+h1(i)-1.5)^2-1.625^2; x0=solve(f1);

y0=[y0 min(eval(x0))]; end y0

v411=[]; v412=[]; v41=[]; v42=[]; v431=[]; v432=[]; v43=[]; v4=[];

for i=1:length(h1) %求左球冠体积

d2=int(f2,t,0);

v411=[v411 eval(d2)];

f6=pi*(R^2-(0.625-y)^2); t=y0(i);

d5=int(f6,-1,t);

v412=[v412 eval(d5)]; v41=v411+v412;

%求圆柱体积

f3=((1.5^2)/2*asin((-y*tan(a*pi/180)+h1(i)-1.5)/1.5)+((-y*tan(a*pi/180)+h1(i)-1.5)/2)*sqrt(1.5^2-(-y*tan(a*pi/180)+h1(i)-1.5)^2)+pi*1.5^2/4)*2;

d3=int(f3,0,8);

v42=[v42 eval(d3)];

%求水平面与右球冠的交点

f4=(y-7.375)^2+(-y*tan(a*pi/180)+h1(i)-1.5)^2-1.625^2; x00=solve(f4);

y00=max(eval(x00));

%求右球冠体积

f5=((R^2-(y-7.375)^2)/2*asin((h1(i)-y*tan(a*pi/180)-1.5)/(sqrt(R^2-(y-7.375)^2)))+((h1(i)-y*tan(a*pi/180)-1.5)/2)*sqrt(R^2-(y-7.375)^2/2-(h1(i)-y*tan(a*pi/180)-1.5)^2)+pi*(R^2-(y-7.375)^2)/4)*2; y00

d4=int(f5,8,y00);

v431=[v431 eval(d4)];

f9=pi*(R^2-(y-7.375)^2); d9=int(f9,y00,9);

v432=[v432 eval(d9)]; end

v43=v431+v432; v4=v41+v42+v43;

disp('第四部分的体积'); v4'

附录八：

用来求纵向和横向同时变位的第五部分的体积：

syms y;

a=3;R=1.625;b=4.3;

%第五部分体积 h=2.9:0.1:3;

h1=h*cos(b*pi/180)+1.5*(1-cos(b*pi/180))+2*tan(a*pi/180); y0=[];