Adsorption Modelling - Solving PDE - Axial Dispersion Model
54 views (last 30 days)
Show older comments
Ashley Gratton
on 1 Mar 2018
Commented: Federico Urbinati
on 19 Oct 2022
Dear all,
I am trying to develop a simple model for an adsorption column. I have attached a file that contains the equations that I am trying to solve.
To solve the system, I am trying to solve the axial dispersion model (neglecting pressure drop, and hence, the dv/dz term) alongside the Linear Driving Force (LDF) model, which gives an expression for dq/dt in the axial dispersion PDE (see attached file). I am also using the Langmuir equation to give an expression for q* (equilibrium concentration) in the LDF equation.
Unknown variables are q and c, all other variables are known (i.e. D, v, epsilon etc.( / are the independent variables (i.e. time, t, and length, z).
I have tried using the Method of Lines to solve this system (again, see attached file), but I am aware that the system I have developed is flawed.
Any suggestions on how to approach this problem would be appreciated. I would preferably like to solve the system using either a PDE solver or ode45 applied to the method of lines.
My code so far is as follows:
Data File:
%%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
bH = 2.5*10^-5; % Langmuir parameter for H2
L = 1; % Column length
%%MoL Coefficients & Mesh Generation
N = 20;
dz = L/N;
z = [0:dz:L];
% MoL Coefficients from Discretisation
alph = ((D/(dz^2)) + ((v)/(2*dz)));
bet = ((-2*D)/(dz^2));
gamm = ((D/(dz^2)) - ((v)/(2*dz)));
%%Initial / Boundary Conditions
% Initial condition for ODE 1:
c0 = 0; % C at t = 0 for all z
% Boundary Condition for ODE 1:
cFeed = 10; % C at z = 0 for all t
qs = (1+bH*cFeed)/(bH*cFeed); % Conc. of adsorbed phase at surface saturation
% Initial Condition for ODE 2:
q0 = qs*((bH*cFeed)/(1+bH*cFeed))-1; % Initial concentration in adsorbed phase
% Time Boundaries
t0 = 0;
tf = 10000;
% Initial equilibrium conc.
qeq0 = qs*((bH*c0)/(1+bH*c0));
dq0 = k*qeq0;
Solving using ODE45:
close all
clear all
clc
% Call Data File
DataFile;
tspan = [0:10000];
[t,c] = ode45('cprime', tspan, c0);
Cprime file
function [ dcdt ] = cprime(~, c)
% Call data file to define constants
DataFile;
t=t0;
q = q0;
[dqdt] = qprime(t, q, k, qs, bH, c);
%%Create Tri-diagonal Matrix A
A = zeros(N,N);
for i = 2 : N-1
A(i, i-1) = alph;
A(i, i) = bet;
A(i, i+1) = gamm;
end
A(1,1) = bet;
A(1,2) = gamm;
A(N,N-1) = alph+gamm;
A(N-1, N) = gamm;
A(N,N) = bet;
b = ones(N,1);
b(1) = cFeed*alph-((1-epsilon)/(epsilon)*(k*(qeq-q)));
dcdt = (A*c + b);
end
qprime file:
function [dqdt] = qprime(~, q, k, qs, bH, c)
DataFile;
qeq = qs*((bH*c)/(1+bH*c));
dqdt = k*(qeq-q);
end
I am really struggling to figure out how to structure my code / approach these equations.
Any help would be muchly appreciated!
Thanks!
46 Comments
Ridhwan Lawal
on 3 May 2022
Hello @Ashley Gratton please Did you find a solution to your problem? I am havin this exact same challenge with my project too
Federico Urbinati
on 19 Oct 2022
hi halok
did you manage to define the matlab code for multiadsorbate system?
I can't go ahead, can you publish it?
if you want I can give you a contact to write to us in private
Thanks in advance
federico
Accepted Answer
Torsten
on 2 Mar 2018
The code I provided here for a very similar problem might help:
https://de.mathworks.com/matlabcentral/answers/371313-error-in-solving-system-of-two-reaction-diffusion-equations
Best wishes
Torsten.
20 Comments
Yash Toshniwal
on 26 Jun 2019
function main
% Initialisation
close all
clear all
clc
%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
b = 2.5*10^-5; % Langmuir parameter
qs = 2*10^4; % Saturation capacity
cFeed = 10; % Feed concentration
L = 1; % Column length
t0 = 0; % Initial Time
tf = 1000; % Final time
dt = 0.5; % Time step
z = [0:0.01: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);
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
b = 2.5*10^-5; % Langmuir parameter
qs = 20000; % Saturation capacity
% 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));
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*( ((qs*b*c(1))/(1 + b*c(1))) - q(1) );
DcDt(1) = 0.0;
% Set time derivatives in remaining points
for i=2:n
%Equation: dq/dt = k(q*-q) where q* = qs * (b*c)/(1+b*c)
DqDt(i) = k*( ((qs*b*c(i))/(1 + b*c(i))) - 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
When i try to run this code in matlab, it does not produce any result, my problem is quite similar with the above code, if anyone can provide me with the working code it will be very helpful.
Zain Shabeeb
on 9 Jul 2019
Hi Ashley, Torsten,
I have been trying to model a system that is similar to this one. So far I have had trouble understanding the system, but after analysing it closely, I realized that it is very similar to your system, except that there is an extra term in the equation to account for transport under electric field. Basically this is a coupled ion exchange - electrodialysis system that I am trying to model. I have attached the model equations in the first snapshot. In the second snapshot I have rearranged the equations to match the model that you have provoided in this post. As you can see, the difference is the extra term containing phi (electric field).
Your model is with respect to time and one direction. In my model, the bulk concentration is varying with time and the z direction, and the solid concentration is varying with time and the x direction. The electric potential gradient is in the x direction only.
Please let me know what changes to make to your MATLAB code so that I can model my system. Is it even possible to use this kind of code to model my system?
Thanks in advance.
More Answers (7)
Abraham Boayue
on 2 Mar 2018
Hey Ashley, referring to Equation 1.1 in your paper, q appeared as a partial differential equation and then as an ODE in Eq 1.5, can you tell why it is so? If it is a PDE, what are the initial and boundary conditions on q (C has these conditions given but not for q, recall that when you solve an ODE without the initial conditions given, you obtain a general solution, with these conditions given, you can obtain a particular solution), these conditions are necessary to obtain reasonable solutions. Do you have an idea of the kind of plot that you are expecting from the solutions of your equations, a plot would be helpful to allow others to work on the problem. I am trying to apply finite difference method to solve the problem, but would need some answers to my questions.
8 Comments
Nelson Guy
on 4 Jul 2019
any clues Ash and Trosten for modeling binary gas mixtures in breakthrough mode? say Propane/Propylene
Abraham Boayue
on 5 Mar 2018
I still have come up with a reasonable solution, but I am making some progress so far. Here is an image that I obtained from using the finite difference method. Take a look at it and tell me what you think.
0 Comments
Abraham Boayue
on 5 Mar 2018
An image-
2 Comments
Abraham Boayue
on 6 Mar 2018
The graphs represent the concentration profile. The x axis is the timeline and the vertical axis shows how the concentration depreciates over time.
Abraham Boayue
on 6 Mar 2018
I am using the following papers as a guide: http://www.cellulosechemtechnol.ro/pdf/CCT7-8(2014)/p.717-726.pdf
0 Comments
Patricia Sáez González
on 10 Sep 2019
function main
% Initialisation
close all
clear all
clc
%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
b = 2.5*10^-5; % Langmuir parameter
qs = 2*10^4; % Saturation capacity
cFeed = 10; % Feed concentration
L = 1; % Column length
t0 = 0; % Initial Time
tf = 1000; % Final time
dt = 0.5; % Time step
z = [0:0.01: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);
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
b = 2.5*10^-5; % Langmuir parameter
qs = 20000; % Saturation capacity
% 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));
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*( ((qs*b*c(1))/(1 + b*c(1))) - q(1) );
DcDt(1) = 0.0;
% Set time derivatives in remaining points
for i=2:n
%Equation: dq/dt = k(q*-q) where q* = qs * (b*c)/(1+b*c)
DqDt(i) = k*( ((qs*b*c(i))/(1 + b*c(i))) - 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, so, I would like to know how to obtain these curves (C/C0 vs t) from this code. I would appreciate any help. Thank you very much!
4 Comments
Luis Gilberto Domínguez López
on 25 May 2020
Hi there!
You have help me a lot with this code thank you so much. I was wondering what would happen if instead of a "step" change (breakthrough analysis) you make a rectangular pulse of concentration. I was trying to implement this change in the code unseccesfully because of my lack of experience. could you give some advice?
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/299983/image.png)
nashwa tarek
on 13 Jul 2020
hi patricia ,
i want to ask you , your code doesnt study the effect of initial concentration, can you have a clue what if i want to study the effect of concentration with the column length and the velocity
nashwa tarek
on 4 Jul 2020
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
0 Comments
Federico Urbinati
on 18 Oct 2022
hi does anyone have the matlab code to solve adsorption when you have multiple components and not just one?
I hope someone can answer, I can not go forward. thank you
0 Comments
See Also
Categories
Find more on Dimensionality Reduction and Feature Extraction in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!