How to solve coupled partial differential equations with method of lines?

I am getting problem in solving the partial differential equations used for modelling of packed bed adsorption column. I have attached the equations which I need to solve.

I solved it using method of lines approach with the help of some code which I got on mathworks ask community. But the graph which I am getting is different from which I need to get.

Can anyone help me in solving these equations?

I am also providing the code which I used to solve.

% parameters
eps = 0.41;             % bed porosity
volflow = 5e-7;         % volumetric flow rate in m3/sec
dia = 0.022;            % bed internal dia in m
bedarea = pi*dia^2/4;
u = volflow/bedarea;    % superficial velocity in m/s
Dl = 0.08e-4;    % axial dipsersion coefficient in m2/sec
rhop = 1228.5;  % particle density 
Kl = 0.226;     % overall mass transfer coeff in sec-1
P = 102000;     % total pressure in pascals
R = 8.314;      % gas constant in jol/mol-k
T = 500;        % gas temp in K
yco2 = 0.2;     % co2 mole fraction
cfeed = yco2*P/(R*T);       % feed conc in mol/m3
K0 = 4.31e-9;           % in pascal-1
delh = -29380;          % heat of adsorption in j/mol
Keq = K0*exp(-delh/(R*T));      % equilibrium constant in pascal-1

% using method of lines
tf = 300;      % End time
nz = 400;       % Number of gateways
np = nz + 1;    % Number of gateposts

% Initial Concentrations
C = zeros(np,1);  
q = zeros(np,1);

c0 = [C; q];   
dt = tf/300;
tspan = 0:dt:tf;

[t, c] = ode15s(@rates, tspan, c0);

% Plot results
plot(t,c(:,np-4)/cfeed),grid
xlabel('time'), ylabel('C/Cfeed')

function dcdt = rates(~,c)
eps = 0.41;             % bed porosity
volflow = 5e-7;         % volumetric flow rate in m3/sec
dia = 0.022;            % bed internal dia in m
bedarea = pi*dia^2/4;
u = volflow/bedarea;    % superficial velocity in m/s
Dl = 0.08e-4;    % axial dipsersion coefficient in m2/sec
rhop = 1228.5;  % particle density 
Kl = 0.226;     % overall mass transfer coeff in sec-1
P = 102000;     % total pressure in pascals
R = 8.314;      % gas constant in jol/mol-k
T = 500;        % gas temp in K
yco2 = 0.2;     % co2 mole fraction
cfeed = yco2*P/(R*T);       % feed conc in mol/m3

K0 = 4.31e-9;           % in pascal-1
delh = -29380;          % heat of adsorption in j/mol
Keq = K0*exp(-delh/(R*T));      % equilibrium constant in pascal-1
qm = 5.09;          % max adsrobed conc in mol/kg
n = 0.429;          % toth isotherm parameter
L = 0.171;
nz = 400;       
np = nz + 1;   
dz = L/nz;

C = c(1:np);
q = c(np+1:2*np);
dCdt = zeros(np,1);
dqdt = zeros(np,1);

%  boundary condition at z=0
 C(1) = (1/(eps*Dl+u*dz))*((eps*Dl*C(2)+u*cfeed*dz));

for i = 2:np-1
       dCdt(i) = (Dl/dz^2)*(C(i+1)-2*C(i)+C(i-1)) - (u/(2*eps*dz))*(C(i+1)-C(i-1)) - (((1-eps)*rhop)/eps)*dqdt(i);
       dqdt(i) = Kl*((qm*Keq*C(i)*R*T/(1+(Keq*C(i)*R*T)^n)^(1/n)) - q(i)); 
end

% boundary conditions at z=L
dCdt(np) = dCdt(np-1); 

% Combine rates for the function to return to the calling routine
dcdt = [real(dCdt); real(dqdt)];

end

Here I used time interms of seconds. So when I plot the graph the x-axis is also in seconds and my breakthrough time which I am getting is around 35 seconds. But in the experimental values for showed that the breakthrough time for the graph is around 300 minutes.

so please someone help whether there is any mistake which I am doing in my code.

[Davide Masiello]

I can see several problems with your code. For instance:

1) you are not actually imposing a boundary condition at z = 0, because what the functions really sees is that dCdt is zero at that point (from the pre-allocation of the dCdt vector). If you need to impose an algebraic boundary condition then your problem becomes a DAE system and you need to specify that to MatLab.

2) Why do you write

dcdt = [real(dCdt); real(dqdt)];

?

There should not be the need to use the real function.

Some help could be devised if you could post the governing equations and b.c.

[Ari Dillep Sai]

Hi Davide,

Thanks for your reply.

I am attaching the governing equations.

[Davide Masiello]

What are z+ and z- in the flux boundary condition?

[Ari Dillep Sai]

C at z+ corresponds to the concentration after entering into the column and C at z- corresponds to the concentration before entering into the column which is inlet concentration

Answer 1

[Davide Masiello]

Hi @Ari Dillep Sai, take a look at the code below.

The breakthrough now is found to be sometimes after 4000 s (or ~67 mins), so not the 300 mins you have suggested, but maybe a bit better than the 35 s you were getting before.

Also note the following improvements to the code:

1 - I have passed all the constants as external parameters, so you don't have to define them again in the function.

2 - I have written the method of lines without use of for loops, but with wise use of logical indexing.

3 - ode15s solves a system of DAEs now, which is necessary given the algebraic nature of the boundary conditions for the diffusion equation.

I am not sure this will work for you but it might be a good starting point to develop your final solution.

I suggest you go throught the equations and the values of the parameters again, and ensure that everything is correct.

clear,clc

% Parameters
p.eps         = 0.41;                           % bed porosity
p.volflow     = 5e-7;                           % volumetric flow rate in m3/sec
p.dia         = 0.022;                          % bed internal dia in m
p.bedarea     = pi*p.dia^2/4;
p.u           = p.volflow/p.bedarea;            % superficial velocity in m/s
p.Dl          = 0.08e-4;                        % axial dipsersion coefficient in m2/sec
p.rhop        = 1228.5;                         % particle density 
p.Kl          = 0.226;                          % overall mass transfer coeff in sec-1
p.P           = 102000;                         % total pressure in pascals
p.R           = 8.314;                          % gas constant in jol/mol-k
p.T           = 500;                            % gas temp in K
p.yco2        = 0.2;                            % co2 mole fraction
p.cfeed       = p.yco2*p.P/(p.R*p.T);           % feed conc in mol/m3
p.K0          = 4.31e-9;                        % in pascal-1
p.delh        = -29380;                         % heat of adsorption in j/mol
p.Keq         = p.K0*exp(-p.delh/(p.R*p.T));    % equilibrium constant in pascal-1
p.qm          = 5.09;                           % max adsrobed conc in mol/kg
p.n           = 0.429;                          % toth isotherm parameter
p.L           = 0.171;                          % bed length

% using method of lines
tf      = 300*60;                               % End time
p.nz    = 400;                                  % Number of nodes
z       = linspace(0,p.L,p.nz);                 % Space domain
dz      = diff(z); p.dz = dz(1);                % Space increment


% Initial conditions
c       = zeros(p.nz,1);
q       = zeros(p.nz,1);
y0      = [c;q];   
tspan   = [0 tf];

% Mass matrix (for DAE system)
M               = eye(2*p.nz);
M(1,1)          = 0;
M(p.nz,p.nz)    = 0;
options         = odeset('Mass',M,'MassSingular','yes');
[t,y]           = ode15s(@(t,y)rates(t,y,p),tspan,y0,options);

% Plot results
plot(t/60,y(:,p.nz)/p.cfeed)
grid
xlabel('time (mins)')
ylabel('C/Cfeed')
title('Dimensionless concentration at ''z = L''')

hold on
plot(t/60,0.05*ones(size(t)),'--r')

legend('C/Cfeed','breakthrough')

function out = rates(~,y,p)
% Variables
c       = y(1:p.nz,1);
q       = y(p.nz+1:2*p.nz,1);

% Space derivatives
dcdz            = zeros(p.nz,1);
d2cdz2          = zeros(p.nz,1);
dcdz(2:end-1)   = (c(3:end)-c(1:end-2))/(2*p.dz);
d2cdz2(2:end-1) = (c(3:end)-2*c(2:end-1)+c(1:end-2))/(p.dz^2);

% Saturation
qstar = p.qm*p.Keq*c*p.R*p.T./((1+(p.Keq*c*p.R*p.T).^(p.n)).^(1/p.n));

% Governing equations
dqdt = p.Kl*(qstar-q); 
dcdt = p.Dl*d2cdz2-p.u*dcdz/p.eps-p.rhop*(1-p.eps)*dqdt/p.eps;
       
% Boundary conditions
dcdt(1)     = p.eps*p.Dl*(c(2)-c(1))/p.dz+p.u*(p.cfeed-c(1));
dcdt(end)   = c(end)-c(end-1);

% Output
out = [dcdt;dqdt];
end

[Ari Dillep Sai]

Thank you Davide!

Can you please explain me how you written mass matrix or can you provide any documentation if you have?

[Davide Masiello]

Hi @Ari Dillep Sai, I would suggest a careful reading of the documentation at the following link

https://ch.mathworks.com/help/matlab/math/solve-differential-algebraic-equations-daes.html

If you think my answer could be helpful for other users with the same problem, consider clicking the Accept button on the top left.

Cheers

[Ari Dillep Sai]

Thank you so much for the support Davide!

[Una]

Dear Davide,

I have the same mass transfer equation as Ari but with different isotherm. I have tried to follow your code but i get an error massage. I dont know what i am doing wrong. I hope you can help me.

clear all
clc

% Design parameters of the adsorbent bed
p.eps         = 0.87;                                    % bed porosity
%p.volflow     = 2.86;                                    % volumetric flow rate in m3/sec
%p.dia         = 5.05;                                    % bed internal dia in m
%p.bedarea     = pi*p.dia^2/4;                            % bed area in m2
p.u           = 3;                                       % superficial velocity in m/s
p.rhop        = 500;                                     % particle density in kg/m3
p.r_1         = 0.000525;                                % particle radius in m
p.r_2         = 0.000585;                                % particle radius in m                          
p.L           = 0.3;                                     % bed length in m

% General parameters
p.Sc          = 2.03e-2;                                 % schmidt number                   
p.De          = 1e-11;                                   % effective diffusivity 
p.Dm          = 1.6e-6;                                  % molecular diffusivity of co2 in air in m2/s
p.Dl          = p.Dm + (p.u*p.r_1*2)^2/(192*p.Dm);       % axial dipsersion coefficient in m2/sec
%p.rohg        = 1.2;                                     % density of air in kg/m3
%p.viskog      = 18.37e-6;                                % dynamic viscosity of air in Ns/m2
%p.Re          = 190;
%p.sh          = 2+1.1*p.Re^0.6+p.Sc^0.33;
%p.kf          = (p.sh*p.Dm)/p.r_1;
p.Kl          = 0.005;
%p.P           = 101325;                                   % total pressure in pascals
p.R           = 8.314;                                    % gas constant in J/mol-k
p.T           = 298;                                      % gas temp in K
%p.yco2        = 0.0004;                                   % co2 mole fraction
p.cfeed       = 16;                     % feed conc in mol/m3

% Parameters of the isotherm (depend on the adsorbent used)
p.qsat1       = 0.16;
p.b1          = 4.89;
p.qsat2       = 3.50;
p.b2          = 22;
p.b3          = 0.00062;

% using method of lines
tf      = 300*60;                                        % End time
p.nz    = 300;                                           % Number of nodes
z       = linspace(0,p.L,p.nz);                          % Space domain
dz      = diff(z); 
p.dz    = dz(1);                                         % Space increment

% Initial conditions
c       = zeros(p.nz,1);
q       = zeros(p.nz,1);
y0      = [c;q];   
tspan   = [0 tf];

% Mass matrix (for DAE system)
M               = eye(2*p.nz);
M(1,1)          = 0;
M(p.nz,p.nz)    = 0;
options         = odeset('Mass',M,'MassSingular','yes');
[t,y]           = ode15s(@(t,y)rates(t,y,p),tspan,y0,options);

ans = 1×2

300 300

Error using odearguments
@(T,Y)RATES(T,Y,P) must return a column vector.

Error in ode15s (line 153)
odearguments(odeIsFuncHandle, odeTreatAsMFile, solver_name, ode, tspan, y0, options, varargin);

% Plot results
plot(t,y(:,p.nz)/p.cfeed)
grid
xlabel('time (mins)')
ylabel('C/Cfeed')
title('Dimensionless concentration at ''z = L''')
hold on

function out = rates(~,y,p)
% Variables
c       = y(1:p.nz,1);
q       = y(p.nz+1:2*p.nz,1);

% Space derivatives
dcdz            = zeros(p.nz,1);
d2cdz2          = zeros(p.nz,1);
dcdz(2:end-1)   = (c(3:end)-c(1:end-2))/(2*p.dz);
d2cdz2(2:end-1) = (c(3:end)-2*c(2:end-1)+c(1:end-2))/(p.dz^2);

% Isotherm

if p.R*p.T*c < 0.16
    qstar=(p.qsat1*p.b1*p.R*p.T*c)/(1+p.b1*p.R*p.T*c);
else 
    qstar=(p.qsat2*p.b2*(p.R*p.T*c-0.16))./(1+p.b2*(p.R*p.T*c-0.16)) + p.b3*p.R*p.T*c;
end
size(qstar)

% Governing equations
dqdt = p.Kl*(qstar-q); 
dcdt = p.Dl*d2cdz2-p.u*dcdz/p.eps-p.rhop*(1-p.eps)*dqdt/p.eps;
       
% Boundary conditions
dcdt(1)     = p.eps*p.Dl*(c(2)-c(1))/p.dz+p.u*(p.cfeed-c(1));
dcdt(end)   = c(end)-c(end-1);

% Output
out = [dcdt;dqdt];
end

[Torsten]

As you can see above, qstar is a 300x300 matrix instead of a 300x1 vector.

[Una]

I am a noob. How can I change that? or is it possible to change it?

[Torsten]

If you are a noob, use a safe loop:

qstar = zeros(p.nz,1);
for i=1:p.nz
    if p.R*p.T*c(i) < 0.16
        qstar(i)=(p.qsat1*p.b1*p.R*p.T*c(i))/(1+p.b1*p.R*p.T*c(i));
    else 
        qstar(i)=(p.qsat2*p.b2*(p.R*p.T*c(i)-0.16))/(1+p.b2*(p.R*p.T*c(i)-0.16)) + p.b3*p.R*p.T*c(i);
    end
end

Is there a continuous transition between the equilibrium loadings (i.e. at R*T*c = 0.16) ?

Otherwise, you will most probably have trouble with the integrator.

[Una]

Thank you veery much Torsten!

I am not sure whether there is a continuous transation. I have taken the isotherm from a peer reviewed paper that works with the same adsorbent as I have to . But now the results seems to be reasonable. Since there is no experimental data available for that isotherm I am trying to validate the code with another isotherm. But now i am am getting another error.

clear all
clc

% Design parameters of the adsorbent bed
p.eps         = 0.4;                                     % bed porosity
%p.volflow     = 2.86;                                    % volumetric flow rate in m3/sec
%p.dia         = 5.05;                                    % bed internal dia in m
%p.bedarea     = 20; %pi*p.dia^2/4;                       % bed area in m2
p.u           = 7.06e-2;                                  % superficial velocity in m/s
p.rhop        = 880;                                     % particle density in kg/m3
p.dp          = 0.00026*2;                              % particle diameter in m
p.Rp          = 0.00026;                                  % particle radius in m
p.L           = 0.01;                                      % bed length in m

% General parameters
%p.Sc          = 2.03e-2;                                 % schmidt number                   
%p.De          = 1e-11;                                   % effective diffusivity 
p.Dm          = 1.3e-5;                                   % molecular diffusivity of co2 in air in m2/s
p.Dl          =(0.45+0.55*p.eps)*1.3e-5+0.35*p.Rp*p.u;    % axial dipsersion coefficient in m2/sec
%p.rohg        = 1.2;                                     % density of air in kg/m3
%p.viskog      = 18.37e-6;                                % dynamic viscosity of air in Ns/m2
%p.Re          = p.rohg*p.u*p.dp/p.viskog;
%p.sh          = 2+1.1*p.Re^0.6+p.Sc^0.33;
%p.kf          = (p.sh*p.Dm)/p.dp;
p.Kl          = 0.003; %1/(((p.Rp/3*p.kf)+(p.Rp^2/15*p.De))*(p.eps/1-p.eps));
%p.P           = 101325;                                  % total pressure in pascals
p.R           = 8.314;                                    % gas constant in J/mol-k
p.T           = 298.15;                                      % gas temp in K
%p.yco2        = 0.0004;                                  % co2 mole fraction
p.cfeed       = 0.016;                                    % feed conc in mol/m3

% Parameters of the isotherm (depend on the adsorbent used)
p.b0          = 93;                                   % Toth parameter in 1/bar
p.delh        = 95.3e3;                               % heat of adsorption in J/mol
p.T0          = 353.15;                               % refernece temperature in K
p.T           = 298.15;                                  % temperature in K
p.b           = p.b0*exp((p.delh/(p.R*p.T0))*((p.T0/p.T)-1));
p.ns0         = 3.40;                                 % Toth perameter in mol/kg
p.exe         = 0;                                    % dimensionless constant 
p.ns          = p.ns0*exp(p.exe*(1-p.T/p.T0));        % max adsrobed conc in mol/kg
p.t0          = 0.37;                                 % Toth parameter
p.alp         = 0.33;                                 % dimensionless constant
p.t           = p.t0+p.alp*(1-p.T0/p.T);


% using method of lines
tf      = 3000;                                        % End time
p.nz    = 400;                                        % Number of nodes
z       = linspace(0,p.L,p.nz);                       % Space domain
dz      = diff(z); 
p.dz = dz(1);                                         % Space increment

% Initial conditions
c       = zeros(p.nz,1);
q       = zeros(p.nz,1);
y0      = [c;q];   
tspan   = [0 tf];

% Mass matrix (for DAE system)
M               = eye(2*p.nz);
M(1,1)          = 0;
M(p.nz,p.nz)    = 0;
options         = odeset('Mass',M,'MassSingular','yes');
[t,y]           = ode15s(@(t,y)rates(t,y,p),tspan,y0,options);
% Plot results
plot(t,y(:,p.nz)/p.cfeed)
grid
xlabel('time (s)')
ylabel('C/Cfeed')
title('concentration at ''z = L''')
hold on

function out = rates(~,y,p)
% Variables
c       = y(1:p.nz,1);
q       = y(p.nz+1:2*p.nz,1);

% Space derivatives
dcdz            = zeros(p.nz,1);
d2cdz2          = zeros(p.nz,1);
dcdz(2:end-1)   = (c(3:end)-c(1:end-2))/(2*p.dz);
d2cdz2(2:end-1) = (c(3:end)-2*c(2:end-1)+c(1:end-2))/(p.dz^2);

% Isotherm

qstar = p.ns*p.b*c*p.R*p.T./((1+(p.b*c*p.R*p.T).^(p.t)).^(1/p.t));

% Governing equations
dqdt = p.Kl*(qstar-q); 
dcdt = p.Dl*d2cdz2-p.u*dcdz-p.rhop*(1-p.eps)*dqdt/p.eps;
       
% Boundary conditions
dcdt(1)     = p.Dl*(c(2)-c(1))/p.dz+p.u*(p.cfeed-c(1));
dcdt(end)   = c(end)-c(end-1);

% Output
out = [dcdt;dqdt];
end

[Torsten]

p.b is Infinity. So check the expression

p.b = p.b0*exp((p.delh/p.R*p.T0)*(p.T0/p.T-1));

for errors.

I am not sure whether there is a continuous transation.

Then I'd plot qstar as a function of c.

[Torsten]

Don't you have inconsistent units since c is in mmol/m^3 and R is in J/(mol*K) ?

And your equation for qstar is quite challenging since the expression b*c(i)*R*T must always remain positive. This is rarely the case in numerical computations if you start with c = 0.

My advice: Plot qstar against c and see if you get reasonable results !!

[Una]

It is supposd to be mol/m3. I have written the unit here incorrectly. I will try to do it. Thank you Torsten.

[Sona]

i am facing a similar problem right now. how can qstar be plotted against c? 😬 I seem to be doing something wrong.

[Torsten]

@Sona Ghulami

I don't know your equation for qstar.

Please include the appropriate code to calculate it.

[Sona]

Thank you Torsten for your fast reply.

function main

cFeed= 0.017;            % Feed concentration in mol/m3
L = 3.2;              % Column length in m
t0 = 0;                % Initial Time in s            
tf= 5000;             % Final time in s
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 and boundary condition %%%%%%%$%
c0 = zeros(n,1);        % c = 0 at t =0 for all z
c0(1)= cFeed;             % c = cFeed at z =0 and t>0
q0 = zeros(n,1);        % t = 0, q = 0 for all z
y0 = [c0 ; q0];         

%%%%%ODE15S Solver%%%%%%
[T, Y] = ode15s(@(t,y) MyFun(t,y,z,n),t,y0);

%%%%%% Plots %%%%%
plot(T,Y(:,n)/cFeed)
legend('c/cFeed')
figure
plot(T,Y(:,2*n))
legend('q')

end

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

epsilon  = 0.5;                         % bed porosity
volflow = 2.86;                        % volumetric flow rate in m3/s
bedarea = 20;                         % bed area in m2
v= volflow/bedarea;                         % superficial velocity in m/s
rho= 55.4;                                    % particle density in kg/m3
Rp = 2.25e-3;                                 % particle radius in m
Dm = 1.3e-5;                                  % molecular diffusivity in m2/s
Dl = (0.45+0.55*epsilon)*Dm+0.35*Rp*v;        % axial dipsersion coefficient in m2/sec
k  = 3.8e-3;
R  = 8.314;                                   % gas constant in J/mol-k
T  = 298;                                     % gas temp in K


% isotherm 
b0  = 2.25e4;                                   % Toth parameter in 1/bar
delh = 6e4;                                      % heat of adsorption in J/mol
T0 = 296;                                      % refernece temperature in K
T  = 298;                                      % temperature in K
b  = b0*exp((delh/(R*T0))*(T0/T-1));
ns0 = 1.97;                                     % Toth perameter in mol/kg
exe = 2.37;                                     % dimensionless constant 
ns  = ns0*exp(exe*(1-T/T0));                    % max adsrobed conc in mol/kg
t0 = 4.22e-1;                                  % Toth parameter
alp= 9.49e-1;                                  % dimensionless constant
t = t0+alp*(1-T0/T);


% Variables being allocated zero vectors
c = zeros(n,1);
q = zeros(n,1);
DcDt = zeros(n,1);
DqDt = zeros(n,1);
zhalf = zeros(n-1,1);
DcDz = zeros(n,1);
DyDt = zeros(2*n,1);
D2cDz2 = zeros(n,1);
c = y(1:n);
q = y(n+1:2*n);

% Interior mesh points
zhalf(1:n-1)=(z(1:n-1)+z(2:n))/2;
for i=2:n-1
 DcDz(i) = (c(i)-c(i-1))/(z(i)-z(i-1));
  D2cDz2(i) = ((c(i+1)-c(i))/(z(i+1)-z(i))-(c(i)-c(i-1))/(z(i)-z(i-1)))/(zhalf(i)-zhalf(i-1));
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
qstar = ns*b*c(1)*R*T./((1+(b*c(1)*R*T).^(t)).^(1/t));             %%%% I have substituted p with R*T*c
DqDt(1) = k*(qstar - q(1) );
DcDt(1) = 0.0;

% Set time derivatives in remaining points
for i=2:n
    % q* 
 qstar = ns*b*c(i)*R*T./((1+((b*c(i)*R*T).^(t))).^(1/t));
    %Equation: dq/dt = k(q*-q)
    DqDt(i) = k*(qstar - q(i) );
      %Equation: dc/dt = D * d2c/dz2 - v*dc/dz - ((1-e)/(e))*roh*dq/dt
      DcDt(i) = Dl*D2cDz2(i) - v*DcDz(i) - ((1-epsilon)/(epsilon))*rho*DqDt(i);
  end
% Concatenate vector of time derivatives
DyDt = [DcDt;DqDt];
end

[Torsten]

Strange curve. Seems it only knows two values: 0 and 1.9...

Maybe a mistake of units in the setting of b0 ? Shouldn't it be in unit 1/Pa ?

epsilon  = 0.5;                         % bed porosity
volflow = 2.86;                        % volumetric flow rate in m3/s
bedarea = 20;                         % bed area in m2
v= volflow/bedarea;                         % superficial velocity in m/s
rho= 55.4;                                    % particle density in kg/m3
Rp = 2.25e-3;                                 % particle radius in m
Dm = 1.3e-5;                                  % molecular diffusivity in m2/s
Dl = (0.45+0.55*epsilon)*Dm+0.35*Rp*v;        % axial dipsersion coefficient in m2/sec
k  = 3.8e-3;
R  = 8.314;                                   % gas constant in J/mol-k
T  = 298;                                     % gas temp in K
% isotherm 
b0  = 2.25e4;                                   % Toth parameter in 1/bar
delh = 6e4;                                      % heat of adsorption in J/mol
T0 = 296;                                      % refernece temperature in K
T  = 298;                                      % temperature in K
b  = b0*exp((delh/(R*T0))*(T0/T-1));
ns0 = 1.97;                                     % Toth perameter in mol/kg
exe = 2.37;                                     % dimensionless constant 
ns  = ns0*exp(exe*(1-T/T0));                    % max adsrobed conc in mol/kg
t0 = 4.22e-1;                                  % Toth parameter
alp= 9.49e-1;                                  % dimensionless constant
t = t0+alp*(1-T0/T);
c = 0:0.01:1;
qstar = ns*b*c*R*T./(1+(b*c*R*T).^t).^(1/t);
plot(c,qstar)

[Sona]

no, i just checked again. It has to be in 1/bar. Do you think it is a problem related to the code and that i should take an other appraoch of solving it? Or maybe my values might have a problem somewhere (Although i have checked them a 100 times now)?

[Torsten]

Sorry, but the above cannot be the correct curve for qstar given c.

There must either be a mistake in the units and/or in the parameters.

Look: The unit of c*R*T is Pa, so in order for b*c*R*T to be dimensionless as required, b has to be in unit 1/Pa, not 1/bar.

[Sona]

that makes sense! I am an idiot 🤦🏻‍♀️ I just used that values that i was given without considering that they might need to be converted. Thank you!

[Torsten]

You should replace your function by a vectorized version.

Note that I changed

Dcdz(n) = 0

to

Dcdz(n) = (c(n)-c(n-1))/(z(n)-z(n-1));

and

b0  = 2.25e4;                                 % Toth parameter in 1/bar

to

b0  = 2.25e4*1e-5;                            % Toth parameter in 1/Pa

function DyDt=MyFun(t, y, z, n) 
epsilon  = 0.5;                               % bed porosity
volflow = 2.86;                               % volumetric flow rate in m3/s
bedarea = 20;                                 % bed area in m2
v= volflow/bedarea;                           % superficial velocity in m/s
rho= 55.4;                                    % particle density in kg/m3
Rp = 2.25e-3;                                 % particle radius in m
Dm = 1.3e-5;                                  % molecular diffusivity in m2/s
Dl = (0.45+0.55*epsilon)*Dm+0.35*Rp*v;        % axial dipsersion coefficient in m2/sec
k  = 3.8e-3;
R  = 8.314;                                   % gas constant in J/mol-k
T  = 298;                                     % gas temp in K
% isotherm
b0  = 2.25e4*1e-5;                            % Toth parameter in 1/Pa
delh = 6e4;                                   % heat of adsorption in J/mol
T0 = 296;                                     % refernece temperature in K
T  = 298;                                     % temperature in K
b  = b0*exp((delh/(R*T0))*(T0/T-1));
ns0 = 1.97;                                   % Toth perameter in mol/kg
exe = 2.37;                                   % dimensionless constant
ns  = ns0*exp(exe*(1-T/T0));                  % max adsrobed conc in mol/kg
t0 = 4.22e-1;                                 % Toth parameter
alp= 9.49e-1;                                 % dimensionless constant
t = t0+alp*(1-T0/T);
% Variables being allocated zero vectors
c = zeros(n,1);
q = zeros(n,1);
DcDt = zeros(n,1);
DqDt = zeros(n,1);
zhalf = zeros(n-1,1);
DcDz = zeros(n,1);
DyDt = zeros(2*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;
% Calculate spatial derivatives
Dcdz(2:n) = (c(2:n)-c(1:n-1))./(z(2:n)-z(1:n-1));
D2cDz2(2:n-1) = ((c(3:n)-c(2:n-1))./(z(3:n)-z(2:n-1))-(c(2:n-1)-c(1:n-2))./(z(2:n-1)-z(1:n-2)))./(zhalf(2:n-1)-zhalf(1:n-2));
% Calculate D2cDz2 at z=L for boundary condition dc/dz = 0
D2cDz2(n) = -1.0/(z(n)-zhalf(n-1))*Dcdz(n);
% Set time derivatives 
qstar = ns*b*c(1:n)*R*T./(1+(b*c(1:n)*R*T).^(t)).^(1/t);             %%%% I have substituted p with R*T*c
DqDt = k*(qstar - q );
DcDt(1) = 0.0;
DcDt(2:n) = Dl*D2cDz2(2:n) - v*DcDz(2:n) - ((1-epsilon)/(epsilon))*rho*DqDt(2:n);
% Concatenate vector of time derivatives
DyDt = [DcDt;DqDt];
end

[Sona]

I get an error:

Unable to perform assignment because the indices on the left side are not compatible with the size of the right side.

Error in main>MyFun (line 71)

Dcdz(2:n) = (c(2:n)-c(1:n-1))./(z(2:n)-z(1:n-1));

function main
cFeed= 0.017;         % Feed concentration in mol/m3
L = 3.2;              % Column length in m
t0 = 0;               % Initial Time in s            
tf= 5000;             % Final time in s
dt= 0.5;              % Time step
z = [0:0.001:L].';    % Mesh generation
t = [t0:dt:tf];       % Time vector
n = numel(z);         % Size of mesh grid


%%%%%%Initial Conditions and boundary condition %%%%%%%$%
c0 = zeros(n,1);        % c = 0 at t =0 for all z
c0(1)= cFeed;           % c = cFeed at z =0 and t>0
q0 = zeros(n,1);        % t = 0, q = 0 for all z
y0 = [c0 ; q0];         

%%%%%ODE15S Solver%%%%%%
[T, Y] = ode15s(@(t,y) MyFun(t,y,z,n),t,y0);

%%%%%% Plots %%%%%
plot(T,Y(:,n)/cFeed)
legend('c/cFeed')
figure
plot(T,Y(:,2*n))
legend('q')
end

function DyDt=MyFun(t, y, z, n) 
epsilon  = 0.5;                               % bed porosity
volflow = 2.86;                               % volumetric flow rate in m3/s
bedarea = 20;                                 % bed area in m2
v= volflow/bedarea;                           % superficial velocity in m/s
rho= 55.4;                                    % particle density in kg/m3
Rp = 2.25e-3;                                 % particle radius in m
Dm = 1.3e-5;                                  % molecular diffusivity in m2/s
Dl = (0.45+0.55*epsilon)*Dm+0.35*Rp*v;        % axial dipsersion coefficient in m2/sec
k  = 3.8e-3;
R  = 8.314;                                   % gas constant in J/mol-k
T  = 298;                                     % gas temp in K

% isotherm
b0  = 2.25e4*1e-5;                            % Toth parameter in 1/Pa
delh = 6e4;                                   % heat of adsorption in J/mol
T0 = 296;                                     % refernece temperature in K
T  = 298;                                     % temperature in K
b  = b0*exp((delh/(R*T0))*(T0/T-1));
ns0 = 1.97;                                   % Toth perameter in mol/kg
exe = 2.37;                                   % dimensionless constant
ns  = ns0*exp(exe*(1-T/T0));                  % max adsrobed conc in mol/kg
t0 = 4.22e-1;                                 % Toth parameter
alp= 9.49e-1;                                 % dimensionless constant
t = t0+alp*(1-T0/T);

% Variables being allocated zero vectors
c = zeros(n,1);
q = zeros(n,1);
DcDt = zeros(n,1);
DqDt = zeros(n,1);
zhalf = zeros(n-1,1);
DcDz = zeros(n,1);
DyDt = zeros(2*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;

% Calculate spatial derivatives
DcDz(2:n) = (c(2:n)-c(1:n-1))./(z(2:n)-z(1:n-1));
D2cDz2(2:n-1) = ((c(3:n)-c(2:n-1))./(z(3:n)-z(2:n-1))-(c(2:n-1)-c(1:n-2))./(z(2:n-1)-z(1:n-2)))./(zhalf(2:n-1)-zhalf(1:n-2));

% Calculate D2cDz2 at z=L for boundary condition dc/dz = 0
D2cDz2(n) = -1.0/(z(n)-zhalf(n-1))*DcDz(n);

% Set time derivatives 
qstar = ns*b*c(1:n)*R*T./(1+(b*c(1:n)*R*T).^(t)).^(1/t);            
DqDt = k*(qstar - q );
DcDt(1) = 0.0;
DcDt(2:n) = Dl*D2cDz2(2:n) - v*DcDz(2:n) - ((1-epsilon)/(epsilon))*rho*DqDt(2:n);

% Concatenate vector of time derivatives
DyDt = [DcDt;DqDt];
end

[Torsten]

I corrected the code.

[Sona]

I cant thank you enough! 🧡

Answer 2

[Torsten]

The sink term

- (((1-eps)*rhop)/eps)*dqdt(i);

for the concentration of your gaseous species is always 0 because you specify dqdt(i) after dCdt(i) in the for-loop.

Thus the breakthrough is at the time instant

t_break = L/(u/eps)

[Ari Dillep Sai]

Hi Torsten.

Thanks for your reply.

I tried by specifying dqdt(i) first in the for loop but the nature of whole graph changed after doing that.

Can you please me in solving the equations which I have given in the attachment?

[Torsten]

I tried by specifying dqdt(i) first in the for loop but the nature of whole graph changed after doing that.

That's for sure, but the old version was wrong.

[Ari Dillep Sai]

Can you please let me know the mistakes in the code?

[Ari Dillep Sai]

Hi Torsten,

Can you help me in finding mistakes in my code with the respect to the boundary conditions and modelling equations?