Help with HSDM model for lithium adsorption: simulated curves do not match experimental data (based on Jiang et al., 2020)

Hello everyone,

I’m trying to replicate the results from the paper by Jiang et al. (2020):

*“Application of concentration-dependent HSDM to the lithium adsorption from brine in fixed bed columns” – Separation and Purification Technology 241, 116682.

The paper models lithium adsorption in a fixed bed packed with Li–Al layered double hydroxide resins. It uses a concentration-dependent Homogeneous Surface Diffusion Model (HSDM), accounting for axial dispersion, film diffusion, surface diffusion, and Langmuir equilibrium. I've implemented a full MATLAB simulation including radial discretization inside the particles and a coupling with the axial profile.

I’ve followed the theoretical development very closely and used the same parameters as those reported by the authors. However, the breakthrough curves generated by my model still don’t fully match the experimental data shown in the paper (especially at intermediate values of C_out/C_in).

I suspect there may be a mistake in my implementation of either:

• the mass balance in the fluid phase,
• the boundary condition at the surface of the particle,
• or how I define the rate of mass transfer using Kf.

I'm sharing my complete code below, and I would appreciate any suggestions or corrections you may have.

function hsdm_column_simulation_v2
% Full HSDM model with radial diffusion dependent on loading (Jiang et al. 2023)
clc; clear; close all

%% ==== PARÁMETROS DEL SISTEMA ====   -->   %% ==== SYSTEM PARAMETERS ====
% Parámetros de la columna              -->   % Column parameters
Q = (15e-6)/60;            % Brine flow rate [m3/s]
D = 0.02;                  % Column diameter [m]
A = pi/4 * D^2;            % Cross-sectional area [m2]
v = Q/A;                   % Superficial or interstitial velocity [m/s]
epsilon = 0.355;           % Bed void fraction
mu = 6.493e-3;             % Fluid dynamic viscosity [Pa·s] 

% Parámetros de la resina               -->   % Resin properties
rho = 1.3787;              % Solid density [g/L] (coherent with q in mg/L)
dp = 0.65/1000;            % Particle diameter [m]
Dm = 1.166e-5 / 10000;     % Lithium diffusivity [m2/s]
R = 0.000325;              % Particle radius [m]
Dax = 0.44*Dm + 0.83*v*dp; % Axial dispersion [m²/s] = 4.2983e-5

%% Isoterma de Langmuir                -->   %% Langmuir Isotherm
qmax = 5.9522;             % Langmuir parameter [mg/g]   
b = 0.03439;               % Langmuir parameter [L/mg]

%% Difusión superficial dependiente de la carga   -->   %% Surface diffusion dependent on loading
Ds0 = 4e-14 ;              % Surface diffusion coef. at 0 coverage [m²/s] (original was 3.2258e-14)
k_exp = 0.505;             % Dimensionless constant

%% Correlaciones empíricas            -->   %% Empirical correlations
Re = (rho * v * dp) / mu;
Sc = mu / (rho * Dm);
Sh = 1.09 + 0.5 * Re^0.5 * Sc^(1/3);
Kf = Sh * Dm / dp;         % 6.5386e-5

%% Discretización                    -->   %% Discretization
L = 0.60;                  % Column height [m]
Nz = 20;                   % Axial nodes
Nr = 5;                    % Radial nodes per particle
dz = L / (Nz - 1);
dr = R / Nr;

%% Condiciones operacionales         -->   %% Operating conditions
cFeedVec = [300, 350, 400];  % mg/L
tf = 36000;                  % Final time [s] (600 min)
tspan = [0 tf];
colores = ['b','g','r'];

%% Figura                           -->   %% Plot
figure; hold on
for j = 1:length(cFeedVec)
    cFeed = cFeedVec(j);

    % Condiciones iniciales para el lecho y la partícula
    c0 = zeros(Nz,1);     % Initial concentration in fluid: C = 0
    q0 = zeros(Nz*Nr,1);  % Initial loading in particles: q = 0
    y0 = [c0; q0];

    % Agrupación de parámetros       -->   % Parameter grouping
    param = struct('Nz',Nz,'Nr',Nr,'dz',dz,'dr',dr,'R',R,...
                   'epsilon',epsilon,'rho',rho,'v',v,'Dax',Dax,...
                   'qmax',qmax,'b',b,'Ds0',Ds0,'k',k_exp,...
                   'cFeed',cFeed,'Kf',Kf);

    %% Resolución del sistema con ode15s   -->   %% System solution using ode15s
    [T, Y] = ode15s(@(t,y) hsdm_rhs(t,y,param), tspan, y0);

    % Ploteo de gráfico con salida normalizada  -->  % Plot normalized outlet
    C_out = Y(:,Nz);
    plot(T/60, C_out / cFeed, colores(j), 'LineWidth', 2, ...
        'DisplayName', ['C_{in} = ' num2str(cFeed) ' mg/L']);
end

% %% Evaluación de carga adsorbida para corroborar modelo
% Nz = param.Nz;
% Nr = param.Nr;
% R = param.R;
% dr = param.dr;
% qmax = param.qmax;
% q_final = reshape(Y(end, Nz+1:end), Nr, Nz);  % q(r,z) at final time
% q_avg = trapz(linspace(0, R, Nr), q_final .* linspace(0, R, Nr)', 1) * 2 / R^2;
% fprintf('\n----- Saturation Analysis for C_in = %d mg/L -----\n', cFeed);
% fprintf('Global average q     : %.4f mg/g\n', mean(q_avg));
% fprintf('Max q in column      : %.4f mg/g\n', max(q_avg));
% fprintf('Theoretical qmax     : %.4f mg/g\n', qmax);

%% Gráfico de Cout/Cin vs tiempo    -->   %% Plot Cout/Cin vs time
xlabel('Tiempo (min)')
ylabel('C_{out} / C_{in}')
title('Curva de ruptura Modelo HSDM')
legend('Location','southeast')
set(gca, 'FontName', 'Palatino Linotype')  % Axes font
box on

%% Puntos experimentales aproximados (visualmente desde el paper)
t_exp = 0:30:600;
Cexp_300 = [0.00 0.22 0.36 0.48 0.57 0.63 0.69 0.73 0.77 0.80 ...
            0.82 0.84 0.85 0.86 0.87 0.88 0.88 0.89 0.89 0.89 0.90];
Cexp_350 = [0.00 0.26 0.40 0.53 0.62 0.69 0.74 0.78 0.81 0.83 ...
            0.85 0.86 0.87 0.88 0.88 0.89 0.89 0.90 0.90 0.90 0.91];
Cexp_400 = [0.00 0.29 0.45 0.59 0.68 0.75 0.79 0.83 0.85 0.87 ...
            0.88 0.89 0.90 0.90 0.91 0.91 0.91 0.91 0.92 0.92 0.92];

% Superponer puntos experimentales    -->   % Plot experimental data
plot(t_exp, Cexp_400, 'r^', 'MarkerFaceColor', 'r', 'DisplayName', 'Exp 400 mg/L')
plot(t_exp, Cexp_350, 'go', 'MarkerFaceColor', 'g', 'DisplayName', 'Exp 350 mg/L') 
plot(t_exp, Cexp_300, 'bs', 'MarkerFaceColor', 'b', 'DisplayName', 'Exp 300 mg/L') 

end

%% FUNCIÓN RHS DEL MODELO HSDM    -->    %% RHS FUNCTION FOR HSDM MODEL
function dydt = hsdm_rhs(~, y, p)
% Extraer parámetros           -->   % Extract parameters
Nz = p.Nz; Nr = p.Nr; dz = p.dz; dr = p.dr;
R = p.R; epsilon = p.epsilon; rho = p.rho;
v = p.v; Dax = p.Dax; qmax = p.qmax; b = p.b;
Ds0 = p.Ds0; k_exp = p.k; cFeed = p.cFeed; Kf = p.Kf;

% Separar variables            -->   % Split variables
c = y(1:Nz);
q = reshape(y(Nz+1:end), Nr, Nz);  % q(r,z)

% Inicializar derivadas        -->   % Initialize derivatives
dc_dt = zeros(Nz,1);
dq_dt = zeros(Nr,Nz);

%% Balance de masa axial (fase fluida)  -->   %% Mass balance (fluid phase)
for i = 2:Nz-1
    dcdz = (c(i+1)-c(i-1))/(2*dz);       
    d2cdz2 = (c(i+1) - 2*c(i) + c(i-1))/dz^2;
    qsurf = q(end,i);                      
    Csurf = c(i);                          
    dqR_dt = (3/R) * Kf * (Csurf - qsurf);
    dc_dt(i) = Dax * d2cdz2 - v * dcdz - ((1 - epsilon)/epsilon) * rho * 1000 * dqR_dt; 
end

%% Condiciones de contorno en la columna
% Entrada z = 0 - tipo Danckwerts
qsurf = q(end,1);
Csurf = c(1);
dqR_dt_in = (3 / R) * Kf * (Csurf - qsurf);
dc_dt(1) = Dax * (c(2) - c(1)) / dz^2 - v * (cFeed - c(1)) - ((1 - epsilon)/epsilon) * rho * 1000 * dqR_dt_in;

% Salida z = L
dc_dt(end) = Dax * (c(end-1) - c(end)) / dz;

%% Difusión radial al interior de la partícula
for iz = 1:Nz
    for ir = 2:Nr-1
        rq = (ir-1)*dr;
        d2q = (q(ir+1,iz) - 2*q(ir,iz) + q(ir-1,iz)) / dr^2;
        Ds = Ds0 * (1 - q(ir,iz)/qmax)^k_exp;
        dq_dt(ir,iz) = Ds * (d2q + (2/rq)*(q(ir+1,iz)-q(ir-1,iz))/(2*dr));
    end

    % Centro de la partícula
    dq_dt(1,iz) = 0;

    % Superficie de la partícula
    q_eq = (qmax * b * c(iz)) / (1 + b * c(iz));
    dq_dt(Nr,iz) = 3 * Kf / R * (q_eq - q(Nr,iz));
end

% Vector final
dydt = [dc_dt; dq_dt(:)];
end

[Torsten]

dc_dt(1) = Dax * (c(2) - c(1)) / dz^2 - v * (cFeed - c(1)) - ((1 - epsilon)/epsilon) * rho * 1000 * dqR_dt_in;

How did you come up with this expression for dc/dt at z = 0 ?

If you plot the concentration curves from your simulation over the reactor height L for certain times, you will see that this setting leads to negative concentrations at the inlet (whereas it should be near to the feed concentration for larger times).

[Ignacia]

Hi Torsten, thank you for your comment, that's a very helpful observation. I initially thought the expression I used was a form of the Danckwerts boundary condition, but I now see that it might not be correctly applied. I'm still trying to understand how to implement the inlet condition properly, especially under the presence of axial dispersion and the coupling with the sink term due to particle uptake.

Would you suggest applying the Danckwerts boundary by directly discretizing the first derivative at the inlet? Or is there a preferred way to formulate it when surface diffusion is also involved? Any clarification or suggestion would be much appreciated. Thanks again for pointing this out!

[Torsten]

Your equation for c reads

dc/dt = -v*dc/dz + Dax*d^2c/dz^2 - (1-eps)/eps*rhos*dq/dt = d/dz (-v*c+Dax*dc/dz) - (1-eps)/eps*rhos*dq/dt

At z = 0, you have the boundary condition

Dax*dc/dz = v*(c-cFeed)

Thus the term

d/dz (-v*c+Dax*dc/dz) at z = 0

becomes approximately

((-v*c+Dax*dc/dz) @z=dz/2 - (-v*c+Dax*dc/dz)@z=0)/(dz/2 - 0) =

((-v*c(1)+Dax*(c(2)-c(1))/dz) + v*cFeed) / (dz/2) =

2*(v*(cFeed-c(1))/dz + Dax*(c(2)-c(1))/dz^2)

Thus you get

dc_dt(1) = 2*(v*(cFeed-c(1))/dz + Dax*(c(2)-c(1))/dz^2) - ((1 - epsilon)/epsilon) * rho * 1000 * dqR_dt;

as time derivative for c at z = 0.

I suggest to discretize v*dc/dz at z = iz*dz (1<=iz<=Nz) as v*(c(iz)-c(iz-1))/dz, not as v*(c(iz+1)-c(iz-1))/(2*dz) as you do because of instabilities that might occur due to a central differences discretization.

Further,

% Centro de la partícula

dq_dt(1,iz) = 0;

is wrong: you have dq/dr = 0 at the center, not dq/dt.

Similarly,

% Salida z = L

dc_dt(end) = Dax * (c(end-1) - c(end)) / dz;

is wrong: you have dc/dz = 0 at z = L, not dc/dt.

And this is not boundary condition (13) from the article which is necessary to implement in order to establish mass conservation together with equation (9) from the article:

% Superficie de la partícula

q_eq = (qmax * b * c(iz)) / (1 + b * c(iz));

dq_dt(Nr,iz) = 3 * Kf / R * (q_eq - q(Nr,iz));

[Ignacia]

Hi Torsten,

Thank you very much for your helpful feedback. Based on your suggestion, I corrected the discretization for the convective and dispersive terms at the inlet, and updated the corresponding section of my code as follows:

%% Balance de masa axial (fase fluida)
for i = 2:Nz-1
    dcdz = (c(i) - c(i-1)) / dz;
    d2cdz2 = (c(i+1) - 2*c(i) + c(i-1)) / dz^2;

    dqR_dt = (3 / R) * Kf * (c(i) - q(end,i));
    dc_dt(i) = Dax * d2cdz2 - v * dcdz - ((1 - epsilon)/epsilon) * rho * 1000 * dqR_dt; 
end

%% Condición de contorno en la entrada z = 0
qsurf = q(end,1);
Csurf = c(1);
dqR_dt_in = (3 / R) * Kf * (Csurf - qsurf);
dc_dt(1) = 2*(v*(cFeed - c(1))/dz + Dax*(c(2) - c(1))/dz^2) - ((1 - epsilon)/epsilon) * rho * 1000 * dqR_dt_in;

Just to confirm — is this what you meant regarding the correction for the inlet boundary condition and convective term?

I’m still unsure how to correctly implement the outlet condition at z = L and the particle surface condition at r = R, particularly with regard to ensuring full mass conservation as per Eq. (13). Any further guidance would be very appreciated.

[Torsten]

Just to confirm — is this what you meant regarding the correction for the inlet boundary condition and convective term?

Yes.

I’m still unsure how to correctly implement the outlet condition at z = L and the particle surface condition at r = R, particularly with regard to ensuring full mass conservation as per Eq. (13). Any further guidance would be very appreciated.

The discretization of equation (10) together with boundary conditions (12) and (13) and of the outlet condition (8) at z = L is treated in Chapter (3) of the enclosed paper for the cases c = 2 (spherical coordinates) and c = 0 (cartesian coordinates). Come back if you need confirmation for your attempt to implement this in your code.

[Ignacia]

Hi Torsten,

Following your guidance regarding the boundary conditions at z=L and r=R, I reviewed Chapter 3 of the Sinkovec & Madsen document and implemented both in my MATLAB code as follows:

1. For the Boundary at z=L (Cartesian coordinates c=0):

Since the domain is in Cartesian coordinates at this boundary, I considered the general second-order derivative formulation (eq. 3.1c) with c=0. However, at z=Lz (last axial node), the centered scheme would require a ghost node cNz+1 which is not defined. As suggested in the document, I used a mirror node assumption cNz+1=cNz−1​, which leads to a second derivative approximated by:

​​

This simplifies the term to:

d2cdz2_L = (c(end-2) - 2*c(end-1) + c(end-2)) / dz^2;
dc_dt(end) = Dax * d2cdz2_L - ((1 - epsilon)/epsilon) * rho * 1000 * dqR_dt_end;

2. Boundary at r=R (Spherical coordinates c=2):

For the radial boundary at the particle surface, I implemented the boundary condition from eq. (13) in Jiang et al., which matches the spherical diffusion equation with surface flux continuity:

Since this is also a boundary, a backward difference was used for the radial derivative, as recommended in Sinkovec paper. The approximation is:

I directly equated this to Jiang’s expression (i think) and solved for the time derivative at the surface:

qR = q(Nr, iz);
qR_dr = q(Nr-1, iz);  

Ds_R = Ds0 * (1 - qR / qmax)^k_exp;  

dq_dt(Nr, iz) = Ds_R * ( (qR - qR_dr) / dr ) / R^2

About the condition in r=0 im not sure.

Let me know if these implementations correctly reflect the methods described, or if you’d suggest any refinements. And again, thanks for the help!

%% FUNCIÓN RHS DEL MODELO HSDM
function dydt = hsdm_rhs(~, y, p)
% Extraer parámetros
Nz = p.Nz; Nr = p.Nr; dz = p.dz; dr = p.dr;
R = p.R; epsilon = p.epsilon; rho = p.rho;
v = p.v; Dax = p.Dax; qmax = p.qmax; b = p.b;
Ds0 = p.Ds0; k_exp = p.k; cFeed = p.cFeed; Kf = p.Kf;


% Separar variables
c = y(1:Nz);
q = reshape(y(Nz+1:end), Nr, Nz);  % q(r,z)

% Inicializar derivadas
dc_dt = zeros(Nz,1);
dq_dt = zeros(Nr,Nz);

%% Balance de masa axial (fase fluida)
for i = 2:Nz-1
    dcdz = (c(i) - c(i-1)) / dz; %dcdz = (c(i+1)-c(i-1))/(2*dz);               % Prmera derivada convectiva, parte de la ecuación 5 (u dC/dx)
    d2cdz2 = (c(i+1) - 2*c(i) + c(i-1))/dz^2;    % Segunda derivada axial de C, parte de la ecuación 5 (Dax d2C/dx2)

    %% Flujo al sólido en la superficie de la partícula
    %qsurf = q(end,i);                            % Carga en la superficie de la partícula
    %Csurf = c(i);                                % Concentración en fluido en ese nodo
    dqR_dt = (3 / R) * Kf * (c(i) - q(end,i));       % Tasa de transferencia al sólido según kf

    %% Balance de masa en el fluido 
    dc_dt(i) = Dax * d2cdz2 - v * dcdz - ((1 - epsilon)/epsilon) * rho * 1000  * dqR_dt; 
end

%% Condiciones de contorno en la columna
%% Condicion de contorno en la entrada z = 0 - tipo Danckwerts
qsurf = q(end,1);
Csurf = c(1);
dqR_dt_in = (3 / R) * Kf * (Csurf - qsurf);
dc_dt(1) = 2*(v*(cFeed - c(1))/dz + Dax*(c(2) - c(1))/dz^2) - ((1 - epsilon)/epsilon) * rho * 1000 * dqR_dt_in;

%% Condiciones de contorno en la salida z = L
dqR_dt_end = (3 / R) * Kf * (c(end) - q(end,end));

% Nodo ficticio c(end+1) = c(end-1) se impone que la derivada espacial sea 0 en el extremo usando esquema de 3 puntos para garantizar la conservación (Sinkovec-Madsen)
d2cdz2_L = (c(end-2) - 2*c(end-1) + c(end-2)) / dz^2;  % en el borde derecho (z=L) no existe u_(Nz+1) por lo que se aplica nodo ficticio (Sinkovec) se asuma du/dz para z=L = 0
dc_dt(end) = Dax * d2cdz2_L - ((1 - epsilon)/epsilon) * rho * 1000 * dqR_dt_end;

%% Difusión radial al interior de la partícula 
for iz = 1:Nz
    for ir = 2:Nr-1
        rq = (ir-1)*dr;
        d2q = (q(ir+1,iz) - 2*q(ir,iz) + q(ir-1,iz)) / dr^2; % Ecuación de la segunda derivada NUMERICO
        Ds = Ds0 * (1 - q(ir,iz)/qmax)^k_exp;  % Coeficiente de difusión superficial
        dq_dt(ir,iz) = Ds * (d2q + (2/rq)*(q(ir+1,iz)-q(ir-1,iz))/(2*dr));  % Ecuacion de la difusión radial al interior de la partícula (Ds) primera derivada
    end

    %% Condición en r=0 (centro de la partícula)
    dq_dt(1, iz) = (q(2, iz) - q(1, iz)) / dr;  % debe tender a cero
    %dq_dt(1, iz) = dq_dt(2, iz);  % this or the other one???

    %% Condición en r=R (flujo hacia partícula = flujo desde el fluido)
    % basada en la ecuación (13) del artículo -Ds * (∂q/∂r)|r=R = Kf * (C - q_surf)
    %Superficie r=R usando diferencias finitas hacia atrás según Sinkovec-Madsen
    
    % Nodo superficie r = R
    qR = q(Nr, iz);
    qR_dr = q(Nr-1, iz);  % nodo anterior

    Ds_R = Ds0 * (1 - qR / qmax)^k_exp;  % Ds evaluado en R

    % Condición de contorno completa en r = R (derivada hacia atrás + flujo igualado)
    dq_dt(Nr, iz) = (Kf * (c(iz) - qR)) / rho;

    
end

% Vector final
dydt = [dc_dt; dq_dt(:)];

end

[Torsten]

dc_dt(end) = Dax * d2cdz2_L - ((1 - epsilon)/epsilon) * rho * 1000 * dqR_dt_end;

This cannot be correct because there is no convective outflow term. All the mass will accululate in the last cell.

And there is nothing in d2cdz2_L that implements dc/dz = 0 at z = L.

    dq_dt(1, iz) = (q(2, iz) - q(1, iz)) / dr;  % debe tender a cero
    % Condición de contorno completa en r = R (derivada hacia atrás + flujo igualado)
    dq_dt(Nr, iz) = (Kf * (c(iz) - qR)) / rho;

You made the same mistake again: The right-hand side is dq/dr, not dq/dt.

Take the time. Read and digest the schemes described in chapter (3) of the article.

I honestly don't know how you can start with such a complicated problem if you have no experience with the discretization and numerical solution of partial differential equations. During my studies I was always told to start simple and add complexity step by step to have control that the results remain senseful and to make it easy to recover potential mistakes in the implementation.

Here is the approximation for dc/dt at z = L. It's not taken from the article.

    dqR_dt_end = (3 / R) * Kf * (c(end) - q(end,end));
    dc_dt(end) = -v*(c(end)-c(end-1))/dz - 2*Dax * (c(end)-c(end-1))/dz^2 - ((1 - epsilon)/epsilon) * rho * 1000 * dqR_dt_end;        

But for the spherical case (c=2) (diffusive transport in the particle), you should take both the discretization in interior points and for the boundary conditions from the article because it's quite complicated to deduce. A backward difference approximation is used only for the derivative in a boundary point if there is no spatial derivative in the given boundary condition. This is not the case for your problem: at both sides (at r = 0 and at r = R), dc_p/dr is contained in the boundary conditions (equations (12) and (13)). Thus your statement

Since this is also a boundary, a backward difference was used for the radial derivative, as recommended in Sinkovec paper.

is not correct.

[Torsten]

Do you know how c_p and q are related ? I don't think your dqR_dt is correct because you subtract c and q which have different units.

Do you know where the Langmuir isotherm is used in the article's model ?