×
INTELLIGENT WORK FORUMS
FOR ENGINEERING PROFESSIONALS

Contact US

Log In

Come Join Us!

Are you an
Engineering professional?
Join Eng-Tips Forums!
  • Talk With Other Members
  • Be Notified Of Responses
    To Your Posts
  • Keyword Search
  • One-Click Access To Your
    Favorite Forums
  • Automated Signatures
    On Your Posts
  • Best Of All, It's Free!

*Eng-Tips's functionality depends on members receiving e-mail. By joining you are opting in to receive e-mail.

Posting Guidelines

Promoting, selling, recruiting, coursework and thesis posting is forbidden.

Students Click Here

Reaction Engineering

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?


RE: Reaction Engineering

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.

RE: Reaction Engineering

(OP)
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!

RE: Reaction Engineering

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.

RE: Reaction Engineering

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.

RE: Reaction Engineering

(OP)
electricpete,

Thanks again. It works!

RE: Reaction Engineering

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.

RE: Reaction Engineering

...and another pause at the end

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

Red Flag This Post

Please let us know here why this post is inappropriate. Reasons such as off-topic, duplicates, flames, illegal, vulgar, or students posting their homework.

Red Flag Submitted

Thank you for helping keep Eng-Tips Forums free from inappropriate posts.
The Eng-Tips staff will check this out and take appropriate action.

Reply To This Thread

Posting in the Eng-Tips forums is a member-only feature.

Click Here to join Eng-Tips and talk with other members! Already a Member? Login



News


Close Box

Join Eng-Tips® Today!

Join your peers on the Internet's largest technical engineering professional community.
It's easy to join and it's free.

Here's Why Members Love Eng-Tips Forums:

Register now while it's still free!

Already a member? Close this window and log in.

Join Us             Close