3 views (last 30 days)

Show older comments

I understand that this code returns a vector length of 303 but the initail vector (y0) only contains 3. I am not sure how to fix this error.

Vspan = [0 10];

y0 = [0 570 (50*101.325)];

[V,y] = ode45(@(~,y) fun(del_h1,del_cp1,cp1),Vspan,y0);

function dydt = fun(~,y,del_h1,del_cp1,cp1)

X = y(1);

T = y(2);

P = y(3);

% inital parameters

T0 = 570; % [K] inlet temp

P0 = 50*101325; % [Pa] inlet pressure

Ft0 = 8000; % inlet total flow rate w mixture [kmol/hr]

%T = 570:2.3:800; % 800 is max temp allowed [K]

Ac = 3; % cross sectional area [m^3]

R = 8.314*1000; % [J/kmol.K]

%X= 0:.01:1;

%T = T';

%X = X';

% mol fractions

% let 1 = H2, 2 = CO, 3 = CH4, 4 = H2O

% 5 = CO2, 6 = C2H4, 7 = N2

y10 = 0.555; y20 = 0.215; y30 = 0.167;

y50 = 0.033; y60 = 0.019; y70 = 0.011; % inert gases

% molar masses [kg/kmol]

m1 = 2.016; m2 = 28.01; m3 = 16.04; m4 = 18.016; m5 = 44.01;

m6 = 30.07; m7 = 28.02;

% finding density of gas mixture (not sure if this is right way lol)

% should inert gases be included?

% found this at https://chemistry.stackexchange.com/questions/91120/finding-density-of-the-gas-mixture

mt = y10*m1 + y20*m2 +y30*m3 + y50*m5 + y60*m6 + y70*m7; % total mass [kg/kmol]

yt = y10 + y20 + y30 + y50 + y60 + y70; % total mole fraction

rho = (P.*mt)./(yt.*R.*T); % density of mixture

% flow rates [kmol/hr]

F10 = Ft0*y10; F20 = Ft0*y20; F30 = Ft0*y30;

% thetas

theta_b = F20/F10; theta_c = F30/F10;

% flow rates [kmol/hr]

delta = -2/3;

F1 = F10.*(1-X);

F2 = F10.*(theta_b - (X./3));

F3 = F10.*(theta_c + (X./3));

F4 = (F10.*X)/3;

Ft = Ft0 + delta.*F10.*X;

% masses

% finding concentrations

% note that: density = (molar flow)/(volumetric flow)

E = y10*delta;

vo = Ft0/Ac; % not the right volumetric flow rate

v = vo*(1+E.*X).*(T./T0).*(P0./P); % volumetric flow rate

C1 = F1./v; C2 = F2./v; C3 = F3./v; C4 = F4./v;

% let 1 = H2, 2 = CO, 3 = CH4, 4 = H2O, 5 = CO2, 6 = C2H4, 7 = N2

% parameters

a = [28.84 28.95 34.31 33.46 36.11 49.37 29]/10^3;

b = [.00765 .411 5.469 .688 4.233 13.92 .2199]/10^5;

c = [.3288 .3548 .3661 .7604 -2.887 -5.816 .5723]/10^8;

d = [-.8698 -2.22 -11 -3.593 7.464 7.28 -2.871]/10^12;

a = a';

b = b';

c = c';

d = d';

rho_cat = (1.14*(100^3))/1000; % [kg/m^3] catalyst density

D = ((2.81E-06).*T + 1.68E-04)/(100^2); % [m^2/s] effective diffusivity

area = 150*1000; % [m^2/kg.catalyst]

volume = (0.395*1000)/(100^3); % pore volume [m^3/kg]

UaTaT = 0; % this is U*a*(Ta-T) which is neglible

% arrhenius equation

A = [0.2687 9.706E-12 4.124E-09]/1000; % [in Pascal]

Ea = [70 -90 -68.5]*1000; % [kPa] may need to multiply by 1000 to get into Pa

K0 = A(1)*exp(-Ea(1)./(R.*T));

K1 = A(2)*exp(-Ea(2)./(R.*T));

K2 = A(3)*exp(-Ea(3)./(R.*T));

% partial pressures

P1 = P.*y10;

P2 = P.*y20;

r1 = -( (K0.*P1) ./ ( (1+sqrt(K1.*P1)+(K2.*P2)).^2) );

r1_prime = r1/rho_cat;

% effectiveness factor

tmod = 0.16.*sqrt(1./((2.81E-06).*T + 1.68E-04));

n = (3.*(tmod.*cot(tmod) - 1))./(tmod.^2); % effectiveness factor

% Erguen equation stuff

y = P./P0;

U = v./Ac;

phi = 0.4;

gc = 1; % for metric system

Dp = 0.32/100; % [m] particle diameter

G = rho_cat.*U;

u = 3.2*10^-4; % viscosity

Bo = ((G.*(1-phi))./(rho_cat.*gc.*Dp.*(phi.^3))).*(((150.*(1-phi)*u)./Dp) + 1.75.*G);

gamma = (2.*Bo)./(Ac.*rho_cat.*(1-phi).*P0);

dydt = [ -r1_prime.*n./F10; % dXdW

(1./rho_cat).*( (UaTaT + n.*r1_prime.*del_h1) ./ (F10.*(cp1+ del_cp1.*X)) ); % dTdW

(-gamma./(2.*y)).*(T./T0).*(1+E.*X) ];% dydW

%dydt = dydt';

end

James Tursa
on 6 May 2021

You have a mismatch in the number of arguments. E.g., this

[V,y] = ode45(@(~,y) fun(del_h1,del_cp1,cp1),Vspan,y0);

should be this

[V,y] = ode45(@(t,y) fun(t,y,del_h1,del_cp1,cp1),Vspan,y0);

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

Start Hunting!