Solving a set of differential equations for a packed bed reactor - chemical reaction engineering

Hello,

I am trying to solve a set of differential equations for conversion and pressure drop with respect to catalyst weight within a PBR, with a rather complex rate law, which I believe is causing MATLAB to throw an error:

"ode_sys

Not enough input arguments.

Error in ode_sys (line 8)

X=var(1);"

Here is my code for the explicit equations and differential equations which for the ode_sys.

% Model solves differential equations for conversion and pressure drop in
% PBR with respect to catalyst weight.

function f=ode_sys(W,var)

% var = [0.5 0.5];

X=var(1);
y=var(2);

%% Parameters (Excplicits)

% Initial Conditions
P0=5; % kPa, Initial Ammonia Pressure CHECK UNITS 50 mbar
T0=500+273; % K, Initial Temperature
Fa0=1.67e-4; % molA/s, Initial Ammonia Molar Flow
ya0 = 1; % Mole Fraction Ammonia at Inlet (Pure Ammonia)

% Catalyst
rhob=1.2e3; % kg/m^3
phi=0.5; % Porosity
rhoc=rhob/(1-phi); % kg/m^3, Catalyst Density
Dp=1.5e-3; % m, Particle Diameter

% Reactor Tube
Dt=21e-3; % m, Tube Internal Diameter 25OD 2mm wall thickness
Ac=((pi)*(Dt)^2)/4; % m^2, Tube Cross Sectional Area

% Gas Properties
R=8.314; % kPa.dm^3/(mol.K), Gas Constant
mu=1.7e-5; % kg/m.s, Gas Viscosity @ To, Po
rhog=0.434; % kg/m^3, Gas Density @ To, Po
m=Fa0*17.01*3600; % kg/hr, Gas Mass Flowrate CHECK FAO MATCHES
G=(m/3600)/Ac; % kg/m^2.s, Superficial Mass Velocity
Ca0=(ya0)*((P0)/(R*T0)); % mol/dm^3, Initial Ammonia Concentration

% Ergun Parameters
beta=((G*(1-phi))/((rhog)*(Dp)*(phi)^3))*(((150*(1-phi)*mu)/Dp)+1.75*G); % kPa/m, Pressure Drop Constant 
alpha=(2*beta)/(Ac*rhoc*(1-phi)*P0); % kPa^-1, Pressure Drop Parameter

% Rate Constants (From Literature)
ka=3.615e-4; % mol/gcatalyst.s, Rate Constant
Kad=6.293e-4*1.01e2; % kPa^-1, Ammonia Adsorption Constant
Khd=3.937e-3*1.01e2; % kPa^-1, Hydrogen Adsorption Constant

% Stoichiometry
Pa=Ca0*R*T0*((1-X)/(1+X))*y; % kPa, Ammonia Partial Pressure as f(X,y) From Stoic  
Ph=Ca0*R*T0*((1.5*X)/(1+X))*y; % kPa, Hydrogen Partial Pressure as f(X,y) From Stoic

% y=(P/P0) Pressure Drop Ratio

% Rate Law
raprime=-(ka*(Kad^2)*(Pa^2))/((Kad*Pa)+((Khd^(3/2))*(Ph^(3/2)))*(1+sqrt(Khd*Ph)))^2;  % molA reacted/gcatalyst.s, Non-Elementary Reaction Rate Law (Literature)

%% Differentials

%dX/dW
dXdW = -raprime/Fa0;
%dY/dW
dydW = -(alpha*(1+X))/(2*y);

f = [dXdW; dydW];

end

And here is the code I am currently using to solve the equations:

clc
Wspan = [0 10]; % Range for independant variable, W
var0 = [0 0];
[w, var]=ode45(@ode_sys,Wspan,var0);

Answer 1

[Star Strider]

I did not run the code, however the problem is likely that using ‘var’ as a variable name is overshadowing the var (variance) function. That would conform to the error message.

The solution may be as easy as naming it something else, perhaps to a more generic ‘v’.

[William Entwistle]

Tried this suggestion, no errors now, solution seems to be NaN....

[Star Strider]

The NaN values are due either to NaN values in data, 0/0 or Inf/Inf somewhere in the code.

When I try to run it in the opnline Run feature, I get a message that simply says ‘Something went wrong. Please try again’. All the initial conditions are 0 and the first value of ‘W’ is zero, so I though those could be the problem. Adding 1E-8 to ‘var0’ that usually solves these, did not this time. I do not see any obvious problems and the error message is not helpful. I defer to you to discover the reason for the NaN values.

Nevertheless,I did manage to solve the posted problem.

.

[William Entwistle]

I have managed to solve the problem, thank you.

[Star Strider]

My pleasure.

I solved the original problem you requested help with.

[Wen Hong Teo]

Hi Mr William. Good day. I am also trying to create the similar model. However, I am facing problem on the codings. May I know how you managed to solve the problem and obtain successful model? I hope you don't mind sharing your codes if possible. Thanks a lot in advance.