Eng-Tips is the largest engineering community on the Internet

Intelligent Work Forums for Engineering Professionals

  • Congratulations waross on being selected by the Tek-Tips community for having the most helpful posts in the forums last week. Way to Go!

Matlab coding for a PBR with pressure drop and adiabatic operation problems

Status
Not open for further replies.

anthony0028

Chemical
Mar 15, 2015
2
0
0
US
I am having issues getting this code to run. I actually received this code back in college when taking reactor design, now it will not run. I am receiving several errors which I have failed to resolve. what is wrong with the coding?

here are the errors when I run the ODE45 script:

>> exampleGas_graphs
Error using odearguments (line 90)
EXAMPLEGASPD_X_2013 must return a column vector.

Error in ode45 (line 113)
[neq, tspan, ntspan, next, t0, tfinal, tdir, y0, f0, odeArgs, odeFcn, ...

Error in exampleGas_graphs (line 1)
[w,y]=ode45('ExampleGasPD_X_2013',[0 70],[1 0]);



The function file is as follows:
--------------------------------------
function dy=ExampleGasPD_X_2013(w,y);
dy=zeros(1,0);
alpha=0.01;
vo=10;
Fao=2; Fbo=5; Fto=Fao+Fbo;
T1=(73000*y(2)+53550);
T2=105+30*y(2);
T=(73000*y(2)+53550)/(105+(30*y(2)));
x=y(2);
Cto=Fto/vo;
ko=0.001; Ea=10000; R1=1.987;
To=300;
kprime=ko*exp((Ea/R1)*((1/To)-(1/T)));
Ca=Cto*(1-y(2))*y(1)/(3.5-y(2))*To/T;
Cb=Cto*(2.5-2*y(2))*y(1)/(3.5-y(2))*To/T;
Cc=Cto*(2*y(2))*y(1)/(3.5-y(2))*To/T;
rate=kprime*(Ca*Cb*Cb);
dy(1)=-(alpha*(3.5-y(2))/2/y(1))*To/T;
dy(2)=rate/Fao;
-------------------------------------

The ODE45 file is as follows:
-------------------------------------
[w,y]=ode45('ExampleGasPD_X_2013',[0 70],[1 0]);
figure; plot(w,y:),1));
xlabel('w');
ylabel('P/Po');
figure; plot(w,y:),2));
xlabel('W');
ylabel('X');

T=T1/T2;
figure; plot(y:),2),T);
xlabel('w');
ylabel('T');
--------------------------------------
 
Status
Not open for further replies.
Back
Top