储油罐变位识别与标注.docx

储油罐变位识别与标注
摘要
本文主要采用微积分为主体思想，建立出主体模型，在模型改进中尝试采用蒙特卡罗的思想进行进一步探索。求解时，模型全部封装成程序实现，方便了应用与推广，而且误差验证表明两个模型都达到了很高的精度。
对于问题一：
1、对原数据分析，对无变位时的数据进行分析，研究对应高度的理论容积与实际装油的量之间的关系，其差值存在说明可能是罐体内装置及替他因素造成，最终用回归拟合的方法得出一个三次函数，通过作图发现发现吻合程度很高。
2、采用微积分的思想计算倾斜后应高度的容积，得出容积与高度的函数。
3、结合前两个步骤得出最终的变位后标注于高度的函数关系。采用MATLAB求解出新的标定表，摘部分结果：
高度\cm
25
50
75
100
标注\L
427.214
1325.557
2360.045
3322.853
4、误差分析，用自己的模型函数算出题中给出数据的高度对应标注值，与实际油量进行对比，得出平均绝对值误差率为，模型达到了很高的精度，有力的证明了模型的正确性。
对于问题二：
1、用积分的方法算出对应油浮子标度下的理论容积。在求解的过程中，采用分块积分的细想，水平分上中下，竖直分左中右，然后结合结合具体的图形算出相关的关系，进行积分的出最后的函数：，然后将其此模型的求解与建立封装成MATLAB程序，方便下面的校正与求解。
2、把上一步的体积模型函数两个变位参数，取0度，此时把得出对应高度的结果与题目中给的各高度对应容量进行对比，这里均匀选取23个代表点计算，差值基本恒定,其均值为，然后加到第一步的体积函数中作校正，得出最终的模型关系，然后修正模型对应的完整封装的MATLAB程序。得出的最终表达式进行误差分析：，可见本模型已经达到了极佳的精度。
3、根据总容积恒定的参照，即为累计的出油加上此时刻的标注值应该等于一个定值，考虑到角度不会太大，所以这里采用计算机遍历搜索的思想，以这23个代表点算出的每个不同参数总容积的方差最小为目标进行搜索，编程序直接调用上面的模型对应封装好的MATLAB源程序，在0到15度之间以0.1为步长遍历搜索出两个变位参数的值： ，
4、将求出的参数带入模型，及出每个10cm的新的标注表，部分结果摘录：
高度/cm
50
100
150
200
标注/L
5576.665
16953.22
30622.84
44426.59
5、在模型改进中，我们又进一步尝试了蒙特卡罗算法，利用约束范围内的随机点占一个容易求解的规则图形随机点的个数之比，得出两者之间体积的对应关系，进而求出容积与高度及变位关系的函数。我们对比结果后发现误差比积分模型大一些，但此模型在形状等更复杂的场合还是更有适应力的。
关键字 ： 微积分 回归拟合 计算机搜索 蒙特卡罗算法
1问题重述
通常加油站都有若干个储存燃油的地下储油罐，并且一般都有与之配套的“油位计量管理系统”，采用流量计和油位计来测量进/出油量与罐内油位高度等数据，通过预先标定的罐容表（即罐内油位高度与储油量的对应关系）进行实时计算，以得到罐内油位高度和储油量的变化情况。
许多储油罐在使用一段时间后，由于地基变形等原因，使罐体的位置会发生纵向倾斜和横向偏转等变化（以下称为变位），从而导致罐容表发生改变。按照有关规定，需要定期对罐容表进行重新标定。图1是一种典型的储油罐尺寸及形状示意图，其主体为圆柱体，两端为球冠体。图2是其罐体纵向倾斜变位的示意图，图3是罐体横向偏转变位的截面示意图。
请你们用数学建模方法研究解决储油罐的变位识别与罐容表标定的问题。
（1）为了掌握罐体变位后对罐容表的影响，利用如图4的小椭圆型储油罐（两端平头的椭圆柱体），分别对罐体无变位和倾斜角为a=4.10的纵向变位两种情况做了实验，实验数据如附件1所示。请建立数学模型研究罐体变位后对罐容表的影响，并给出罐体变位后油位高度间隔为1cm的罐容表标定值。
（2）对于图1所示的实际储油罐，试建立罐体变位后标定罐容表的数学模型，即罐内储油量与油位高度及变位参数（纵向倾斜角度a和横向偏转角度b ）之间的一般关系。请利用罐体变位后在进/出油过程中的实际检测数据（附件2），根据你们所建立的数学模型确定变位参数，并给出罐体变位后油位高度间隔为10cm的罐容表标定值。进一步利用附件2中的实际检测数据来分析检验你们模型的正确性与方法的可靠性。

2模型假设
1） 假设温度的变化不会引起较大的容积变化；
2） 假设罐体壁使用一段时间不会发生较大形变；
3） 假设探针的刻度显示受外界影响因素不大；
4） 假设罐体内不会积累固体杂物；

3符号说明
符号
符号意义
变量类型

油罐液面高度
问题一中
变量

油罐液高为的储油容积

油罐纵向倾斜角度

油罐内一些附件所占的体积

罐体长度
问题二中
变量

罐体油浮子左边部分长度

罐体油浮子右边部分长度

标注高度

柱体高度

正面两端的球体部分弦高
其他的函数等变量在模型求解中给予说明。
4问题分析
4.1 问题一：
首先，需要算出对应高度液面以下，罐体内设施部分所占的体积，最终将其拟合成一个关于高度的函数，为修正标注做准备。考虑到设施都是固定的柱体且横截面积并不是很大，故这里采用无变位放置时的数据进行计算时，在修正标注时可以直接使用。这里分两步进行：
1、先用积分的方法算出对应注油高度时的罐体容积；
2、通过对应高度的罐体容积减去此高度对应的油的体积得出罐体内装置等造成的差值体积；
3、拟合出罐体装置等其他体积随高度的变化函数
然后，进行第一次标注修正，分为两步骤进行：
1、用积分的方法算出变位后每个高度对应容积；
2、用变位后对应各高度的容积减去此高度的装置等造成的差值体积，然后得出变位后的标注。
4.2 问题二：
对于此问变位后的标注修正，然后计算变位后对应油浮子高度液面以下的容积，而且这是一个规则形状的物体，用微积分的方法最终可以得到一个精度很高的标注于不同变位角度参数和高度的关系模型。
对于积分我们采用分块讨论的细想，将球罐体分为三步分进行讨论，左边的球体部分，中间的柱体部分，右边的球体部分，垂直方向分成上中下三个部分，从倾斜后中间部分非常容易计算，两端部分很小且相等，所以采用整体减去中部再除以二得出。
对于横向变位参数考虑，由于油罐两端的顶板为圆缺，所以油罐分别沿横直径和纵直径呈轴对称，由于前面已经对纵向倾斜变位进行了较为深入的研究。所以在此对横行偏转倾斜进行研究：

图1横向倾斜截面图
油罐如果发生横向变位（如图1(b)所示），油罐绕其中心轴发生旋转，油罐的空间位置发生变化，但是由于油罐体的顶板为球冠体，呈轴对称，所以罐内的液体相对于水平面是没有变化的，即实际液面高度是不变化的，但横向变位后，浮子的读数发生偏差。
我们考虑油罐体未发生变位的状态，通过改变油位探针的位置模拟横向变位所带来的误差。即若发生横向变位转动角为，则我们使未变位状态下油位探针在横截面内沿转动的反方向倾斜角（如图6(c)所示）
然后结合微积分可以得出左边球体，中间主题，右边球体，随油浮子显示高度的函数，对于上中下的三部分，算出中间部分的总体后，然后用整个罐体的体积减去其得出上下端的，然后从下端与中部的分界面出发，下端部分为常数，加从此分界面以上到制定高度的容积，即可算出每个高度对应的总容积。
得出变位后油浮子标注于高度，两个变位参数的关系后，将零度带入两个参数与实际数据进行对比验证，修正函数，然后算出误差，检验模型精度。
4.3模型优化
考虑到用积分的方法计算溶液的体积需要考虑的一度比较复杂，而其计算困难，所以本文又给出了用蒙特卡洛方法计算油罐内储油容积的方法。运用蒙特卡洛方法的思想，将计算容积的过程转化为随机点落在液面与油罐体在坐标系中围成的区域中的概率。
5 模型建立与求解
经过以上的分析与准备，我们将逐步建立以下的数学模型，进一步阐述模型的实际建立过程。
问题一：
1计算无变位储油容积
首先研究油罐罐体横截面，建立如
图2坐标系

图2坐标系中油罐横截面
横截面椭圆的方程为：
液面高度为时，，点坐标为，所以，其中，容易求出
(1.1)
其中为油罐内液高。
在无变位情况下油罐液面高度为时的容积即可用下面定积分积分公式计算得到：
(1.2)
其中为椭圆筒的长度。
在无变位情况下油罐液面高度为时的容积即可用下面定积分积分公式计算得到：
(1.3)
其中为椭圆筒的长度。
2 变位后椭圆筒的部分容积计算：

图3纵向倾斜罐体剖面图
1． 为水平线；
2． 为罐中油面；
3． 罐体向左下方倾斜角度；
4． 为油位探针，为浮子的位置；
5． 为从油罐右边底部发出的水平线；
的长度反映了油位探针在直圆筒上的位置，为油位探针的长度，的长度是直圆筒的径向长度，这些变量可以通过测量得到。的长度是可以直接读出来（油位探针读出的液面高度）。
由于，所以，进而可以计算，
(1.4)
在中，
(1.5)

图4油罐液面的水平剖面
注：为罐体中液面，由于油罐倾斜角度较小故计算时将其近似为梯形。
此时储油液面为等腰梯形，且横截面为弓形。记端和端液面高度分别为。由式可以分别计算出端和端的页面宽度，分别记为。
研究油罐的部分容积，即计算液面下罐体的容积。
将液面下的容积分两部分分别计算：一部分为水平面下的罐体容积，另一部分为液面和水平面之间的罐体容积。
1） 液面和水平面之间的罐体容积
(1.6)
2） 水平面下的罐体容积
(1.7)

3模型的求解
用MATLAB解出每一个题中所给数据高度时的装置体积，为拟合做准备。
下面来考虑罐体内管道等设施所占体积对高度的函数，以毫米为单位，由于开始部分与结束部分非线性化程度较高，经过筛选采用三次拟合，拟合函数为

拟合函数曲线与原始数据拟合如下图：

图 5装置体积拟合对比曲线图
由此可见拟合的精度还是相当良好的。
综上，可以的出装置体积随装置的变化函数：

下面计算变位后的的每个标注高度对应水平面以下的容积，这将容器分为三部分，上端部分，中间部分，下端部分，中间部分截面近似为一个梯形，以此为基准来确定这三个部分的分界面。其中上端与下端部分所占的空间相等。
对于标注，只有中间部分实际意义最大，上端下端此时的游标尺寸作图即可以很容易计算出来上端下端对应的油浮子位置显示，题中要求以每厘米为间隔作为标注，所以这里采用厘米为此问求解的统一单位：

中间部分采用横向切片微元，由于倾斜角度很小，所以界面可以近似等效梯形。
利用此思想结合图形的角度变换，整理的到如下中间部分积分公式：

用总容积减去中间部分容积得出下部的容积：

对应高度的容积为：

用变位后的各个高度对应的容积减去此高度的容器内其他装置的到的修正后标注，标注从下部与中部的分界面开始，到上面分界面结束，及对应的油浮子的刻度应该从到，对应公式为:

然后，用MATLAB编程算出每1cm的变位后修正的标注值，变位后的结果见附录。

问题二：
由于油罐两端的顶板为圆缺，所以油罐分别沿横直径和纵直径呈轴对称，由于前面已经对纵向倾斜变位进行了较为深入的研究。所以在此对横行偏转倾斜进行研究：

图6横向倾斜截面图
油罐如果发生横向变位（如图6(b)所示），油罐绕其中心轴发生旋转，油罐的空间位置发生变化，但是由于油罐体的顶板为球冠体，呈轴对称，所以罐内的液体相对于水平面是没有变化的，即实际液高是不变化的，但横向变位后，浮子的读数发生误差。
我们考虑油罐体未发生变位的状态，通过改变油位探针的位置模拟横向变位所带来的误差。即若发生横向变位转动角为，则我们使未变位状态下油位探针在横截面内沿转动的反方向倾斜角（如图6(c)所示）
可以看出，用图6(c)的方法模拟油罐横向变位时，浮子的读数值与油罐发生偏转时（图6(b)）的读数是不同的。分别记图6(b)与图6(c)中浮子的读数为。

在图6(b)中，为截面圆心，均为直径，截面圆的半径为,，
在中，，则

在中，，所以，
(1.8)
1 计算罐体的总体积
中间部分主体体积：

两端球体的体积，采用积分的进行求解：

总体积为：

计算时统一采用分米为单位，结果为：
dm
2 对横向变位参数进行折算
由于是圆筒形的光体，随意横向旋转可以等效为探针的的变位，所以经过几何分析，当探针上的油浮子显示为时，等效为没有横向变位的高度为：

再推到出此函数的反函数，用来确定上中下三层分界面时的对应油浮子高度所显示的值：

下层分界面浮子对应显示高度为：

上层分界面浮子对应显示高度为；

然后将此函数编成MATLAB函数文件代码，为以下计算做铺垫。
3 计算中间柱体部分的中间层对应高度函数
中间柱体部分的函数，以横向切片为微元，切片近似为梯形
侧面圆的方程，取中心为原点建坐标：

带入上一步折算的，可以算出此时对应的近似梯形截面积的两个上下底，
最终整理的出积分函数为：

再将此文件编程MATLAB函数文件，为下面计算铺垫。
4 计算左端球体和右端球体中间层对应油浮子高度的容积函数
在计算左右两球端的液面高度时，我们做了一个等效，在液面的高度作出一个垂直于此球体与柱体结合面的平面，两端以此平面一下的体积作为此高度时两端球体所装的液面，等效前后差别很小，且两端成互补的作用，所以这样的精度还是不会有明显的降低的。
同样这里先对面微元，然后积分求出体积即可：
先用勾股定理列方程算出两端球体所对应的半径
dm
然后得出积分函数与对应高度的函数表达式：

其中有：

然后可以得出左右球体部分中层所高度变化的函数：
左端球为：

右端球为；

再将此文件编程MATLAB函数文件，为下面计算铺垫。

5 计算出下层的容积
先用前几部算出的式子得出中部总体积：

然后可以的下部的体积：

6 得出标注与高度即变为关系的参数
这里求出下层体积，加上中层的关于高度的表达式，即为对应高度的最后标注，式子如下：

由于最后展开式子过长，表述不清楚，而且也不方便程序的实现，所以本文不与展开最后的表达式，只是将其总结如下：

然后总结所有前面自己的MATLAB文件函数，将整个模型程序化，代码见附录。
7 将变位参数置零度带入模型与实际参数对比改进得出最终模型
这里选取一次注油前的数据中，每隔1dm晒出一个点，取出23个点进行对比校正，用上一步的模型的出这些点对应高度的计算容积后，在Excel中比较了它们之间的差值，发现基本为一个固定的数值。
这些点的对应油浮子高度为：
这些点的对应题中标注量为：
这些点带入模型求出的标注量：

差值为：

观察发现，基本为一个恒定值，求助其平均值为：

其值为：
故最后修正后模型确立如下：

8 利用累积出油加罐体油总体积不变的参照遍历搜索变位参数
首先在Excel中算出每个高度时的累积出油量，如果没有倾斜的画那么标度所示的容量加上累积的出油量应该是一个常值，我们取第一次加油前的数据做了一个累积出油量与此高度所对应的标注值的图像如下：

表1题中数据油量与累计出油和的变化规律
所以图形反映出原来的标注已经不够准确，这里必须筛选出更加准确的变位参数。
这里对于变位参数的确定，我们采用遍历求解法，以筛选出来的数据的累计出油量与此高度的标注量之和方差最小为目标进行变为参数的搜索。
搜索的MATLAB程序见附录，结果为

9 利用最终的模型计算最终修正变位后每隔10厘米的标注
带入上个步骤的求出变位参数值，带入模型求出每个10cm的变位后的修正标注。

带入前面编号的程序即可得出新的标注，结果为：

浮子对应高度/cm
变位后标注/L
浮子对应高度/cm
变位后标注/L
10
659.5428
160
33435.88
20
1392.841
170
36238.41
30
2386.95
180
39014.87
40
3838.127
190
41749.58
50
5576.665
200
44426.59
60
7538.602
210
47029.48
70
9684.613
220
49541.18
80
11984.98
230
51943.64
90
14415.03
240
54217.46
100
16953.22
250
56341.2
110
19579.99
260
58290.29
120
22277.18
270
60034.76
130
25027.63
280
61533.69
140
27814.86
290
62709.78
150
30622.84
300
63449.31

6模型改进与推广
1 蒙特卡罗方法的原理

蒙特卡罗方法又被称为统计模拟法、随机抽样法，这种方法是通过抓住事物运动的几何数量和几何特征，将其同一定的概率模型联系起来，用电子计算机实现统计模拟或抽样，以获得问题的近似解。
蒙特卡罗方法的基本原理及思想如下：当所求问题的解是某个事件的概率，或者是某个随机变量的数学期望，或者是与概率、数学期望有关的量时，通过某种试验的方法，得出该事件发生的频率，或者该随机变量若干个具体观察值的算术平均值，通过它得到问题的解。
对于本题，在计算油面高度为的油罐部分容积时涉及到多重积分。建立适当的三维坐标系，计算出卧式油罐表面的方程和油罐内液面的方程，记油罐表面与液面围成的封闭区域为。在三维坐标系中一个体积已知的区域()内进行均匀随机的取点，同时判断该点是否在区域中。当抽样点足够多时，样本点落在油罐表面和油面所围成的封闭区域内的个数与的比值与区域和区域的体积比是一致的，即，从而可得到区域体积的无偏估计。具体操作中，用计算机产生均匀分布在区域内的随机数序列进行模拟。
2 利用蒙特卡罗方法求油罐部分容积新模型
变位前,直椭圆筒处于水平状态,此时可以以油罐体的中心为原点，以直椭圆筒体的截面椭圆长轴所在直线为坐标轴，截面椭圆短轴所在直线为轴，以直椭圆筒的中心轴所在直线为轴，建立三维直角坐标系。
此时容易获得：
直椭圆筒的方程为：

平面方程为：

其中为水平状态下直椭圆筒内液面的高
则油罐表面与液面围成的封闭区域为可表示为：

对于变位后的油罐体，以油罐体中心为原点，仍以直椭圆筒体的截面椭圆长轴所在直线为坐标轴，截面椭圆短轴所在直线为轴，以直椭圆筒的中心轴所在直线为轴，建立三维直角坐标系。
此时，液面在此坐标系中的剖面图如下：

图7变位后新坐标系下油面剖面图
我们很容易在上述坐标系中给出直椭圆筒的方程：

根据本题中数据，则有

由变位后浮子的读数已知，则椭圆纵截面中

在上述所建坐标系中可计算得到液面上点

显然，液面所在的平面平行于坐标轴轴，所以可以取向量，可知
不平行，则可以求出液面所在平面的一个法向量：

再有已知点的坐标，且根据题目中实际数据，则可求出平面一般方程为
则油罐表面与液面围成的封闭区域为可表示为：
(1.9)

取一个长方体区域，使得个椭圆顶卧式油罐包含在中。然后通过计算机程序【见附录】产生均与分布在区域中的随机数列进行模拟。并统计出，取样点的个数和落在区域中的取样点个数，根据蒙特卡罗原理，我们可以得到液高为油罐中油的体积：

7 模型检验
问题一：油罐标定模型的误差分析
利用题目所给的数据对此模型，进行误差检验分析：
带入每一个题中所给高度，用MATLAB对应算出相应的值，然后与实际的值进行比较，作出比较图如下：

图2 模型算出数值与题目数据对比图
由于题目数据太多，说以这里选取十组数据的对比表格如下：

高度
计算值
实际值
误差
680.63
2071.599925
2062.73
8.86992526
693.03
2123.357344
2112.73
10.62734376
704.67
2171.890086
2162.73
9.160085555
716.45
2220.934854
2212.73
8.204854144
727.66
2267.521361
2262.73
4.791361112
739.39
2316.160876
2312.73
3.430875812
750.90
2363.761758
2362.73
1.031757963
761.55
2407.677855
2412.73
-5.052144666
773.43
2456.501431
2462.73
-6.228569189
785.39
2505.457576
2512.73
-7.272423632
表1 一组数据对比
下面计算此模型的平均误差率：
对所有的数据误差先取绝对值，在求平直：

求出结果为：

由此得误差满足要求。

问题二：对变位参数为零时模型值做误差分析检测精度
在这里仍然用模型建立与求解所选取的23个点，非常具有代表性，做出数据如下
显示油高/dm
显示油量/L
计算油量
相对的差值
差值的绝对值
26.32
60448.88
60447.88
1.00
1.00
25.28
58534.01
58536.79
-2.78
2.78
24.32
56572.86
56585.36
-12.50
12.50
23.28
54297.75
54297.76
-0.01
0.01
22.33
52066.65
52072.51
-5.86
5.86
21.38
49744.05
49736.48
7.57
7.57
20.36
47119.35
47123.9
-4.55
4.55
19.33
44382.33
44394.99
-12.66
12.66
18.31
41615.55
41621
-5.45
5.45
17.29
38778.25
38793.24
-14.99
14.99
16.36
36192.94
36182.4
10.54
10.54
15.26
33077.38
33071.86
5.52
5.52
14.32
30416.13
30409.15
6.98
6.98
13.38
27766.55
27755.69
10.86
10.86
12.31
24764.66
24762.87
1.79
1.79
11.31
22006.33
22008.71
-2.38
2.38
10.95
21037.34
21030.52
6.82
6.82
9.94
18338.22
18332.72
5.50
5.50
8.90
15648.56
15641.08
7.48
7.48
7.95
13280.49
13275.62
4.87
4.87
6.92
10830.14
10830.99
-0.85
085
5.98
8714.85
8727.536
-12.69
12.69
利用这23组数据得出平均误差率为：

7模型评价
1优点
1) 利用积分方法求液高为的油罐储油容积的思路比较清晰，结果精确。
2) 所有模型都matlab程序实现，方便模型的应用于推广。
3) 以积分为主导的模型的精确度很高。
4) 蒙特卡罗方法有很强的适应性，问题的几何形状的复杂性对它的影响不大。该方法的收敛性是指概率意义下的收敛，因此问题维数的增加不会影响它的收敛速度，而且存贮单元也很省。

2缺点
1) 积分模型思路虽然比较清晰，容易想到，但是具体积分过程比较复杂，受维度和具体几何形状的限制；
2) 得出积分模型后，具体求解计算过程比较复杂；
3) 蒙特卡罗方法方法主要弱点是收敛速度较慢和误差的概率性质，其概率误差正比于，如果单纯以增大抽样粒子个数N来减小误差，就要增加很大的计算量。

8参考文献：
[1] 姜启源 邢文训，谢金星，杨顶辉，大学数学实验，北京：清华大学出版社，2005；
[2] 薛定宇 陈阳泉，高等应用数学问题的MATLAB求解，北京：清华大学出版社，2004；
[3] 韩中庚 数学建模方法及其应用 高等教育出版社 2005
[4] 姜启源 数学模型 高教出版社 2002
[5]董辰挥， 彭雪峰， MATLAB2008全程指南 北京：电子工业出版社 2009

附录
附录一：第一问的变位后的罐容标定表
标号\cm
修正标注\L
标号\cm
修正标注\L
15
181.5694014
66
1985.424191
16
199.8288983
67
2027.201679
17
220.2072046
68
2068.969184
18
242.1845401
69
2110.714649
19
265.51268
70
2152.42599
20
290.0328182
71
2194.091088
21
315.629865
72
2235.697774
22
342.214217
73
2277.233816
23
369.7127743
74
2318.68691
24
398.0639056
75
2360.044663
25
427.2143663
76
2401.29458
26
457.1172858
77
2442.424051
27
487.7307867
78
2483.420336
28
519.0169982
79
2524.270548
29
550.9413287
80
2564.961637
30
583.4719157
81
2605.480373
31
616.5791984
82
2645.813326
32
650.2355817
83
2685.946851
33
684.4151646
84
2725.867057
34
719.0935207
85
2765.559795
35
754.2475154
86
2805.010624
36
789.8551543
87
2844.20479
37
825.8954546
88
2883.127192
38
862.3483361
89
2921.762353
39
899.1945271
90
2960.094379
40
936.415483
91
2998.106925
41
973.9933158
92
3035.783147
42
1011.910732
93
3073.105652
43
1050.150976
94
3110.056444
44
1088.697783
95
3146.616858
45
1127.535337
96
3182.767488
46
1166.648225
97
3218.488106
47
1206.021411
98
3253.757563
48
1245.640194
99
3288.553678
49
1285.490188
100
3322.853103
50
1325.557289
101
3356.631172
51
1365.827651
102
3389.861709
52
1406.287665
103
3422.516799
53
1446.92394
104
3454.566512
54
1487.723276
105
3485.978547
55
1528.672655
106
3516.717789
56
1569.759215
107
3546.745733
57
1610.970239
108
3576.019721
58
1652.293138
109
3604.491892
59
1693.715437
110
3632.107721
60
1735.224758
111
3658.803857
61
1776.80881
112
3684.504784
62
1818.455374
113
3709.117302
63
1860.152291
114
3732.520479
64
1901.887447
115
3754.544682
65
1943.648766
116

3774.916186

117
3792.995826

附录二：第一问平稳放置时对应注油高度的罐容积
load h; %题目所给注油的对应高度值
fun=@(x)(sqrt(1-x.^2/60^2)*2*89*245);
rongji=zeros(78,1);
for k=1:78
rongji(k)=quad(fun,-60,(-60+h(k)/10))/1000-262;
end

附录三：第一问罐体内各装置所占体积随高度变化的函数拟合及图形对比
load h;
load cha; % 每个高度对应的罐容积与油量的体积差
p1=polyfit(h,cha,3);
x1=1:1200;
y1=polyval(p1,x1);
plot(h,cha,'o')
hold on
plot（x1,y1);

附录四：第一问变位后的油浮对应水平面以下的容积
fun=@(x)(sqrt(1-x.^2/60^2)*2*89);
s=quad(fun,(-60+17.5618),60);
vDown=(pi*60*89-s)*245*0.5/1000; %油浮下面的14.7cm对应的容积
%------------------------------
biaozhu=zeros(103,1);
hbiao=[15:117]'; %标注从15cm到117cm
for k=1:103
biaozhu(k)=quad(@fun,14.69458,hbiao(k))/1000+ vDown；%
end
final=real(biaozhu);
%--------------------------------
function y=fun(x) %积分函数
a=4.1/180*pi;
high=245/cos(a);
shangdi=x+tan(a)*40-60;
xiadi=x-tan(a)*205-60;
y=(fun1(shangdi)+fun1(xiadi))*high*0.5;
end
function y = fun1(x) %此函数的子函数
y=sqrt(1-x.^2/60^2)*2*89;
end

附录五：第一问修正变位后的标注和实际数据图形对比
load ygao %第一列是题中数据高度，第二列是油量
biaozhu=zeros(53,1);
hbiao=ygao(:,1)/10;
for k=1:53
biaozhu(k)=quad(@fun,14.69458,hbiao(k))/1000+186.5114;
end
biaozhu=real(biaozhu);
load p1
zhuangzhi=polyval(p1, ygao(:,1));
biaozhu=biaozhu-zhuangahi;
shiji=ygao(:,2);
plot(ygao(:,1),biaozhu);
hold on
plot(ygao(:,1),shiji,’o’);

附录六：第一问变位后修正的每隔1cm的标注
clc
biao=zeros(103,1);
h=15:117；
for k=1:103
biao(k)=quad(@fun,14.69458,h(k))/1000+186.5114;
end
biao=real(biao);
load p1
zhuang=polyval(p1, h);
biao=biao-zhuang;
xlswrite('biao.xlsx',biao) %将结果输出到Excel中

附录七： 第二问模型的求解导出及对参数封装的模型源程序
global L1 L2 A B D L
L1=20;
L2=60;
D=30;
L=80;
format long g
jiaoA=input(' 输入纵向变位参数');
jiaoB=input('输入横向变位参数');
disp('hBiao向量为修正的标注值');
A= jiaoA/180*pi;
B=jiaoB/180*pi;
hDown=chu(L2*tan(A));
hUp=chu(D-L1*tan(A));
vTotal=pi*(D/2)^2*L+3500/3*pi;
vCLeftBall=quad(@vleft,hDown,hUp);
vCRightBall=quad(@vleft,hDown,hUp);
vCMiddleZhu=quad(@middle,hDown,hUp);
vCMiddleZhu=real(vCMiddleZhu);
vCDown=(vTotal-vCMiddleZhu-vCLeftBall-vCRightBall)/2;
hBiao=zeros23,1);
load hyuan
for k=1:23
vLeftBall=quad(@vleft,hDown,hyuan(k)); vRightBall=quad(@vleft,hDown,hyuan(k));
vMiddleZhu=quad(@middle,hDown,hyuan(k));
vMiddleZhu=real(vMiddleZhu);
hBiao(k)=vCDown+vLeftBall+vRightBall+vMiddleZhu;
end

%-----------------------------
function y=ball(x)
global L1 A
hh=L1*tan(A)+x;
y=8*sin(A)/(3*pi)*sqrt(hh/30.*(1-hh/30));
end

%-----------------------------
function y=chu(x)
global B
y=(x-30*(1-cos(B)))/cos(B);
end

%-------------------------------
function y=di(x)
global A
y=sqrt(225-(x-15).^2);
end

%------------------------------