simulation of adsorption with heat transfer

10 次查看（过去 30 天）
Saradsorption
Saradsorption 2023-7-10
编辑： Torsten 2023-7-10
Ran in:
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

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

回答（1 个）

ProblemSolver
ProblemSolver 2023-7-10
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.
  2 个评论
Saradsorption
Saradsorption 2023-7-10
移动：Torsten 2023-7-10
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
Torsten 2023-7-10
编辑：Torsten 2023-7-10
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.

请先登录，再进行评论。

类别

在 Help Center 和 File Exchange 中查找有关 Programming 的更多信息

标签

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by