Info

This question is closed. Reopen it to edit or answer.

Index exceeds matrix dimensions.

12 views (last 30 days)
Riccardo Rinaldi
Riccardo Rinaldi on 2 Nov 2017
Closed: Walter Roberson on 2 Nov 2017
clearvars
close all
clc
global eta eps rext rint delta rhoc cpi PinSweep
global DH0 Trif mui ki MW S R Ac Dp p ls Tw
% Gas Constant [J/(kmol*K)]
R = 8314.45;
% Catalyst Density [kgcat/m^3]
rhoc = 2400;
% Membrane Thickness [m]
delta = 5e-6;
% Reaction efficiency [-]
eta = 0.28;
%External Radius [m]
rext = 0.08;
% Internal Radius [m]
rint = 0.045;
% Cross Sectional Area [m^2]
Ac = pi*rint^2;
% Reactor Length [m]
L = 10;
% Void Fraction [-]
eps = 0.6;
% Stoichiometric Matrix [-]
nu = [-1 -1 1 1 0 0 0];
%CO %H2O %CO2 %H2 %Inert %H2perm %sweep
% Stoichiometric Matrix at Permeate side [-]
nup = [0 0 0 -1 0 -1 0];
%CO %H2O %CO2 %H2 %Inert %H2perm %sweep
% Specific Heats matrix [kJ/kmol*K]
cpi = 1e-3*[2.9108e+04 3.3363e+04 2.9370e+04 2.7617e+04 2.9105e+04 2.7617e+04 3.3363e+04 %C1
8.7730e+03 2.6790e+04 3.4540e+04 9.5600e+04 8.6149e+03 9.5600e+04 2.6790e+04 %C2
3.0851e+03 2.6105e+03 1.4280e+03 2.4660e+03 1.7016e+03 2.4660e+03 2.6105e+03 %C3
8.4553e+03 8.8960e+03 2.6400e+04 3.7600e+03 1.0347e+02 3.7600e+03 8.8960e+03 %C4
1.5382e+03 1.1690e+03 5.8800e+02 5.6760e+02 9.0979e+02 5.6760e+02 1.1690e+03]; %C5
%CO %H2O %CO2 %H2 %Inert %H2perm %sweep
% Viscosity [Pa*s]
mui = [ 1.1127e-06 1.7096e-08 2.148e-06 1.797e-07 6.5592e-07 1.797e-07 1.7096e-08 %C1
0.5338 1.1146 0.46 0.685 0.6081 0.685 1.1146 %C2
94.7 0 290 -0.59 54.714 -0.59 0 %C3
0 0 0 140 0 140 0 ]; %C4
%CO %H2O %CO2 %H2 %Inert %H2perm %sweep
% Thermal conductivity [W/m*K]
ki = [ 5.9882e-04 6.2041e-06 3.69 2.653e-03 3.3143e-04 2.653e-03 6.2041e-06 %C1
0.6863 1.3973 -0.3838 0.7452 0.7722 0.745 1.3973 %C2
57.13 0 964 12 16.323 12 0 %C3
501.92 0 1.86e+06 0 373.72 0 0 ]; %C4
%CO %H2O %CO2 %H2 %Inert %H2perm %sweep
% Ideal gas enthalpy of formation [kJ/kmol] Trif =298.15 K
DH0 = 1e-3*[-1.1053e+08 -2.41814e+08 -3.9351e+08 0 0 0 -2.41814e+08 ];
%CO %H2O %CO2 %H2 %Inert %H2perm %sweep
% Reference Temperature [K]
Trif = 298.15;
% Wall Temperature [K]
Tw = 628.15;
% Initial Sweep Temperature [K]
TinSweep = 628.15;
% Initial Temperature [K]
Tin = 628;
% Initial Sweep Pressure [kPa]
PinSweep = 101;
% Initial Pressure [kPa]
Pin = 2376.53;
% Initial Sweep Flow Rate [kmol/h]
FinSweep = 8;
% Initial Flow Rate [kmol/h]
Ftot0 = 14;
% Inlet Compositions at Permeation Side [-]
F0Perm = 0;
% Inlet Compositions of reacting Mixture [-]
x0 = [0.0521 0.4453 0.0272 0.4649 0.0105];
%CO %H2O %CO2 %H2 %Inert
% MW [kg/kmol]
MW = [ 28.010 18.015 44.010 2.016 28.013 2.016 18.015 ];
%CO %H2O %CO2 %H2 %Inert %H2perm %sweep
% Normal Boiling Temperature and corresponding Parameters [K]
Tb = [81.15 373.15 194.65 79 77.35 79 373.15 ];
%CO %H2O %CO2 %H2 %Inert %H2perm %sweep
S = 1.5.*Tb;
% Global Heat Exchange coefficient for permeate side [kJ/(m^2*h*K)]
Um = 8.64;
% Particle Diameter [m]
Dp = 0.0053;
% Thermal Conductivity of packing Material [kJ/(m*h*K)]
ls = 1.256;
% Solid Surface Emissivity [-]
p = 0.8;
% Inlet Conditions
In = [Ftot0*x0, F0Perm, FinSweep, Pin, Tin, TinSweep - 100];
% Sistem Solution
zspan = [0 L]; % [m]
options = odeset('Reltol',1e-11,'Abstol',1e-12);
[Z,X] = ode15s(@Tronky,zspan,In,options);
function f = Tronky(z,x)
% Variable Definition
global eta eps rext rint delta rhoc cpi PinSweep
global DH0 Trif mui ki MW S R Ac Dp p ls Tw
Fi = x(1:5);
FiPS = x(6:7);
P = x(8);
T = x(9);
Tperm = x(10);
% Other Variables
Ftot = sum(Fi); % [kmol/h]
xi = Fi./Ftot; % [-]
Pi = P*xi; % [kPa]
FtotPS = sum(FiPS); % [kmol/h]
xiPS = FiPS./FtotPS; % [-]
PiPS = PinSweep*xiPS; % [kPa]
Tmem = (T + Tperm)/2; % [K]
% Reaction, Adsorption, Equilibrium Constants
k = exp(-29364/(1.987*T) + 40.32/1.987); % [kmol/(kgcat*h)]
Kco = exp(3064/(1.987*T) - 6.74/1.987)*1/(101325); % [1/Pa]
Kh2o = exp(-6216/(1.987*T) + 12.77/1.987); % [1/Pa]
Kco2 = exp(12542/(1.987*T) - 18.45/1.987); % [1/Pa]
Keq = exp(4400/T - 4.063); % [-]
% Reaction Rate [kmol/m^3*h]
rco = k*Kco*Kh2o*(Pi(1)*Pi(2) - Pi(3)*Pi(4)/Keq)*(1 + Kco*Pi(1) + Kh2o*Pi(2) + Kco2*Pi(4))^-2*rhoc;
Rco = rco*eta*(1 - eps)*pi*(rext^2 - rint^2);
% Membrane Permeability [kmol/(m*h*Pa^0.5)]
Bh = 2.95*1e-4*3600/1e3*exp(-5833.5/Tmem);
% Hydrogen Flux [kmol/(h*m^2)]
Jh2perm = Bh/delta*(Pi(4)^0.5 - PiPS(2)^0.5);
Rh2 = Jh2perm*2*pi*rint;
% Specific Heats, Enthalpies and molecular Weights of reacting Mixture
for i = 1:5
cp(i) = cpi(1,i) + cpi(2,i)*((cpi(3,i)/T)/sinh(cpi(3,i)/T))^2 +...
cpi(4,i)*((cpi(5,i)/T)/cosh(cpi(5,i)/T))^2; % [kJ/(kmol*K)]
DHi(i) = DH0(i) + (2*cpi(2,i)*cpi(3,i))*((1/(exp(2*cpi(3,i)/T) - 1) - 1/(exp(2*cpi(3,i)/Trif) - 1)) + ...
(1/(exp(2*cpi(3,i)/T) + 1) - 1/(exp(2*cpi(3,i)/Trif) + 1))); % [kJ/kmol]
cpmixi(i) = cp(i)*xi(i); % [kJ/(mol*K)]
MWi(i) = MW(i)*xi(i);
end
% Mixture Molecular Weight [kg/kmol]
MWm = sum(MWi);
% Mixture Density [kg/m^3]
rhom = P/(R*T)*MWm;
% Mixture Specific Heat [kJ/(mol*K)]
cpm = sum(cpmixi(i));
% Reaction Heat [kJ/mol]
DHr = -DHi(1) - DHi(2) + DHi(3) + DHi(4);
% Viscosities and Conductivities
for i = 1:5
mu = mui(1,i)*(T^(mui(2,i)))/(1 + mui(3,i)/T + mui(4,i)/(T^2)); % [Pa*s]
k = 3600/1e3*(ki(1,i)*(T^(ki(2,i)))/(1 + ki(3,i)/T + ki(4,i)/(T^2))); % [kJ/(m*h*K)]
end
% Mixing Rules for Viscosities and Conductivities
for i = 1:5
for j= 1:5
phi(i,j) = (1 + ((mu(i)/mu(j))^(1/2))*((MW(j)/MW(i))^(1/4)))^2/...
(8*(1 + MW(i)/MW(j)))^(1/2); % [-]
Ss(i,j) = 0.735*sqrt(S(i)*S(j)); % [K]
Ai(i,j) = 0.25*(1 + (mu(i)/mu(j)*(MW(j)/MW(i))^(3/4)*((T + S(i))/(T + S(j))))^(1/2))^2*...
((T + Ss(i,j))/(T + S(i))); % [-]
DENm(i) = phi(i,j)*xi(j); % [-]
DENk(i) = Ai(i,j)*xi(j); % [-]
end
DENmt(i) = sum(DENm); % [-]
DENkt(i) = sum(DENk); %
NOMm(i) = mu(i)*xi(i); % [Pa*s]
NOMk(i) = k(i)*xi(i); % [kJ/(m*h*K)]
end
mum = sum(NOMm)/sum(DENmt);
km = sum(NOMk)/sum(DENkt);
% Mixture Velocity [m/s]
v = Ftot*MWm*1/(rhom*Ac*3600);
% Adimensional Numbers [-]
Re = rhom*v*Dp/mum;
Pr = mum*cpm/km;
% Parameters for global Heat Exchange Coefficient
alpharu = (0.8171*(T/100)^3)/(1 + eps/(2*(1 - eps))*((1 - p)/p));
alphars = 0.8171*(p/(2 - p))*(T/100)^3;
ler0 = eps*(km + 0.95*alpharu*Dp) + ((0.95*(1 - eps))/(2/(3*ls) + 1/(10*km + alphars*Dp))); % [kJ/(m*K*h)]
ler = ler0 + 0.111*km*(Re*Pr^(1/3))/(1 + 46*(Dp/(2*rext))^2); % [kJ/(m*K*h)]
alphaw = 8.694*(ler0/(2*rext)^(4/3)) + 0.512*km*2*rext*Re*Pr^(1/3)/Dp;
Bi = alphaw*rext/ler; % [-]
U = (1/alphaw + rext/(3*ler)*(Bi + 3)/(Bi + 4))^(-1); % [kJ/(m^2*K*h)]
% Specific Heats and Enthalpies of Permeation Side
for i = 6:7
cp(i) = cpi(1,i) + cpi(2,i)*((cpi(3,i)/T)/sinh(cpi(3,i)/T))^2 +...
cpi(4,i)*((cpi(5,i)/T)/cosh(cpi(5,i)/T))^2; % [kJ/(kmol*K)]
DHi(i) = DH0(i) + (2*cpi(2,i)*cpi(3,i))*((1/(exp(2*cpi(3,i)/T) - 1) - 1/(exp(2*cpi(3,i)/Trif) - 1)) + ...
(1/(exp(2*cpi(3,i)/T) + 1) - 1/(exp(2*cpi(3,i)/Trif) + 1))); % [kJ/kmol]
DHTpermi(i) = DH0(i) + (2*cpi(2,i)*cpi(3,i))*((1/(exp(2*cpi(3,i)/Tperm) - 1) - 1/(exp(2*cpi(3,i)/Trif) - 1)) + ...
(1/(exp(2*cpi(3,i)/Tperm) + 1) - 1/(exp(2*cpi(3,i)/Trif) + 1))); % [kJ/kmol]
end
% Ergun Expression
fv = 150 + 1.75*Re/(1 - eps);
Pressure = fv*v*mum/Dp^2*((1 - eps)^2/eps^3);
% Balance Equations
f(1) = -Rco;
f(2) = -Rco;
f(3) = Rco;
f(4) = Rco - Rh2;
f(5) = 0;
f(6) = -Rh2;
f(7) = 0;
f(8) = Pressure;
f(9) = 1/(Ftot*cpm)*(Rco*DHr + U*2*pi*rext*(Tw - T) - Um*2*pi*rint*(T - Tperm));
f(10) = 1/(FiPS(1)*cp(6) + FiSweep*cp(7))*(Um*2*pi*rint*(T - Tperm) + Rh2*(DHTpermi(6) - DHi(6)));
f = f(:);
end
  2 Comments
Geoff Hayes
Geoff Hayes on 2 Nov 2017
Edited: Geoff Hayes on 2 Nov 2017
Riccardo - please copy and paste the full error message so that we can see which line is generating this error message. The error message is clear - the code is trying to access (get) an element in a matrix using an index that is greater than the dimensions of the matrix.
Walter Roberson
Walter Roberson on 2 Nov 2017
Duplicates https://www.mathworks.com/matlabcentral/answers/364780-index-exceeds-matrix-dimensions

This question is closed.

Answers (0)

This question is closed.

Tags

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!