Plot-Parameter study- Adsorption

2 次查看（过去 30 天）
Una
Una 2023-1-12
编辑： Torsten 2023-1-18
Ran in:
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

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

回答（1 个）

Torsten
Torsten 2023-1-12
编辑：Torsten 2023-1-12
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.
  5 个评论
显示 3更早的评论
Una
Una 2023-1-18
Ran in:
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
Torsten 2023-1-18
编辑：Torsten 2023-1-18
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 );

请先登录，再进行评论。

类别

在 Help Center 和 File Exchange 中查找有关 Solver Outputs and Iterative Display 的更多信息

标签

产品


版本

R2021b

Community Treasure Hunt

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

Start Hunting!

Translated by