储油罐问题Word格式.docx
- 文档编号:21808161
- 上传时间:2023-02-01
- 格式:DOCX
- 页数:14
- 大小:215.73KB
储油罐问题Word格式.docx
《储油罐问题Word格式.docx》由会员分享,可在线阅读,更多相关《储油罐问题Word格式.docx(14页珍藏版)》请在冰豆网上搜索。
椭圆中阴影部分面积为
油液容积为
。
根据如上的公式,整理得到容积V与高度h之间的函数关系式为:
其中a=0.89,b=0.6,l=2.45(m)分别为罐体截面椭圆长半轴、短半轴和罐体长度,h为罐内的油位高度
油量/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)
alpha=4.1*3.1416/180;
fun=@(y)a*sqrt(1-((y-b)/b).^2);
f1=@(u)quadl(fun,0,hx(u));
f2=@(x)arrayfun(f1,x);
v=2*quadl(f2,-0.4,2.05);
问题二
对于实际的储油罐,坐标系的建立类似于问题一,将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
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;
s=rx^2*acos((rx-H)/rx)-sqrt(rx^2-(rx-H)^2)*(rx-H);
计算体积的matlab程序
functionv=vol(h,alpha,beta)
f=@(x)arrayfun(@(x)surf(x,h,alpha,beta),x);
v=1000*quadl(f,-3,7);
计算误差平方和的函数
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);
dv1=-diff(v);
y=sum((dv-dv1).^2);
直接搜索参数的程序,步长较长,计算时步长应逐渐变小。
errmin=inf;
foralpha=0.5:
0.5:
5
forbeta=0.5:
errs=errsum(alpha,beta);
iferrs<
errmin
errmin=errs
alphamin=alpha
betamin=beta
计算结果为
在这个结果基础上,可进一步搜索。
加快计算速度的方法
上述程序在实际计算的时候,速度慢,如果实际数据取得较多的时候,计算速度将会更慢,为了加快速度,我们可以采用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升左右。
程序如下
doublesurf(doublex,doubleh,doublealpha,doublebeta)
{doublerx,HH,hx,s;
alpha=alpha/180*3.1416;
beta=beta/180*3.1416;
if((x<
=-2.0)&
&
-3.0))
rx=sqrt(2.6406-(x+11.0/8.0)*(x+11.0/8.0));
elseif((x>
=6.0)&
7.0))
rx=sqrt(2.6406-(x-43.0/8.0)*(x-43.0/8.0));
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;
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文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 储油罐 问题
![提示](https://static.bdocx.com/images/bang_tan.gif)