MATLAB Answers

0

Error when using ode15

Asked by Can Cakiroglu on 23 May 2019
Latest activity Commented on by Can Cakiroglu on 24 May 2019
Hello,
I'm a chemical engineering student, so these equations might come difficult to analyze. However, I am receiving following error on my equations
This is the main script
close all
clear
clc
global Fao Tr z a b1 b2 R K30 E2 E3 omega alpha k20 Bo rho_c Phi Ac T0 P0
Fao = 1487.351; % kmol/h
Tr = 273; % K
z = 1;
a = (3 * 1.5) / ((3 * 1.5) + 1);
b1 = (0.5 + (0.5 / 1.5)) / 1;
b2 = (-0.5+(1.5*1.5))/1;
R = 8.314 / 1000; % kJ/mol*K
K30 = 0.0307;
E2 = 72100.82; % kJ/mol
E3 = -80989.98; % kJ/mol
omega = 1.564;
alpha = 0.640;
k20 = 2.12 * (10 ^ 13);
Bo = 0.254626; %atm/m
rho_c = 1600; %kg/m^3
Phi = 0.5;
D = 0.05; %Diameter of tube, m
Ac = pi * (D^2)/4; %m^2
X0 = 0;
T0 = 673; % K
P0 = 150; %atm
x0 = [X0 T0 P0];
Wspan = [0 5];
[W, x] = ode15s(@fun, Wspan, x0);
And this is the function used
function func = fun(x)
%% Identification of global variables
global Fao Tr z a b1 b2 R K30 E2 E3 omega alpha k20 Bo rho_c Phi Ac T0 P0
%% Decleration of the non-linear equations
func = [((k20*(-E2/R*x(2))*(x(3)*(a/3)*(1-b2*z)*((10^(-5.519265*(10^(-5))*x(2)+...
1.848863*(10^(-7))*(x(2)^2)+(2001.6/x(2))+2.6899)/(x(2)^(5.519265)))^2)-...
(((x(3)*z)^2)/((x(3)*a*(1-b1*z))^3))))/((1+(K30*exp(-E3/R*x(2)))*(x(3)*z)/...
((x(3)*a*(1-b1*z))^omega))^(2*alpha)))/Fao; % dX/dW, x(1) = X
((2.12*(10^13)*(-E2/R*x(2))*(x(3)*(a/3)*(1-b2*z)*((10^(-5.519265*(10^(-5))*x(2)...
+ 1.848863*(10^(-7))*(x(2)^2)+(2001.6/x(2))+2.6899)/(x(2)^(5.519265)))^2)-...
(((x(3)*z)^2)/((x(3)*a*(1-b1*z))^3))))/((1+(K30*exp(-E3/R*x(2)))*(x(3)*z)/...
((x(3)*a*(1-b1*z))^omega))^(2*alpha)))/(-(-91.63+(-12.96*(x(2)-Tr)+0.00917*((x(2)-Tr)^2)/2)))/...
(Fao*(6.50+0.00100*x(2))+x(1)*2*(6.70+ 0.00630*x(2))-...
3*(6.62+0.00081*x(2))-(6.50+0.00100*x(2))); % dT/dW, x(2) = T
-(Bo/(Ac*(1-Phi)*rho_c))*(P0/x(3))*(x(2)/T0)*(1-0.5*x(1)) %dp/dW, x(3) = P
];
end
and I receive following error when i start the script:
Error in ode15s (line 150)
odearguments(FcnHandlesUsed, solver_name, ode, tspan, y0, options, varargin);
Error in Reactor_Project (line 40)
[W, x] = ode15s(@fun, Wspan, x0);
I would be very happy if someone could identify this error. I thank everybody who spends time on this matter.

  0 Comments

Sign in to comment.

Tags

Products


Release

R2019a

1 Answer

Answer by David Wilson on 24 May 2019
 Accepted Answer

It is good practice to try and avoid globals, but rather pass them as auxilary variables into the function.
Fao = 1487.351; % kmol/h
Tr = 273; % K
z = 1;
a = (3 * 1.5) / ((3 * 1.5) + 1);
b1 = (0.5 + (0.5 / 1.5)) / 1;
b2 = (-0.5+(1.5*1.5))/1;
R = 8.314 / 1000; % kJ/mol*K
K30 = 0.0307;
E2 = 72100.82; % kJ/mol
E3 = -80989.98; % kJ/mol
omega = 1.564;
alpha = 0.640;
k20 = 2.12 * (10 ^ 13); % should be 2.12e13 I suspect
Bo = 0.254626; %atm/m
rho_c = 1600; %kg/m^3
Phi = 0.5;
D = 0.05; %Diameter of tube, m
Ac = pi * (D^2)/4; %m^2
X0 = 0;
T0 = 673; % K
P0 = 150; %atm
x0 = [X0 T0 P0];
Wspan = [0 5];
[W, x] = ode15s(@(t,x) fchem(t,x,Fao, Tr, z, a, b1, b2, R, K30, E2, E3, omega, alpha, k20, Bo, rho_c, Phi, Ac, T0, P0) ...
, Wspan, x0);
plot(W,x)
Now you function is (which I've renamed to make it more user-friendly)
function func = fchem(t,x,Fao, Tr, z, a, b1, b2, R, K30, E2, E3, omega, alpha, k20, Bo, rho_c, Phi, Ac, T0, P0)
%% Identification of global variables
%global Fao Tr z a b1 b2 R K30 E2 E3 omega alpha k20 Bo rho_c Phi Ac T0 P0
%% Decleration of the non-linear equations
func = [((k20*(-E2/R*x(2))*(x(3)*(a/3)*(1-b2*z)*((10^(-5.519265*(10^(-5))*x(2)+...
1.848863*(10^(-7))*(x(2)^2)+(2001.6/x(2))+2.6899)/(x(2)^(5.519265)))^2)-...
(((x(3)*z)^2)/((x(3)*a*(1-b1*z))^3))))/((1+(K30*exp(-E3/R*x(2)))*(x(3)*z)/...
((x(3)*a*(1-b1*z))^omega))^(2*alpha)))/Fao; % dX/dW, x(1) = X
((2.12*(10^13)*(-E2/R*x(2))*(x(3)*(a/3)*(1-b2*z)*((10^(-5.519265*(10^(-5))*x(2)...
+ 1.848863*(10^(-7))*(x(2)^2)+(2001.6/x(2))+2.6899)/(x(2)^(5.519265)))^2)-...
(((x(3)*z)^2)/((x(3)*a*(1-b1*z))^3))))/((1+(K30*exp(-E3/R*x(2)))*(x(3)*z)/...
((x(3)*a*(1-b1*z))^omega))^(2*alpha)))/(-(-91.63+(-12.96*(x(2)-Tr)+0.00917*((x(2)-Tr)^2)/2)))/...
(Fao*(6.50+0.00100*x(2))+x(1)*2*(6.70+ 0.00630*x(2))-...
3*(6.62+0.00081*x(2))-(6.50+0.00100*x(2))); % dT/dW, x(2) = T
-(Bo/(Ac*(1-Phi)*rho_c))*(P0/x(3))*(x(2)/T0)*(1-0.5*x(1)) %dp/dW, x(3) = P
];
end
When I plot the solution, it looks pretty boring, but that's a different problem.

  1 Comment

Thank you very much for your hard work. The pretty boring part is indeed so, but hopefuly will be fixed since this is our final project calculations. :)

Sign in to comment.