Adiabatic Flame Temperature (C2H2) [MATLAB Code] (Extensible, by adding pressure effect..)

MATLAB CODE AUTO: clear all close all clc % c=input('Enter number of C in fuel: '); % h=input('Enter number of H in fuel: '); c=2; h=2; o=c+h/4; %yH yH2 yO yO2 yOH yH2O yCO yCO2 y= [10.44 16 10.53 17.8 16.23 22.17 17.05 27.67;... 11.44 17.71 11.52 19.66 17.95 24.94 18.81 30.61;... 12.43 19.45 12.52 21.51 19.7 27.55 20.58 33.56;... 13.42 21.21 13.52 23.45 21.46 30.2 22.36 36.53;... 14.42 22.99 14.52 25.36 23.23 32.47 24.14 39.51;... 15.41 24.79 15.53 27.3 25.03 35.57 25.93 42.50;... 16.4 26.62 16.54 29.25 26.83 38.3 27.72 45.51;... 17.4 28.16 17.55 31.22 28.65 41.04 29.51 48.52;... 18.39 30.32 18.57 33.2 30.49 43.8 31.32 51.51]; disp(['C',num2str(c),'H',num2str(h),' +(',num2str(o),')O2==>(',num2str(c),')CO2+(',num2str(h/2),')H2O']) T0=input('Enter Initial Temperature[K] of combustion: '); e=input('Enter maximum error: '); mov=input('Enter 1 to see errors plot, else enter 0: '); q2=0; for T=T0:200:4000 q1=q2; % if T==2300 % K1=0.01155; % K2=0.00205; % K3=0.00157; % K4=0.00963; % K5=0.00493; if T==2400 K1=0.0211; K2=0.00353; K3=0.00294; K4=0.0158; K5=0.00861; yH=y(1,1); yH2=y(1,2); yO=y(1,3); yO2=y(1,4);

yOH=y(1,5); yH2O=y(1,6); yCO=y(1,7); yCO2=y(1,8); % elseif T==2500 % K1=0.0366; % K2=0.00585; % K3=0.00523; % K4=0.0250; % K5=0.0144; elseif T==2600 K1=0.0610; K2=0.00933; K3=0.00890; K4=0.0382; K5=0.0231; yH=y(2,1); yH2=y(2,2); yO=y(2,3); yO2=y(2,4); yOH=y(2,5); yH2O=y(2,6); yCO=y(2,7); yCO2=y(2,8); % elseif T==2700 % K1=0.0974; % K2=0.01434; % K3=0.01467; % K4=0.0565; % K5=0.0359; elseif T==2800 K1=0.151; K2=0.0214; K3=0.0230; K4=0.0814; K5=0.0539; yH=y(3,1); yH2=y(3,2); yO=y(3,3); yO2=y(3,4); yOH=y(3,5); yH2O=y(3,6); yCO=y(3,7); yCO2=y(3,8); % elseif T==2900

% K1=0.226; % K2=0.0311; % K3=0.0353; % K4=0.114; % K5=0.0788; elseif T==3000 K1=0.330; K2=0.0441; K3=0.0525; K4=0.157; K5=0.1125; yH=y(4,1); yH2=y(4,2); yO=y(4,3); yO2=y(4,4); yOH=y(4,5); yH2O=y(4,6); yCO=y(4,7); yCO2=y(4,8); % elseif T==3100 % K1=0.469; % K2=0.0610; % K3=0.0721; % K4=0.223; % K5=0.1567; elseif T==3200 K1=0.653; K2=0.0829; K3=0.1078; K4=0.279; K5=0.2138; yH=y(5,1); yH2=y(5,2); yO=y(5,3); yO2=y(5,4); yOH=y(5,5); yH2O=y(5,6); yCO=y(5,7); yCO2=y(5,8); elseif T==3400 K1=1.191; K2=0.1447; K3=0.203; K4=0.465; K5=0.377;

yH=y(6,1); yH2=y(6,2); yO=y(6,3); yO2=y(6,4); yOH=y(6,5); yH2O=y(6,6); yCO=y(6,7); yCO2=y(6,8); elseif T==3600 K1=2.025; K2=0.237; K3=0.356; K4=0.732; K5=0.625; yH=y(7,1); yH2=y(7,2); yO=y(7,3); yO2=y(7,4); yOH=y(7,5); yH2O=y(7,6); yCO=y(7,7); yCO2=y(7,8); elseif T==3800 K1=3.25; K2=0.369; K3=0.588; K4=1.099; K5=0.984; yH=y(8,1); yH2=y(8,2); yO=y(8,3); yO2=y(8,4); yOH=y(8,5); yH2O=y(8,6); yCO=y(8,7); yCO2=y(8,8); elseif T==4000 K1=4.97; K2=0.549; K3=0.924; K4=1.585; K5=1.477; yH=y(9,1); yH2=y(9,2); yO=y(9,3);

yO2=y(9,4); yOH=y(9,5); yH2O=y(9,6); yCO=y(9,7); yCO2=y(9,8); end % f=PCO2/PCO; errorP=0; errorn=0; if mov==1 M=moviein(12); z=0; end for i=0.01:0.01:4 f=i; if mov==1 clf end for j=0.001:0.001:4 PH2O=j; PO2=double((K1*f)^2); PO=double(K5*sqrt(PO2)); PH2=double((K2*PH2O)/sqrt(PO2)); POH=double((K3*PH2O)/sqrt(PH2)); PH=double(K4*sqrt(PH2)); nH=2*PH2+PH+POH+2*PH2O; nC=nH*(c/h); PCO=double(nC/(1+f)); PCO2=nC-PCO; nO=PH2O+2*PO2+PO+POH+PCO+2*PCO2; P=PH2O+PO2+PO+PH2+PH+POH+PCO+PCO2; errorP=P-1; errorn=nO/nH-(2*o)/h; if mov==1 z=z+1; plot(j,abs(errorP),'b-o') hold on plot(j,abs(errorn),'r-o') plot(j,e,'g-') title(['PCO2/PCO=',num2str(i)]) xlabel(['PH2O=',num2str(PH2O)]) ylabel(['Errors | Max-Error=',num2str(e)]) legend('P-error','nO/nH-error','maximum error line','Location','NorthEastOutside') M(:,z)=getframe;

end if mov==1 if (abs(errorP)>10*e || abs(errorn)>10*e) && PH2O>=2 break end end if abs(errorP)