Reaction Engineering
Reaction Engineering
(OP)
Hi
I am new with Matlab and trying to solve this problem.
The problem is a mutli-reaction system of the form
dc1/dt=kc1+kc2
dc2/dt=kc2-kc3
... etc.
where k=reaction constant
Here is a code I wrote in Matlab:
function dc=bio(t,c)
T=298
E=13145 %constant
R=1.98 % constant
A=1.4e5 %constant
k=A*e^(-E/RT) %(note K is a function of T)
dc=[dc1;dc2]
Then I made another m.file
t=0
tf=20
initial=[1 0]
[t,c]=ode45(@bio,[t0 tf],initial);
So far everything is fine and the model works.
Now I like to change the T (over a wide range)and
run the model again. I tried a few things but none worked. Can someone please tell me how I can do that?
I am new with Matlab and trying to solve this problem.
The problem is a mutli-reaction system of the form
dc1/dt=kc1+kc2
dc2/dt=kc2-kc3
... etc.
where k=reaction constant
Here is a code I wrote in Matlab:
function dc=bio(t,c)
T=298
E=13145 %constant
R=1.98 % constant
A=1.4e5 %constant
k=A*e^(-E/RT) %(note K is a function of T)
dc=[dc1;dc2]
Then I made another m.file
t=0
tf=20
initial=[1 0]
[t,c]=ode45(@bio,[t0 tf],initial);
So far everything is fine and the model works.
Now I like to change the T (over a wide range)and
run the model again. I tried a few things but none worked. Can someone please tell me how I can do that?





RE: Reaction Engineering
=====================================
Eng-tips forums: The best place on the web for engineering discussions.
RE: Reaction Engineering
The problem was when I put the loop in the the "solver" file (I mean the one which I call the ode45), it takes only the first value of T and repeats the calculation. As I said I am new with Matlab so I may be making a stupid mistake!
RE: Reaction Engineering
T=298
in your bio file (derivative file).
You have to do it in your main file. Something like:
global T
t=0
tf=20
initial=[1 0]
T=298
for index = 1:100
T=T+0.1*index
[t,c]=ode45(@bio,[t0 tf],initial);
end
Then your bio file doesn't change T but instead pulls it in as a global variable. Bio.m/ Something lik
dc=bio(t,c)
global T
E=13145 %constant
R=1.98 % constant
A=1.4e5 %constant
k=A*e^(-E/RT) %(note K is a function of T)
dc=[dc1;dc2]
Note it looks like you may have left out a few steps in defining dc1 and dc2 but I assume they can be calculated in the dc=bio() derivative procedure as a funciton of T.
=====================================
Eng-tips forums: The best place on the web for engineering discussions.
RE: Reaction Engineering
% file Main.m
global T
t0=0
tf=20
initial=[1 1]
T=298
for index = 1:3
T=T+10*index
[t,c]=ode45('bio',[t0 tf],initial);
plot(t,c)
pause
end
%file bio.m
function dc=bio(t,c)
global T
E=13145 %constant
R=1.98 % constant
A=1.e5 %constant
k=A*exp(-E/R/T) %(note K is a function of T)
d1=k*c(1)+k*c(2)
d2= k*c(1)-k*c(2)
dc=[d1;d2]
=====================================
Eng-tips forums: The best place on the web for engineering discussions.
RE: Reaction Engineering
Thanks again. It works!
RE: Reaction Engineering
plot(t,c(:,1))
title(['C1, T=' num2str(T)])
pause
plot(t,c(:,2))
title(['C2, T=' num2str(T)])
=====================================
Eng-tips forums: The best place on the web for engineering discussions.
RE: Reaction Engineering
=====================================
Eng-tips forums: The best place on the web for engineering discussions.