simulation of adsorption with heat transfer

Hi, I try to simulate a TSA of 2 components with Matlab, for that I firstly try to simulate the adsorption step with the equations in the fields ''eq'' and ''eq2''. I have a problem with the quantity adsorbed, so if you have any time, could you help me please?

I use a main text onlyads.m that is going to change to add the others step (desorption, cooling and drying):

% Temperature Swing Adsorption Simulation
clc
clear
close all

%%%%%%%%%%%%%%%%%%%%%%%%%%%   GENERAL IMPUT %%%%%%%%%%%%%%%%%%%%%%%%%%%


% COLUMN IMPUT

D_column = 1.8;  % Column diameter (m)
L=6 ; %L bed (m)
S=L*pi()*D_column ; %S bed (m2)

% GAS IMPUT

cpmethane=2220; %(J/kgK)
cpco2=819; %(J/kgK)
cp_biogas=(0.5*cpmethane+0.5*cpco2); %(J/kgK)
lambdag=0.408; %(W/mK)
rho_gas=1.133; %Gas density (kg/m3)
R=8.3145 ; %J/molK
% butanol proprieties

c10=73.628*10^-3; % butanol concentration (kg/m3)
M1=74.12*10^-3; % masse molaire butanol (kg/mol)
carbon1 = 4; % carbon moles in the compounds (mol)
oxygen1 =1; % oxygen moles in the compounds (mol)
hydrogen1=5;% hydrogen moles in the compounds (mol)
V1=(15.9*carbon1+6.11*oxygen1+2.31*hydrogen1)*10^-6; %m3/mol

% toluene proprieties

c20=8.425*10^-3;  % toluene concentration (kg/m3)
M2=92.14*10^-3; % masse molaire toluene (kg/mol)
carbon2 = 6;% carbon moles in the compounds (mol)
oxygen2 =1;% oxygen moles in the compounds (mol)
hydrogen2=10;% hydrogen moles in the compounds (mol)
V2=(16.5*carbon2+5.48*oxygen2+1.98*hydrogen2)*10^-6;


%BPL activated carbon

eps_bed = 0.44;  % Interparticle voidage (m3voidage/m3bed)
eps_bead = 0.35;  % Intraparticle voidage (m3voidage/m3bead)
rho_bed = 480 ;%Bulk density (kg/m3)
d_particule = 2*10^-3; %Particule diameter (m)
r_microporous = 1*10^-6 ;%Microprous radius (m)
A_spe=1680; %Specific area (m-1)
lambda = 0.63; %Thermal conductivity (W/m.K)
tor=1.325; %Tortuosity
lambdas=555;  %(W/mK)
cp_bed=1560; % specific heat capacity (J/K kg)


Heat=2*lambdag/d_particule; % heat transfert coefficient fluid solid (W/m2/K)


%%%%%%%%%%%%%%%%%%%%%%%%%%%   ADSORPTION IMPUT %%%%%%%%%%%%%%%%%%%%%%%%%%%  


T_inlet = 293.15;  % Inlet temperature (K)
P_inlet = 101325;   %Inlet pressure (Pa)    
Q_inlet=1200/3600; % Feed flow (Nm3/s)
u=Q_inlet/S ; % vitesse superficielle (m/s)

%isotherm parameters
b01=9.98*10^-7;  %  Langmuir-Freundlich constants (mol/kg)
b02=1.88*10^-5; % Langmuir-Freundlich constants (mol/kg)
qm1=5.38; % mol/kg
qm2=5.52; % mol/kg
DeltaH1=58601; % J/mol
DeltaH2=50085; %J/mol
t1=0.5;
t2=0.33;


%increament for length bed
Nz=101;
z=linspace(0,L,Nz);
dz=z(2)-z(1);


%time
t=0:43200;

A_parameters = zeros(1,32);
A_parameters(1)=M1;
A_parameters(2)=M2;
A_parameters(3)=c10;
A_parameters(4)=c20;
A_parameters(5)=Nz;
A_parameters(6)=P_inlet;
A_parameters(7)=T_inlet;
A_parameters(8)=V1;
A_parameters(9)=V2;
A_parameters(10)=u;
A_parameters(11)=d_particule;
A_parameters(12)=eps_bead;
A_parameters(13)=R;
A_parameters(14)=tor;
A_parameters(15)=eps_bed;
A_parameters(16)=b01;
A_parameters(17)=b02;
A_parameters(18)=qm1;
A_parameters(19)=qm2;
A_parameters(20)=dz;
A_parameters(21)=lambdag;
A_parameters(22)=lambdag;
A_parameters(23)=cp_biogas;
A_parameters(24)=rho_gas;
A_parameters(25)=A_spe;
A_parameters(26)=Heat;
A_parameters(27)=DeltaH1;
A_parameters(28)=DeltaH2;
A_parameters(29)=cp_bed;
A_parameters(30)=rho_bed;
A_parameters(31)=t1;
A_parameters(32)=t2;



[c1ads,c2ads,q1ads,q2ads,Tgads,Tsads,y]=Adsorptionstep(A_parameters);

Diff = 0.9230

Diff = 0.0078

%In this code I use the function Adsorptionstep like you can see: 
function [c1ads,c2ads,q1ads,q2ads,Tgads,Tsads,y]=Adsorptionstep(A_parameters)

M1=A_parameters(1);
M2=A_parameters(2);
c10=A_parameters(3);
c20=A_parameters(4);
Nz=A_parameters(5);
P_inlet=A_parameters(6);
T_inlet=A_parameters(7);
V1=A_parameters(8);
V2=A_parameters(9);
u=A_parameters(10);
d_particule=A_parameters(11);
eps_bead=A_parameters(12);
R=A_parameters(13);
tor=A_parameters(14);
eps_bed=A_parameters(15);
b01=A_parameters(16);
b02=A_parameters(17);
qm1=A_parameters(18);
qm2=A_parameters(19);
dz=A_parameters(20);
lambdag=A_parameters(21);
lambdas=A_parameters(22);
cp_biogas=A_parameters(23);
rho_gas=A_parameters(24);
A_spe=A_parameters(25);
Heat=A_parameters(26);
DeltaH1=A_parameters(27);
DeltaH2=A_parameters(28);
cp_bed=A_parameters(29);
rho_bed=A_parameters(30);
t1=A_parameters(31);
t2=A_parameters(32);


t=0:43200;
%initial conditions
ICA=zeros(1,Nz); %for PDE-1
ICB=zeros(1,Nz);%for PDE-2
ICC=zeros(1,Nz);%for ODE-1
ICD=zeros(1,Nz);%for ODE-2
ICE=293.15*ones(1,Nz);%for ODE-2
ICF=400*ones(1,Nz);%for ODE-2
IC=[ICA,ICB,ICC,ICD,ICE,ICF];

[y1,y2]=fracmolaire(M1,M2,c10,c20,R,T_inlet,P_inlet);
p1=y1*P_inlet;
p2=y2*P_inlet;

[y1bis,y2bis]=fracmolairebis(y1,y2);
D1= Diffusivite(T_inlet,y2bis,M1,M2,P_inlet,V1,V2,u,d_particule,eps_bead);
D2=Diffusivite(T_inlet,y1bis,M1,M2,P_inlet,V1,V2,u,d_particule,eps_bead);

LDF1=Mass_transfert_coeff(d_particule,R,T_inlet,M1,D1,tor,eps_bead);
LDF2=Mass_transfert_coeff(d_particule,R,T_inlet,M1,D2,tor,eps_bead);
% ODE SOLVER
parameters = zeros(1,34);
parameters(1)=eps_bed ;
parameters(2)=u ;
parameters(3)=LDF1;
parameters(4)= LDF2;
parameters(5)=b01 ;
parameters(7)=b02;
parameters(6)=qm2 ;
parameters(8)=qm1;
   parameters(9)=c10;
    parameters(10)=c20;
    parameters(11)=rho_bed;
    parameters(12)=Nz;
    parameters(13)=dz;
    parameters(14)=D1;
    parameters(15)=D2;
    parameters(16)=lambdag;
    parameters(17)=lambdas;
    parameters(18)=rho_gas;
    parameters(19)=cp_biogas;
    parameters(20)=A_spe;
    parameters(21)=Heat;
    parameters(22)=DeltaH1;
    parameters(23)=DeltaH2;
    parameters(24)=R;
    parameters(25)=p1;
    parameters(26)=p2;
    parameters(27)=T_inlet;
    parameters(28)=cp_bed; 
    parameters(29)=t1;
    parameters(30)=t2;
    parameters(31)=eps_bead;
    parameters(32)=M1;
    parameters(33)=M2;
    parameters(34)=P_inlet;
    
[t_sol,y]=ode15s(@f2,t,IC,[],parameters);
    

% extract value
c1ads=y(:,1:Nz);
c2ads=y(:,Nz+1:2*Nz);
q1ads=y(:,2*Nz+1:3*Nz);
q2ads=y(:,3*Nz+1:4*Nz);
Tgads=y(:,4*Nz+1:5*Nz);
Tsads=y(:,5*Nz+1:6*Nz);

%reinput BC
c1ads(:,1)=c10;
c2ads(:,1)=c20;
Tgads(:,1)=T_inlet;
c1ads(:,end)=(4*c1ads(:,end-1)-c1ads(:,end-2))./3;
c2ads(:,end)=(4*c2ads(:,end-1)-c2ads(:,end-2))./3;
q1ads(:,end)=(4*q1ads(:,end-1)-q1ads(:,end-2))./3;
q2ads(:,end)=(4*q2ads(:,end-1)-q2ads(:,end-2))./3;
Tgads(:,end)=(4*Tgads(:,end-1)-Tgads(:,end-2))./3;
Tsads(:,end)=(4*Tsads(:,end-1)-Tsads(:,end-2))./3;
end 

%that use fracmolaire fracmolairebis diffusivite and mass_transfer_coeff to calculate parameters useful for the ODE solver :
function [frac1,frac2] = fracmolaire(M1,M2,c10,c20,R,T_inlet,P_inlet)
frac1=c10*R*T_inlet/(M1*P_inlet); % Assumptions P and T constant to simplify calculs
frac2=c20*R*T_inlet/(M2*P_inlet);
end
function [frac1bis,frac2bis] = fracmolairebis(y1,y2)
frac1bis= y1/(y2);
frac2bis= y2/(y1);
end
function Diff = Diffusivite(T_inlet,ybis,M1,M2,P_inlet,V1,V2,u,d_particule,eps_bead)
Dmol12=10^(-3)*T_inlet^(1.75)*(1/M1+1/M2)^(1/2)/(P_inlet*(V1^(1/3)+V2^(1/3))^2);
Dmol12bis=1/(ybis/Dmol12);
Diff=0.73*Dmol12bis+u*d_particule/2/(eps_bead*(1+9.49*eps_bead*Dmol12bis/(2*u*d_particule/2)))
end
function kLDF=Mass_transfert_coeff(d_particule,R,T0,M,Diff,tort,eps_bead)
       Dk=0.97*d_particule/2*(T0/M)^0.5;
       Dp=1/(1/Dk+1/Diff) ;
       Deff=Dp*eps_bead/tort ;
       Delta=(2+1)*(2+3);
       kLDF=Delta*Deff/(d_particule/2)^2;
end 
%The ODE solver is made with the function f2:
% fonction dy/dt
function dydt=f2(t,y,parameters)

    eps_bed = parameters(1);
    u = parameters(2);
    LDF1 = parameters(3);
    LDF2 = parameters(4);
    b01 = parameters(5);
    b02 = parameters(6);
    qm1 = parameters(7);
    qm2 = parameters(8);
    c10 = parameters(9);
    c20 = parameters(10);
    rho_bed = parameters(11);
    Nz = parameters(12);
    dz = parameters(13);
    D1 = parameters(14);
    D2 = parameters(15);
    lambdag = parameters(16);
    lambdas = parameters(17);
    rho_gas = parameters(18);
    cp_biogas = parameters(19);
    A_spe = parameters(20);
    Heat = parameters(21);
    DeltaH1 = parameters(22);
    DeltaH2 = parameters(23);
    R = parameters(24);
    p1 = parameters(25);
    p2 = parameters(26);
    T_inlet = parameters(27);
    cp_bed=parameters(28); 
    t1=parameters(29);
    t2=parameters(30);
    eps_bead=parameters(31);
    M1=parameters(32);
    M2=parameters(33);
    P_inlet=parameters(34);
    

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

dTgdt=zeros(Nz,1);
dTsdt=zeros(Nz,1);


%define value
c1=y(1:Nz); %concentration of the 1st component in mol/m3
c2=y(Nz+1:2*Nz); %concentration of the second component in mol/m3
q1=y(2*Nz+1:3*Nz);%quantity of component 1 adsorbed in mol/kg
q2=y(3*Nz+1:4*Nz);%quantity of component 2 adsorbed in mol/kg
Tg=y(4*Nz+1:5*Nz);%Gas temperature in K
Ts=y(5*Nz+1:6*Nz); %Solid temperature in K




%BC
c1(1)=c10;
c2(1)=c20;
c1(end)=(4*c1(end-1)-c1(end-2))./3;
c2(end)=(4*c2(end-1)-c2(end-2))./3;
Tg(1)=T_inlet;
Ts(1)=400;
Tg(end)=(4*Tg(end-1)-Tg(end-2))./3;
Ts(end)=(4*Ts(end-1)-Ts(end-2))./3;
b1(1)=b01*exp(-DeltaH1/R*Ts(1));
b2(1)=b02*exp(-DeltaH1/R*Ts(1));
p1(1)=c10*R*Tg(1)/(M1*P_inlet)*P_inlet;
p2(1)=c20*R*Tg(1)/(M2*P_inlet)*P_inlet;

%interior

for i=2:Nz-1
    p1(i)=c1(i)*R.*Tg(i)/(M1*P_inlet)*P_inlet; % Assumptions P  constant to simplify calculs
    p2(i)=c2(i)*R.*Tg(i)/(M2*P_inlet)*P_inlet;
    b1(i)=b01*exp(-DeltaH1/R*Ts(i));
    b2(i)=b02*exp(-DeltaH2/R*Ts(i));
    q1star(i)=qm1.*(b1(i).*p1(i))^t1./(1+(b1(i).*p1(i)+b2(i).*p2(i))^t1);
    q2star(i)=qm2.*(b2(i).*p2(i))^t2./(1+(b1(i).*p1(i)+b2(i).*p2(i))^t2);
    dq1dt(i)=LDF1.*(q1star(i)-q1(i));
    dq2dt(i)=LDF2.*(q2star(i)-q2(i));
    q1ads(i, :) = q1;
    q2ads(i, :) = q2;
    dc1dz(i)=(c1(i+1)-c1(i-1))./2./dz;
    dc2dz(i)=(c2(i+1)-c2(i-1))./2./dz;
    dcc1dz(i)=(c1(i+1)-2*c1(i)+c1(i-1))./2./dz;
    dcc2dz(i)=(c2(i+1)-2*c2(i)+c2(i-1))./2./dz;
    dc1dt(i)= D1*dcc1dz(i)-u*dc1dz(i)-rho_bed*((1-eps_bead)./eps_bead).*dq1dt(i);
    dc2dt(i)=D2*dcc2dz(i)-u*dc2dz(i)-rho_bed*((1-eps_bead)./eps_bead).*dq2dt(i);
    dTgdz(i)=(Tg(i+1)-Tg(i-1))./2./dz;
    dTsdz(i)=(Ts(i+1)-Ts(i-1))./2./dz;
    dTTgdz(i)=(Tg(i+1)-2*Tg(i)+Tg(i-1))./2./dz;
    dTTsdz(i)=(Ts(i+1)-2*Ts(i)+Ts(i-1))./2./dz;
    dTgdt(i)= lambdag/(rho_gas*cp_biogas*eps_bead)*dTTgdz(i)-u/eps_bead*dTgdz(i)-(1-eps_bed)/(rho_gas*cp_biogas*eps_bead)*A_spe*Heat*(Tg(i)-Ts(i));
    dTsdt(i)=lambdas/(rho_bed*cp_bed)*dTTsdz(i)-(1-eps_bed)*A_spe*Heat*(Tg(i)-Ts(i))/(rho_bed*cp_bed)-(DeltaH1*dq1dt(i)+DeltaH2*dq2dt(i))/(rho_bed*cp_bed);
end
dydt=[dc1dt;dc2dt;dq1dt;dq2dt;dTgdt;dTsdt];
end

So my problem is that i obtain q1ads and q2ads as matrix of zeros like if there is nothing that is adsorbed. Is there a problem with my initial conditions?

Sorry for the english errors i am not an english speaker

[ProblemSolver]

@Saradsorption- Please use the code format, commenting and proper Indentation about the functions. Use these guidelines:

https://www.mathworks.com/matlabcentral/answers/6200-tutorial-how-to-ask-a-question-on-answers-and-get-a-fast-answer

[Saradsorption]

Thank you, I didn't know

Answer 1

[ProblemSolver]

Hello @Saradsorption,

From the code that you have provided shows that the variables q1ads and q2ads are not assigned properly. I also see that you have then overwritten these two variables in the function "f2". However, when you are using ODE Solver, you should assign the correct initial conditions before calling the ODE solver. Therefore you should do something like this:

function [c1ads,c2ads,q1ads,q2ads,Tgads,Tsads,y]=Adsorptionstep(A_parameters)
    % ...
    
    % ODE SOLVER
    parameters = zeros(1,34);
    % ...
    [t_sol,y]=ode15s(@f2,t,IC,[],parameters);
    
    % extract value
    c1ads=y(:,1:Nz);
    c2ads=y(:,Nz+1:2*Nz);
    q1ads=y(:,2*Nz+1:3*Nz);
    q2ads=y(:,3*Nz+1:4*Nz);
    Tgads=y(:,4*Nz+1:5*Nz);
    Tsads=y(:,5*Nz+1:6*Nz);
    
    % ...
end

I hope this helps.

[Saradsorption]

Thank you for youre quick response, by changing these errors and by changing an error on the calculation on the adsorption constants b01 and b02 with a positive parameters inside the exp (it was negatif), i have this message and I don't know why?

Warning: Failure at t=0.000000e+00. Unable to meet integration tolerances without reducing the step size below the

smallest value allowed (7.905050e-323) at time t.

> In ode15s (line 653)

In Adsorptionstep (line 96)

In onlyads (line 123)

[Torsten]

    dTTgdz(i)=(Tg(i+1)-2*Tg(i)+Tg(i-1))./2./dz;
    dTTsdz(i)=(Ts(i+1)-2*Ts(i)+Ts(i-1))./2./dz;

These approximations for the 2nd spatial derivatives of Tg and Ts are not correct. The division should be by dz^2, not by 2*dz.

And my guess is that there are some more errors in your code that make the integrator fail right at the beginning (that's what the error message says).

You should look at what you return to the integrator. Maybe there are some Inf or NaN values in the array dydt.