Plot-Parameter study- Adsorption

Dear all,

I am trying to do a parameter study (velocity and length) on an adsorbent. I had the following idea on how to do it but it doesnt work and I cant figure out how to solve the error.

I am relatively new in matlab so for any suggestion I am grateful.

velocity = [0.01 0.1 0.5 1 1.5 2 2.5 3 3.5];   % Diffferent velocity values
for ii=1:numel(velocity)
    v=velocity(ii);
    [T, Y] = Main(v);
    plot(T,Y(:,n), 'linewidth', 1),
    hold all
end

Unrecognized function or variable 'v'.

Error in solution>MyFun (line 44)
Dl = Dm +(v*r_1*2)^2/(192*Dm); % axial dipsersion coefficient in m2/sec

Error in solution>@(t,Y)MyFun(t,Y,z,n) (line 27)
[T, Y] = ode15s(@(t,Y) MyFun(t,Y,z,n),t,y0);

Error in odearguments (line 92)
f0 = ode(t0,y0,args{:}); % ODE15I sets args{1} to yp0.

Error in ode15s (line 153)
odearguments(odeIsFuncHandle, odeTreatAsMFile, solver_name, ode, tspan, y0, options, varargin);

Error in solution>Main (line 27)
[T, Y] = ode15s(@(t,Y) MyFun(t,Y,z,n),t,y0);

function [T, Y] = Main(~)

cFeed = 0.016;            % Feed concentration in mol/m3
L = 0.3;                  % Column length in m
t0 = 0;                   % Initial Time in s            
tf = 8000;                % Final time in s
dt = 0.5;                 % Time step
z = 0:0.005:L;            % Mesh generation
t = t0:dt:tf;             % Time vector
n = numel(z);             % Size of mesh grid

% Initial Conditions / Vector Creation
c0 = zeros(n,1);          % c = 0 at t = 0 for all z
c0(1) = cFeed;            % c = cFeed at z= 0 for all t
q0 = zeros(n,1);          %  q = 0 at t = 0 for all z
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)

% Design parameters of the adsorbent bed
epsilon         = 0.87;                                % bed porosity
%v               = 3;                                  % superficial velocity in m/s
rho             = 500;                                 % particle density in kg/m3
r_1             = 0.000525;                            % Channel inner radius in m                         
L               = 0.3;                                 % bed length in m                                 

% General parameters                
De          = 1e-11;                                   % effective diffusivity 
Dm          = 1.381e-5;                                % molecular diffusivity of co2 in air in m2/s
Dl          = Dm +(v*r_1*2)^2/(192*Dm);                % axial dipsersion coefficient in m2/sec
k           = 0.005;                                   % mass transfer coefficient in 1/sec
R           = 8.314;                                   % gas constant in J/mol-K
T           = 298;                                     % gas temp in K

% Parameters of the isotherm 
qsat1       = 0.16;                                    % mol/kg
b1          = 4.89e-2;                                 % 1/pascal
qsat2       = 3.5;                                     % mol/kg
b2          = 22e-2;                                   % 1/pascal
b3          = 0.00062e-2;                              % mol/kg-pascal

% 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 
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) = ((c(i+1)-c(i))/(z(i+1)-z(i))-(c(i)-c(i-1))/(z(i)-z(i-1)))/(zhalf(i)-zhalf(i-1));
end

% 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));

% time derivatives at z=0 and q* 
 if R*T*c(1) < 16
    qstar=(qsat1*b1*R*T*c(1))/(1+b1*R*T*c(1));
else 
    qstar=(qsat2*b2*(R*T*c(1)-16))/(1+b2*(R*T*c(1)-16)) + b3*R*T*c(1);
end
DqDt(1) = k*(qstar - q(1) );
DcDt(1) = 0.0;

% Set time derivatives in remaining points
for i=2:n
    % q* for RTc < or > 16
 if R*T*c(i) <16
    qstar=(qsat1*b1*R*T*c(i))/(1+b1*R*T*c(i));
else 
    qstar=(qsat2*b2*(R*T*c(i)-16))/(1+b2*(R*T*c(i)-16)) + b3*R*T*c(i);
 end
    % Equation: dq/dt = k(q*-q)
    DqDt(i) = k*(qstar - q(i) );
      % Equation: dc/dt = Dl * d2c/dz2 - v*dc/dz - ((1-e)/(e))*roh*dq/dt
      DcDt(i) = Dl*D2cDz2(i) - v*DcDz(i) - ((1-epsilon)/(epsilon))*rho*DqDt(i);
  end
% Concatenate vector of time derivatives
DyDt = [DcDt;DqDt];
end

Answer 1

[Torsten]

I don't know whether you and Sona Ghulami are the same person, but i suggest you start with the last code under

https://de.mathworks.com/matlabcentral/answers/1730075-how-to-solve-coupled-partial-differential-equations-with-method-of-lines?s_tid=srchtitle

because several mistakes in the code have been eliminated.

Especially change

DcDz(n) = 0;

to

DcDz(n) = (c(n)-c(n-1))/(z(n)-z(n-1));

Concerning your question:

Change

function [T, Y] = Main(~)

to

function [T, Y] = Main(v)

and

[T, Y] = ode15s(@(t,Y) MyFun(t,Y,z,n),t,y0);

to

[T, Y] = ode15s(@(t,Y) MyFun(t,Y,z,n,v),t,y0);

and

function DyDt=MyFun(~, y, z, n)

to

function DyDt=MyFun(~, y, z, n, v)

And I explicitly suggest that you plot qstar over c before you continue.

[Una]

Thank you. I have followed that post. Why does it need to be changed? I, of course, followed your suggestion.The plot doesnt seem right but I cant figure out what the problem is yet, so I am going with what I have for now (I am attaching the isotherm equations and parmaters).

R           = 8.314;              % gas constant in J/mol-K
T           = 298;                % gas temp in K

% Parameters of the isotherm 
qsat1       = 0.16;               % mol/kg (changed from 0.16 mmol/g)
b1          = 4.89e-2;            % 1/pascal (changed from 4.89 1/mbar)
qsat2       = 3.5;                % mol/kg  (changed from 3.50 mmol/g)
b2          = 22.00e-2;           % 1/pascal (changed from 22.00 1/mbar)
b3          = 0.00062e-2;         % mol/kg-pascal (changed from 0.00062 mmol/g-mbar)
c=[0:0.01:1];
 if R*T*c < 16
    qstar=(qsat1*b1*R*T*c)/(1+b1*R*T*c);
else 
    qstar=(qsat2*b2*(R*T*c-16))/(1+b2*(R*T*c-16)) + b3*R*T*c;
 end
 plot (c,qstar)

[Torsten]

Here is the plot of qstar against c. I doubt that the equation for equilibrium loading qstar is correct.

R           = 8.314;              % gas constant in J/mol-K
T           = 298;                % gas temp in K
% Parameters of the isotherm 
qsat1       = 0.16;               % mol/kg (changed from 0.16 mmol/g)
b1          = 4.89e-2;            % 1/pascal (changed from 4.89 1/mbar)
qsat2       = 3.5;                % mol/kg  (changed from 3.50 mmol/g)
b2          = 22.00e-2;           % 1/pascal (changed from 22.00 1/mbar)
b3          = 0.00062e-2;         % mol/kg-pascal (changed from 0.00062 mmol/g-mbar)
c=[0:0.1:1];
for i=1:numel(c)
    if R*T*c(i) < 16
        qstar(i)=(qsat1*b1*R*T*c(i))/(1+b1*R*T*c(i));
    else 
        qstar(i)=(qsat2*b2*(R*T*c(i)-16))/(1+b2*(R*T*c(i)-16)) + b3*R*T*c(i);
    end
end
plot (c,qstar)

[Torsten]

Could you change the parameters such that one can see the overshoot c/d_feed = 1.2 ?

cFeed = 0.016;            % Feed concentration in mol/m3
L = 0.3;                  % Column length in m
t0 = 0;                   % Initial Time in s            
tf = 8000;                % Final time in s
dt = 0.5;                 % Time step
z = 0:0.005:L;            % Mesh generation
t = t0:dt:tf;             % Time vector
n = numel(z);             % Size of mesh grid

% Initial Conditions / Vector Creation
c0 = zeros(n,1);          % c = 0 at t = 0 for all z
c0(1) = cFeed;            % c = cFeed at z= 0 for all t
q0 = zeros(n,1);          %  q = 0 at t = 0 for all z
y0 = [c0 ; q0];           % Appends conditions together

%ODE15S Solver
[T, Y] = ode15s(@(t,Y) MyFun(t,Y,z,n),t,y0);
figure(1)
plot(z,[Y(100,1:n);Y(250,1:n);Y(500,1:n);Y(1000,1:n);Y(2000,1:n);Y(4000,1:n);Y(8000,1:n)], 'linewidth', 1)

figure(2)
plot(z,[Y(100,n+1:2*n);Y(250,n+1:2*n);Y(500,n+1:2*n);Y(1000,n+1:2*n);Y(2000,n+1:2*n);Y(4000,n+1:2*n);Y(8000,n+1:2*n)], 'linewidth', 1)

function DyDt=MyFun(~, y, z, n)

% Design parameters of the adsorbent bed
epsilon         = 0.87;                                % bed porosity
v               = 1;                                  % superficial velocity in m/s
rho             = 500;                                 % particle density in kg/m3
r_1             = 0.000525;                            % Channel inner radius in m                         
L               = 0.3;                                 % bed length in m                                 

% General parameters                
De          = 1e-11;                                   % effective diffusivity 
Dm          = 1.381e-5;                                % molecular diffusivity of co2 in air in m2/s
Dl          = Dm +(v*r_1*2)^2/(192*Dm);                % axial dipsersion coefficient in m2/sec
k           = 0.005;                                   % mass transfer coefficient in 1/sec
R           = 8.314;                                   % gas constant in J/mol-K
T           = 298;                                     % gas temp in K

% Parameters of the isotherm 
qsat1       = 0.16;                                    % mol/kg
b1          = 4.89e-2;                                 % 1/pascal
qsat2       = 3.5;                                     % mol/kg
b2          = 22e-2;                                   % 1/pascal
b3          = 0.00062e-2;                              % mol/kg-pascal

% 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 
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) = ((c(i+1)-c(i))/(z(i+1)-z(i))-(c(i)-c(i-1))/(z(i)-z(i-1)))/(zhalf(i)-zhalf(i-1));
end

% DcDz and D2cDz2 at z=L for boundary condition dc/dz = 0
DcDz(n) = (c(n)-c(n-1))/(z(n)-z(n-1)); 
D2cDz2(n) = -1.0/(z(n)-zhalf(n-1))*DcDz(n);

% time derivatives at z=0 and q* 
 if R*T*c(1) < 16
    qstar=(qsat1*b1*R*T*c(1))/(1+b1*R*T*c(1));
else 
    qstar=(qsat2*b2*(R*T*c(1)-16))/(1+b2*(R*T*c(1)-16)) + b3*R*T*c(1);
end
DqDt(1) = k*(qstar - q(1) );
DcDt(1) = 0.0;

% Set time derivatives in remaining points
for i=2:n
    % q* for RTc < or > 16
 if R*T*c(i) <16
    qstar=(qsat1*b1*R*T*c(i))/(1+b1*R*T*c(i));
else 
    qstar=(qsat2*b2*(R*T*c(i)-16))/(1+b2*(R*T*c(i)-16)) + b3*R*T*c(i);
 end
    % Equation: dq/dt = k(q*-q)
    DqDt(i) = k*(qstar - q(i) );
      % Equation: dc/dt = Dl * d2c/dz2 - v*dc/dz - ((1-e)/(e))*roh*dq/dt
      DcDt(i) = Dl*D2cDz2(i) - v*DcDz(i) - ((1-epsilon)/(epsilon))*rho*DqDt(i);
  end
% Concatenate vector of time derivatives
DyDt = [DcDt;DqDt];
end

[Una]

Hallo Torston, I am sorry to bother you again. I have made a few changes according to your suggestions but now i get the following error. I tried using the "which function" but I couldnt find anything odd.

velocity = [0.01  0.5 1.5 2.5 3];   % Diffferent velocity values
for ii=1:numel(velocity)
  v=velocity(ii);
  [T, Y] = Mainmmen(v);
  plot(T,Y(:,n), 'linewidth', 1),
   hold all
end

Error using solution>Mainmmen
Too many output arguments.

function Mainmmen(v)
cFeed= 0.016;         % Feed concentration in mol/m3
L = 0.3;              % Column length in m
t0 = 0;               % Initial Time in s            
tf= 15000;            % Final time in s
dt= 0.05;                % Time step
z = [0:0.005:L].';    % Mesh generation
t = [t0:dt:tf];       % Time vector
n = numel(z);         % Size of mesh grid


%Initial Conditions and boundary condition
c0 = zeros(n,1);        % c = 0 at t =0 for all z
c0(1)= cFeed;           % c = cFeed at z =0 and t>0
q0 = zeros(n,1);        % t = 0, q = 0 for all z
y0 = [c0 ; q0];         

%ODE15S Solver
[T, Y] = ode15s(@(t,y) MyFun(t,y,z,n,v),t,y0);

end 



function DyDt=MyFun(t, y, z, n,v)

% Design parameters of the adsorbent bed
epsilon         = 0.80;                                % bed porosity
%vs              = 3;                                  % superficial velocity in m/s
%v               = vs/epsilon;
rho             = 500;                                 % particle density in kg/m3
r_1             = 0.000525;                            % Channel inner radius in m
r_2             = 0.000585;                                % Channel radius in m 
L               = 0.3;                                 % bed length in m                                 

% General parameters                
De          = 1e-11;                                       % effective diffusivity 
Dm          = 1.6e-5;                                      % molecular diffusivity of co2 in air in m2/s
Dl          = Dm +(v*r_1*2)^2/(192*Dm);                    % axial dipsersion coefficient in m2/sec
rohg        = 1.18;                                        % density of air in kg/m3
viskog      = 18.37e-6;                                    % dynamic viscosity of air in Ns/m2
Sc          = viskog/Dm;                                   % schmidt number
Re          = (rohg*v*2*r_1)/viskog;
sh          = 3.6*(1+0.139*(r_1*2/L)*Re*Sc).^0.81;
kf          = (sh*Dm)/(2*r_1);
k           = 1/((r_1/(2*kf))+(r_2^2-r_1^2)/(8*De));
R           = 8.314;                                    % gas constant in J/mol-k
T           = 298;                                      % gas temp in K

% Parameters of the isotherm 
qsat1       = 0.16;                                    % mol/kg
b1          = 4.89e-2;                                 % 1/pascal
qsat2       = 3.5;                                     % mol/kg
b2          = 22e-2;                                   % 1/pascal
b3          = 0.00062e-2;                              % mol/kg-pascal
% Variables being allocated zero vectors
c = zeros(n,1);
q = zeros(n,1);
DcDt = zeros(n,1);
DqDt = zeros(n,1);
zhalf = zeros(n-1,1);
DcDz = zeros(n,1);
DyDt = zeros(2*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;

% Calculate spatial derivatives
DcDz(2:n) = (c(2:n)-c(1:n-1))./(z(2:n)-z(1:n-1));
D2cDz2(2:n-1) = ((c(3:n)-c(2:n-1))./(z(3:n)-z(2:n-1))-(c(2:n-1)-c(1:n-2))./(z(2:n-1)-z(1:n-2)))./(zhalf(2:n-1)-zhalf(1:n-2));

% Calculate D2cDz2 at z=L for boundary condition dc/dz = 0
D2cDz2(n) = -1.0/(z(n)-zhalf(n-1))*DcDz(n);

% Set time derivatives

for i=1:numel(c)
    if R*T*c(i) < 16
        qstar(i)=(qsat1*b1*R*T*c(i))/(1+b1*R*T*c(i));
    else 
        qstar(i)=(qsat2*b2*(R*T*c(i)-16))/(1+b2*(R*T*c(i)-16)) + b3*R*T*c(i);
    end
end           
DqDt = k*(qstar(i) - q );
DcDt(1) = 0.0;
DcDt(2:n) = Dl*D2cDz2(2:n) - v*DcDz(2:n) - ((1-epsilon)/(epsilon))*rho*DqDt(2:n);

% Concatenate vector of time derivatives
DyDt = [DcDt;DqDt];
end

[Torsten]

Read again my suggested changes in the answer above.

If you want to have T and Y as output from "Mainmmen", you must define them as output arguments within the function:

function [T,Y] = mainmmen(v)

And change

DqDt = k*(qstar(i) - q );

to

DqDt = k*(qstar - q );