fixed-bed adsorption column-modeling solving PDE
11 Ansichten (letzte 30 Tage)
Ältere Kommentare anzeigen
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 Kommentar
Antworten (1)
Sulaymon Eshkabilov
am 30 Dez. 2022
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 Kommentare
Torsten
am 30 Dez. 2022
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.
Siehe auch
Kategorien
Mehr zu Plasma Physics finden Sie in Help Center und File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!