柴油機燃燒MATLAB計算方法

1、%-從文件讀取數(shù)據(jù)filename,pathname=uigetfile('*.txt''*.dat','選擇示功圖數(shù)據(jù)文件');str_filename=pathname,filename;Deg,A = textread(str_filename,'%s %s','headerlines',2);%讀取數(shù)據(jù)Deg= str2num(char(Deg); %讀取角度數(shù)據(jù) 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;%-預(yù)設(shè)初始計算參數(shù)%-若需修改模型參數(shù)，可在以下修改gc=0.870; gh=0.126; go=0.004; %柴油平均質(zhì)量成分（碳、氫、氧） tao=4; fai=1.06; %沖程、掃氣系數(shù)、 n=2500; i=6; %轉(zhuǎn)速，汽缸數(shù) bc=235

4、.6; %燃油消耗率 Pe=161.51; %有效功率 Gma=0.02; %殘余廢氣系數(shù) 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; %進氣參數(shù)（壓力、空氣相對質(zhì)量常數(shù)、環(huán)境溫度、環(huán)境壓力、比熱容比） Sx=1.9*10(-3); Twt=600; Twb=600; Tww=600; %余隙高度、壁面溫度（缸蓋、活塞、缸套）b=0.25; %sitkei公式經(jīng)驗常數(shù)（取值參照該公式）ps=277; %壓縮始點角（選取示功圖壓力變

5、化開始點）Camst=300; Camend=450; %計算始點角，計算終點角,360度為上止點 Alf=1.9; %過量空氣系數(shù) %-非改動參數(shù) Mua=28.97; %空氣相對分子質(zhì)量R0=8314; %氣體常數(shù) Hu=42287; %低熱值(柴油) ao=4.678; bo=6.8723*10(-4); co=-6.0683*10(-8); %F.schmidt氣體回歸比熱容參數(shù)（空氣）ar=4.7513; br=1.199*10(-3); cr=-1.4232*10(-7); %F.schmidt氣體回歸比熱容參數(shù)（純?nèi)紵a(chǎn)物）%-計算過程穩(wěn)定參數(shù)deg=Deg+360; %計算角度

6、平移gf=Pe*bc*tao*10(-3)/(120*n*i); %循環(huán)噴油量 namd=2*cor/s; %連桿曲柄比L0=(1/0.21)*(gc/12+gh/4-go/32); %燃料燃燒理論空氣量 b0=1+(gh/4+go/32)/L0; %理論分分子變化系數(shù)Ta=tn*(p0/pn)(k-1)/k); %進氣初始溫度%-計算參數(shù)初始化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(d

7、eg);Afa=zeros(size(deg);afa=10000+Afa;T=zeros(size(deg);T=T+Ta;%-活塞行程參數(shù)計算for j=1:N caudeg(j)=asin(sqrt(sin(pi/180*deg(j)/namd)2); %計算角度轉(zhuǎn)換 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); %瞬時體積%-傳熱計算參數(shù) Aw=pi*d*(Sx+x(j); %缸壁面積 At=0.25*pi*d2; %缸

8、蓋面積 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); %缸內(nèi)非燃燒期間溫度計算（燃燒期間溫度將在下面重新賦值） endend %-內(nèi)能變化計算for j=(Camst+1):(Camend+1) X(j)=Q(j-1)/(gf*Hu*1e3); %燃料燃燒百分數(shù) Ma=Alf*gf*L0*(1+Gma); %壓縮始點物質(zhì)的量 M(j)=Ma*(1+0.065*X(j)/(1+Gma)*Alf); %某

9、時刻物質(zhì)的量 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); %純?nèi)紵a(chǎn)物所占比例 %-cva=4.1868*(ao+bo*T(j)+co*T(j)2); %F.Schimidt（純?nèi)紵a(chǎn)物定容比熱）cve=4.1868*(ar+br*T(j)+cr*T(j)2); %F.Schimidt（空氣定容比熱）cv(j)=kr(j)*cve+(1-kr(j)*cva; %定容比熱 %- U(j)=103*(M(j)*cv(j)*

10、T(j)-Ma*cv(j)*Ta); %內(nèi)能變化 DU(j)=U(j)-U(j-1); %瞬態(tài)內(nèi)能變化%-做功計算 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)=DQ(j); %總放熱量 % els

11、e 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');gri