Error using / Matrix dimensions must agree

1 view (last 30 days)
Riccardo Rinaldi
Riccardo Rinaldi on 18 May 2018
Answered: Walter Roberson on 18 May 2018
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 Comment
Geoff Hayes
Geoff Hayes on 18 May 2018
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?

Sign in to comment.

Sign in to answer this question.

Answers (1)

Walter Roberson
Walter Roberson on 18 May 2018
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.

Categories

Find more on Loops and Conditional Statements in Help Center and File Exchange

Tags

Community Treasure Hunt

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

Start Hunting!