储油罐问题.doc

问题重述通常加油站都采用与储油罐配套的“油位计量管理系统”。采用流量计和油位计来测量进/出油量与罐内油位高度等数据，通过预先标定的罐容表（即罐内油位高度与储油量的对应关系）进行实时计算，以得到罐内油位高度和储油量的变化情况。使用一段时间后，由于地基变形等原因，使罐体的位置会发生纵向倾斜和横向偏转等变化（以下称为变位），油箱液面不再平行于油箱底部，从而使得容积的计算公式不再准确，导致罐容表发生改变，需要重新定位。图1是一种典型的储油罐尺寸及形状示意图，其主体为圆柱体，两端为球冠体。图2是其罐体纵向倾斜变位的示意图，图3是罐体横向偏转变位的截面示意图。要求用数学建模方法研究解决储油罐的变位识别与罐容表标定的问题。 （1）为了掌握罐体变位后对罐容表的影响，利用如图4的小椭圆型储油罐（两端平头的椭圆柱体），分别对罐体无变位和倾斜角为a=4.10的纵向变位两种情况做了实验，实验数据如附件1所示。请建立数学模型研究罐体变位后对罐容表的影响，并给出罐体变位后油位高度间隔为1cm的罐容表标定值。（2）对于图1所示的实际储油罐，试建立罐体变位后标定罐容表的数学模型，即罐内储油量与油位高度及变位参数（纵向倾斜角度a和横向偏转角度b ）之间的一般关系。请利用罐体变位后在进/出油过程中的实际检测数据（附件2），根据你们所建立的数学模型确定变位参数，并给出罐体变位后油位高度间隔为10cm的罐容表标定值。进一步利用附件2中的实际检测数据来分析检验你们模型的正确性与方法的可靠性。油油浮子出油管油位探测装置注油口检查口地平线2m6m1m1m3 m油位高度图1 储油罐正面示意图油位探针油位探针 地平线图2 储油罐纵向倾斜变位后示意图油油浮子出油管油位探测装置注油口检查口水平线(b) 小椭圆油罐截面示意图 油油浮子出油管油位探针注油口水平线2.05mcm0.4m1.2m1.2m1.78m(a) 小椭圆油罐正面示意图图4 小椭圆型油罐形状及尺寸示意图图3 储油罐截面示意图（b）横向偏转倾斜后正截面图地平线垂直线油位探针（a）无偏转倾斜的正截面图油位探针油位探测装置3m14问题一无变位情形1. 建立坐标系将椭圆柱体的储油箱水平放置，箱底接触地面的只有一条直线，以该条直线为x轴，油位探针所在直线为y轴，两轴所在直线相交与O点，过O点做z轴垂直于x,y轴所确定的平面。x2.05mzO-0.4m,Ohyy2公式推导在z-y平面内列写椭圆方程： ；椭圆中阴影部分面积为；油液容积为。根据如上的公式，整理得到容积V与高度h之间的函数关系式为： 其中 a= 0.89, b=0.6, l=2.45 (m)分别为罐体截面椭圆长半轴、短半轴和罐体长度， h 为罐内的油位高度油量/L3659.883946.554110.15油量/L3659.883946.554110.15有变位情形 坐标系同上，无变位时，油平面法向为（0,1,0），由于发生纵向变位，相当于坐标系绕z轴旋转了一个角度，因此，油平面法向变为（，油平面方程为无变位时，油平面高度不随x坐标变化，当有变位时，油平面高度是x的函数，表达式为计算体积的Matlab程序function v= getv(h)%要求matlab版本较高a=0.89;b=0.6;alpha=4.1*3.1416/18;fun=(x,y)a*sqrt(1-(y-b)/b).2);hx=(x)min(max(0,h-x*tan(alpha),1.2);v=2*quad2d(fun,-0.4,2.05,0,hx);end或者function v= getv2(h)%低版本atlab也可以运行a=0.89;b=0.6;alpha=4.1*3.1416/180; fun=(y)a*sqrt(1-(y-b)/b).2);hx=(x)min(max(0,h-x*tan(alpha),1.2);f1=(u)quadl(fun,0,hx(u);f2=(x)arrayfun(f1,x);v=2*quadl(f2,-0.4,2.05);end问题二对于实际的储油罐，坐标系的建立类似于问题一，将y轴建立在油位探针上，原点取油位探针的最低点，x轴为圆柱体水平放置时和地面的接触线，因此，储油罐左边球缺的x坐标范围为（-3，-2），右边球缺的x 坐标范围是（6，7），中间圆柱体部分的x坐标范围是（-2,6），圆柱体底圆半径为1.5米，容易计算出球缺的半径为1.625米当储油罐发生纵向旋转角度，相当于坐标系关于z轴旋转了，旋转矩阵为油平面法向由（0,1,0）变为，而当储油罐横向旋转时，相当于坐标系关于x轴旋转了，旋转矩阵为此时方向由变为。不管罐体如何旋转，油平面总是经过点（0,h,0）,因此油平面方程为圆柱体油量的计算用垂直于x轴的平面截储油罐，截得是圆，圆半径为 根据油平面方程，油平面和xy面交线为解出为，考虑到球缺上的点和 xz面的距离，记. 当储油罐横向偏转时，油体和垂直于x轴的截面为弓形，其高度为由于我们只在储油罐所在区域考虑问题，应作修正，修正后的表达式为由弓形的面积公式，我们得到油体和垂直于x轴的截面的截面积为于是，油量为变位参数的确定数据给定了各时间段的进出油量值D及显示油位高度D，据此可以计算参数。将油位高度值代入体积表达式中，计算得到=，于是得到,通过最小二乘准则，利用数据编程计算得到变位参数与的估计值（结果为）。计算截面积的matlab程序function s=surf(x,h,alpha,beta)alpha=alpha/180*pi;beta=beta/180*pi; if(x=-3) rx=sqrt(13/8)2-(x+11/8)2);else if (x6)&(x=7) rx=sqrt(13/8)2-(x-43/8)2);else rx=3/2; endendhx=h-x*tan(alpha)/cos(beta)-(3/2-rx); H=min(2*rx,max(rx-(rx-hx)*cos(beta),0); if rx=0 s=0;elses=rx2*acos(rx-H)/rx)-sqrt(rx2-(rx-H)2)*(rx-H) ;endend计算体积的matlab程序function v=vol(h,alpha,beta)f=(x)arrayfun(x)surf(x,h,alpha,beta),x); v=1000*quadl(f,-3,7);end计算误差平方和的函数function y=errsum(alpha,beta) % dv 和h来自于题中少量数据，实际计算时应多取点数据dv=149.09 68.45 199.27 70.05 136.36 232.74 107.97 49.24 80.65 120.29;h=2632.23 2624.30 2620.67 2610.29 2606.61 2599.59 2587.60 2582.05 2579.57 2575.44 2569.46 ; for i=1:length(h) hi=h(i)/1000;v(i)=vol(hi,alpha,beta);enddv1=-diff(v);y=sum(dv-dv1).2);end直接搜索参数的程序，步长较长，计算时步长应逐渐变小。errmin=inf;for alpha=0.5:0.5:5 for beta=0.5:0.5:5 errs=errsum(alpha,beta); if errserrmin errmin=errs alphamin=alpha betamin=beta end endend计算结果为。在这个结果基础上，可进一步搜索。加快计算速度的方法上述程序在实际计算的时候，速度慢，如果实际数据取得较多的时候，计算速度将会更慢，为了加快速度，我们可以采用matlab和c语言混合编程的方法来解决。计算参数时，影响速度的主要因素为函数的计算和积分的计算，我们可以用C语言实现，将其编译为mex文件，供matlab 调用。下面我们先用一个简单例子来说明，1. mex 配置在命令行输入 mex setup 然后根据提示操作即可。2. 编写c语言文件/* mymex.c*/#include mex.hvoid mexFunction (int nlhs,mxArray* plhs,int nrhs,const mxArray* prhs)mexPrintf(c语言编译为mex文件例子);说明mexfunction 是c语言编写mex文件必须要有的函数，相当于main()函数。参数说明：nlhs:输出参数个数 plhs: 输出参数的mxArray 数组 nrhs:输入参数个数 prhs:输入参数的mxArray数组 3将c程序编译为mex文件mex mymex.c编译成功后，将生成一个mymex.mexw32或mymex.mexw64文件 4.在命令行输入mymex,将会 输出c语言编译为mex文件例子对于本题，我们用C语言写求数值积分的程序，可以将速度提高100多倍，不过误差有0.2升左右。程序如下#include mex.h double surf(double x,double h,double alpha,double beta)double rx,HH,hx,s; alpha=alpha/180*3.1416;beta=beta/180*3.1416; if (x-3.0) )rx=sqrt(2.6406-(x+11.0/8.0)*(x+11.0/8.0); else if (x=6.0)&(x2.0*rx ) HH=2.0*rx; if (rx=0.0|HH=0.0) s=0.0;elses=rx*rx*acos(rx-HH)/rx)-sqrt(rx*rx-(rx-HH)*(rx-HH)*(rx-HH) ;return s; double Simpson(double(*f)(double,double,double,double),double a,double b,int n,double hh,double alpha,double beta) int k; double s,s1,s2; double h=(b-a)/n; s1=f(a+h/2,hh,alpha,beta); s2=0.0; for(k=1;k=n-1;k+) s1+=f(a+k*h+h/2,hh,alpha,beta); s2+=f(a+k*h,hh,alpha,beta); s=h/6*(f(a,hh,alpha,beta)+4*s1+2*s2+f(b,hh,alpha,beta); return s; voidmexFunction(int nlhs,mx