Simulation of a PFR catalyst reactor with impact of heat and mass transfer

I've tried to solve the mass balance and energy balance in bulk and on the surface but i encountered the error below. please if someone can help me to do it.Thanks alot in advance

((Error using *

Incorrect dimensions for matrix multiplication. Check that the number of columns in the first matrix

matches the number of rows in the second matrix. To perform elementwise multiplication, use '.*'.

Error in daeic12 (line 63)

F = UM' * f;

Error in ode15s (line 304)

[y,yp,f0,dfdy,nFE,nPD,Jfac] = daeic12(odeFcn,odeArgs,t,ICtype,Mt,y,yp0,f0,...

Error in creapp1c (line 56)

[V, sol] = ode15s(@BMiBE,Vspan,Y_in,opt, nu, Dreact, Dp, P, Rg);

Related documentation ))

the error lne has been mension in the inseted matlab code .I know the meaning of the error but I couldnt figure out in which line it occur.

clc, clear, close all

%Process data

P = 5; %[bar]

Tb0 = 850 + 273.15; %[K]

Tint0=Tb0;

Rg=8.31447;

rhocat = 1900;

rhoskel = 2900;

Dp = 0.02;

Dreact = 0.125;

eps = 1 - (rhocat/rhoskel);

aspec=6*(1-eps)/Dp;

VfracS=0.6;

N0 = [50 150 0 2.63 0];

Cguess = N0/sum(N0)*P*10^5/(Rg*Tb0)/2;

Vfin = 0.0107; %[m^3]

Vspan = linspace(0,0.0107,200);

% R1 R2 R3

nu = [-1 , 0 , -1 %CH4

-1 , -1 , -2 %H2O

1 , -1 , 0 %CO

3 , 1 , 4 %H2

0 , 1 , 1]; %CO2

[hm, h, Vdot]=transfcoeff (Tb0, N0', Dreact, Dp, Rg, P);

Cb = N0/Vdot;

guess0 = [Cguess];

disp(guess0)

options = optimset('Display','off','TolFun',1e-15);

[guess, ~] = fsolve(@BMintf,guess0,options,hm, h, Cb ,Tb0,P ,Rg, aspec, nu,N0);

Cint0 = guess(1:end);

disp(Cint0)

Y_in=[N0, Cint0,Tb0,Tb0];

% Mass matrix 10x10: MBs at bulk, MBs at interface, EBs=0

m= [1 1 1 1 1 0 0 0 0 0];

M = diag(m);

opt= odeset('Mass',M,'Abstol',1e-3);

[V, sol] = ode15s(@BMiBE,Vspan,Y_in,opt, nu, Dreact, Dp, P, Rg); error line (line 56)

Nb = sol(:, 1:5);

Cint = sol(:, 6:10);

disp(Nb)

C_bulk = [];

for i = 1:length(V)

[y, cb] = BMiBE(V(i),sol(i,:)', nu, Dreact, Dp, P, Rg);

C_bulk = [C_bulk cb];

end

disp(C_bulk)

%PLOTS

figure(1)

plot(V,Cint(:,1))

hold on

plot(V,C_bulk(1,:))

legend('CH4 int','CH4 bulk');xlabel('Volume [m^3]');ylabel('C [mol/m^3]')

hold off

figure(2)

plot(V,Cint(:,2))

hold on

plot(V,C_bulk(2,:))

legend('H2O int','H2O bulk');xlabel('Volume [m^3]');ylabel('C [mol/m^3]')

hold off

figure(3)

plot(V,Cint(:,3))

hold on

plot(V,C_bulk(3,:))

legend('CO int','CO bulk');xlabel('Volume [m^3]');ylabel('C [mol/m^3]')

hold off

figure(4)

plot(V,Cint(:,4))

hold on

plot(V,C_bulk(4,:))

legend('H2 int','H2 bulk');xlabel('Volume [m^3]');ylabel('C [mol/m^3]')

hold off

figure(5)

plot(V,Cint(:,5))

hold on

plot(V,C_bulk(5,:))

legend('CO2 int','CO2 bulk');xlabel('Volume [m^3]');ylabel('C [mol/m^3]')

hold off

figure (6)

XCH4 = 1 - (Nb(:,1)./N0(1));

plot(V, XCH4)

xlabel('Volume [m^3]');ylabel('XCH4');

% MB for 5 species + EB

function [Yprimo, Cb] = BMiBE(V,y, nu, D_react, Dp, P, Rg)

Tb = y(12);

Tint = y(11);

Nb = y(1:5); %molar flowrates

Cint= y(6:10); %concentration

eps = 0.345;

aspec = 6*(1-eps)/Dp; %specific area [m^2]

%specific heats

acp=[5.529e-3 1.558e-3 0.716e-3 0.370e-3 1.593e-3];

bcp =[3.292 3.409 3.269 3.266 4.925];

cp = (acp.*Tb+bcp).*Rg;

Ncp = (cp(1)*Nb(1)+cp(2)*Nb(2)+cp(3)*Nb(3)+cp(4)*Nb(4)+cp(5)*Nb(5));

[rsurf, Rsurf] = r(y, P, Rg, aspec, nu);

[hm, h, Vdot] = transfcoeff (Tb,Nb, D_react, Dp, Rg, P);

%Entalpy of reaction for each reaction [kJ/mol]

H=[-4.47e13.*Tint^(-4.459)+226.9;

-271.4.*Tint^(-0.2977) ;

99.52.*Tint^(0.0937) ].*10^3;

%switch caso

Q = -h.*aspec.*(Tint-Tb); %isoT

% Q=0; %adb

% Q=600e6; %uniform heat w/m3

Cb = Nb./Vdot;

y1 = hm.*(Cint-Cb).*aspec;

y2 = rsurf - hm.*(Cint-Cb);

t1 = (Q + h.*aspec.*(Tint-Tb))./Ncp;

t2 = H'*Rsurf+h.*(Tint-Tb);

Yprimo = [y1; y2;t1;t2];

end

function [hm, h, Vdot]=transfcoeff (Tb,Nb, D_react, Dp, Rg, P)

Surf = pi*D_react.^2/4;

MW = [16.04 18.015 28.01 2.016 44.01]*10^(-3); %kg/mol

MWmix = (MW*Nb)/sum(Nb');

rhomix = ((P*10^5)/(Rg*Tb))*MWmix;

Vdot = sum(Nb)*Rg*Tb/(P*10^5);

v_0 = Vdot./Surf;

%diffusion parameters

a = 10^(-10).*[1.932, 1.718, 2.454, 11.275, 1.248];

b = [1.7888, 1.7884, 1.7370, 1.6936, 1.8033];

Dimix=a.*Tb.^b; %m^2/s

% for i=1:5

% Dimix(1,i) = a(i).*Tb^(b(i));

% end

mu = 3.144e-08*Tb + 2.474e-06;

cpmix = 1.386*Tb + 1648; %[J/kg K]

k = 0.0002781*Tb - 0.107;

Re = (rhomix*v_0*Dp)/mu;

Sc = mu./(rhomix.*Dimix);

Pr=cpmix*mu/k;

Jd = 2.19./Re.^(2/3) + 0.78./Re.^0.381;

% Sh = Jd.*Sc.^(1/3).*Re;

hm=(Jd*v_0)./(Sc.^(2/3))';

h=Jd./(Pr.^(2/3))*cpmix.*rhomix.*v_0;

end

function f = BMintf(Guess,hm, h, Cb,Tb0, P, Rg, aspec, nu,N0) % =========================================

Tint0=Tb0;

Cint = Guess(1:end);

H=[-4.47e13*Tint0^(-4.459)+226.9, -271.4.*Tint0^(-0.2977), 99.52.*Tint0^(0.0937)]*10^3;

y0=[N0 Cint Tint0 Tb0];

[rsurf, Rsurf] = r(y0, P, Rg, aspec, nu);

f1 = hm .* (Cb - Cint)' + rsurf;

f2 = H*Rsurf + h*(Tint0-Tb0);

f = [f1' f2];

end

function [rsurf, Rsurf] = r(y, P, Rg, aspec, nu)

rhocat = 1900;

p = y(6:10)./sum(y(6:10)).*P;

A = [4.225e15 , 1.955e6 , 1.020e15];

Aeq = [ 8.23e-5 , 6.12e-9 , 6.65e-4 , 1.77e5]; %A, Aeq

E = [240.1 , 67.13 , 243.9]*10^3;

deltaH = [-70.65 , -82.90 , -38.28 , 88.68 ]*10^3; %E, deltaH

Kp1 = exp(30.481-27.187e3/y(11));

Kp2 = exp(-3.924+4.291e3/y(11));

Kp3 = exp(26.891-23.258e3/y(11));

Kp=[Kp1,Kp2,Kp3];

k = A.*exp(-E/(Rg*y(11)));

Keq = Aeq.*exp(-deltaH/(Rg*y(11)));

DEN = 1+Keq(1)*p(3)+Keq(2)*p(4)+Keq(3)*p(1)+Keq(4)*p(2)/p(4);

Ri = [k(1)*(p(1)*p(2)-(p(4))^3*p(3)/Kp(1))./((p(4)^2.5).*DEN^2)*1000/3600.*rhocat

k(2)*(p(3)*p(2)-p(4)*p(5)/Kp(2))./(p(4).*DEN^2)*1000/3600.*rhocat

k(3)*(p(1)*(p(2))^2-(p(4))^4*p(5)/Kp(3))./((p(4)^3.5).*DEN^2)*1000/3600.*rhocat]; %[mol/s/m^3]

Rsurf = Ri./aspec;

rsurf = (nu*Rsurf);

end