|
|
chemical equilibrium equation
Posted:
May 11, 2012 4:05 AM
|
|
Dear sirs,
I cannt run this programme. can any one tell me where is the problem
masum
function [alpha,beta,gamma,delta,Afuel]=fueldata(fuel);
%
% [alpha,beta,gamma,delta,Afuel]=fueldata(fuel)
%
% Routine to specify the thermodynamic properties of a fuel.
% Data taken from:
% 1. Ferguson, C.R., 1986, "Internal Combustion Engines", Wiley;
% 2. Heywood, J.B., 1988, "Internal Combustion Engine Fundamentals",
% McGraw-Hill; and
% 3. Raine, R. R., 2000, "ISIS_319 User Manual", Oxford Engine Group.
% ********************************************************************
% input:
% fuel switch
% from Ferguson: 'gasoline', 'diesel', 'methane', 'methanol',
% 'nitromethane', 'benzene';
% from Heywood: 'methane_h', 'propane', 'hexane', 'isooctane_h',
% 'methanol_h', 'ethanol', 'gasoline_h1', gasoline_h2', 'diesel_h';
% from Raine: 'toluene', 'isooctane'.
% output:
% alpha, beta, gamma, delta - number of C, H, 0, and N atoms
% Afuel - vector of polynomial coefficients for cp/R, h/RT, and s/R
% of the form h/RT=a1+a2*T/2+a3*T^2/3+a4*T^3/4-a5/T^2+a6/T (for
% example) where T is expressed in K.
% ********************************************************************
% Set values for conversion of Heywood data to nondimensional format
% with T expressed in K
SVal=4.184e3/8.31434;
SVec=SVal*[1e-3 1e-6 1e-9 1e-12 1e3 1 1];
switch fuel
case 'gasoline' % Ferguson
alpha=7; beta=17; gamma=0; delta=0;
Afuel=[4.0652 6.0977E-02 -1.8801E-05 0 0 -3.5880E+04 15.45];
case 'diesel' % Ferguson
alpha=14.4; beta=24.9; gamma=0; delta=0;
Afuel=[7.9710 1.1954E-01 -3.6858E-05 0 0 -1.9385E+04 -1.7879];
case 'methane' % Ferguson
alpha=1; beta=4; gamma=0; delta=0;
Afuel=[1.971324 7.871586E-03 -1.048592E-06 0 0 -9.930422E+03 8.873728];
case 'methanol' % Ferguson
alpha=1; beta=4; gamma=1; delta=0;
Afuel=[1.779819 1.262503E-02 -3.624890E-06 0 0 -2.525420E+04 1.50884E+01];
case 'nitromethane' % Ferguson
alpha=1; beta=3; gamma=2; delta=1;
Afuel=[1.412633 2.087101E-02 -8.142134E-06 0 0 -1.026351E+04 1.917126E+01];
case 'benzene' % Ferguson
alpha=6; beta=6; gamma=0; delta=0;
Afuel=[-2.545087 4.79554E-02 -2.030765E-05 0 0 8.782234E+03 3.348825E+01];
case 'toluene' % Raine
alpha=7; beta=8; gamma=0; delta=0;
Afuel=[-2.09053 5.654331e-2 -2.350992e-5 0 0 4331.441411 34.55418257];
case 'isooctane' % Raine
alpha=8; beta=18; gamma=0; delta=0;
Afuel=[6.678E-1 8.398E-2 -3.334E-5 0 0 -3.058E+4 2.351E+1];
case 'methane_h' % Heywood
alpha=1; beta=4; gamma=0; delta=0;
Afuel=[-0.29149 26.327 -10.610 1.5656 0.16573 -18.331 19.9887/SVal].*SVec;
case 'propane' % Heywood
alpha=3; beta=8; gamma=0; delta=0;
Afuel=[-1.4867 74.339 -39.065 8.0543 0.01219 -27.313 26.4796/SVal].*SVec;
case 'hexane' % Heywood
alpha=6; beta=14; gamma=0; delta=0;
Afuel=[-20.777 210.48 -164.125 52.832 0.56635 -39.836 79.5542/SVal].*SVec;
case 'isooctane_h' % Heywood
alpha=8; beta=18; gamma=0; delta=0;
Afuel=[-0.55313 181.62 -97.787 20.402 -0.03095 -60.751 27.2162/SVal].*SVec;
case 'methanol_h' % Heywood
alpha=1; beta=4; gamma=1; delta=0;
Afuel=[-2.7059 44.168 -27.501 7.2193 0.20299 -48.288 31.1406/SVal].*SVec;
case 'ethanol' % Heywood
alpha=2; beta=6; gamma=1; delta=0;
Afuel=[6.990 39.741 -11.926 0 0 -60.214 8.01623/SVal].*SVec;
case 'gasoline_h1' % Heywood
alpha=8.26; beta=15.5; gamma=0; delta=0;
C farg.m 21
Afuel=[-24.078 256.63 -201.68 64.750 0.5808 -27.562 NaN].*SVec;
case 'gasoline_h2' % Heywood
alpha=7.76; beta=13.1; gamma=0; delta=0;
Afuel=[-22.501 227.99 -177.26 56.048 0.4845 -17.578 NaN].*SVec;
case 'diesel_h' % Heywood
alpha=10.8; beta=18.7; gamma=0; delta=0;
Afuel=[-9.1063 246.97 -143.74 32.329 0.0518 -50.128 NaN].*SVec;
end
function [h,u,v,s,Y,cp,dlvlT,dlvlp]=farg(p,T,phi,f,fueltype,airscheme);
%
% [h,u,v,s,Y,cp,dlvlT,dlvlp]=farg(p,T,phi,f,fueltype,airscheme)
%
% Routine to determine the state of mixtures of fuel, air
% and residual combustion products at low temperatures.
% Method closely follows that of:
% 1. Ferguson, C.R., 1986, "Internal Combustion Engines", Wiley, p108;
% who uses the results of:
% 2. Hires, S.D., Ekchian, A., Heywood, J.B., Tabaczynski, R.J., and
% Wall, J.C., 1976, "Performance and NOx Emissions Modeling of a Jet
% Ignition Pre-Chamber Stratified Charge Engine", SAE Trans., Vol 85,
% Paper 760161.
% ********************************************************************
% input:
% p,T,phi - pressure (Pa), temperature (K), and equivalence ratio
% f - residual mass fraction; set f=0 if no combustion products
% are present and f=1 if only combustion products are present
% fueltype - 'gasoline', 'diesel', etc - see fueldata.m for full list
% airscheme - 'GMcB' (Gordon and McBride) or 'Chemkin'
% output:
% h - enthalpy (J/kg), u - internal energy (J/kg),
% v - specific volume (m^3/kg), s - entropy (J/kgK),
% Y - mole fractions of 6 species: CO2, H2O, N2, O2, CO, and H2,
% cp - specific heat (J/kgK),
% dlvlT - partial derivative of log(v) wrt log(T)
% dlvlp - partial derivative of log(v) wrt log(p)
% ********************************************************************
[alpha,beta,gamma,delta,Afuel]=fueldata(fueltype);
switch airscheme
case 'GMcB'
A=airdata('GMcB_low');
case 'Chemkin'
A=airdata('Chemkin_low');
end
Ru=8314.34; % J/kmolK
table=[-1 1 0 0 1 -1]';
M=[44.01 18.02 28.008 32.000 28.01 2.018]'; % kg/kmol
MinMol=1e-25;
dlvlT=1; dlvlp=-1;
eps=0.210/(alpha+0.25*beta-0.5*gamma);
if phi <= 1.0 % stoichiometric or lean
nu=[alpha*phi*eps beta*phi*eps/2 0.79+delta*phi*eps/2 ...
0.21*(1-phi) 0 0]';
dcdT=0;
else % rich
z=1000/T;
K=exp(2.743+z*(-1.761+z*(-1.611+z*0.2803)));
dKdT=-K*(-1.761+z*(-3.222+z*0.8409))/1000;
a=1-K;
b=0.42-phi*eps*(2*alpha-gamma)+K*(0.42*(phi-1)+alpha*phi*eps);
c=-0.42*alpha*phi*eps*(phi-1)*K;
nu5=(-b+sqrt(b^2-4*a*c))/2/a;
dcdT=dKdT*(nu5^2-nu5*(0.42*(phi-1)+alpha*phi*eps)+ ...
0.42*alpha*phi*eps*(phi-1))/(2*nu5*a+b);
nu=[alpha*phi*eps-nu5 0.42-phi*eps*(2*alpha-gamma)+nu5 ...
0.79+delta*phi*eps/2 0 nu5 0.42*(phi-1)-nu5]';
end
% mole fractions and molecular weight of residual
tmoles=sum(nu);
Y=nu/tmoles;
Mres=sum(Y.*M);
% mole fractions and molecular weight of fuel-air
fuel=eps*phi/(1+eps*phi);
o2=0.21/(1+eps*phi);
n2=0.79/(1+eps*phi);
Mfa=fuel*(12.01*alpha+1.008*beta+16*gamma+14.01*delta)+ ...
32*o2+28.02*n2;
% mole fractions of fuel-air-residual gas
Yres=f/(f+Mres/Mfa*(1-f));
Y=Y*Yres;
Yfuel=fuel*(1-Yres);
Y(3)=Y(3)+n2*(1-Yres);
Y(4)=Y(4)+o2*(1-Yres);
% component properties
Tcp0=[1 T T^2 T^3 T^4]';
Th0=[1 T/2 T^2/3 T^3/4 T^4/5 1/T]';
Ts0=[log(T) T T^2/2 T^3/3 T^4/4 1]';
cp0=A(1:6,1:5)*Tcp0;
h0=A(1:6,1:6)*Th0;
s0=A(1:6,[1:5 7])*Ts0;
Mfuel=12.01*alpha+1.008*beta+16.000*gamma+14.01*delta;
a0=Afuel(1); b0=Afuel(2); c0=Afuel(3); d0=Afuel(6); e0=Afuel(7);
cpfuel=Afuel(1:5)*[1 T T^2 T^3 1/T^2]';
hfuel=Afuel(1:6)*[1 T/2 T^2/3 T^3/4 -1/T^2 1/T]';
s0fuel=Afuel([1:5 7])*[log(T) T T^2/2 T^3/3 -1/T^2/2 1]';
% set min value of composition so log calculations work
if Yfuel<MinMol
Yfuel=MinMol;
end
i=find(Y<MinMol);
Y(i)=ones(length(i),1)*MinMol;
% properties of mixture
h=hfuel*Yfuel+sum(h0.*Y);
s=(s0fuel-log(Yfuel))*Yfuel+sum((s0-log(Y)).*Y);
cp=cpfuel*Yfuel+sum(cp0.*Y)+sum(h0.*table*T*dcdT*Yres/tmoles);
MW=Mfuel*Yfuel+sum(Y.*M);
R=Ru/MW;
h=R*T*h;
u=h-R*T;
v=R*T/p;
s=R*(-log(p/101.325e3)+s);
cp=R*cp;
function [h,u,v,s,Y,cp,dlvlT,dlvlp]=ecp(p,T,phi,fueltype,airscheme,Yguess);
%
% [h,u,v,s,Y,cp,dlvlT,dlvlp]=ecp(p,T,phi,fueltype,airscheme,Yguess)
%
% Routine to determine the equilibrium state of combustion products.
% Method closely follows that of:
% 1. Ferguson, C.R., 1986, "Internal Combustion Engines", Wiley, p122;
% which uses the method described by:
% 2. Olikara, C., and Borman, G.L., 1975, "A Computer Program for
% Calculating Properties of Equilibrium Combustion Products with
% Some Applications to I.C. Engines", SAE Paper 750468.
% ********************************************************************
% input:
% p,T,phi - pressure (Pa), temperature (K), and equivalence ratio
% fueltype - 'gasoline', 'diesel', etc - see fueldata.m for full list
% airscheme - 'GMcB' (Gordon and McBride) or 'Chemkin'
% Yguess - (optional) initial estimate for mole fractions of the
% species CO2 H2O N2 O2 CO H2 H O OH and NO
% output:
% h - enthalpy (J/kg), u - internal energy (J/kg),
% v - specific volume (m^3/kg), s - entropy (J/kgK),
% Y - mole fractions of 10 species, cp - specific heat (J/kgK),
% dlvlT - partial derivative of log(v) wrt log(T)
% dlvlp - partial derivative of log(v) wrt log(p)
% ********************************************************************
[alpha,beta,gamma,delta,Afuel]=fueldata(fueltype);
switch airscheme
case 'GMcB'
A0=airdata('GMcB_hi');
case 'Chemkin'
A0=airdata('Chemkin_hi');
end
% Equilibrium constant data from Olikara and Borman via Ferguson
Kp=[ 0.432168E+00 -0.112464E+05 0.267269E+01 -0.745744E-04 0.242484E-08
0.310805E+00 -0.129540E+05 0.321779E+01 -0.738336E-04 0.344645E-08
-0.141784E+00 -0.213308E+04 0.853461E+00 0.355015E-04 -0.310227E-08
0.150879E-01 -0.470959E+04 0.646096E+00 0.272805E-05 -0.154444E-08
-0.752364E+00 0.124210E+05 -0.260286E+01 0.259556E-03 -0.162687E-07
-0.415302E-02 0.148627E+05 -0.475746E+01 0.124699E-03 -0.900227E-08];
MinMol=1e-25;
tol=3e-12;
Ru=8314.34; % J/kmol.K
M=[44.01 18.02 28.008 32.000 28.01 2.018 1.009 16 17.009 30.004]'; % kg/kmol
dcdT=zeros(4,1);
dcdp=zeros(4,1);
dfdT=zeros(4,1);
dfdp=zeros(4,1);
dYdT=zeros(10,1);
dYdp=zeros(10,1);
B=zeros(4,1);
% check if solid carbon will form
eps=0.210/(alpha+0.25*beta-0.5*gamma);
if phi>(0.210/eps/(0.5*alpha-0.5*gamma))
error('phi too high - c(s) and other species will form');
end
if nargin==5 % no Yguess so estimate the composition using farg
[h,u,v,s,Y,cp,dlvlT,dlvlp]=farg(p,T,phi,1,fueltype,airscheme);
Y(7:10)=ones(4,1)*MinMol; % since farg only returns first 6 species
else
Y=Yguess;
D ecp.m 25
end
% evaluate constants
patm=p/101.325e3; % convert Pa to atmospheres
TKp=[log(T/1000) 1/T 1 T T^2]';
K=10.^(Kp*TKp);
c=K.*[1/sqrt(patm) 1/sqrt(patm) 1 1 sqrt(patm) sqrt(patm)]';
d=[beta/alpha (gamma+0.42/eps/phi)/alpha (delta+1.58/eps/phi)/alpha]';
if abs(phi-1)<tol
phi=phi*(1+tol*sign(phi-1));
end
i=find(Y<MinMol);
Y(i)=ones(length(i),1)*MinMol;
DY3to6=2*tol*ones(4,1);
MaxIter=500;
MaxVal=max(abs(DY3to6));
Iter=0;
DoneSome=0;
while (Iter<MaxIter)&((MaxVal>tol)|(DoneSome<1))
Iter=Iter+1;
if Iter>2,
DoneSome=1;
end
D76=0.5*c(1)/sqrt(Y(6));
D84=0.5*c(2)/sqrt(Y(4));
D94=0.5*c(3)*sqrt(Y(6)/Y(4));
D96=0.5*c(3)*sqrt(Y(4)/Y(6));
D103=0.5*c(4)*sqrt(Y(4)/Y(3));
D104=0.5*c(4)*sqrt(Y(3)/Y(4));
D24=0.5*c(5)*Y(6)/sqrt(Y(4));
D26=c(5)*sqrt(Y(4));
D14=0.5*c(6)*Y(5)/sqrt(Y(4));
D15=c(6)*sqrt(Y(4));
A(1,1)=1+D103;
A(1,2)=D14+D24+1+D84+D104+D94;
A(1,3)=D15+1;
A(1,4)=D26+1+D76+D96;
A(2,1)=0;
A(2,2)=2*D24+D94-d(1)*D14;
A(2,3)=-d(1)*D15-d(1);
A(2,4)=2*D26+2+D76+D96;
A(3,1)=D103;
A(3,2)=2*D14+D24+2+D84+D94+D104-d(2)*D14;
A(3,3)=2*D15+1-d(2)*D15-d(2);
A(3,4)=D26+D96;
D ecp.m 26
A(4,1)=2+D103;
A(4,2)=D104-d(3)*D14;
A(4,3)=-d(3)*D15-d(3);
A(4,4)=0;
B(1)=-(sum(Y)-1);
B(2)=-(2*Y(2)+2*Y(6)+Y(7)+Y(9)-d(1)*Y(1)-d(1)*Y(5));
B(3)=-(2*Y(1)+Y(2)+2*Y(4)+Y(5)+Y(8)+Y(9)+Y(10)-d(2)*Y(1)-d(2)*Y(5));
B(4)=-(2*Y(3)+Y(10)-d(3)*Y(1)-d(3)*Y(5));
invA=inv(A);
DY3to6=invA*B;
MaxVal=max(abs(DY3to6));
Y(3:6)=Y(3:6)+DY3to6/10;
i=find(Y<MinMol);
Y(i)=ones(length(i),1)*MinMol;
Y(7)=c(1)*sqrt(Y(6));
Y(8)=c(2)*sqrt(Y(4));
Y(9)=c(3)*sqrt(Y(4)*Y(6));
Y(10)=c(4)*sqrt(Y(4)*Y(3));
Y(2)=c(5)*sqrt(Y(4))*Y(6);
Y(1)=c(6)*sqrt(Y(4))*Y(5);
end
if Iter>=MaxIter
warning('convergence failure in composition loop');
end
TdKdT=[1/T -1/T^2 1 2*T]';
dKdT=2.302585*K.*(Kp(:,[1 2 4 5])*TdKdT);
dcdT(1)=dKdT(1)/sqrt(patm);
dcdT(2)=dKdT(2)/sqrt(patm);
dcdT(3)=dKdT(3);
dcdT(4)=dKdT(4);
dcdT(5)=dKdT(5)*sqrt(patm);
dcdT(6)=dKdT(6)*sqrt(patm);
dcdp(1)=-0.5*c(1)/p;
dcdp(2)=-0.5*c(2)/p;
dcdp(5)=0.5*c(5)/p;
dcdp(6)=0.5*c(6)/p;
x1=Y(1)/c(6);
x2=Y(2)/c(5);
x7=Y(7)/c(1);
x8=Y(8)/c(2);
x9=Y(9)/c(3);
x10=Y(10)/c(4);
dfdT(1)=dcdT(6)*x1+dcdT(5)*x2+dcdT(1)*x7+dcdT(2)*x8+dcdT(3)*x9+dcdT(4)*x10;
dfdT(2)=2*dcdT(5)*x2+dcdT(1)*x7+dcdT(3)*x9-d(1)*dcdT(6)*x1;
dfdT(3)=2*dcdT(6)*x1+dcdT(5)*x2+dcdT(2)*x8+dcdT(3)*x9+dcdT(4)*x10-d(2)*dcdT(6)*x1;
dfdT(4)=dcdT(4)*x10-d(3)*dcdT(6)*x1;
dfdp(1)=dcdp(6)*x1+dcdp(5)*x2+dcdp(1)*x7+dcdp(2)*x8;
dfdp(2)=2*dcdp(5)*x2+dcdp(1)*x7-d(1)*dcdp(6)*x1;
dfdp(3)=2*dcdp(6)*x1+dcdp(5)*x2+dcdp(2)*x8-d(2)*dcdp(6)*x1;
dfdp(4)=-d(3)*dcdp(6)*x1;
B=-dfdT;
dYdT(3:6)=invA*B;
dYdT(1)=sqrt(Y(4))*Y(5)*dcdT(6)+D14*dYdT(4)+D15*dYdT(5);
dYdT(2)=sqrt(Y(4))*Y(6)*dcdT(5)+D24*dYdT(4)+D26*dYdT(6);
dYdT(7)=sqrt(Y(6))*dcdT(1)+D76*dYdT(6);
dYdT(8)=sqrt(Y(4))*dcdT(2)+D84*dYdT(4);
dYdT(9)=sqrt(Y(4)*Y(6))*dcdT(3)+D94*dYdT(4)+D96*dYdT(6);
dYdT(10)=sqrt(Y(4)*Y(3))*dcdT(4)+D104*dYdT(4)+D103*dYdT(3);
B=-dfdp;
dYdp(3:6)=invA*B;
dYdp(1)=sqrt(Y(4))*Y(5)*dcdp(6)+D14*dYdp(4)+D15*dYdp(5);
dYdp(2)=sqrt(Y(4))*Y(6)*dcdp(5)+D24*dYdp(4)+D26*dYdp(6);
dYdp(7)=sqrt(Y(6))*dcdp(1)+D76*dYdp(6);
dYdp(8)=sqrt(Y(4))*dcdp(2)+D84*dYdp(4);
dYdp(9)=D94*dYdp(4)+D96*dYdp(6);
dYdp(10)=D104*dYdp(4)+D103*dYdp(3);
% calculate thermodynamic properties
Tcp0=[1 T T^2 T^3 T^4]';
Th0=[1 T/2 T^2/3 T^3/4 T^4/5 1/T]';
Ts0=[log(T) T T^2/2 T^3/3 T^4/4 1]';
cp0=A0(:,1:5)*Tcp0;
h0=A0(:,1:6)*Th0;
s0=A0(:,[1:5 7])*Ts0;
% Y(1) and Y(2) reevaluated
Y(1)=(2*Y(3)+Y(10))/d(3)-Y(5);
Y(2)=(d(1)/d(3)*(2*Y(3)+Y(10))-2*Y(6)-Y(7)-Y(9))/2;
i=find(Y<MinMol);
Y(i)=ones(length(i),1)*MinMol;
% properties of mixture
h=sum(h0.*Y);
s=sum((s0-log(Y)).*Y);
cp=sum(Y.*cp0+h0.*dYdT*T);
MW=sum(Y.*M);
MT=sum(dYdT.*M);
Mp=sum(dYdp.*M);
R=Ru/MW;
v=R*T/p;
cp=R*(cp-h*T*MT/MW);
dlvlT=1+max(-T*MT/MW,0);
dlvlp=-1-max(p*Mp/MW,0);
h=R*T*h;
s=R*(-log(patm)+s);
u=h-R*T
function Tb=Tadiabatic(p,Tu,phi,f,fueltype,airscheme);
%
% Tb=Tadiabatic(p,Tu,phi,f,fueltype,airscheme)
%
% Routine for calculating the adiabatic flame temperature.
% Method involves iteratively selecting flame temperatures until
% the enthalpy of the combustion products (in equilibrium) matches
% the enthalpy of the initial gas mixture.
% farg.m is used to determine the enthalpy of the unburned mixture,
% and ecp.m is used to determine the enthalpy of the burned gas.
% ********************************************************************
% input:
% p - pressure (Pa)
% Tu - temperature of the unburned mixture (K)
% phi - equivalence ratio
% f - residual mass fraction; set f=0 if no combustion products
% are present and f=1 if only combustion products are present
% fueltype - 'gasoline', 'diesel', etc - see fueldata.m for full list
% airscheme - 'GMcB' (Gordon and McBride) or 'Chemkin'
% output:
% Tb - temperature of the burned gas (K) - adiabatic flame temperature
% ********************************************************************
MaxIter=50;
Tol=0.00001; % 0.001% allowable error in temperature calculation
Tb=2000; % initial estimate
DeltaT=2*Tol*Tb; % something big
Iter=0;
[hu,u,v,s,Y,cp,dlvlT,dlvlp]=farg(p,Tu,phi,f,fueltype,airscheme);
while (Iter<MaxIter)&(abs(DeltaT/Tb)>Tol)
Iter=Iter+1;
[hb,u,v,s,Y,cp,dlvlT,dlvlp]=ecp(p,Tb,phi,fueltype,airscheme);
DeltaT=(hu-hb)/cp;
Tb=Tb+DeltaT;
end
if Iter>=MaxIter
warning('convergence failure in adiabatic flame temperature loop');
end
% enginedata.m
%
% Script file used by the function ahrind.m to
% define the engine properties and initial conditions
% ***** engine geometry **********************************************
b=0.1; % engine bore (m)
stroke=0.08; % engine stroke (m)
eps=0.25; % half stroke to rod ratio, s/2l
r=10; % compression ratio
Vtdc=pi/4*b^2*stroke/(r-1); % volume at TDC
Vbdc=pi/4*b^2*stroke+Vtdc; % volume at BDC
% ***** engine thermofluids parameters *******************************
Cblowby=0.8; % piston blowby constant (s^-1)
f=0.1; % residual fraction
fueltype='gasoline';
airscheme='GMcB';
phi=0.8; % equivalence ratio
thetas=-35*pi/180; % start of burning
thetab=60*pi/180; % burn duration angle
RPM=2000;
omega=RPM*pi/30; % engine speed in rad/s
heattransferlaw='constant'; % 'constant', or 'Woschni'
hcu=500; % unburned zone heat transfer coefficient/weighting
hcb=500; % burned zone heat transfer coefficient/weighting
Tw=420; % engine surface temperature
% ***** initial conditions *******************************************
p1=100e3;
T1=350;
theta1=-pi;
V1=Vbdc;
[h1,u1,v1,s1,Y1,cp1,dlvlT1,dlvlp1]=farg(p1,T1,phi,f,fueltype,airscheme);
mass1=Vbdc/v1;
U1=u1*mass1;
function yprime=RatesComp(theta,y,flag);
%
% yprime=RatesComp(theta,y,flag)
%
% Function that returns the drivatives of the following 5 variables
% w.r.t. crank angle (theta) for the compression phase:
G RatesComp.m & RatesComb.m & RatesExp.m 30
% 1) pressure; 2) unburned temperature;
% 3) work; 4) heat transfer; and 5) heat leakage.
% See Ferguson, C.R., 1986, "Internal Combustion Engines", Wiley,
% p174.
global b stroke eps r Cblowby f fueltype airscheme phi ...
thetas thetab omega ...
heattransferlaw hcu ...
Tw theta1 Vtdc mass1 ...
%p=y(1);
Tu=y(2);
yprime=zeros(5,1);
% mass in cylinder accounting for blowby:
mass=mass1*exp(-Cblowby*(theta-theta1)/omega);
% volume of cylinder:
V=Vtdc*(1+(r-1)/2*(1-cos(theta)+ ...
1/eps*(1-(1-eps^2*sin(theta).^2).^0.5)));
% derivate of volume:
dVdtheta=Vtdc*(r-1)/2*(sin(theta)+ ...
eps/2*sin(2*theta)./sqrt(1-eps^2*sin(theta).^2));
switch heattransferlaw
case 'constant'
hcoeff=hcu;
case 'Woschni'
upmean=omega*stroke/pi; % mean piston velocity
C1=2.28;
hcoeff=hcu*130*b^(-0.2)*Tu^(-0.53)*(p/100e3)^(0.8)*C1*upmean;
end
A=1/mass*(dVdtheta+V*Cblowby/omega);
Qconv=hcoeff*(pi*b^2/2+4*V/b)*(Tu-Tw);
Const1=1/omega/mass;
[h,u,v,s,Y,cp,dlvlT,dlvlp]=farg(p,Tu,phi,f,fueltype,airscheme);
B=Const1*v/cp*dlvlT*Qconv/Tu; % note typo on p174, eq. 4.76
C=0;
D=0;
E=v^2/cp/Tu*dlvlT^2+v/p*dlvlp;
yprime(1)=(A+B+C)/(D+E);
yprime(2)=-Const1/cp*Qconv+v/cp*dlvlT*yprime(1);
yprime(3)=p*dVdtheta;
yprime(4)=Qconv/omega;
yprime(5)=Cblowby*mass/omega*h;
function yprime=RatesComb(theta,y,flag);
%
% yprime=RatesComb(theta,y,flag)
%
% Function that returns the drivatives of the following 6 variables
% w.r.t. crank angle (theta) for the combustion phase:
% 1) pressure; 2) unburned temperature; 3) burned temperature;
% 4) work; 5) heat transfer; and 6) heat leakage.
% See Ferguson, C.R., 1986, "Internal Combustion Engines", Wiley,
% p174.
global b stroke eps r Cblowby f fueltype airscheme phi ...
thetas thetab RPM omega ...
heattransferlaw hcu hcb ...
p1 T1 V1 Tw theta1 Vtdc Vbdc mass1
p=y(1);
Tb=y(2);
Tu=y(3);
yprime=zeros(6,1);
% mass in cylinder accounting for blowby:
mass=mass1*exp(-Cblowby*(theta-theta1)/omega);
% volume of cylinder:
V=Vtdc*(1+(r-1)/2*(1-cos(theta)+ ...
1/eps*(1-(1-eps^2*sin(theta).^2).^0.5)));
% derivate of volume:
dVdtheta=Vtdc*(r-1)/2*(sin(theta)+ ...
eps/2*sin(2*theta)./sqrt(1-eps^2*sin(theta).^2));
% mass fraction burned and derivative:
x=0.5*(1-cos(pi*(theta-thetas)/thetab));
dxdtheta=pi/2/thetab*sin(pi*(theta-thetas)/thetab);
if x<0.0001, x=0.0001; end;
if x>0.9999, x=0.9999; end;
switch heattransferlaw
case 'constant'
hcoeffu=hcu;
hcoeffb=hcb;
case 'Woschni'
upmean=omega*stroke/pi; % mean piston velocity
C1=2.28;
C2=3.24e-3;
Vs=Vbdc-Vtdc;
k=1.3;
pm=p1*(V1/V)^k; % motoring pressure
hcoeffu=hcu*130*b^(-0.2)*Tu^(-0.53)*(p/100e3)^(0.8)* ...
(C1*upmean+C2*Vs*T1/p1/V1*(p-pm))^(0.8);
hcoeffb=hcb*130*b^(-0.2)*Tb^(-0.53)*(p/100e3)^(0.8)* ...
(C1*upmean+C2*Vs*T1/p1/V1*(p-pm))^(0.8);
end
A=1/mass*(dVdtheta+V*Cblowby/omega);
Qconvu=hcoeffu*(pi*b^2/2+4*V/b)*(1-sqrt(x))*(Tu-Tw);
Qconvb=hcoeffb*(pi*b^2/2+4*V/b)*sqrt(x)*(Tb-Tw);
Const1=1/omega/mass;
[hu,uu,vu,s,Y,cpu,dlvlTu,dlvlpu]= ...
farg(p,Tu,phi,f,fueltype,airscheme);
[hb,ub,vb,s,Y,cpb,dlvlTb,dlvlpb]= ...
ecp(p,Tb,phi,fueltype,airscheme);
B=Const1*(vb/cpb*dlvlTb*Qconvb/Tb+ ...
vu/cpu*dlvlTu*Qconvu/Tu);
C=-(vb-vu)*dxdtheta-vb*dlvlTb*(hu-hb)/cpb/Tb*(dxdtheta- ...
(x-x^2)*Cblowby/omega);
D=x*(vb^2/cpb/Tb*dlvlTb^2+vb/p*dlvlpb);
E=(1-x)*(vu^2/cpu/Tu*dlvlTu^2+vu/p*dlvlpu);
yprime(1)=(A+B+C)/(D+E);
yprime(2)=-Const1/cpb/x*Qconvb+vb/cpb*dlvlTb*yprime(1)+ ...
(hu-hb)/cpb*(dxdtheta/x-(1-x)*Cblowby/omega);
yprime(3)=-Const1/cpu/(1-x)*Qconvu+vu/cpu*dlvlTu*yprime(1);
yprime(4)=p*dVdtheta;
yprime(5)=Const1*mass*(Qconvb+Qconvu);
yprime(6)=Cblowby*mass/omega*((1-x^2)*hu+x^2*hb);
function yprime=RatesExp(theta,y,flag);
%
% yprime=RatesExp(theta,y,flag)
%
% Function that returns the drivatives of the following 5 variables
% w.r.t. crank angle (theta) for the expansion phase:
% 1) pressure; 2) unburned temperature;
% 3) work; 4) heat transfer; and 5) heat leakage.
% See Ferguson, C.R., 1986, "Internal Combustion Engines", Wiley,
% p174.
global b stroke eps r Cblowby f fueltype airscheme phi ...
thetas thetab RPM omega ...
heattransferlaw hcb ...
p1 T1 V1 Tw theta1 Vtdc Vbdc mass1
p=y(1);
Tb=y(2);
yprime=zeros(5,1);
% mass in cylinder accounting for blowby:
mass=mass1*exp(-Cblowby*(theta-theta1)/omega);
% volume of cylinder:
V=Vtdc*(1+(r-1)/2*(1-cos(theta)+ ...
1/eps*(1-(1-eps^2*sin(theta).^2).^0.5)));
% derivate of volume:
dVdtheta=Vtdc*(r-1)/2*(sin(theta)+ ...
eps/2*sin(2*theta)./sqrt(1-eps^2*sin(theta).^2));
switch heattransferlaw
case 'constant'
hcoeff=hcb;
case 'Woschni'
upmean=omega*stroke/pi; % mean piston velocity
C1=2.28;
C2=3.24e-3;
Vs=Vbdc-Vtdc;
k=1.3;
pm=p1*(V1/V)^k; % motoring pressure
hcoeff=hcb*130*b^(-0.2)*Tb^(-0.53)*(p/100e3)^(0.8)* ...
(C1*upmean+C2*Vs*T1/p1/V1*(p-pm))^(0.8);
end
A=1/mass*(dVdtheta+V*Cblowby/omega);
Qconv=hcoeff*(pi*b^2/2+4*V/b)*(Tb-Tw);
Const1=1/omega/mass;
if Tb<1000
[h,u,v,s,Y,cp,dlvlT,dlvlp]=farg(p,Tb,phi,1,fueltype,airscheme);
else
[h,u,v,s,Y,cp,dlvlT,dlvlp]=ecp(p,Tb,phi,fueltype,airscheme);
end
B=Const1*v/cp*dlvlT*Qconv/Tb;
C=0;
D=v^2/cp/Tb*dlvlT^2+v/p*dlvlp;
E=0;
yprime(1)=(A+B+C)/(D+E);
yprime(2)=-Const1/cp*Qconv+v/cp*dlvlT*yprime(1);
yprime(3)=p*dVdtheta;
yprime(4)=Qconv/omega;
yprime(5)=Cblowby*mass/omega*h;
% ahrind.m
%
% Script file to determine the performance of a fuel inducted engine
% based on a (user-specified) arbitrary heat release profile as a
% function of crank angle.
% Method closely follows that of:
% Ferguson, C.R., 1986, "Internal Combustion Engines", Wiley.
% ********************************************************************
% input:
% enginedata.m - this is another script file that defines all of the
% relevant engine parameters and operating conditions.
% output:
% ahrind.mat - this file contains all of the variables. For plotting
% the results, see the example script file plotresults.m
% ********************************************************************
timestart=cputime;
global b stroke eps r Cblowby f fueltype airscheme phi ...
thetas thetab omega ...
heattransferlaw hcu hcb ...
Tw theta1 Vtdc Vbdc mass1 ...
p1 T1 V1
% load the engine parameters and initial conditions
enginedata
switch heattransferlaw
case 'Woschni'
if (abs(hcu) > 10)|(abs(hcb) > 10),
warning('Woschni model with weighting factor > 10')
end
end
% integration parameters
dtheta=1*pi/180;
options=odeset('RelTol',1e-3);
% integration during compression phase
disp(['integrating over the compression phase']);
[thetacomp,pTuWQlHl]=ode45('RatesComp', ...
[-pi:dtheta:thetas],[p1 T1 0 0 0],options);
% specification of initial conditions at start of combustion phase
% b - beginning of combustion
pb=interp1(thetacomp,pTuWQlHl(:,1),thetas);
Tub=interp1(thetacomp,pTuWQlHl(:,2),thetas);
Tbb=Tadiabatic(pb,Tub,phi,f,fueltype,airscheme);
Wb=interp1(thetacomp,pTuWQlHl(:,3),thetas);
Qlb=interp1(thetacomp,pTuWQlHl(:,4),thetas);
Hlb=interp1(thetacomp,pTuWQlHl(:,5),thetas);
% integration during combustion phase
disp(['integrating over the combustion phase']);
[thetacomb,pTbTuWQlHl]=ode45('RatesComb', ...
[thetas:dtheta:thetas+thetab],[pb Tbb Tub Wb Qlb Hlb],options);
% specification of initial conditions at start of expansion phase
% e - end of combustion / start of expansion
pe=interp1(thetacomb,pTbTuWQlHl(:,1),thetas+thetab);
Tbe=interp1(thetacomb,pTbTuWQlHl(:,2),thetas+thetab);
We=interp1(thetacomb,pTbTuWQlHl(:,4),thetas+thetab);
Qle=interp1(thetacomb,pTbTuWQlHl(:,5),thetas+thetab);
Hle=interp1(thetacomb,pTbTuWQlHl(:,6),thetas+thetab);
% integration during expansion phase
disp(['integrating over the expansion phase']);
[thetaexp,pTbWQlHl]=ode45('RatesExp', ...
[thetas+thetab:dtheta:pi],[pe Tbe We Qle Hle],options);
% error checks
mass4=mass1*exp(-Cblowby*2*pi/omega);
p4=interp1(thetaexp,pTbWQlHl(:,1),pi);
T4=interp1(thetaexp,pTbWQlHl(:,2),pi);
W4=interp1(thetaexp,pTbWQlHl(:,3),pi);
Ql4=interp1(thetaexp,pTbWQlHl(:,4),pi);
Hl4=interp1(thetaexp,pTbWQlHl(:,5),pi);
[h4,u4,v4,s4,Y4,cp4,dlvlT4,dlvlp4]= ...
farg(p4,T4,phi,1,fueltype,airscheme);
U4=u4*mass4;
error1=1-v4*mass4/Vbdc;
error2=1+W4/(U4-U1+Ql4+Hl4);
% indicated mean effective pressure and thermal efficiency
imep=W4/(pi*b^2/4*stroke);
eta=W4/mass1*(1+phi*0.06548*(1-f))/phi/0.06548/(1-f)/47870/1e3;
% calcuate the heat flux in W/m^2
qcomp=calcq(thetacomp,pTuWQlHl,'comp'); % compression
qcombu=calcq(thetacomb,pTbTuWQlHl,'combu'); % combustion-unburned zone
qcombb=calcq(thetacomb,pTbTuWQlHl,'combb'); % combustion-burned zone
qexp=calcq(thetaexp,pTbWQlHl,'exp'); % expansion
timefinish=cputime;
timetaken=timefinish-timestart;
% save all data
save ahrind.mat
clear
I plotresults.m
% plotresults.m
%
I plotresults.m 36
% Script file to plot output from ahrind.m and compare results
% with output listed in:
% Ferguson, C.R., 1986, "Internal Combustion Engines", Wiley, p178;
load ahrind.mat; % results from ahrind.mat
load ferguson.txt; % tabulated output from ferguson p178-179
% Set some parameters to make figures look attractive
NLW=1; % normal line width
NFS=18; % normal font size
NMS=1; % normal marker size
close all;
figure(1);
plot(thetacomp*180/pi,pTuWQlHl(:,1)/1e6); hold on;
plot(thetacomb*180/pi,pTbTuWQlHl(:,1)/1e6);
plot(thetaexp*180/pi,pTbWQlHl(:,1)/1e6);
plot(ferguson(:,1),ferguson(:,4)*1e5/1e6,'o');
axis([-180 180 0 7]);
set(gca,'FontSize',NFS)
set(gca,'LineWidth',NLW)
set(gca,'XTick',[-180 -90 0 90 180]);
set(gca,'XTickLabel',[-180 -90 0 90 180]);
xlabel('crank angle (degrees ATC)')
ylabel('pressure (MPa)')
print -deps p_ahr.eps
figure(2);
plot(ferguson(:,1),ferguson(:,5),'o'); hold on;
plot(ferguson(:,1),ferguson(:,6),'s');
plot(thetacomp*180/pi,pTuWQlHl(:,2));
plot(thetacomb*180/pi,pTbTuWQlHl(:,3));
plot(thetacomb*180/pi,pTbTuWQlHl(:,2));
plot(thetaexp*180/pi,pTbWQlHl(:,2));
axis([-180 180 0 3000]);
set(gca,'FontSize',NFS)
set(gca,'LineWidth',NLW)
set(gca,'XTick',[-180 -90 0 90 180]);
set(gca,'XTickLabel',[-180 -90 0 90 180]);
xlabel('crank angle (degrees ATC)')
ylabel('temperature (K)')
legend('burned gas','unburned gas',2);
print -deps T_ahr.eps
figure(3);
nn=4; % for work plot
plot(thetacomp*180/pi,pTuWQlHl(:,nn-1)); hold on;
plot(thetacomb*180/pi,pTbTuWQlHl(:,nn));
plot(thetaexp*180/pi,pTbWQlHl(:,nn-1));
I plotresults.m 37
plot(ferguson(:,1),ferguson(:,7),'o');
axis([-180 180 -300 600]);
set(gca,'FontSize',NFS)
set(gca,'LineWidth',NLW)
set(gca,'XTick',[-180 -90 0 90 180]);
set(gca,'XTickLabel',[-180 -90 0 90 180]);
set(gca,'YTick',[-300 -150 0 150 300 450 600]);
set(gca,'YTickLabel',[-300 -150 0 150 300 450 600]);
xlabel('crank angle (degrees ATC)')
ylabel('work (J)')
print -deps W_ahr.eps
figure(4);
nn=5; % for heat transfer plot
plot(thetacomp*180/pi,pTuWQlHl(:,nn-1)); hold on;
plot(thetacomb*180/pi,pTbTuWQlHl(:,nn));
plot(thetaexp*180/pi,pTbWQlHl(:,nn-1));
plot(ferguson(:,1),ferguson(:,8),'o');
axis([-180 180 -50 300]);
set(gca,'FontSize',NFS)
set(gca,'LineWidth',NLW)
set(gca,'XTick',[-180 -90 0 90 180]);
set(gca,'XTickLabel',[-180 -90 0 90 180]);
xlabel('crank angle (degrees ATC)')
ylabel('heat transfer (J)')
print -deps Ql_ahr.eps
figure(5);
nn=6; % for heat leakage plot
plot(thetacomp*180/pi,pTuWQlHl(:,nn-1)); hold on;
plot(thetacomb*180/pi,pTbTuWQlHl(:,nn));
plot(thetaexp*180/pi,pTbWQlHl(:,nn-1));
plot(ferguson(:,1),ferguson(:,10),'o');
axis([-180 180 -8 0]);
set(gca,'FontSize',NFS)
set(gca,'LineWidth',NLW)
set(gca,'XTick',[-180 -90 0 90 180]);
set(gca,'XTickLabel',[-180 -90 0 90 180]);
xlabel('crank angle (degrees ATC)')
ylabel('heat leakage (J)')
print -deps Hl_ahr.eps
figure(6);
doflamearrivalplot=0;
thetaflame=-5*pi/180;
lengththeta=length(thetacomb);
plot(thetacomp*180/pi,qcomp); hold on;
if doflamearrivalplot==1,
nflameb=min(find(thetacomb>thetaflame));
nflameu=nflameb-1;
plot(thetacomb(nflameu:nflameb)*180/pi, ...
[qcombu(nflameu) qcombb(nflameb)]);
else
nflameb=1;
nflameu=length(thetacomb);
end
Nunburned=1:nflameu;
Nburned=nflameb:lengththeta;
plot(thetacomb(Nunburned)*180/pi,qcombu(Nunburned));
plot(thetacomb(Nburned)*180/pi,qcombb(Nburned));
plot(thetaexp*180/pi,qexp);
axis([-180 180 -0.2e6 1.2e6]);
set(gca,'FontSize',NFS)
set(gca,'LineWidth',NLW)
set(gca,'XTick',[-180 -90 0 90 180]);
set(gca,'XTickLabel',[-180 -90 0 90 180]);
xlabel('crank angle (degrees ATC)')
ylabel('heat flux (W/m^2)')
print -deps qflux_ahr.eps
function q=calcq(theta,pTarray,phase);
%
% calculation of the heat flux (W/m^2) from the data generated
% by ahrind.
% theta is an array of crank angles
% pTarray is the corresponding array of pressure, Temperature,
% Work, etc data as generated by running arhind.m
% phase is a switch indicating the part of the cycle:
% 'comp' - compression phase
% 'combu' - combustion phase, unburned gas zone
% 'combb' - combustion phase, burned gas zone
% 'exp' - expansion phase
global b stroke eps r ...
omega ...
heattransferlaw hcu hcb ...
Tw Vtdc Vbdc ...
p1 T1 V1
switch phase
case 'comp'
p=pTarray(:,1);
T=pTarray(:,2);
hc=hcu;
C2=0;
case 'combu'
K ferguson.txt 39
p=pTarray(:,1);
T=pTarray(:,3);
hc=hcu;
C2=3.24e-3;
case 'combb'
p=pTarray(:,1);
T=pTarray(:,2);
hc=hcb;
C2=3.24e-3;
case 'exp'
p=pTarray(:,1);
T=pTarray(:,2);
hc=hcb;
C2=3.24e-3;
end
switch heattransferlaw
case 'constant'
hcoeff=hc;
case 'Woschni'
V=Vtdc*(1+(r-1)/2*(1-cos(theta)+ ...
1/eps*(1-(1-eps^2*sin(theta).^2).^0.5)));
upmean=omega*stroke/pi; % mean piston velocity
Vs=Vbdc-Vtdc;
k=1.3;
C1=2.28;
pm=p1*(V1./V).^k; % motoring pressure
hcoeff=hc*130*b^(-0.2)*T.^(-0.53).*(p/100e3).^(0.8).* ...
(C1*upmean+C2*Vs*T1/p1/V1*(p-pm)).^(0.8);
end
q=hcoeff.*(T-Tw);
|
|
