Fixed-bed adsorption

19 次查看（过去 30 天）
Francesco Caputo
Francesco Caputo 2020-9-5
评论： Federico Urbinati 2022-10-24
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
  1 个评论
shubham chauhan
shubham chauhan 2020-11-17
Can u tell me ur model equations and plzz tell units of parameters u have taken

请先登录，再进行评论。

请先登录，再回答此问题。

采纳的回答

Alan Stevens
Alan Stevens 2020-9-5
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
  24 个评论
显示 22更早的评论
Federico Urbinati
Federico Urbinati 2022-10-18
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

请先登录，再进行评论。

更多回答（0 个）

类别

在 Help Center 和 File Exchange 中查找有关 Mathematics 的更多信息

标签

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by