fixed-bed adsorption column-modeling solving PDE

hi all i am trying obtainbreakthrough curve from that model (c/co vs time) i want to study the effect of initial concentration, flow rate and column length. my problem now that by changing the initial concentration ther is no change in the curve, i am not good at matlab i just want to use it to finish my thesis. can anyone help me please to figure out the problem

function main

%System Set-up %

%Define Variables

%D = 1.85*10^-5; % Axial Dispersion coefficient

%v = 3.38*10^-3; % Superficial velocity

%epsilon = 0.47; % Voidage fraction

%k = 1.9*10^-6; % Mass Transfer Coefficient

%ap= 0.003; % radius of pellet

%b = 0.05; % Langmuir parameter

%qm = 14.5; % saturation capacity

% rho= 1970; % paricle denisty

cFeed = 20; % Feed concentration

L = 0.02; % Column length

t0 = 0; % Initial Time

%tf = 20; % Final time

tf = 20; % 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(:,n)/cFeed)

end

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

% Defining Constants

D = 1.85*10^-5; % Axial Dispersion coefficient

v =3.38*10^-3 ; % Superficial velocity

epsilon = 0.47; % Voidage fraction

k = 1.9*10^-6; % Mass Transfer Coefficient

b = 0.05; % Langmuir parameter

qm = 14.5; % saturation capacity

rho= 1970; %particle denisty

% 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*( q(1)- ((qm*b*c(1))/(1 + b*c(1))) );

DcDt(1) = 0.0;

% Set time derivatives in remaining points

for i=2:n

%Equation: dq/dt = k(q-q*) where q* = qm * (b*c)/(1+b*c)

DqDt(i) = k*( q(i)- ((qm*b*c(i))/(1 + b*c(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))*rho *DqDt(i);

end

% Concatenate vector of time derivatives

DyDt = [DcDt;DqDt];

end

[Sona]

have you managed to solve it? I am facing the same problem for my thesis ^^

Answer 1

[Sulaymon Eshkabilov]

Here is an edited and corrected version of your code (smaller time step 'cause it is a stiff DE) that shows that the concentration change makes an impact on the simulation results:

L = 0.02;         % Column length
t0 = 0;        % Initial Time
tf = 20;      % Final time
dt = 0.05;      % Time step        
z = 0:0.0005:L; %Mesh generation
t = t0:dt:tf;% Time vector
n = numel(z);  % Size of mesh grid

cF = [0.1 1.5 5 10 50 100 250 500 1000];   % Diffferent feed concentration values
for ii=1:numel(cF)
    cFeed=cF(ii);
    [T, Y] = Main(cFeed);
    plot(T,Y(:,n), 'linewidth', 2), hold all
end

function [T,Y]=Main(cFeed)
%System Set-up %
%Define Variables
%D = 1.85*10^-5;    % Axial Dispersion coefficient
%v = 3.38*10^-3;    % Superficial velocity
%epsilon = 0.47;  % Voidage fraction
%k = 1.9*10^-6;    % Mass Transfer Coefficient
%ap= 0.003;    % radius of pellet
%b = 0.05;  % Langmuir parameter
%qm = 14.5; %  saturation capacity
 
% rho= 1970;  % paricle denisty
%cFeed = 10;     % Feed concentration
L = 0.02;         % Column length
t0 = 0;        % Initial Time
%tf = 20;      % Final time
tf = 20;      % Final time
dt = 0.05;      % 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);
 
end 
function DyDt=MyFun(~, y, z, n)
% Defining Constants
D = 1.85*10^-5;    % Axial Dispersion coefficient
v =3.38*10^-3 ;    % Superficial velocity
epsilon = 0.47;  % Voidage fraction
k = 1.9*10^-6;    % Mass Transfer Coefficient
b = 0.05;  % Langmuir parameter
qm = 14.5; % saturation capacity
 
rho= 1970;   %particle denisty
% 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) = (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*( q(1)- ((qm*b*c(1))/(1 + b*c(1))) );
DcDt(1) = 0.0;
% Set time derivatives in remaining points
for i=2:n
    %Equation: dq/dt = k(q-q*) where q* = qm * (b*c)/(1+b*c)
    DqDt(i) = k*( q(i)- ((qm*b*c(i))/(1 + b*c(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))*rho *DqDt(i);
  end
% Concatenate vector of time derivatives
DyDt = [DcDt;DqDt];
end

[Torsten]

The discretization of the second-order derivatives D2cDz2 in the inner grid points is incorrect.

[Sona]

How is it supposed to be?

[Torsten]

In cartesian coordinates:

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

In cylindrical or spherical coordinates, it will be different, of course.