Got Questions? Get Answers.
Discover MakerZone

MATLAB and Simulink resources for Arduino, LEGO, and Raspberry Pi

Learn more

Discover what MATLAB® can do for your career.

Opportunities for recent engineering grads.

Apply Today

Thread Subject:
chemical equilibrium equation

Subject: chemical equilibrium equation

From: masum billah

Date: 11 May, 2012 08:05:10

Message: 1 of 2

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);

Subject: chemical equilibrium equation

From: Steven_Lord

Date: 11 May, 2012 14:10:43

Message: 2 of 2



"masum billah" <masum05me@gmail.com> wrote in message
news:joih7m$k6j$1@newscl01ah.mathworks.com...
> Dear sirs,
> I cannt run this programme. can any one tell me where is the problem
> masum

*snip nearly 900 lines of code*

I think it unlikely that anyone is going to sit and read through pages and
pages of source code to help you locate the problem when you haven't showed
the group how you're calling the function or what you mean by "cannt[sic]
run this programme" -- does it error? Throw a warning? Give an answer you
didn't expect?

But even if you provided that information, I think you're going to be the
best suited to determine what's wrong, since you know best what you're
trying to do. So I'm going to point you to the tools you'll need to
determine where the problem is on your own. If you're receiving an error or
warning, set an error breakpoint as described in the last section of this
page from the documentation:

http://www.mathworks.com/help/techdoc/matlab_env/brqxeeu-175.html

Then run your code. MATLAB will stop execution on the line where the
error/warning occurred, so you can investigate what went wrong in the
function call on that line. You can also use regular breakpoints to step
through your code if you're receiving an unexpected answer to determine
where the answer MATLAB computes starts to differ from the one you expected.

I would also review the Code Analyzer messages for your code:

http://www.mathworks.com/help/techdoc/matlab_env/brqxeeu-151.html#brqxeeu-155

Problems that will probably prevent your code from running to completion are
generally indicated by red bars, while areas where you may be able to
improve the efficiency, performance, robustness, or memory consumption of
your code are indicated by orange bars.

If you're able to isolate the problem down to a small section of code (no
more than a dozen lines) but can't narrow it down any further, feel free to
post that small section of code along with a description of the problem
you're experiencing to the newsgroup. When you do that, someone may be able
to help diagnose the problem and determine how to fix it.

--
Steve Lord
slord@mathworks.com
To contact Technical Support use the Contact Us link on
http://www.mathworks.com

Tags for this Thread

No tags are associated with this thread.

What are tags?

A tag is like a keyword or category label associated with each thread. Tags make it easier for you to find threads of interest.

Anyone can tag a thread. Tags are public and visible to everyone.

Contact us