数模等奖A题储油罐变位识别与罐容表标定

1、储油罐的变位识别与罐容表的标定摘要本文研究储油罐的变位识别与罐容表的标定。分别以小椭圆型油罐和实际卧式储油罐为研究对象，运用高等数学的积分的知识，分别建立罐体变位前后罐内油体积与油高读数之间的积分模型，使用Matlab软件得出结论。对于问题一，以小椭圆型储油罐为研究对象，在无变位时，小椭圆型储油罐为规则的椭球柱体，可利用解析几何与高等数学的知识建立油罐内体积与油高读数之间的积分模型，得出罐体无变位时的理论值。当罐体发生纵向变位时，小椭圆型储油罐的截面不再是规则的几何形体，但根据倾角a及所给小椭圆型罐体的尺寸，可得其截面面积的表达式，利用高等数学中积分的方法，根据不同油高，建立了模型一，得到了储

2、油量和油高的关系公式。最后，根据实验数据的处理，用拟合的方法，修正了某些系统误差的影响，计算出罐体变位后油位高度间隔1cm的罐容表的标定值。对于问题二，由于实际储油罐内没油的高度不同，我们将其分为五种情况分别讨论，并对每种情况建立积分公式，得出罐内油体积与油位高度及变位参数（纵向倾斜角a和横向偏转角B）之间的函数关系式，利用所给的实验数据，运用最小二乘法，建立非线性规划模型nminerrorC二QilData）八（vcj-vc,Hj-oiiDataj2argaP2用Matlab非线性规划求解得出使得总体误差最小的a与B值：a=2.12，B=4.06。通过a与B的数值计算出出油量理论值与实测值的

3、平均相对误差小于0.5%。对模型进行了较为充分的正确性验证和稳定性验证：在a与B的值为0时，其计算出来的罐容值与理论值完全吻合，说明模型在体积计算上是正确的；当对油高进行0.1%的扰动时，a的值变化也在0.1%左右，说明a的稳定性很好，但是B的值从4.06。变成了3.75。，变化了大约8%所以我们详细分析了B的数学表达式，从理论上分析了影响其稳定性的因素。根据得到的变位参数计算出实际罐体变位后油位高度间隔为10cm的罐容表的标定值。最后，本文对模型的优缺点进行了评价，并讨论模型的推广。关键字:储油罐；变位识别；罐容表标定；非线性规划问题重述通常加油站都有若干个储存燃油的地下储油罐，并且一般都有

4、与之配套的“油位计量管理系统”，采用流量计和油位计来测量进/出油量与罐内油位高度等数据，通过预先标定的罐容表（即罐内油位高度与储油量的对应关系）进行实时计算，以得到罐内油位高度和储油量的变化情况。许多储油罐在使用一段时间后，由于地基变形等原因，使罐体的位置会发生纵向倾斜和横向偏转等变化（以下称为变位），从而导致罐容表发生改变。按照有关规定，需要定期对罐容表进行重新标定。根据上述所述，求解下列问题：（1）为了掌握罐体变位后对罐容表的影响，利用小椭圆型储油罐（两端平头的椭圆柱体），分别对罐体无变位和倾斜角为：=4.1的纵向变位两种情况做了实验。请建立数学模型研究罐体变位后对罐容表的影响，并给出罐体

5、变位后油位高度间隔为1cnm勺罐容表标定值。（2）对于实际储油罐，试建立罐体变位后标定罐容表的数学模型，即罐内储油量与油位高度及变位参数（纵向倾斜角度和横向偏转角度一：）之间的一般关系。请利用罐体变位后在进/出油过程中的实际检测数据，根据你们所建立的数学模型确定变位参数，并给出罐体变位后油位高度间隔为10cnm勺罐容表标定值。进一步利用实际检测数据来分析检验你们模型的正确性与方法的可靠性。问题分析本文研究罐容表的读数与储油罐的变位的关系。借助高等数学积分的方法，求出储油量与油高读数的函数关系式，并对倾斜的储油罐进行容量标定。对问题一的分析问题一中用小椭圆储油罐分别对罐体无变位和纵向倾斜进行实验

6、，研究变位对罐容表的影响，因此我们分别建立变位前和变位后的罐容表读数与罐内油体积的函数关系式，通过函数关系式计算出理论值，再与所给的实际值相比较，得出其相对误差，然后通过分析系统误差进行修正，出罐体变位后油位高度间隔为1cm的罐容表的标定值。对问题二的分析问题二中是以实际储油罐为研究对象，不仅考虑了储油罐的纵向倾斜，而且还考虑了横向偏转，为了使问题简化，我们先只考虑纵向倾斜，由于储油罐的形体不规则，所以我们将它分成如图1所示的三部分，分别算出每部分的体积与罐容表读数的函数关系式，然后对其求和。再考虑横向偏转，建立它与所给的油高的函数关系式。然后将二者进行综合考虑得出变位后罐容表读数与储油罐内油

7、体积的函数关系式，通过关系式和所给数据，运用最小二乘法，通过MATLAB?序,搜索出a和B的最小误差解，再对模型的稳定性和正确性进行评定，最后给出高度间隔10cm的罐容表的标定值。模型假设假设一：数据是储油罐的内壁参数。假设二：忽略温度、压力对汽油的密度的影响。假设三：储油罐在偏移的过程中，油位探针始终与油罐底面垂直。假设四：对卧式储油罐来说，不考虑其长期埋在地下所发生的蠕变假设五：累加进出油量数据是准确可靠的。符号说明H:对应于罐容表读数的液面实际高度。Hi：球冠中与油罐圆柱左侧底面距离为x处的油高球冠中与油罐左侧底面相距为x处的小圆半径。H2:球冠中与油罐圆柱右侧底面距离为x处的油高F3:

8、球冠中与油罐右侧底面相距为x处的小圆半径。R：储油罐圆柱部分的底面半径。R:球冠所在球体的大圆半径。Hi：第i条数据所对应的罐容表读数。OilData:用于分析的油量进出数据。a:椭圆长半轴长。b:椭圆短半轴长。n:用于分析的进出油测量数据个数。ho:罐容表读数。模型的建立与求解5.1模型一的建立与求解问题一要求研究罐体变位后对罐容表的影响，并给出罐体变位后油位高度间隔为1cm的罐容表标定值。计算未变位和变位的理论罐内油位高度与储油量的关系利用高等数学中微元法求体积的方法建立罐容表读数与罐内油体积的函数关系式的模型（1）在无变位的情况下，储油罐内的油所占空间为柱体，其体积为v=SLL其中S为柱

9、体底面面积，L为柱体的长度。hSb2XJdy底面椭圆方程为22a2b2a2b2将（4）代入（2），得到其积分解析表达式为aI222h12S(h、b-hbarcsin-b)bb其中，h=Hb图2微元法求椭圆切面面积图2微元法求椭圆切面面积（1）如图TOCo1-5hza2H12S(H-b).、H(2b-H)barcsin(-1)二bbba12H12(9)(9)VL(Hb).H(2b一H)b2arcsin(1)b2bb2图3油罐无倾斜时示意图(2)当油罐发生纵向偏转时，油罐中油所占空间为一倾斜柱体，如图4所示:h=H(X0.4)tan：(10)又根据油面的高度不同，可分为以下三种情况：若油面位于图5

10、所示位置，则:TOCo1-5hz(H+OAtan/tanaar22(11)Vi(H0.4tan：-xtan.-b)b-(H0.4tan；:-b),2.H0.4tan：-xtan：-b1,2,barcsinbdxb2若油面位于图6所示位置，则:2.45ab(H0.4tan；.；-xtan-b).b2-(H0.4tan：-b)2-b2arcsinH0.4tan：-xtan二一b12bdxb2(12)图7情况3高油位若油面位于图7位置，贝二二abL-2.45ay220心2山)伽:.b(xtan：b).b-(xtan：b)2xtan；.-b12barcsinbdxb2由上述公式知，油罐的变位会对罐内油

11、高与储油量的对应关系(罐容表)，产生较大的影响。综合式(11)-(13)，可以得到模型1如下：(H+0.4tan-)/tan-已0b(H0-4tan:-xtan：-b).b2_(H20.4tan：-b)b2arcsinH0.4tan；:-xtanb1,2,二bdx，b2当0H2.05*tan:a(H0.4tan：-xtan；:-b).b2(H0.4tan：-b)2b【0V(H)=5十b2arcsinM+4b2dx,(当2.05*tanH.2-0.4ta切b2a_b),b2_(xtan_b)2：abL_0心心/tan/(xtan：b2arcsinxtanJb丄二b2dx，2(当a-xtan:)R

12、_(R-H-2tan：亠xtan_：:)(22).(H2tan：-xtan：)(2R-H-2tant*xtan：)dx对于第一部分体积可用图17求解油高：H=GT图中：R=OAR二AlR2=MNCI=1R2_(R_1)2=r2二R=1.625R22=R2-(R-1x)2二R2R2-(R-1x)2HMT-MGMG-BQ-PQBQ=H2tan:PQ二xtan:-(23)所以,巴=R2-(R-H-2tan：-xtan：)然后再利用下图确定积分区间0,xi:则：dvi=R22cosjR2-山_(R2_H)H12R2-H1dxR2图18图中CM即为所要求解的x1。CI=1OI=R-1IQ=R-H2tan

13、：KI=(R-1)tan：，”KN=KGCosa=(R-1)tana+R-(H+2tana)cos。=OPPF=R2-OP2EQ=PE-PN+NQNQ=KQn:OK=PN=!cos。联立以上各式可得：x=R2-R-1tan:亠R-H2tanJcos2-_-cos：+(R-1)tana+R-(H+2tana)sinacosa(25)(25)V=oR22cosj电土!_(R2_H),H12R2-H1dxR2对于第三部分的体积，方法与第一部分的体积求解相似。其体积为:式中:x2=、R2-|R-R-1tan：-H6tancos1cos-|R-R-1tan：-H-6tan：|sin：cos：R3-R2-

14、R-1x2H2二R3-Ri-H6tan二xtan：所以:二:R/cosR-H)已2只2-已dxR282R-(H2tan：-xtan：)亠1Rcos一(R-H-2tan一xtan二)Ri,(H2tanJxta门一：)(2民-H-2tan1xtan-：)dx0!|R32COSR3_H2_侃十2血网_H2JdxILR3X20R3其他几种情况均可用类似的方法进行求解。对于图13所示情况：V=0，&2込，R2一已一(&HJH12R-H1dxR2H2tan：tan:RicosR,-(H2tan：-xtan：)R(26)(27)(28)(29)(30)(R-H-2tan：xtan：)、(H2tan：-xta

15、n：)(2R_H_2tan：xtan-)dx(R-H-2tan：xtan：)、(H2tan：-xtan：)(2R_H_2tan：xtan-)dx(31)对于图14所示的情况:V=0R22对于图14所示的情况:V=0R22x3JRHiR2+丿iH+2tanaxtanaR黛一cos-二-cos(HiR2).Hi2R2-H,dx十(H+2tanaxtana-R,).(H2tan：-xtan：JQRUH匸2tanx，xtan：)82iR-H2tan.：-xtan：亠1尺cosR-H-2tan：亠xtan：x3R,H2tanxtan二2R-H-2tantxtan：dx:Ricos只3口出_出.H22R-

16、H2dxR31二R22dxXi积分限中出现的表达式H+2tanaRX3*1tana(33)对于图15所示情况,川，2V0R22-4蔥-cosHi-R2+R丿82JR1H+2tan。-xtan。:-cos-R,(Hi-R2),Hi2R2二Hdx-R、+(H+2tan。-xtan-R)x,Rs2dx(34)R(R-H-2tan：xtan：).(H2tan：-xtan：)(2R-H-2tan：xtan：)dx(35)对于图16所示情况:V=2；二R2dx8-RiPR，*f,x1-：R22cos4民_Hi_H)Hi2R-HidxR2七曲R2COSRTH2taFan：)(36)(37)(38)2_1fv

17、ol6(x)二&(蔥一cos罟)HT”2厂叮(41)因此，我们总结出7个分段被积函数如下：fvol1(x)=R2cos21_(R2-HJ,Hj2R2-已R2gRW严严fVOl3(X)二R32COSR3_H2_R3_H2、H22R3-出R3/H_RIfvol4(x)=&2兀-cos12+(巴-&)已(2&-出)(39)R2丿，“JH+2tana-xtana-R)_fvol(x)=R兀cos+(H+2tanaxtanaRJR丿.(H2tan：-xtan：)(2R,-H-2tan=xtan：)(40)(42)2fvol7(x)=:R?F面再对横向偏移角B进行分析研究，如图19所示:图中ho为罐容表的

18、读数，所以真实液面高度为：H=Rcos：（h0一R）cos：R一R1cos-即：Hf_R）cos一：R（43）所以最终所得的体积关系式只需将上述体积关系式中的H换为式（30）,即可得出。经以上分析，我们得到模型二如下：V（G,ho）才H2tan:-fVOh(X)dX亠|tan：fvol2(x)dX,$08xfvoh(x)dx亠ifvol2xdxofvol3xdx,1fvol4xdxfvol7X1当6tan：_H_1.5-3tan：8dxfvol2xdxfvol3(x)dx,LX3L0当1.5-3tan：HX212fvol5xdx亠ifvol6xdx亠iR3dx,101x?当1.57tan：_H

19、_3-2tanH亠2tan：Iafvol?(x)dxX3dxofvol5xX2X3_1.57tan18fvol4(x)dx+Xfvol7(x)dx+J02fvol7xdx&R2_、ol1(x)dx：当3-2tan：一H乞3(44)nminerror（：,IQilData）（VChj-vc：屮2）-08込）2argap7（45）式（44）中H由式（43）确定，其他变量的定义参见符号说明：H=（hR）cosP+R使得式（45）取得最小值的/即是待定的变位参数模型的求解我们利用MATLAB软件，将模型二中的式（44）和式（45）编程（相关程序见附录），利用非线性规划的方法求解a和B值，a=2.12,

20、B=4.06。根据得到的变位参数值，得到罐体变位后油位间隔为10cm的罐容表标定值。罐容表标定值H(mm)V(L)H(mm)V(L)046.40662150030232.31100354.5192160033046.072001061.754170035855.973002214.65180038646.714003691.607190041402.875005419.525200044108.696007356.6210046747.987009471.698220049303.8780011739.27230051758.569001413797100016645

21、.525005628685260058314.42120021920.95270060149.79130024654.29280061756.69140027429.93290063083模型的分析6.1模型的正确性如图所示，横坐标为题目所作实验的出油量，竖坐标为经过所建模型求解的两次体积差，即模型的出油量，二者相对误差较小，线性拟合较好。6.2模型的灵敏度分析当对油高进行0.1%的随机扰动,当对油高进行0.1%的随机扰动,即让H=H10.001*rand-1,1时，模型求得a=2.119。，比较无扰动的a=2.12。，变化幅度也在0.1%左右，说明a的稳定性较好，但

22、是B的值从4.06。变成了3.75。，变化了大约8%所以模型在B方向上抗干扰能力较低，因此我们对B进行理论分析。由H=(h0-R)cos一：R可知：iH-R1二cos-所以当h。在R附近变化时对B的影响较大。用Matlab对其进行数据模拟,结果如下图所示：表明B在Ri附近比较敏感7.1模型的优缺点七.模型的评价与推广模型的优点本文借助高等数学微积分的思想，建立罐体在变位前后标定罐容表的数学模型，得出罐内储油量与油位高度及变位参数的函数关系式，理论基础成熟，可信度较高。该模型以微积分为基础，简单易懂，又有相应的软件(Matlab软件)支持，算法简单，容易推广。运用多种拟合方法使结果更加精确，通过

23、灵敏度和误差分析使模型更具有实际意义，增加了应用价值。模型二的误差分析运用最小二乘法，使得模型得出的结果更加准确。模型的缺点模型二中的B值的稳定性不是很好在用拟合法处理数据的时候，由于模型假设具有一定的主观性，导致拟合的曲线不是十分精确。7.2模型的推广本文所建立的微分模型不仅适用于储油罐的顶部为球冠的情况，还可以推广到顶部为弧形顶、平顶、椭球顶、锥顶的情况，易于计算，在实际应用中具有延伸和推广的价值。八.模型的改进本文所建立的积分模型中a与B是独立无相关的，而在模型二的求解中用B的值调整了油高读数，然后再去求解a,因此使得这两个方向的变化独立了，影响了模型的精度。因此我们必须考虑a方向上对油

24、高读数的影响。不妨假设H=(h0cosa-R)cosB+R,，则模型二的形式如下：V(a,B，ho)“argminapH：2tanfvoh(x)dxotan：fvol2(x)dx,当0乞H8xfvoh(x)dx亠Ifvol2xdx亠ifvol3xdx,当6tan：-H_1.5-3tan:1X38fvol4xdx%fvol7xdx0fvol5xdx电fvol2fvol2(x)dx.6tan:X2电_、x)dx+J。fvol3(x)dx,当1.5-3tan：_H_1.57tanr8X212fvol4xdx亠Ifvol7xdx亠Ifvol5xdx亠ifvol6xdx亠！R3dx,v00 x?当1.5

25、7tan：_H乞3-2tan：fvol2(x)dx12X20fvob(x)dx十8応Rifvok(x)dx+H亠2tatana0当3-2tan：_H乞3(47)式中:H=(hocos：-RJcos：Rnminerror(：,QilData)八(v(：,hj-vg,：屮匸)-心氏內)2arg:,-i4用matlab软件搜索求得：a=2.0925B=3.5055误差的平方和为2.4574，较前个模型二的误差平方和2.4570大了一点点，但这微小的差别几乎可以忽略不计，但是我们得到的a与B更小，从实际意义来说，更加符合实际，我们认为这个两个值更加准确。当然还有很多的假设，他们都满足一定的数学条件，有

26、自己的物理意义。如可以假设H=ho(1tana)RcosP+R，或者H=(h0RJcosBcos十R等等。但是a与B真正意义上的关系式还需要很复杂的数学推导和耐心的去解决，在此由于时间有限，我们不予解答出来。等以后有时间再仔细推敲。模型二数据的误差检校，由于模型二以实际卧式储油罐为研究对象，则不能忽略注油管和出油罐等管道对罐内油体积的影响。因此我们必须对模型二的数据进行检校，排除系统误差的影响。罐内油体积的相对误差为17%左右，用油高读数和实际油的体积与理论算出来的体积之差拟合，如下图所示：2500200030001000其相对误差大约为0.6%,较前面的相对误差减小了100多倍，数据更加精确

27、使模型更加优化。九.参考文献【1】姜启源，数学模型（第三版）高等教育出版社，2003【2】石辛民，基于MATLAB勺实用数值计算.清华大学出版社，2006【3】黄永建，MATLAB!言在运筹学中的应用，湖南大学出版社，2005附录附表无变位进油无变位进油实际油量修正后油量相对误差实际油量修正后油量相对误差312311.96390.00011617621761.7960.000116362361.96439.87E-0518121811.7810.000121412411.94680.00012918621861.7920.000112462461.96318E-0519121911.7670.

28、000122512511.93420.00012919621961.7720.000116562561.93790.0001120122011.7850.000107612611.93430.00010720622061.7890.000102662661.91410.0001321122111.7620.000113712711.92929.94E-0521622161.7710.000106762761.91080.00011722122211.7510.000112812811.91190.00010922622261.7260.000121862861.890.000128231223

29、11.7590.000104912911.88630.0001252315.832315.560.000116962961.87660.0001282365.832365.5760.00010710121011.8810.0001172367.062366.7840.00011710621061.8890.0001052417.062416.8070.00010511121111.8929.72E-052467.062466.7780.00011411621161.8889.63E-052517.062516.7960.00010512121211.8780.0001012567.062566

30、.7550.00011912621261.8660.0001062617.062616.7510.00011813121311.860.0001072666.982666.6760.00011413621361.8270.0001272668.832668.5410.00010814121411.8220.0001262718.832718.5160.00011514621461.8150.00012615121511.8240.00011715621561.8220.00011416121611.8290.00010616621661.820.00010817121711.8160.0001

31、07附表二变位进油变位进油变位进油实际油量修正后油量相对误差实际油量修正后油量相对误差实际油量修正后油量相对误差962.86962.00467910.0008891912.731913.4811240.0003932862.732861.2431860.0005196391012.861007.2547970.0055651962.731966.211970.0017712912.732913.1836870.0001557361062.861063.4190020.0005262012.732016.142670.0016932962.732966.4459250.0012526521112

32、.861110.1227830.0024662062.732060.9123560.0008823012.733013.4231760.0002300291162.861161.827970.0008882112.732114.3150550.000753062.733066.5006350.0012296211212.861216.0182890.0025972162.732164.4988790.0008173112.733117.344870.0014803851262.861262.0214390.0006642212.732215.3194070.0011693162.733163.

33、5720920.0002661841312.791314.680580.0014382262.732263.6933880.0004263212.733215.9466540.001000221362.791363.5369920.0005482312.732314.3049570.0006813262.733268.9423470.0019004151412.731412.8611729.28E-052362.732363.941710.0005133312.733315.7196880.0009016711462.731462.097550.0004332412.732409.830290

34、.0012033362.733362.7681831.13545E-051512.731515.9647940.0021342462.732460.9545680.0007213412.733415.4184920.0007871631562.731567.5356790.0030662512.732512.3336690.0001583462.733462.7728041.23611E-051612.731618.2414780.0034062562.732557.991510.0018523512.733518.1828140.0015498951662.731669.3006630.00

35、39362612.732610.2927070.0009343514.743522.2274290.213%1712.731714.8669640.0012462662.732663.7077850.0003671762.731765.8661320.0017762712.732714.6784420.0007181812.731816.0458360.0018262762.732764.0505670.0004781862.731867.7639010.0026952812.732813.8361240.000393MATLAB程序functionvolumes=totalEval(x)ti

36、c;oilTank=getOil;globalHalphaalpha=x(1);beta=x(2);n=size(oilTank,1);volumes=zeros(n,1);fori=1:31HO=(i-1)/10;H=(H0-1.5)*cos(beta)+1.5;x1x2=calcL1(H,alpha);ifH6*tan(alpha)volumes(i)=real(quad(intV2,0,x1)+quad(intV1,0,(H+2*tan(alpha)/tan(alpha);elseifH1.5-3*tan(alpha)volumes(i)=real(quad(intV1,0,8)+qua

37、d(intV2,0,x1)+quad(intV3,0,x2);elseifH1.5+7*tan(alpha)lmax=(H+2*tan(alpha)-1.5)/tan(alpha);volumes(i)=real(quad(intV4,0,1)+quad(intV5,0,lmax)+quad(intV1,lmax,8)+quad(intV3,0,x2);elseifH3-2*tan(alpha)volumes(i)=real(quad(intV4,0,1)+quad(intV1,0,8)+quad(intV6,0,1);elsevolumes(i)=real(quad(intV7,0,1)+8

38、*pi*2.25-real(quad(intV2,0,x1)+quad(intV1,0,(H+2*tan(alpha)/tan(alpha);endendtoc;functionxmax1xmax2=calcL1(H,alpha)R=1.625;R1=1.5;talpha=tan(alpha);2*cos(alpha)2*cos(alpha)*cos(alpha);Rc=(R-1)./cos(alpha);xmaxl=(RA2-(R-1)*talpha+R1-(H+2*talpha)A2).a0.5-Rc+.(R-1)*talpha+R1-(H+2*talpha)*sin(alpha)*cos

39、(alpha);xmax2=(RA2-(-(R-1)*talpha+R1-(H-6*talpha).A2*cos(alpha)A2).A0.5-Rc-.(-(R-1)*talpha+R1-(H-6*talpha)*sin(alpha)%functionV1=intV1(H,x,alpha)functionV1=intV1(x)R1=1.5;globalHalpha%H=1.29;%alpha=0.07;HH=H+2*tan(alpha);V1=R1A2.*acos(R1-(HH-x.*tan(alpha)/R1)-(R1-HH+x.*tan(alpha).*(HH-x.*tan(alpha).*(2*R1-HH+x.*tan(alpha).A0.5;functionV=intV2(x)globalHalphaR=1.