zeldovich
function [ 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
ve
run ecp
clear;
phi = 0.9; % enter equivalence ratio input
T = 800; % enter temperature (K) input
P = 100; % enter pressure (kPa) input
fuel_id =1;
% fuel_id - 1=Methane, 2=Gasoline, 3=Diesel, 4=Methanol, 5=Nitromethane
% call ecp function
[ierr, Y, h, u, s, v, R, Cp, MW, dvdT, dvdP] = ecp( T, P, phi, fuel_id );
%echo input
fprintf(' \n Equilibrium Combustion Solver \n' );
fprintf(' Pressure (kPa) = \t \t %6.1f \n', P );
fprintf(' Temperature (K) = \t \t %6.1f \n', T); ...
fprintf(' Fuel Air Equivalence ratio = \t% 4.2f \n ', phi);
%print output mole fractions and properties
fprintf(' \n Mole Fractions \n' );
fprintf(' CO2 = \t %1.4e \n', Y(1) );
fprintf(' H2O = \t %1.4e \n', Y(2) );
fprintf(' N2 = \t %1.4e \n', Y(3) );
fprintf(' O2 = \t %1.4e \n', Y(4) );
fprintf(' CO = \t %1.4e \n', Y(5) );
fprintf(' H2 = \t %1.4e \n', Y(6) );
fprintf(' H = \t %1.4e \n', Y(7) );
fprintf(' O = \t %1.4e \n', Y(8) );
fprintf(' OH = \t %1.4e \n', Y(9) );
fprintf(' NO = \t %1.4e \n', Y(10) );
fprintf(' \n Mixture Properties \n' );
fprintf(' h(kJ/kg) = \t %6.1f \n', h );
fprintf(' u(kJ/kg) = \t %6.1f \n', u );
fprintf(' s (kJ/Kg K) = \t %6.4f \n', s );
fprintf(' v (m3/kg) = \t %6.4f \n', v );
fprintf(' cp (kJ/Kg K) =\t %6.3f \n', Cp );
fprintf(' Molecular Mass = %5.2f \n', MW );
fprintf(' dvdt = %1.4e \n', dvdT );
fprintf(' dvdp = %1.4e \n', dvdP );
şu iki kodu birleştirirmem lazım fakat run ecp yi zeldovic e entegre etmeliyim çıktılarda farklı thetha açılarında CO C Nox olmalı fakat bi türlü yapamadım hocam yapamazsan tezin yarım kalır dedi kaç gündür uğraşıyorum bir türlü olmadı yardımcı olabilecek en azından yol gösterebilcek olan biri varsa çok minnettar olurum şimdiden çok teşekkürler saygılarımla…