MATLAB: Fixed-bed adsorption

#fxed-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

Best Answer

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
Related Question