Adsorption Modelling - Solving PDE - Axial Dispersion Model

Question

Ashley Gratton on 1 Mar 2018

0
Link

Direct link to this question

https://nl.mathworks.com/matlabcentral/answers/385756-adsorption-modelling-solving-pde-axial-dispersion-model

Commented: Federico Urbinati on 19 Oct 2022

AdsorptionEquations.pdf

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
Show 44 older commentsHide 44 older comments

Ashley Gratton on 22 Jun 2018

I think you are misunderstanding. When solving a multivariate equation (PDE) such as the axial dispersion model (ADM), that varies both spatially and with time, you only need to discretise the PDE in the spatial domain. Then, you are effectively left with a system of ODEs, which can be solved using any iterative ODE method you like, be it ode45/ode23, Euler's method, RK4 etc. So the process looks something like this:

Disctretise the PDE in the spatial domain using the centralised difference method (a step in the method of lines). You will have to choose how many mesh points to have along your spatial domain (N) - i.e. how many slices to split your adsorption column up into. Your mesh grid step size (\delta x) is then defined as the column length (L) divided by the number of mesh points (N) ==> \delta x = L/N. Now you have a column that is discretised spatially (i.e. split up into a tangible number of sections - not a continuum).
Now, at each of these mesh points along the column, the ADM can be applied. For example, this means that for 20 mesh points (N = 20), you will have an ODE (Equation 3.23 in my previously attached document) at each point that varies ONLY with time - giving you 20 ODEs (i.e. a system of ODEs) to solve however you would usually solve an ODE.
Set up a loop to solve this system of ODEs between your initial (t0) and final (tf) times. I used a function file that defines the ADM and calls said function file in the main script using ode45.
Plot the outlet concentration (Cout) against time (t) for the breakthrough curve. Plot a given column of the concentration vector (C) against length (z) to obtain a concentration profile at a given time (determined by which column of the C matrix you plot). Note: I created the z vector as z = [0 : \delta x : L].

Hope this is clear, I've tried to explain as best I can, but I understand it may all be a bit confusing - this project was effectively my entire dissertation, so trying to explain in a few paragraphs is tricky! If you want any more info, let me know, I will try to help.

Best,

Ash

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

@ALOK PRASAD

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

Sign in to comment.

Sign in to answer this question.

Answer 1

Torsten on 2 Mar 2018

3
Link

Direct link to this answer

https://nl.mathworks.com/matlabcentral/answers/385756-adsorption-modelling-solving-pde-axial-dispersion-model#answer_308000

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
Show 18 older commentsHide 18 older comments

Ashley Gratton on 2 Mar 2018

Hello,

Thank you for your answer.

I have had a look at your example, but still cannot seem to get this to work. The problem that I am getting stuck on is that I have an ODE within a PDE and need to solve the two simultaneously, but the PDE's first solution in the method of lines matrix depends on the first solution of the ODE, which depends on the first solution of the PDE... It seems like a circular argument...

How would you discretise the PDE when I have two variables (c and q, where q* = fn(c) [see attached file])?

Essentially I am unsure on:

a) How to discretise the PDE using MoL due to the complexity of the system

b) How to solve the system of ODEs generated by the MoL and the dq/dt ODE simultaneously. (Usually, I would know how to solve a system of ODEs, but in this case, I have a system of ODEs from the MoL plus another ODE, which throws me off).

Any further suggestions on how to tackle this problem would be great! If you would like me to provide any further info, let me know, as I understand what I am trying to communicate may not be entirely clear!

Thanks again!

Ashley Gratton on 2 Mar 2018

Thank you again for your answer, I really appreciate your help.

I have modified your provided code for my problem, however, I am not sure where I have made errors. I will provide the code below. Please note: I have commented the code to make my interpretation of your code clear to you, alongside any questions.

% 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.1:L]; % Mesh generation
t = [t0:dt:tf];% Time vector
n = numel(z);  % Size of mesh grid
%Initial Conditions / Vector Creation
c0 = ones(n,1);     % t = 0, c = 0 for all z, so shouldn't this be zeros too?
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);
% Expected results:
% c is the concentration in the gas phase, so should eventually reach the
% feed concentration (cFeed) when the adsorbent in the column is fully saturated.
% So effectively, the column should reach steady-state whereby
% for all z, c = cfeed.
% Note q = concentration of gas adsorbed (i.e. on the adsorbent surface)
% Problems:
% I am unsure where to code in the cFeed condition so that the above result
% is achieved.
% Y matrix just displays loads of ones for 1/2 of the matrix, but then does
% give something for the latter half?

Function:

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
cFeed = 10;     % Feed concentration
% Variables being allocated empty 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,1);
D2cDz2 = zeros(n-1,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;
% Boundary Condition at z = L (dc/dz = 0 at z = L)
%DuDx(1) = 0.0;
DcDz(n) = 0; % Changed from commented-out part above, can I do this?
% Since my BC is dc/dz = 0 at z = L (i.e. at the Nth point)
% Discretisation process
D2cDz2(1) = zhalf(1)/(zhalf(1)-z(1)) * (c(2)-c(1))/(z(2)-z(1));
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
for i=1:n-1
    %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);
    %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) );
end
% Unsure of the purpose of these? %%%%%%%%%%%
DcDt(n) = 0.0;                              %
                                            %
DqDt(n) = k*(((qs*b*c(n))/(1+b*c(n)))-q(n));%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Regardless, I have changed them to fit my equations.
DyDt = [DcDt;DqDt];
end

I understand how you have set up the system (i.e. discretised the PDE) for the most part, I have commented parts that I do not understand if I have modified them correctly for my problem. You have discretised the system slightly differently to how I learned to do it, hence my trouble understanding!

I have tried but am not certain on how to enter my boundary / initial conditions, which are as follows:

t = 0, c = 0 & q = 0, for all z

z = 0, c = cfeed, for t>0

z = L, dc/dz = 0, for t>0

Again, thank you so much for your help, I appreciate it.

Torsten on 5 Mar 2018

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

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.

Sign in to comment.

Answer 2

Abraham Boayue on 2 Mar 2018

1
Link

Direct link to this answer

https://nl.mathworks.com/matlabcentral/answers/385756-adsorption-modelling-solving-pde-axial-dispersion-model#answer_308108

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
Show 6 older commentsHide 6 older comments

Nelson Guy on 24 Apr 2019

Edited: Nelson Guy on 24 Apr 2019

Oscillation.pdf

Hello

I used the same code for mimicking breakthrough curve for single gas adsorption on solid zeolite- fixed bed.

%Assumed plug flow, langmuir, Dubinin-A, LDF

% Initialization

close all

clear all

clc

%System Set-up %

%Propane Zeolite

D = 3.37*10^-13; % Axial Dispersion coefficient [m^2/s]

v =1*10^-3; % Superficial velocity [m/s]

epsilon = 0.4; % Voidage fraction []

k = 0.02516; % Mass Transfer Coefficient [1/s]

cFeed = 10; % Feed concentration [mol/m^3]

d=1285720; % density zeolite [g/m^3]Y-DP

L = 0.1; % Column length [m]

%isotherm parameters

E=6224;

W0=41.72*4.4615;

n=0.8646;

b=10^(-4);

qs=2.5;

%Gas parameters

Temp=273.15+60; %Temperature

ps=869200; %Saturation pressure of gas

R=8.3144621; %Gas constant

Cs=ps/(R*Temp);

qe=@(C) qs*b*C/(1+b*C);

%qe=@(C) W0*exp(-(R*Temp*log(Cs./C)/E).^n);

t0 = 0; % Initial Time

tf = 10000; % Final time

dt = 100; % Time step

z = [0:L/20: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] = ode23s(@(t,y) user(t,y,z,n,D,v,epsilon,k,d,qe),t,y0);

plot(T,Y);

%filename = 'testdata.xlsx';

%xlswrite(filename,[t.' T Y])

function DyDt=user(~, y, z, n,D,v,epsilon,k,d,qe)

% 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*( qe(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*( qe(c(i)) - q(i) );

%Equation: dc/dt = D * d2c/dz2 - v/e*dc/dz - (d*(1-e)/(e))*dq/dt

DcDt(i) = D*D2cDz2(i) - (v/epsilon)*DcDz(i) - (d*(1-epsilon)/(epsilon))*DqDt(i);

end

% Concatenate vector of time derivatives

DyDt = [DcDt;DqDt];

end

MathWorks

But aftter sometimes the curve starts oscillating. AGAIN Same question above Why please? is this a bounday condition, mesh or time step issue etc?

Regards

Nelson

Ashley Gratton on 24 Apr 2019

Edited: Ashley Gratton on 24 Apr 2019

Hi Nelson,

Perhaps you've tried this already but I would suggest incrementally reducing the time step (

) to try and eliminate the oscillations in your concentration-time plot. Just be aware that as you decrease your time step, you will increase the computation time. If that doesn't work, then I'm out of ideas - personally I would try cycling between different ODE solvers (i.e. ode45/23/15 etc.) to see if the oscillations stop.

I also noticed that your density value is given in units of

not in standard SI units (

), so that could possibly be a source of instability if your other parameters are in SI units and you haven't accounted for the conversion. I would also do a sensitivity analysis on your system (i.e. observe the model response to changing different parameters, whilst keeping others constant). For example, I varied my axial dispersion coefficient across several orders of magnitude and ran the simulation with each value, to get an idea of how the system would respond - what I found was that higher values of the axial dispersion coefficient caused instability in the model. Performing a sensitivity analysis may highlight parameters that are particularly prone to causing instability in your model.

Hope this helps...

Ash

Nelson Guy on 25 Apr 2019

Edited: Nelson Guy on 28 Apr 2019

Thanks, it worked this way:

z = [0:L/500:L]; % Mesh generation

plot(T,Y);

DcDz(i) = (c(i)-c(i-1))/(z(i)-z(i-1));

I would like to perform a binary mixture 50% Propane and 50% Propylene breakthroughs for selectivity purposes?

Any ideas for the code below

%Assumed plug flow, langmuir, Dubinin-A, LDF

% Initialization

close all

clear all

clc

%System Set-up %

%Propane/Propylene zeolite

D1 =3.37E-13; % Axial Dispersion coefficient [m^2/s] of propane

D2= 1.61E-12; % Axial Dispersion coefficient [m^2/s] of propylene

v =1*10^-3; % Superficial velocity [m/s]

epsilon = 0.4; % Voidage fraction []

k1 =0.02516; % Mass Transfer Coefficient [1/s] of propane

k2 =0.2074; % Mass Transfer Coefficient [1/s] of propylene

c1Feed=5; % Feed concentration [mol/m^3]of propane

c2Feed=5; % Feed concentration [mol/m^3]of propylene

d=1285720; % density acetylftw [g/m^3]Y-DP

L = 0.1; % Column length [m]

%isotherm parameters

E=6224;

W0=41.72*4.4615;

n=0.8646;

b1=10^(-4);

qs1=53;

b2= 10^(-3);

qs2=56;

%Gas parameters

Temp=273.15+60; %Temperature

ps=869200; %Saturation pressure of gas

R=8.3144621; %Gas constant

Cs=ps/(R*Temp);

qe=@(C1,C2) (qs1*b1*C1/(1+b1*C1))+(qs2*b2*C2/(1+b2*C2));

%qe=@(C1,C2) W01*exp(-(R*Temp*log(Cs1./C1)/E1).^n1)+W02*exp(-(R*Temp*log(Cs2./C2)/E2).^n2);

t0 = 0; % Initial Time

tf = 10000; % Final time

dt = 100; % Time step

z = [0:L/500:L]; % Mesh generation

t = [t0:dt:tf]; % Time vector

n = numel(z); % Size of mesh grid

%Initial Conditions / Vector Creation

c01 = zeros(n,1);

C02=zeros(n,1);

c01(1) = c1Feed;

c02(1)=c2Feed;

q01 = zeros(n,1); % t = 0, q = 0 for all z

q02=zeros(n,1);

y01 = [c01 ; q01]; % Appends conditions together

y02=[c02;q02]

%ODE15S Solver

[T1, Y1] = ode23s(@(t,y1) user(t,y1,z,n,D1,v,epsilon,k1,d1,qe1),t,y01);

[T2, Y2] = ode23s(@(t,y2) user(t,y2,z,n,D2,v,epsilon,k2,d2,qe2),t,y02);

plot(T1,Y1);

plot(T2,Y2);

filename = 'testdata.xlsx';

xlswrite(filename,[t.' T1 Y1 T2 Y2])

function DyDt=user(~, y, z, n,D,v,epsilon,k,d,qe)

% Variables being allocated zero vectors

c1 = zeros(n,1);

c2= zeros(n,1);

q1 = zeros(n,1);

q2= zeros(n,1);

Dc1Dt = zeros(n,1);

Dc2Dt = zeros(n,1);

Dq1Dt = zeros(n,1);

Dq2Dt = zeros(n,1);

Dy1Dt = zeros(2*n,1);

Dy2Dt = zeros(2*n,1);

zhalf = zeros(n-1,1);

Dc1Dz = zeros(n,1);

Dc2Dz = zeros(n,1);

D2c1Dz2 = zeros(n,1);

D2c2Dz2 = zeros(n,1);

c1 = y1(1:n);

c2 = y2(1:n);

q1= y1(n+1:2*n);

q2= y2(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

Dq1Dt(1) = k1*( qe1(c1(1)) - q1(1) );

Dq2Dt(1) = k2*( qe2(c2(1)) - q2(1) );

Dc1Dt(1) = 0.0;

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

Dq1Dt(i) = k1*( qe1(c1(i)) - q1(i) );

Dq2Dt(i) = k2*( qe2(c2(i)) - q2(i) );

%Equation: dc/dt = D * d2c/dz2 - v/e*dc/dz - (d*(1-e)/(e))*dq/dt

Dc1Dt(i) = D1*D2c1Dz2(i) - (v/epsilon)*Dc1Dz(i) - (d1*(1-epsilon)/(epsilon))*Dq1Dt(i);

Dc2Dt(i) = D2*D2c2Dz2(i) - (v/epsilon)*Dc2Dz(i) - (d2*(1-epsilon)/(epsilon))*Dq2Dt(i);

end

% Concatenate vector of time derivatives

Dy1Dt = [Dc1Dt;Dq1Dt];

Dy2Dt = [Dc2Dt;Dq2Dt];

end

Regards

N

Nelson Guy on 27 Apr 2019

and Thanks Ash and Trosten.

Nelson

Nelson Guy on 4 Jul 2019

any clues Ash and Trosten for modeling binary gas mixtures in breakthrough mode? say Propane/Propylene

Sign in to comment.

Answer 3

Abraham Boayue on 5 Mar 2018

1
Link

Direct link to this answer

https://nl.mathworks.com/matlabcentral/answers/385756-adsorption-modelling-solving-pde-axial-dispersion-model#answer_308520

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
Show -2 older commentsHide -2 older comments

Sign in to comment.

Answer 4

Abraham Boayue on 5 Mar 2018

1
Link

Direct link to this answer

https://nl.mathworks.com/matlabcentral/answers/385756-adsorption-modelling-solving-pde-axial-dispersion-model#answer_308521

AxialDisp.PNG

An image-

2 Comments
Show NoneHide None

Ashley Gratton on 6 Mar 2018

Hello Abraham,

What do each of the axes represent?

Thank you,

Ash

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.

Sign in to comment.

Answer 5

Abraham Boayue on 6 Mar 2018

0
Link

Direct link to this answer

https://nl.mathworks.com/matlabcentral/answers/385756-adsorption-modelling-solving-pde-axial-dispersion-model#answer_308524

I am using the following papers as a guide: http://www.cellulosechemtechnol.ro/pdf/CCT7-8(2014)/p.717-726.pdf

http://www.imrfjournals.in/pdf/MATHS/MSIRJ-Volume-3-Issue-1-2014/28.pdf

0 Comments
Show -2 older commentsHide -2 older comments

Sign in to comment.

Answer 6

Patricia Sáez González on 10 Sep 2019

0
Link

Direct link to this answer

https://nl.mathworks.com/matlabcentral/answers/385756-adsorption-modelling-solving-pde-axial-dispersion-model#answer_391182

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
Show 2 older commentsHide 2 older comments

Torsten on 10 Sep 2019

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
%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
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
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));
  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*( ((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

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?

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

Sign in to comment.

Answer 7

nashwa tarek on 4 Jul 2020

0
Link

Direct link to this answer

https://nl.mathworks.com/matlabcentral/answers/385756-adsorption-modelling-solving-pde-axial-dispersion-model#answer_460805

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