Paxton asked . 2020-12-08

Index exceeds number of array elements (1)

%% 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

matlab , array index , matrix