I'm trying to run a code to plot conversion vs catalyst mass in a packed bed reactor, but something is wrong.

[Walter Roberson]

How does the output differ from what you expect?

Answer 1

[Torsten]

You tell the integrator that you have four additional inputs to your function

[W,F]=ode45(@(t,F)ODEfun(t,F,T,P_total,F0_total,R),W,F0); 

but you only have two:

function dFdW = ODEfun(W,F,T,P_total) 

[Luna Naddi]

[W,F]=ode45(@(W,F)ODEfun(W,F,T,P_total));

I changed it to this and still got an error

[Torsten]

I don't get an error (see above). If the plots are not as expected, I think we will not be able to help.

[Luna Naddi]

this plot is as expected but why does it not display it for me?

[Torsten]

Since you edited away your question, I assume the issue has been solved :-)

Will such a harmless question be thus strictly punished that you feel you have to do that ?

[Luna Naddi]

I managed to get one of the plots required of me. I just now have to figure out how to plot conversion vs catalyst mass.

[Walter Roberson]

Since you deleted your code, we are unable to assist you in that difficulty.

[Luna Naddi]

clc 
clear all
% Components Directory for number indexed 
% 1 : B : H2
% 2 : D : CO2
% 3 : A : CO 
% 4 : C : CH3OH
% 5 : E : H2O
% 6 : F : N2
% R is gas constant  
% T is temperature 
x_0 = 0; 
x_end = 6.75; %g 
%Variables
T = 473; %degree Kelvin  
P_total = 50; %bar
y1_0 = 0.698;     
y2_0 = 0.178;
y3_0 = 0;
y4_0 = 0;
y5_0 = 0;
y6_0 = 0.124;
Ftotal = 4.6064*(10^-4); %mol/s
R = 8.314; %J/mol*K
F1_0 = y1_0 * Ftotal;      
F2_0 = y2_0 * Ftotal;        
F3_0 = 0;       
F4_0 = 0; 
F5_0 = 0;
F6_0 = y6_0 * Ftotal;
% Partial pressure 
P1 = P_total * y1_0; 
P2 = P_total * y2_0;
P3 = P_total * y3_0;
P4 = P_total * y4_0;
P5 = P_total * y5_0;
P6 = P_total * y6_0;
x = [x_0 x_end]; 
y0 = [y1_0 y2_0 y3_0 y4_0 y5_0 y6_0];
% Define Kinetic Constants: % mol s-1 bar-1 kg-1
R = 8.314; 
k1 = 4.89*10^7*exp(-113000/(R*T)); 
k2 = 9.64*10^11*exp(-152900/(R*T));
k3 = 1.09*10^5*exp(-87500/(R*T));
% Define Adsorption constants: 
bCO  = 2.16*10^(-5)*exp(46800/(R*T)); 
bCO2 = 7.05*10^(-7)*exp(61700/(R*T));
bfrac = 6.37*10^(-9)*exp(84000/(R*T)); 
% Define Equilibrium Constants: 
K1 = 10^((5139/T)-12.621);
K2 = 10^((-2073/T)+2.029); 
K3 = 10^((3066/T)-10.592);
% Base Rate Expressions
r1 = k1*bCO*((P3*P1^(3/2)-P4/(P1^(1/2)*K1))/((1+bCO*P3+bCO2*P2)*(P1^(1/2)+bfrac*P5)));
r2 = k2*bCO2*((P2*P1-P3*P5/K2)/((1+bCO*P3+bCO2*P2)*(P1^(1/2)+bfrac*P5)));
r3 = k3*bCO2*((P2*P1^(3/2)-P4*P5/(P1^(3/2)*K3))/((1+bCO*P3+bCO2*P2)*(P1^(1/2)+bfrac*P5)));
% Calling the ODE solver
[x,y]=ode45(@(t,Y)ODEfun(t,Y,r1,r2,r3,y1_0,y2_0,y3_0,y4_0,y5_0,y6_0),x,y0);
%Catalyst mass vs conversion
figure(1)
hold on
plot(x,y(:,1), x,y(:,2), x,y(:,3), x,y(:,4), x,y(:,5), x,y(:,6), 'LineWidth',2) 
leg=legend('H2','CO2','CO','CH3OH','H2O','N2');
leg.ItemTokenSize=[0,10];
xlabel('catalyst mass, kg');
ylabel('Molar fraction, -');
title('Molar fraction vs. Catalyst Mass');
grid on 
ylim([0, max(max(y))+0.05])
function dFdW = ODEfun(x,Y,r1,r2,r3,y1_0,y2_0,y3_0,y4_0,y5_0,y6_0);
F1 = Y(1);
F2 = Y(2);
F3 = Y(3);
F4 = Y(4);
F5 = Y(5);
F6 = Y(6);
y1 = F1 / (F1 + F2 + F3 + F4 + F5 + F6);
y2 = F2 / (F1 + F2 + F3 + F4 + F5 + F6);
y3 = F3 / (F1 + F2 + F3 + F4 + F5 + F6);
y4 = F4 / (F1 + F2 + F3 + F4 + F5 + F6);
y5 = F5 / (F1 + F2 + F3 + F4 + F5 + F6);
y6 = F6 / (F1 + F2 + F3 + F4 + F5 + F6);
rA_net = y1*(-2*r1 -r2 -3*r3);
rB_net = y2*(-1*r1 - r3);
rC_net = y3*(-1*r1 + r2);
rD_net = y4*(-1*r1 + -1*r3);
rE_net = y5*(r2 + r3);
rF_net = 0;
% These the values that input in the equations above
if y1 <0 | y2 <0 | y3 <0 | y4 <0
    Z1 = 0;
    Z2 = 0;
    Z3 = 0;
    Z4 = 0;
    Z5 = 0;
    Z6 = 0;
end
% The code that equate the value of Z to the ODE's.
dy1dx = rA_net; 
dy2dx = rB_net; 
dy3dx = rC_net; 
dy4dx = rD_net; 
dy5dx = rE_net;
dy6dx = rF_net;
dFdW = [dy1dx; dy2dx; dy3dx; dy4dx; dy5dx; dy6dx];
end 

[Torsten]

Why do you use the initial mole fractions to calculate the reaction rates ? They have to be changed according to the actual composition.

And what is your question about the plot ? Technically, I had no problems with your above code and the plot.

[Luna Naddi]

I need to obtain a plot of conversion vs catalyst mass and I need to alter the code above to do so.

[Luna Naddi]

this is the plot given to us, I have to obtain the same one with my code

[Torsten]

I'm not an expert in chemical reaction engineering and can't help you in this respect, but I think your code is far off to answer the question behind the graphics.