Adsorption Modelling - Solving PDE - Axial Dispersion Model

Question

Ashley Gratton 2018-3-1

0
链接

此问题的直接链接

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

评论： Federico Urbinati 2022-10-19

采纳的回答： Torsten

AdsorptionEquations.pdf

在 MATLAB Online 中打开

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 个评论
显示 44更早的评论隐藏 44更早的评论

Ashley Gratton 2018-6-22

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 2022-5-3

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 2022-10-19

@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

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

Torsten 2018-3-2

3
链接

此回答的直接链接

https://ww2.mathworks.cn/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 个评论
显示 18更早的评论隐藏 18更早的评论

Ashley Gratton 2018-3-2

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 2018-3-2

在 MATLAB Online 中打开

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 2018-3-5

在 MATLAB Online 中打开

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 2019-6-26

在 MATLAB Online 中打开

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 2019-7-9

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.

请先登录，再进行评论。

Answer 2

Abraham Boayue 2018-3-2

1
链接

此回答的直接链接

https://ww2.mathworks.cn/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 个评论
显示 6更早的评论隐藏 6更早的评论

Nelson Guy 2019-4-24

编辑：Nelson Guy 2019-4-24

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 2019-4-24

编辑：Ashley Gratton 2019-4-24

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 2019-4-25

编辑：Nelson Guy 2019-4-28

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 2019-4-27

and Thanks Ash and Trosten.

Nelson

Nelson Guy 2019-7-4

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

请先登录，再进行评论。

Answer 3

Abraham Boayue 2018-3-5

1
链接

此回答的直接链接

https://ww2.mathworks.cn/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 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

Answer 4

Abraham Boayue 2018-3-5

1
链接

此回答的直接链接

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

AxialDisp.PNG

An image-

2 个评论
显示无隐藏无

Ashley Gratton 2018-3-6

Hello Abraham,

What do each of the axes represent?

Thank you,

Ash

Abraham Boayue 2018-3-6

The graphs represent the concentration profile. The x axis is the timeline and the vertical axis shows how the concentration depreciates over time.

请先登录，再进行评论。

Answer 5

Abraham Boayue 2018-3-6

0
链接

此回答的直接链接

https://ww2.mathworks.cn/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 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

Answer 6

Patricia Sáez González 2019-9-10

0
链接

此回答的直接链接

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

在 MATLAB Online 中打开

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 个评论
显示 2更早的评论隐藏 2更早的评论

Torsten 2019-9-10

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 2020-5-25

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 2020-7-13

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

请先登录，再进行评论。

Answer 7

nashwa tarek 2020-7-4

0
链接

此回答的直接链接

https://ww2.mathworks.cn/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