Index exceeds matrix dimensions.

1 次查看（过去 30 天）
Riccardo Rinaldi
Riccardo Rinaldi 2017-11-2
编辑： Walter Roberson 2017-11-2
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

请先登录，再回答此问题。

回答（1 个）

Rik
Rik 2017-11-2
Let's start with the obvious problem: You didn't post a question, just a wall of unformatted code. Ask an actual question instead of letting us guess what you want us to do.
Another issue is the use of global variables. Global variables are a bad idea. In some very specific cases they are unavoidable, but under almost every circumstance you will be able to solve your problem without them. Use function arguments to pass variables from one workspace to another. If you really, absolutely, undoubtedly need global variables, use as few as possible, and with long names that extensively describe what they are and what they should do. This minimizes the risk of using the same global variable name in two scripts/functions.
After spending some time pasting your code into Matlab and adding back some returns, I have many m-lint warnings. Sometimes they can be ignored, but it is always a good place to start if you encounter problems. No warnings is no guarantee for good code, but at least you have a solid foundation.
This looks like a lot of manual labour, where the chance of a typo is huge. Finding out which index is wrong might be a needle in a haystack. In this case however, the solution is with my previous remark: there is an m-lint warning about S not being unused or unset. This means that it is an empty matrix when you try to access elements 1 through 5.
Next time, have a read here and here. It will greatly improve your chances of getting an answer.

类别

在 Help Center 和 File Exchange 中查找有关 Matrix Indexing 的更多信息

标签

Community Treasure Hunt

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

Start Hunting!

Translated by