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

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

[Davide Masiello]

Hi @Samridh Sharma, I;ve had experience with similar problems before, it would be helpful if you could post the governing equations.

[Samridh Sharma]

@Davide Masiello thanks for your remark, I have updated my question with more details.

[Davide Masiello]

Much better. Would you mind adding you BCs too please?

Answer 1

[Davide Masiello]

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)
    sol(k)  = adsorptionSolver(t,z,cF(:,k));
    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
function out = adsorptionSolver(t,z,cFeed)
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
out = ode15s(@(t,y) adsorptionModel(t,y,h,n,cFeed),t,y0,opts);
 
end

%% Adsorption model
function dydt = adsorptionModel(~,y,h,n,cFeed)

% 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

[Samridh Sharma]

Thanks for your promptness @Davide Masiello, this solves the problem when adsorption bed has only one feed gas. If the feed comprises of 2 gas in unequal ratio, then the Langmuir term has the feed concentration of the two gases couple, as shown in equation (2).

I wish to know how to solve equation (3) when two gases with their different langmuir parameter are getting adsorbed in single bed filter in a pressure swing adsorption.

[Torsten]

I wish to know how to solve equation (3) when two gases with their different langmuir parameter are getting adsorbed in single bed filter in a pressure swing adsorption.

Now the solution vector is

y0 = [c0 ; q0];

In case of two components, it will be

y0 = [c10 ; q10 ; c20 ; q20]

Now you rewrite the solution in local variables as

c = y(1:n);
q = y(n+1:2*n);

In case of two components, it will be

c1 = y(1:n);
q1 = y(n+1:2*n);
c2 = y(2*n+1:3*n);
q2 = y(3*n+1:4*n);

Modify everything else accordingly.

[Davide Masiello]

I agree with @Torsten, you can add as many components as you want in that way.

With some more work you could even generalize the code to accept a variable number of gas components which you pass to the solver every time you solve the system.

[Davide Masiello]

And please, @Samridh Sharma do let us know if it worked!

[Samridh Sharma]

@Davide Masiello and @Torsten, thanks for patiently helping me with my problem.

For the function adsorptionModel, the langmuir term should be:

Here, when I have added the second gas with concentration c2, then for determining the concentration at every point of defined length, I need to solve with the c1 and c2 values together. But the array formed on RHS for solving LT is not matching and dcdt also gives error. Please help me with these corrections.

function dydt = adsorptionModel(~,y,h,n,cFeed)
% Constants
D   = 1.85e-5;                                                              % Axial Dispersion coefficient
v   = 3.38e-3 ;                                                             % Superficial velocity
epl = 0.47;                                                                 % Void fraction
k   = 1.9e-6;                                                               % Mass Transfer Coefficient
b   = [0.05 0.1];                                                           % Langmuir parameter
qm  = 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);
c = [c1 ; c2];
q = [q1 , q2];
% Langmuir term
LT = (c*b)./(1+c*b)-q;  
% LT = [b(1).*c1./(1+b(1).*c1+b(2).*c2)-q1; b(2).*c2./(1+b(1).*c1+b(2).*c2)-q2];
% LT = b*c./(1+b*c)-q;
% Concentration equation
dcdz            = zeros(n,1);
dcdz(2:n-1)     = (c(3:n)-c(2:n-1))/h;
d2cdz2          = zeros(n,1);
d2cdz2(2:n-1)   = (c(3:n)-2*c(2:n-1)+c(1:n-2))/h^2;
dcdt            = -v*dcdz+D*d2cdz2-k*qm*(1/epl-1)*LT;
dcdt(1)         = c(1)-cFeed;
dcdt(n)         = c(n)-c(n-1);
% Saturation equation
dqdt = k*qm*LT;
% Concatenate vectors of time derivatives
dydt = [dcdt;dqdt];
end

[Davide Masiello]

@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]

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.