I" m trying to build function that solves pyrolysis problem. But, when I run it, it keep saying "Index exceeds number of array elements (1)." I need help…
Error is in this funcation on line 53.
Error in nonstqsol (line 53)
%% Function of non-stoichiometric thermodynamic equilibrium model
function fun=nonstqsol(x)global n M_biomass T M_gas M_char n_charR=8.314;a=0.2;b=0.13;x=1.44;y=0.66;if nargin==0 x=ones(1,8);end% To calculate the Gibbs formation Value of parameter a-g taken from
% Synthetic fuels, by Probstein, R F; Hicks, R E, 1982.
% H2 CO, CO2, H2O, CH4, C, O2
prm=[0 5.619e-3 -1.949e-2 -8.950e-3 -4.620e-2 0 0;... %a 0 -1.190e-5 3.122e-5 -3.672e-6 1.130-5 0 0;... %b 0 6.383e-9 -2.448e-8 5.209e-9 1.319e-8 0 0;... %c 0 -1.846e-12 6.946e-12 -1.478e-12 -6.647e-12 0 0;... %d 0 -4.891e2 -4.891e2 0.0 -4.891e2 0 0;... %e 0 8.684e-1 5.270 2.868 1.411e1 0 0;... %f 0 -6.131e-2 -1.207e-1 -1.722e-2 -2.234e-1 0 0];... %g%% Enthalpy of formation at temperature of 298K
% [H2 CO CO2 H2O CH4 C O2]
h_f298K=[0 -110.5 -393.5 -241.8 -74.8 0 0];%% Standard Gibbs function of formation at given temperature of the gas species i.
%Standard Gibbs of CO
del_gCO=h_f298K(1,2)-(prm(1,2).*T(n).*log(T(n))-(prm(2,2).*(T(n))^2)-(prm(3,2)/2.*(T(n))^3)... -(prm(4,2)/3.*(T(n))^4)+(prm(5,2)/(2.*T(n))+prm(6,2)+(prm(7,2).*T(n))));%Standard Gibbs of CO2
del_gCO2=h_f298K(1,3)-(prm(1,3).*T(n).*log(T(n))-(prm(2,3).*(T(n))^2)-(prm(3,3)/2.*(T(n))^3)... -(prm(4,3)/3.*(T(n))^4)+(prm(5,3)/(2.*T(n))+prm(6,3)+(prm(7,3).*T(n))));%Standard Gibbs of H2O
del_gH2O=h_f298K(1,4)-(prm(1,4).*T(n).*log(T(n))-(prm(2,4).*(T(n))^2)-(prm(3,4)/2.*(T(n))^3)... -(prm(4,4)/3.*(T(n))^4)+(prm(5,4)/(2.*T(n))+prm(6,4)+(prm(7,4).*T(n))));%Standard Gibbs of CH4
del_gCH4=h_f298K(1,5)-(prm(1,5).*T(n).*log(T(n))-(prm(2,5).*(T(n))^2)-(prm(3,5)/2.*(T(n))^3)... -(prm(4,5)/3.*(T(n))^4)+(prm(5,5)/(2.*T(n))+prm(6,5)+(prm(7,5).*T(n))));del_gspecies=[0; del_gCO; del_gCO2; del_gH2O; del_gCH4; 0; 0];% Stoichiometric Co-efficient Reactions
% H2 CO CO2 H2O CH4 C O2
vi = [1, 1, 0, -1, 0, -1, 0;... % C + H2O --> CO + H21,-1, 1, -1, 0, 0, 0;... % CO + H2O --> CO2+ H2-2, 0, 0, 0, 1, -1, 0;... % C + 2H2 --> CH40, 2,-1, 0, 0, -1, 0;... % C + CO2 --> 2CO3, 1, 0, -1, -1, 0, 0;... % CH4+ H2O --> CO + 3H20, 2, 0, 0, 0, -2,-1;... % 2C + O2 --> 2CO2, 0, 1, -2, 0, -1, 0]; % C + 2H2O--> CO2+ 2H2
%% Gibbs enthalpy of formation of species i
%G_Ft=v_i*g_(f,T,i)
GFt(n,:)=vi*del_gspecies;%% mass of components H2, CO, CO2, H2O, CH4, CO2
nH2 =x(1); %Mole of H2O
nCO =x(2); %Mole of carbon monoxide
nCO2 =x(3); %Mole of carbon dioxide
nH2O =x(4); %Mole of hydrogen
nCH4 =x(5); %Mole of methane
Lg_C =x(6); %La grangian for Carbon
Lg_O =x(7); %La grangian for Oxygen
Lg_H =x(8); %La grangian for Hydrogen
% The Composition of Wood in weight percent.
% C - 51.6% H - 6.2% O - 43.2% CH(1.44)O(0.66)
Mw_Char = 14.2; %Molecular weight of char
Mw_biomass = 24; %Molecular weight of biomass
n_biomass= M_biomass/Mw_biomass; %Mole of biomass
n_char =M_char/Mw_Char; %Mole of char
n_gas=n_biomass-n_char; %Mole of gass
Mw_gas= M_gas/n_gas; %Molecular weight of gas
%%
%Carbon mass balance
f(1)=n_biomass-(n_gas*(nCO+nCO2+nCH4+n_char));%Hydrogen Mass Balance
f(2)=(x*n_biomass)-(n_gas*(2*nH2+4*nCH4+2*nH2O)+a*n_char);%Oxygen Mass Balance
f(3)=(y*n_biomass)-(n_gas*(nCO+2*nCO2+nH2O)+b*n_char);%H2 Gibbs function
f(4)=GFt(n,1)+R*T(n)*log(nH2/sum(x(1:5)+n_char))-2*Lg_H;%CO Gibbs
f(5)=GFt(n,2)+R*T(n)*log(nCO/(sum(x(1:5)+n_char)))-Lg_C-Lg_O;%CO2 Gibbs
f(6)=GFt(n,3)+R*T(n)*log(nCO2/(sum(x(1:5)+n_char)))-Lg_C-2*Lg_O;%H2O gibbs
f(7)=GFt(n,4)+R*T(n)*log(nH2O/(sum(x(1:5)+n_char)))-Lg_H-2*Lg_O;%CH4 Gibbs
f(8)=GFt(n,5)+R*T(n)*log(nCH4/(sum(x(1:5)+n_char)))-Lg_C-4*Lg_H;fun = [f(1),f(2),f(3),f(4),f(5),f(6),f(7),f(8)];Solver function
% Calling function
function nlinear_eq_solverclose allclear allclcglobal n M_biomass T M_char M_gas n_charT = 953:5:1053;%input('Enter the Reactor temperature in the Range between 700-1165 [K] :');
M_biomass = 15;%input('Enter the Biomass feed rate in [Kg/hr] :')*1000;%kg
Matrix =[0.41621 0.4973 0.1801 0.0656 0.0051 11508 202070 2820];for n=1:length(T)% Pyrolysis from "Modeling circulating fluidized bed biomass gasifiers.
[t1,z] =ode113(@Tinitialpyrolysis,[0 1],[1 0 0 0]);% Material balance
M_tar = M_biomass*z(end,3); % Rate of tar into the reactor by devolatilization, [kg/s]
M_char = M_biomass*z(end,4); % Rate of char into the reactor by devolatilization, [kg/s]
M_gas=M_biomass*z(end,2); % [Kg/s]
InitialGuess = Matrix;options=optimset('Display','iter','TolFun',1e-32,'MaxFunEvals',10000,'MaxIter',10000);[x,fval,exitflag] = fsolve(@nonstqsol,InitialGuess,options);Eflag(n)=[exitflag];yh2=x(1)/sum(x(1:5)); %mol frac hydrogen
yco=x(2)/sum(x(1:5)); %mol frac carbon monoxide,exitflag
yco2=x(3)/sum(x(1:5)); %mol frac carbon dioxide
yh2o=x(4)/sum(x(1:5)); %mol frac h2o
ych4=x(5)/sum(x(1:5)); %mol frac of methane
Lg_C=x(6) %La grangian for CarbonLg_O=x(7); %La grangian for OxygenLg_H=x(8); %La grangian for Hydrogenfn(n,:)=fval;mol(n,:) =[x(1:5) n_char];molfr(n,:) =[yh2 yco yco2 yh2o ych4];Matrix=x;endEflagfnmolfigure(1)hold ongrid onbox onH=plot(T-273,molfr(:,1),'.r');K=plot(T-273,molfr(:,2),'.g');L=plot(T-273,molfr(:,3),'.b');M=plot(T-273,molfr(:,4),'.c');N=plot(T-273,molfr(:,5),'.k');legend([H,K,L,M,N],'H2','CO','CO2','H2O','CH4');title('Molfraction of Pyrolysis products in Biomass Gasification');xlabel('Temperature in [C]');ylabel('massfraction');figure (2)hold ongrid onbox onplot(T-273,molfr(:,7),'.r');plot(T-273,molfr(:,8),'*b');plot(T-273,molfr(:,9),'+g');%% Pyrolysis of wood (pine Radiata)
function [dpyro] = Tinitialpyrolysis(z,x)Te=T(n);R=8.314;% Taken from "Modeling chemical and physical processes of wood and
% biomass pyrolysis ",Progress in Energy and Combustion
% Science 34 (2008) 47–90, Colomba Di Blasi.
% Reaction Schema for devolatilization of wood by proposed by
% Colomba Di Blasi
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% k_1
% |---------------> gas
% | ^
% | |k_4
% | k_2 |
% Biomass ----------------> tar
% | |
% | |k_5
% | k_3 |
% |------------->char
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Rxn 1,2,3,4,5 based on Colomba Di Blasi.
% k_j=k_0j exp(-E_Aj/((RT) ))
% k_0 is pre-exponential factor of the reaction
% E_A is activation energy of the reaction in (J(mol)^(-1) )
% R is the universal gas constant, in (J(mol)^(-1) K^(-1) )
% T is the operating temperature, in (K)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% k_oj(1/s) E_Aj(J/mole) deltH_rxnj(J/kg)
ko1=1.30e8; E_A1=140000; deltH_rxn1=64000; ko2=2.00e8; E_A2=133000; deltH_rxn2=64000; ko3=1.08e7; E_A3=121000; deltH_rxn3=64000; ko4=1.00e5; E_A4=93300; deltH_rxn5=-42000; ko5=1.48e6; E_A5=144000; deltH_rxn6=-42000;%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
k1=ko1*exp(-E_A1/(R*Te));k2=ko2*exp(-E_A2/(R*Te));k3=ko3*exp(-E_A3/(R*Te));k4=ko4*exp(-E_A4/(R*Te));k5=ko5*exp(-E_A5/(R*Te));%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%dM_biomass = -(k1+k2+k3)*x(1);dM_char = (k3*x(1))+(k5*x(3));dM_tar = (k2*x(1))-((k4*x(3))+(k5*x(3)));dM_gas = (k1*x(1))+(k4*x(3));dpyro = [dM_biomass; dM_gas; dM_tar; dM_char];endend
Best Answer