Fixed-bed adsorption

Hello everybody, I am trying to obtain breakthrough curves from my matlab code, but probably there are some errors that I am not able to find. In particular, I have few doubts on time course. I have attached the matlab code and I hope that someone wants to help me. Thanks a lot.

epsilon = 0.62; % Voidage fraction

Kc = 4.5542e-5; % Mass Transfer Coefficient

Kl = 2.37; % Langmuir parameter

qs = 0.27; % Saturation capacity

cFeed = 10; % Feed concentration

L = 0.01; % Column length

R = 4.52e-10;

rho = 400;

t0 = 0; % Initial Time

tf = 900; % Final time

dt = 10; % Time step

z = [0:0.01:L]; % Mesh generation

t = [t0:dt:tf];% Time vector

n = numel(z); % Size of mesh grid

%Initial Conditions / Vector Creation

c0 = zeros(n,1);% t = 0, c = 0

c0(1)=cFeed;

q0 = zeros(n,1); % t = 0, q = 0 for all z,

q0(1) = zeros;

y0 = [c0 ; q0]; % Appends conditions together

%ODE15S Solver

[T, Y] = ode15s(@(t,y) Myfun(t,y,z,n),[t0 tf],y0);

plot(Y(:,n)/cFeed,'b')

xlabel('time')

ylabel('exit conc.')

function DyDt=Myfun(t,y,z,n)

% Defining Constants

D = 6.25e-4; % Axial Dispersion coefficient

v = 1*10^-3; % Superficial velocity

epsilon = 0.62; % Voidage fraction

Kc = 4.5542e-5; % Mass Transfer Coefficient

Kl = 2.37; % Langmuir parameter

qs = 0.27;% Saturation capacity

R = 4.52e-10;

rho = 400;

%Tc = zeros(n,1);

%q = zeros(n,1);

DcDt = zeros(n,1);

DqDt = zeros(n,1);

DyDt = zeros(2*n,1);

zhalf = zeros(n-1,1);

DcDz = zeros(n,1);

D2cDz2 = zeros(n,1);

c = y(1:n);

q = y(n+1:2*n);

% Interior mesh points

zhalf(1:n-1)=(z(1:n-1)+z(2:n))/2;

for i=2:n-1

DcDz(i) = ((z(i)-z(i-1))/(z(i+1)-z(i))*(c(i+1)-c(i))+(z(i+1)-z(i))/(z(i)-z(i-1))*(c(i)-c(i-1)))/(z(i+1)-z(i-1));

D2cDz2(i) = (zhalf(i)*(c(i+1)-c(i))/(z(i+1)-z(i))-zhalf(i-1)*(c(i)-c(i-1))/(z(i)-z(i-1)))/(zhalf(i)-zhalf(i-1));

end

% Calculate DcDz and D2cDz2 at z=L for boundary condition dc/dz = 0

DcDz(n) = 0;

D2cDz2(n) = -1.0/(z(n)-zhalf(n-1))*(c(n)-c(n-1))/(z(n)-z(n-1));

% Set time derivatives at z=0

% DcDt = 0 since c=cFeed remains constant over time and is already set as initial condition

% Standard setting for q

DqDt(1) = 3/R*Kc*( ((qs*c(1))/(Kl + c(1))) - q(1) );

DcDt(1) = 0.0;

% Set time derivatives in remaining points

for i=2:n

%Equation: dq/dt = k(q*-q) where q* = qs * (b*c)/(1+b*c)

DqDt(i) = 3/R*Kc*( ((qs*c(i))/(Kl + c(i))) - q(i) );

%Equation: dc/dt = D * d2c/dz2 - v*dc/dz - ((1-e)/(e))*dq/dt

DcDt(i) = D/epsilon*D2cDz2(i) - v/epsilon*DcDz(i) - ((1-epsilon)/(epsilon))*rho*DqDt(i);

end

% Concatenate vector of time derivatives

DyDt = [DcDt;DqDt];

end

[shubham chauhan]

Can u tell me ur model equations and plzz tell units of parameters u have taken

Answer 1

[Alan Stevens]

Do you just need a finer mesh? See:

epsilon = 0.62;  % Voidage fraction
Kc = 4.5542e-5;    % Mass Transfer Coefficient
Kl = 2.37;  % Langmuir parameter
qs = 0.27;    % Saturation capacity
cFeed = 10;     % Feed concentration
L = 0.01;         % Column length
R = 4.52e-10; 
rho = 400;
t0 = 0;        % Initial Time
tf = 90;      % Final time                             %%%%%%%%%%% Shorter time
dt = 1;      % Time step                               %%%%%%%%%%% Smaller timestep
z = 0:0.001:L; % Mesh generation                       %%%%%%%%%%% Finer mesh
t = t0:dt:tf;% Time vector
n = numel(z);  % Size of mesh grid
%Initial Conditions / Vector Creation
c0 = zeros(n,1);% t = 0, c = 0 
c0(1)=cFeed;
q0 = zeros(n,1); % t = 0, q = 0 for all z,
q0(1) = zeros; 
y0 = [c0 ; q0];     % Appends conditions together
%ODE15S Solver
[T, Y] = ode15s(@(t,y) Myfun(t,y,z,n),[t0 tf],y0);
plot(T, Y(:,n)/cFeed,'b')
xlabel('time')
ylabel('exit conc.')
function DyDt=Myfun(~,y,z,n)
% Defining Constants
D = 6.25e-4;    % Axial Dispersion coefficient
v = 1*10^-3;    % Superficial velocity
epsilon = 0.62;  % Voidage fraction
Kc = 4.5542e-5;    % Mass Transfer Coefficient
Kl = 2.37;  % Langmuir parameter
qs = 0.27;% Saturation capacity
R = 4.52e-10;
rho = 400;
%Tc = zeros(n,1);
%q = zeros(n,1);
DcDt = zeros(n,1);
DqDt = zeros(n,1);
%DyDt = zeros(2*n,1);
zhalf = zeros(n-1,1);
DcDz = zeros(n,1);
D2cDz2 = zeros(n,1);
c = y(1:n);
q = y(n+1:2*n);
% Interior mesh points
zhalf(1:n-1)=(z(1:n-1)+z(2:n))/2;
for i=2:n-1
  DcDz(i) = ((z(i)-z(i-1))/(z(i+1)-z(i))*(c(i+1)-c(i))+(z(i+1)-z(i))/(z(i)-z(i-1))*(c(i)-c(i-1)))/(z(i+1)-z(i-1));
  D2cDz2(i) = (zhalf(i)*(c(i+1)-c(i))/(z(i+1)-z(i))-zhalf(i-1)*(c(i)-c(i-1))/(z(i)-z(i-1)))/(zhalf(i)-zhalf(i-1));
end
% Calculate DcDz and D2cDz2 at z=L for boundary condition dc/dz = 0
DcDz(n) = 0; 
D2cDz2(n) = -1.0/(z(n)-zhalf(n-1))*(c(n)-c(n-1))/(z(n)-z(n-1));
% Set time derivatives at z=0
% DcDt = 0 since c=cFeed remains constant over time and is already set as initial condition
% Standard setting for q 
DqDt(1) = 3/R*Kc*( ((qs*c(1))/(Kl + c(1))) - q(1) );
DcDt(1) = 0.0;
% Set time derivatives in remaining points
for i=2:n
    %Equation: dq/dt = k(q*-q) where q* = qs * (b*c)/(1+b*c)
    DqDt(i) = 3/R*Kc*( ((qs*c(i))/(Kl + c(i))) - q(i) );
      %Equation: dc/dt = D * d2c/dz2 - v*dc/dz - ((1-e)/(e))*dq/dt
      DcDt(i) = D/epsilon*D2cDz2(i) - v/epsilon*DcDz(i) - ((1-epsilon)/(epsilon))*rho*DqDt(i);
  
end
% Concatenate vector of time derivatives
DyDt = [DcDt;DqDt];
end 

[Francesco Caputo]

Thankssssssssssssssssssssssss

[Francesco Caputo]

Hi Alan, I have another question...Why, if I change the parameters (qs or Kl, for example),does the breakthrough curve not change? Thank you.

[Alan Stevens]

Why not try changing them and see!

[Francesco Caputo]

I'm trying, but the curve doesn't change...

[Alan Stevens]

Ok, I'll have a look, though I'm not familiar with the physics here!

[Alan Stevens]

It does change! Your problem must be that you are changing the values at the top of the program, but these aren't used! You have redefined them within Myfun. These are the ones you need to change.

[Francesco Caputo]

It's my first time with matlab, I have to improve my knowledge. Thank you very much. I have really appreciated your help.

[Francesco Caputo]

Hi Alan, I need to visualize the area under the graph (which represents the adsorption efficiency), but I have some problems with the command "trapz". Could you help me? Thank you

[Alan Stevens]

Because you wrote "visualize" I assume you want to plot area under the curve as a function of time. If so, use cumtrapz (the cumulative area). Use the following in your code:

figure
subplot(2,1,1)
plot(T, Y(:,n)/cFeed,'b')
grid
xlabel('time')
ylabel('exit conc.')
title(['Kl = ', num2str(Kl), ' qs = ', num2str(qs)])

% Cumulative area

A = cumtrapz(T,Y(:,n)/cFeed);
subplot(2,1,2)
plot(T,A),grid
xlabel('time'), ylabel('Area')

This results in something like the following:

[Francesco Caputo]

Thanks very very much! Your help is foundamental for me!

[Francesco Caputo]

Dear Alan, I have one more doubt....I am not sure if I have well implemented the boundary conditions...Could you have a lot? I have attached the matlab code and the BCs. Thanks for your help.

epsilon = 0.62; % Voidage fraction

Ko = 17.3075e-6; % Mass Transfer Coefficient

Kl = 0.27; % Langmuir parameter

qs = 0.038; % Saturation capacity

cFeed = 10; % Feed concentration

L = 0.01; % Column length

R = 0.0885/100;

rho = 52500;

t0 = 0; % Initial Time

tf = 500; % Final time %%%%%%%%%%% Shorter time

dt = 1; % Time step %%%%%%%%%%% Smaller timestep

z = [0:0.001:L]; % Mesh generation %%%%%%%%%%% Finer mesh

t = [t0:dt:tf];% Time vector

n = numel(z); % Size of mesh grid

%Initial Conditions / Vector Creation

c0 = zeros(n,1);% t = 0, c = 0

c0(1)=cFeed;

q0 = zeros(n,1); % t = 0, q = 0 for all z,

q0(1) = zeros;

y0 = [c0 ; q0]; % Appends conditions together

%ODE15S Solver

[T, Y] = ode15s(@(t,y) Myfun(t,y,z,n),[t0 tf],y0);

plot(T, Y(:,n)/cFeed,'b')

xlabel('time [min]')

ylabel('Effluent concentration [mg/dl]')

title('Breakthrough curve')

Q = trapz(T,Y(:,n));

Efficiency_ads = Q/(tf*cFeed);

function DyDt=Myfun(~,y,z,n)

% Defining Constants

D = 6.25e-4; % Axial Dispersion coefficient

v = 0.01; % Superficial velocity

epsilon = 0.62; % Voidage fraction

Ko = 17.3075e-6; % Mass Transfer Coefficient

Kl = 0.027; % Langmuir parameter

qs = 0.038;% Saturation capacity

R = 0.0885/100;

rho = 52500;

%q = zeros(n,1);

DcDt = zeros(n,1);

DqDt = zeros(n,1);

%DyDt = zeros(2*n,1);

zhalf = zeros(n-1,1);

DcDz = zeros(n,1);

D2cDz2 = zeros(n,1);

c = y(1:n);

q = y(n+1:2*n);

% Interior mesh points

zhalf(1:n-1)=(z(1:n-1)+z(2:n))/2;

for i=2:n-1

DcDz(i) = ((z(i)-z(i-1))/(z(i+1)-z(i))*(c(i+1)-c(i))+(z(i+1)-z(i))/(z(i)-z(i-1))*(c(i)-c(i-1)))/(z(i+1)-z(i-1));

D2cDz2(i) = (zhalf(i)*(c(i+1)-c(i))/(z(i+1)-z(i))-zhalf(i-1)*(c(i)-c(i-1))/(z(i)-z(i-1)))/(zhalf(i)-zhalf(i-1));

end

% Calculate DcDz and D2cDz2 at z=L for boundary condition dc/dz = 0

DcDz(n) = 0;

D2cDz2(n) = -1.0/(z(n)-zhalf(n-1))*(c(n)-c(n-1))/(z(n)-z(n-1));

% Set time derivatives at z=0

% DcDt = 0 since c=cFeed remains constant over time and is already set as initial condition

% Standard setting for q

DqDt(1) = 3/R*Ko*( ((qs*c(1))/(Kl + c(1))) - q(1) );

DcDt(1) = 0.0;

% Set time derivatives in remaining points

for i=2:n

DqDt(i) = 3/R*Ko*( ((qs*c(i))/(Kl + c(i))) - q(i) );

DcDt(i) = D/epsilon*D2cDz2(i) - v/epsilon*DcDz(i) - ((1-epsilon)/(epsilon))*rho*DqDt(i);

end

% Concatenate vector of time derivatives

DyDt = [DcDt;DqDt];

end

[Alan Stevens]

I can't really tell. Possibly include

DcDz(1) = v*(c(1)-cFeed)/D;

though this doesn't seem to make any noticeable difference.

[Francesco Caputo]

Perfect, thank you very much!

[shubham chauhan]

Mr Alan steavns the code u solve for fixed bed absorption I need modelling equations and the parameters values

[shubham chauhan]

Mr Alan Stevens also tell which mathematical technique you used to solve the code for partial differential equation

[Alan Stevens]

@ shubham chauhan All the coding, values etc are given in the replies above.

[shubham chauhan]

I Need breakthrough curve at the exit of the bed help me developing code for this

[Alan Stevens]

@shubham chauhan The coding listed by Francesco on 10 Sep 2020 above, together with my reply just below it produces the breakthrough curve (whatever that might be!).

[shubham chauhan]

@alan stevens i have different modelling equations to obtain breakthrough curve. If u are interested in making code for that i will share my modelling equations. I have two PDEs to solve

[Aïli Bourbah]

@shubham chauhan Hi. I am a student working for a company about tetrachloro-ethene adsorption and desorption into activated carbon as a school project. Could I have your email to ask you some things about making a code for simulation ?

[Aniruddha Jain]

Can you send the research paper from where you solved the above equations?

[Anshul Suri]

Hi can someone please refer me to theoretical background over this problem?

[Federico Urbinati]

hi does anyone have the matlab code to solve adsorption when you have multiple components and not just one?

I hope someone can answer, I can not go forward. thank you

[Federico Urbinati]

Hi @Alan Stevens can you help me?