Fixed bed adsorption using pore diffusion model

Hello!

As the title suggests. Almost all of my code is from this thread, specifically given by one Alan Stevens (Freundlich4.m): fixed bed adsorption column model-solving PDE-freundlish isotherm - MATLAB Answers - MATLAB Central (mathworks.com)

I am trying to use the pore diffusion model which is given by:

(equation 1)

where q is the amount adsorbed, Ds is the diffusion coefficient of adsorbate in the solid phase, r is the particle radius coordinate. This model is usually applied on the an adsorbent particle.

The initial condition is:

and the boundary conditions are:

where, rp is the particle radius, kf is the external mass transfer coefficient, rho_p is the particle density, Ct is the concentration in the bulk fluid at time t, same as "C" in the code in the forum linked above and Cs is the concentration of adsorbate on the surface of the particle.

The dqdt is replaces the dqdt in the following equation:

(equation 2)

and has the same boundary and initial conditions as given in the link.

I modified the code but I am confused. Equation 1 and 2 have to be solved simultaneosuly, but this introduces an additional spatial parameter r. So one of the ways this maybe possible is by averaging the values of q at all r at a particular time instant to calculate dqdt and then replacing it in the equations 2 to calculate dCdt. However, this requires equation 1 to be fully solved before starting the calculations for dCdt.

I need help implementing this code. I could really use some.

Answer 1

[Torsten]

The "dq/dt" in equation 2 is the change of mass of q in the particle with respect to time. If you set

q_in_particle = integral_{r=0}^{r=r_p} 4*pi*r^2*q dr

and differentiate q_in_particle with respect to time, you arrive at

d(q_in_particle)/dt = 4* pi*r_p^2*D*(dq/dr)_{r=r_p}.

and (dq/dr)_{r=r_p} comes from the second boundary condition in your PDE for q.

[James99]

Excellent. I'll try that. Thank you :)

[James99]

Torsten, where does the q_in_particle equation come from? is it the same as

source: https://doi.org/10.1021/acs.iecr.0c05032 (equation 12)

[Torsten]

No, this is the mean concentration in the particle:

q_in_particle_mean = 1/(4/3*pi*r_p^3) * integral_{r=0}^{r=r_p} 4*pi*r^2 * q dr.

Thus your q_t from above equals q_in_particle/volume of particle.

[James99]

Noted. Thank you!

[James99]

Hi Torsten,

Just an additional question. How do I fit dq/dt with my own q vs t data? I am stuck on the boundary condition at r = rp since I have neither Ct nor Cs. I suppose they are calculated using eq. 2.

[Torsten]

dq/dt ~ gradient(q)/gradient(t) 

where q and t are your measurement vectors ?

[James99]

I thought about that but wouldn't it be better to use lsqnonlin? My objective is to find the diffusion coefficient

[Torsten]

Ct and Cs follow during the simulation of your equations. But you will never be able to estimate Ds because the mass transfer coefficient k_f is also unknown.

[James99]

I have already calculated k_f, or atleast a good approximation

[Torsten]

Ok, if you think that Ds is the only unknown parameter in your problem, you can start fitting it using lsqcurvefit, e.g..

[James99]

Hello Torsten,

Thank you for your previous reply.

I wanted to try modifying the original code to try and work with the original model (equation 1) and it seems to work. Although I get a curve I want for q vs t graph with just the current parameters, I would like to see a similar curve for Cb/Cin vs t graph. I've tried changing the Dz but to no avail. I am still in process of finding Dz but I just wanted to see how things change with changing parameters. Basically, Cb doesn't achieve equilibrium gradually but rather instantly. Do you mind taking a look at my code?

Thank you

Cfeed = 43.7;     % Inlet concentration ug/L

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

dp = 0.8; % mm
rp = dp/2;
r = linspace(0, rp, nz);
% q = zeros(np,1);

% 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/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('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')

% qm = 23660;
% KL = 0.00051;
k = 2.5*10^-5; nf = 1.45;
q = k*c(:,2*np).^(1/nf); % Freundlich equation 
% q = (qm*KL.*c(:,2*np))/(1 + KL.*c(:,2*np)); % Langmuir equation
figure
plot(t,q), grid
xlabel('time'), ylabel('q')
title('Figure 4: Exit q vs time')






function dcdt = rates(~, c)
Dz = 3*10^-10;   % axial Diffusion coefficient
Ds = 3*10^-8;   % surface Diffusion coefficient
u = 0.036;      % Velocity m/s 
% u = 2*10^-4;    % Velocity
e = 0.71;       % Bed porosity
Kf = 2.973E-06;
rho_p = 2100;   %kg/m3
% ap = 1.1/1000;  % m
ap = 0.0011;   % Particle diameter m
ep = 0.71;       % Particle porosity
L = 0.04;       % Bed length in m
nz = 200;       % Number of gateways
np = nz + 1;    % Number of gateposts
dz = L/nz;
dr = ap/nz;
r = linspace(0, ap, nz);

C = c(1:np);
Cs = c(np+1:end);
dCdt = zeros(np,1);
dCsdt = zeros(np,1);
dqdt = zeros(np,1);
q = zeros(np,1);

for i = 2:np-1
    q(np) = (-2*ap*Kf*(C(i)-Cs(i)))/(rho_p*Ds) + q(np-1);
    q(1) = q(2);
    dqdr = q(i+1)/2/dr - q(i-1)/2/dr;
    d2qdr2 = q(i+1)/dr^2 - 2*q(i)/dr^2 + q(i-1)/dr^2;
    dqdt(i) = Ds/ap^2*(2*r(i)*dqdr + r(i)*r(i)*d2qdr2);

    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)*dqdt(i);
    dCsdt(i) = 3*Kf/(ep*ap)*(C(i)-Cs(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
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

[Torsten]

I'm confused. What do you think the program that you posted solves ?

[James99]

The only modification I made from the original code is the dq/dt equation. The program still solves for C, Cs and q. My apologies if I am not being clear.

[Torsten]

Sorry, but this code is completely wrong.

If you want to solve equation (1) and equation (2), you have to solve nz*(np+1) differential equations where nz is the number of discretization points in axial direction and np is the number of discretization points "into the particle" in each axial discretization point. The extra +1 in (np+1) is for the concentration of the adsorbens in the gas phase. You should imagine that in each cell in axial direction of the reactor, there exists a "typical" state of particle loading which is simulated by solving equation (1). The concentration of the adsorbat in the particle (equation (1)) and the concentration of the adsorbens in the gas phase (equation (2)) are coupled by the mass transfer term at the particle surface in each discretization point in axial direction.

I only know about this book in German

https://onlinelibrary.wiley.com/doi/abs/10.1002/cite.330530304

but you should really first read about the physical background before you continue coding.

[James99]

Understood, thank you!

[James99]

Hello Torsten,

I thought about what you said and I made some modifications to the code. The code I attached is incomplete. I can't get my head around how I should relate Cs between the two equations and generally just implement that into the code. Do you mind taking a look at the code?

[Torsten]

As far as I remember, Cs is the equilibrium concentration corresponding to the concentration in the bulk fluid. You must supply Cs as a function of C and maybe of temperature. So I don't know how Csdt comes into play: C (bulk concentration of the adsorbens in the gas phase) and q (adsorbat concentration in the particle) are the variables to be solved for.

[James99]

Hello Torsten,

I am still trying to develop the model however I have a couple of queries. Firstly, let me clarify the equations

Fluid phase mass balance:

The pore diffusion model is still the same with a slight change in the second boundary condition

I have already written the code but my worry is that the code isn't exactly sharing the variable Cs among the two equations very well. And I am not quite sure where exactly I should put the above mentioned boundary condition. I put it outside the loop done for the calculations for particle.

Do you mind taking a look at the code and giving me some feedback? I really want to try and make this work. The code doesn't exactly give any errors but it does give warnings about integration tolerances and whatever it does calculate before the warning is quite abnormal. Furthermore, the parameters (Dp, Dz and kf) may not be right, but nevertheless, they are enough to give me a breakthough curve.

Thank you!

(The equations are taken from this source https://doi.org/10.1016/j.psep.2018.04.027 should you wish to see it. The fluid phase equation above isn't exactly the same as the paper but that shouldn't be a problem)

Answer 2

[James99]

For anyone interested, the attached code kinda works, or atleast serve you as a starting point. But make sure that you have your diffusion coefficients right because mine aren't.