# Solve for pressure swing adsorption cycle for two gases in one bed.

115 vues (au cours des 30 derniers jours)
Samridh Sharma le 24 Juil 2023
Commenté : Samridh Sharma le 27 Juil 2023
Hello,
I am trying to solve for pressure swing adsorption cycle to see the final concentration achieved in filtered fluents. I have been able to find various resources shared by matlab user but I need to use the governing partial differential equation that has coupled parameters with concentration of other gas in the Langmuir adsorption isotherm (equation 2).
Equations:
Bondary conditions:
Initial conditions:
So, the solution from method of lines approach uses ode15s but it does not let me use array with column indices 2 or more. Therefore I need help to have my fixed single bed adsorption filter operate for two different gases and solve equation 3 for both spatial and temporal solutions.
'C' is concentration in fluid phase.
'q' is concentration in adsorbed phase.
clc
close all
clear all
L = 0.02; % Column length
t0 = 0; % Initial Time
tf = 20; % Final time
dt = 0.05; % Time step
z = 0:0.0005:L; %Mesh generation
t = t0:dt:tf;% Time vector
n = numel(z); % Size of mesh grid
cF = [0.1 1.5 5 10 50 100 250 500 1000]; % Diffferent feed concentration values
for ii=1:numel(cF)
cFeed=cF(ii);
[T, Y] = Main(cFeed);
plot(T,Y(:,n)/cF(9), 'linewidth', 2), hold all
end
function [T,Y]=Main(cFeed)
%System Set-up %
%Define Variables
%D = 1.85*10^-5; % Axial Dispersion coefficient
%v = 3.38*10^-3; % Superficial velocity
%epsilon = 0.47; % Voidage fraction
%k = 1.9*10^-6; % Mass Transfer Coefficient
%ap= 0.003; % radius of pellet
%b = 0.05; % Langmuir parameter
%qm = 14.5; % saturation capacity
% rho= 1970; % paricle denisty
%cFeed = 10; % Feed concentration
L = 0.02; % Column length
t0 = 0; % Initial Time
%tf = 20; % Final time
tf = 20; % Final time
dt = 0.05; % Time step
% z=[0:0.01:L]; %Mesh generation
z = [0:0.0005:L]; %Mesh generation
t = [t0:dt:tf];% Time vector
n = numel(z); % Size of mesh grid
%Initial Conditions / Vector Creation
c0 = zeros(n,1);
c0(1) = cFeed;
q0 = zeros(n,1); % t = 0, q = 0 for all z, this makes sense to me
y0 = [c0 ; q0]; % Appends conditions together
%ODE15S Solver
[T, Y] = ode15s(@(t,y) MyFun(t,y,z,n),t,y0);
end
function DyDt=MyFun(~, y, z, n)
% Defining Constants
D = 1.85*10^-5; % Axial Dispersion coefficient
v =3.38*10^-3 ; % Superficial velocity
epsilon = 0.47; % Voidage fraction
k = 1.9*10^-6; % Mass Transfer Coefficient
b = 0.05; % Langmuir parameter
qm = 14.5; % saturation capacity
rho= 1970; %particle denisty
% Variables being allocated zero vectors
c = zeros(n,1);
q = zeros(n,1);
DcDt = zeros(n,1);
DqDt = zeros(n,1);
DyDt = zeros(2*n,1);
zhalf = zeros(n-1,1);
DcDz = zeros(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) = (zhalf(i)*(c(i+1)-c(i))/(z(i+1)-z(i))-zhalf(i-1)*(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
% DcDt = 0 since c=cFeed remains constant over time and is already set as initial condition
% Standard setting for q
DqDt(1) = k*( q(1)- ((qm*b*c(1))/(1 + b*c(1))) );
DcDt(1) = 0.0;
% Set time derivatives in remaining points
for i=2:n
%Equation: dq/dt = k(q-q*) where q* = qm * (b*c)/(1+b*c)
DqDt(i) = k*( q(i)- ((qm*b*c(i))/(1 + b*c(i))) );
%Equation: dc/dt = D * d2c/dz2 - v*dc/dz - ((1-e)/(e))*dq/dt
DcDt(i) = D*D2cDz2(i) -v* DcDz(i) - ((1-epsilon)/(epsilon))*rho *DqDt(i);
end
% Concatenate vector of time derivatives
DyDt = [DcDt;DqDt];
end
##### 3 commentairesAfficher 1 commentaire plus ancienMasquer 1 commentaire plus ancien
Samridh Sharma le 24 Juil 2023
@Davide Masiello thanks for your remark, I have updated my question with more details.
Davide Masiello le 24 Juil 2023

Connectez-vous pour commenter.

### Réponse acceptée

Davide Masiello le 25 Juil 2023
Modifié(e) : Davide Masiello le 27 Juil 2023
Hi @Samridh Sharma, please see the modifications I made to the code and the resulting output.
Please let me know if you have further questions.
clear,clc
%% MAIN
% Parameters
L = 0.02; % Column length
t0 = 0; % Initial Time
tf = 20; % Final time
dt = 0.05; % Time step
t = t0:dt:tf; % Time vector
dz = 5e-5; % Mesh size
z = 0:dz:L; % Mesh generation
n = numel(z);
cF_1 = [0.1 1.5 5 10 50 100 250 500 1000]; % Different feed concentration values
cF_2 = 0.1*ones(size(cF_1));
cF = [cF_1;cF_2];
% Solution and plotting
for k = 1:size(cF,2)
name = sprintf('Feed concentration of compound #1 cFeed = %f.\n',cF(1,k));
figure
subplot(4,2,1)
plot(z,sol(k).y(1:n,1:5:end))
xlabel('z')
ylabel('c1')
title('c1(z) at different t')
subplot(4,2,2)
plot(z,sol(k).y(n+1:2*n,1:5:end))
xlabel('z')
ylabel('q1')
title('q1(z) at different t')
subplot(4,2,3)
plot(sol(k).x,sol(k).y(1:40:n,:))
xlabel('t')
ylabel('c1')
title('c1(t) at different z')
subplot(4,2,4)
plot(sol(k).x,sol(k).y(n+1:40:2*n,:))
xlabel('t')
ylabel('q1')
title('q1(t) at different z')
subplot(4,2,5)
plot(z,sol(k).y(2*n+1:3*n,1:5:end))
xlabel('z')
ylabel('c2')
title('c2(z) at different t')
subplot(4,2,6)
plot(z,sol(k).y(3*n+1:4*n,1:5:end))
xlabel('z')
ylabel('q2')
title('q2(z) at different t')
subplot(4,2,7)
plot(sol(k).x,sol(k).y(2*n+1:40:3*n,:))
xlabel('t')
ylabel('c2')
title('c2(t) at different z')
subplot(4,2,8)
plot(sol(k).x,sol(k).y(3*n+1:40:4*n,:))
xlabel('t')
ylabel('q2')
title('q2(t) at different z')
sgtitle(name)
disp('-------------------------------------------------------------------------------------------------')
end
-------------------------------------------------------------------------------------------------
-------------------------------------------------------------------------------------------------
-------------------------------------------------------------------------------------------------
-------------------------------------------------------------------------------------------------
-------------------------------------------------------------------------------------------------
-------------------------------------------------------------------------------------------------
-------------------------------------------------------------------------------------------------
-------------------------------------------------------------------------------------------------
-------------------------------------------------------------------------------------------------
%% Solver function
n = numel(z); % Size of mesh grid
h = diff(z); h = h(1);
% Initial Conditions
c1_0 = zeros(n,1);
c1_0(1) = cFeed(1);
q1_0 = zeros(n,1);
c2_0 = zeros(n,1);
c2_0(1) = cFeed(2);
q2_0 = zeros(n,1);
y0 = [c1_0 ; q1_0; c2_0 ; q2_0];
% Creating mass matrix
M = eye(4*n); % Mass matrix
M(1,1) = 0; % Algebraic equation for c1 at z = 0;
M(n,n) = 0; % Algebraic equation for c1 at z = L
M(2*n+1,2*n+1) = 0; % Algebraic equation for c2 at z = 0;
M(3*n,3*n) = 0; % Algebraic equation for c2 at z = L
opts = odeset('Mass',M,'MassSingular','yes');
% Solving
end
% Constants
D = [1.85e-5;1.85e-5]; % Axial Dispersion coefficient
v = 3.38e-3 ; % Superficial velocity
epl = 0.47; % Void fraction
k = [1.9e-6;1.9e-6]; % Mass Transfer Coefficient
b = [0.05;0.1]; % Langmuir parameter
qm = [14.5;14.5]; % Saturation capacity
% Variables allocation
c1 = y(1:n);
q1 = y(n+1:2*n);
c2 = y(2*n+1:3*n);
q2 = y(3*n+1:4*n);
% Langmuir term
c = [c1,c2];
q = [q1,q2];
LT = (b'.*c)./(1+c*b)-q;
% Concentration equations
dc1dz = zeros(n,1);
dc1dz(2:n-1) = (c1(3:n)-c1(2:n-1))/h;
d2c1dz2 = zeros(n,1);
d2c1dz2(2:n-1) = (c1(3:n)-2*c1(2:n-1)+c1(1:n-2))/h^2;
dc1dt = -v*dc1dz+D(1)*d2c1dz2-k(1)*qm(1)*(1/epl-1)*LT(:,1);
dc1dt(1) = c1(1)-cFeed(1);
dc1dt(n) = c1(n)-c1(n-1);
dc2dz = zeros(n,1);
dc2dz(2:n-1) = (c2(3:n)-c2(2:n-1))/h;
d2c2dz2 = zeros(n,1);
d2c2dz2(2:n-1) = (c2(3:n)-2*c2(2:n-1)+c2(1:n-2))/h^2;
dc2dt = -v*dc2dz+D(2)*d2c2dz2-k(1)*qm(1)*(1/epl-1)*LT(:,2);
dc2dt(1) = c2(1)-cFeed(2);
dc2dt(n) = c2(n)-c2(n-1);
% Saturation equations
dq1dt = k(1)*qm(1)*LT(:,1);
dq2dt = k(2)*qm(2)*LT(:,2);
% Concatenate vectors of time derivatives
dydt = [dc1dt;dq1dt;dc2dt;dq2dt];
end
##### 7 commentairesAfficher 5 commentaires plus anciensMasquer 5 commentaires plus anciens
Davide Masiello le 27 Juil 2023
Modifié(e) : Davide Masiello le 27 Juil 2023
@Samridh Sharma I have edited my answer, now it solves the system for 2 gases.
I think the problem in your case was how you were concatenating the c and q vectors, and possibly other modifications needed in the rest of the code to accommodate the different number of species.
I also assume that the two different compounds will have different values of D, k, and qm, but not knowing them I just used the same values for both.
Moreover, they would have different cFeed I guess. I just used 0.1 for the second species for all the simulations.
You can easily change these values in the code.
Please let us know if this works and don't forget to accept the answer if you think it helped!
Cheers.
Samridh Sharma le 27 Juil 2023
Thank you @Davide Masiello, this works.
Yes, you are right that various parameters will be different for different species of gases. I will edit those changes with correct values.
Thank you, this was very helpful.

Connectez-vous pour commenter.

### Catégories

En savoir plus sur Mathematics and Optimization dans Help Center et File Exchange

R2023a

### Community Treasure Hunt

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

Start Hunting!

Translated by