题目储油罐的变位识别与罐容表标定1.docx

1、题目 储油罐的变位识别与罐容表标定1题目 储油罐的变位识别与罐容表标定摘要 本文分别就题目提出的两个问题建立储油罐的椭圆柱体模型实际罐体模型，来解决储油罐的的变为识别和罐容标定问题。第一问中，用积分的方法得到正常和倾斜变位两种积分函数，通过matlab控制积分上下限得到一组罐容标定值，再用图像的吻合性来说明模型的可靠性。对于第二个问题，分正常罐体和横纵向变位罐体两种情况，分别用二重积分的办法得到实际罐体两边半球体和中间圆柱体中油体积与油高的函数关系式，从而得到总油体积和油高的函数关系，再利用附表二中的出油量和油高变化关系进行最小二乘法遍历求解，即得到a，b向量空间的最小二乘法中的最小值，求出变

2、位参数，再回代总油体积和油高函数关系式中，从而得到实际油罐的罐容定表。关键字 积分法 图像检验 最小乘法 罐容表标定 一、问题的重述通常加油站都有若干个储存燃油的地下储油罐，并且一般都有与之配套的“油位计量管理系统”，采用流量计和油位计来测量进/出油量与罐内油位高度等数据，通过预先标定的罐容表（即罐内油位高度与储油量的对应关系）进行实时计算，以得到罐内油位高度和储油量的变化情况。许多储油罐在使用一段时间后，由于地基变形等原因，使罐体的位置会发生纵向倾斜和横向偏转等变化（以下称为变位），从而导致罐容表发生改变。按照有关规定，需要定期对罐容表进行重新标定。图1是一种典型的储油罐尺寸及形状示意图，其

3、主体为圆柱体，两端为球冠体。图2是其罐体纵向倾斜变位的示意图，图3是罐体横向偏转变位的截面示意图。请你们用数学建模方法研究解决储油罐的变位识别与罐容表标定的问题。 （1）为了掌握罐体变位后对罐容表的影响，利用如图4的小椭圆型储油罐（两端平头的椭圆柱体），分别对罐体无变位和倾斜角为 =4.10的纵向变位两种情况做了实验，实验数据如附件1所示。请建立数学模型研究罐体变位后对罐容表的影响，并给出罐体变位后油位高度间隔为1cm的罐容表标定值。（2）对于图1所示的实际储油罐，试建立罐体变位后标定罐容表的数学模型，即罐内储油量与油位高度及变位参数（纵向倾斜角度 和横向偏转角度 ）之间的一般关系。请利用罐体

4、变位后在进/出油过程中的实际检测数据（附件2），根据你们所建立的数学模型确定变位参数，并给出罐体变位后油位高度间隔为10cm的罐容表标定值。进一步利用附件2中的实际检测数据来分析检验你们模型的正确性与方法的可靠性。二、问题的假设（1）附表中所给的数据真实可靠，统计时没有偶然误差和系统误差；（2）油罐与所给尺寸吻合，无变形和破坏现象；（3）不计罐内油位探针、进出油管对罐内油量体积的影响。三、问题的分析 对于第一个问题，要研究罐体倾斜变位对罐内油体体积的影响，就要分别建立油体积与油高度的函数关系，得到图像后通过对照分析变位后的影响，并利用附表一中倾斜变位出油量和高度的数值关系和罐容标定值画一个图中

5、，分析检验倾斜变位油体罐容标定表的可靠性。 对于第二歌问题，要研究变位后油体积和高度及变位参数之间的函数关系，并确定参数大小，利用附表二中的数据说明模型的可靠性就要先用积分求出不规则半球体和中间不规则半柱体体积之和与高度的函数关系，再利用附表二中出油量和高度对函数进行最小二乘法拟合。从而求出变位参数。在回代函数关系式求出实际储油罐变位后的罐容标定表。 四、符号说明表1： 符号说明表符号 符号含义h探针指示的油高 S(h)高度为h对应的椭圆（或圆）面积V(h）正常罐体高度为h时对应的油体体积L倾斜时h*cot(4.1)对应的长度其他符号说明：（1、其他符号在问题中说明。2、相同符号的含义由所在部

6、分说明决定）五、模型的建立与求解第一个问题：建立正常和倾斜变位体积和油高的函数关系 椭圆高度与面积坐标图 图（1）椭圆的方程为：；椭圆弦长AB=高度为h时对应的面积S(h）；则S(h)= 计算机运行结果 （一）、正常罐体油高为h时对应的体积V（h）； V（h）=S(h)*（0.4+2.05） 运用matlab求出每隔0.01m的标定值（附录1）。将此函数的理论值与附表一中累积进油量加上初始262L和油高关系图线对照，对照图如下： 程序见附件Tu2 正常椭圆体储油罐油量油高理论实际图线对照图 图（2） 显然正常油罐油量和高度关系的理论值和实际值绝大多数吻合得很好，只有在油高接近1200时才出现偏

7、差，而且偏差不是很大，图线的吻合性说明正常罐体油量与高度的函数关系是可靠的。（二）、倾斜变位4.1度时h与L的关系式为h=tan（4.1）*L带入函数S(h)得到 图（3）运用matlab软件通过for语句控制积分上限（附录2）从而得到一组罐容标定表，并画出其图线与附件一中倾斜变位进油量和高度的关系程序见附件Tu4 正常椭圆体储油罐油量油高理论实际图线对照图 图（4） 结果很明显理论值图线和实际值图线吻合地非常好，误差很小，说明倾斜变位模型可靠。表2： 椭圆柱体倾斜变位标定表 （数值为程序2的运行结果） 高度/m容积/L高度/m容 积/L高度/m容 积/L高度/m容 积/L0.011.6744

8、0.31595.24540.611885.13130.913228.61220.023.53110.32630.14600.621928.51320.923266.72190.036.26350.33665.58090.631971.93100.933304.42140.049.97460.34701.95840.642015.37200.943341.69120.0514.75610.35774.85780.652058.82400.953378.51130.0620.69100.36812.20320.662102.27520.963414.86140.0727.85410.37849.9

9、7480.672189.12540.973450.72000.0836.31630.38888.15380.682231.50010.983486.06420.0946.14230.39926.72180.692275.82440.993520.87040.1057.39330.40965.66060.702319.08631.003555.11440.1170.12690.411004.95400.712405.37241.013588.76930.1284.39670.421044.58400.722448.37181.023621.80780.13100.25410.431124.790

10、70.732491.25871.033654.20030.14117.74750.441165.33610.742576.64321.043685.91480.15136.92280.451206.15530.752619.11511.053716.91720.16157.81850.461247.23400.762661.42261.063747.17140.17180.25910.471288.55770.772703.55241.073776.63560.18204.9950.481330.11130.782745.49111.0838.5.26620.19228.90650.49137

11、1.88110.792787.22491.093259.81900.20254.88520.501413.85370.802828.73991.103776.63560.21281.85740.511456.01520.812870.02191.113805.26620.22309.76110.521498.35210.822911.82991.123833.01310.23338.53880.531540.85110.832911.05701.133859.81900.24368.14270.541583.49920.842951.82991.143885.61830.25398.52870

12、.551626.28320.852992.32641.153910.33150.26429.65680.561669.19020.863032.53071.163933.85850.27461.49080.571712.20770.873072.42741.173956.05500.28493.99680.581755.32320.883112.00041.183976.65550.29527.14380.591798.52390.893151.23381.193995.53730.30560.90270.601841.79730.903190.10991.204012.7448再将上述两图形

13、中的正常值和理论值在一个图中画出，来观察变位对罐体油量的影响 程序见附件Tu5 正常和变位的理论值对照图 图（5）由图（2）和图（4）说明这两条曲线是可靠的，通过对照发现差距比较大，说明倾斜变位时相同高度下油量的变化很大，倾斜变位对灌容影响很大。第二个问题：1、用二重积分法求解正常实际罐体油量和油高的函数关系 2、先用积分法求解油量和油高纵向变为角之间的函数关系，再将函数里的高度用显示高度与横向变位的关系式替换，从而得到油量与油高、纵向变位角、横向变为角之间的一个函数关系式，再用最小二乘法求出变位参数。(一)、正常实际罐体油量和油高的关系 先考虑弓形的面积与对应高度的函数关系式:再考虑体积:对

14、于中间圆柱体r=1.5;所以中间圆柱体部分高度为h时对应的体积： 弓形对应的高度为h的坐标图 图（6）先用二重积分的办法求两端不规则部分球体的体积：坐标图（7）通过已知数据容易求得球体半径R0=1.625m；运用matlab控制高度h得到一组标定值（附录3） 实际储油罐两端不规则球体的几何立体图形 图（7） 程序见附件Tu10 实际储油罐理论值图线和实际值图线对照 图（8） 由图像可知，实际值和理论值几乎完全重合，可见计算误差不大，（二）、横向纵向变位时油体积与油位高度的函数关系 总体思路：想考虑纵向变位油量与油高的函数关系，再用实际油高与显示油高、横向变位的关系式代换函数中的油高h，从而得到

15、油体积与显示油高，纵向变位参数，横向变为参数的函数关系式。 对于纵向变位，油位高度不同时，组成油体积的几何形状不同，如下图：由图知要分三种情况计算油体体积。 对于左右不规则球体体积：由于纵向变位角度不会很大，而且左右两边可以部分抵消，所以其体积的计算可以和正常罐体一样积分。但左右两边的油高不同. 左半球油高h1=h+2tan（a），h为油位探针位置实际油高； 右半球油高h2=h-6tan（a），下面分三种情况油体体积于油位高度的函数关系（式中h为探针位置实际油位高度，式中S(L)为图（6）中弓形的面积函数，L=h*tan（a） 实际储油罐倾斜变位图 图（9） 对于横向变位，图如下，要寻找显示油

16、高和实际油高、变位参数之间的关系式，由下图中间直角三角形可得： 横向变位时实际油位和显示油位高度之间的关系图 图（10） （三）、油体体积与油高、变位参数a，b，之间的函数关系将上式代人上面三个积分表达式得 再用matlab积分就可得到V1、V2、V3关于油高h，纵向变位角a，横向变位角b的函数表达式（程序见附录4），由于表达式很长，占用很大篇幅，在此不方便书写。再将V1、V2、V3保存为关于h，a，b的函数。再将V1、V2、V3保存为一个分段函数fv=fv(h,a,b)(程序如附录5)。（四）、 下面利用已建立的储油量与油位高度及变位参数的函数关系，并结合给出的实际数据，反过来对、进行求解。

17、因为实际储油罐的储油量初值未知，因此附件2中的D列所给储油量数据不准确，所以用非线性最小二乘法求参数，即在参数解空间中找到参数，使得最小。即。其中，表示了不同高度之间对应的理论储油量差，而表示了附件2中出油量值。在这个范围内以步长为0.01用遍历搜索算法求出、的局部最优解，求解过程通过Matlab编程实现，程序见附录（6）。最后得到变位参数的局部最优解为=，=。将、的值代入到前面建立的模型中，得到实际体积与显示高度的关系，通过这个关系建立出变位后的罐容表标定，如表3所示。表3：变位后的罐容表高度/实际体积/高度/实际体积/高度/实际体积/0.10.35341.119.20332.146.719

18、80.21.05241.221.87862.249.27760.32.19521.324.61262.351.73430.43.66541.427.38932.454.07070.55.38821.530.19292.556.26550.67.32151.633.00822.658.29530.79.43371.735.81992.760.13180.811.69911.838.61242.861.73980.914.09571.941.37042.963.06681.016.60332.044.07843.064.8804六、模型的检验问题一： 1、在无变位情况下：对给出的数据做三次拟合，得

19、到h=0:0.01:1.2对应的120组数据，和理论值对比分析，由图像可知误差很小，说明模型可靠2、在有变位情况下： 同样对进油的数据做三次拟合，得到h=0:0.01:1.2对应的120组数据，和理论值对比分析，问题二：1、在无变位情况下：给出的数据中，显示的油量容积不是真实值，而是假设没有变位情况下的值，我们可以根据这个值验证无变位情况的数学模型。2、在有变位情况下：、利用有变位下的数学模型和求得的，值，验证第二次出油时的相对误差。先求出对应每个显示高度的实际值V0，让后用上一项减去下一项，得到实际的相对误差，和注入油量作比较 ，发现误差很小。七、模型评价与优化 我们建立的数学模型是基于实际

20、图形的积分计算，但是对于第二问有位偏时候积分比较复杂，我们采用近似积分和数值积分的方法，存在一定的误差。考虑到外界环境已经内部因素如储油罐壁厚等因素，理论值和实际测量值难免有一定的误差，但是误差都在可接受的范围内。我们做出的罐容表的值都是对不同情况积分的理论值，没有添加调谐系数，这样更能直观的反应理论和实际测量之间的误差关系。为了更好的优化，还可以采用拉格朗日近似积分法，差值算法等方法，为了和给出的测量值吻合，也可以适当的在理论值前边加上一个系数。对于两个方向的位偏立体不规则图形的近似积分法，还有很多的研究空间。参考文献：1刘来福 杨淳 黄海洋等，数学建模方法与分析，北京，机械工业出版社，20

21、09.52张传义 包革军 张彪，工科数学分析（下册），北京，科学出版社，2001.93王文波，数学建模及其基础知识详解，武汉，武汉大学出版社，2006.54郑阿奇 曹弋，MATLAB实用教程，北京，电子工业出版社，2007.85姜健飞 胡良剑 唐俭，数值分析及其MATLAB实验，北京，科学出版社，2004.66蒲延炳，常见卧式油罐罐表，四川石油财经学院7高恩强（山东昌邑麻纺厂）丰培云（莱芜钢铁总厂供销处），卧式倾斜安装圆柱体油罐不同液面高度时贮油量的计算 附录：附录1：syms L h k=2*0.89/0.6*sqrt(2*0.6*h-h2);H=0:0.01:1.2;for i=1:121

22、 v(i)=vpa(1000*2.45*int(k,h,0,H(i),6); end附录2：syms L h k=2*0.89/0.6*sqrt(2*0.6*h-h2);s1=int(k,h,0,L*tand(4.1);h=0:0.01:1.2;for i=1:15v(i)=int(s1,L,0,0.4+h(i)*cotd(4.1);endfor i=16:117 v(i)=int(s1,L,0,0.4+h(i)*cotd(4.1)-int(s1,L,0,h(i)*cotd(4.1)-2.05);end for i=118:121 v(i)=pi*0.6*0.89*2.45-int(s1,L,

23、0,(1.2-h(i)*cotd(4.1)+2.05);endvpa(1000*v,8)附录3：syms y x hH=0:0.1:3;for i=1:31v(i)=vpa(2*int(int(sqrt(1.6252-x2-(y-1.5)2)-0.625,x,-sqrt(3*y-y2),sqrt(3*y-y2),y,0,H(i)+8*int(2*sqrt(3*h-h2),h,0,H(i),10);End附录4：syms h y x L a b hsL=int(2*sqrt(3*h-h2),h,0,tan(a)*L);v1=int(int(sqrt(1.6252-x2-(y-1.5)2)-0.6

24、25,x,-sqrt(3*y-y2),sqrt(3*y-y2),y,0,(h-1.5)*cos(b)+1.5+2*tan(a);v2=int(sL,L,0,2+(h-1.5)*cos(b)+1.5)/tan(a);v3=int(sL,L,0,(h-1.5)*cos(b)+1.5)/tan(a)-6);v4=int(int(sqrt(1.6252-x2-(y-1.5)2)-0.625,x,-sqrt(3*y-y2),sqrt(3*y-y2),y,0,(h-1.5)*cos(b)+1.5)-6*tan(a);v5=int(sL,L,0,6+(3-(h-1.5)*cos(b)+1.5)/tan(a)

25、;V1=simple(v1+v2);V2=simple(v1+v2-v3+v4);V3=simple(31/24*pi+8*pi*1.52-v5+v4);附录5：function fv=fv(h,a,b)if 1.5*(1-1/cos(b)=h=6*tan(a)./cos(b)+1.5*(1-1/cos(b) fv=VV1(h,a,b); end if 6*tan(a)/cos(b)+1.5*(1-1/cos(b)h=(3-2*tan(a)/cos(b)+1.5*(1-1/cos(b)fv=VV2(h,a,b);end if (3-2*tan(a)/cos(b)=h=3/cos(b)+1.5*

26、(1-1/cos(b)fv=VV3(h,a,b); end附录6：k=size(h);c=0,0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1.0,1.1,1.2,1.3,1.4,1.5,1.6,1.7,1.8,1.9,2.0,2.1,2.2,2.3,2.4,2.5,2.6,2.7,2.8,2.9,3.0,3.1,3.2,3.3,3.4,3.5,3.6,3.7,3.8,3.9,4.0;d=0,0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1.0,1.1,1.2,1.3,1.4,1.5,1.6,1.7,1.8,1.9,2.0,2.1,2.2,2

27、.3,2.4,2.5,2.6,2.7,2.8,2.9,3.0,3.1,3.2,3.3,3.4,3.5,3.6,3.7,3.8,3.9,4.0; e=0;f=0;y1(1,1)=0;for j=1:41 for m=1:41 for i=1:k v(i)=f3(h(i),c(j),d(m); end sum=0; for n=2:k v1(n-1)=v(n-1)-v(n); sum=sum+(v1(n-1)-v0(n-1)2; end y1(j,m)=sum; endend smin=y1(1,1);for j=1:41 for m=1:41 if (y1(j,m)-smin)0 smin=y1

28、(j,m); e=j; f=m; end end endef保存函数下列函数为m文件function f3=f3(h,c,d)a=c*pi/180;b=d*pi/180; h1=1.5-(1.5-h)*cos(b);h2=(1.5+h)*cos(b)-1.5; if 0h=6*tan(a) y=0:0.00001:(h1/tan(a)+2); m=1.5-h1+(y-2)*tan(a); f=2.25.*acos(m./1.5)-m.*sqrt(2.25-m.2); t1=trapz(y,f); y1=0:0.00001:(h1+2*tan(a); f1=(0.390625-y1.2+3.*y1).*asin(sqrt(3*y1-y1.2)./(0.390625+3.*y1-y1.2)-0.625.*sqrt(3.*y1-y1.2); t2=trapz(y1,f1); f3=t1+t2;endif 6*tan(a)h=1.5 y=0:0.0001:8; m=1.5-h1+(y-2)*tan(a); f=2.25.*acos(m./1.5)-m.*sqrt(2.25-m.2); t3=trapz(y,f); y1=0:0.0001:(h1+2*tan(