Error using / Matrix dimensions must agree

1 次查看（过去 30 天）
Riccardo Rinaldi
Riccardo Rinaldi 2018-5-18
回答： Walter Roberson 2018-5-18
clearvars
close all
clc
Ea = 98299.57711; %[-]
P = 101325; % [Pa]
R = 8.314; %[J/mol K]
d_t = 7*10^-3; %[m]
d_p = 2*10^-3; %[m]
rho_cat = 885.06692; %[kg/m3]
epsi = 0.4512; % [-]
% Fume Hood Temperature [K]
T_hood = 26 + 273.15;
% Normalized initial volumetric Flow [NmL/min]
Flux_norm0 = 209.930;
% Reactor Cross Section [m2]
Area = pi*d_t^2;
% Specific Heats matrix [J/(mol*K)]
cpM = [ 2.9108e+04 2.7617e+04 3.3298e+04 3.3363e+04 2.9370e+04 2.9105e+04 %C1
8.7730e+03 9.5600e+04 7.9933e+04 2.6790e+04 3.4540e+04 8.6149e+03 %C2
3.0851e+03 2.4660e+03 2.0869e+03 2.6105e+03 1.4280e+03 1.7016e+03 %C3
8.4553e+03 3.7600e+03 4.1602e+04 8.8960e+03 2.6400e+04 1.0347e+02 %C4
1.5382e+03 5.6760e+02 9.9196e+02 1.1690e+03 5.8800e+02 9.0979e+02 ]*1e-3; %C5
%CO %H2 %CH4 %H2O %CO2 %N2
%Viscosity [Pa*s]
muM = [ 1.1127e-06 1.797e-07 5.2546e-07 1.7096e-08 2.148e-06 6.5592e-07 %C1
0.5338 0.685 0.59006 1.1146 0.46 0.6081 %C2
94.7 -0.59 105.67 0 290 54.714 %C3
0 140 0 0 0 0 ]; %C4
%CO %H2 %CH4 %H2O %CO2 %N2
%Ideal gas enthalpy of formation [J/mol]
DHf = [-1.1053e+08 0 -7.481e+04 -2.4182e+05 -3.9351e+05 0 ]*1e-3;
%CO %H2 %CH4 %H2O %CO2 %N2
% Thermal conductivity [W/(m*K)]
kM = [ 5.9882e-04 2.653e-03 8.3983e-06 6.2041e-06 3.69 3.3143e-04 %C1
0.6863 0.7452 1.4268 1.3973 -0.3838 0.7722 %C2
57.13 12 0 0 964 16.323 %C3
501.92 0 0 0 1.86e+06 373.72 ]; %C4
%CO %H2 %CH4 %H2O %CO2 %N2
% Molecular Weight [kg/mol]
MW = [ 28.010 2.016 16.04 18.015 44.010 28.013 ]*1e-3;
%CO %H2 %CH4 %H2O %CO2 %N2
% Boiling Temperature [K]
Tb = 273.15+[ -192 -252.7 -161.4 100 -78.5 -195.8 ];
%CO %H2 %CH4 %H2O %CO2 %N2
S = 1.5.*Tb;
% Reference Temperature [K]
Tref = 298.15;
% Composition (Experimental Data) [-]
xi = [0.009786795 0.198499487 1.16216E-05 0.000487182 0 0.791214915 % T = 200 °C
0.00978733 0.198960411 3.15295E-05 0.000372155 0 0.790848575 % T = 220 °C
0.009596056 0.199206182 9.66385E-05 0.000398412 0 0.790702711 % T = 240 °C
0.009577825 0.198737454 0.00023838 0.000490008 1.17116E-05 0.79094462 % T = 260 °C
0.009197159 0.198126343 0.000611563 0.000825504 1.77622E-05 0.791221668 % T = 280 °C
0.008538797 0.196208971 0.001282088 0.001516948 3.31114E-05 0.792420085 % T = 300 °C
0.00713546 0.192796089 0.002666397 0.002838614 4.21538E-05 0.794521286 % T = 320 °C
0.004761065 0.186357537 0.005154589 0.005456679 3.2777E-05 0.798237352 % T = 340 °C
0.002403765 0.179820036 0.007523381 0.007810173 1.3767E-05 0.802428878 % T = 360 °C
0.001063646 0.176333059 0.008834865 0.009277417 2.06217E-05 0.804470392 % T = 380 °C
0.000623009 0.175212648 0.009234694 0.009890363 3.66325E-05 0.805002655 % T = 400 °C
0.000455837 0.174995752 0.009508677 0.009927986 4.47736E-05 0.805066975 % T = 420 °C
0.000422779 0.175275613 0.009556078 0.009932302 5.01422E-05 0.804763086 % T = 440 °C
0.000407729 0.174741413 0.009672711 0.009889148 7.15368E-05 0.805217461 % T = 460 °C
0.00050211 0.175579558 0.009462071 0.009827122 9.52404E-05 0.804533898 % T = 480 °C
0.000606778 0.176285873 0.009381173 0.009593589 0.000128888 0.804003699 ]; % T = 500 °C
%CO %H2 %CH4 %H2O %CO2 %N2
% Observed Reaction Rates [mol/(kg_CAT*s)]
r_MET = [1.23395E-05 3.34928E-05 0.000102675 0.000253193 0.000649339 0.001359222 0.002819341 0.005424882 0.007876526...
0.009226098 0.009637257 0.00992239 0.009975619 0.010091675 0.009880299 0.009802285];
r_WGS = [0 0 0 1.24394E-05 1.88594E-05 3.51035E-05 4.45717E-05 3.44958E-05 1.44132E-05 2.15349E-05 3.82294E-05...
4.67216E-05 5.23436E-05 7.46353E-05 9.945E-05 1.34674e-04];
% Solid Conductivity [W/(m*K)]
k_s = 35; %[W/m*K]
% Parameters
beta = 0.895;
gamma = 2/3;
b = 0;
% Vectors Initialization
T_vect = zeros(18,1);
rhomix_vect = zeros(18,1);
k_vect = zeros(18,6);
cp_vect = zeros(18,6);
Dh_vect = zeros(18,6);
mu_vect = zeros(18,6);
MW_xi_vect = zeros(18,6);
cp_xi_vect = zeros(18,6);
phi_vect = zeros(18,6,6); %?? PERCHé non è raccomandata la preallocazione?
Ss_vect = zeros(18,6,6);
Ai_vect = zeros(18,6,6);
DENm_vect = zeros(18,6,6);
DENk_vect = zeros(18,6,6);
DENmt_vect = zeros(18,6,1);
DENkt_vect = zeros(18,6,1);
NOMm = zeros(18,6,1);
NOMk = zeros(18,6,1);
mu_mix_vect = zeros(18,1,1);
k_mix_vect = zeros(18,1,1);
cp_mix_vect = zeros(18,1,1);
MW_mix_vect = zeros(18,1,1);
rho_mix_vect = zeros(18,1,1);
DHrMET_vect = zeros(18,1,1);
DHrWGS_vect = zeros(18,1,1);
SUM_DHr_r_obs_vect = zeros(18,1,1);
K_vect = zeros(18,1);
phi1_vect = zeros(18,1);
phi2_vect = zeros(18,1);
phi_vect = zeros(18,1);
lambdab0_lambdag_vect = zeros(18,1);
Re_vect = zeros(18,1);
Pr_vect = zeros(18,1);
Pe_rf_vect = zeros(18,1);
lambdaconv_lambdag_vect = zeros(18,1);
lambdaeff_lambdag_vect = zeros(18,1);
lambdaeff_vect = zeros(18,1);
for i = 1:16
T = 180 + i*20 + 273.15;
% Initial volumetric Flow [m3/s]
Flux0 = Flux_norm0/60 * T_hood/T * 10^-6;
% Velocity [m/s]
v = Flux0/Area;
% Viscosities [Pa*s] and Conductivities [W/(m*K)] of species
for j = 1:6
k = (kM(1,j)*(T^kM(2,j)))/(1 + kM(3,j)/T + kM(4,j)/(T^2));
cp = cpM(1,j) + cpM(2,j)*((cpM(3,j)/T)/sinh(cpM(3,j)/T))^2 + cpM(4,j)*((cpM(5,j)/T)/cosh(cpM(5,j)/T))^2;
Dh = DHf(j) + (2*cpM(2,j)*cpM(3,j))*((1/(exp(2*cpM(3,j)/T) - 1) - 1/(exp(2*cpM(3,j)/Tref) - 1)) + ...
(1/(exp(2*cpM(3,j)/T) + 1) - 1/(exp(2*cpM(3,j)/Tref) + 1)));
mu = (muM(1,j)*(T^muM(2,j)))/(1 + muM(3,j)/T + muM(4,j)/(T^2));
MW_xi = (MW(j)*xi(i,j));
cp_xi = (cp*xi(i,j));
T_vect(i,1) = T;
k_vect(i,j) = k;
cp_vect(i,j) = cp;
Dh_vect(i,j) = Dh;
mu_vect(i,j) = mu;
MW_xi_vect(i,j) = MW_xi;
cp_xi_vect(i,j) = cp_xi;
% Mixture Viscosity and Conductivity (for each Temperature)
for z = 1:6
phi = (1 + ((mu_vect(i,j)/mu_vect(i,z))^(1/2))*((MW(z)/MW(j))^(1/4)))^2/...
(8*(1 + MW(j)/MW(z)))^(1/2); % [-]
Ss = 0.735*sqrt(S(j)*S(z)); % [K]
Ai = 0.25*(1 + (mu_vect(i,j)/mu_vect(i,z)*(MW(z)/MW(j))^(3/4)*((T + S(j))/(T + S(z))))^(1/2))^2*...
((T + Ss)/(T + S(j))); % [%]
DENm = phi*xi(z); % [-]
DENk = Ai*xi(z); % [-]
phi_vect(i,j,z) = phi;
Ss_vect(i,j,z) = Ss;
Ai_vect(i,j,z) = Ai;
DENm_vect(i,j,z) = DENm;
DENk_vect(i,j,z) = DENk;
end
DENmt = sum(DENm,3); % [-]
DENkt = sum(DENk,3); % [-]
NOMm = mu_vect(i,j)*xi(j); % [Pa*s]
NOMk = k_vect(i,j)*xi(j); % [W/(m*K)]
DENmt_vect(i,j) = DENmt;
DENkt_vect(i,j) = DENkt;
NOMm(i,j) = NOMm;
NOMk(i,j) = NOMk;
end
mu_mix = sum(NOMm,2)/sum(DENmt,2);
k_mix = sum(NOMk,2)/sum(DENkt,2);
cp_mix = sum(cp_xi,2);
MW_mix = sum(MW_xi,2);
rho_mix = P*MW_mix/(R*T);
DHrMET = -Dh_vect(i,1) - 3*Dh_vect(i,2) + Dh_vect(i,3) + Dh_vect(i,4);
DHrWGS = -Dh_vect(i,1) - Dh_vect(i,4) + Dh_vect(i,5) + Dh_vect(i,2);
SUM_DHr_r_obs = DHrMET*r_MET(i) + DHrWGS*r_WGS(i);
% Estimation of lambdaeff
K = k_s/k_mix;
phi1 =(0.333 * (1 - 1/K)^2)/(log(K - 0.577 * (K - 1)) - 0.423 * (1 - 1/K)) - 2/(3 * K);
phi2 =( 0.072*(1-(1 - 1/K))^2)/(log(K - 0.925*(K - 1))-0.075*(1 -(1 - 1/K))) - 2/(3*K);
phi = phi2 + (phi1 - phi2)*((epsi - 0.260)/0.216);
lambdab0_lambdag = epsi + beta *(1 - epsi)/(phi + gamma*(k_mix/k_s));
Re = rho_mix * v * d_p/mu_mix;
Pr = mu_mix/(rho_mix * k_mix);
Pe_rf = 8.65 *(1 + 19.4 *(d_p/d_t)^2);
lambdaconv_lambdag = Re * Pr / Pe_rf;
lambdaeff_lambdag = lambdab0_lambdag + lambdaconv_lambdag;
lambdaeff = lambdaeff_lambdag*k_mix;
% Biot Number [-]
phi_w = exp(-0.004*(log(K))^4 + 0.0068*(log(K))^3 - 0.0306*(log(K))^2 ...
+ 0.1986*(log(K)) + 1.0698);
alfa_w = 0.041;
k_ew0 = epsi * k_mix +(1 - epsi)*k_s*(1/(phi_w*k_s/k_mix+1/3));
h_w = 2*k_ew0/d_p + alfa_w*(cp_mix*rho_mix*v);
Bi_w = h_w*d_p/lambdaeff;
%T gradient check %OBIETTIVO: FARSI RESTITUIRE QUESTO VALORE CHECK PER
%OGNI TEMPERATURA
Check = Ea/(R*T)*((SUM_DHr_r_obs*(1 - epsi)*(1 - b)*(d_t/2)^2)/(lambdaeff*T)*(1/8 + d_p/(Bi_w*d_t)));
mu_mix_vect(i) = mu_mix;
k_mix_vect(i) = k_mix;
cp_mix_vect(i) = cp_mix;
MW_mix_vect(i) = MW_mix;
rho_mix_vect(i) = rho_mix;
DHrMET_vect(i) = DHrMET;
DHrWGS_vect(i)= DHrWGS;
SUM_DHr_r_obs_vect(i) = SUM_DHr_r_obs;
K_vect(i) = K;
phi1_vect(i) = phi1;
phi2_vect(i) = phi2;
phi_vect(i) = phi;
lambdab0_lambdag_vect(i) = lambdab0_lambdag;
Re_vect(i) = Re;
Pr_vect(i) = Pr;
Pe_rf_vect(i) = Pe_rf;
lambdaconv_lambdag_vect(i) = lambdaconv_lambdag;
lambdaeff_lambdag_vect(i) = lambdaeff_lambdag;
lambdaeff_vect(i) = lambdaeff;
end
  1 个评论
Geoff Hayes
Geoff Hayes 2018-5-18
Riccardo - which line of code is throwing the error? Try putting a breakpoint at that line and re-run your code. When the debugger pauses at that line, take a look at the division - do the dimensions of the matrix make sense for what you are trying to do?

请先登录，再进行评论。

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

回答（1 个）

Walter Roberson
Walter Roberson 2018-5-18
When i becomes 2, mu_mix and k_mix are 2 x 1 vectors, so K = k_s/k_mix is matrix division giving a 1 x 2 result; the calculation is similar to k_s * pinv(k_mix) and is finding a least-squared solution. So your K is a vector, but the code below that point clearly expects it to be a scalar. If it is expected to be a scalar then you need to figure out which scalar it should be; if it is appropriate that it is a vector then you need to start converting * to .* and / to ./ and ^ to .^ .
I notice that further on you have
mu_mix_vect(i) = mu_mix
which indicates that you expect mu_mix to be a scalar rather than a vector.
So how are you constructing mu_mix? Well, you have
mu_mix = sum(NOMm,2)/sum(DENmt,2);
which is going to be a vector if NOMm has more than one row. So does it? Look above:
NOMm = mu_vect(i,j)*xi(j); % [Pa*s]
NOMk = k_vect(i,j)*xi(j); % [W/(m*K)]
DENmt_vect(i,j) = DENmt;
DENkt_vect(i,j) = DENkt;
NOMm(i,j) = NOMm;
NOMk(i,j) = NOMk;
The first two of those statements assign NOMm and NOMk as being scalars, but the last two lines there extend NOMm and NOMk into 2d arrays with i rows, and except for the case of i=1, j=1, two identical copies of the scalars that were calculated. I have to say that I have my doubts about whether that extension to 2D and the way that mu_mix and k_mix are calculated.

类别

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

标签

Community Treasure Hunt

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

Start Hunting!

Translated by