柴油机燃烧MATLAB计算方法

1、%-从文件读取数据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

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;endfor 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;endY(N-2)=(A(N-6)-4*(A(N-5)+A(N)+A(N-4)+12

3、*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

4、.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; %压缩始点角（选取示功图压力变

5、化开始点）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; %计算角度

6、平移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(

7、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*d2*(s/(skm-1)+x(j); %瞬时体积%-传热计算参数 Aw=pi*d*(Sx+x(j); %缸壁面

8、积 At=0.25*pi*d2; %缸盖面积 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); %缸内非燃烧期间温度计算（燃烧期间温度将在下面重新赋值） endend %-内能变化计算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(

9、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(

10、j)=103*(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*Cm0.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)

11、=DQ(j); %总放热量 % else DW(j)=106*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*Cm0.7; %sitkei公式 DQW(j)=103*(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); %总放热量 endend%-图形输出figuresubplot 121plot(Camst+1:Camend+1),T(Camst+1:Camend+1);xlabel('crank-deg');ylabel('Temp.K');gr