Unable to meet integration tolerance
1 view (last 30 days)
Show older comments
Hi! I'm trying to design an steady-state packed bed reactor, with a nonlinear second-order differential equation, but i'm getting this "Warning: Failure at t=2.455844e-06. Unable to meet integration tolerances without reducing the step size below the smallest value allowed (6.776264e-21) at time t."
I would really appreciate your time and help, thanks.
The script is:
clear all; clc
global U Deff ro_b k
%butanotiol+heptano
%T=25ºC
Deff=0.046E-8 %m2/s
ro_b=167.9400;
L=4.27;
U=0.0024;
k=4.68E-2
Ca0=6.05;
%ODE
z=0:0.1:L;
z1=z';
y0=[Ca0; 0]
[z,y]=ode15s(@react,z1,y0)
plot(z,y(:,1))
ylabel('Concentración Ca [mol/m3]')
xlabel('z [m]')
and the function is:
function F=react(z,y)
global U Deff ro_b k
A=y(2);
B=(U./Deff).*y(2)-k.*ro_b.*(y(1).^2)./Deff;
F=[A; B];
end
2 Comments
Ced
on 21 Oct 2014
Edited: Ced
on 21 Oct 2014
I'm not a chemist, but are you sure your equations are correct? Looking at the change of B, it is initialized at zero, and then pulled down to minus infinity. The reason is that no matter the sign of y(1), B becomes smaller. And being initialized at zero, y(2) itself will always be smaller or equal to zero. There is no way this can be stable, hence the integration fails. Also, the difference in scale between the different terms is bound to bring numerical issues, even if the equations are correct.
On a sidenote: It's easier to read if you use e.g. dy and y than changing the names from y to A,B to F. Also, global variables are really not pretty in my opinion. A better option would be:
% script
clear all; clc
%butanotiol+heptano
%T=25ºC
Deff=0.046E-8 %m2/s
ro_b=167.9400;
L=4.27;
U=0.0024;
k=4.68E-2
Ca0=6.05e8;
%ODE
z=0:0.1:L;
z1=z';
y0=[Ca0; 0]
ode_react = @(z,y)react(z,y,U,Deff,ro_b,k);
[z,y]=ode15s(ode_react,z1,y0)
plot(z,y(:,1))
ylabel('Concentración Ca [mol/m3]')
xlabel('z [m]')
end
% and your function
function dy=react(z,y,U,Deff,ro_b,k)
dy = zeros(2,1);
dy(1)=y(2);
dy(2)=(U./Deff).*y(2)-k.*ro_b.*(y(1).^2)./Deff;
end
Answers (0)
See Also
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!