fixed bed adsorption column model-solving PDE-Langmiur isotherm
4 views (last 30 days)
Show older comments
Dear all,
I am trying to get a breakthrogh curve from the equations(attached) but am not getting proper plot. I need to model a fixed bed to obtain the breakthrough curves using langmiur isotherm. I am having problem with this code but dont know where is the problem. I would appreciate any help. Thank you very much.
Cfeed = 10; % Inlet concentration
tf = 86400; % End time
nz = 200; % Number of gateways
np = nz + 1; % Number of gateposts
% Initial Concentrations
C = zeros(np,1); C(1) = Cfeed;
Cs = zeros(np,1);
% Prepare data to pass to ode solver
c0 = [C; Cs];
dt = tf/10000;
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('Cs')
% title('Figure 2: Exit Cs vs time')
% figure
% plot(t, c(:,np)-c(:,2*np)), grid
% xlabel('time'), ylabel('Cb - Cs')
% title('Figure 3: Exit Cb-Cs vs time')
%
% k = 2.5*10^-5; nf = 1.45;
% q = k*c(:,2*np).^(1/nf); % 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)
Dz = 3*10^-8; % Diffusion coefficient
u = 2*10^-5; % Velocity
e = 0.4;
Kf = 3*10^-5;
ap = 2*10^-3; % Particle diameter
ep = 0.53;
L = 0.1; % Length
nz = 200; % Number of gateways
np = nz + 1; % Number of gateposts
dz = L/nz;
dr = ap/nz;
rho = 400;
Dp = 1.2*10^-10;
b=0.03;
qm=118;
C = c(1:np);
Cs = c(np+1:end);
Cs(1) = 0;
dCdt = zeros(np,1);
dCsdt = zeros(np,1);
dCsdr = zeros(np,1);
% rhalf = zeros(np-1,1);
% DCsDr = zeros(np,1);
D2CsDr2 = zeros(np,1);
%
% rhalf(1:np-1)=(z(1:np-1)+z(2:np))/2;
for i = 2:np-1
dCsdr(i) = (2*dr)*(Cs(i+1)-Cs(i-1));
D2CsDr2(i) = (1/dr^2)*(Cs(i+1)-2*Cs(i)+Cs(i-1));
% dCsdt(i) = 3*Kf/(ep*ap)*(C(i)-Cs(i)); % Needed to change the sign from -ve to +ve here.
end % i loop
for i = 2:np-1
dCdt(i) = (Dz/dz^2)*(C(i+1)-2*C(i)+C(i-1)) - u/(2*dz)*(C(i+1)-C(i-1)) - (1-e)*3*Kf/(e*ap)*(C(i)-Cs(i));
% dCsdt(i) = 3*Kf/(ep*ap)*(C(i)-Cs(i)); % Needed to change the sign from -ve to +ve here.
dCsdt(i) = 1/(1+rho*((1-ep)/ep)*(qm*b/((1+b*Cs(i))^2)))*Dp*(D2CsDr2(i)+((2/ap)*dCsdr(i)));
end % i loop
% For simplicity, set exit rates to those of the immediately previous nodes
dCdt(np) = dCdt(np-1);
dCsdt(np) = dCsdt(np-1);
% Combine rates for the function to return to the calling routine
dcdt = [dCdt; dCsdt];
end % of rates function
0 Comments
Answers (0)
See Also
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!