Info

この質問は閉じられています。 編集または回答するには再度開いてください。

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
  2 件のコメント
Geoff Hayes
Geoff Hayes 2017 年 11 月 2 日
編集済み: Geoff Hayes 2017 年 11 月 2 日
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 2017 年 11 月 2 日
Duplicates https://www.mathworks.com/matlabcentral/answers/364780-index-exceeds-matrix-dimensions

この質問は閉じられています。

回答 (0 件)

この質問は閉じられています。

タグ

Community Treasure Hunt

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

Start Hunting!

Translated by