清空筛选条件
清空筛选条件

fixed-bed adsorption column-modeling solving PDE

25 次查看（过去 30 天）
nashwa tarek
nashwa tarek 2021-4-8
评论： Torsten 2022-12-30
hi all i am trying obtainbreakthrough curve from that model (c/co vs time) i want to study the effect of initial concentration, flow rate and column length. my problem now that by changing the initial concentration ther is no change in the curve, i am not good at matlab i just want to use it to finish my thesis. can anyone help me please to figure out the problem
function main
%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 = 20; % Feed concentration
L = 0.02; % Column length
t0 = 0; % Initial Time
%tf = 20; % Final time
tf = 20; % Final time
dt = 0.5; % 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);
plot(T, Y(:,n)/cFeed)
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) = ((z(i)-z(i-1))/(z(i+1)-z(i))*(c(i+1)-c(i))+(z(i+1)-z(i))/(z(i)-z(i-1))*(c(i)-c(i-1)))/(z(i+1)-z(i-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

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

回答（1 个）

Sulaymon Eshkabilov
Sulaymon Eshkabilov 2022-12-30
Ran in:
Here is an edited and corrected version of your code (smaller time step 'cause it is a stiff DE) that shows that the concentration change makes an impact on the simulation results:
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), '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 个评论
显示 1更早的评论
Torsten
Torsten 2022-12-30
In cartesian coordinates:
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))
In cylindrical or spherical coordinates, it will be different, of course.

请先登录，再进行评论。

类别

在 Help Center 和 File Exchange 中查找有关 Plasma Physics 的更多信息

标签

Community Treasure Hunt

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

Start Hunting!

Translated by