储油罐的变位识别与罐容表标定大学生数学建模一等奖—-毕业论文设计.doc

储油罐的变位识别与罐容表标定
摘要
本文对储油罐的变位识别和罐容表的标定问题进行了深入探讨，建立了储油量和油位高度以及变位参数之间关系的数学模型，主要应用了mtalab进行求解。
针对问题一，我们利用积分的方法推导出小椭圆储油罐在无变位和发生纵向倾斜变位时的一般公式。讨论了在储油罐发生纵向倾斜变位后对罐容表的影响，定义了平均影响率(变位前后储油量之差绝对值的平均值占总罐体容积的比例)作为评价罐体变位对罐容表的影响程度的大小的指标，求出。并分别给出了小椭圆储油罐在无变位和在纵向倾斜变位角取的罐容表。
表1 小椭圆储油罐罐容表(纵向变位)
油位高度
0.00
0.01
0.02
0.03
0.04
0.05
…
储油体积
00.0017
0.0035
0.0063
0.0100
0.0148
0.0207
…
…
1.15
1.16
1.17
1.1713
1.18
1.19
1.20
…
3.9103
3.9338
3.9560
3.9588
3.9766
3.9955
4.10174.1101
针对问题二，将储油罐分为5个区域分别进行讨论，考虑到在球冠处的体积表达式过于复杂，我们省略了球冠处的一小部分体积，进行了近似求解，得出了罐内储油量与油位高度以及变位参数之间的一般关系的数学模型。
在利用储油罐的实际测量值估计变位参数时，我们建立了最小二乘拟合模型，得到了最佳的变位参数为：纵向倾斜变位，横向偏转变位。并据此对储油罐的罐容表进行了标定(见表3)。
在模型验证中，我们又采用蒙特卡洛模拟的方法对在问题二的模型中忽略的部分球冠体积进行了模拟计算。得到用问题二模型中求出的总储油量与模拟得出的总储油量一致度达到了99%，误差非常小，验证了我们所建立的模型的合理性和准确性。

关键词 平均影响率 最小二乘参数估计法 蒙特卡洛模拟

一 问题重述
通常加油站都有若干个储存燃油的地下储油罐，并且一般都有与之配套的“油位计量管理系统”，采用流量计和油位计来测量进/出油量与罐内油位高度等数据，通过预先标定的罐容表（即罐内油位高度与储油量的对应关系）进行实时计算，以得到罐内油位高度和储油量的变化情况。
许多储油罐在使用一段时间后，由于地基变形等原因，使罐体的位置会发生纵向倾斜和横向偏转等变化（以下称为变位），从而导致罐容表发生改变。按照有关规定，需要定期对罐容表进行重新标定。
问题一
为了掌握罐体变位后对罐容表的影响，利用小椭圆型储油罐（两端平头的椭圆柱体），分别对罐体无变位和倾斜角为a=4.1的纵向变位两种情况做了实验。请建立数学模型研究罐体变位后对罐容表的影响，并给出罐体变位后油位高度间隔为1cm的罐容表标定值。
问题二
对于实际的储油罐，试建立罐体变位后罐内储油量与油位高度及变位参数（纵向倾斜角度a和横向偏转角度b ）之间的一般关系。利用罐体变位后在进/出油过程中的实际检测数据，根据你们所建立的数学模型确定变位参数，并给出罐体变位后油位高度间隔为10cm的罐容表标定值。然后进一步用实际检测数据来分析检验模型的正确性与方法的可靠性。
二 模型假设
1、假设题中所给数据均为储油罐内壁测量值；
2、不考虑由于温度、压强变化等原因而引起储油罐的体积变化；
3、油位探针被固定在储油罐上，其上油浮子能够准确测量油位高度；

三 符号说明
椭圆的半长轴长
椭圆的半短轴长
储油罐的总长
油位探针到油罐底部左侧的距离
储油罐的总体积
储油罐的纵向倾斜角度
储油罐的横向偏转角度
油位高度
球冠体的半径
圆柱体的底面半径
注：未说明符号在文中用到时注明
四 问题一的解答
小椭圆储油罐罐体变位前后都可以应用积分的方法求出罐体的储油量和油位高度之间的关系。对于纵向倾斜的小椭圆储油罐，考虑分段求出其储油量和油位高度之间的关系，从而得到重新标定后的罐容表。
4.1 小椭圆储油罐无变位时的模型
由于此时的椭圆无变位，考虑先对二维椭圆进行积分。为方便表示油位高度，建立如图所示的坐标系，椭圆的半长轴长为，半短轴长为，则椭圆方程为

图1 对椭圆的积分示意图
在方向上取椭圆面中的一微元积分得到油的侧面积

储油罐内油的体积为

查积分表得到
(1)
利用matlab计算得到

经验证两种方法得到的体积公式完全等价，(1)式即为小椭圆储油罐无变位时的储油量和油位高度关系的模型。
根据此模型，我们可以求出小椭圆储油罐无变位时罐容表标定值(油位高度间隔取，结果见附录一)。
4.2 小椭圆储油罐纵向倾斜变位时的模型
储油罐纵向倾斜之后，油位计在油位过高或者过低时将不起作用(如图2所示的和区域)，考虑到倾斜角变化一般不会很大，所以我们可以将储油罐按液面高低分成五个部分，来求其储油量和油位高度之间的关系。我们讨论的是小椭圆储油罐纵向倾斜变位为逆时针旋转，如图2。对于储油罐顺时针旋转变位(即为负值)时的情况与此非常类似，在此不再详细讨论。

图2 储油罐分区示意图
4.2.1 对区域的讨论
在区域，其油位低于油位探针的油浮子，所以油位计量系统中显示油位高度为零。
当油位计刚开始有示数时，计算其储油体积。将区域放大得到图3

图3 区域的放大图
图中，从原点纸面向里为轴，利用三重积分可以得到

其中为油位探针到储油罐左侧的距离
积分得到
(2)
4.2.2 对区域的讨论
由区域很容易得到区域的储油量和油位高度的变化关系，直接给出结论：

所以
(3)
4.2.3 对区域的讨论

图4 区域示意图
在小椭圆储油罐无变位模型中我们已经求出了的计算公式，同区域中的积分原理可以计算出，我们就可以得到此时的油量体积为
(4)
其中

4.2.4 对区域的讨论
由区域4和区域2的相似性，将(3)式中的换为，将换为，并用总体积减去即为区域4的储油体积和油位高度的变化关系。

其中为小椭圆储油罐的总体积
化简并积分可得
(5)
4.2.5 对区域的讨论
在此区域中油浮子到达油位探针顶点，无法进一步测量油位高度。无法测量的总体积为：
(6)
4.2.6 综合各区域的罐容表标定的数学模型
综上所述，我们得到了储油量和油位高度、纵向倾斜角之间的分段函数关系式：
表2 分段函数关系
区域
油位高度
储油量
1

2

3

4

5

根据储油量和油位高度的分段函数关系我们可以得到罐体纵向倾斜变位()后油位高度间隔为的罐容表标定值(见附录一)。
当时，各区域油位高度及体积变化范围为：
表3 各区域油位高度和储油量变化范围
区域
油位高度
储油量
1

2

3

4

5

4.3 罐体变位后对罐容表的影响
为了能更加准确地刻画罐体的纵向倾斜变位对罐容表的影响，我们分别对罐体变位前后的理论值和测量值进行多方面的比较。
4.3.1 罐体变位前理论值与测量值比较
根据附件一中所给数据，我们计算出在附件中所给的油位高度下理论值和测量值，并画出其曲线。

图5 罐体变位前的曲线对比
通过对比我们发现，对于任意，储油量的理论值和实际值始终成如(7)式的比例关系。
(7)
4.3.2 罐体变位后理论值与测量值比较
同样根据附件一中所给数据，我们计算出在附件中所给的油位高度下理论值和测量值，并画出其曲线，如图6。

图6 罐体变位后的曲线对比
从图6中可以看出测量值仍然始终小于理论值，进一步求得理论值与测量值之差的变化范围为[0.0454,0.0910]，测量值的相对误差范围为[1.56%，5.18%]。
4.3.3 罐体纵向倾斜变位前后理论值比较

图7 罐体变位前后的曲线对比

图8 同一高度下储油量的理论值与测量值之差变化关系

图9 储油量的测量值的相对误差随油位高度的变化关系
由以上各图可以清晰地看出纵向倾斜变位后，使得在同一个油位高度下，变位后比变位前的储油量减小。但是这样仍不够直观，我们需要找到一个指标来定量刻画罐体变位后对罐容表的影响。从图9中可以看出，当油位高度较小时(0.1附近)，变位后相对于变位前的相对误差几乎达到了60%以上，但是此时的储油量的差别并不大，鉴于此，我们定义平均影响率：

来刻画罐体变位后对罐容表的影响。可以求出在纵向倾斜变位时，4.87%。
五 问题二的解答
如图8实际的储油罐示意图，其主体为圆柱体，两端为球冠体，在储油罐无变位时我们计算其各部分体积。

图10 实际储油罐示意图
圆柱体积计算公式

带入数据得到：主体圆柱体积为56.5487
一端球缺体积计算公式

带入数据得到：两端球缺总体积为，则储油罐的总体积为64.6645
5.1考虑不发生变位时储油量和油位高度的关系
圆柱内的油体积随油位高度变化关系：
(8)
一端球缺内油体积随油位高度变化关系：
(9)

其中
为圆柱体底面半径
为球缺对应的半径
为球内小圆半径
总的储油罐内的油量对油位高度的变化关系为：

5.2只考虑横向偏转变位时的储油量和油位高度的关系
只考虑储油罐横向偏转为时：对实际的油位高度没有影响，但此时的油位探针已经随储油罐发生偏转(如下图)

图11 只考虑横向偏转示意图

由油浮子测量得到的油位高度与实际油位高度的关系为：
(10)
所以对罐容表影响转换公式：

5.3只考虑纵向倾斜变位时的储油量和油位高度的关系
只考虑纵向倾斜变位时，我们利用问题一中的思想，将储油罐分成5个区域，分别计算储油量和油位高度的关系(如图12所示)。

图12 实际储油罐分区示意图
5.3.1对区域1的讨论
油位探针测得的油位高度始终为零，分成三部分来计算这部分的体积，如区域1的放大图(13)

图13 实际储油罐分区示意图

区域1的总体积应为：

对圆柱体部分进行三重积分得到：

由球缺部分的体积随油位高度的计算公式(9)得：

由于部分的体积精确计算非常复杂，而且储油罐的纵向倾斜角度一般不会高于所我们考虑将这部分体积省略，进行近似计算。由于倾斜角较小，所以区域3占据了储油罐的绝大部分，而在区域3中这种近似计算的误差将由于左右两个球冠的省略体积一正一负而有所减小，所以

由此得到区域1的总体积公式：
(11)
5.3.2 对区域2的讨论
在区域2中，油位高度的变化范围。
各部分储油体积和油位高度的变化关系:

],
总储油量和油位高度变化关系：
(12)
5.3.3 对区域3的讨论
在区域3中，油位高度的变化范围。
圆柱体部分的储油量：

其中

球缺部分的近似储油量：
，
总储油量和油位高度变化关系：
(13)
5.3.4 对区域4的讨论
在区域4中，油位高度的范围。
圆柱体部分的储油量：

球缺部分的总体积：
，
总储油量和油位高度变化关系：
(14)
5.3.5 对区域5的讨论
在区域4中，油位高度始终等于。
各部分体积分别为：

，
总储油量和油位高度变化关系：

5.4综合考虑储油罐纵向倾斜和横向偏转
根据如上讨论，我们可以得出结论，可以直接把5.3中各区域的公式用式(10)进行变换即可得到综合考虑了储油罐纵向倾斜和横向偏转的一般关系式。

我们依然考虑将储油罐分成五个区域分别求解。
5.4.1 对区域1的讨论
当储油罐未发生偏转时，油位高度为，而当储油罐发生横向偏转后，就可能使得油位探针测得的示数变为零。

图14 对区域一的讨论示意图
即当

时，储油罐发生横向偏转后，油位探针测得的油位高度恰好为零。本属于第二区域的部分横向偏转变位后转为了变位后的第一区域。在积分时只需要将替换在第一区域的积分即可得到的一般关系式。
在此区域依然恒等于零，的(因此区域，所以只是求其最大体积)一般关系式为

(15)
5.4.2 对区域2的讨论
由式(22)

可以得到

在第二区域须满足条件

计算可得

再由(12)式：
(16)
然后只需要做下式的变换即可得到一般关系式：

5.4.3 对区域3的讨论
同区域2中讨论，应该满足

得到：

由(13)式：
(16)
然后只需要做下式的变换即可得到一般关系式：

5.4.4 对区域4的讨论
应该满足关系：

5.4.5 对区域5的讨论
在区域5中油位高度始终保持为。

5.5 变位参数的确定和罐容表的计算
5.5.1用最小二乘参数估计法确定参数
最小二乘参数估计法基本思想：根据的关系表达式求得几组油量高度，计算出相邻高度油量的体积之差 通过与附件的实际储油量进行比较，通过对进行等间距的穷举最终求得理论值与实际值的差值的平方和，当取得最小值，此时的即为所求的最佳值。即求解如下最小二乘拟合模型

算法描述：
输入： 组显示油高。
输出： 纵横向偏角的值
Step1: 根据的高度值以及关系式，求得关于的表达式；
Step2: 根据附录找出实际对应的出油量；
Step3: 对进行等间距穷举，同时计算出，当取得最小值时，求得的值。
算法结束
用最小二乘参数估计法得到的变位参数为：，，角度都符合实际情况。
5.5.2 实际储油罐罐容表的制定
估计出实际储油罐的纵向倾斜变位参数和横向旋转变位参数后，我们就可以根据所建立的储油量和油位高度以及变位参数的一般模型计算得到罐容表。
表4 实际储油罐罐容表(纵向变位，横向偏转)
油位高度
0
0.1
0.2
0.2224
0.3
0.4
0.5
储油体积
00.0478
0.3558
1.0582
1.2789
2.2043
3.6764
5.4005
0.6
0.7
0.8
0.9
1
1.1
1.2
1.3
7.3345
9.4471
11.7128
14.1092
16.6165
19.2162
21.891
24.6243
1.4
1.5
1.6
1.7
1.8
1.9
2
2.1
27.4001
30.2029
33.0172
35.8279
38.6194
41.3765
44.0834
46.7237
2.2
2.3
2.4
2.5
2.6
2.7
2.8
2.9
49.2807
51.7364
54.0719
56.266
58.2952
60.1314
61.7392
63.0665
2.929
3

63.3862
64.003464.6645

5.6利用所给数据进行进一步的检验
我们根据推导的到的关系式对附件2中的第一次补充进油后的300组数据进行了验证，得到在相同油位高度差下，理论出油值值与实际出油值的关系图：

图15 理论储油量和实际储油量的差别对比图
从图可以看出所得的理论值与实际值基本重合，表示所算出的和的值是较为合理的。
六 模型的检验与评价
6.1 模型的检验（蒙特卡洛模拟方法）
在实际储油罐罐容表模型的建立和求解过程中，我们对球冠体内倾斜的部分燃油的体积进行了近似计算，忽略了一小部分球冠体体积。鉴于此，我们考虑通过用计算机模拟对该部分的体积进行模拟计算，观察近似计算值与精确模拟数值的吻合情况，同时也多我们建立的模型进行验证。
模拟过程中的主要步骤：
Step1: 划分空间，确定被忽略区域Q的空间限制范围，建立空间限制函数表达式。并寻求一包含该区域Q的最小长方体。建立坐标系，确定Q所在的区域范围；
Step2: 均匀做点，在长方体内分别从三个坐标轴依次等间距的产点，记录落入该区域的点以及生成的点的总数，计算该长方体区域的总体积；
Step3：统计落在该区域的点的个数，求该部分体积，计算公式为：。

在模拟中对不同区域分别进行求解所忽略部分的体积,再与所得到的理论值相加可得实际测量的精确值。根据蒙特卡洛模拟得到的被省略部分的体积，我们可以画出实际体积和近似体积的差别图。

图16蒙特卡洛模拟验证图

由图中可以看到，模型所推导出的结果与实际精确的数值吻合的非常好，模型对数据求解的精度很高，验证了模型的稳定性与可行性。
6.2 模型的评价
优点：
1、模型是由简单到复杂一步步建立的，增强可读性，逻辑性强；
2、在建模过程中将一些给定的数据参数化，使得模型更具一般性，可应用范围变广；
3、通过计算机模拟的方法较好地对近似计算所得结果做出了验证。
缺点：
在对附件2中实际储油罐的变位参数进行参数估计时，由于储油量和油位高度以及变位参数的函数关系非常复杂，用最小二乘参数估计实现很困难，我们仅采用了附件2中一些数据进行了参数估计，使得得到的变为参数与实际值可能相差较远，值得进一步改进的函数关系以使得更好地利用最小二乘参数估计法。
八 模型的改进和推广
7.1 模型的改进
建立模型积分出的函数关系很复杂，用最小二乘参数估计法进行参数估计时，程序运行时间较长。对此改进方案有：
1、将模型积分出的函数关系式通过一定的变换，近似化简函数关系，这样会减少一定的程序运行时间；
2、利用其它更为简便的参数估计方法进行求解。
7.2 模型的推广
对于不同于题中所给的油罐也可用此方法和理论进行分析罐内油位高度与储油量的对应关系，并且进行变位识别和罐容表的重新标定；我们用到的参数估计方法也可用到其他的实际生活中参数估计问题的求解。
参考文献
[1] 林雪松，周婧，林德新，MATLAB7.0应用锦集，北京：机械工业出版社，2006
[2] 王郑耀 ，卧式加油灌剩余油料体积的计算，西安交通大学 理学院，2004
[3] 王连群，李莉，石油油罐体积计算方法的探讨，吉林化工学院学报，第6期：45-50,1989
附录
附录一：部分结果
附表1 小椭圆储油罐罐容表(无变位)
油位高度
0
0.01
0.02
0.03
0.04
0.05
0.06
储油体积
0
0.0053
0.0149
0.0274
0.042
0.0586
0.0768
0.07
0.08
0.09
0.1
0.11
0.12
0.13
0.14
0.0966
0.1177
0.14
0.1636
0.1882
0.2139
0.2405
0.2681
0.15
0.16
0.17
0.18
0.19
0.2
0.21
0.22
0.2965
0.3258
0.3558
0.3866
0.4181
0.4503
0.4831
0.5165
0.23
0.24
0.25
0.26
0.27
0.28
0.29
0.3
0.5506
0.5852
0.6204
0.656
0.6922
0.7289
0.766
0.8035
0.31
0.32
0.33
0.34
0.35
0.36
0.37
0.38
0.8415
0.8799
0.9186
0.9578
0.9972
1.0371
1.0772
1.1176
0.39
0.4
0.41
0.42
0.43
0.44
0.45
0.46
1.1583
1.1993
1.2405
1.282
1.3237
1.3657
1.4078
1.4501
0.47
0.48
0.49
0.5
0.51
0.52
0.53
0.54
1.4926
1.5353
1.5781
1.621
1.6641
1.7072
1.7505
1.7938
0.55
0.56
0.57
0.58
0.59
0.6
0.61
0.62
1.8373
1.8808
1.9243
1.9679
2.0115
2.0551
2.0987
2.1423
0.63
0.64
0.65
0.66
0.67
0.68
0.69
0.7
2.1858
2.2294
2.2729
2.3163
2.3596
2.4029
2.4461
2.4891
0.71
0.72
0.73
0.74
0.75
0.76
0.77
0.78
2.5321
2.5749
2.6175
2.66
2.7023
2.7445
2.7864
2.8281
0.79
0.8
0.81
0.82
0.83
0.84
0.85
0.86
2.8696
2.9108
2.9518
2.9925
3.033
3.0731
3.1129
3.1524
0.87
0.88
0.89
0.9
0.91
0.92
0.93
0.94
3.1915
3.2303
3.2686
3.3066
3.3442
3.3813
3.4179
3.4541
0.95
0.96
0.97
0.98
0.99
1
1.01
1.02
3.4898
3.5249
3.5596
3.5936
3.627
3.6599
3.692
3.7235
1.03
1.04
1.05
1.06
1.07
1.08
1.09
1.1
3.7543
3.7844
3.8136
3.842
3.8696
3.8962
3.9219
3.9466
1.11
1.12
1.13
1.14
1.15
1.16
1.17
1.18
3.9701
3.9925
4.0136
4.0333
4.0515
4.0681
4.0828
4.0952
1.19
1.2

4.1049
4.1101

附表2 小椭圆储油罐罐容表(纵向变位)
油位高度
0.00
0.01
0.02
0.03
0.04
0.05
0.06
储油体积
00.0017
0.0035
0.0063
0.0100
0.0148
0.0207
0.0279
0.07
0.08
0.09
0.10
0.11
0.12
0.13
0.14
0.0363
0.0461
0.0574
0.0701
0.0844
0.1002
0.1177
0.1369
0.15
0.15
0.16
0.17
0.18
0.19
0.20
0.21
0.1513
0.1578
0.1802
0.2040
0.2289
0.2548
0.2818
0.3097
0.22
0.23
0.24
0.25
0.26
0.27
0.28
0.29
0.3385
0.3681
0.3985
0.4296
0.4614
0.4939
0.5271
0.5609
0.30
0.31
0.32
0.33
0.34
0.35
0.36
0.37
0.5952
0.6301
0.6655
0.7015
0.7379
0.7748
0.8121
0.8499
0.38
0.39
0.40
0.41
0.42
0.43
0.44
0.45
0.8881
0.9267
0.9656
1.0049
1.0445
1.0845
1.1247
1.1653
0.46
0.47
0.48
0.49
0.50
0.51
0.52
0.53
1.2061
1.2472
1.2885
1.3300
1.3718
1.4138
1.4559
1.4983
0.54
0.55
0.56
0.57
0.58
0.59
0.60
0.61
1.5408
1.5834
1.6262
1.6691
1.7121
1.7553
1.7985
1.8417
0.62
0.63
0.64
0.65
0.66
0.67
0.68
0.69
1.8851
1.9284
1.9719
2.0153
2.0588
2.1022
2.1456
2.1891
0.70
0.71
0.72
0.73
0.74
0.75
0.76
0.77
2.2324
2.2758
2.3190
2.3622
2.4053
2.4483
2.4912
2.5340
0.78
0.79
0.80
0.81
0.82
0.83
0.84
0.85
2.5766
2.6190
2.6614
2.7035
2.7454
2.7872
2.8287
2.8700
0.86
0.87
0.88
0.89
0.90
0.91
0.92
0.93
2.9110
2.9518
2.9923
3.0325
3.0724
3.1119
3.1512
3.1900
0.94
0.95
0.96
0.97
0.98
0.99
1.00
1.01
3.2285
3.2667
3.3044
3.3416
3.3785
3.4148
3.4507
3.4860
1.02
1.03
1.04
1.05
1.06
1.07
1.08
1.09
3.5208
3.5551
3.5887
3.6218
3.6541
3.6859
3.7169
3.7471
1.10
1.11
1.12
1.13
1.14
1.15
1.16
1.17
3.7766
3.8052
3.8330
3.8598
3.8856
3.9103
3.9338
3.9560
1.1713
1.18
1.19
1.20

3.9588
3.9766
3.9955
4.10174.1101

附录二：文中所涉及的程序
2.1求解小椭圆储油罐变位前的公式及罐容表
One1.m
syms y a b h L
v0=L*int(2*a*sqrt(1-(y-b)^2/b^2),y,0,h);
v=simple(v0)
pretty(v)

one2.m
h=0:0.01:1.2;
v=5/9*(-13083/10000*(18/25-3/5*h).^(1/2).*h+4361/2000*(18/25-3/5*h).^(1/2).*h.^2+39249/125000*15^(1/2)*asin(1/6*30^(1/2)*h.^(1/2)).*h.^(1/2))./h.^(1/2)*15^(1/2)
plot(h,v)
grid on
axis([0 1.2 0 4.5])
xlabel('油位高度/m')
ylabel('储油量/L')

one3.m 测量值和理论值数据对比
v1=[50:50:2050 2053.83 2103.83 2105.06 2155.06 2205.06 2255.06 2305.06 2355.06 2404.98 2406.83:50:2906.83 2906.91:50:3706.91];
v1=(262+v1)/1000;
h=[159.02 176.14 192.59 208.50 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.40 462.62 474.78 486.89 498.95 510.97 522.95 534.90 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.70
764.16 776.53 788.99 801.54 814.19 826.95 839.83 852.84 866.00 879.32
892.82 892.84 906.53 920.45 934.61 949.05 963.80 978.91 994.43 1010.43
1026.99 1044.25 1062.37 1081.59 1102.33 1125.32 1152.36 1193.49]';
h=h/1000;
v2=5/9*(-13083/10000*(18/25-3/5*h).^(1/2).*h+4361/2000*(18/25-3/5*h).^(1/2).*h.^2+39249/125000*15^(1/2)*asin(1/6*30^(1/2)*h.^(1/2)).*h.^(1/2))./h.^(1/2)*15^(1/2);
plot(h,v1,'r',h,v2,'k','markersize',3)
xlabel('油位高度h/m')
ylabel('储油量V/m3')
v1./v2

%变位后理论值和实际值相比较
v1=[
1.01 1.0583 1.118 1.1675 1.2221 1.2791 1.3274 1.3825 1.4335 1.4849 1.536 1.5918 1.6451 1.6973 1.7498 1.7966 1.8487 1.9 1.9527 1.9991 2.0526 2.1032 2.1484 2.2022 2.2527 2.3037 2.3521 2.4027 2.4522 2.4978 2.5486 2.5995 2.6446 2.6962 2.7488 2.7988 2.8472 2.8958 2.942 2.9925 3.0441 3.0896 3.1408 3.1897 3.234 3.2841 3.3345 3.379 3.4235 3.4731 3.5176 3.5694 3.5732];
v2=[0.9629 1.0129 1.0629 1.1129 1.1629 1.2129 1.2629 1.3128 1.3628 1.4127 1.4627 1.5127 1.5627 1.6127 1.6627 1.7127 1.7627 1.8127 1.8627 1.9127 1.9627 2.0127 2.0627 2.1127 2.1627 2.2127 2.2627 2.3127 2.3627 2.4127 2.4627 2.5127 2.5627 2.6127 2.6627 2.7127 2.7627 2.8127 2.8627 2.9127 2.9627 3.0127 3.0627 3.1127 3.1627 3.2127 3.2627 3.3127 3.3627 3.4127 3.4627 3.5127 3.5147 ];
s1=v1-v2;
s=(v1-v2)./v1;
h=[411.29 423.45 438.33 450.54 463.90 477.74 489.37 502.56 514.69 526.84 538.88 551.96 564.40 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.90 761.55 773.43 785.39 796.04 808.27 820.80 832.80 844.47 856.29 867.60 880.06 892.92 904.34 917.34 929.90 941.42 954.60 968.09 980.14 992.41 1006.34 1019.07 1034.24 1035.36];
h=h/1000;
plot(h,s)
grid on
xlabel('油位高度/m')
ylabel('储油量V/立方米')
figure(2)
plot(h,s1)
grid on
xlabel('油位高度/m')
ylabel('储油量/m³')
figure(3)
plot(h,v1,'k*',h,v2,'r')
grid on
xlabel('油位高度/m')
ylabel('储油量/m³')

clear all
syms y h
va=2*int(0.89*sqrt(1-y^2/0.6^2)*((h+0.4*0.0717-0.6+y)/0.0717),y,(0.6-(h/0.0717+0.4)*0.0717),0.6);
for h=0:0.0005:2.05*tan(4.1*pi/180)
r=eval(va);
hold on
plot(h,r)
grid on
end
axis([0 1.4 0 4.5])
syms h
v1=2*int(((h+0.4*0.0717)-y)/0.0717*(0.89*(1-((y-0.6)/0.6)^2)^0.5),y,h,(h+0.4*0.0717));
v=0.89*0.6*2.45*((h-0.6)/0.36*(h*(1.2-h))^0.5+asin(h/0.6-1)+pi/2);
v2=2*int((2.45-((h+0.4*0.0717)-y)/0.0717)*(0.89*(1-((y-0.6)/0.6)^2)^0.5),y,(h-2.05*0.0717),h);
v=v1+v-v2;
for h=2.05*tan(4.1*pi/180):0.005:1.2-0.4*tan(4.1*pi/180)
r=eval(v);
hold on
plot(h,r,'r.')
grid on
end
syms y h
va=4.1101-2*int(0.89*sqrt(1-y^2/0.6^2)*((1.2-h+2.05*0.0717-0.6+y)/0.0717),y,(0.6-((1.2-h)/0.0717+2.05)*0.0717),0.6);
for h=1.2-0.4*tan(4.1*pi/180):0.0001:1.2
r=eval(va);
hold on
plot(h,r,'g.')