unction [ ETA, IMEP, NOX_ppm ] = Homogeneous(varargin)
% Homogeneous Two Zone Combustion Cycle
% This program computes the pressure and temperature vs crank angle
% and the Work, Indicated thermal efficiency - ETA
% and the Indicated mean effective pressure - IMEP (kPa)
%
%CFR Engine geometry and operating conditions
R = 12; % Compression ratio - R
B = .0825; % Bore - B (m)
S = .1143; % Stroke - S (m)
EPS = 0.25; % Half stroke to rod ratio - EPS
RPM = 1200; % Engine speed - RPM
HEAT = 500; % Heat transfer coefficient
BLOWBY = 0.8; % Blowby coefficient
THETAS = -25; % Start of heat release (deg ATDC)
THETAB = 60; % Burn angle (deg)
PHI = 0.88; % Equivalence ratio - PHI
F = 0.05; % Residual fraction - F
TW = 400; % Wall temperature - TW
fuel_type = 2; % gasoline
FS = 0.06548; % stoichiometric fuel-air ratio for gasoline
A0 = 47870; %
T1 = 350;
P1 = 100; % kPa
if ( nargin == 3 )
PHI = varargin{1};
F = varargin{2};
RPM = varargin{3};
end
OMEGA = RPM*pi/30;
to_ppm = 10^6; % convert from mass fraction to ppm
MW_NO = 30; % molecular weight of NO, g/mol
THETA = -180;
DTHETA = 1;
THETAE = THETA+DTHETA;
[ VOL, X, EM ] = auxiliary( THETA );
NNOX = THETAB/DTHETA;
NY = 6+NNOX;
Y = zeros(NY,1);
Y(1) = P1;
Y(2) = nan;
Y(3) = T1;
[~, ~, ~, ~, vU, ~, ~, ~, ~, ~] = farg( Y(3), Y(1), PHI, F, fuel_type );
MNOT = VOL/vU;
M = EM*MNOT;
NN = 36*10;
SAVE.THETA = zeros( NN, 1 );
SAVE.VOL = zeros( NN, 1 );
SAVE.T = zeros(NN, 1 );
SAVE.P = zeros( NN, 1 );
SAVE.MDOTFI = zeros( NN, 1 );
SAVE.NOX = zeros(NN,5);
SAVE.NOe = zeros(NN,1);
fprintf( 'THETA VOLUME BURN FRAC PRESSURE BURN TEMP UNBURNED T WORK HEAT LOSS MASS H-LEAK NOx\n' );
fprintf( ' deg cm^3 -- kPa K K J J g J ppm\n' );
fprintf('%7.1f %6.1f %3.3f %6.1f %6.1f %6.1f %5.0f %5.0f %5.3f %5.2f %6.1f\n', ...
THETA, VOL*1000000, X, Y(1), Y(2), Y(3), Y(4)*1000, Y(5)*1000, M*1000, Y(6)*1000, 0.0 );
II = 1;
for III=1:36,
for JJJ=1:10,
NOe_save = 0;
[ Y ] = integrate( THETA, THETAE, Y );
[ VOL, X, EM ] = auxiliary( THETA );
M = EM*MNOT;
THETA=THETAE;
THETAE=THETA+DTHETA;
% save data for plotting later
SAVE.THETA(II) = THETA;
SAVE.VOL(II) = VOL;
SAVE.P(II) = Y(1);
SAVE.TB(II) = Y(2);
SAVE.TU(II) = Y(3);
SAVE.X(II) = X;
SAVE.NOX(II,:) = [ Y(6+1), Y(round(6+0.25*NNOX)),...
Y(round(6+0.5*NNOX)), Y(round(6+0.75*NNOX)), Y(6+NNOX) ]*to_ppm;
SAVE.NOe(II) = NOe_save*to_ppm; % save equilibrium NOx mass fraction for plotting
II=II+1;
if ( THETAS >= THETA && THETAS < THETAE )
Y(2) = tinitial( Y(1), Y(3), PHI, F );
end
if ( THETA > THETAS + THETAB )
Y(3) = nan;
end
end
fprintf('%7.1f %6.1f %3.3f %6.1f %6.1f %6.1f %5.0f %5.0f %5.3f %5.2f %6.1f\n', ...
THETA, VOL*1000000, X, Y(1), Y(2), Y(3), Y(4)*1000, Y(5)*1000, M*1000, Y(6)*1000, Y(7)*to_ppm );
end
% integrate total NOx value
NOX_ppm = 0;
for nn=1:NNOX;
THETA = THETAS + (nn-1)/(NNOX-1)*THETAB;
dxbdtheta = 0.5*sin(pi*(THETA-THETAS)/THETAB)*pi/THETAB;
dxb = dxbdtheta*DTHETA;
NOX_ppm = NOX_ppm + Y(6+nn)*dxb*to_ppm;
end
ETA = Y(4)/MNOT*(1+PHI*FS*(1-F))/PHI/FS/(1-F)/A0;
IMEP = Y(4)/(pi/4*B^2*S);
fprintf('ETA = %1.4f IMEP = %7.3f kPa NOx = %6.1f ppm\n', ETA, IMEP, NOX_ppm );
if ( nargin == 0 )
% if not called externally with custom PHI, F, and RPM parameters,
% generate some plots
sTitle = sprintf('Homogenous 2 zone, methane, PHI=%.2f F=%.2f RPM=%.1f\nETA=%.3f IMEP=%.2f kPa NOx=%.1f ppm ', PHI, F, RPM, ETA, IMEP, NOX_ppm );
figure;
plot( SAVE.THETA, SAVE.X, 'linewidth',2 );
set(gca,'fontsize',18,'linewidth',2,'Xlim',[-100 100]);
xlabel( '\theta','fontsize',18);
ylabel('burn fraction','fontsize',18);
print -deps2 NOxfrac
title( sTitle );
figure;
plot( SAVE.THETA, SAVE.P,'linewidth',2 );
set(gca,'fontsize',18,'linewidth',2,'Xlim',[-100 100]);
xlabel( '\theta','fontsize',18);
ylabel('pressure (kPa)','fontsize',18);
print -deps2 NOxpress
figure;
plot( SAVE.THETA, SAVE.TU, '-',SAVE.THETA, SAVE.TB,'--','linewidth',2 );
set(gca,'fontsize',18,'linewidth',2,'Xlim',[-100 100]);
xlabel( '\theta','fontsize',18);
ylabel( 'temperature (K)', 'fontsize',18);
legend('Unburned','Burned', 'Location', 'SouthEast');
figure;
plot( SAVE.THETA, SAVE.NOX, SAVE.THETA, SAVE.NOe,'linewidth',2 );
set(gca,'fontsize',18,'linewidth',2,'Xlim',[-100 100]);
xlabel('\theta','fontsize',18);
ylabel('NOx (ppm)','fontsize',18);
axis( [ THETAS, 110, 0, max(SAVE.NOe)*1.1 ] );
print -deps2 NOxNOx
legend( 'X=0', 'X=0.25', 'X=0.5', 'X=0.75', 'X=1', 'Equilibrium', 'Location', 'NorthEast' );
%title( sTitle );
end
function [ TB ] = tinitial( P, TU, PHI, F )
TB = 2000;
[~, HU,~, ~, ~, ~, ~, ~, ~, ~] = farg( TU, P, PHI, F, fuel_type );
for ITER=1:50,
[ierr, ~, HB,~, ~, ~, ~, CP, ~, ~, ~] = ecp( TB, P, PHI, fuel_type );
if ( ierr ~= 0 )
fprintf('Error in ECP(%g, %g, %g): %d\n', TB, P, PHI, ierr )
end
DELT = +(HU-HB)/CP;
TB = TB + DELT;
if ( abs(DELT/TB) < 0.001 )
break;
end
end
end
function [ VOL, X, EM ] = auxiliary( THETA )
VTDC = pi/4*B^2*S/(R-1); % m3
VOL = VTDC*(1 + (R-1)/2*(1-cosd(THETA) + 1/EPS*(1-sqrt(1-(EPS*sind(THETA))^2))));
X = 0.5*(1-cos(pi*(THETA-THETAS)/THETAB));
if ( THETA <= THETAS )
X = 0.;
end
if ( THETA >= THETAS+THETAB )
X = 1.;
end
EM = exp(-BLOWBY*(THETA*pi/180 + pi)/OMEGA);
end
function [Y] = integrate( THETA, THETAE, Y )
[TT, YY ] = ode23( @rates, [ THETA, THETAE ], Y );
for J=1:NY,
Y(J) = YY(length(TT),J);
end
function [ YPRIME ] = rates( THETA, Y )
YPRIME = zeros(NY,1);
[ VOL, X, EM ] = auxiliary( THETA );
M = EM*MNOT;
DUMB = sqrt(1-(EPS*sind(THETA))^2);
DV = pi/8*B^2*S*sind(THETA)*(1+EPS*cosd(THETA)/DUMB);
AA = (DV + VOL*BLOWBY/OMEGA)/M;
C1 = HEAT*(pi*B^2/2 + 4*VOL/B)/OMEGA/M/1000;
C0 = sqrt(X);
P = Y(1);
TB = Y(2);
TU = Y(3);
% three different computations are required depending upon the size
% of the mass fraction burned
if ( X > 0.999 )
% EXPANSION
[ierr, YB, HL, ~, ~, VB, ~, CP, ~, DVDT, DVDP] = ecp( TB, P, PHI, fuel_type );
if ( ierr ~= 0 )
fprintf('Error in ECP(%g, %g, %g): %d\n', TB, P, PHI, ierr );
end
BB = C1/CP*DVDT*TB*(1-TW/TB);
CC = 0;
DD = 1/CP*TB*DVDT^2 + DVDP;
EE = 0;
YPRIME(1) = (AA + BB + CC)/(DD + EE);
YPRIME(2) = -C1/CP*(TB-TW) + 1/CP*DVDT*TB*YPRIME(1);
YPRIME(3) = 0;
elseif ( X > 0.001 )
% COMBUSTION
[~, HU, ~, ~, VU, ~, CPU, ~, DVDTU, DVDPU] = farg( TU, P, PHI, F, fuel_type );
[ierr, YB, HB, ~, ~, VB, ~, CPB, ~, DVDTB, DVDPB] = ecp( TB, P, PHI, fuel_type );
if ( ierr ~= 0 )
fprintf('Error in ECP(%g, %g, %g): %d\n', TB, P, PHI, ierr );
end
BB = C1*(1/CPB*TB*DVDTB*C0*(1-TW/TB) + 1/CPU*TU*DVDTU*(1-C0)*(1-TW/TU));
DX = 0.5*sin( pi*(THETA-THETAS)/THETAB )*180/THETAB;
CC = -(VB-VU)*DX - DVDTB*(HU-HB)/CPB*(DX-(X-X^2)*BLOWBY/OMEGA);
DD = X*(VB^2/CPB/TB*(TB/VB*DVDTB)^2 + DVDPB);
EE = (1-X)*(1/CPU*TU*DVDTU^2 + DVDPU);
HL = (1-X^2)*HU + X^2*HB;
YPRIME(1) = (AA + BB + CC)/(DD + EE);
YPRIME(2) = -C1/CPB/C0*(TB-TW) + 1/CPB*TB*DVDTB*YPRIME(1) + (HU-HB)/CPB*(DX/X - (1-X)*BLOWBY/OMEGA);
YPRIME(3) = -C1/CPU/(1+C0)*(TU-TW) + 1/CPU*TU*DVDTU*YPRIME(1);
else
% COMPRESSION
[~, HL, ~, ~, ~, ~, CP, ~, DVDT, DVDP] = farg( TU, P, PHI, F, fuel_type );
BB = C1*1/CP*TU*DVDT*(1-TW/Y(3));
CC = 0;
DD = 0;
EE = 1/CP*TU*DVDT^2 + DVDP;
YPRIME(1) = ( AA + BB + CC )/(DD + EE);
YPRIME(2) = 0;
YPRIME(3) = -C1/CP*(Y(3)-TW) + 1/CP*Y(3)*DVDT*YPRIME(1);
end
% save equilibrium NO concentration
NOe_save = 0;
if ( X > 0.001 )
% save the equilibrium NO concentration as mass fraction
N_V = (1000*P)/(8.314*TB)*(1/100)^3; % calculate total molar concentration [mol/cm^3]
NOe_save = YB(10)*N_V*MW_NO*VB*1000; % NO mass fraction
end
% common to all cases
YPRIME(4) = Y(1)*DV;
YPRIME(5) = 0;
if ( ~isnan(TB) )
YPRIME(5) = YPRIME(5) + C1*M*C0*(TB-TW);
end
if ( ~isnan(TU) )
YPRIME(5) = YPRIME(5) + C1*M*(1-C0)*(TU-TW);
end
YPRIME(6) = BLOWBY*M/OMEGA*HL;
% perform NOx integration for each element burned
if ( X > 0.001 )
% COMBUSTION OR EXPANSION
for k=1:NNOX,
if ( THETA >= THETAS + (k-1)/(NNOX-1)*THETAB )
% convert Y(6+k) to [NO] mol/cm^3 from mass fraction
% and then back
YPRIME(6+k) = zeldovich( TB, P/100, YB, Y(6+k)/(MW_NO*VB*1000) )*MW_NO*VB*1000/OMEGA;
end
end
end
% 1/omega is s/rad, so convert to s/deg
for JJ=1:NY,
YPRIME(JJ) = YPRIME(JJ)*pi/180;
end
end
end
function [ dNOdt ] = zeldovich( T, P, y, NO )
% calculate rate of NO formation d[NO]/dt given
% inputs:
% T [K] : gas mixture temperature, kelvin
% P [bar] : cylinder pressure, bar
% y [...] : equilibrium mole fraction of constituents
% NO [mol/cm^3] : current NOx concentration
% outputs:
% dNOdt [ (mol/cm^3) / sec ] : rate of NO formation
% extended zeldovich rate constants from Heywood Table 11.1 (cm^3/mol-s)
k1 = 7.6*10^13*exp(-38000/T);
k2r = 1.5*10^9*T*exp(-19500/T);
k3r = 2*10^14*exp(-23650/T);
N_V = (100000*P)/(8.314*T)*(1/100)^3; % calculate molar concentration [mol/cm^3]
N2e = y(3)*N_V;
He = y(7)*N_V;
Oe = y(8)*N_V;
NOe = y(10)*N_V;
R1 = k1*Oe*N2e;
R2 = k2r*NOe*Oe;
R3 = k3r*NOe*He;
% ratios to check against Heywood Table 11.2
% at @ 2600 K, 10 bar, phi={0.8, 1, 1.2}
% R1
% R1R2 = R1/R2
% R1R2R3 = R1/(R2+R3)
alpha = NO/NOe;
dNOdt = 2*R1*(1-alpha*alpha)/(1+alpha*R1/(R2+R3));
end
end
bu kodu çalıştırmak istediğimde
Homogeneous(varargin)
Execution of script varargin as a function is not supported:
D:\MATLAB\toolbox\matlab\lang\varargin.m
bu hatayı alıyorum napacagımı çözemedim yardımcı olurmusunuz. teşekkürelr, saygılarımla
