fixed bed adsorption column model-solving PDE-freundlish isotherm

function main
%System Set-up %
%Define Variables
%D = 3*10^-8;    % Axial Dispersion coefficient
%v = 1*10^-3;    % Superficial velocity
%epsilon = 0.4;  % Voidage fraction
%k = 3*10^-5;    % Mass Transfer Coefficient
%Kf=2.5*10^-5; %freundlish parameter
%nf= 1.45; %freundlish constant
cFeed = 10;     % Feed concentration
L = 1;         % Column length
t0 = 0;        % Initial Time
%tf = 1000;      % Final time
tf = 2000;      % Final time
dt = 0.5;      % Time step
% z=[0:0.01:L]; %Mesh generation
z = [0:0.0005:L]; %Mesh generation
t = [t0:dt:tf];% Time vector
n = numel(z);  % Size of mesh grid
%Initial Conditions / Vector Creation
c0 = zeros(n,1);     
c0(1) = cFeed;
q0 = zeros(n,1);    % t = 0, q = 0 for all z, this makes sense to me
y0 = [c0 ; q0];     % Appends conditions together
%ODE15S Solver
[T, Y] = ode15s(@(t,y) MyFun(t,y,z,n),t,y0);
%plot(T,Y);
plot(T,Y(:,n)/cFeed)
end 
function DyDt=MyFun(~, y, z, n)
% Defining Constants
D = 3*10^-8;    % Axial Dispersion coefficient
v = 1*10^-3;    % Superficial velocity
epsilon = 0.4;  % Voidage fraction
k = 3*10^-5;    % Mass Transfer Coefficient
Kf=2.5*10^-5; %freundlish parameter
nf= 1.45; %freundlish constant
% Variables being allocated zero vectors
c = zeros(n,1);
q = zeros(n,1);
DcDt = zeros(n,1);
DqDt = zeros(n,1);
DyDt = zeros(2*n,1);
zhalf = zeros(n-1,1);
DcDz = zeros(n,1);
D2cDz2 = zeros(n,1);
c = y(1:n);
q = y(n+1:2*n);
% Interior mesh points
zhalf(1:n-1)=(z(1:n-1)+z(2:n))/2;
for i=2:n-1
  %DcDz(i) = ((z(i)-z(i-1))/(z(i+1)-z(i))*(c(i+1)-c(i))+(z(i+1)-z(i))/(z(i)-z(i-1))*(c(i)-c(i-1)))/(z(i+1)-z(i-1));
  DcDz(i) = (c(i)-c(i-1))/(z(i)-z(i-1));
  D2cDz2(i) = (zhalf(i)*(c(i+1)-c(i))/(z(i+1)-z(i))-zhalf(i-1)*(c(i)-c(i-1))/(z(i)-z(i-1)))/(zhalf(i)-zhalf(i-1));
end
% Calculate DcDz and D2cDz2 at z=L for boundary condition dc/dz = 0
DcDz(n) = 0; 
D2cDz2(n) = -1.0/(z(n)-zhalf(n-1))*(c(n)-c(n-1))/(z(n)-z(n-1));
% Set time derivatives at z=0
% DcDt = 0 since c=cFeed remains constant over time and is already set as initial condition
% Standard setting for q 
DqDt(1) = k*(kf*(c(1)^(1/nf))-q(1));
DcDt(1) = 0.0;
% Set time derivatives in remaining points
for i=2:n
    %Equation: dq/dt = K(q*-q) where q*=kf*(c^(1/nf))
    DqDt(i) = k*(kf*(c(i)^(1/nf))-q(i));
      %Equation: dc/dt = D * d2c/dz2 - v*dc/dz - ((1-e)/(e))*dq/dt
      DcDt(i) = D*D2cDz2(i) - v*DcDz(i) - ((1-epsilon)/(epsilon))*DqDt(i);
  end
% Concatenate vector of time derivatives
DyDt = [DcDt;DqDt];
end

This is the code which I am using. I don't have much knowledge about Matlab but I need to finish the code for my thesis. My question is, I need to model a fixed bed to obtain the breakthrough curves using freundlish isotherm, so, I would like to know how to obtain these curves (C/C0 vs t) from this code. i cant run this code i dont know where is the problem. I would appreciate any help. Thank you very much

[darova]

can you attach original equations?

[nashwa tarek]

[nashwa tarek]

darova thanks for your help, waiting your feedback

Answer 1

[Alan Stevens]

The following works. I can't say if the result is sensible or not - you will need to judge that. I'm not convinced that the coded equations match those in your hand-written attachment!

cFeed = 10;      % Feed concentration
L = 1;           % Column length
t0 = 0;          % Initial Time
tf = 2000;       % Final time
dt = 0.5;        % Time step
z = 0:0.005:L;   % Mesh generation
t = t0:dt:tf;    % Time vector
n = numel(z);    % Size of mesh grid

%Initial Conditions / Vector Creation
c0 = zeros(n,1);     
c0(1) = cFeed;
q0 = zeros(n,1);    % t = 0, q = 0 for all z, this makes sense to me
y0 = [c0 ; q0];     % Appends conditions together

%ODE15S Solver
[T, Y] = ode15s(@(t,y) MyFun(t,y,z,n),t,y0);

%plot(T,Y);
plot(T,Y(:,n)/cFeed)

function DyDt=MyFun(~, y, z, n)
% Defining Constants
D = 3*10^-8;    % Axial Dispersion coefficient
v = 1*10^-3;    % Superficial velocity
epsilon = 0.4;  % Voidage fraction
k = 3*10^-5;    % Mass Transfer Coefficient
Kf=2.5*10^-5;   %freundlish parameter
nf= 1.45;       %freundlish constant

% Variables being allocated zero vectors
c = zeros(n,1);
q = zeros(n,1);
DcDt = zeros(n,1);
DqDt = zeros(n,1);
DyDt = zeros(2*n,1);
zhalf = zeros(n-1,1);
DcDz = zeros(n,1);
D2cDz2 = zeros(n,1);
c = y(1:n);
q = y(n+1:2*n);

% Interior mesh points
zhalf(1:n-1)=(z(1:n-1)+z(2:n))/2;
for i=2:n-1
  %DcDz(i) = ((z(i)-z(i-1))/(z(i+1)-z(i))*(c(i+1)-c(i))+(z(i+1)-z(i))/(z(i)-z(i-1))*(c(i)-c(i-1)))/(z(i+1)-z(i-1));
  DcDz(i) = (c(i)-c(i-1))/(z(i)-z(i-1));
  D2cDz2(i) = (zhalf(i)*(c(i+1)-c(i))/(z(i+1)-z(i))-zhalf(i-1)*(c(i)-c(i-1))/(z(i)-z(i-1)))/(zhalf(i)-zhalf(i-1));
end
% Calculate DcDz and D2cDz2 at z=L for boundary condition dc/dz = 0
DcDz(n) = 0; 
D2cDz2(n) = -1.0/(z(n)-zhalf(n-1))*(c(n)-c(n-1))/(z(n)-z(n-1));
% Set time derivatives at z=0
% DcDt = 0 since c=cFeed remains constant over time and is already set as initial condition
% Standard setting for q 
DqDt(1) = k*(Kf*(c(1)^(1/nf))-q(1));
DcDt(1) = 0.0;


% Set time derivatives in remaining points
for i=2:n
    %Equation: dq/dt = K(q*-q) where q*=kf*(c^(1/nf))
    DqDt(i) = k*(Kf*(c(i)^(1/nf))-q(i));
     %Equation: dc/dt = D * d2c/dz2 - v*dc/dz - ((1-e)/(e))*dq/dt
    DcDt(i) = D*D2cDz2(i) - v*DcDz(i) - ((1-epsilon)/(epsilon))*DqDt(i);
end

% Concatenate vector of time derivatives
% Very small imaginary components need to be removed!
DyDt = [real(DcDt);real(DqDt)];
end

[nashwa tarek]

thank you so much Alan, i am not that good in matlab i am trying to write a matlab code which fits my equations to get the breakthrough curve for a fixed bed adsorption coumn(c/c0 vs time). i did some batch experimental work and i found that the freundlish isotherm fits better

i dont know if i applied my equations right or not in matlab, can you help

[Alan Stevens]

Well, if you extend tf (the end time) to 4000, then it reaches 1.

As for getting the equations right, I would need the original equations written clearly, with each parameter defined, and the data for each supplied. It's possible you have implemented them correctly already, but there seemed some strange mismatches. For example your code has

DqDt(i) = k*(Kf*(c(i)^(1/nf))-q(i));

but in your .jpg file the RHS has only C's not q's. Perhaps you did some manipulation before coding to get it in this form, but it isn't obvious to me.

[nashwa tarek]

hi Alan

first of all thank you so much really i cant thank you more, you helped me alot after i lost hope:)))

the answer of your question is that DqDt(i)= K*(q*-q) and q*= Kf(c(i)^1/nf) so i just substituted .

now i still need you help , i need to add this equation to my code (pore diffusion equation) attached where Ep is particle porosity =0.53 and sp = particle denisty= 1970

i want to add this equation and get the breakthrough curve t vs c/co

thanks alot

[Alan Stevens]

Like the following perhaps. Because dqdt(i) and dcdt(i) have to be solved simultaneously, I've added some extrra variables to make the equations easier to write. You will need to check what I've done very carefully!

%System Set-up %
%Define Variables
%D = 3*10^-8;    % Axial Dispersion coefficient
%v = 1*10^-3;    % Superficial velocity
%epsilon = 0.4;  % Voidage fraction
%k = 3*10^-5;    % Mass Transfer Coefficient
%Kf=2.5*10^-5; %freundlish parameter
%nf= 1.45; %freundlish constant
% Ep = 0.53; pore porosity
% sp = 1970; particle density

cFeed = 10;      % Feed concentration
L = 1;           % Column length
t0 = 0;          % Initial Time
tf = 4000;       % Final time
dt = 0.5;        % Time step
z = 0:0.005:L;   % Mesh generation
t = t0:dt:tf;    % Time vector
n = numel(z);    % Size of mesh grid

%Initial Conditions / Vector Creation
c0 = zeros(n,1);     
c0(1) = cFeed;
q0 = zeros(n,1);    % t = 0, q = 0 for all z, this makes sense to me
y0 = [c0 ; q0];     % Appends conditions together

%ODE15S Solver
[T, Y] = ode15s(@(t,y) MyFun(t,y,z,n),t,y0);

%plot(T,Y);
plot(T,Y(:,n)/cFeed)
legend('c')
figure
plot(T,Y(:,2*n))
legend('q')

function DyDt=MyFun(~, y, z, n)
% Defining Constants
D = 3*10^-8;    % Axial Dispersion coefficient
v = 1*10^-3;    % Superficial velocity
epsilon = 0.4;  % Voidage fraction
k = 3*10^-5;    % Mass Transfer Coefficient
Kf=2.5*10^-5;   % freundlish parameter
nf= 1.45;       % freundlish constant
Ep = 0.53;      % pore porosity
sp = 1970;      % particle density

% Variables being allocated zero vectors
c = zeros(n,1);
q = zeros(n,1);
DcDt = zeros(n,1);
DqDt = zeros(n,1);
DyDt = zeros(2*n,1);
zhalf = zeros(n-1,1);
DcDz = zeros(n,1);
D2cDz2 = zeros(n,1);
c = y(1:n);
q = y(n+1:2*n);

% Interior mesh points
zhalf(1:n-1)=(z(1:n-1)+z(2:n))/2;
for i=2:n-1
  %DcDz(i) = ((z(i)-z(i-1))/(z(i+1)-z(i))*(c(i+1)-c(i))+(z(i+1)-z(i))/(z(i)-z(i-1))*(c(i)-c(i-1)))/(z(i+1)-z(i-1));
  DcDz(i) = (c(i)-c(i-1))/(z(i)-z(i-1));
  D2cDz2(i) = (zhalf(i)*(c(i+1)-c(i))/(z(i+1)-z(i))-zhalf(i-1)*(c(i)-c(i-1))/(z(i)-z(i-1)))/(zhalf(i)-zhalf(i-1));
end
% Calculate DcDz and D2cDz2 at z=L for boundary condition dc/dz = 0
DcDz(n) = 0; 
D2cDz2(n) = -1.0/(z(n)-zhalf(n-1))*(c(n)-c(n-1))/(z(n)-z(n-1));
% Set time derivatives at z=0
% DcDt = 0 since c=cFeed remains constant over time and is already set as initial condition
% Standard setting for q 
DqDt(1) = k*(Kf*(c(1)^(1/nf))-q(1));
DcDt(1) = 0.0;


% Set time derivatives in remaining points
for i=2:n
    % revised equations
    alpha = k*(Kf*(c(i)^(1/nf))-q(i)); beta = - Ep/((1-Ep)*sp);
    gamma = D*D2cDz2(i) - v*DcDz(i);  delta = - ((1-epsilon)/(epsilon));
    DqDt(i) = (alpha + beta.*gamma)./(1 - beta.*delta);
    DcDt(i) = (gamma + alpha.*delta)./(1 - beta.*delta);
    
    % original eqns 
    %Equation: dq/dt = K(q*-q) where q*=kf*(c^(1/nf))
    % DqDt(i) = k*(Kf*(c(i)^(1/nf))-q(i));
    %Equation: dc/dt = D * d2c/dz2 - v*dc/dz - ((1-e)/(e))*dq/dt
    % DcDt(i) = D*D2cDz2(i) - v*DcDz(i) - ((1-epsilon)/(epsilon))*DqDt(i);
end

% Concatenate vector of time derivatives
% Very small imaginary components need to be removed!
DyDt = [real(DcDt); real(DqDt)];
end

[Alan Stevens]

Ok, but there are a couple of apparent anomalies in your equations. I've typed them out to ensure I understand them - see the attachment. I've included questions which need to be answered before I spend any time implementing them.

[nashwa tarek]

thank alot Alan

[Alan Stevens]

See attached. Am I getting any closer to what you want? Check before I try to implement!

Also, I don't have data for particle diameter.

[nashwa tarek]

thanks alot Alan

i hope this can work, and the breakthrough curve must be s shape like the one i attched

ireally cant thank you more really appricated

[Alan Stevens]

Hmm! You have now supplied a value (for ap) with units. All the other constants need to have units too, so that we can be sure the equations are all consistent. For example, you give L = 1, without units. I suspect this isn't meant to be 1 mm! Please provide units so I can check for consistency.

[nashwa tarek]

let's say L=0.1m , ap= 0.002m by this they will be consistant with other variables with Kg, m and sec

[Alan Stevens]

Do you have an equation for Ce?

[nashwa tarek]

no but we can assume Ce=Cs

[Alan Stevens]

Hmm. If you look carefully at your last hand written equation (number 9), you will see that it implies dCs/dt is always zero. I've therefore written a program using my previous equations - see attached. In order to generate sensible looking values for C and Cs and q I had to change the sign of the right hand side of the dCsdt equation.

Actually, I don't know if the results are sensible or not! You will have to be the judge of that.

Notice that the attached program produces four graphs: one for the exit values of Cb/Cfeed vs t, one for Cs vs t, one for the difference between Cb and Cs vs time and one for q vs time.

I've probably helped (or hindered!) as much as I can here. Good luck!

[nashwa tarek]

thank you Alan so much really appriated

[nashwa tarek]

thanks alot Alan you helped me alot, really cant thank you more

i have one last question, when i change the value of Cfeed the curve of cb/cfeed vs time should change but it doesnt? do you have any clue

[Alan Stevens]

The curve of Cb should (and does) change, but the equations are linear in C, so the normalised curve (i.e. Cb/Cfeed) should be (and is) independent of Cfeed. As Cfeed goes up(or down) so to does C, proportionately.

The curve of Cb/Cfeed should be the same for all values of Cfeed as long as nothing else changes. (Cb itself does change, of course.)

[Rohan Singh]

@Alan Stevens How will the code change for equation of Cm change if spherical coordinates are used for particle continuity equations ?

[Alan Stevens]

Cm?

Have a look at http://pleasemakeanote.blogspot.com/2010/02/9-derivation-of-continuity-equation-in.html for continuity equation in spherical coordinates.

[Rohan Singh]

Thanks a lot @Alan Stevens

I am doing something similar to nashwa but second equation is also radial variations to it (equations attached). I tried to modify the code but I am not getting a proper curve. Can you please help me out with this ?

Cfeed = 15; % Inlet concentration

tf = 8640000; % End time

nz = 200; % Number of gateways

np = nz + 1; % Number of gateposts

% 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/10000;

tspan = 0:dt:tf;

% Call ode solver

[t, c] = ode15s(@rates, tspan, c0);

% Plot results

hold on

plot(t, c(:,np)/Cfeed),grid

xlabel('time'), ylabel('Cb/Cfeed')

title('Figure:F(Treated laterite) 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')

%k = 2.5*10^-5; nf = 1.45;

%q = k*c(:,2*np).^(1/nf); % equation

% figure

% plot(t,q), grid

% xlabel('time'), ylabel('q')

% title('Figure 4: Exit q vs time')

% Function to calculate rate of change of concentrations with time

function dcdt = rates(~, c)

Dz = 4.3*10^-9; % Diffusion coefficient

u = 0.001; % Velocity

e = 0.4;

Kf = 2.5*10^-5;

ap = 0.9*10^-3; % Particle diameter

ep = 0.3;

L = 1; % Length

nz = 200; % Number of gateways

np = nz + 1; % Number of gateposts

dz = L/nz;

R=0.00001;

dr=R/nz;

rho = 480;

qm = 118;

b= 0.03;

Dp = 1.2*10^-10;

C = c(1:np);

Cs = c(np+1:end);

dCdt = zeros(np,1);

dCsdt = zeros(np,1);

%dCsdr = zeros(np,1);

%dCsdr2 = zeros(np,1);

for i = 2:np-1

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)*3*Kf/(e*ap)*(C(i)-Cs(i));

%dCsdt(i) = 3*Kf/(ep*ap)*(C(i)-Cs(i)); % Needed to change the sign from -ve to +ve here.

%dCsdr =

dCsdt(i) = (1/(1+(rho*((1-ep)/ep)*(qm*b)/((1+b*Cs(i))^2)))*Dp*((Cs(i+1)-2*Cs(i)+Cs(i-1))/(dr^2)+(2/R*(Cs(i+1)-Cs(i-1))/(2*dr))));

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

[Alan Stevens]

Hmm. From the equations in your .jpg files it looks like Cp is independent of Z, so can be solved as a function of time and R (you don’t have R varying anywhere as far as I can see!).

Cb seems to be independent of R, except at R=Rp (wherever that is – the outer surface perhaps?).

I suggest you try to solve Cp as a function of t and R first. Then, completely separately, solve Cb as a function of t and Z, using the values of Cp previously obtained at t and Rp.

For Cp you need initial conditions for Cp and dCp/dR (the latter is probably zero); and for Cb you need initial conditions for Cb and dCb/dt.

[Rohan Singh]

@Alan Stevens Thank you so much

[Rohan Singh]

Sorry to bother you again @Alan Stevens. I am still learning this stuff. I tried to implement your suggested way but for some reason, i am getting a strange curve. It would really help me if you can guide me further what to do. Thanks in advance.(I have attached the equations and boundary conditions for your reference)

Cfeed = 10; % Inlet concentration

tf = 86400; % End time

nz = 200; % Number of gateways

np = nz + 1; % Number of gateposts

% 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/10000;

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')

%

% k = 2.5*10^-5; nf = 1.45;

% q = k*c(:,2*np).^(1/nf); % Freundlcih equation

% figure

% plot(t,q), grid

% xlabel('time'), ylabel('q')

% title('Figure 4: Exit q vs time')

% Function to calculate rate of change of concentrations with time

function dcdt = rates(~, c)

Dz = 3*10^-8; % Diffusion coefficient

u = 2*10^-5; % Velocity

e = 0.4;

Kf = 3*10^-5;

ap = 2*10^-3; % Particle diameter

ep = 0.53;

L = 0.1; % Length

nz = 200; % Number of gateways

np = nz + 1; % Number of gateposts

dz = L/nz;

dr = ap/nz;

rho = 400;

Dp = 1.2*10^-10;

b=0.03;

qm=118;

C = c(1:np);

Cs = c(np+1:end);

Cs(1) = 0;

dCdt = zeros(np,1);

dCsdt = zeros(np,1);

dCsdr = zeros(np,1);

% rhalf = zeros(np-1,1);

% DCsDr = zeros(np,1);

D2CsDr2 = zeros(np,1);

%

% rhalf(1:np-1)=(z(1:np-1)+z(2:np))/2;

for i = 2:np-1

dCsdr(i) = (2*dr)*(Cs(i+1)-Cs(i-1));

D2CsDr2(i) = (1/dr^2)*(Cs(i+1)-2*Cs(i)+Cs(i-1));

% dCsdt(i) = 3*Kf/(ep*ap)*(C(i)-Cs(i)); % Needed to change the sign from -ve to +ve here.

end % i loop

for i = 2:np-1

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)*3*Kf/(e*ap)*(C(i)-Cs(i));

% dCsdt(i) = 3*Kf/(ep*ap)*(C(i)-Cs(i)); % Needed to change the sign from -ve to +ve here.

dCsdt(i) = 1/(1+rho*((1-ep)/ep)*(qm*b/((1+b*Cs(i))^2)))*Dp*(D2CsDr2(i)+((2/ap)*dCsdr(i)));

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

Answer 2

[Anthony Ozzello]

Hi @Alan Stevens,

I'm sorry to follow up on an old thread. I'm a relatively new user to Matlab and I'm trying to model a fixed bed reactor like the folks here, but unlike them, I'd like to be able to visualize the values of C and q throughout the depth of the bed as a function of t, so I need the function to return C and q as a matrix of both t and z, I believe. Is there a simple way to do this? I originally thought I would have to use the pdesolver, but this is an interesting way to handle this. Please let me know if there is a simple way to look at the results as both a function of t and z so that I can do cumtraps and surface plots throughout the depth of the fixed bed. Thanks, Tony

[Alan Stevens]

Look at the Y output from:

[T, Y] = ode15s(@(t,y) MyFun(t,y,z,n),t,y0);

The rows of Y represent values at times given by T.

The columns of Y represent values of c and q at heights, z. The columns from 1 to n are values of c, while the columns from n+1 to 2n are values of q.

So, for example, plotting

plot(z,Y(100,1:n))

will show the values of c as a function of z at time T(100); and

plot(z,Y(150,n+1:2*n))

will show the values of q as a function of z at time T(150).

[Anthony Ozzello]

Thanks @Alan Stevens,

I thought there was something like this to it, I also figured out how to do the surface plot by extracting the values into matrices, maybe there's a simpler way.

Another question for you too. I don't know the values for Kc and D in the function call, but I should be able to fit them from some known effluent curves from existing column data. I would have to wrap all this into some kind of fitting routine with a least squares or something like that against my known effluent curve. Can you point me to a good reference to take the output of the ode and do some minimization with it?

Finally, have you ever looked at multicomponent adsorption in a fixed bed like this? I assume the code would be pretty similar?

Thanks,

Tony

[Alan Stevens]

For fitting two unknowns:

help fminsearch

I've never looked at adsorption of any kind in a fixed bed like this! I simply looked at the equations the original poster submitted!

[Anthony Ozzello]

Thanks Alan,

I've gotten everything working, will look into fminsearch. Really appreciate the help.

Tony

[Anthony Ozzello]

Hi @Alan Stevens,

There is something interesting in the code that I need your help understanding. Likely due to me being a newer user to MATLAB.

In this code, we define the time vector initially in the call to ode15s like so:

[T, Y] = ode15s(@(t,y) Myfun(t,y,z,n),[t0 tf],y0);

I do note that in the function definition that it looks like the time vector is ignored and maybe it only cares about t0 abd tf?

function DyDt=Myfun(~,y,z,n)

But the return variable 'T' is very different from 't'. It is much larger, but it never gets modified in my function.

I have a set of experimental data points for (t,y) that I want to fit to some constants that go into Myfun and I have a pretty good idea how to use fminsearch to do that, but the return of [T Y] is on a different set of time values than my original data so I don't have a simple way to compare my experimental and theoretical data.

For example, I have 12 experimental data points covering a range of 380 minutes, with intervals of 10 minutes. so I'd like to set my time vector to 380 minutes with a delta t of 10, but 'T' comes back with a lot more values than this at different intervals. What to do?

Thanks,

Tony

[Torsten]

tspan = 0:10:380;
[T, Y] = ode15s(@(t,y) Myfun(t,y,z,n),tspan,y0);

[Anthony Ozzello]

That got it, thanks!

[Anthony Ozzello]

Hey again,

I'm working on fitting some of the parameters in the ODE using the problem-based approach outline in the MATLAB documentation here: https://www.mathworks.com/help/optim/ug/fit-ode-problem-based-least-squares.html. I've got everything working properly except for plotting out the results of the fitted data points against the original data. My data exists in time space at non-evenly spaced intervals from 5 to 380 minutes and to get a smooth function, I'm integrating every 0.1 minute. I modified the RtoODE function to only return the solution points that match the original dataset but using this code, where exp_t contains the timestamps of the original data and tspan(2) is the time interval in the integration. There can be multiple components in the analysis all concatented together into y.

function solpts = RtoODE(r,z,n,nComponents,tspan,y0, exp_t)

[T, Y] = ode15s(@(t,y) Myfun(t,y,z,n,nComponents,r),tspan,y0);

for i = 1:length(exp_t)

time = exp_t(i);

j= time/tspan(2)+1;

for k = 1: nComponents

solpts(i,k)=Y(j,k*n);

end

end

end

But when I got to plot the data that comes back, I get an error, where exp_t and exp_y are my raw data, there are 2 components for y, but I'm only trying to plot the first one right now.

myfcn = fcn2optimexpr(@RtoODE,r,z,n,nComponents,tspan,y0, exp_t);

figure

plot (exp_t,exp_y(:,1),'g')

hold on

plot(exp_t,myfcn(:,1),'b')

hold off

Error using plot

Data must be numeric, datetime, duration or an array convertible to double.

It should work based on the example in the documentation, but it doesn't and I can't find anything in any of the previous question and answers or anything elsewhere online that gives me some ideas on how to plot optimization variables if this isn't it.

Any clues?

Thanks,

Tony

[Anthony Ozzello]

nvm, I found the error.

Tony

Answer 3

[Anthony Ozzello]

Hi yet again,

I followed Torsten's response, which was very helpful to setting the time span, which enabled me to compare the results of the calculated integration to my experimental data points very well, but I encountered a problem. It seems that my problem requires significant calculation at a lot of points near the origin and then can space out away from it. When I specify the tspan to ensure I capture my experimental points, I have to use more of a linspace to ensure I have the right timestamps and, even if I go with pretty small time intervals, it's nowhere near tight enough near the origin and I'm getting negative concentrations. So it's clear that I have to let MATLAB chose the timespan properly, I think, by just specifying the end points. But then my question is how do I do some kind of fuzzy matching between the timestamps of my experimental data and the timestamps that come back from ode15s for calculating the least squares fit? Can someone point me in a direction or some examples? Thanks, Tony

[Torsten]

The time steps don't matter, it's the number of discretization points you must enlarge (number of elements of the z-vector).

[Anthony Ozzello]

Thanks Torsten, I will try that.

[Anthony Ozzello]

Hi Torsten, I changed the z discretization by a factor of 5 and, in one case it didn't change the t-value at which the concentration went negative and in another case it pushed it a little more to the right. I did see how the time steps don't matter though. It will stay positive for a variety of Langmuir feed parameters, but never does for any Freundlich ones, once the exponent for the Freundlich isotherm deviates from 1, the concentration goes negative and then things are a mess.

[Torsten]

If you include the code, I'll have a look at it.

[Anthony Ozzello]

Hi Torsten, thanks for the offer. I'm attaching two files, the MATLAB script and the .xlsx file with the experimental data I am reading in. If you run it, as it, it works fine. If you change the second curve from Langmuir to Freundlich in cell A3 of the spreadsheet, that's where we start running into problems. Thanks, Tony

[Anthony Ozzello]

you can also change the Langmuir constants outside a reasonable range and it'll also give negative concentration values in the middle of the calculation, but it doesn't fail and it gives reasonable effluent values. You can try a combination like (10, 0.4). Tony

[Torsten]

Changing the line

DcDz(i+k*n) = ((z(i)-z(i-1))/(z(i+1)-z(i))*(c(i+1+k*n)-c(i+k*n))+(z(i+1)-z(i))/(z(i)-z(i-1))*(c(i+k*n)-c(i-1+k*n)))/(z(i+1)-z(i-1));

to

DcDz(i+k*n) = (c(i+k*n)-c(i-1+k*n))/(z(i)-z(i-1))

together with more elements for the z-vector to avoid numerical diffusion usually helps.

[Anthony Ozzello]

I made the change and the concentration still goes complex with z=100 steps and it stops with stepssize smaller than minimum about 2.3min. I changed to z=200 steps and the time where it stops changes to 1.76mins. I changed it to 250 steps and it's still going to complex concentrations, but it is taking forever to run now and it seems to be hung at t=1.68mins. I think it's going in the wrong direction as I increase the amount of elements in the z-vector. Tony

[Anthony Ozzello]

Hi Torsten, I've got it working much better but I cannot escape the negative concentration and complex results with the Freundlich isotherm. I've increased the z discretization from 100 up to 3000 and it still does not resolve the issue and it ends up stopping after shorter and shorter times as I decrease these step sizes. Is it possible that the method of lines is just not compatible with the Freundlich isotherm?

[Torsten]

I think in your case concentrations in the gas phase get negative because the driving force between gaseous and adsorbed phase is very strong. And as soon as one concentration in the gas phase becomes negative, the ^1/n in the Freundlich isotherm gives complex results for equilibrium q*.

Try to use smaller RelTol and/or AbsTol for ODE15s.

Try to replace c(...)^1/n in the Freundlich isotherm by abs(c(...))^1/n.

[Anthony Ozzello]

Hi Torsten, I did all these things and I'm still getting the issue of a complex result. I took some data and fitted it to both a Langmuir and a Freundlich isotherm so I knew I would have good constants and the Freundlich still goes complex on me. I tried stepping through the iterations of the loop as it works through the packed bed step by step and I found something interesting and I'm not sure if it makes sense or not. Every time the script calls the function to integrate it, I can watch a small pulse of concentration moving from z=0 to z=L at t=0, it's only one z-slice wide and it always has the same value for q and c. All the other values for q and c in the vector are zero except for the value at z=0, which is the feed. It doesn't make sense to me that this single pulse of concentration is existing through the bed at t=0 with each iteration of the function call. This happens for both Langmuir and Freundlich isotherms, is it an artifact or expected behavior? If I use the absolute value of the concentration in the isotherm calculation, it gives me complete garbage.

[Torsten]

Maybe you could post a ready-to-use version only of the reactor model with the Freundlich isotherm (not coupled with optimization) where it is possible for me to reproduce the problem you encounter.

[Anthony Ozzello]

Hi Torsten, Here is the code, stripped of the optimization. Let me know what you might find. thanks, Tony

[Torsten]

I changed a few things in the output section, but I can't reproduce the behaviour you describe concerning the Freundlich isotherm.

[Anthony Ozzello]

Hi Torsten,

Thanks for looking at it so rapidly.

You have to change the commenting in two places, swap line 198 with line 202 as well as swapping line 184 with line 188. Then you get the issue.

Tony

[Torsten]

This code works, but I'm not sure if this is a stable solution in every instance.

[Anthony Ozzello]

Hi Torsten, I see how the code completes, but the concentrations and masses on the bed are still complex. It hasn't resolved that issue. I'm hoping to find a solution that doesn't result in complex concentrations. I'm not sure the math is trustworthy if it's giving me complex solutions. I don't know why the concentration is going negative in the first place which is driving all this. With the Langmuir isotherm, the concentration never drops below zero. So we have to be driving dq/dt into a weird space that isn't allowing c to assympote to zero from above. But I can't figure out how this is happening. Do you have any more ideas?

[Torsten]

I must admit that I used "Octave" to run your problem, and the ode15s solver therein did not give negative or complex concentrations with the settings I've done.

Since I can't test the MATLAB ode15s solver, I think I can't give you further advice.

[Anthony Ozzello]

ahh, bummer. so can you run a quick check on you solution while it's running an verify that c(i), for all i, never goes negative when you use the Freundlich isotherm? Is that's true, then there must be something I need to do differently with the MATLAB ode15s setup I would assume. Thanks, Tony

[Torsten]

Negative concentrations appear in the solution for c, but with smallest value at about -5e-9.

This will influence neither convergence nor reliability of the results in my opinion.

[Anthony Ozzello]

Hi Torsten, I agree that it's a really small number, but it's what drives the solution into complex numbers since the isotherm is qe=qs*c^(1/n) and n is not 1, I'm not sure why you are not getting complex numbers? It still converges for me as well, but the complex numbers are present.

[Torsten]

For negative concentrations c, DyDt becomes partially complex, but the ode15s integrator in Octave only accounts for the real part. The imaginary part is ignored.

Maybe you could test my settings by inserting the line

DyDt = real(DyDt)

at the end of the function where you define the ODEs.

[Anthony Ozzello]

Hi Torsten, I tried this as well as putting the real statement right at the isotherm calculation and now the integral comes back with an error that it cannot complete at the time step required after more than an hour of processing where before it would complete in about 0.2 seconds. This won't work. Any idea on how to set ode15s or any other stiff solver in Matlab to function that way yours is in Octave? Thanks, Tony

[Torsten]

You could use the octave solver under MATLAB or try a different stiff solver, e.g. RADAU5, for which a MATLAB version exists, from

http://www.unige.ch/~hairer/software.html

All other options for the MATLAB - ODE15S are trial-and-error and had to be tested under the MATLAB environment. But I think RelTol, AbsTol and BDF-order are the relevant settings from which you already reported as not being successful.

[Anthony Ozzello]

Hi Torsten, I may look into that. I can get it to complete for a limited range of values, and if I ignore the complex portion, I get reasonable curves, so maybe I'm okay? The optimization always fails due to it stepping outside the range of acceptable values, maybe I can deal with that using constraints. I'll try to figure out how to use other solvers under MATLAB, there's probably some posts I can find. Thanks for all your help, it's been invaluable. Tony

[Anthony Ozzello]

Hi Torsten or anyone, coming back to this topic. I have a question about how the second derivative is calculated in this code. I'm searching literature and having a hard time finding a reference to the approach used here which is looking like (z*dc/dz(forwards)-z*dc/dz(backwards))/dz. Why is this used vs central differences for the second derivative? Are there advantages? If I change to central differences, I get different results, so I would like to get a feel for which is correct. Thanks, Tony:

[Torsten]

The approximation of the second derivative used above is incorrect.

The correct discretization is

D2cDz2(i) = ((c(i+1)-c(i))/(z(i+1)-z(i))-(c(i)-c(i-1))/(z(i)-z(i-1)))/(zhalf(i)-zhalf(i-1));

[Anthony Ozzello]

thanks Torsten

[Anthony Ozzello]

How would that change the boundary condition: D2cDz2(n) = -1.0/(z(n)-zhalf(n-1))*(c(n)-c(n-1))/(z(n)-z(n-1)); then?

[Torsten]

This is the correct approximation of the second derivative at the boundary if the boundary condition is set to dc/dz = 0.

d^2c/dz^2 
~ (dc/dz(n)-dc/dz(n-1/2))/(z(n)-z(n-1/2)) 
= -dc/dz(n-1/2)/(z(n)-z(n-1/2)) (assuming dc/dz(n) = 0)
~ -[(c(n)-c(n-1))/(z(n)-z(n-1))]/(z(n)-z(n-1/2))

where

z(n-1/2)     = zhalf(n-1), 
dc/dz(n-1/2) = dc/dz at zhalf(n-1)

[Anthony Ozzello]

Thanks again Torsten

[Federico Urbinati]

@Anthony Ozzello

Hi Anthony, have you completed the code related to the fixed bed adsorption of multicomponents?

i'm really stuck, you would do me a favor if i posted it here.

thanks in advance,

Federico

Answer 4

[Puteri Ellyana Mohmad Latfi]

Hi @Alan Stevens and @Torsten,

I am not really good in Matlab but I have to model a fixed bed adsorption to obtain the breakthrough curve of the adsorption. I tried to incorporate my model in Alan's code for Nashwa but I am not really sure if the initial and boundary condition in the code are correct though. Beside that, I got this error when i ran this code.

Warning: Failure at t=2.006182e+02. Unable to meet integration tolerances without reducing the step size below the

smallest value allowed (4.547474e-13) at time t.

> In ode15s (line 655)

In fyp (line 19)

I would really appreciate any of your help. Thank you very much. I have also attached the model's equations in the file.

Cfeed = 0.0224;     % Inlet concentration

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

% Initial Concentrations
Cb = zeros(np,1); Cb(1) = Cfeed;  
Cp = zeros(np,1);

% Prepare data to pass to ode solver
c0 = [Cb; Cp];   
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('Cp')
title('Figure 2: Exit Cp vs time')
figure
plot(t, c(:,np)-c(:,2*np)), grid
xlabel('time'), ylabel('Cb - Cp')
title('Figure 3: Exit Cb-Cp vs time')

k = 598;
q = k*c(:,2*np);  % Freundlcih equation
figure
plot(t,q), grid
xlabel('time'), ylabel('q')
title('Figure 4: Exit q vs time')

% Function to calculate rate of change of concentrations with time
function dcdt = rates(~, c)
    a=0.00000000261;
    b=0.2102;
    d=1705.4719;
    e=-38.8576;
    L=0.08;
    nz = 200;       % Number of gateways
    np = nz + 1;    % Number of gateposts
    dz = L/nz;

    Cb = c(1:np);
    Cp = c(np+1:end);
    dCbdt = zeros(np,1);
    dCpdt = zeros(np,1);
    
   for i = 2:np-1
       dCbdt(i) = (a/dz^2)*(Cb(i+1)-2*Cb(i)+Cb(i-1)) - b/(2*dz)*(Cb(i+1)-Cb(i-1)) - d*(Cb(i)-Cp(i));
       dCpdt(i) = e*(Cb(i)-Cp(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
    dCbdt(np) = dCbdt(np-1); 
    dCpdt(np) = dCpdt(np-1);
    
    % Combine rates for the function to return to the calling routine
    dcdt = [dCbdt; dCpdt];

end % of rates function 

[Alan Stevens]

By making e positive it works. However, I don't know if the results are sensible or not!

Cfeed = 0.0224;     % Inlet concentration

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

% Initial Concentrations
Cb = zeros(np,1); Cb(1) = Cfeed;  
Cp = zeros(np,1);

% Prepare data to pass to ode solver
c0 = [Cb; Cp];   
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('Cp')
title('Figure 2: Exit Cp vs time')

figure
plot(t, c(:,np)-c(:,2*np)), grid
xlabel('time'), ylabel('Cb - Cp')
title('Figure 3: Exit Cb-Cp vs time')

k = 598;
q = k*c(:,2*np);  % Freundlcih equation
figure
plot(t,q), grid
xlabel('time'), ylabel('q')
title('Figure 4: Exit q vs time')

% Function to calculate rate of change of concentrations with time
function dcdt = rates(~, c)
    a=0.00000000261;
    b=0.2102;
    d=1705.4719;
    e=38.8576;      % By making e positive this works.  
    L=0.08;
    nz = 200;       % Number of gateways
    np = nz + 1;    % Number of gateposts
    dz = L/nz;

    Cb = c(1:np);
    Cp = c(np+1:end);
    dCbdt = zeros(np,1);
    dCpdt = zeros(np,1);
    
   for i = 2:np-1
       dCbdt(i) = (a/dz^2)*(Cb(i+1)-2*Cb(i)+Cb(i-1)) - b/(2*dz)*(Cb(i+1)-Cb(i-1)) - d*(Cb(i)-Cp(i));
       dCpdt(i) = e*(Cb(i)-Cp(i));  
   end % i loop
   
   % For simplicity, set exit rates to those of the immediately previous nodes
    dCbdt(np) = dCbdt(np-1); 
    dCpdt(np) = dCpdt(np-1);
    
    % Combine rates for the function to return to the calling routine
    dcdt = [dCbdt; dCpdt];

end % of rates function 

[Puteri Ellyana Mohmad Latfi]

I see. Thanks Alan. Is the initial and boundary conditions in the code correspond to the ones in the file attached? and how do i edit the code to get a nicer S-shaped graph for the Cb/Cfeed vs t graph?

[Alan Stevens]

"Is the initial and boundary conditions in the code correspond to the ones in the file attached?"

I don't know - I'll leave that for you to check.

"how do i edit the code to get a nicer S-shaped graph for the Cb/Cfeed vs t graph"

Set tf = 30, say, rather than 2000.

[Puteri Ellyana Mohmad Latfi]

Hi @Alan Stevens. I tried using the coding for different set of experimental values, however, i didnt get the same graph as the experimental breakthrough graph, which is shown in the image attached (10000ppm).

Cfeed = 0.0227;     % Inlet concentration

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

% Initial Concentrations
Cb = zeros(np,1); Cb(1) = Cfeed;  
Cp = zeros(np,1);

% Prepare data to pass to ode solver
c0 = [Cb; Cp];   
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('Cp')
% title('Figure 2: Exit Cp vs time')
% 
% figure
% plot(t, c(:,np)-c(:,2*np)), grid
% xlabel('time'), ylabel('Cb - Cp')
% title('Figure 3: Exit Cb-Cp vs time')
% 
% k = 56.5;
% q = k*c(:,2*np);  % Freundlcih equation
% figure
% plot(t,q), grid
% xlabel('time'), ylabel('q')
% title('Figure 4: Exit q vs time')

% Function to calculate rate of change of concentrations with time
function dcdt = rates(~, c)
    a=0.00000001278;
    b=0.160556;
    d=244.942875;
    e=-1.3964;
    L=0.065;
    
    nz = 200;       % Number of gateways
    np = nz + 1;    % Number of gateposts
    dz = L/nz;
    Cb = c(1:np);
    Cp = c(np+1:end);
    dCbdt = zeros(np,1);
    dCpdt = zeros(np,1);
    
   for i = 2:np-1
       dCbdt(i) = (a/dz^2)*(Cb(i+1)-2*Cb(i)+Cb(i-1)) - b/(dz)*(Cb(i+1)-Cb(i-1)) - d*(Cb(i)-Cp(i));
       dCpdt(i) = -e*(Cb(i)-Cp(i)); 
   end % i loop
   
   % For simplicity, set exit rates to those of the immediately previous nodes
    dCbdt(np) = dCbdt(np-1); 
    dCpdt(np) = dCpdt(np-1);
    
    % Combine rates for the function to return to the calling routine
    dcdt = [dCbdt; dCpdt];
end % of rates function 

[Alan Stevens]

All I can suggest is check your data carefully. I don't know enough about your system to connect the data values in your graph.png to the constants in the program.