柴油机燃烧MATLAB计算方法.doc

%------------------------------------------------------------------------------------------从文件读取数据
[filename,pathname]=uigetfile({'*.txt';'*.dat'},'选择示功图数据文件');
str_filename=[pathname,filename];
[Deg,A] = textread(str_filename,'%s %s','headerlines',2);%读取数据
Deg= str2num(char(Deg)); %读取角度数据
A = str2num(char(A)); %读取未光顺压力值 Mpa
%------------------------------------------------------------------------------------------七点五次光顺
N=numel(A);
Y(1)=(39*A(1)+8*A(2)-4*(A(3)+A(4)-A(6))+A(5)-2*A(7))/42;
Y(2)=(8*A(1)+19*A(2)+16*A(3)+6*A(4)-4*A(5)-7*A(6)+4*A(7))/42;
Y(3)=(-4*(A(1)+A(6))+16*A(2)+19*A(3)+12*A(4)+2*A(5)+7*A(7))/42;
if N<7
return;
end
for k=4:(N-3)
Y(k)=(-2*(A(k-3)+A(k+3))+3*(A(k-2)+A(k+2))+6*(A(k-1)+A(k+1))+7*A(k))/21;
end
Y(N-2)=(A(N-6)-4*(A(N-5)+A(N))+A(N-4)+12*A(N-3)+19*A(N-2)+16*A(N-1))/42;
Y(N-1)=(4*A(N-6)-7*A(N-5)+4*A(N-4)+6*A(N-3)+16*A(N-2)+19*A(N-1)+8*A(N))/42;
Y(N)=(-2*A(N-6)+4*(A(N-5)-A(N-3)-A(N-2))+A(N-4)+8*A(N-1)+39*A(N))/42;
%-----------------------------------------------------预设初始计算参数
%-----------------------------------------------------若需修改模型参数，可在以下修改
gc=0.870; gh=0.126; go=0.004; %柴油平均质量成分（碳、氢、氧）
tao=4; fai=1.06; %冲程、扫气系数、
n=2500; i=6; %转速，汽缸数
bc=235.6; %燃油消耗率
Pe=161.51; %有效功率
Gma=0.02; %残余废气系数
d=0.102; s=0.12; skm=17.3; cor=0.192; %缸径、冲程、压缩比、连杆长度
p0=0.261; Ra=287.06; tn=298.15; pn=0.10133; k=1.4; %进气参数（压力、空气相对质量常数、环境温度、环境压力、比热容比）
Sx=1.9*10^(-3); Twt=600; Twb=600; Tww=600; %余隙高度、壁面温度（缸盖、活塞、缸套）
b=0.25; %sitkei公式经验常数（取值参照该公式）
ps=277; %压缩始点角（选取示功图压力变化开始点）
Camst=300; Camend=450; %计算始点角，计算终点角,360度为上止点
Alf=1.9; %过量空气系数
%-------------------------------------------------------------------------------------非改动参数
Mua=28.97; %空气相对分子质量
R0=8314; %气体常数
Hu=42287; %低热值(柴油)
ao=4.678; bo=6.8723*10^(-4); co=-6.0683*10^(-8); %F.schmidt气体回归比热容参数（空气）
ar=4.7513; br=1.199*10^(-3); cr=-1.4232*10^(-7); %F.schmidt气体回归比热容参数（纯燃烧产物）
%-------------------------------------------------------------------------------------计算过程稳定参数
deg=Deg+360; %计算角度平移
p=0.1.*Y; %压力值单位转换
gf=Pe*bc*tao*10^(-3)/(120*n*i); %循环喷油量
namd=2*cor/s; %连杆曲柄比
L0=(1/0.21)*(gc/12+gh/4-go/32); %燃料燃烧理论空气量
b0=1+(gh/4+go/32)/L0; %理论分分子变化系数
Ta=tn*(p0/pn)^((k-1)/k); %进气初始温度
%--------------------------------------------------------------------------------------计算参数初始化
Q=zeros(size(deg));
X=zeros(size(deg));
M=zeros(size(deg));
DW=zeros(size(deg));
DQW=zeros(size(deg));
DU=zeros(size(deg));
U=zeros(size(deg));
Afa=zeros(size(deg));
afa=10000+Afa;
T=zeros(size(deg));
T=T+Ta;
%------------------------------------------------－--－－－－－－－－－－－－------------活塞行程参数计算
for j=1:N
caudeg(j)=asin(sqrt((sin(pi/180*deg(j))/namd)^2)); %计算角度转换
x(j)=(s/2)*((1+namd)-namd*sqrt((cos(caudeg(j)))^2)-cos(pi/180*deg(j))); %活塞离缸盖距离
v(j)=0.25*pi*d^2*(s/(skm-1)+x(j)); %瞬时体积
%--------------------------------------------------------------传热计算参数
Aw=pi*d*(Sx+x(j)); %缸壁面积
At=0.25*pi*d^2; %缸盖面积
Ab=1.25*At; %活塞顶面积
de=2*d*(Sx+x(j))/(d+2*(Sx+x(j))); %当量直径
Cm=n*2*0.12/60; %活塞平均速度
if p(j)>ps
T(j)=tn*(p(j)/pn)^((k-1)/k); %缸内非燃烧期间温度计算（燃烧期间温度将在下面重新赋值）
end
end
%------------------------------------------------------------－－－－－－－－－－－－--内能变化计算
for j=(Camst+1):(Camend+1)
X(j)=Q(j-1)/(gf*Hu*1e3); %燃料燃烧百分数
Ma=Alf*gf*L0*(1+Gma); %压缩始点物质的量
M(j)=Ma*(1+0.065*X(j)/((1+Gma)*Alf)); %某时刻物质的量
T(j) = 1e6*p(j)*v(j)/(M(j)*R0); %瞬时温度
kr(j)=((Alf-1+b0)*b0*X(j)+b0*Alf*Gma)/((Alf-1+b0)*((1+Gma)*Alf+X(j)*(b0-1))); %纯燃烧产物所占比例
%-------------------------------------------------------
cva=4.1868*(ao+bo*T(j)+co*T(j)^2); %F.Schimidt（纯燃烧产物定容比热）
cve=4.1868*(ar+br*T(j)+cr*T(j)^2); %F.Schimidt（空气定容比热）
cv(j)=kr(j)*cve+(1-kr(j))*cva; %定容比热
%--------------------------------------------------------
U(j)=10^3*(M(j)*cv(j)*T(j)-Ma*cv(j)*Ta); %内能变化
DU(j)=U(j)-U(j-1); %瞬态内能变化
%--------------------------------------------------------－－－－－－－－－－－－－－---做功计算
if j==1
DW(j)=0.5*(p(j)+p(N-1))*(v(j)-v(N-1)); %瞬时做功
ag(j)=0.205*(1+b)*de^(-0.3)*T(j)^(-0.2)*p(j)^0.7*Cm^0.7; %sitkei公式
DQW(j)=(1/(6*n))*ag(j)*(At*(T(j)-Twt)+Ab*(T(j)-Twb)+Aw*(T(j)-Tww)); %瞬时传热量
DQ(j)=DU(j)+DW(j)+DQW(j); %瞬时放热量
Q(j)=DQ(j); %总放热量
%－－－－－－－－－－－－－－－－－－－－－－－－－
else
DW(j)=10^6*0.5*(p(j)+p(j-1))*(v(j)-v(j-1)); %瞬时做功
ag(j)=0.205*(1+b)*de^(-0.3)*T(j)^(-0.2)*p(j)^0.7*Cm^0.7; %sitkei公式
DQW(j)=10^3*(1/(6*n))*ag(j)*(At*(T(j)-Twt)+Ab*(T(j)-Twb)+Aw*(T(j)-Tww)); %瞬时传热量
DQ(j)=DU(j)+DW(j)+DQW(j); %单位放热量
Q(j)=Q(j-1)+DQ(j); %总放热量
end
end
%-------------------------------------------------------------------------－－－－－－图形输出
figure
subplot 121
plot((Camst+1:Camend+1),T(Camst+1:Camend+1));
xlabel('crank-deg');ylabel('Temp.K');
grid
subplot 122
plot((Camst+1:Camend+1),DQ(Camst+1:Camend+1));
xlabel('crank-deg');ylabel('ROHR.1/deg');
grid
hold on;
figure
subplot 121
plot(A);
grid
subplot 122
plot(p);
grid