储油罐的变位识别.doc

储油罐的变位识别与罐容表标定模型
摘要
石油作为工业中的重要原料，加油站行业同样维系着国民经济命脉，因此加油站的地下储油罐的储油量问题也越来越引起人们的关注与重视，而由储油罐罐底引起的计量误差是目前最普遍也是从事油品储存加工企业最需解决的问题。

本文就针对罐体发生变位后对罐容表的影响进行研究。

通过对题目中提出的储油罐是否变位，以及变位之后罐内储油量与油位高度的关系和罐容表如何标定的问题进行分析，忽略实际操作过程中不可控的因素，通过对罐内油的横截面面积求积分的方法，建立了求罐内油量容积的模型，然后根据实际采集的数据来验证与矫正模型，但在积分的过程中，)(arcsin )(1)(2x g x g x g +-积分非常困难，为了简化计算，用契比切夫多项式逼近法将其近似处理为：
()
)(48)(80)(61515751653x g x g x g --π. 在问题一中针对小椭圆油罐的情况建立了变位前后储油量与油位高度关系的模型，利用给定的无变位与变位后的油位高度数据代入模型，在MATLAB 中编程计算得到了罐体内油容量的数据，比较发现变位前后计算得到的理论值与实际的油容量值之间的最大相对误差分别为：3.4917%，5.4695% ，说明了模型的准确性较高，然后根据模型获得了变位后油位高度间隔1cm 的罐容表标定值。

对于问题二，针对实际储油罐，罐体发生纵向倾斜α角度及横向偏转β角度的变位，选择适当的坐标系，建立了变位情况下的模型，通过对题目中给出的显示油量容积与出油量的数据的比较分析，发现显示的油量容积不是罐体内油量的准确值。

通过建立一个无变位情况下的新模型检验了模型的正确性较高，设定,101︒≤≤︒α,101︒≤≤︒β以︒1为间隔得到100组变位参数，再将油位高度代入模型求解出对应的100组理论的容油量准确值, 将求解出的容油量相邻之间的差值与对应出油量进行比较。

在不考虑异常点的情况下，当这个比值越小，说明给定的位变参数越准确，否则反之。

在这100组数据中就得到比值最小的一组()︒︒3,2，将其定为初始变位参数，然后采用逐步调试算法求得最终变位参数为︒=︒=2.3,9.1βα。

将最优变位参数代回模型求出最大相对误差为4.0613%，更加说明了模型的准确性。

最后将这组变位参数代入模型中得到罐体变位后油位高度间隔10cm 的罐容表标定值。

关键词：变位 积分 契比切夫多项式逼近法 罐容表标定
一、问题提出
加油站通常都有若干个与之配套的“油位计量管理系统”的储存燃油的地下储油罐，采用流量计和油位计来截量进/出油量与罐内油位高度等数据，通过预先标定的罐容表（即罐内油位高度与储油量的对应关系）进行实时计算，以得到罐内油位高度和储油量的变化情况。

但许多储油罐在使用一段时间后，由于地基变形等原因，使罐体的位置会发生纵向倾斜和横向偏转等变化，从而导致罐容表发生改变。

按照有关规定，需要定期对罐容表进行重新标定。

要求通过对罐体的纵，横截面的分析用数学建模方法研究解决储油罐的变位识别与罐容表标定的问题。

现根据题目给出的图形与附表解决以下问题：
1.根据对小椭圆型储油罐（两端平头的椭圆柱体）变位前后进/出油过程中采集的实验数据，建立数学模型研究罐体变位后对罐容表的影响，同时给出变位后油位高度间隔为1cm的罐容表标定值。

2.对实际储油罐，建立储油量与油位高度及变位参数（纵向倾斜角度α和横向偏转角度β）之间的一般关系，再根据附表数据来检验模型和确定变位参数，并给出罐体变位后油位高度间隔为10cm的罐容表标定值。

二、问题分析
为了掌握储油罐变位后对罐容表的影响，就需要讨论变位前后罐内油位高度与储油量的变化情况。

对问题一，小椭圆型储油罐罐体无变位时，罐内储油的横截面积是椭圆的一部分，可以计算出面积，且面积还是关于油位高度的函数，而计算储油量只需对侧面面积进行积分即可。

但当小椭圆型储油罐罐体发生倾斜角为︒1.4的纵向变位时，此时计算储油量要分三种情况考虑，第一种是油位高度在0的位置时，容积
215，因此油位高度在0以下的部分可以不讨论；=，而罐内容量初值为L
L
.1
v6744
第二种是油位高度在0到m
2.1间，储油量还是通过对罐内储蓄的油的横侧面积进行积分求得，只是此时油位高度是关于位置变化的函数；第三种是油位高度高于m
2.1时，油位高度不会再随着进油量的增加而变化，这种情况也不在谈论的范围内。

建立油位高度与储油量的变化关系模型后，再根据题目给出附件1的实验采集数据表，将表中油位高度带入模型后得出的储油量与表中的储油量求误差，通过误差值来检验模型的优良。

确定了模型就易得罐体变位后油位高度间隔1cm的罐容表标定值。

问题二中的实际储油罐是由圆柱体和两个球冠体构成。

要建立罐内储油量与油位高度及变位参数（纵向倾斜角度α和横向偏转角度β）之间的一般关系式，首先仅建立纵向倾斜α角度时的罐内储油量的模型，其圆柱体部分储油量的模型与问题一中变位后的模型相同，而球冠体部分的储油量通过对罐内储蓄的油的横侧面积求积分获得，横截面积的计算以罐体半高度所在的面为分界面，油位高度在分界面以下时，横截面是圆的一部分；而在分界面以上时，横截面积有整个圆和圆的一部分两种，对两部分积分再求和即可。

再建立仅横向偏转β角度的模型，横向偏转只改变了油位高度的真实值，因β非常小，我们可以近似看成油位实际高度与油位截量高度之间成三角函数关系。

当罐体横向和纵向同时发生变位时，只需把横向偏转后的真实油位高度代
入纵向变位罐内储油量的模型就可获得横纵向同时变位时罐内储油量的模型。

α,的一般关系式，先此时建立的罐内储油量的模型是关于油位高度及变位参数()β
要通过题目给出的附件2的实际采集数据来检验我们的模型，再任意代入一组变位参数值，通过表中给出的油位高度求出储油量，计算出理论的出油量，再与实际的出油量进行比较，然后不断的调试变位参数使理论与实际出油量之间的误差尽可能小，最终确定一组变位参数。

三、模型假设
1.忽略温度变化对储油罐的容量影响
2.假设油位探针底部固定，且一定过圆心；
3.在进油结束到出油开始的时间段储油罐无变化；
4.大气压对储油罐的容量没有影响；
5.不考虑管道里面的油量对罐容表的影响；
6.忽略储油罐底部和罐壁污垢对油位探针的高度的影响；
四、符号规定
p':表示油罐体的高度；
l:表示储油罐的长度；
h:表示油位高度；
α:表示倾斜角度；
v:表示油罐装油的体积；
R:表示油罐体高度的一半；
p:表示小椭圆型储油罐的高度；
a:表示小椭圆型储油罐横截面的长半轴；
b:表示小椭圆型储油罐横截面的短半轴；
t:表示储油罐油位探针距罐体左端最远距离；
1
t:表示储油罐油位探针距罐体右端最远距离；
2
r:表示储油罐两端球冠体横侧面的半径；
d:表示储油罐两端球冠体内油位实际高度；
L:表示储油罐弓形所确定的圆的半径；
m:表示小椭圆型储油罐油位探针距罐体左壁的距离；
n:表示小椭圆型储油罐油位探针距罐体右壁的距离；
m':表示储油罐油位探针距罐体的圆柱体部分左壁的距离；
n ':表示储油罐油位探针距罐体的圆柱体部分右壁的距离；
),(R Q O :表示储油罐左弓形所确定圆的圆心坐标；
),'('R Q O :表示储油罐右弓形所确定圆的圆心坐标；
五、模型的建立
5.1小椭圆型油罐储油量与油位高度的关系模型
5.1.1罐体无变位时储油量与油位高度关系的模型
对问题一，要分别计算出油罐没有变位和发生纵向变位（变位角度为4.1°）时小椭圆型储油罐中储存油量与高度的关系，首先当没有发生变位的时候，油罐里面的油构成了一个规则体，规则体的横截面如图所示：
图 1
根据柱体的体积公式，要计算储油量就需要求出图中阴影部分的面积，再乘上油罐的长度。

根据相关的数学知识建立如图 1的直角坐标系，在知道椭圆的长轴和短轴的数据之后，阴影部分的面积与高度之间的函数关系为：
()⎰---=b h b dy y b b
a h S 222 ⎪⎭
⎫ ⎝⎛+-+=)(arcsin )(1)(22h f h f h f ab π 其中 b
b h h f -=)( 由符号定义可知m ,n ,l 之间的关系为：
n m l += 于是体积与油位高度之间的函数关系为：
()h lS h v =)(
即罐体未变位的模型为：
⎥⎦
⎤⎢⎣⎡+-+=)(arcsin )(1)(2)(2h f h f h f abl h v π 其中 b
b h h f -=)( 5.1.2 变位罐体储油量的模型建立
使用一段时间后的储油罐因地基变形等原因导致罐体的位置发生纵向倾斜或横向偏转等变化，从而引起罐容表发生变化。

当小椭圆型储油罐纵向变位α角度时，油罐的横截面弓形的高度d 与坐标x 成一次函数关系（如图 2）。

h xtan x d +-=α)(
图 2
将d 的高度带入前面的椭圆的弓形面积公式，得侧面面积为：
()dy y b b
a S
b x d b ⎰---=222 ⎪⎭
⎫ ⎝⎛+-+=)(a r c s i n )(1)(22x g x g x g ab π 其中 b
b x d x g -=)()( 然后对弓形的面积进行积分就是体积，但是对于不同的情况要分别讨论，首先通过matlab 程序编程(程序参见附件1)对一些特殊油位高度计算油罐的容油量，结果如表1：
表 1
油面特殊高度如图3所示：
图 3
经过对表中的数据比较发现，要讨论变位后对罐容表的影响，需分为以下几种情况考虑：
1.当h=0时即图2区域I 以下部分，此时油面高度不会随着加入油量的多少而改变，即此时的油位高度为0，由表知在油位高度为0时储油量为1.6744L ，而罐内油量初值为215L ，远大于1.6744L ，故此情况无需研究。

2.当油位高度大于0且小于等于1h 时，即为图（3）中的I 区域：
dx S h v h
m ⎰-=α
tan )(
3.当油位高度大于1h 且小于等于2h 时，即为图（3）中的II 区域：
dx S h v m ⎰-=n
)( 4.当油位高度大于2h ，且小于等于p 时，即为图（3）中的II I 区域：
()()h g m ab dx S h v h -+=⎰-πn
g .)(）（ 其中 ()α
tan h p h g -= 5.在油位高度达到m 2.1以后继续往里面注入油，此时油位高度已经达到最大值 ，也不会随着注入量的增加而增加，即此时的油位高度始终是1.2m ，无研究意义 综上所述，罐体纵向倾斜 1.4=α对罐容表影响的模型为：
()()()⎪⎪⎪⎩⎪⎪⎪⎨⎧≤≤-=-+≤≤≤≤=⎰⎰⎰---p h h h h h h h 0h v 2n g 21n 1tan απαtan h p h g h g m ab dx S d S dx S h m h m ，）（）（
5.2变位后的油罐储油量与油位高度的关系模型
5.2.1倾斜油罐储油量与油位高度的关系模型
储油罐体纵向倾斜α角度时（如图 3）
图 3
储油罐罐体纵向倾斜α角度时，球罐圆柱体部分的油面的高度d 是关于x 的函数，即
h x +-=αtan d(x)
1.首先计算油罐较高一端凸头部分的储油量1v ，经分析发现，凸头部分的储油量分两种情况考虑：
(1) 当R y ≤1时，易求得油面所在直线h x x f +-=αtan )(和凸头部分曲线()2
2)(Q x L R x y ---=的交点()11,y x A ，凸头部分储蓄的油的侧面积为： ()dy y x r S x r x d x r ⎰---=)
()(22)(2
⎪⎭
⎫ ⎝⎛+-+=)(a r c s i n )(1)(2)(22x g x g x g x r π 其中)
()()()(x r x r x d x g -=，)()(),()()(x y R x r x y x f x d -=-=，然后再对面积积分即得球冠体部分的储油量；
()⎰'
-=m x dx S h v 11
(2) 当R y >1时，（如图 3）在1t -到1x -之间的球冠体部分储蓄的油的侧面是整个圆面积，而在1x -到m '-之间侧面只是圆的一部分，所以球冠体部分的储油量为两部分之和，即为：
()()⎰⎰'--+-=m x x t dx S dx x y R h v 11121)(π
()()
()⎰'-+⎪⎭⎫ ⎝⎛--+---+=m x dx S x t Q x x t t Q L x t 111211*********π 2.同理油罐较低一端凸头部分的储油量为曲线方程为()2
2)(Q x L R x y '---='，且储油量3v 为： 当()R n f ≤'时
()dx S h v x n ⎰'=2
3 当()R n f >'时
()()dx x y R dx S h v t x x n ⎰⎰'-+='2
222
3)(π 即：
()()
()⎪⎭⎫ ⎝⎛+'+++---+=⎰'222222222222331)(2x t Q x x t t Q L x t dx S h v x n π 通过d 我们易求的圆柱体内油的侧面积为：
()dy y R x S R
x d R ⎰---=222)(
3.中间圆柱体部分里面的容油量的计算模型如下：
由模型一可知分三种情况讨论，第一种情况是在油位高度为0以下部分，第二种情况是油位高度h 在0到m 3之间，第三种情况为油位高度为3m 时，仍然往储油的面加油。

显然第一与第三种情况无实际研究意义，第二种情况飞模型为：
当αtan 0n h '<<时
()()dx x s h v h
m ⎰'-=α
tan 2
当ααtan tan m P h n '-'≤≤'时
()()dx x s h v n m ⎰'
'-=2 当P h m P '<<'-'αtan 时
()()()()()x t m R dx x s h v n x t -'+=⎰'
-22π
其中 ()αt a n h p x t -'=
.
综上所述，倾斜油罐储油量的模型为：
)()()()(321h v h v h v h v ++=
5.2.2横向偏转油罐储油量与油位高度的关系模型
储油罐罐体横向偏转β角度时，截出的油位高度为h ，因为偏转β角度非常小，所以实际的油位高度h '可以近似的取值，即m h 5.1≤时βcos /h h ≈'；m h 5.1>时.cos βh h ≈'下面就分两种情况讨论：
1.当R h ≤'时，（如图 4），油浮子所在的直线和圆弧的交点()11,y x A 易求，故横侧面的面积即为：
dy y x r S x r x d x r ⎰---=)
()()(22)(2 ⎪⎭
⎫ ⎝⎛+-+=)(arcsin )(1)(2)(22x g x g x g x r π 其中)
()()()(x r x r x d x g -=
，)()(x y h x d -'=，)()(x y R x r -=，然后再对面积积分即得球罐凸头部分的储油量； ()dx S h v m x ⎰
'-='11 此时储油罐中圆柱体部分的储油量为：
()dy y R l h v R h R ⎰-'--'='2222
⎪⎭
⎫ ⎝⎛'+'-'+'=)(a r c s i n )(1)(222h p h p h p R l π 其中R
R h h p -'=
')(，n m l '+'='，当球体只发生横向偏转的特殊情况时，球罐两边的凸头部分油量是对称的，因此该情况的模型为： ()()h v h v v '+'=212
图 4
2．当R h >'时
()()dx S dx x y R h v m x x t ⎰⎰'
--+-='11121)(π
()()
()⎰'-+⎪⎭⎫ ⎝⎛--+---+=m x dx S x t Q x x t t Q L x t 1
11211121221131π 此时油罐中圆柱体部分的储油量的模型不变，仍为： ()dy y R l h v R
h R ⎰-'--'='2222
⎪⎭
⎫ ⎝⎛'+'-'+'=)(a r c s i n )(1)(222h p h p h p R l π 故此时总油量的模型仍为：
()()h v h v v '+'=212
5.2.3纵横向同时变位时油罐储油量与油位高度的关系模型
储油罐同时发生纵向倾斜和横向偏转时，横向偏转β角度引起油位高度发生变化，但角度β非常小，可以近似的处理为实际油位高度h '，纵向横向同时变位只需把实际油位高度h '带入5.2.1的模型中（即h h '=）即可。

此时的模型为
)()()()(321h v h v h v h v ++=
其中R h ≤时βcos /h h =；R h >时βcos h h =.
六、模型求解
6.1问题一的模型求解
小椭圆型储油罐纵向倾斜︒=1.4α时，根据建立的储油量与油位高度关系的模型，求解出罐体变位后油位高度间隔1cm 的罐容表标定值如表2（程序参见附件2）：
表 2
6.2问题二的模型求解
对实际的储油罐建立了储油量与油位高度及变位参数（纵向倾斜角度α和横向偏转角度β）之间的一般关系的模型，根据题目给出的附件2中集数据表通过逐步调试
算法来确定变位参数。

逐步调试算法思想如下：
首先考虑变位参数,101︒≤≤︒α,101︒≤≤︒β以︒1为间隔得到100组变位参数，再将油位高度代入模型求解出对应的100组理论的容油量准确值, 将求解出的容油量相邻之间的差值与对应出油量进行比较。

在不考虑异常点的情况下，当这个比值越小，说明给定的位变参数越准确，否则反之。

在这100组数据中就得到比值最小的一组
()︒︒3,2，
再将()︒︒3,2设置为初始值代入模型中，根据表中的油位高度算出理论储油量，相邻理论储油量的差即是理论出油量，通过理论出油量与实际出油量之间的相对误差来调试变位参数，此时的变位参数所对应的最大误差为w 。

我们固定α值，以︒1.0的间隔来不断的调试β，得到一组变位参数()1,βα使得最大误差1w 较小（w w <1），然后再固定1β的值，不断调试α，又得到一组变位参数()11,βα使得最大误差2w 较小（12w w <），这样不断的重复此过程，最终得到一组的变位参数()βα'',使得最大误差相对最小。

此时的()βα'',即是要确定的最优变位参数。

将题中附件2中的数据分为前300个数据和后300个数据两组，再通过此算法思想调试变位参数，如表3：（程序参见附件3）
表3：
α
β
前300 后300 α
β
前300 后300 1.5 3.2 5.2503 4.0736 2 3.1 3.989 4.2492 1.5 3.5 5.291 4.2492 2 3.2 3.6732 5.1767 1.6 3.2 4.8772 4.14 2 3.3 3.6121 4.2342 1.6 3.5 4.926 4.1277 2 3.4 3.5549 4.2188 1.7 3.2 4.4996 4.1157 2 3.5 3.3594 4.2029 1.7 3.5 4.5484 4.1041 2.1 3.2 3.9215 4.0772 1.8 3.2 4.1175 4.0502 2.1 3.3 3.907 4.1482 1.8 3.5 4.2177 4.0639 2.1 3.4 3.8706 4.0901 1.9 3.2 3.7794 4.0613 2.1 3.5 3.892 4.2917 1.9 3.3 3.7626 4.2638 2.1 3.5 4.1662 5.1266 1.9 3.4 3.7664 4.2492 2.1 3.6 3.8607 4.6217 1.9 3.5 3.7306 4.4801 2.1 3.7 3.8443 4.2296 2 2.5 3.5699 4.2039 2.1 3.8 3.8276 4.1681 2
3
3.5844
4.2653
2.2
3.3
4.2481
4.278
最终确定的变位参数为︒=︒=2.3,9.1βα。

将此时的变位参数代回模型中，得到储油量与油位高度的关系的关系模型，易得油位高度间隔为10cm 的罐容表标定值，如表4（程序参见附件4）
表4 ：
罐体变位后油位高度间隔为10cm的罐容表标定值
七、结果分析和检验
7.1问题一的模型检验与结果分析
对问题一中无变位罐体储油量的模型，罐体是一个标准的椭圆柱体，用罐内油的横截面积乘上罐体的长度求储油量是精确算法，将题目中附表1内的无变位进油的油位高度代入模型求得相应的储油量，与实际的储油量进行定性比较（如图 5）（程序参见附件5）
图 5
从图中可以看出理论储油量与实际储油量的曲线相差不大，我们再通过误差定量比较，误差的计算方式为：
|理论储油量—实际储油量|
相对误差= —————————————⨯100%
实际储油量
由此公式求得最大误差为3.4917%，在误差允许范围内，故此模型是符合实际的。

α，罐体的储油量是通过对罐内油的横截面积积分小椭圆型储油罐纵向倾斜︒
=1.4
建立的模型，现在仍然还是通过定性和定量的方法来检验模型。

图6是变位后的理论油量与实际储油量随油位高度的变化曲线（程序参见附件6）。

图 6
通过对相对误差的计算，最大误差为：5.4659%,通过对曲线发现，在曲线的中间段，理论值与实际值几乎平行，产生这个现象的原因在模型求解计算的过程中的近似替代以及忽略的一些不可控因素。

但是他们之间的差异很小，所以在误差允许范围内，变位后建立的罐内储油量模型是符合实际的。

7.2问题二的模型检验与结果分析
对问题二中建立的球罐内储油量与油位高度以及变位参数之间的关系的模型，
在建立模型之前通过对数据的比较发现油罐此时显示的容油量与出油量之间的关系不正确，而在实际过程中统计的出油量数据相对显示容油量要准确，所以说明显示的
α(即罐体不变位)时，将题目中给出容油量不是罐内容油量的准确值，而当0
=β
=
的附件2中油位高度代入模型求得的理论储油量与实际储油量几乎相等，而题目给出的是变位后的实际检截数据，进一步说明表中的油量容积不是罐体内油量的准确值，同时说明油量容积是罐体在没有发生变位的情况下，通过油位的高度建立的模型计算得到的。

现在就通过建立罐体不变位时的模型验证这一现象：
罐体不变位时罐内储油量的模型与5.2.2中偏转油罐储油量的模型相同，只是其中
'。

再将表中的油位高度代入此时的模实际油位高度就是截量出的油位高度，即h
h=
型，得出的理论储油量与表中的油量容积也几乎相等，由此说明模型建立的准确性。

通过对截试出来中的数据（部分数据在表三中）比较发现，当α固定时，β的变化对整个模型影响相对较小，而当β固定时，α的变化对整个模型相对较大，这一现象比较符合实际情况，同时说明该模型在实际应用的过程中要更多的考虑α的值。

该模型通过逐步调试求解出来的最终变位参数(α,β)为)
︒时，此时计算出来的相
2.3,
9.1(︒
对误差为4.0613%，在允许的误差范围内，说明模型准确。

八、模型的评价和推广
优点：
1.在对圆或者椭圆面里面的弓形进行积分时，遇到不能够直接积分的时侯，采用了契比切夫多项式逼近法来进行近似化处理，使积分能够算出并且使得出的结果符合实际。

2.对于本文中两问题来确定的模型，根据题目所给的数据运用Mtlab软件求得结果。

经过对模型求解的数据仔细分析，不难发现各项数据与实际情况相符合。

因此模型在满足一定的条件下是准确的，同时模型又很直观，容易理解，很利于应用和推广。

3.本文用到的数学方法（积分）都比较简单，适合推广。

4.由于使用到的数据都是显示中的数据，对数据通过分析建立模型，因此得出的模型具有可操作性，该模型对满足这些相关条件加油站具有较强的指导意义。

缺点：
1.对罐容表的表示表定值的影响，除考虑到罐体的倾斜和横向偏转外，还有其他一些因素对罐容表的表示表定值的影响，如果将这些因素全部都考虑，这将增加需要处理的数据量。

鉴于时间关系，本文不可能做出更加精细的描述。

2.本文在对问题求解时，均因条件限制，不能够准确的计算出积分，采用逼近法来进似代替求积分，因此结果均存在一定误差。

推广：
当目前环境下，客户的需求量快速增长，加油站如何在这样的社会大背景下满足客户的需求量，以及同时使客户的满意度更高，这是对每个加油站都需要知道储油量的多少，我们建立的这类数学模型在实际情况下能使自己知道自己的储油量，对类似的服务业也适用，只要给出的数据足够，实际和精确，则模型得出的结果将具有很强的实际意义。

九、参考文献
[1] 姜启源，数学模型，北京：高等教育出版社，2004年。

[2] 孙祥徐流美吴清，MATLAB7.0基础教程，清华大学出版社，2005-5。

[3] 苏金明阮沈勇王永利，matlab 工程数学，北京：电子工业出版社，2005年。

[4] 王沫然，matlab与科学计算，北京：电子工业出版社，2005，年。

[5] 孙宏达关进波，用逼近法计算横截面为椭圆形(圆形)储油罐的储油体积，
/view/74604c7931b765ce05081474.html，2010-9-10。

十、附件
附件1：
%体积计算的函数 tj.m
function dd=tj(d)
a=0.89;
b=0.6;
c=4.1*pi/180;
syms x y h
h=-tan(c)*x+d;
v=(2*a)/b*sqrt(b^2-y^2);
if d==0
dd=int(int(v,y,-b,tan(pi-c)*x+d-b),x,-0.4,0);
elseif d>0 & d<=0.1649
dx=d/tan(c);
dd=int(int(v,y,-b,tan(pi-c)*x+d-b),x,-0.4,dx);
elseif d>0.1649 & d<=1.1713
dd=int(int(v,y,-b,tan(pi-c)*x+d-b),x,-0.4,2.05);
elseif d>1.1713 & d<=1.2
m=(d-1.2)/tan(c);
dd=int(int(v,y,-b,tan(pi-c)*x+d-b),x,m,2.05)+(m+0.4)*pi*a*b; else
end
%主程序ywzcx1.m
clear all;
r(1)=tj(0);
r(2)=tj(0.1649);
r(3)=tj(1.1713);
r(4)=tj(1.2);
v=double(r)
附件2
%问题1体积计算函数1 ywtj1.m
function v=ywtj1(d)
syms x m h dd y
c=4.1*pi/180;
h=-tan(c)*x+d;
dd=2*0.89/0.6*sqrt(0.6^2-y^2);
v=int(int(dd,y,-0.6,h-0.6),x,-0.4,d/tan(c));
%问题1体积计算函数2 ywtj2.m
function v=ywtj2(d)
syms x m h
c=4.1*pi/180;
h=-tan(c)*x+d;
m=0.6*0.89*(pi/2+(h-0.6)/0.6*sqrt(1-((h-0.6)/0.6)^2)+asin((h-0.6)/0.6)); v=int(m,x,-0.4,2.05);
%问题1体积计算函数3 ywtj3.m
function v=ywtj3(d)
syms x m h
c=4.1*pi/180;
h=-tan(c)*x+d;
m=0.6*0.89*(pi/2+(h-0.6)/0.6*sqrt(1-((h-0.6)/0.6)^2)+asin((h-0.6)/0.6)); dx1=(d-1.2)/tan(c);
v=int(m,x,dx1,2.05)+pi*0.6*0.89*(0.4-(1.2-d)/tan(c));
主程序ywzcx2.m
clear;
clc;
r=[];
for i=1:14
d1=i/100;
w1=ywtj1(d1);
r=[r,w1];
end
for i=15:116
d2=i/100;
w2=ywtj2(d2);
r=[r,w2];
end
for i=117:120
d3=i/100;
w3=ywtj3(d3);
r=[r,w3];
end
附表3：
%圆柱体积函数 yztj.m
function v=yztj(k,b)
if k==0
syms ss x y
ss=2*sqrt(1.5^2-y*y);
v=int(int(ss,y,-1.5,k*x+b-1.5),x,-2,6);
else
m=-(b/k);
syms ss x y
ss=2*sqrt(1.5^2-y*y);
if b<=(-6*k)
v=int(int(ss,y,-1.5,k*x+b-1.5),x,-2,m);
elseif b<(3+2*k)
v=int(int(ss,y,-1.5,k*x+b-1.5),x,-2,6);
else
v=int(int(ss,y,-1.5,k*x+b-1.5),x,(3-b)/k,6)+pi*1.5^2*(2+(3-b)/k);
end
end
%左凸球冠体积函数zqbtj.m
function v=zqbtj(k,b,dx)
if dx>=-2
v=0;
else
dy=k*dx+b;
syms x y r h dd ss r1 s
dd=k*x+b;
h=dd-1.5+sqrt(1.625^2-(x+1.375)^2);
r=sqrt(1.625^2-(x+1.375)^2);
y=(h-r)/r;
ss=r^2*(pi/2+16/(1575*pi)*(615*y-80*y^3-48*y^5));
r1=sqrt(1.625^2-(x+1.375)^2);
s=pi*r1^2;
if dy>1.5
v=int(s,x,-3,dx)+int(ss,x,dx,-2);
else
v=int(ss,x,dx,-2);
end
end
%右凸球冠体积函数yqbtj.m
function v=yqbtj(k,b,dx)
if dx<=6
v=0;
else
dy=k*dx+b;
syms x y r h dd ss r1 s
dd=k*x+b;
h=dd-1.5+sqrt(1.625^2-(x-5.375)^2);
r=sqrt(1.625^2-(x-5.375)^2);
y=(h-r)/r;
ss=r^2*(pi/2+16/(1575*pi)*(615*y-80*y^3-48*y^5));
r1=sqrt(1.625^2-(x-5.375)^2);
s=pi*r1^2;
if dy>1.5
v=int(s,x,dx,7)+int(ss,x,6,dx);
else
v=int(ss,x,6,dx);
end
end
%球冠与油平面交点横坐标函数hzb.m function x=hzb(k,b,dx,dy)
d=abs(k*dx-dy+b)/sqrt(k*k+1);
x0=(k*dy+dx-k*b)/(k*k+1);
y0=(k*k*dy+k*dx+b)/(k*k+1);
r=1.625;
xuch=r*r-d*d;
if dx==-1.375
x=x0-sqrt(xuch/(k*k+1));
end
if dx==5.375
x=x0+sqrt(xuch/(k*k+1));
end
%第二问确定偏转角度主程序ew.m clear all;
close all;
clc;
a=xlsread('ew.xls');
h=a(:,1);
h=h./1000;
cy=a(:,2);
cy=cy./1000;
a=2;
beta=10*pi/180;
c=a*pi/180;
k=-tan(c);
m=304; % m 至少取2
f=h(m-1:end);
tt=cy(m:end);
% tt 表示实际出油量
for i=1:length(f)
if f(i)<1.5
f(i)=f(i)/cos(beta);
else
f(i)=f(i)*cos(beta);
end
x1=hzb(k,f(i),-1.375,1.5);
x2=hzb(k,f(i),5.375,1.5);
v1(i)=zqbtj(k,f(i),x1);
v2(i)=yqbtj(k,f(i),x2);
v(i)=yztj(k,f(i))+v1(i)+v2(i);
end
v=double(v);
% wc=[];
% zsc=[];
for i=2:length(f)
zsc(i-1)=-(v(i)-v(i-1));
% zsc表示理论出油量
end
for i=1:length(f)-1
wc(i)=abs(zsc(i)-tt(i))/tt(i)*100; end
% 抛出异常的点
[zdwc,n]=max(wc);
wc(n)=[];
zdwc=max(wc);
附件4
%第二问求标度的主程序 ewer.m clear all;
clc;
f=10:10:300;
f=f./100;
a=1.9;
beta=3.2*pi/180;
c=a*pi/180;
k=-tan(c);
syms v1 v2 v
for i=1:30
if f(i)<1.5
f(i)=f(i)/cos(beta);
else
f(i)=f(i)*cos(beta);
end
x1=hzb(k,f(i),-1.375,1.5);
x2=hzb(k,f(i),5.375,1.5);
v1(i)=zqbtj(k,f(i),x1);
v2(i)=yqbtj(k,f(i),x2);
v(i)=yztj(k,f(i))+v1(i)+v2(i); end
v=double(v);
v=v*1000;
v=v';
%plot(f,v,'r-');
%fenxi.m
clear;
clc;
b=xlsread('untitled3.xls');
a=2.0;
b=b/1000;
beta=3.5*pi/180;
c=a*pi/180;
k=-tan(c);
m=500;
f=b(m:m+9);
for i=1:10
if f(i)<1.5
f(i)=f(i)/cos(beta);
else
f(i)=f(i)*cos(beta);
end
x1=hzb(k,f(i),-1.375,1.5);
x2=hzb(k,f(i),5.375,1.5);
v1(i)=zqbtj(k,f(i),x1);
v2(i)=yqbtj(k,f(i),x2);
v(i)=yztj(k,f(i))+v1(i)+v2(i);
end
v=double(v);
plot(f,v,'r-');
附件5：dddf.m
function v=dddf(h)
v=89/2301893700960125939166572010851114631246643200000*(-62725908771887790 02551381099945+47343316306301763082307218964480*h-31691265005705735037417 580134400*h^2)^(3/2)-267/201763772073550000*(-627259087718877900255138109994 5+47343316306301763082307218964480*h-31691265005705735037417580134400*h^2 )^(1/2)+13083/20000*pi+1122716626669857009/201763772073550000*asin(-420493118 6029427/3377699720527872+5/3*h)-18788454695436288/2522047150919375*asin(-420 4931186029427/3377699720527872+5/3*h)*h-89/44958861346877459749347109586935 83264153600000*(16630349716423503158784128615+565815492334254236105829253 120*h-495176015714152109959649689600*h^2)^(3/2)+267/25220471509193750*(1663 0349716423503158784128615+565815492334254236105829253120*h-49517601571415 2109959649689600*h^2)^(1/2)-107343635458253943/25220471509193750*asin(-402036 087858629/422212465065984+5/3*h)+18788454695436288/2522047150919375*asin(-40 2036087858629/422212465065984+5/3*h)*h;
%untitled2.m
clear;
clc;
a=xlsread('untitled2.xls');
m=a(:,1);
h=a(:,2);
for i=1:53
m(i)=(m(i)+215)/1000;
h(i)=h(i)/1000;
c(i)=dddf(h(i));
f(i)=abs(c(i)-m(i))/m(i);
end
plot(h,m,'ro')
hold on;
fplot('dddf',[0,1.2])
r=max(f);
附件6：
%函数tj.m
function dd=tj(d)
a=0.89;
b=0.6;
c=4.1*pi/180;
syms x y h
h=-tan(c)*x+d;
v=(2*a)/b*sqrt(b^2-y^2);
if d==0
dd=int(int(v,y,-b,tan(pi-c)*x+d-b),x,-0.4,0);
elseif d>0 & d<=0.1649
dx=d/tan(c);
dd=int(int(v,y,-b,tan(pi-c)*x+d-b),x,-0.4,dx);
elseif d>0.1649 & d<=1.1713
dd=int(int(v,y,-b,tan(pi-c)*x+d-b),x,-0.4,2.05);
elseif d>1.1713 & d<=1.2
m=(d-1.2)/tan(c);
dd=int(int(v,y,-b,tan(pi-c)*x+d-b),x,m,2.05)+(m+0.4)*pi*a*b; else
end
%jifen.m
clear;
clc;
a=0.89;
b=0.6;
c=4.1*pi/180;
syms x y h
v=(2*a)/b*sqrt(b^2-y^2);
dd=int(int(v,y,-b,tan(pi-c)*x+h-b),x,-0.4,2.05)。