Reaction Engineering

[enaynad]

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?

 

[electricpete]

It's not clear what your problem is. I think you need to declare T as a global variable in both the main program and the derivative program (you probably already did that). And then create a loop which calls the procedure repeatedly. Then you can either printout the results or write them to a file during the process. Perhaps you could also store the t and c arrays for each loop in a multidimensional array where one index is the index of the loop. But likely the arrays may not all have the same size so you'll have to pay attention to that (may need a longer array with some zero-padding on the shorter results)

=====================================
Eng-tips forums: The best place on the web for engineering discussions.

 

[enaynad]

Thank you much for the reply. I guess I was not clear about my problem. As you mentioned, I wanted to make a loop to change T (and subsequently k) and then solve the ode and finally come up with the profile of c based on the new T.

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!

 

[electricpete]

You can't assign
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.

 

[electricpete]

Try this:

% 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.

 

[enaynad]

electricpete,

Thanks again. It works!

 

[electricpete]

you're welcome. I realized my plot didn't plot both curves. A better set of plots would be:
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.

 

[electricpete]

...and another pause at the end

=====================================
Eng-tips forums: The best place on the web for engineering discussions.