Fixed Bed Adsorption Model

I have to model a fixed bed adsorption to obtain the breakthrough curve of the adsorption. I tried to code the model but i think the initial and boundary condition in the code are not correct. Beside that, I got this error when i ran this code.

Warning: Failure at t=2.006182e+02. Unable to meet integration tolerances without reducing the step size below the

smallest value allowed (4.547474e-13) at time t.

> In ode15s (line 655)

In fyp (line 19)

I am not really good with MATLAB, so I would really appreciate any of your help. Thank you very much. I have also attached the model's equations in the file.

Cfeed = 0.0224;     % Inlet concentration

tf = 2000;      % End time
nz = 200;       % Number of gateways
np = nz + 1;    % Number of gateposts

% Initial Concentrations
Cb = zeros(np,1); Cb(1) = Cfeed;  
Cp = zeros(np,1);

% Prepare data to pass to ode solver
c0 = [Cb; Cp];   
dt = tf/100;
tspan = 0:dt:tf;

% Call ode solver
[t, c] = ode15s(@rates, tspan, c0);

% Plot results
plot(t, c(:,np)/Cfeed),grid
xlabel('time'), ylabel('Cb/Cfeed')
title('Figure 1: Exit Cb/Cfeed vs time')
figure
plot(t, c(:,2*np)),grid
xlabel('time'), ylabel('Cp')
title('Figure 2: Exit Cp vs time')
figure
plot(t, c(:,np)-c(:,2*np)), grid
xlabel('time'), ylabel('Cb - Cp')
title('Figure 3: Exit Cb-Cp vs time')

k = 598;
q = k*c(:,2*np);  % Freundlcih equation
figure
plot(t,q), grid
xlabel('time'), ylabel('q')
title('Figure 4: Exit q vs time')

% Function to calculate rate of change of concentrations with time
function dcdt = rates(~, c)
    a=0.00000000261;
    b=0.2102;
    d=1705.4719;
    e=-38.8576;
    L=0.08;
    nz = 200;       % Number of gateways
    np = nz + 1;    % Number of gateposts
    dz = L/nz;

    Cb = c(1:np);
    Cp = c(np+1:end);
    dCbdt = zeros(np,1);
    dCpdt = zeros(np,1);
    
   for i = 2:np-1
       dCbdt(i) = (a/dz^2)*(Cb(i+1)-2*Cb(i)+Cb(i-1)) - b/(2*dz)*(Cb(i+1)-Cb(i-1)) - d*(Cb(i)-Cp(i));
       dCpdt(i) = e*(Cb(i)-Cp(i));  % Needed to change the sign from -ve to +ve here.
   end % i loop
   
   % For simplicity, set exit rates to those of the immediately previous nodes
    dCbdt(np) = dCbdt(np-1); 
    dCpdt(np) = dCpdt(np-1);
    
    % Combine rates for the function to return to the calling routine
    dcdt = [dCbdt; dCpdt];

end % of rates function 

[Brandon]

Can I ask you, how did you derive the differential equations that you are using in your model?

Answer 1

[Torsten]

Cfeed = 0.0224;     % Inlet concentration
tf = 3000;      % End time
nz = 400;       % Number of gateways
np = nz + 1;    % Number of gateposts
% Initial Concentrations
Cb = zeros(np,1); Cb(1) = Cfeed;
Cp = zeros(np,1);
% Prepare data to pass to ode solver
c0 = [Cb; Cp];
dt = tf/200;
tspan = 0:dt:tf;
% Call ode solver
[t, c] = ode15s(@rates, tspan, c0);
% Plot results
figure(1)
plot(t, c(:,np)/Cfeed),grid
xlabel('time'), ylabel('Cb/Cfeed')
title('Figure 1: Exit Cb/Cfeed vs time')
figure(2)
plot(t, c(:,2*np)),grid
xlabel('time'), ylabel('Cp')
title('Figure 2: Exit Cp vs time')
figure(3)
plot(t, c(:,np)-c(:,2*np)), grid
xlabel('time'), ylabel('Cb - Cp')
title('Figure 3: Exit Cb-Cp vs time')
k = 598;
q = k*c(:,2*np);  % Freundlcih equation
figure(4)
plot(t,q), grid
xlabel('time'), ylabel('q')
title('Figure 4: Exit q vs time')
L=0.08;
dz = L/nz;
z=0:dz:L;
figure(5)
plot(z,[c(1,1:np);c(2,1:np);c(3,1:np);c(4,1:np);c(5,1:np)])
figure(6)
plot(z,[c(1,np+1:2*np);c(2,np+1:2*np);c(3,np+1:2*np);c(4,np+1:2*np);c(5,np+1:2*np)])

% Function to calculate rate of change of concentrations with time
function dcdt = rates(t, c)
a=0.00000000261;
b=0.2102;
d=1705.4719;
e=8.8576;
L=0.08;
nz = 400;       % Number of gateways
np = nz + 1;    % Number of gateposts
dz = L/nz;
Cb = c(1:np);
Cp = c(np+1:end);
dCbdt = zeros(np,1);
dCpdt = zeros(np,1);

for i = 2:np-1
    %dCbdt(i) = (a/dz^2)*(Cb(i+1)-2*Cb(i)+Cb(i-1)) - b/(2*dz)*(Cb(i+1)-Cb(i-1)) - d*(Cb(i)-Cp(i));
    dCbdt(i) = (a/dz^2)*(Cb(i+1)-2*Cb(i)+Cb(i-1)) - b/dz*(Cb(i)-Cb(i-1)) - d*(Cb(i)-Cp(i));
    dCpdt(i) = e*(Cb(i)-Cp(i));  % Needed to change the sign from -ve to +ve here.
end % i loop
% For simplicity, set exit rates to those of the immediately previous nodes
%dCbdt(np) = dCbdt(np-1);
%dCpdt(np) = dCpdt(np-1);
dCbdt(np) = -2*a*(Cb(np)-Cb(np-1))/dz^2 - d*(Cb(np)-Cp(np));
dCpdt(np) = e*(Cb(np)-Cp(np));

% Combine rates for the function to return to the calling routine
dcdt = [dCbdt; dCpdt];
t
end % of rates function

Concerning your equations: Since you don't solve equations for Cp within the particle, but only define an overall mass transfer coefficient at r=R, you don't need the conditions at r=0 and r=R.

[Puteri Ellyana Mohmad Latfi]

Hi @Torsten, I tried running your code but somehow the Cb/Cp vs time graph shows that it takes a long time to reach 1. Not sure if that's correct tho

[Torsten]

This is not a question of the code, but of your constants a,b,d and e.

[Puteri Ellyana Mohmad Latfi]

But it kinda look different from the experimental result as shown in the attached picture (10,000ppm). the breakthrough time is similar but the time taken for it to reach 1 is different.

Cfeed = 0.0227;     % Inlet concentration
tf = 250;      % End time
nz = 400;       % Number of gateways
np = nz + 1;    % Number of gateposts

% Initial Concentrations
Cb = zeros(np,1); Cb(1) = Cfeed;
Cp = zeros(np,1);

% Prepare data to pass to ode solver
c0 = [Cb; Cp];
dt = tf/200;
tspan = 0:dt:tf;

% Call ode solver
[t, c] = ode15s(@rates, tspan, c0);

% Plot results
figure(1)
plot(t, c(:,np)/Cfeed),grid
xlabel('time'), ylabel('Cb/Cfeed')
title('Figure 1: Exit Cb/Cfeed vs time')

% figure(2)
% plot(t, c(:,2*np)),grid
% xlabel('time'), ylabel('Cp')
% title('Figure 2: Exit Cp vs time')
% 
% figure(3)
% plot(t, c(:,np)-c(:,2*np)), grid
% xlabel('time'), ylabel('Cb - Cp')
% title('Figure 3: Exit Cb-Cp vs time')
% 
% k = 56.5;
% q = k*c(:,2*np);  % Freundlcih equation
% figure(4)
% plot(t,q), grid
% xlabel('time'), ylabel('q')
% title('Figure 4: Exit q vs time')
% 
% L=0.065;
% dz = L/nz;
% z=0:dz:L;
% figure(5)
% plot(z,[c(1,1:np);c(2,1:np);c(3,1:np);c(4,1:np);c(5,1:np)])
% figure(6)
% plot(z,[c(1,np+1:2*np);c(2,np+1:2*np);c(3,np+1:2*np);c(4,np+1:2*np);c(5,np+1:2*np)])

% Function to calculate rate of change of concentrations with time
function dcdt = rates(t, c)
a=0.00000001278;
b=0.16056;
d=244.9429;
e=1.3964;
L=0.065;

nz = 400;       % Number of gateways
np = nz + 1;    % Number of gateposts
dz = L/nz;
Cb = c(1:np);
Cp = c(np+1:end);
dCbdt = zeros(np,1);
dCpdt = zeros(np,1);

for i = 2:np-1
    %dCbdt(i) = (a/dz^2)*(Cb(i+1)-2*Cb(i)+Cb(i-1)) - b/(2*dz)*(Cb(i+1)-Cb(i-1)) - d*(Cb(i)-Cp(i));
    dCbdt(i) = (a/dz^2)*(Cb(i+1)-2*Cb(i)+Cb(i-1)) - b/dz*(Cb(i)-Cb(i-1)) - d*(Cb(i)-Cp(i));
    dCpdt(i) = e*(Cb(i)-Cp(i));  % Needed to change the sign from -ve to +ve here.
end % i loop
% For simplicity, set exit rates to those of the immediately previous nodes
%dCbdt(np) = dCbdt(np-1);
%dCpdt(np) = dCpdt(np-1);
dCbdt(np) = -2*a*(Cb(np)-Cb(np-1))/dz^2 - d*(Cb(np)-Cp(np));
dCpdt(np) = e*(Cb(np)-Cp(np));

% Combine rates for the function to return to the calling routine
dcdt = [dCbdt; dCpdt];
end % of rates function

[Torsten]

Where are your adsorption equilibrium curves ?

[Puteri Ellyana Mohmad Latfi]

from the article, there is no adsorption equilibrium curves, only the adsorption breakthrough curve, which is supposed to be an S-shaped curve.

[Torsten]

And how did you come up with the values for d and e ? Maybe you should run an optimization to fit them to the experimental breakthrough curves of the article.

[Torsten]

You might replace

dCbdt(np) = -2*a*(Cb(np)-Cb(np-1))/dz^2 - d*(Cb(np)-Cp(np));

by

dCbdt(np) = -b/dz*(Cb(np)-Cb(np-1)) - 2*a/dz^2*(Cb(np)-Cb(np-1)) - d*(Cb(np)-Cp(np));

Although the usual exit condition is dC/dz = 0 and -b/dz*(Cb(np)-Cb(np-1)) should be set to 0 in this case, the revised setting gives a better profile for Cb at the exit. I don't know exactly why this is the case and whether the inclusion of the term is justified.

[Puteri Ellyana Mohmad Latfi]

I think it's probably because the original code tries to use Freundlich equation, which is a thermodynamically inconsistent isotherm as the concentration, Cp goes to zero the gradient goes to infinity, so it creates a numerical stiffness at the leading edge of the breakthrough curve. However, I want to use linear isotherm. so maybe that's why. but im not sure how to recorrect it though.

Btw, d is a calculated constant while e is an adjustable parameter, i've already optimized the e value though.

[Torsten]

And are you satisfied with the breakthrough curve now after the modification of dCbdt(np) ?

[Puteri Ellyana Mohmad Latfi]

it does give a nice S-graph and the correct breakthrough time but the time for it to reach 1 is different from the graph in the article as shown in the attached file.

[Torsten]

So you manually changed the coefficient e depending on temperature, I guess ?

To better reproduce the breakthrough curves, you might want to optimize e using lsqcurvefit, e.g.

[Puteri Ellyana Mohmad Latfi]

the other constant changes as well but e can be adjusted to get the right breakthrough curve. do you mind showing me how to do that?

[Torsten]

time = [0 40 50 65 75 90 110];
cnorm = [0 0 0.01 0.2 0.6 0.95 1.0];
p0 = 2.0;
info = 0;
sol = lsqnonlin(@(p)fit(p,time,cnorm,info),p0);
info = 1;
res = fit(sol,time,cnorm,info)

function res = fit(p,time,cnorm,info)
  format long
  p
  Cfeed = 0.0227;     % Inlet concentration
  tf = 250;      % End time
  nz = 100;       % Number of gateways
  np = nz + 1;    % Number of gateposts
  % Initial Concentrations
  Cb = zeros(np,1); Cb(1) = Cfeed;
  Cp = zeros(np,1);
  % Prepare data to pass to ode solver
  c0 = [Cb; Cp];
  if info==0
    tspan = time;
  elseif info==1
    tf = 250;      % End time
    dt = tf/200;
    tspan = 0:dt:tf;
  end
  % Call ode solver
  [t, c] = ode15s(@(t,c)rates(t,c,p), tspan, c0);
  if info==0
    res = cnorm.' - c(:,np)/Cfeed
  elseif info==1 
    res = cnorm - interp1(t,c(:,np),time)/Cfeed
  end
  if info == 1
    % Plot results
    figure(1)
    plot(t, c(:,np)/Cfeed),grid
    hold on
    plot(time,cnorm,'*')
    xlabel('time'), ylabel('Cb/Cfeed')
    title('Figure 1: Exit Cb/Cfeed vs time')
    L=0.065;
    dz = L/nz;
    z=0:dz:L;
    figure(2)
    plot(z,[c(40,1:np);c(80,1:np);c(120,1:np);c(160,1:np);c(200,1:np)])
    figure(3)
    plot(z,[c(40,np+1:2*np);c(80,np+1:2*np);c(120,np+1:2*np);c(160,np+1:2*np);c(200,np+1:2*np)])
  end
end
% Function to calculate rate of change of concentrations with time
function dcdt = rates(t, c, e)
  persistent iter
  if isempty(iter)
    iter = 0;
  end
  iter = iter + 1;
  a=0.00000001278;
  b=0.16056;
  d=244.9429;
  %e=1.3964;
  L=0.065;
  nz = 100;       % Number of gateways
  np = nz + 1;    % Number of gateposts
  dz = L/nz;
  Cb = c(1:np);
  Cp = c(np+1:end);
  dCbdt = zeros(np,1);
  dCpdt = zeros(np,1);
  for i = 2:np-1
    %dCbdt(i) = (a/dz^2)*(Cb(i+1)-2*Cb(i)+Cb(i-1)) - b/(2*dz)*(Cb(i+1)-Cb(i-1)) - d*(Cb(i)-Cp(i));
    dCbdt(i) = (a/dz^2)*(Cb(i+1)-2*Cb(i)+Cb(i-1)) - b/(2*dz)*(Cb(i+1)-Cb(i-1)) - d*(Cb(i)-Cp(i));
    dCpdt(i) = e*(Cb(i)-Cp(i));  % Needed to change the sign from -ve to +ve here.
  end % i loop
  % For simplicity, set exit rates to those of the immediately previous nodes
  %dCbdt(np) = dCbdt(np-1);
  %dCpdt(np) = dCpdt(np-1);
  dCbdt(np) = -b/dz*(Cb(np)-Cb(np-1)) - 2*a/dz^2*(Cb(np)-Cb(np-1)) - d*(Cb(np)-Cp(np));
  dCpdt(np) = e*(Cb(np)-Cp(np));
  % Combine rates for the function to return to the calling routine
  dcdt = [dCbdt; dCpdt];
  if mod(iter,100)==0
    disp(t)
    iter = 0;
  end
end % of rates function