hello, everybody
I would like to make the rounded rectangle and I would like them to translated and rotated.
First I tried with rectangle function using ‘Curvature’ of 1. then I can make the rounded rectangle.
However, it is not possible for them to translated and rotated.
Therefore, I made the rectangle with polyshape function. and it is okay for them to translated and rotated.
However, i do not know how I can make them to have rounded edges.
Could you please helpe me how to make them to rounded edges with polyshape function?
clear; close all; clc;
figure(1); hold on; axis equal
% Before translating and rotating
rectangle(‘Position’,[-14.1276, 226.1976, 6.5929, 9.4184],’FaceColor’,’r’,’Curvature’,1)

% After translating and rotating
% translate and rotate polygon
x = [-14.1276; -14.1276+6.5929; -14.1276+6.5929; -14.1276];
y = [226.1976; 226.1976; 226.1976+9.4184; 226.1976+9.4184];
p = [x,y];
pgon = polyshape(p);
pgon = translate(pgon,50,-75);
pgon = rotate(pgon,18,[57, 349]);
plot(pgon,’FaceColor’,’red’,’FaceAlpha’,1)
% object by rectangle function can not translated and rotated.hello, everybody
I would like to make the rounded rectangle and I would like them to translated and rotated.
First I tried with rectangle function using ‘Curvature’ of 1. then I can make the rounded rectangle.
However, it is not possible for them to translated and rotated.
Therefore, I made the rectangle with polyshape function. and it is okay for them to translated and rotated.
However, i do not know how I can make them to have rounded edges.
Could you please helpe me how to make them to rounded edges with polyshape function?
clear; close all; clc;
figure(1); hold on; axis equal
% Before translating and rotating
rectangle(‘Position’,[-14.1276, 226.1976, 6.5929, 9.4184],’FaceColor’,’r’,’Curvature’,1)

% After translating and rotating
% translate and rotate polygon
x = [-14.1276; -14.1276+6.5929; -14.1276+6.5929; -14.1276];
y = [226.1976; 226.1976; 226.1976+9.4184; 226.1976+9.4184];
p = [x,y];
pgon = polyshape(p);
pgon = translate(pgon,50,-75);
pgon = rotate(pgon,18,[57, 349]);
plot(pgon,’FaceColor’,’red’,’FaceAlpha’,1)
% object by rectangle function can not translated and rotated. hello, everybody
I would like to make the rounded rectangle and I would like them to translated and rotated.
First I tried with rectangle function using ‘Curvature’ of 1. then I can make the rounded rectangle.
However, it is not possible for them to translated and rotated.
Therefore, I made the rectangle with polyshape function. and it is okay for them to translated and rotated.
However, i do not know how I can make them to have rounded edges.
Could you please helpe me how to make them to rounded edges with polyshape function?
clear; close all; clc;
figure(1); hold on; axis equal
% Before translating and rotating
rectangle(‘Position’,[-14.1276, 226.1976, 6.5929, 9.4184],’FaceColor’,’r’,’Curvature’,1)

% After translating and rotating
% translate and rotate polygon
x = [-14.1276; -14.1276+6.5929; -14.1276+6.5929; -14.1276];
y = [226.1976; 226.1976; 226.1976+9.4184; 226.1976+9.4184];
p = [x,y];
pgon = polyshape(p);
pgon = translate(pgon,50,-75);
pgon = rotate(pgon,18,[57, 349]);
plot(pgon,’FaceColor’,’red’,’FaceAlpha’,1)
% object by rectangle function can not translated and rotated. curvature MATLAB Answers — New Questions

​

THIS IS A FILE DESCRIBE THE EQUATION AND THE , the Lagrangian -derived equations of motionTHIS IS A FILE DESCRIBE THE EQUATION AND THE , the Lagrangian -derived equations of motion THIS IS A FILE DESCRIBE THE EQUATION AND THE , the Lagrangian -derived equations of motion matlab, mathematics, code MATLAB Answers — New Questions

​

I’d like to know where the mistake in this code, and the solution

h = 0.00172; rho = 1620;
L = 100*h; c=L*L/4;
M=10; N=10;
syms xi eta
shape_w=sym(zeros(M,N)); x_m=sym(zeros(M,1)); y_n=sym(zeros(N,1));

for m=1:M
for n=1:N
x_m(m)=(xi^m)*(1+xi);
y_n(n)=(eta^n)*(1-eta)*(1+eta)^2;
shape_w(m,n)=x_m(m)*y_n(n);
end
end
fun=(shape_w*shape_w’);
M_w =c*(rho*h)*integral2(@(xi,eta)fun,-1,1,-1,1);I’d like to know where the mistake in this code, and the solution

h = 0.00172; rho = 1620;
L = 100*h; c=L*L/4;
M=10; N=10;
syms xi eta
shape_w=sym(zeros(M,N)); x_m=sym(zeros(M,1)); y_n=sym(zeros(N,1));

for m=1:M
for n=1:N
x_m(m)=(xi^m)*(1+xi);
y_n(n)=(eta^n)*(1-eta)*(1+eta)^2;
shape_w(m,n)=x_m(m)*y_n(n);
end
end
fun=(shape_w*shape_w’);
M_w =c*(rho*h)*integral2(@(xi,eta)fun,-1,1,-1,1); I’d like to know where the mistake in this code, and the solution

h = 0.00172; rho = 1620;
L = 100*h; c=L*L/4;
M=10; N=10;
syms xi eta
shape_w=sym(zeros(M,N)); x_m=sym(zeros(M,1)); y_n=sym(zeros(N,1));

for m=1:M
for n=1:N
x_m(m)=(xi^m)*(1+xi);
y_n(n)=(eta^n)*(1-eta)*(1+eta)^2;
shape_w(m,n)=x_m(m)*y_n(n);
end
end
fun=(shape_w*shape_w’);
M_w =c*(rho*h)*integral2(@(xi,eta)fun,-1,1,-1,1); how i can fix error using integral2calc integral2t/tensor (line 231)? MATLAB Answers — New Questions

​

Hi everyone,
I created a neural network that uses a functionLayer that I defined by :
functionLayer(@Feature_wise_LM,Name="FiLM_64",Formattable=1,NumInputs=2,NumOutputs=1)
However, when I save the trained neural network with :
save("./trained_networks/FiLM", "trained_network");
and reload it in the same script with :
trained_network = load("./trained_networks/FiLM.mat");
trained_network = trained_network.trained_network;
I get the error :
Warning: While loading an object of class ‘dlnetwork’:
Error using nnet.internal.cnn.layer.GraphExecutor/propagate (line 354)
Execution failed during layer(s) ‘FiLM_64’.

Error in deep.internal.network.ExecutableNetwork/configureForInputsAndForwardOnLayer (line 347)
propagate(this, fcn, Xs, outputLayerIdx, outputLayerPortIdx);

Error in deep.internal.network.EditableNetwork/convertToDlnetwork (line 101)
[executableNetwork, layerOutputSizes] = configureForInputsAndForwardOnLayer(…

Error in dlnetwork.loadobj (line 741)
net = convertToDlnetwork(privateNet, exampleInputs, initializeNetworkWeights);

Caused by:
Undefined function handle.

Error in nnet.cnn.layer.FunctionLayer/predict (line 61)
[varargout{1:layer.NumOutputs}] = layer.PredictFcn(varargin{:});

I tried to save the function in a dedicated file "Feature_wise_LM.m" , but it didn’t workHi everyone,
I created a neural network that uses a functionLayer that I defined by :
functionLayer(@Feature_wise_LM,Name="FiLM_64",Formattable=1,NumInputs=2,NumOutputs=1)
However, when I save the trained neural network with :
save("./trained_networks/FiLM", "trained_network");
and reload it in the same script with :
trained_network = load("./trained_networks/FiLM.mat");
trained_network = trained_network.trained_network;
I get the error :
Warning: While loading an object of class ‘dlnetwork’:
Error using nnet.internal.cnn.layer.GraphExecutor/propagate (line 354)
Execution failed during layer(s) ‘FiLM_64’.

Error in deep.internal.network.ExecutableNetwork/configureForInputsAndForwardOnLayer (line 347)
propagate(this, fcn, Xs, outputLayerIdx, outputLayerPortIdx);

Error in deep.internal.network.EditableNetwork/convertToDlnetwork (line 101)
[executableNetwork, layerOutputSizes] = configureForInputsAndForwardOnLayer(…

Error in dlnetwork.loadobj (line 741)
net = convertToDlnetwork(privateNet, exampleInputs, initializeNetworkWeights);

Caused by:
Undefined function handle.

Error in nnet.cnn.layer.FunctionLayer/predict (line 61)
[varargout{1:layer.NumOutputs}] = layer.PredictFcn(varargin{:});

I tried to save the function in a dedicated file "Feature_wise_LM.m" , but it didn’t work Hi everyone,
I created a neural network that uses a functionLayer that I defined by :
functionLayer(@Feature_wise_LM,Name="FiLM_64",Formattable=1,NumInputs=2,NumOutputs=1)
However, when I save the trained neural network with :
save("./trained_networks/FiLM", "trained_network");
and reload it in the same script with :
trained_network = load("./trained_networks/FiLM.mat");
trained_network = trained_network.trained_network;
I get the error :
Warning: While loading an object of class ‘dlnetwork’:
Error using nnet.internal.cnn.layer.GraphExecutor/propagate (line 354)
Execution failed during layer(s) ‘FiLM_64’.

Error in deep.internal.network.ExecutableNetwork/configureForInputsAndForwardOnLayer (line 347)
propagate(this, fcn, Xs, outputLayerIdx, outputLayerPortIdx);

Error in deep.internal.network.EditableNetwork/convertToDlnetwork (line 101)
[executableNetwork, layerOutputSizes] = configureForInputsAndForwardOnLayer(…

Error in dlnetwork.loadobj (line 741)
net = convertToDlnetwork(privateNet, exampleInputs, initializeNetworkWeights);

Caused by:
Undefined function handle.

Error in nnet.cnn.layer.FunctionLayer/predict (line 61)
[varargout{1:layer.NumOutputs}] = layer.PredictFcn(varargin{:});

I tried to save the function in a dedicated file "Feature_wise_LM.m" , but it didn’t work deep learning, neural networks, function, matlab function MATLAB Answers — New Questions

​

I need programme format for simulating rotor inter-turn fault in synchronous generator.I need programme format for simulating rotor inter-turn fault in synchronous generator. I need programme format for simulating rotor inter-turn fault in synchronous generator. rotor inter-turn fault simulation. MATLAB Answers — New Questions

​

in R2025a I can not find anymore the "plot browser" menu. I used a lot the hide/show options for lines.
Has this function disappeared?
Property Inspector is not really helpfull if many lines are present in the plot.in R2025a I can not find anymore the "plot browser" menu. I used a lot the hide/show options for lines.
Has this function disappeared?
Property Inspector is not really helpfull if many lines are present in the plot. in R2025a I can not find anymore the "plot browser" menu. I used a lot the hide/show options for lines.
Has this function disappeared?
Property Inspector is not really helpfull if many lines are present in the plot. plotedit MATLAB Answers — New Questions

​

I have tried to simulate a PV module with Simulink for a change in radiation and I am unable to get the characteristics.The P-V and I-V characteristics I am expecting are <http://ecee.colorado.edu/~ecen2060/materials/simulink/PV/PV_module_model.pdf here> . I used the repeating sequence for my voltage source and the repeating sequence stair for my insolation. I get incomplete characteristics such as the figures below.

<</matlabcentral/answers/uploaded_files/192/P-V%20and%20I-V%20characteristicspage-0.jpg>>

_ . What could be the problem? I have changed the time values and output values for the repeating sequence source and the simulation time but without success. Anyone with an explanation?I have tried to simulate a PV module with Simulink for a change in radiation and I am unable to get the characteristics.The P-V and I-V characteristics I am expecting are <http://ecee.colorado.edu/~ecen2060/materials/simulink/PV/PV_module_model.pdf here> . I used the repeating sequence for my voltage source and the repeating sequence stair for my insolation. I get incomplete characteristics such as the figures below.

<</matlabcentral/answers/uploaded_files/192/P-V%20and%20I-V%20characteristicspage-0.jpg>>

_ . What could be the problem? I have changed the time values and output values for the repeating sequence source and the simulation time but without success. Anyone with an explanation? I have tried to simulate a PV module with Simulink for a change in radiation and I am unable to get the characteristics.The P-V and I-V characteristics I am expecting are <http://ecee.colorado.edu/~ecen2060/materials/simulink/PV/PV_module_model.pdf here> . I used the repeating sequence for my voltage source and the repeating sequence stair for my insolation. I get incomplete characteristics such as the figures below.

<</matlabcentral/answers/uploaded_files/192/P-V%20and%20I-V%20characteristicspage-0.jpg>>

_ . What could be the problem? I have changed the time values and output values for the repeating sequence source and the simulation time but without success. Anyone with an explanation? pv module, simulink, simulation. MATLAB Answers — New Questions

​

Hello
I am working on control design using ucover and musyn functions with a colleague, and we are sharing scripts to generate controllers. We see that the same script does not produce the same controller on our computers, even though all inputs and Matlab versions are identical (24.2.0.2923080 (R2024b) Update 6).
Basically, we input FRD objects into ucover, which outputs an uncertainty model, which is processed and in turn fed into musyn, whcih outputs the controller.
I’ve found that the output from musyn is different, even if I fetch all inputs to musyn from my colleagues workspace. So the difference is introduced in the function itself, likely related to some optimization or floating point processing.
I’ve seen some other posts dealing with this, for example this one:
https://se.mathworks.com/matlabcentral/answers/130493-how-come-i-get-different-output-answers-with-the-same-matlab-version-the-same-code-installed-on-two?s_tid=ta_ans_results
We both get the same output from this prompt, so I assume differing BLAS versions is not a problem (both apparently use AVX2):
>> version(‘-blas’)
ans =
‘Intel(R) oneAPI Math Kernel Library Version 2024.1-Product Build 20240215 for Intel(R) 64 architecture applications (CNR branch AVX2)’
I’ve also found that, if we restrict MATALB to 1 CPU by running
maxNumCompThreads(1)
the output from the ucover function becomes identical (at least for our test case). But the output from musyn is still different on our computers, even with identical input. The difference is large enough to be very significant from a control design perspective, so this is a bit frustrating since we want to be able to reproduce the same controllers in the future, using the script as a recipe (for tracability).
Any tips are greatly appreciated, thanks!Hello
I am working on control design using ucover and musyn functions with a colleague, and we are sharing scripts to generate controllers. We see that the same script does not produce the same controller on our computers, even though all inputs and Matlab versions are identical (24.2.0.2923080 (R2024b) Update 6).
Basically, we input FRD objects into ucover, which outputs an uncertainty model, which is processed and in turn fed into musyn, whcih outputs the controller.
I’ve found that the output from musyn is different, even if I fetch all inputs to musyn from my colleagues workspace. So the difference is introduced in the function itself, likely related to some optimization or floating point processing.
I’ve seen some other posts dealing with this, for example this one:
https://se.mathworks.com/matlabcentral/answers/130493-how-come-i-get-different-output-answers-with-the-same-matlab-version-the-same-code-installed-on-two?s_tid=ta_ans_results
We both get the same output from this prompt, so I assume differing BLAS versions is not a problem (both apparently use AVX2):
>> version(‘-blas’)
ans =
‘Intel(R) oneAPI Math Kernel Library Version 2024.1-Product Build 20240215 for Intel(R) 64 architecture applications (CNR branch AVX2)’
I’ve also found that, if we restrict MATALB to 1 CPU by running
maxNumCompThreads(1)
the output from the ucover function becomes identical (at least for our test case). But the output from musyn is still different on our computers, even with identical input. The difference is large enough to be very significant from a control design perspective, so this is a bit frustrating since we want to be able to reproduce the same controllers in the future, using the script as a recipe (for tracability).
Any tips are greatly appreciated, thanks! Hello
I am working on control design using ucover and musyn functions with a colleague, and we are sharing scripts to generate controllers. We see that the same script does not produce the same controller on our computers, even though all inputs and Matlab versions are identical (24.2.0.2923080 (R2024b) Update 6).
Basically, we input FRD objects into ucover, which outputs an uncertainty model, which is processed and in turn fed into musyn, whcih outputs the controller.
I’ve found that the output from musyn is different, even if I fetch all inputs to musyn from my colleagues workspace. So the difference is introduced in the function itself, likely related to some optimization or floating point processing.
I’ve seen some other posts dealing with this, for example this one:
https://se.mathworks.com/matlabcentral/answers/130493-how-come-i-get-different-output-answers-with-the-same-matlab-version-the-same-code-installed-on-two?s_tid=ta_ans_results
We both get the same output from this prompt, so I assume differing BLAS versions is not a problem (both apparently use AVX2):
>> version(‘-blas’)
ans =
‘Intel(R) oneAPI Math Kernel Library Version 2024.1-Product Build 20240215 for Intel(R) 64 architecture applications (CNR branch AVX2)’
I’ve also found that, if we restrict MATALB to 1 CPU by running
maxNumCompThreads(1)
the output from the ucover function becomes identical (at least for our test case). But the output from musyn is still different on our computers, even with identical input. The difference is large enough to be very significant from a control design perspective, so this is a bit frustrating since we want to be able to reproduce the same controllers in the future, using the script as a recipe (for tracability).
Any tips are greatly appreciated, thanks! ucover, musyn, avx2 MATLAB Answers — New Questions

​

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.

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(:)];
end 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(:)];
end adsorption, hsdm, resin, adsorption column, lithium MATLAB Answers — New Questions

​