储油罐的变位识别与罐容表标定

1、 WORD 储油罐的变位识别与罐容表标定薛申芳（学院数学系， ：054001）摘要 该文就储油罐的变位识别与罐容表标定问题1，用四次多项式灌油体离散变率进行回归得到连续变化率，建立了连续型微分方程模型，通过该微分方程初值问题，得到倾斜灌油体体积随高度的变化规律，然后把无变位与倾斜情况的灌油体随测量高度的变化引起的灌油量差异，以与当灌油体一样时，而油标高度测量值差差异来判断分析影响。该文就问题2的倾斜旋转油罐，通过坐标系旋转方法，建立了计算灌油体与油位关系的精确积分模型；由于灌油体不是油位、倾斜角以与旋转角的初等函数，利用所给测量数据和所建立的数学模型进行变位参数的确定时采用了最小二乘方法，且对

2、部分被积函数的二元二次多项式展开近似处理，以与对倾斜角和旋转角的离散化（把倾斜角和旋转角的围化为一些离散小区间）搜索，计算结果显示倾斜角为，旋转角为；最后利用所给出的模型和所求出的倾斜角和旋转角给出灌容标定值。关键词 数学建模；罐容表标定；多项式回归；坐标旋转变换；积分1 问题 通常加油站都有若干个储存燃油的地下储油罐，并且一般都有与之配套的“油位计量管理系统”，采用流量计和油位计来测量进/出油量与罐油位高度等数据，通过预先标定的罐容表（即罐油位高度与储油量的对应关系）进行实时计算，以得到罐油位高度和储油量的变化情况。许多储油罐在使用一段时间后，由于地基变形等原因，使罐体的位置会发生纵向倾斜和

3、横向偏转等变化（以下称为变位），从而导致罐容表发生改变。按照有关规定，需要定期对罐容表进行重新标定。现用数学建模方法研究解决储油罐的变位识别与罐容表标定的问题。 （1）对小椭圆型储油罐（两端平头的椭圆柱体），分别对罐体无变位和倾斜角为a=4.10的纵向变位两种情况做了实验(实验数据参看CMCM2010A附件1),去建立数学模型研究罐体变位后对罐容表的影响，并给出罐体变位后油位高度间隔为1cm的罐容表标定值。（2）对中间部分为柱面，两端为球冠面的实际储油罐，建立罐体变位后标定罐容表的数学模型，即罐储油量与油位高度与变位参数（纵向倾斜角度a和横向偏转角度b）之间的一般关系、利用罐体变位后在进/出油

4、过程中的实际检测数据（参看CMCM2010A附件2）确定变位参数，并给出罐体变位后油位高度间隔为10cm的罐容表标定值。2 问题12.1 问题1分析图2 出油变化率图1 进油变化率从附表1中，可以得到油位以与油体的改变量，从而可以建立微分方程模型。图1和图2分别给出了无变位时进油和出油情况的罐油体增加量（出油时罐油体增加量取负）与油位增量之比（油体变化率）随油位变的离散图形。在图1中有几个不正常的点，可能由于测量误差造成的；而从图2中的点比较正常，所以在无变位情况建模时采用附表1中的出油数据。对倾斜情况，通过对附件1中的数据分析，采用了附件1中进油测量数据，且在计算变化率时把一组点：累计进油量

5、为0.0017l（该值由积分计算得到）时，油位高度为0。上述两种情况，通过对灌油体的离散变化率进行拟合，得到连续变化率，进而可以求解连续变量的微分方程去得到灌油体体积随高度的变化规律。2.2 问题1变量说明记：罐油体的体积（）；: 罐油位（）；：油位增量（）；: 罐油体的增加量（）；：，即为油体关于的变化率的拟合函数；其中，为1时表示无变位情况，为2时表示倾斜情况。3 问题1数学模型与其与解下面通过建立模型，得到无变位和倾斜两种情况下罐油体的体积与油位高度的关系。3.1 无变位情况由2.1的分析以与2.2定义的变量知，现选用出油数据对随着油位高度变化的数据，用四次多项式进行回归得1： （1）这

6、里显著性水平取为，其决定系数为,，可知回归方程非常显著，回归模型成立。通过解微分方程初值问题： （2） 得到，其解为 （3）上式即为无变位情况下罐油体的体积与油位高度的关系。3.2 倾斜情况为了让结果更符合实际，利用附件1中倾斜进油变位数据，对随着油位高度变化的数据用四次多项式进行回归得： （4）显著性水平也取为，决定系数为,，回归模型成立（经验证比用出油数据拟合效果好）。通过解微分方程初值问题： （5） 得到，其解为 （6）无变位与倾斜情况的灌油体随测量高度的变化以与灌油量差参看图3。当灌油体一样时，而油标高度测量值差差异的部分数据参看表1。倾斜灌容标定值参看附表1。图3 无变位与倾斜比较油

7、量差表1 一样油体对应的油位V(m3)0.00170.50171.00171.50172.00172.50173.00173.5017h1(m)0.00150.21900.35940.48430.60390.72350.84790.9870h2 (m)00.18540.32700.45820.58220.70.81440.9333注：上表中为无变位油位，为倾斜油位。4问题2下面建立中间部分为柱面，两端为球冠面的实际储油罐变位后标定罐容表的数学模型。4.1问题2分析油罐变位包括倾斜和旋转两种变位方式，这里先对倾斜进行讨论，然后再加上旋转。首先建立适当的坐标系，对于倾斜变位，相当于实行坐标系旋转，

8、在新的坐标系下，计算罐体中的油体与坐标高度的关系；对于旋转变位，根据罐体的轴对称性，便可以得到油位的坐标高度与测量浮标高度的关系，从而得到灌油体与测油浮标高度的关系的数学模型。然后利用最小二乘法得到参数和。4.2问题2变量记号记：油罐纵向倾斜角(弧度),假设不超过10o ×/180)，且假定油罐右端上斜；：油罐旋转角（弧度，假设不超过10o ×/180）；：油罐柱面半径（1.5）；：油罐柱体的长度（8）；：油罐油面铅垂坐标()；：变位后罐油浮子的高度,即浸泡在罐油中的浮杆长度（m）；V：罐油体积()；：球冠半径(1.625)。5问题2模型现在建立问题2的数学模型。5.1 无

9、变位油罐的坐标系对无变位油罐而言，以柱面左侧圆的圆心为坐标中心，指向正上方，轴穿过柱面中心轴指向右方，建立右手坐标系（参看图4）.图4 罐体无变位坐标系 记分别为平面与罐体柱面所交矩形的四个顶点（参看图4）。则在坐标系下，油罐柱面方程为：(7)左、右球冠面满足的方程分别为： (8)和 (9)柱面左、右侧所在的平面方程分别为： (10)和 (11)5.2 倾斜油罐的坐标系假设油罐右端向上倾斜，倾斜角度为，这等价于将坐标系绕轴顺时针旋转角度，假设得到的新坐标系为,则这两个坐标系之关系为2：或 (12)其中为旋转矩阵（正交矩阵），且 (13)这里旋转矩阵中的取旋转角度的负值（绕轴逆时针旋转取正，顺时

10、旋转针取负）。于是在坐标系下，油罐柱面方程为： (14)左、右球冠面满足的方程分别为： (15)和 (16)柱面左、右侧所在的平面方程分别为： (17)和 (18)记分别为点在坐标系下的坐标，则 (19)且在假定最大倾斜角度不会超过10º的情况下，易得： (20)5.3倾斜油罐灌容模型下面在坐标系下讨论倾斜油罐油体与坐标高度和倾斜角的关系： (21)利用微元法知，在方向取微元，用过轴上的点作垂直于轴的平面去截取油罐体，得到的截面面积与有关，记之为,则微元体的体积为，则罐油体积可用积分 (22)得到。由于罐面是的分片函数，故对的不同区间围进行分别积分，即分为三个区间围分别进行积分3，即

11、 (23)其中为当在不同区间上用过轴上的点作垂直于轴的平面去截取油罐体，得到的截面面积。分别表示用过轴上的点作垂直于轴的平面去截取油罐体左、右球冠部分所得到的截面面积。它们的表达式分别为： (24) (25) (26) (25) (26) 其中将和 、的表达式代入(23)式便得倾斜油罐油体与坐标高度和倾斜角的关系的数学模型。5.4倾斜旋转油罐灌容模型下面考虑加上旋转发生（旋转角度为，参看图9），为罐体倾斜和旋转后的斜，且记：油罐最低点到油面的高度（参看图5）；：油罐只发生倾斜时，游标测得的高度（参看图5和图6）；图5 油罐倾斜后坐标图6 油罐旋转断面图由几何关系以与所建的坐标系，得：由于假定油

12、罐右端往上倾斜，取负值，从而得到与的关系为 (27)将(27)式代入(23)式，即得到了罐油体积与、与的关系的一般模型： (28)至此，就建立了倾斜旋转油罐灌容的精确数学模型。6 问题2变位参数确定与灌容标定（1）首先讨论、的确定。利用附件2中E列（显示油高）和F列（显示油量容积）的的数据，即由测量值与得： (29)为确定和，理论上讲课采用最小二乘法，令误差为，即通过求极值问题： (30)s.t. (31)去得出和。遇到问题是并不是的初等函数。故采用离散方法和MATLAB软件去解决。方法是把的区间-10/180,0以与的区间0, 10/180均分成一些离散个小区间（每个小区间长度适当小），对每

13、个去得到，去求 (32)由于 、的积分表示不是初等函数，故采用近似处理，即把 、积分表达式中的被积函数对、分别在（0，-0.6）、（0，8.32）进行二元二次多项式展开（带来的体积误差），然后再进行积分求得 、的近似表示式，由于其解析式繁杂，用软件进行处理，不再写出相应的展开表达式。计算结果为（参看程序附录2Vprog2_1）。根据所得角度，把部分测量值相应的理论值进行了比较，绝对误差（参看图7）。图7 结果比较（2）罐容标定。利用上述求得的角度值对罐容进行标定，计算结果参看图8、计算程序参看附录2Vprog2_24。图8 标定结果图7总结问题1用四次多项式灌油体离散变率进行回归得到连续变化率

14、，建立了连续型微分方程模型，通过该微分方程初值问题，得到倾斜灌油体体积随高度的变化规律，也可以采用问题2的方法去建立精确地积分模型得到倾斜灌油体体积随高度的变化规律。问题2的精确地积分模型中，变位油罐油体体积随高度的变化规律用积分表示的，由于不是油位、倾斜角以与旋转角的初等函数，在变为参数确定是带来困难，计算机编程时只能近似处理，二元被积函数在何处进行多项式展开是导致误差大小的关键问题；为了克服灌油体体积的积分表达式难以表示成油位、倾斜角以与旋转角的初等函数，采用了把倾斜角和旋转角的围化为一些离散小区间，在区间各分点得到测量数据与模型数据误差平方和，然后搜索到使得该平方和最小的倾斜角以与旋转角

15、，该方法一会有一定的误差，但倾斜角和旋转角的离散化越细密误差将会越小。参考文献1 茆诗松，程依明，濮晓龙 编。概率论与数理统计教程.：高等教育，2004.72 林 编。航天器轨道理论。：国防工业，2000.63 树嫄 编。微积分。：中国人民大学，2007.34 王沫然 编。MATLAB与科学计算（第二版）。：电子工业，2006.7附录1（问题1程序）% %无变位出油模型(油标底高：0.0286m,对应体积：0.0017l立方米,当高度1.2时体积为2.6523,灌体：4.1101m3)% load fujian1% h=xc(:,1);% v=xc(:,2);% dh=diff(h);% dv

16、=diff(v);% dvh=-dv./dh;% % plot(h(1:end-1),dvh,'k.')% b,bint,r,rint,sarts=regress(dvh,h(1:end-1).4,h(1:end-1).3,h(1:end-1).2,h(1:end-1),ones(size(h(1:end-1);% % hold on% f=poly2sym(b);% % ezplot(f,0,1.2);% vh1=dsolve('Dvh=-05335/1312*h4+80293/664*h3-2223/832*h2+86795/1312*h+89097/85248&#

17、39;,'vh(0)=0','h');% % figure(2)% ezplot(vh1,0,1.2);hold on;% q1=subs(vh1,1.2)% %变位进油模型(油标底高：0.0286m,对应体积：0.0017l立方米,当高度1.2时体积为2.6523,灌体：4.1101m3)% load fujian1% h=0;xxj(:,1);% v=0.00171;xxj(:,2);% dh=diff(h);% dv=diff(v);% dvh=dv./dh;% % plot(h(1:end-1),dvh,'k:.')% b,bint,r,

18、rint,sarts=regress(dvh,h(1:end-1).4,h(1:end-1).3,h(1:end-1).2,h(1:end-1),ones(size(h(1:end-1);% % hold on% f=poly2sym(b);% % ezplot(f,0,1.2);% vh2=dsolve('Dvh=-47049/0656*h4+97013/5328*h3-01607/5328*h2+9697/832*h+80339/70496','vh(0)=0.00171','h');% % figure(2)% q2=subs(vh2,1.

19、2)% ezplot(vh2,0,1.2)% ezplot(vh2-vh1)% axis(0,1.2,-1,5)syms h% 0.0017 0.5017 1.0017 1.5017 2.0017 2.5017 3.0017 3.5017v0=input('v0=') h01=double(solve(vh1-v0) h02=double(solve(vh2-v0)表1 倾斜灌融表（标定值）hvhvhvhv00.00170.30000.90230.60002.07570.91003.40820.01000.02040.31000.93890.61002.11750.92003.

20、44870.02000.04040.32000.97570.62002.15950.93003.48870.03000.06160.33001.01270.63002.20170.94003.52810.04000.08390.34001.04980.64002.24410.95003.56670.05000.10720.35001.08720.65002.28660.96003.60470.06000.13160.36001.12470.66002.32930.97003.64180.07000.15700.37001.16240.67002.37220.98003.67800.08000.

21、18320.38001.20020.69002.45840.99003.71330.09000.21020.39001.23820.70002.50171.00003.74760.10000.23810.40001.27640.71002.54511.01003.78080.11000.26670.41001.31480.72002.58861.02003.81270.12000.29590.42001.35330.73002.63221.03003.84340.13000.32590.43001.39200.74002.67591.04003.87280.14000.35640.44001.

22、43080.75002.71971.05003.90070.15000.38750.45001.46980.76002.76351.06003.92700.16000.41920.46001.50900.77002.80731.07003.95170.17000.45130.47001.54830.78002.85111.08003.97470.18000.48390.48001.58780.79002.89501.09003.99580.19000.51690.49001.62750.80002.93871.10004.01490.20000.55040.50001.66730.81002.

23、98241.11004.03200.21000.58420.51001.70730.82003.02601.12004.04690.22000.61840.52001.74750.83003.06951.13004.05950.23000.65290.53001.78790.84003.11281.14004.06960.24000.68770.54001.82850.85003.15591.15004.07720.25000.72290.55001.86920.86003.19881.16004.08210.26000.75830.56001.91010.87003.24141.17004.

24、08410.27000.79390.57001.95120.88003.28371.18004.08320.28000.82980.58001.99250.89003.32561.19004.07910.29000.86600.59002.03400.90003.36711.20004.0718附录2（问题2程序）%Vprog2_1(求、的MATLAB程序)a=1.5;%柱面半径1.5z0=8;R=1.625;%球体半径t=-7*pi/180:0.01:0;t1=0:0.01:6*pi/180;load alfahk=x(1:40:end,1);v=x(1:40:end,2);%容积测量值fo

25、r i=1:length(t) for j=1:length(t1) c=t(i);%sin(t)近似地用t代替; d=1-t(i)2/2;%cos(t)近似地用1-t2/2代替; XA=-3/2*d; XB=-3/2*d-8*c; XC=3/2*d-8*c; XD=-XA; for k=1:length(hk) h=XA-2*c+1/d*(1.5+(hk(k)-1.5)*(1-t1(j)2/2);syms X Z;f=sqrt(R2-(d*X+c*Z)2-(-c*X+d*Z-R+1)2);ft=tay(f,X,Z,0,-0.6);%f的泰勒展开式（2阶）,经试验在z=-0.6处结果好（），球

26、冠理论体积为4.0579，近似计算结果为4.0577Zm=d*(R-1)-sqrt(d2*R2-2*d2*R+d2-X2-2*c*X*R+2*c*X+2*R-1);SL=int(2*ft,Z,Zm,c/d*X);%VL=int(SL,X,XA,h);g=sqrt(R2-(d*X+c*Z)2-(-c*X+d*Z-z0+R-1)2);gt=tay(g,X,Z,0,z0+0.32);ZM=d*(1+z0-R)+sqrt(d2*(1+2*z0-2*R+z02-2*z0*R+R2)-X2+2*c*X*(R-1-z0)-2*(z0-R-z0*R)-z02-1);SR=int(2*gt,Z,z0/d+c/d

27、*X,ZM);%VR=int(SR,X,XB,h);f1=sqrt(a2-(d*X+c*Z)2+1e-3);%由于计算误差，避免根式下出现负值，加了一个小的正数。S1=int(2*f1,Z,c/d*X,-a/c-d/c*X);%V1=int(S1,X,XA,h)+VL;S2=int(2*f1,Z,c/d*X,z0/d+c/d*X);%V2=int(S1,X,XA,XB)+int(S2,X,XB,h)+VL+VR;%柱体体积为三部分分别为：2.98，53.5682，56.5487。S3=int(2*f1,Z,a/c-d/c*X,z0/d+c/d*X);% V3=int(S1,X,XA,XB)+i

28、nt(S2,X,XB,XD)+int(S3,X,XD,h)+int(SL,X,XA,XD)+VR; if h<XA V(k)=real(int(S1,X,XA,h)+int(SL,X,XA,h); elseif h>XB & h<=XD V(k)=real(int(S1,X,XA,XB)+int(S2,X,XB,h)+int(SL,X,XA,h)+int(SR,X,XB,h); else V(k)=real(int(S1,X,XA,XB)+int(S2,X,XB,XD)+int(S3,X,XD,h)+int(SL,X,XA,XD)+int(SR,X,XB,h); en

29、d end S(i,j)=sum(v'-double(V).2); i,j endendm=min(min(S);U1,U2=find(S=m);alfa=t(U1)beita=t1(U2)%Vprog2_2 (罐容标定计算)a=1.5;%柱面半径1.5z0=8;R=1.625;%球体半径t=-0.001;t1=0.171;load alfahk=0:0.1:3; c=t;%sin(t)近似地用t代替; d=1-t2/2;%cos(t)近似地用1-t2/2代替; XA=-3/2*d; XB=-3/2*d-8*c; XC=3/2*d-8*c; XD=-XA; for k=1:length(hk) h=XA-2*c+1/d*(1.5+(hk(k)-1.5)*(1-t12/2);syms X Z;f=sqrt(R2-(d*X+c*Z)2-(-c*X+