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/gn’, mean(q_avg));
% fprintf(‘Max q in column : %.4f mg/gn’, max(q_avg));
% fprintf(‘Theoretical qmax : %.4f mg/gn’, 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(:)];
endHello 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.