simulation of adsorption with heat transfer
18 次查看(过去 30 天)
显示 更早的评论
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);
%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
2 个评论
ProblemSolver
2023-7-10
@Saradsorption- Please use the code format, commenting and proper Indentation about the functions. Use these guidelines:
回答(1 个)
ProblemSolver
2023-7-10
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 个评论
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.
另请参阅
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!