数模优秀论文：储油罐的变位识别与罐容表标定.doc

储油罐的变位识别与罐容表标定
专家点评:
本文基于所给数据准确、罐体几何形状因有附属构件而含有误差进而导致推导的罐容与油位高度之间函数关系的理论公式含有较大偏差的理解下，通过对理论公式计算结果与实测数据的偏差的曲线拟合，对小椭圆型储油罐给出了修正的罐容表。文中分析研究了无变位和有纵向变位的小椭圆型储油罐的罐容与油位高度的函数表达式、有纵向变位和横向变位的实际储油罐罐容与油位高度的函数表达式、以储油罐中油量随高度的变化率为依据识别纵向倾斜角度和横向偏转角度，由此给出了罐容表的标定、检验了所给出的数学模型的正确性和可靠性，思路正确、方法有效、所得结果合理，但是，对问题一利用祖暅原理将有变位近似转换为无变位的方法略欠妥当。

中国海洋大学 曹圣山教授

摘要

对于两端平头的小椭圆型储油罐与实际球冠封口的储油罐，本文分别建立了相应的数学模型，解决了储油罐变位后的识别和罐容表的标定的相关问题。在建立两个模型的过程中充分的运用了MATLAB和EXCEL两个软件，利用祖暅原理将变位容积计算转换为未变位时的计算，在保证精度情况下，避免了复杂的积分运算。
对于模型1，首先，我们通过积分，得出无变位时的储油量与油位高度关系，此时，所得理论容积与实测容积出现由罐内附属构件占有一定体积造成的偏差，及时的运用曲线拟合的方法获得了其偏差函数，对模型1的公式进行了修正，获得了很好的结果。在变位条件下，依据油位高度，对变位后的小椭圆形储油罐划分了三种高度条件来讨论了其罐容标定，然后利用几何关系将高度转化为无变位条件下的高度来计算容积。
对于模型2，无变位时，同样，我们先积分，积出无变位情况下实际油罐的储油量与油位高度表达式；变位时，我们依然依据油位高度，对实际的球冠封口的储油罐划分为三种情况来讨论，同样采用一些转化将高度转化为无变位条件下的高度来计算容积；在求解α,β的过程中，利用导数间的关系建立了油位高度的关系，编写了导数返查的MATLAB程序以及依据数值逼近思想所利用的迭代公式和最小二乘法的线性拟合，精确地计算出了α，β的值 ，进而促成模型2的正确建立，然后利用模型计算出罐容标定表并利用给定数据分析检验。

关键词： 祖暅原理 数值逼近思想 最小二乘法 迭代 模型

1．问题重述
通常加油站都有若干个储存燃油的地下储油罐，并且一般都有与之配套的“油位计量管理系统”，采用流量计和油位计来测量进/出油量与罐内油位高度等数据，通过预先标定的罐容表（即罐内油位高度与储油量的对应关系）进行实时计算，以得到罐内油位高度和储油量的变化情况。
许多储油罐在使用一段时间后，由于地基变形等原因，使罐体的位置会发生纵向倾斜和横向偏转等变化（以下称为变位），从而导致罐容表发生改变。按照有关规定，需要定期对罐容表进行重新标定。原文图1是一种典型的储油罐尺寸及形状示意图，其主体为圆柱体，两端为球冠体。原文图2是其罐体纵向倾斜变位的示意图，原文图3是罐体横向偏转变位的截面示意图。
请你们用数学建模方法研究解决储油罐的变位识别与罐容表标定的问题。
（1）为了掌握罐体变位后对罐容表的影响，利用如图4的小椭圆型储油罐（两端平头的椭圆柱体），分别对罐体无变位和倾斜角为a=4.10的纵向变位两种情况做了实验，实验数据如附件1所示。请建立数学模型研究罐体变位后对罐容表的影响，并给出罐体变位后油位高度间隔为1cm的罐容表标定值。
（2）对于图1所示的实际储油罐，试建立罐体变位后标定罐容表的数学模型，即罐内储油量与油位高度及变位参数（纵向倾斜角度a和横向偏转角度b ）之间的一般关系。请利用罐体变位后在进/出油过程中的实际检测数据（原文附件2），根据你们所建立的数学模型确定变位参数，并给出罐体变位后油位高度间隔为10cm的罐容表标定值。进一步利用附件2中的实际检测数据来分析检验你们模型的正确性与方法的可靠性。

2、符号说明
长度，m；
面积,m2
体积,m3
椭圆管道短半长轴长,m
椭圆管道长半长轴长,m
液面高度,m
误差；
圆筒管道半径,m
球冠封头所对应的球的半径,m

3、基本假设
(1)实际储油罐的偏转角度，角很小；
(2)假设油位探针固定
(3)假设实际油罐内的油位探针装置，注油管，出油管体积不计；
(4)假设实验数据都是可信的

4、问题分析与模型建立

4.1 问题一的求解
4.1.1 无变位储油量与油位高度关系计算

小椭圆油罐模型采用的是两短平头的椭圆柱体，其三维坐标示意图如下所示：
其油位高度H与储油量V的数学关系式推到如下：

; ;
; 则;
又因为，由参考文献[1]可知

； （1）

通过此公式计算出的储油量为空油罐理想值，当实际上该油罐内部还有油位探针、注油管、出油管，是实际值比公式计算值偏小，并且误差将随油位高度增加而增加，经过EXCEL对数据的处理计算，如下图，发现误差确实随油位高度增加而增加

因此须对公式进行修正，于是我们作出了误差随油位高度变化的图像，并对该图像进行了拟合，得出偏差公式：； （2）
由此可得的修正值为：
（3）

4.1.2 油罐纵向变位后储油量与油位高度关系计算
首先下图给出了纵向变位后油罐的剖视图：

由以上两图来计算倾斜油罐的油位高度与储油量的关系，我们设法将倾斜油罐的油位高度，转化为相同储油量下无变位油罐的油位高度，再利用无变位的体积公式计算，由祖暅定理：“等高处横截面积常相等的两个立体，其体积也必然相等”，显然我们将变位时的等腰梯形转化为面积相等的矩形，其体积也必然相等，也就是我们的方法，如下：
(1)液面AB处于与之间时：
；
；
则储油量的纵截面梯形面积计算公式为:

；
即：；[2]
；
(2)当液面AB处于以下时，
用矩形面积与直角三角形面积相等的方法导出H2与H1的关系，这时矩形底边长小于L，矩形和三角形的底边长为，矩形面积，直角三角形面积；
所以； 同时还要将体积计算式中L改为

(3)当液面AB处于以上时，容积采用总容积减去空余部分容积，而空余部分容积可由2计算，其中，空余部分的等效高度为， 则可算出容积。
综上所述:
（4）
将(4)式代入(3)式,可以得到修正后的发生纵向变位的储罐的储油量与油位高度的关系式利用MATLAB进行计算（程序见附录一），可得出标定结果数据见附录二图像如下：

4.2 问题二求解
4.2.1模型的建立：
对于未变位储油罐，我们建立了储油量与油位高度的一般关系式，其简略推导过程如下：
把总体积分成圆筒部分体积与封头部分（两边）体积之和
即 ； （5）

；
；
则；

对于
;
即：； （6）

对于，同样积分可得：[3]
(7)考虑纵形倾斜角度α与横向偏转角度β后，我们将油浮子所测高度h,
等效转化为未变位液面的高度，转化过程如下：
首先将考虑β角的偏转，油罐发生β偏转后液面水平高度并未变化，有几何关系
可以得到液面水平高度与油浮子所测高度h之间的关系为：
； (8)

然后将等效转化为油罐未发生纵向倾斜的高度，分三段讨论如下：
(1)当液面位于以下时，近似液体截面面积 ；
即：;
;
();可得到容积
(2)当液面位于与之间时，；可得到容积
(3)当液面高于时，可利用1中所讨论的结论计算空余容积，再用总容积减去空余容积，直接得到储油量与高度的关系式，其中空余容积截面中 ，从而得出空余容积等效高度，可得到容积.

4.2.2变位参数的确定
在实验数据高度h给定的范围内，我们得出等效半径后

; (9)
即：；
两边求导得：; (10)
因为实际工程中β角一般比较小，这里取，即, 因为，以为步长，得出了由0.25到2.699的一系列导数值，导数表见附录3；由原题附件2的出油量与显示油高，可以计算出高度高度为的近似导数值，利用此的值在附录三导数表中反查出的值（自己编的MATLAB反差程序见附录四），从而得出一系列的实数对 ，在平面直角坐标系中以绘制散点图, 然后用EXCEL进行最小二乘线性拟合，得到=1.0065，显然不合理。根据逼近思想概括出该题的适用公式：，即当取时的最小二乘拟合结果，显然该公式在时收敛。我们又验证其它的值，如=0.998，由拟合曲线得到的值为0.999，当=0.997，有拟合曲线得到的值为0.995。当选取的值为0.9975时发现得到的曲线的拟合公式为 ，图像如下

计算出的值为0.9975，此时迭代收敛。说明此时假设合理，此时 ， 。

4.2.3 罐容标定
由4.2.1与4.2.2中所得的最终储油量与油位高度的关系上式，选取油位高度间隔为10cm对实际油罐灌容进行标定，计算中所采用的MATLAB程序见附录五，实际储油罐罐容标定图如下，（其中标定的数值表见附录六）

4.2.4 模型检测
(1)评价方式一
因为附件2中的显示油量容积是根据油罐未变位时的体积计算公式得来的，因此我们选用了出油量和显示油高两组数据来检验我们模型的正确性和方法的可靠性。因为从数据表无法得到储油罐的初始油量体积，即真实的储油量的体积无法直接读取，因此可选择累积出油量来分析。
由显示油高，利用以上建立的，的储油量计算公式，算出储油量，即我们模型的结果。由数据算出各个油位高度实际出油量的累加值，以油位高度为横坐标，分别以我们模型的结果和出油量的累加值为纵坐标，作出图形如下，显然将出油量累加值曲线向上平移V升，V是初始储油量并且未知。

要以实际检测数据分析检验模型，可以考虑两段曲线的相似程度，为此可以可以在相同油位高度下做差，差值随油位高度变化如下图：

若差值曲线为水平曲线则完全相似，由上图看出我们模型存在误差。这一系列差值的平均值为59167.86，由图看出偏差在[-59,59]，即最大占总体积的0.09%,并可求出这一系列差值的标准差为43.35。说明我们的模型满足要求，建模方法可靠性很好。

(2)评价方式二
该方法是利用模型算出的容积值，得出各高度间的出油量，与真实出油量进行比较，发现差值区间在[-4.7L，3.7L]内，并算出，我们算出ε平均值=0.721%，由此看出我们的模型还是比较精确的。
5模型的评价与该进
5.1模型的优点
(1)综合运用了MATLAB和EXCEL两个软件，求解准确，在运用MATLAB进行储油量变位体积的计算时，精确地得出了函数公式，与实验数据吻合很好。
(2)论文中建立了两端平头的椭圆柱体模型和实际储油罐两个体积计算模型，模型的稳定性较高。
(3)此模型经过实际测量数据的检验，精度较高，对于实际问题具有较高的实用性。

5.2模型的缺点
(1)本文在计算储油罐储油量的过程当中，等效装换高度会给计算结果带来一些误差。

5.3模型的改进
两个模型中采用的方法都是依据油位高度来分段对储油量体积进行计算，但是这样不可避免会造成一定的误差，建模时可以采取整个体积完全积分的方法，这样得到的公式将是极为准确的，当然实施起来也是相当困难的。

6模型的推广
我们建立模型的思想是可以推广到其它类似问题上的，例如对平封头圆筒形卧式容器、椭圆形封头圆筒形卧式容器、碟形封头圆筒形卧式容器等形状的储油罐都可以通过变换封头的积分公式来得到相应的体积计算方法。

7参考文献
[1] 管冀年 赵海， 卧式储油罐罐内油品体积标定的实用方法，计量与测试技术，3:21，2004.
[2] 田铁军, 倾斜卧式罐直圆筒部分的容积计算，现代计量测试，5：33-36，1999.
[3] 蒋心亚 宗光, 各种形状封头的圆筒形卧式容器在不同液面高度时液体体积计算的统一表达式，化工设备与管道，395)：30-34，2002.
[4] 姜启源 谢金星，数学模型 北京 高等教育出版社 2003.

附录一：MATLAB程序
function y=vjs(x)%简单无修正计算公式
y=2.45*1.78/1.2*((x-0.6).*sqrt(x.*(1.2-x))+0.6^2*asin(x/0.6-1)+pi*0.5*0.6^2)

function y=vjs1(x)%有修正公式
y=2.45*1.78/1.2*((x-0.6).*sqrt(x.*(1.2-x))+0.6^2*asin(x/0.6-1)+pi*0.5*0.6^2)-(-79.596*x.^4+142.47*x.^3-75.57*x.^2+150.72*x-14.503)/1000;
function y=vjs2(x)%液面低于CM1有修正公式
y=2*x/0.07168*1.78/1.2.*((x-0.6).*sqrt(x.*(1.2-x))+0.6^2*asin(x/0.6-1)+pi*0.5*0.6^2)-(-79.596*x.^4+142.47*x.^3-75.57*x.^2+150.72*x-14.503)/1000;

x1=0:0.01:0.14;
x2=0.15:0.01:1.17;
x3=1.18:0.01:1.2;
n1=x1+002867;
y1=vjs2(n1/2);
y2=vjs1(x2-1.65/2*0.07168);
y3=3.96-vjs2((0.1469+1.2-x3)/2);
plot(x1,y1,x2,y2,x3,y3)

附录二：罐容表标定值
0
0.0138
0.41
0.96665
0.81
2.6149
0.01
0.0147
0.42
1.0051
0.82
2.6556
0.02
0.0164
0.43
1.0439
0.83
2.6961
0.03
0.019
0.44
1.083
0.84
2.7364
0.04
0.0225
0.45
1.1223
0.85
2.7764
0.05
0.0271
0.46
1.162
0.86
2.8162
0.06
0.0327
0.47
1.2018
0.87
2.8558
0.07
0.0395
0.48
1.242
0.88
2.8951
0.08
0.0475
0.49
1.2823
0.89
2.9341
0.09
0.0569
0.5
1.3228
0.9
2.9728
0.1
0.0675
0.51
1.3636
0.91
3.0113
0.11
0.0796
052
1.4045
0.92
3.0493
0.12
0.093
0.53
1.4456
0.93
3.0871
0.13
0.108
0.54
1.4869
0.94
3.1245
0.14
0.1244
0.55
1.5283
0.95
3.1615
0.15
0.1434
0.56
1.5698
0.96
3.1982
0.16
0.1656
0.57
1.6114
0.97
3.2344
0.17
0.189
0.58
1.6532
0.98
3.2702
0.18
0.2133
0.59
1.695
0.99
3.3056
0.19
0.2387
0.6
1.737
1
3.3405
0.2
0.2649
0.61
1.779
1.01
3.3749
0.21
0.2921
0.62
1.821
1.02
3.4089
0.22
0.32
0.63
1.8631
1.03
3.4423
0.23
0.3488
0.64
1.9052
1.04
3.4751
0.24
0.3783
0.65
1.9474
1.05
3.5074
0.25
0.4085
0.66
1.9895
1.06
3.5391

0.26
0.4393
0.67
2.0317
1.07
3.5701
0.27
0.4709
0.68
2.0738
1.08
3.6005
0.28
0.503
0.69
2.1159
1.09
3.6302
0.29
0.5358
0.7
2.158
1.1
3.6592
0.3
0.5691
0.71
2.2
1.11
3.6874
0.31
0.6029
0.72
2.242
1.12
3.7148
0.32
0.6373
0.73
2.2839
1.13
3.7414
0.33
0.6722
0.74
2.3257
1.14
3.7671
0.34
0.7075
0.75
2.3674
1.15
3.7919
0.35
0.7433
0.76
2.409
1.16
3.8156
0.36
0.7796
0.77
2.4505
1.17
3.8383
0.37
0.8162
0.78
2.4918
1.18
3.8566
0.38
0.8533
0.79
2.533
1.19
3.8728
0.39
0.8907
0.8
2.5741
1.2
3.8875
0.4
0.9285

附录三（导数表）
0.001
16.65
0.014
16.94805
0.027
17.24041
0.002
16.67362
0.015
16.97072
0.028
17.26269
0.003
16.69666
0.016
16.99336
0.029
17.28494
0.004
16.71966
0.017
17.01597
0.03
17.30716
0.005
16.74264
0.018
17.03855
0.031
17.32935
0.006
16.76559
0.019
17.0611
0.032
17.35151
0.007
16.7885
0.02
17.08362
0.033
17.37364
0.008
16.81139
0.021
17.10611
0.034
17.39574
0.009
16.83424
0.022
17.12857
0.035
17.41781
0.01
16.85706
0.023
17.151
0.036
17.43985
0.011
16.87986
0.024
17.17339
0.037
17.46186
0.012
16.90262
0.025
17.19576
0.038
17.48384
0.013
16.92535
0.026
17.2181
0.039
17.5058
……………………………………………………………………
2.158
19.37924
2.172
19.10495
2.186
18.82537
2.159
19.35982
2.173
19.08515
2.187
18.8052
2.16
19.34038
2.174
19.06533
2.188
18.785
2.161
19.3209
2.175
19.04548
2.189
18.76477
2.162
193014
2.176
19.02561
2.19
18.74452
2.163
19.28188
2.177
19.00571
2.191
18.72423
2.164
19.26233
2.178
18.98578
2.192
18.70392
2.165
19.24275
2.179
18.96582
2.193
18.68358
2.166
19.22314
2.18
18.94584
2.194
18.66322
2.167
19.20351
2.181
18.92583
2.195
18.64283
2.168

19.18385
2.182
18.90579
2.196
18.62241
2.169
19.16416
2.183
18.88573
2.197
18.60196
2.17
19.14445
2.184
18.86564
2.198
18.58148
2.171
19.12471
2.185
18.84552
2.199
18.56098

附录四（导数返查MATLAB程序）
load '问题A附件2：实际采集数据表.xls'
load ' 导数表.xsl'
d=Sheet3(:,6);d(1)=[];
b=Sheet3(:,9)b(1)=[];
a=data(:,10);a(1)=[];
for i=1:140,
for k=1:1054,
if (b(2199-k)<=a(i))&(b(2199-k+1)>=a(i))|(b(2199-k)>=a(i))&(b(2199-k+1)<=a(i))
c(i)=d(2199-k);
end
end
end
for i=141:301,
for k=1055:2198,
if (b(2199-k)<=a(i))&(b(2199-k+1)>=a(i))|(b(2199-k)>=a(i))&(b(2199-k+1)<=a(i))
c(i)=d(2199-k);
end
end
end
for i=303:440,
for k=1:1054,
if (b(2199-k)<=a(i))&(b(2199-k+1)>=a(i))|(b(2199-k)>=a(i))&(b(2199-k+1)<=a(i))
c(i)=d(2199-k);
end
end
end
for i=441:602,
for k=1055:2198,
if (b(2199-k)<=a(i))&(b(2199-k+1)>=a(i))|(b(2199-k)>=a(i))&(b(2199-k+1)<=a(i))
c(i)=d(2199-k);
end
end
end

附录五
%将油浮子高度转化为液面高度
function y=datahx(x)
y=(x-1.5)*cos(4.0523*pi/180)+1.5;
end
%将高度转化为等效高度
function y=datazh(x)
rx/tan(1.9244*pi/180)+sqrt(1.625^2-(1.5-x).^2)-0.625+x/sin(1.9244*pi/180))/2
y=r.^2*1.9244/180*pi/8/2;
end
%计算两端球冠容积与高度的关系
function y=vjs(x)
if x>1.5
x1=sqrt(1.625^2-(x-1.5).^2);
y=2*pi/3*(2*1.625+0.625)-pi*1.625^2*(x1-0.625)+pi/3*(x1.^3-0.625^3)...
+(x-1.5).*(x1.*sqrt(1.625^2-x1.^2)-0.625*sqrt(1.625^2-0.625^2))+(x-1.5).*(1.625^2+(x-1.5).^2/3).*(asin(x1/1.625)-asin(0.625/1.625))+(x-1.5).*(pi*x1.^2/2-0.625*sqrt(x1.^2-0.625^2)-x1.^2.*asin(0.625/x1));
else y=8.1158-vjs(3-x);
end
%计算中间圆柱容积与高度的关系
function y=vjs3(x)
y=8*1.5^2*(pi-acos((x-1.5)/1.5)+(x-1.5)/1.5^2*sqrt(x*(3-x)));

附录六（实际储油罐罐容表定表）
0
0.0533
1.1
19.39
2.1
46.989
0.1
0.387
1.2
22.078
2.2
49.541
0.2
1.1235
1.3
24.824
2.3
51.985
0.3
2.3389
1.4
27.614
2.4
54.302
0.4
3.8256
1.5
30.429
2.5
56.47
0.5
5.5492
1.6
33.254
2.6
58.467
0.6
7.4798
1.7
36.075
2.7
60.265
0.7
9.5911
1.8
38.875
2.8
61.838
0.8
11.859
1.9
41.639
2.9
63.152
0.9
14.262
2
44.349
3
63.4352
1
16.778

15