Why is my code not working how I want it to?

Good News Everyone!

My matlab code is malfunctioning and I have no idea why... So I am trying to convert this Polymath code into a Matlab code and it is not working out. I attach the Polymath code and what is supposed to come out but I am just not getting it. Any help will be greatly appreciated. :D

This is my code:

function nonisothermal
vspan=[0 1]; %volume of the reactor dm3
x0=[0 1035]; %Conversion and Temperature(in K)
[v x]=ode45(@thermalfun,vspan,x0);
subplot(2,1,1)
plot(v,x(:,1),'m--','linewidth',1)
xlabel('Volume of the Reactor')
title('Conversion')
grid on
subplot(2,1,2)
plot(v,x(:,2),'k--','linewidth',1)
xlabel('Volume of the Reactor')
title('Temperature')
grid on
end
function dXdV=thermalfun(v,x)
Cao=18.8;
Fao=0.0376;
Cpa=163;
delCp=-9;
To=1035;
delH=80770+delCp*(x(2)-298);
ra=-Cao*3.58*exp(34222*(1/To-1/x(2)))*(1-x(1))*(To/x(2))/(1+x(1));
Ua=165000;
Ta=1150;
dXdV=[-ra/Fao; (Ua*(Ta-x(2))+ra*delH)/(Fao*(Cpa+x(1)*delCp))]
end

[Star Strider]

I do not see any significant problems with your code, with respect to the published equations. There is an insignificant problem in that ‘Ua’ should be 16500, and that the second plot is ‘-ra’ as a function of ‘v’, so you need to calculate ‘ra’ outside the ode45 call with ‘x(:,1)’, ‘x(:,2)’ and the relevant constants.

Check to see the ODE solver the authors used. The MATLAB numerical solvers are the best I’ve encountered, however some papers use much less sophisticated solvers, leading to results that are difficult to reproduce with the MATLAB solvers. (I invariably trust the MATLAB results.) Also experiment with a much shorter upper limit for ‘vspan’.

Answer 1

[MHZ]

I tried out your code and ran it fine. First make sure to save 2 files in a location in matlab path. Just creat a folder to save in and then go to Matlab home> set path> add folder then hit save and close. Although you do not need to save the first one as a function and can do a normal m-script, save the following as nonisothermal

function nonisothermal
vspan=[0 1]; %volume of the reactor dm3
x0=[0 1035]; %Conversion and Temperature(in K)
[v x]=ode45(@thermalfun,vspan,x0);
subplot(2,1,1)
plot(v,x(:,1),'m--','linewidth',1)
xlabel('Volume of the Reactor')
title('Conversion')
grid on
subplot(2,1,2)
plot(v,x(:,2),'k--','linewidth',1)
xlabel('Volume of the Reactor')
title('Temperature')
grid on
end

Now save the following as thermalfun

function dXdV=thermalfun(v,x)
Cao=18.8;
Fao=0.0376;
Cpa=163;
delCp=-9;
To=1035;
delH=80770+delCp*(x(2)-298);
ra=-Cao*3.58*exp(34222*(1/To-1/x(2)))*(1-x(1))*(To/x(2))/(1+x(1));
Ua=165000;
Ta=1150;
dXdV=[-ra/Fao; (Ua*(Ta-x(2))+ra*delH)/(Fao*(Cpa+x(1)*delCp))]
end

When done, go to the command window and type nonisothermal and hit enter.

[Marlon Brutus]

I don't think you understand... The code is supposed to come with results that are pretty similar to what is attached above.

[MHZ]

Your problem is in the implementation of the integration equation. You are using x(2) instead of T. According to Matlab, x(2) will be 1035 which is your final and not the value of T. You need to separate dT/dV and actually use T.