Adsorption Modelling - Solving PDE - Axial Dispersion Model
54 vues (au cours des 30 derniers jours)
Afficher commentaires plus anciens
Ashley Gratton
le 1 Mar 2018
Commenté : Federico Urbinati
le 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 commentaires
Ridhwan Lawal
le 3 Mai 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
le 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
Réponse acceptée
Torsten
le 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 commentaires
Yash Toshniwal
le 26 Juin 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
le 9 Juil 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.
Plus de réponses (7)
Abraham Boayue
le 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 commentaires
Nelson Guy
le 4 Juil 2019
any clues Ash and Trosten for modeling binary gas mixtures in breakthrough mode? say Propane/Propylene
Abraham Boayue
le 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 commentaires
Abraham Boayue
le 5 Mar 2018
An image-
2 commentaires
Abraham Boayue
le 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
le 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 commentaires
Patricia Sáez González
le 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 commentaires
Luis Gilberto Domínguez López
le 25 Mai 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
le 13 Juil 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
le 4 Juil 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 commentaires
Federico Urbinati
le 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 commentaires
Voir également
Catégories
En savoir plus sur Dimensionality Reduction and Feature Extraction dans Help Center et File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!