Error using / Matrix dimensions must agree

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

[Geoff Hayes]

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?

Answer 1

[Walter Roberson]

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.