Index exceeds number of array elements (1)
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_char R=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 + H2 1,-1, 1, -1, 0, 0, 0;... % CO + H2O --> CO2+ H2 -2, 0, 0, 0, 1, -1, 0;... % C + 2H2 --> CH4 0, 2,-1, 0, 0, -1, 0;... % C + CO2 --> 2CO 3, 1, 0, -1, -1, 0, 0;... % CH4+ H2O --> CO + 3H2 0, 2, 0, 0, 0, -2,-1;... % 2C + O2 --> 2CO 2, 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_solver close all clear all clc global n M_biomass T M_char M_gas n_char T = 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 Carbon Lg_O=x(7); %La grangian for Oxygen Lg_H=x(8); %La grangian for Hydrogen fn(n,:)=fval; mol(n,:) =[x(1:5) n_char]; molfr(n,:) =[yh2 yco yco2 yh2o ych4]; Matrix=x; end Eflag fn mol figure(1) hold on grid on box on H=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 on grid on box on plot(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]; end end
Prashant Kumar answered .
2025-11-20
function fun=nonstqsol(x) global n M_biomass T M_gas M_char n_char R=8.314; a=0.2; b=0.13; x=1.44;M
This means that, no matter what your input for x is, it immediately gets replaced by x=1.44. A little later in your function, you try to define nCO:
nCO =x(2); %Mole of carbon monoxide
However, this results in an error because x only contains a single value. Hence, your index (2) exceeds the number of elements (1) in x.