储油罐问题.docx
- 文档编号:8671965
- 上传时间:2023-02-01
- 格式:DOCX
- 页数:14
- 大小:215.73KB
储油罐问题.docx
《储油罐问题.docx》由会员分享,可在线阅读,更多相关《储油罐问题.docx(14页珍藏版)》请在冰豆网上搜索。
储油罐问题
问题重述
通常加油站都采用与储油罐配套的“油位计量管理系统”。
采用流量计和油位计来测量进/出油量与罐内油位高度等数据,通过预先标定的罐容表(即罐内油位高度与储油量的对应关系)进行实时计算,以得到罐内油位高度和储油量的变化情况。
使用一段时间后,由于地基变形等原因,使罐体的位置会发生纵向倾斜和横向偏转等变化(以下称为变位),油箱液面不再平行于油箱底部,从而使得容积的计算公式不再准确,导致罐容表发生改变,需要重新定位。
图1是一种典型的储油罐尺寸及形状示意图,其主体为圆柱体,两端为球冠体。
图2是其罐体纵向倾斜变位的示意图,图3是罐体横向偏转变位的截面示意图。
要求用数学建模方法研究解决储油罐的变位识别与罐容表标定的问题。
(1)为了掌握罐体变位后对罐容表的影响,利用如图4的小椭圆型储油罐(两端平头的椭圆柱体),分别对罐体无变位和倾斜角为=4.10的纵向变位两种情况做了实验,实验数据如附件1所示。
请建立数学模型研究罐体变位后对罐容表的影响,并给出罐体变位后油位高度间隔为1cm的罐容表标定值。
(2)对于图1所示的实际储油罐,试建立罐体变位后标定罐容表的数学模型,即罐内储油量与油位高度及变位参数(纵向倾斜角度和横向偏转角度)之间的一般关系。
请利用罐体变位后在进/出油过程中的实际检测数据(附件2),根据你们所建立的数学模型确定变位参数,并给出罐体变位后油位高度间隔为10cm的罐容表标定值。
进一步利用附件2中的实际检测数据来分析检验你们模型的正确性与方法的可靠性。
问题一
无变位情形
1.建立坐标系
将椭圆柱体的储油箱水平放置,箱底接触地面的只有一条直线,以该条直线为x轴,油位探针所在直线为y轴,两轴所在直线相交与O点,过O点做z轴垂直于x,y轴所确定的平面。
2公式推导
在z-y平面内列写椭圆方程:
;
椭圆中阴影部分面积为
;
油液容积为
。
根据如上的公式,整理得到容积V与高度h之间的函数关系式为:
其中a=0.89,b=0.6,l=2.45(m)分别为罐体截面椭圆长半轴、短半轴和罐体长度,h为罐内的油位高度
油量/L
3659.88
3946.55
4110.15
油量/L
3659.88
3946.55
4110.15
有变位情形
坐标系同上,无变位时,油平面法向为(0,1,0),由于发生纵向变位,相当于坐标系绕z轴旋转了一个角度
,因此,油平面法向变为(
,油平面方程为
无变位时,油平面高度不随x坐标变化,当有变位时,油平面高度是x的函数,表达式为
计算体积的Matlab程序
functionv=getv(h)
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
或者
functionv=getv2(h)
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面交线为
解出y为
,考虑到球缺上的点和xz面的距离,记
.
当储油罐横向偏转时,油体和垂直于x轴的截面为弓形,其高度为
由于我们只在储油罐所在区域考虑问题,
应作修正,修正后的表达式为
由弓形的面积公式,我们得到油体和垂直于x轴的截面的截面积为
于是,油量为
变位参数的确定
数据给定了各时间段的进出油量值
及显示油位高度
,据此可以计算参数。
将油位高度值
代入体积表达式
中,计算得到
,于是得到
通过最小二乘准则
,利用数据编程计算得到变位参数
与
的估计值(结果为
)。
计算截面积的matlab程序
functions=surf(x,h,alpha,beta)
alpha=alpha/180*pi;beta=beta/180*pi;
if(x<=-2)&(x>=-3)
rx=sqrt((13/8)^2-(x+11/8)^2);
elseif(x>6)&(x<=7)
rx=sqrt((13/8)^2-(x-43/8)^2);
else
rx=3/2;
end
end
hx=h-x*tan(alpha)/cos(beta)-(3/2-rx);
H=min(2*rx,max(rx-(rx-hx)*cos(beta),0));
ifrx==0
s=0;
else
s=rx^2*acos((rx-H)/rx)-sqrt(rx^2-(rx-H)^2)*(rx-H);
end
end
计算体积的matlab程序
functionv=vol(h,alpha,beta)
f=@(x)arrayfun(@(x)surf(x,h,alpha,beta),x);
v=1000*quadl(f,-3,7);
end
计算误差平方和的函数
functiony=errsum(alpha,beta)
%dv和h来自于题中少量数据,实际计算时应多取点数据
dv=[149.0968.45199.2770.05136.36232.74107.9749.2480.65120.29];
h=[2632.232624.302620.672610.292606.612599.592587.602582.052579.572575.442569.46];
fori=1:
length(h)
hi=h(i)/1000;
v(i)=vol(hi,alpha,beta);
end
dv1=-diff(v);
y=sum((dv-dv1).^2);
end
直接搜索参数的程序,步长较长,计算时步长应逐渐变小。
errmin=inf;
foralpha=0.5:
0.5:
5
forbeta=0.5:
0.5:
5
errs=errsum(alpha,beta);
iferrs errmin=errs alphamin=alpha betamin=beta end end end 计算结果为 。 在这个结果基础上,可进一步搜索。 加快计算速度的方法 上述程序在实际计算的时候,速度慢,如果实际数据取得较多的时候,计算速度将会更慢,为了加快速度,我们可以采用matlab和c语言混合编程的方法来解决。 计算参数时,影响速度的主要因素为函数的计算和积分的计算,我们可以用C语言实现,将其编译为mex文件,供matlab调用。 下面我们先用一个简单例子来说明, 1.mex配置 在命令行输入mex–setup然后根据提示操作即可。 2.编写c语言文件 /*mymex.c*/ #include"mex.h" voidmexFunction(intnlhs,mxArray*plhs[],intnrhs,constmxArray*prhs[]) {mexPrintf("c语言编译为mex文件例子"); } 说明 mexfunction是c语言编写mex文件必须要有的函数,相当于main()函数。 参数说明: nlhs: 输出参数个数 plhs: 输出参数的mxArray数组 nrhs: 输入参数个数 prhs: 输入参数的mxArray数组 3.将c程序编译为mex文件 mexmymex.c 编译成功后,将生成一个mymex.mexw32或mymex.mexw64文件 4.在命令行输入mymex,将会输出c语言编译为mex文件例子 对于本题,我们用C语言写求数值积分的程序,可以将速度提高100多倍,,不过误差有0.2升左右。 程序如下 #include"mex.h" doublesurf(doublex,doubleh,doublealpha,doublebeta) {doublerx,HH,hx,s; alpha=alpha/180*3.1416;beta=beta/180*3.1416; if((x<=-2.0)&&(x>-3.0)) rx=sqrt(2.6406-(x+11.0/8.0)*(x+11.0/8.0)); elseif((x>=6.0)&&(x<7.0)) rx=sqrt(2.6406-(x-43.0/8.0)*(x-43.0/8.0)); else rx=1.5; if(x==-3.0||x==7)rx=0; hx=h-x*tan(alpha)/cos(beta)-(1.5-rx); HH=rx-(rx-hx)*cos(beta); if(HH>2.0*rx) HH=2.0*rx; if(rx==0.0||HH<=0.0) s=0.0; else s=rx*rx*acos((rx-HH)/rx)-sqrt(rx*rx-(rx-HH)*(rx-HH))*(rx-HH); returns; } doubleSimpson(double(*f)(double,double,double,double),doublea,doubleb,intn,doublehh,doublealpha,doublebeta) { intk; doubles,s1,s2; doubleh=(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)); returns; } void mexFunction(intnlhs,mxArray*plhs[],intnrhs,constmxArray*prhs[]) { inti; doublex,h,alpha,beta,v; if(nrhs<3)return; h=*mxGetPr(prhs[0]); alpha=*mxGetPr(prhs[1]); beta=*mxGetPr(prhs[2]); v=1000*Simpson(surf,-3,7,200,h,alpha,beta); plhs[0]=mxCreateDoubleScalar(v); } Forpersonaluseonlyinstudyandresearch;notforcommercialuse
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 储油罐 问题