# Chromatography PDE MOL

• MATLAB
Hello all

I am using the method of lines to solve the following PDE:

## \frac {\partial C} {\partial t} + F\frac {\partial q} {\partial t} + u \frac {dC} {dz} = D_{ax} \frac{\partial^2 C} {\partial z^2} ##

## \frac {\partial q} {\partial t} = k (q^{*}-q) ##

With these initial conditions:
## C(0,z) = 0 ##
## q(0,z) = 0 ##

and boundary conditions:
## C(t,0) = C_{in} (1) ##
## \frac {\partial C} {\partial z}(t,L) = 0 (2) ##

I discretized the equations above in space and obtained a system of ODE's which I solved with ode15s. To deal with the boundary condition (2) I said:

## \frac{C(t,L+1) - C(t,L-1)}{2\Delta z}= 0 ## or ## C(t,L+1) = C(t,L-1) ##

The problem is that my system gives a solution which tends to explode and makes no sense. My guess is that the way I am discretizing the system is not correct or consistent. I am using between 50-150 nodes.

My gradient matrix for the 1st derivative (GM1) looks like:

\begin{bmatrix}
0 & \frac{1}{2} & 0 & 0 & 0 & 0 & 0 & ... & 0\\
\frac{-8}{12} & 0 & \frac{8}{12} & \frac{-1}{12} & 0 & 0 & 0 & ... & 0\\
\frac{1}{12} & \frac{-8}{12} & 0 & \frac{8}{12} & \frac{-1}{12} & 0 & 0 & ... & 0 \\
0 & \frac{1}{12} & \frac{-8}{12} & 0 & \frac{8}{12} & \frac{-1}{12} & 0 & ... & 0 \\
... & ... & ... & ... & ... & ... & ... & ... & ... \\
0 & ... & 0 & 0 & 0 & \frac{1}{12} & \frac{-8}{12} & \frac{-1}{12} & \frac{8}{12} \\
0 & ... & 0 & 0 & 0 & 0 & 0 & 0 & 1
\end{bmatrix}

and b1 = \begin{bmatrix}
-\frac{C_{in}}{2} \\
\frac{C_{in}}{12} \\
0 \\
... \\
0 \\
-C_{in}
\end{bmatrix}

For the 2nd derivative (GM2):

\begin{bmatrix}
-2 & 1 & 0 & 0 & 0 & 0 & 0 & ... & 0\\
1 & -2 & 1 & 0 & 0 & 0 & 0 & ... & 0\\
0 & 1 & -2 & 1 & 0 & 0 & 0 & ... & 0 \\
... & ... & ... & ... & ... & ... & ... & ... & ... \\
0 & ... & 0 & 0 & 0 & 0 & 0 & 2 & -2
\end{bmatrix}

and b2 = \begin{bmatrix}
C_{in} \\
0 \\
... \\
0
\end{bmatrix}

In the end I say that:

## \frac {\partial C} {\partial z} = \frac{1}{h}GM_{1} \cdot C + \frac{1}{h}b_{1}## and ## \frac {\partial^2 C} {\partial z^2} = \frac{1}{h^2}GM_{2} \cdot C + \frac{1}{h^2}b_{2}##

Last edited:

My code:

FUNCTION 1:

Code:
function [sol, t, x, C] = LDF_mol(isoType, feedProf, parameter, L, Di, epsb, Q, Cfeed, KLDF, Dax, tpulse, tfinal, opt)
% Transport-Dispersive Model (TDM) considering mass transfer resistence in the solid to be dominant and
% using the Linear Driving Force Model (LDF) approach (Glauckauf and Coates, 1947)
% Change the isotherm model by modifying the isotherm function
% Change feed profile between pulse (e.g.: chromatografic peak) and step (e.g.: breakthrough experiment)
% Multicomponent
% Uses pdepe functin to solve the system of partial differential equations

%% Default arguments (Example)
if nargin == 0
isoType =   'mod-langmuir';                                                     % isotherm type, can be 'linear' or 'linear-langmuir'
feedProf =  'pulse';                                                            % feed profile, can be 'pulse' (e.g.: chromatografic peak) or 'step' (e.g.: breakthrough experiment)
parameter = [1e3 0.3/1e3 0.05; 1e3 0.1/1e3 0.045];                               % isotherm parameters (depends on the isotherm model chosen)
L =         10;                                                                 % cm, column length
Di =        1;                                                                  % cm, column internal diameter
epsb =      0.4;                                                                % column bulk porosity
Q =         4;                                                                  % mL/min, flow rate
Cfeed =    [0.6*300/50 ; 0.6*300/50];                             % g/L, feed concentration
KLDF =      [60; 60];                                                    % min-1, linear driving force (LDF) mass transfer coefficient
Dax =       [6.57e-3; 7.57e-3];                                     % cm2/min, axial dispersion coefficient
tpulse =    0.1;                                                                % min, feed pulse duration. For a step injection set tpulse = tfinal
tfinal =    5;                                                                  % min, final time for calculation
opt.npz =   50;                                                                % number of discretization points in z
opt.npt =   100;                                                                % number of discretization points in t
opt.fig =   true;                                                               % true - show figures; false - do not show figures

Cs_wash = 0; % Salt concentration in mM
twash = 1;

t0_elu = tpulse+twash; telu = 8; Cf = 80;
Cs_elu = @(t) 0+(Cf-0)/telu*(t-t0_elu);

Cs_reg = 80;
treg = 3;

Cs_eq = 0;
teq = 3;
end

%% Calculations
data = struct('isoType',isoType,'feedProf',feedProf,'parameter',parameter,'Cfeed',Cfeed,'Cs_wash',Cs_wash,'Cs_elu',Cs_elu,'Cs_reg',Cs_reg,'Cs_eq',Cs_eq,'npz',opt.npz,'npt',opt.npt,'Q',Q,'L',L,'epsb',epsb,'tpulse',tpulse,'twash',twash,'t0_elu',t0_elu,'telu',telu,'treg',treg,'teq',teq,'KLDF',KLDF,'Dax',Dax);
nc = length(Cfeed);
data.nc = nc;
A = pi()*Di^2/4; % cm2
data.F = (1-epsb)/epsb;
data.ui = Q/A; % cm/min

% Run MOL
t = linspace(0,tfinal,opt.npt);
z = linspace(0,L,opt.npz);

u0 = zeros(2*nc*opt.npz,1);

[t,uout] = ode15s(@odefun, t, u0);
for CompNo = 1:data.nc
idx1 = 1+(CompNo-1)*opt.npz; idx2 = CompNo*opt.npz;

C{CompNo} = uout(:,idx1:idx2);
end
res.t = t'; res.z = z'; res.C = C;

Cs_col = zeros(opt.npt,opt.npz);

Cs_col(:,1) = Cs_profile(t);
for i = 1:eek:pt.npt
for j = 2:eek:pt.npz
if t(i)-z(j)/data.ui >= 0
Cs_col(i,j) = Cs_profile(t(i)-z(j)/data.ui);
else
Cs_col(i,j) = 0;
end
end
end

ymat = [];
for CompNo = 1:data.nc
ymat = [ymat C{CompNo}(:,end)];
end
[ax,h1,h2] = plotyy(t,ymat,t,Cs_col(:,1));
axis([0 tfinal  0 inf]) % fix the axes
set(ax(1),'Box','off'); set(ax(1),'YTick',0:0.2:2); set(ax(2),'YTick',0:20:100,'YLim',[0 100]);
xlabel('t')
ylabel('C')

% Concentration over time with slider figure
figure
hold all
slidermin = 1; % t = 0
slidermax = opt.npt; % t = tfinal
ymax_aux = [];
for CompNo = 1:data.nc-1
res.hslider(CompNo) = plot(z',C{CompNo}(end,:));
ymax_aux = [ymax_aux max(C{CompNo}(2,:))];
end
ymax_aux = [ymax_aux max(C{data.nc}(2,:))];
res.ymax = str2double(num2str(max(ymax_aux),1));
axis([0 L 0 2]);

[ax,h1,h2] = plotyy(z',C{data.nc}(end,:),z',Cs_col(end,:));
set(ax(1),'Box','off'); set(ax(2),'YTick',0:20:100,'YLim',[0 100]);
hold off
res.hslider(data.nc) = h1;
res.hslider(data.nc+1) = h2;

xlabel('Position, z')
ylabel('Concentration, C')
anno_text = sprintf('t = %.1f', 0);
res.time_box = annotation('textbox',[.55 .48 .4 .5],'String',anno_text,'EdgeColor', 'none','FitBoxToText','on');
sld = uicontrol('Style', 'slider',...
'Min',slidermin,'Max',slidermax,'Value',slidermin,...
'SliderStep',[1 1]./(slidermax-slidermin),...
'Position', [390 390 120 20],...
'Callback', @sliderplot);

%% Auxiliary functions
function sliderplot(source,event)
% Handles the UI slider in the Concentration vs Position plot
global data res Cs_col

val = round(get(source,'Value'));
for CompNo = 1:data.nc
set(res.hslider(CompNo),'YData',res.C{CompNo}(val,:));
end
set(res.hslider(data.nc+1),'YData',Cs_col(val,:));

anno_text = sprintf('t = %.2f',res.t(val));
set(res.time_box,'String',anno_text);

function udot = odefun(t,u)

nc = data.nc; npz = data.npz; L = data.L;
z = linspace(0,L,npz); h = L/(npz-1);

Cin = setFeedProfile(data.feedProf, t, data.tpulse, data.Cfeed);

C = u(1:nc*npz);
for CompNo = 1:nc
idx1 = 1+(CompNo-1)*npz;
C(idx1) = Cin(CompNo);
end

q = u(nc*npz+1:end);
Cs = zeros(npz,1);

udot = zeros(2*nc*npz,1);
for j = 1:npz
if t-z(j)/data.ui >= 0
Cs(j) = Cs_profile(t-z(j)/data.ui);
else
Cs(j) = 0;
end
end

qeq = isotherm1(data.isoType, C, Cs, nc, data.parameter);

dCdz = zeros(nc*npz,1); d2Cdz2 = zeros(nc*npz,1);
for CompNo = 1:nc
idx1 = 1+(CompNo-1)*npz; idx2 = CompNo*npz;
idx3 = npz*(nc+CompNo-1)+1; idx4 = (nc+CompNo)*npz;

b1 = zeros(npz,1); b1(2) = -Cin(CompNo)/2; b1(3) = Cin(CompNo)/12;
b1(end) = -C(idx2);
b2 = zeros(npz,1); b2(2) = Cin(CompNo);

udot(idx3:idx4) = data.KLDF(CompNo)*(qeq(:,CompNo)-q(idx1:idx2));
udot(idx1:idx2) = -data.ui*dCdz(idx1:idx2) - data.F*udot(idx3:idx4) + data.Dax(CompNo)*d2Cdz2(idx1:idx2);

udot(idx1)=0; udot(idx3) = 0;
end

%% Feed profile
function Cinj=setFeedProfile(feedProf, t, tpulse, Cfeed)
% Defines the feed profile. Can be 'pulse' (e.g.: chromatografic peak) or 'step' (e.g.: breakthrough experiment)

if strcmp(feedProf,'pulse')
if t<=tpulse
Cinj = Cfeed;
else
Cinj = zeros(1,length(Cfeed));
end

elseif strcmp(feedProf,'step')
Cinj = Cfeed;

else
error('Invalid feed profile. feedProf must be "step" or "pulse"')
end

function Cs_in = Cs_profile(t)
global data

Cs_in = zeros(length(t),1);

for i = 1:length(t)

if t(i) <= data.tpulse
Cs_in(i) = 0;
elseif t(i) <= data.tpulse+data.twash
Cs_in(i) = data.Cs_wash;
elseif t(i) <= data.tpulse+data.twash+data.telu
Cs_in(i) = data.Cs_elu(t(i));
elseif t(i) <= data.tpulse+data.twash+data.telu+data.treg
Cs_in(i) = data.Cs_reg;
else
Cs_in(i) = data.Cs_eq;
end
end

for i = 4:N-2
Grad1(i,i-3:i+1) = [1/12, -8/12, 0, 8/12, -1/12];
end

for i = 3:N-1
end

FUNCTION 2

Code:
%% Isotherm Library

function q=isotherm1(isoType, C, Cs, nc, parameter)
% Defines the isotherm
parameter = transpose(parameter); N = length(Cs);

q = zeros(N,nc);
Cmat = zeros(N,nc);
for i = 1:nc
idx1 = 1+(i-1)*N; idx2 = i*N;
Cmat(:,i) = C(idx1:idx2);
end

if strcmp(isoType,'mod-langmuir')
% Linear-Langmuir isotherm:   q1(C1,C2) = a1*b1*exp(-m1*Cs)*C1/(1+b1*exp(-m1*Cs)*C1) , parameter = [a1 ; b1 ; m1]
for i = 1:N
langmuirDenominator = 1;
for j = 1:nc
langmuirDenominator = langmuirDenominator +  parameter(2,j)*exp(-parameter(3,j)*Cs(i))*Cmat(i,j);
end
q(i,:) = parameter(1,:).*parameter(2,:).*exp(-parameter(3,:)*Cs(i)).*Cmat(i,:)/langmuirDenominator;
end

else
error('Invalid isotherm type. isoType must be "linear" or "langmuir" or "linear-langmuir"')
end

I think now I kind of solved it. However when the wave reaches the right boundary there is some noise in the wave. Do you know what could be the reason?

Code:
function [sol, t, x, C] = LDF_mol(isoType, feedProf, parameter, L, Di, epsb, Q, Cfeed, KLDF, Dax, tpulse, tfinal, opt)
% Transport-Dispersive Model (TDM) considering mass transfer resistence in the solid to be dominant and
% using the Linear Driving Force Model (LDF) approach (Glauckauf and Coates, 1947)
% Change the isotherm model by modifying the isotherm function
% Change feed profile between pulse (e.g.: chromatografic peak) and step (e.g.: breakthrough experiment)
% Multicomponent
% Uses pdepe functin to solve the system of partial differential equations

%% Default arguments (Example)
if nargin == 0
isoType =   'mod-langmuir';                                                     % isotherm type, can be 'linear' or 'linear-langmuir'
feedProf =  'pulse';                                                            % feed profile, can be 'pulse' (e.g.: chromatografic peak) or 'step' (e.g.: breakthrough experiment)
parameter = [1e3 1.5/1e3 0.15; 1e3 2/1e3 0.045];                               % isotherm parameters (depends on the isotherm model chosen)
L =         10;                                                                 % cm, column length
Di =        1;                                                                  % cm, column internal diameter
epsb =      0.4;                                                                % column bulk porosity
Q =         4;                                                                  % mL/min, flow rate
Cfeed =    [0.6*300/50 ; 0.6*300/50];                             % g/L, feed concentration
KLDF =      [60; 60];                                                    % min-1, linear driving force (LDF) mass transfer coefficient
Dax =       [6.57e-3; 7.57e-3];                                     % cm2/min, axial dispersion coefficient
tpulse =    0.1;                                                                % min, feed pulse duration. For a step injection set tpulse = tfinal
tfinal =    8;                                                                  % min, final time for calculation
opt.npz =   300;                                                                % number of discretization points in z
opt.npt =   100;                                                                % number of discretization points in t
opt.fig =   true;                                                               % true - show figures; false - do not show figures

Cs_wash = 0; % Salt concentration in mM
twash = 1;

t0_elu = tpulse+twash; telu = 8; Cf = 80;
Cs_elu = @(t) 0+(Cf-0)/telu*(t-t0_elu);

Cs_reg = 80;
treg = 3;

Cs_eq = 0;
teq = 3;
end

%% Calculations
data = struct('isoType',isoType,'feedProf',feedProf,'parameter',parameter,'Cfeed',Cfeed,'Cs_wash',Cs_wash,'Cs_elu',Cs_elu,'Cs_reg',Cs_reg,'Cs_eq',Cs_eq,'npz',opt.npz,'npt',opt.npt,'Q',Q,'L',L,'epsb',epsb,'tpulse',tpulse,'twash',twash,'t0_elu',t0_elu,'telu',telu,'treg',treg,'teq',teq,'KLDF',KLDF,'Dax',Dax);
nc = length(Cfeed);
data.nc = nc;
A = pi()*Di^2/4; % cm2
data.F = (1-epsb)/epsb;
data.ui = Q/A; % cm/min

% Run MOL
t = linspace(0,tfinal,opt.npt);
z = linspace(0,L,opt.npz);

u0 = zeros(2*nc*opt.npz,1);

[t,uout] = ode15s(@odefun, t, u0);
for CompNo = 1:data.nc
idx1 = 1+(CompNo-1)*opt.npz; idx2 = CompNo*opt.npz;

C{CompNo} = uout(:,idx1:idx2);
end
res.t = t'; res.z = z'; res.C = C;

Cs_col = zeros(opt.npt,opt.npz);

Cs_col(:,1) = Cs_profile(t);
for i = 1:opt.npt
for j = 2:opt.npz
if t(i)-z(j)/data.ui >= 0
Cs_col(i,j) = Cs_profile(t(i)-z(j)/data.ui);
else
Cs_col(i,j) = 0;
end
end
end

ymat = [];
for CompNo = 1:data.nc
ymat = [ymat C{CompNo}(:,end)];
end
[ax,h1,h2] = plotyy(t,ymat,t,Cs_col(:,1));
axis([0 tfinal  0 inf]) % fix the axes
set(ax(1),'Box','off'); set(ax(1),'YTick',0:0.2:2); set(ax(2),'YTick',0:20:100,'YLim',[0 100]);
xlabel('t')
ylabel('C')

% Concentration over time with slider figure
figure
hold all
slidermin = 1; % t = 0
slidermax = opt.npt; % t = tfinal
ymax_aux = [];
for CompNo = 1:data.nc-1
res.hslider(CompNo) = plot(z',C{CompNo}(end,:));
ymax_aux = [ymax_aux max(C{CompNo}(2,:))];
end
ymax_aux = [ymax_aux max(C{data.nc}(2,:))];
res.ymax = str2double(num2str(max(ymax_aux),1));
axis([0 L 0 2]);

[ax,h1,h2] = plotyy(z',C{data.nc}(end,:),z',Cs_col(end,:));
set(ax(1),'Box','off'); set(ax(2),'YTick',0:20:100,'YLim',[0 100]);
hold off
res.hslider(data.nc) = h1;
res.hslider(data.nc+1) = h2;

xlabel('Position, z')
ylabel('Concentration, C')
anno_text = sprintf('t = %.1f', 0);
res.time_box = annotation('textbox',[.55 .48 .4 .5],'String',anno_text,'EdgeColor', 'none','FitBoxToText','on');
sld = uicontrol('Style', 'slider',...
'Min',slidermin,'Max',slidermax,'Value',slidermin,...
'SliderStep',[1 1]./(slidermax-slidermin),...
'Position', [390 390 120 20],...
'Callback', @sliderplot);

%% Auxiliary functions
function sliderplot(source,event)
% Handles the UI slider in the Concentration vs Position plot
global data res Cs_col

val = round(get(source,'Value'));
for CompNo = 1:data.nc
set(res.hslider(CompNo),'YData',res.C{CompNo}(val,:));
end
set(res.hslider(data.nc+1),'YData',Cs_col(val,:));

anno_text = sprintf('t = %.2f',res.t(val));
set(res.time_box,'String',anno_text);

function udot = odefun(t,u)

nc = data.nc; npz = data.npz; L = data.L;
z = linspace(0,L,npz); h = L/(npz-1);

Cin = setFeedProfile(data.feedProf, t, data.tpulse, data.Cfeed);

C = u(1:nc*npz);
for CompNo = 1:nc
idx1 = 1+(CompNo-1)*npz;
C(idx1) = Cin(CompNo);
end

q = u(nc*npz+1:end);
Cs = zeros(npz,1);

udot = zeros(2*nc*npz,1);
for j = 1:npz
if t-z(j)/data.ui >= 0
Cs(j) = Cs_profile(t-z(j)/data.ui);
else
Cs(j) = 0;
end
end

qeq = isotherm1(data.isoType, C, Cs, nc, data.parameter);

dCdz = zeros(nc*npz,1); d2Cdz2 = zeros(nc*npz,1);
for CompNo = 1:nc
idx1 = 1+(CompNo-1)*npz; idx2 = CompNo*npz;
idx3 = npz*(nc+CompNo-1)+1; idx4 = (nc+CompNo)*npz;

b1 = zeros(npz-1,1); b1(1) = -Cin(CompNo)/2; b1(2) = Cin(CompNo)/12;
b1(end) = -C(idx2);
b2 = zeros(npz-1,1); b2(1) = Cin(CompNo);

dCdz(idx1) = 3.14; % Value does not matter
d2Cdz2(idx1) = 3.14; % Value does not matter

udot(idx3:idx4) = data.KLDF(CompNo)*(qeq(:,CompNo)-q(idx1:idx2));
udot(idx1:idx2) = -data.ui*dCdz(idx1:idx2) - data.F*udot(idx3:idx4) + data.Dax(CompNo)*d2Cdz2(idx1:idx2);

udot(idx1)=0; udot(idx3) = 0;
end

%% Feed profile
function Cinj=setFeedProfile(feedProf, t, tpulse, Cfeed)
% Defines the feed profile. Can be 'pulse' (e.g.: chromatografic peak) or 'step' (e.g.: breakthrough experiment)

if strcmp(feedProf,'pulse')
if t<=tpulse
Cinj = Cfeed;
else
Cinj = zeros(1,length(Cfeed));
end

elseif strcmp(feedProf,'step')
Cinj = Cfeed;

else
error('Invalid feed profile. feedProf must be "step" or "pulse"')
end

function Cs_in = Cs_profile(t)
global data

Cs_in = zeros(length(t),1);

for i = 1:length(t)

if t(i) <= data.tpulse
Cs_in(i) = 0;
elseif t(i) <= data.tpulse+data.twash
Cs_in(i) = data.Cs_wash;
elseif t(i) <= data.tpulse+data.twash+data.telu
Cs_in(i) = data.Cs_elu(t(i));
elseif t(i) <= data.tpulse+data.twash+data.telu+data.treg
Cs_in(i) = data.Cs_reg;
else
Cs_in(i) = data.Cs_eq;
end
end

%
% for i = 4:N-2
%     Grad1(i,i-3:i+1) = [1/12, -8/12, 0, 8/12, -1/12];
% end

for i = 3:N-3
Grad1(i,i-2:i+2) = [1/12, -8/12, 0, 8/12, -1/12];
end

end