Problem solving coupled PDE / Adsorption

Hello, I'm trying to model a packed bed to adsorps H2O on Activated Alumina. I've worked a lot researching but I'm not getting proper results.

I've watched many papers and codes, but they doesn't use mole fraction as unit of concentration.

The result should be around 2.59 mol/kg adsorbed

Thank you in advance.

The code that I'm working at is this:

L=0.4                                 %Large m
R=8.314;                              %J/molK
P0=493000;                            %Pascal
T0=298;                               %K
yWater=0.3/100;                       % mole fraction, air composed by 0.3% of Water 
pco2=0.1972                           %To use on function, alumina doesnt adsorb co2 on this layer
psWater=3.169                         %kPa

%% Activated Alumina Isotherm temperature dependent from Zhang 2013
nh2o=40548*exp(-2987/T0);
bco2=5.762*10^-2 *exp(191.57/T0);    %1/kPa
bh2o=7.504*exp(490.85/T0);           %1/kPa
qsat=0.0934*exp(1030.8/T0)          %qsaturation mol/kg
epss=0.26
rhos=820                             %Kg/m3

%% Velocity
rof=1.184;                           %kg/m3 %density at 298K
m=10;                                %inlet mass flowrate kg/s    Qdot=V*A  ROF=m/Qdot
A=pi/4 *1.5^2;                       %area m^2
Q=m/rof;                             %flow m^3
v=Q/A;                               %superficial velocity m/s
u=v/epss;                            %intersticial velocity m/s

%% Mass transfer
LDF1=1E-3;                           %1/s

%% Main code
Nz=101;
z=linspace(0,L,Nz);
dz=z(2)-z(1);

% Time step
t=0:14400;

% Initial condition
ICA=zeros(1,Nz);                    %Inlet gas concentration
ICB=zeros(1,Nz);                    %Initial bed concentration
IC=[ICA,ICB];                       %Vector

% Solver
[t y]=ode15s(@f,t,IC,[],Nz,yWater,P0,qsat,bh2o,bco2,pco2,LDF1,u,rhos,R,T0,epss,dz,nh2o,psWater);

% Extract value
c1=y(:,1:Nz);
q1=y(:,Nz+1:2*Nz);

% Re BC
c1(:,1)=yWater;
c1(:,end)=(4*c1(:,end-1)-c1(:,end-2))./3;
q1(:,end)=(4*q1(:,end-1)-q1(:,end-2))./3;

function dydt=f(t,y,Nz,yWater,P0,qsat,bh2o,bco2,pco2,LDF1,u,rhos,R,T0,epss,dz,nh2o,psWater);

dydt=zeros(length(y),1);
dc1dt=zeros(Nz,1);
dq1dt=zeros(Nz,1);

%Assign values
c1=y(1:Nz);
q1=y(Nz+1:2*Nz);

%BC
c1(1)=yWater; 
c1(end)=(4*c1(end-1)-c1(end-2))./3;
q1(end)=(4*q1(end-1)-q1(end-2))./3;

%interior
for i=2:Nz-1
    dc1dz(i)=(c1(i+1)-c1(i-1)) ./2./dz;      %Centered
    
    ph2o(i)=c1(i)*P0/1000;                   %Mole fraction to kPa
    
    q1star(i)=qsat*bh2o.*ph2o(i)*exp(nh2o*ph2o(i)/psWater)/(1+bh2o.*ph2o(i)+bco2*pco2); %Extended Langmuir
    
    dq1dt(i)=LDF1.*(q1star(i)-q1(i));        %LDF
    
    dc1dt(i)=-u.*dc1dz(i)-((rhos*R*T0./P0).*((1-epss)./epss).*dq1dt(i));
end
dydt=[dc1dt;dq1dt];
end

[Torsten]

Are you sure about the velocity u ? Never seen an adsorption with a flow velocity for the adsorbent of 18.4 m/s .

[Jose Adones]

Hi @Torsten, I'm trying to recreate the paper "The effect of air purification on liquid air energy storage – An analysis from molecular to systematic modelling", they use a bed of 1.5 m of diameter and a inlet massflow of 10 kg/s.

It bothers me too, but it seems to be that way.

If you have any different idea, let me know, I'll be very greatful.

Thank you.

Answer 1

[Torsten]

Your code works if you change the centered scheme for the first derivatives to the usual upwind scheme:

dc1dz(i)=(c1(i)-c1(i-1))/dz;      %Upwind

instead of

dc1dz(i)=(c1(i+1)-c1(i-1)) ./2./dz;      %Centered

[Jose Adones]

Yes it does, the code run, but the adsorbed amount on activated alumina is wrong, should be 2.6 mol/kg not 6.75 mol/kg.

I don't know where the error can be, I think that is on the isotherm and the partial pressures.

These are the equations, pressure its on Mpa

[Torsten]

Which variable are you referring to when you say it is 6.75 mol/kg instead of 2.6 mol/kg ?

[Jose Adones]

"q1" as adsorbed amount on solid

[Torsten]

Then I suggest you calculate q1star for pH2O = 3e-3*P0/1000.

If you don't get 2.6 mol/kg, but 6.75 mol/kg, your equation for q1star is wrong.

T0=298;                               %K
P0=493000;                            %Pascal
pco2=0.1972;                          %To use on function, alumina doesnt adsorb co2 on this layer
psWater=3.169;                        %kPa
%% Activated Alumina Isotherm temperature dependent from Zhang 2013
nh2o=40548*exp(-2987/T0);
bco2=5.762*10^-2 *exp(191.57/T0);    %1/kPa
bh2o=7.504*exp(490.85/T0);           %1/kPa
qsat=0.0934*exp(1030.8/T0);          %qsaturation mol/kg

y_water = 0:0.0001:0.005;                 % mole fraction, air composed by 0.3% of Water
ph2o=y_water*P0/1000;
q1star=qsat*bh2o.*ph2o.*exp(nh2o*ph2o/psWater)./(1+bh2o.*ph2o+bco2*pco2); %Extended Langmuir
plot(y_water,q1star)
grid on

[Jose Adones]

It seems like the isotherm is correct, these are from doi/10.1021/je0201267.

I don't know what else could be wrong, "The ambient air is composed of N2 (77.74%), O2 (21.92%), Ar (1%), H2O (0.3%) and CO2 (0.04%);" and the adsorption pressure and temperature are 493 kPa and 298 K, I will check again, I don't think mr. Wang miscalculated.

[Jose Adones]

Hi Torsten, what if I want to try a variable velocity like this?

Which variable should I solve to discretize on eqn 5?

Thanks in advance.

[Torsten]

Is it really necessary to consider the change of velocity due to the mass adsorbed ? How much percent of the inflow is adsorbed so that you think this could be important ?

[Jose Adones]

Less than 2%, do you think that I should consider velocity as a constant?

[Torsten]

Try a linear increase of velocity over the length of the reactor that reflects this decrease in mass. I don't think you will get an effect on the results. Maybe if temperature changes drastically because of the adsorption, the effect on fluid velocity might become more important.

Answer 2

[Jose Adones]

Hi @Torsten, after sometime I finally made a functional adsorption code with heat transfer, I've followed a lot of your answers here.

But I have a question, I'm using langmuir for calculate qstar, but langmuir is useless when the RH exceeds 54%, so I have an excess surface work isotherm as follows:

T=298;
R=8.314;
ps=1.279;
Mw=18;                          %g/mol
%Langmuir parameters
k1 = 0.002082;                  % mol/g
k2 = 3.851e-5;                  % 1/K
k3 = 0.083;                     % 1/(mol/m3)
k4 = 1659;                      % K
k5 = 1.073;
k6 = 1.049;                     % K
qm = k1 + k2 * T;
b = k3 * exp(k4 / T);
n1 = k5 + k6 / T;

p=linspace(0,ps,1000);      %Concentration mol/m3
qstar = (qm * b * p.^n1) ./ (1 + b * p.^n1);
figure(1)
plot(p,qstar)

hold on
%ESW == qstar >54% RH
esw = -0.0386*(log(abs(R*T*log(p/ps)))-13.3)/Mw;
plot(p/ps,esw)
xlim([0, 1]);
ylim([0, 0.020]);
xlabel('Relative humidity')
ylabel('Adsorbed amount mol/g')
hold off

How can use this "double" isotherm on my code?

Thank you in advance

cFeed = 0.523;          % Concentration at the inlet node in mol/m3
Taire = 293;            % Air Temperature at the inlet in K
L = 0.3;                % Column length in m
t0 = 0;                 % Initial time in s            
tf = 50000;             % Final time in s
dt = 1;                 % Time step
z = [0:0.001:L].';      % Mesh generation
t = [t0:dt:tf];         % Time vector
n = numel(z);           % Mesh size

%%%%%% Initial conditions and boundary conditions %%%%%%%$%
c0 = zeros(n,1);        % c = 0 at t = 0 for all z
c0(1) = cFeed;          % c = cFeed at z = 0 and t > 0
q0 = 1E-5*ones(n,1);    % t = 0, q = 0 for all z
T0 = Taire*ones(n,1);
T0(1) = Taire;
Tw0 = Taire*ones(n,1);
y0 = [c0 ; q0 ; T0 ; Tw0];         

%%%%% Solving using the ODE15S solver %%%%%%
[Tiempo, Y] = ode15s(@(t,y) MyFun(t,y,z), t, y0);

%%%%%% Plots %%%%%
figure(2)
plot(Tiempo/60, Y(:,n))
legend('c/cFeed')
xlabel('Time min')
ylabel('Concentration mol/m3')

figure(3)
plot(Tiempo/60, Y(:,2*n))
legend('q')
xlabel('Time min')
ylabel('Adsorbed amount mol/g')

figure(4)
plot(Tiempo/60, Y(:,3*n))
legend('Air temperature')
xlabel('Time in minutes')
ylabel('Temperature in Kelvin')

figure(5)
plot(Tiempo/60, Y(:,4*n))
title('Wall temperature with atmospheric air at 293K')
legend('Wall temperature')
xlabel('Time in minutes')
ylabel('Temperature in Kelvin')

function DyDt = MyFun(t, y, z)
    n=numel(z);
    epsilon = 0.37;                 % Bed porosity
    volflow = 1.5e-4;               % Flow rate in m^3/s
    bedarea = pi/4 * 0.033^2;       % Bed area in m^2
    v = volflow / bedarea;          % Superficial velocity in m/s
    rho = 1100;                     % Particle density in kg/m^3
    Dl = 2.316e-3;                  % Axial dispersion coefficient in m^2/s
    k = 5e-4;                       % Mass transfer coefficient
    T = 293;                        % Temperature in K
    R = 8.314;                      %J/mol K
    
    
    % Langmuir isotherm parameters
    k1 = 0.002082;                  % mol/g
    k2 = 3.851e-5;                  % 1/K
    k3 = 0.083;                     % 1/(mol/m3)
    k4 = 1659;                      % K
    k5 = 1.073;
    k6 = 1.049;                     % K
    
    qm = k1 + k2 * T;
    b = k3 * exp(k4 / T);
    n1 = k5 + k6 / T;
%% Heat Transfer Calculation
    Tatm = T;
    Kl = 2.766;                     %W/m K
    alfa = 0.52;                    %Dimensionless
    rog = 1.159;                    %kg/m3
    rob = 690;                      %kg/m3
    dh = 62700;                     %j/mol
    hi = 1.5;                       %W/m2 s k
    ho = 4.2;                       %W/m2 s k
    Rbi = 0.033;                    %m
    Rbo = 0.03323;                  %m
    cps = 920;                      %j/kg K
    cpa = 2392.6;                   %j/kg K
    Ma = 27.8732/1000;              %kg/mol
    row = 7700;                     
    cpw = 500;
    cpg= 1047;
    
%%
    % Creating vectors for the variables
    c = zeros(n,1);
    q = zeros(n,1);
    T= zeros(n,1);
    Tw= zeros(n,1);
    DcDt = zeros(n,1);
    DqDt = zeros(n,1);
    DTDt = zeros(n,1);
    DTwDt = zeros(n,1);
    zhalf = zeros(n-1,1);
    DcDz = zeros(n,1);
    D2cDz2 = zeros(n,1);
    DTDz = zeros(n,1);
    D2TDz2 = zeros(n,1);
    DyDt = zeros(4*n,1);
    
    c = y(1:n);
    q = y(n+1:2*n);
    T = y(2*n+1:3*n);
    Tw = y(3*n+1:4*n);
    

    % Interior points of the mesh
    zhalf(1:n-1) = (z(1:n-1) + z(2:n)) / 2;

    % Calculation of 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));
    
    DTDz(2:n) = (T(2:n) - T(1:n-1)) ./ (z(2:n) - z(1:n-1));
    D2TDz2(2:n-1) = ((T(3:n) - T(2:n-1)) ./ (z(3:n) - z(2:n-1)) - (T(2:n-1) - T(1:n-2)) ./ (z(2:n-1) - z(1:n-2))) ./ (zhalf(2:n-1) - zhalf(1:n-2));
    
    % Calculation of D2cDz2 at z=L for the boundary condition dc/dz = 0
    D2cDz2(n) = -1.0 / (z(n) - zhalf(n-1)) * DcDz(n);
    
    D2TDz2(n) = -1.0 / (z(n) - zhalf(n-1)) * DTDz(n);

    % Calculation of qstar (adsorption capacity)
    qstar = (qm * b * (c(1:n)).^n1) ./ (1 + b * (c(1:n)).^n1);
    DqDt = (k * (qstar - q));
    DcDt(1) = 0.0;
    DcDt(2:n) = (Dl * D2cDz2(2:n) - v/epsilon * DcDz(2:n)) - ((1 - epsilon) / epsilon) * rho * 1000*DqDt(2:n);

    % Calculation of temperature
    DTDt(1)=0.0;
    DTDt(2:n) = (Kl*D2TDz2(2:n) - rog*cpg*epsilon*v*DTDz(2:n) - rob*dh*DqDt(2:n)- 2*hi/Rbi *(T(2:n)-Tw(2:n)))...
         ./ (alfa*rog*cpg + rob*cps + rob*q(2:n)*cpa*Ma);
    DTwDt(1)=0;
    DTwDt(2:n) = (2*pi*Rbi*hi*(T(2:n)-Tw(2:n)) -2*pi*Rbo*ho*(Tw(2:n)-Tatm))./ (row*cpw*(pi*(Rbo^2 - Rbi^2)));
    % Concatenation of the temporal derivative vectors
    DyDt = [DcDt; DqDt; DTDt; DTwDt];        
end

[Torsten]

A simple if-statement to make qstar depend on relative humidity does not work ?

[Jose Adones]

I tried this, but doesn't work

 qstar = zeros(n, 1);
 for i = 1:n
    if c(i) > 0.6912   %  54% of RH 
        ps=1.279;          %saturation mol/m3
        qstar(i) = -0.0386 * (log(abs(R * T * log(c(i) / ps))) - 13.3) / 18;
    else
        qstar(i) = (qm * b * (c(i))^n1) / (1 + b * (c(i))^n1);
    end
 end
 
 

Expected result for

 cFeed= 0.8456; %mol/m3

[Torsten]

Be careful with the variable "T" in "MyFun". First you use it as a fixed scalar, then you turn it into the array for the fluid temperature. I corrected it, but because of the if-statment, the runtime was too long for MATLAB online to see if it works. Test it with your own MATLAB version.

My guess is that qm, b and n1 don't have to be calculated with the fixed temperature, but with the actual temperature in the packed bed. You should correct this in your code (although temperature differences seem to be tiny).

cFeed = 0.8456;         % Concentration at the inlet node in mol/m3
Taire = 293;            % Air Temperature at the inlet in K
L = 0.3;                % Column length in m
t0 = 0;                 % Initial time in s            
tf = 50000;             % Final time in s
dt = 1;                 % Time step
z = [0:0.001:L].';      % Mesh generation
t = [t0:dt:tf];         % Time vector
n = numel(z);           % Mesh size
%%%%%% Initial conditions and boundary conditions %%%%%%%$%
c0 = zeros(n,1);        % c = 0 at t = 0 for all z
c0(1) = cFeed;          % c = cFeed at z = 0 and t > 0
q0 = 1E-5*ones(n,1);    % t = 0, q = 0 for all z
T0 = Taire*ones(n,1);
T0(1) = Taire;
Tw0 = Taire*ones(n,1);
y0 = [c0 ; q0 ; T0 ; Tw0];         
%%%%% Solving using the ODE15S solver %%%%%%
[Tiempo, Y] = ode15s(@(t,y) MyFun(t,y,z), t, y0);
%%%%%% Plots %%%%%
figure(1)
plot(Tiempo/60, Y(:,n))
legend('c/cFeed')
xlabel('Time min')
ylabel('Concentration mol/m3')

figure(2)
plot(Tiempo/60, Y(:,2*n))
legend('q')
xlabel('Time min')
ylabel('Adsorbed amount mol/g')

figure(3)
plot(Tiempo/60, Y(:,3*n))
legend('Air temperature')
xlabel('Time in minutes')
ylabel('Temperature in Kelvin')

figure(4)
plot(Tiempo/60, Y(:,4*n))
title('Wall temperature with atmospheric air at 293K')
legend('Wall temperature')
xlabel('Time in minutes')
ylabel('Temperature in Kelvin')

function DyDt = MyFun(t, y, z)
    n=numel(z);
    epsilon = 0.37;                 % Bed porosity
    volflow = 1.5e-4;               % Flow rate in m^3/s
    bedarea = pi/4 * 0.033^2;       % Bed area in m^2
    v = volflow / bedarea;          % Superficial velocity in m/s
    rho = 1100;                     % Particle density in kg/m^3
    Dl = 2.316e-3;                  % Axial dispersion coefficient in m^2/s
    k = 5e-4;                       % Mass transfer coefficient
    Temp = 293;                     % Temperature in K
    R = 8.314;                      %J/mol K
    
    
    % Langmuir isotherm parameters
    k1 = 0.002082;                  % mol/g
    k2 = 3.851e-5;                  % 1/K
    k3 = 0.083;                     % 1/(mol/m3)
    k4 = 1659;                      % K
    k5 = 1.073;
    k6 = 1.049;                     % K
    
    qm = k1 + k2 * Temp;
    b = k3 * exp(k4 / Temp);
    n1 = k5 + k6 / Temp;
%% Heat Transfer Calculation
    Tatm = Temp;
    Kl = 2.766;                     %W/m K
    alfa = 0.52;                    %Dimensionless
    rog = 1.159;                    %kg/m3
    rob = 690;                      %kg/m3
    dh = 62700;                     %j/mol
    hi = 1.5;                       %W/m2 s k
    ho = 4.2;                       %W/m2 s k
    Rbi = 0.033;                    %m
    Rbo = 0.03323;                  %m
    cps = 920;                      %j/kg K
    cpa = 2392.6;                   %j/kg K
    Ma = 27.8732/1000;              %kg/mol
    row = 7700;                     
    cpw = 500;
    cpg= 1047;
    
%%
    % Creating vectors for the variables
    c = zeros(n,1);
    q = zeros(n,1);
    T= zeros(n,1);
    Tw= zeros(n,1);
    DcDt = zeros(n,1);
    DqDt = zeros(n,1);
    DTDt = zeros(n,1);
    DTwDt = zeros(n,1);
    zhalf = zeros(n-1,1);
    DcDz = zeros(n,1);
    D2cDz2 = zeros(n,1);
    DTDz = zeros(n,1);
    D2TDz2 = zeros(n,1);
    DyDt = zeros(4*n,1);
    
    c = y(1:n);
    q = y(n+1:2*n);
    T = y(2*n+1:3*n);
    Tw = y(3*n+1:4*n);
    
    % Interior points of the mesh
    zhalf(1:n-1) = (z(1:n-1) + z(2:n)) / 2;
    % Calculation of 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));
    
    DTDz(2:n) = (T(2:n) - T(1:n-1)) ./ (z(2:n) - z(1:n-1));
    D2TDz2(2:n-1) = ((T(3:n) - T(2:n-1)) ./ (z(3:n) - z(2:n-1)) - (T(2:n-1) - T(1:n-2)) ./ (z(2:n-1) - z(1:n-2))) ./ (zhalf(2:n-1) - zhalf(1:n-2));
    
    % Calculation of D2cDz2 at z=L for the boundary condition dc/dz = 0
    D2cDz2(n) = -1.0 / (z(n) - zhalf(n-1)) * DcDz(n);
    
    D2TDz2(n) = -1.0 / (z(n) - zhalf(n-1)) * DTDz(n);
    % Calculation of qstar (adsorption capacity)
    %qstar = (qm * b * (c(1:n)).^n1) ./ (1 + b * (c(1:n)).^n1);
    qstar = zeros(n, 1);
    for i = 1:n
      if c(i) > 0.6912   %  54% of RH 
        ps = 1.279;          %saturation mol/m3
        qstar(i) = -0.0386 * (log(abs(R * T(i) * log(c(i) / ps))) - 13.3) / 18;
      else
        qstar(i) = (qm * b * (c(i))^n1) / (1 + b * (c(i))^n1);
      end
    end
    DqDt = (k * (qstar - q));
    DcDt(1) = 0.0;
    DcDt(2:n) = (Dl * D2cDz2(2:n) - v/epsilon * DcDz(2:n)) - ((1 - epsilon) / epsilon) * rho * 1000*DqDt(2:n);
    % Calculation of temperature
    DTDt(1)=0.0;
    DTDt(2:n) = (Kl*D2TDz2(2:n) - rog*cpg*epsilon*v*DTDz(2:n) - rob*dh*DqDt(2:n)- 2*hi/Rbi *(T(2:n)-Tw(2:n)))...
         ./ (alfa*rog*cpg + rob*cps + rob*q(2:n)*cpa*Ma);
    DTwDt(1)=0;
    DTwDt(2:n) = (2*pi*Rbi*hi*(T(2:n)-Tw(2:n)) -2*pi*Rbo*ho*(Tw(2:n)-Tatm))./ (row*cpw*(pi*(Rbo^2 - Rbi^2)));
    % Concatenation of the temporal derivative vectors
    DyDt = [DcDt; DqDt; DTDt; DTwDt];        
end

[Jose Adones]

Maybe that's why I can't get proper results for desorption, I will check it, thank you Torsten

[Jose Adones]

It worked, Torsten. But I had to add a line or it crashes

options = odeset('RelTol', 1e-3, 'AbsTol', 1e-4);
[Tiempo, Y] = ode15s(@(t,y) MyFun(t,y,z), t, y0,options);

I will work on desorption in the next few days, I hope you can help me as before, if I run into any problem, thank you very much

[Jose Adones]

Hello @Torsten, I hope you are well. I corrected the code, it had an error in the air-bed temperature. I will leave it attached for those who need it, use doi 10.1016/j.cherd.2019.08.018 as guide.

However, I was never able to get proper results for desorption.

Hope you can help me, in case it doesn't work I'll not count desorption on my thesis. Thank you in advance.

Expected results:

tic
    cFeed = 1E-6;         % Concentration at the inlet node in mol/m3
    Taire = 423; 
    L = 0.3;              % Column length in m
    t0 = 0;               % Initial time in s            
    tf = 10000;            % Final time in s
    dt = 1;               % Time step
    z = [0:0.001:L].';     % Mesh generation
    t = [t0:dt:tf];       % Time vector
    n = numel(z);         % Mesh size

    % Initial conditions and boundary conditions
    c0 = 0.523*ones(n,1);   % c = 0 at t = 0 for all z
    c0(1) = cFeed;          % c = cFeed at z = 0 and t > 0
    q0 = 0.0123*ones(n,1);  % t = 0, q = 0 for all z
    T0 = 293*ones(n,1);
    T0(1) = 423;
    Tw0 = 298*ones(n,1);
    y0 = [c0 ; q0 ; T0 ; Tw0];         

    % Solving using the ODE15S solver
    [Time, Y] = ode15s(@(t,y) MyFun(t,y,z), t, y0);

figure(1)
    % Subplot of the first graph
    subplot(2, 1, 1);
    plot(Time, Y(:,n))
    hold on
    title('Outlet concentration')
    xlabel('Time in seconds')
    ylabel('Concentration')
    xlim([0 8000])
   
    % Subplot of the second graph
    subplot(2, 1, 2);
    plot(Time, Y(:,3*n), 'r')
    hold on
    title('Temperature during the Process')
    xlabel('Time in seconds')
    ylabel('Temperature in Kelvin')
    xlim([0 8000])

    toc

Elapsed time is 3.840592 seconds.

function DyDt = MyFun(t, y, z) 
    n=numel(z);
    epsilon = 0.37;                 % Bed porosity
    volflow = 4.3e-4;               % Flow rate in m^3/s
    bedarea = pi/4 * 0.033^2;       % Bed area in m^2
    v = volflow / bedarea;          % Superficial velocity in m/s
    rho = 1100;                     % Particle density in kg/m^3
    Dl = 3.995e-3;                  % Axial dispersion coefficient in m^2/s
    k = 5e-3;                       % Mass transfer coefficient
    Tfeed = 423;                    % Temperature in K
    R = 8.314;                      % J/mol K
    
    % Langmuir isotherm parameters
    k1 = 0.002082;                  % mol/g
    k2 = 3.851e-5;                  % 1/K
    k3 = 0.083;                     % 1/(mol/m3)
    k4 = 1659;                      % K
    k5 = 1.073;
    k6 = 1.049;                     % K
    
    qm = k1 + k2 * Tfeed;
    b = k3 * exp(k4 / Tfeed);
    n1 = k5 + k6 / Tfeed;
    
%% Heat Transfer Calculation
    Tatm = 293;
    Kl = 3.374;                     % W/m K
    alfa = 0.52;                    % Dimensionless
    rog = 1.159;                    % kg/m3
    rob = 690;                      % kg/m3
    dh = 62700;                     % J/mol
    hi = 1.5;                       % W/m2 s K
    ho = 4.2;                       % W/m2 s K
    Rbi = 0.033;                    % m
    Rbo = 0.03323;                  % m
    cps = 920;                      % J/kg K
    cpa = 2392.6;                   % J/kg K
    Ma = 18/1000;                   % kg/mol
    row = 7700;                     
    cpw = 500;
    cpg = 1047;
    
%%    
    % Creating vectors for variables
    c = zeros(n,1);
    q = zeros(n,1);
    T = zeros(n,1);
    Tw = zeros(n,1);
    DcDt = zeros(n,1);
    DqDt = zeros(n,1);
    DTDt = zeros(n,1);
    DTwDt = zeros(n,1);
    zhalf = zeros(n-1,1);
    DcDz = zeros(n,1);
    D2cDz2 = zeros(n,1);
    DTDz = zeros(n,1);
    D2TDz2 = zeros(n,1);
    DyDt = zeros(4*n,1);
    
    c = y(1:n);
    q = y(n+1:2*n);
    T = y(2*n+1:3*n);
    Tw = y(3*n+1:4*n);
    
    % Interior mesh points
    zhalf(1:n-1) = (z(1:n-1) + z(2:n)) / 2;

    % Calculation of 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));
    
    DTDz(2:n) = (T(2:n) - T(1:n-1)) ./ (z(2:n) - z(1:n-1));
    D2TDz2(2:n-1) = ((T(3:n) - T(2:n-1)) ./ (z(3:n) - z(2:n-1)) - (T(2:n-1) - T(1:n-2)) ./ (z(2:n-1) - z(1:n-2))) ./ (zhalf(2:n-1) - zhalf(1:n-2));
    
    % Calculation of D2cDz2 at z=L for the boundary condition dc/dz = 0
    D2cDz2(n) = -1.0 / (z(n) - zhalf(n-1)) * DcDz(n);
    
    D2TDz2(n) = -1.0 / (z(n) - zhalf(n-1)) * DTDz(n);

    % Calculation of qstar (adsorption capacity)
    ps=(0.61094*exp(17.625*(Tatm-273.15)/(Tatm-30.11))*1000/(R*Tfeed)); % mol/m3
    qstar = (qm * b * (c(1:n)).^n1) ./ (1 + b * ((c(1:n)).^n1).*(1-(c(1:n)./ps)).^0.03);
    DqDt = (k * (qstar - q));
    DcDt(1) = 0.0;
    DcDt(2:n) = (Dl * D2cDz2(2:n) - v/epsilon * DcDz(2:n)) - ((1 - epsilon) / epsilon) * rho * 1000*DqDt(2:n);

    % Temperature calculation
    DTDt(1)=0.0;
    DTDt(2:n) = (Kl*D2TDz2(2:n) - rog*cpg*v/epsilon *DTDz(2:n) + ((1-epsilon)/epsilon)*dh*rho*1000*DqDt(2:n)- 2*hi/Rbi *(T(2:n)-Tw(2:n)))...
    ./ (alfa*rog*cpg + ((1-epsilon)/epsilon)*rob*cps); % Luberti 2015
  
     
    DTwDt(1)=0;
    DTwDt(2:n) = (2*pi*Rbi*hi*(T(2:n)-Tw(2:n)) -2*pi*Rbo*ho*(Tw(2:n)-Tatm))./ (row*cpw*(pi*(Rbo^2 - Rbi^2)));
    % Concatenation of temporal derivatives vectors
    DyDt = [DcDt; DqDt; DTDt; DTwDt];        
end

[Torsten]

The water vapor concentration in the gas phase increases strongly according to the article while yours permanently decreases towards 0 because of the inlet concentration of 0. Are you sure your qstar adequately reflects the influence of temperature on the adsorption equilibrium ?

[Jose Adones]

Sorry for the delay in responding.

Is the same isotherm of the adsorption process, the same parameters, should I change something besides the initial conditions?

[Torsten]

As you can see in the last plot from your code, the water vapor concentration decreases in the reactor right from the beginning whereas in Fig. 5 of the article, it rises to about 1.5 mol/m^3 for ~1500 s.

I don't have access to the article, but you should check the design variables of the desorption again (volume flow and temperature of flushing gas (why do you compute ps with Tfeed and not with T ?), your adsorption isotherm, the factor 1000 in the expression "((1 - epsilon) / epsilon) * rho * 1000*DqDt(2:n)",...)

     cFeed = 1E-6;         % Concentration at the inlet node in mol/m3
    Taire = 423;
    L = 0.3;              % Column length in m
    t0 = 0;               % Initial time in s
    tf = 10000;            % Final time in s
    dt = 1;               % Time step
    z = [0:0.001:L].';     % Mesh generation
    t = [t0:dt:tf];       % Time vector
    n = numel(z);         % Mesh size

    % Initial conditions and boundary conditions
    c0 = 0.523*ones(n,1);   % c = 0 at t = 0 for all z
    c0(1) = cFeed;          % c = cFeed at z = 0 and t > 0
    q0 = 0.0123*ones(n,1);  % t = 0, q = 0 for all z
    T0 = 293*ones(n,1);
    T0(1) = 423;
    Tw0 = 298*ones(n,1);
    y0 = [c0 ; q0 ; T0 ; Tw0];

    % Solving using the ODE15S solver
    [Time, Y] = ode15s(@(t,y) MyFun(t,y,z), t, y0);
c = Y(:,1:n);
q = Y(:,n+1:2*n);
T = Y(:,2*n+1:3*n);
Tw = Y(:,3*n+1:4*n);
figure(1)
    % Subplot of the first graph
    subplot(2, 1, 1);
    plot(Time, c(:,end))
    hold on
    title('Outlet concentration')
    xlabel('Time in seconds')
    ylabel('Concentration')
    xlim([0 8000])

    % Subplot of the second graph
    subplot(2, 1, 2);
    plot(Time, T(:,end), 'r')
    hold on
    title('Temperature during the Process')
    xlabel('Time in seconds')
    ylabel('Temperature in Kelvin')
    xlim([0 8000])

    %figure(2)
    %plot(z,[q(1,:);q(50,:);q(end,:)])
    figure(3)
     plot(z,[c(1,:);c(50,:);c(end,:)])

    toc

Elapsed time is 9.826570 seconds.

function DyDt = MyFun(t, y, z)
    n=numel(z);
    epsilon = 0.37;                 % Bed porosity
    volflow = 4.3e-4;               % Flow rate in m^3/s
    bedarea = pi/4 * 0.033^2;       % Bed area in m^2
    v = volflow / bedarea;          % Superficial velocity in m/s
    rho = 1100;                     % Particle density in kg/m^3
    Dl = 3.995e-3;                  % Axial dispersion coefficient in m^2/s
    k = 5e-3;                       % Mass transfer coefficient
    Tfeed = 423;                    % Temperature in K
    R = 8.314;                      % J/mol K

    % Langmuir isotherm parameters
    k1 = 0.002082;                  % mol/g
    k2 = 3.851e-5;                  % 1/K
    k3 = 0.083;                     % 1/(mol/m3)
    k4 = 1659;                      % K
    k5 = 1.073;
    k6 = 1.049;                     % K

    qm = k1 + k2 * Tfeed;
    b = k3 * exp(k4 / Tfeed);
    n1 = k5 + k6 / Tfeed;

%% Heat Transfer Calculation
    Tatm = 293;
    Kl = 3.374;                     % W/m K
    alfa = 0.52;                    % Dimensionless
    rog = 1.159;                    % kg/m3
    rob = 690;                      % kg/m3
    dh = 62700;                     % J/mol
    hi = 1.5;                       % W/m2 s K
    ho = 4.2;                       % W/m2 s K
    Rbi = 0.033;                    % m
    Rbo = 0.03323;                  % m
    cps = 920;                      % J/kg K
    cpa = 2392.6;                   % J/kg K
    Ma = 18/1000;                   % kg/mol
    row = 7700;
    cpw = 500;
    cpg = 1047;

%%
    % Creating vectors for variables
    c = zeros(n,1);
    q = zeros(n,1);
    T = zeros(n,1);
    Tw = zeros(n,1);
    DcDt = zeros(n,1);
    DqDt = zeros(n,1);
    DTDt = zeros(n,1);
    DTwDt = zeros(n,1);
    zhalf = zeros(n-1,1);
    DcDz = zeros(n,1);
    D2cDz2 = zeros(n,1);
    DTDz = zeros(n,1);
    D2TDz2 = zeros(n,1);
    DyDt = zeros(4*n,1);

    c = y(1:n);
    q = y(n+1:2*n);
    T = y(2*n+1:3*n);
    Tw = y(3*n+1:4*n);

    % Interior mesh points
    zhalf(1:n-1) = (z(1:n-1) + z(2:n)) / 2;

    % Calculation of 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));

    DTDz(2:n) = (T(2:n) - T(1:n-1)) ./ (z(2:n) - z(1:n-1));
    D2TDz2(2:n-1) = ((T(3:n) - T(2:n-1)) ./ (z(3:n) - z(2:n-1)) - (T(2:n-1) - T(1:n-2)) ./ (z(2:n-1) - z(1:n-2))) ./ (zhalf(2:n-1) - zhalf(1:n-2));

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

    D2TDz2(n) = -1.0 / (z(n) - zhalf(n-1)) * DTDz(n);

    % Calculation of qstar (adsorption capacity)
    ps=(0.61094*exp(17.625*(Tatm-273.15)./(Tatm-30.11)).*1000./(R*T(1:n))); % mol/m3
    qstar = (qm * b * (c(1:n)).^n1) ./ (1 + b * ((c(1:n)).^n1).*(1-(c(1:n)./ps)).^0.03);
    DqDt = (k * (qstar - q));
    DcDt(1) = 0.0;
    DcDt(2:n) = (Dl * D2cDz2(2:n) - v/epsilon * DcDz(2:n)) - ((1 - epsilon) / epsilon) * rho * 1000*DqDt(2:n);

    % Temperature calculation
    DTDt(1)=0.0;
    DTDt(2:n) = (Kl*D2TDz2(2:n) - rog*cpg*v/epsilon *DTDz(2:n) + ((1-epsilon)/epsilon)*dh*rho*1000*DqDt(2:n)- 2*hi/Rbi *(T(2:n)-Tw(2:n)))...
    ./ (alfa*rog*cpg + ((1-epsilon)/epsilon)*rob*cps); % Luberti 2015


    DTwDt(1)=0;
    DTwDt(2:n) = (2*pi*Rbi*hi*(T(2:n)-Tw(2:n)) -2*pi*Rbo*ho*(Tw(2:n)-Tatm))./ (row*cpw*(pi*(Rbo^2 - Rbi^2)));
    % Concatenation of temporal derivatives vectors
    DyDt = [DcDt; DqDt; DTDt; DTwDt];
end

[Jose Adones]

1. I will check the desing variables
2. "why do you compute ps with Tfeed and not with T" Idk, I believe thats the problem
3. "your adsorption isotherm, the factor 1000 in the expression "((1 - epsilon) / epsilon) * rho * 1000*DqDt(2:n)",...)" That is because the isotherm is in mol/g and must be used in kg

Thank you Torsten