3 replies to this topic

[sharna]

I am trying to model a Packed bed reactor for ethylbenzene dehydrogenation reaction for styrene production on Matlab but the code does not work, not sure why. I am posting what I did on Matlab. 

 

function dydt= @DehyR1 (w,x)

 

 

Eb=x(1); %ethylbenzene

H2=x(2); %hydrogen

Bz=x(3);%benzene

tol=x(4); %toluene

T=x(5);%temperature

P=x(6);%pressure

 

%initial condition: kmol/hr

eb0=478; 

st0=0.4842;

h20=0;

bz0=19.66;

tol0=0.4587;

c2h40=0;

ch40=0;

steam0=5258;

T0=620+273; %K

p0=1; %bar

 

total0=5756.6029;

 

%catalyst information

e_B=0.4312;

rho_B=1422; %catalyst bulk weight

 

%rate constants, kmol/kgcat.hr) (kJ/kmol)

kc1=4594000000*exp(-21095/T);

kc2=(1.60*10^15)*exp(-35637/T);

kc3=(1.246*10^26)*exp(-57104/T);

keb=(1.014*10^-5)*exp(12295/T);

kst=(2.678*10^-5)*exp(12576/T);

kh2=(4.519*10^-7)*exp(14187/T);

keq=exp((S1/8.314)-(H1/(8.314*T)));

kt1=(2.2215*10^16)*exp(-32744/T);

kt2=(2.4217*10^20)*exp(-42433/T);

kt3=(3.8224*10^17)*exp(-37655/T);

 

X=0.44; %conversion

Tr=25+273; %reference temp

 

%selectivities

s1=0.96; %styrene

s2=0.005; %benzene

s3=0.035; %toluene

 

%mole out

feb=eb0(1-X);

fst=st0+eb0*X*s1;

fh2=h20+eb0*X*s1-eb0*X*s3;

fbz=bz0+eb0*X*s2;

ftol=tlo0+eb0*X*s3;

fch4=ch40+eb0*X*s3;

fc2h4=c2h40+eb0*X*s2;

fsteam=steam0;

ft=feb+fest+ftol+fbz+fh2+fch4+fc2h4+fsteam;

 

%partial Pressure

 

PPeb=(feb*p0)/ft;

PPst=(fst*p0)/ft;

PPh2=(fh2*p0)/ft;

 

N=(1+keb*PPeb+kst*PPst+kh2*PPh2)^2;

 

rc1=((kc1*keb)*(PPeb-(PPst*PPh2/keq))/N;

rc2=(kc2*keb*PPeb)/N;

rc3=(kc3*keb*PPeb*kh2*PPh2)/N;

rt1=kt1*(PPeb-(PPst*PPh2/Keq));

rt2=kt2*PPeb;

rt3=kt3*PPeb;

reb=rc1+rc2+rc3+(e_B/rho_B )*(rt1+rt2+rt3;

 

%Specific Heat Capacity

delA1=41.99;

delB1=-0.082026;

delC1=0.00006499;

delD1=-2.311E-08;

S1ref=113.42;

H1ref=117690;

 

delA2=1.2986E+01;

delB2=-7.6700E-02;

delC2=9.5920E-05;

delD2=-4.1250E-08;

H2ref=1.0551E+05;

 

delA3=10.86;

delB3=-0.151844;

delC3=0.0002304;

delD3=-9.9955E-08;

H3ref=-54680;

 

CPn=343503.206;

C=(65.838*T)-(0.3105*T^2)+(39.13*(10^-5)*T^3)-(16.43*(10^-8)*T^4;

%Entropy and Heat of Reaction

S1=S1ref+(delA1*(ln(T/Tr))+(delB1*(T-Tr))+(delC1*(T^2-Tr^2)/2)+(delD1*(T^3-Tr^3)/3);

 

H1=H1ref+(*(T-Tr))+(delB1*(T^2-Tr^2)/2)+(delC1*(T^3-T^3)/3)+(delD1*(T^4-Tr^4)/4);

H2=H2ref+(delA2*(T-Tr))+(delB2*(T^2-Tr^2)/2)+(delC2*(T^3-T^3)/3)+(delD2*(T^4-Tr^4)/4);

H3=H3ref+(delA3*(T-Tr))+(delB3*(T^2-Tr^2)/2)+(delC3*(T^3-T^3)/3)+(delD3*(T^4-Tr^4)/4);

 

H=H1+H2+H3;

 

dydt=zeroes(5,1);

dydt(1)=reb;

dydt(2)=rc1-rc3+(e_B/rho_B )*(rt1-rt3);

dydt(3)=rc2+rt2*(e_B/rho_B );

dydt(4)=rc3+rt3*(e_B/rho_B );

dydt(5)=(reb*(-H))/(total0*(CPn+X*C));

end

 

 

Not sure what do I put on the command window. I need a graph of conversion(1) against space time, and temperature against space time to work out the catalyst weight. Thanks in advance! 

 

 

[Francisco Angel]

Dear Sharna:

 

From what I see, you have defined a function to calculate the derivatives of the composition (x's):

 

function dydt= @DehyR1 (w,x).. (Remove the @ on the function definition, or it will give an error).

 

Regarding your question: "Not sure what do I put on the command window", you need to use a function like ode45, the usage is like this:

 

[w,x]=ode45(@DehyR1,[wIni wEnd],x0)

 

Note that the first parameter is a function that outputs the derivative values, given the independent variable (w) and the values of the dependent variables (x), the second parameter is the interval of the independent variable used for the calculation, and the third parameter is a vector containing the initial values of the dependent variables.

 

As a result of ode45 you will get

w: a vector containing values of the independent variable,

x: the corresponding values of the dependent variables.

 

This can be implement both on the command window or a separate script. After you have w, x from ode45, you can do any postprocessing, like the graphs you mentioned.

 

Best regards.

[sharna]

I have put that,

 

function [w,x] = please_work()

 

%initial condition: kmol/hr

eb0=478; 

st0=0.4842;

h20=0;

bz0=19.66;

tol0=0.4587;

c2h40=0;

ch40=0;

steam0=5258;

T0=620+273; 

p0=1; 

 

[w, x] = ode45(@DehyR1, [0 50], [eb0]);

%plot(w, x(:,14) , 'r'); hold on; plot(w, x(:,15) , 'y');

hold on; plot(w, x(:,1) , 'b'); hold off;

end

 

doesn't work!

[Francisco Angel]

Dear Sharna:

I tried check your code, but the version of the function "DehyR1" have a number of issues, missing parenthesis, undefined variables, etc.

Please be more specific about the error matlab gives to you, and attach the version of the code you are currently working on.

Best regards.