fixed bed adsorption column model-solving PDE-freundlish isotherm
34 Ansichten (letzte 30 Tage)
Ältere Kommentare anzeigen
function main
%System Set-up %
%Define Variables
%D = 3*10^-8; % Axial Dispersion coefficient
%v = 1*10^-3; % Superficial velocity
%epsilon = 0.4; % Voidage fraction
%k = 3*10^-5; % Mass Transfer Coefficient
%Kf=2.5*10^-5; %freundlish parameter
%nf= 1.45; %freundlish constant
cFeed = 10; % Feed concentration
L = 1; % Column length
t0 = 0; % Initial Time
%tf = 1000; % Final time
tf = 2000; % 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);
plot(T,Y(:,n)/cFeed)
end
function DyDt=MyFun(~, y, z, n)
% Defining Constants
D = 3*10^-8; % Axial Dispersion coefficient
v = 1*10^-3; % Superficial velocity
epsilon = 0.4; % Voidage fraction
k = 3*10^-5; % Mass Transfer Coefficient
Kf=2.5*10^-5; %freundlish parameter
nf= 1.45; %freundlish constant
% 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*(kf*(c(1)^(1/nf))-q(1));
DcDt(1) = 0.0;
% Set time derivatives in remaining points
for i=2:n
%Equation: dq/dt = K(q*-q) where q*=kf*(c^(1/nf))
DqDt(i) = k*(kf*(c(i)^(1/nf))-q(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))*DqDt(i);
end
% Concatenate vector of time derivatives
DyDt = [DcDt;DqDt];
end
This is the code which I am using. I don't have much knowledge about Matlab but I need to finish the code for my thesis. My question is, I need to model a fixed bed to obtain the breakthrough curves using freundlish isotherm, so, I would like to know how to obtain these curves (C/C0 vs t) from this code. i cant run this code i dont know where is the problem. I would appreciate any help. Thank you very much
3 Kommentare
Antworten (4)
Alan Stevens
am 9 Jul. 2020
The following works. I can't say if the result is sensible or not - you will need to judge that. I'm not convinced that the coded equations match those in your hand-written attachment!
cFeed = 10; % Feed concentration
L = 1; % Column length
t0 = 0; % Initial Time
tf = 2000; % Final time
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);
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);
plot(T,Y(:,n)/cFeed)
function DyDt=MyFun(~, y, z, n)
% Defining Constants
D = 3*10^-8; % Axial Dispersion coefficient
v = 1*10^-3; % Superficial velocity
epsilon = 0.4; % Voidage fraction
k = 3*10^-5; % Mass Transfer Coefficient
Kf=2.5*10^-5; %freundlish parameter
nf= 1.45; %freundlish constant
% 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*(Kf*(c(1)^(1/nf))-q(1));
DcDt(1) = 0.0;
% Set time derivatives in remaining points
for i=2:n
%Equation: dq/dt = K(q*-q) where q*=kf*(c^(1/nf))
DqDt(i) = k*(Kf*(c(i)^(1/nf))-q(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))*DqDt(i);
end
% Concatenate vector of time derivatives
% Very small imaginary components need to be removed!
DyDt = [real(DcDt);real(DqDt)];
end
22 Kommentare
Rohan Singh
am 27 Jan. 2021
Sorry to bother you again @Alan Stevens. I am still learning this stuff. I tried to implement your suggested way but for some reason, i am getting a strange curve. It would really help me if you can guide me further what to do. Thanks in advance.(I have attached the equations and boundary conditions for your reference)
Cfeed = 10; % Inlet concentration
tf = 86400; % End time
nz = 200; % Number of gateways
np = nz + 1; % Number of gateposts
% Initial Concentrations
C = zeros(np,1); C(1) = Cfeed;
Cs = zeros(np,1);
% Prepare data to pass to ode solver
c0 = [C; Cs];
dt = tf/10000;
tspan = 0:dt:tf;
% Call ode solver
[t, c] = ode15s(@rates, tspan, c0);
% Plot results
plot(t, c(:,np)/Cfeed),grid
xlabel('time'), ylabel('Cb/Cfeed')
title('Figure 1: Exit Cb/Cfeed vs time')
% figure
% plot(t, c(:,2*np)),grid
% xlabel('time'), ylabel('Cs')
% title('Figure 2: Exit Cs vs time')
% figure
% plot(t, c(:,np)-c(:,2*np)), grid
% xlabel('time'), ylabel('Cb - Cs')
% title('Figure 3: Exit Cb-Cs vs time')
%
% k = 2.5*10^-5; nf = 1.45;
% q = k*c(:,2*np).^(1/nf); % Freundlcih equation
% figure
% plot(t,q), grid
% xlabel('time'), ylabel('q')
% title('Figure 4: Exit q vs time')
% Function to calculate rate of change of concentrations with time
function dcdt = rates(~, c)
Dz = 3*10^-8; % Diffusion coefficient
u = 2*10^-5; % Velocity
e = 0.4;
Kf = 3*10^-5;
ap = 2*10^-3; % Particle diameter
ep = 0.53;
L = 0.1; % Length
nz = 200; % Number of gateways
np = nz + 1; % Number of gateposts
dz = L/nz;
dr = ap/nz;
rho = 400;
Dp = 1.2*10^-10;
b=0.03;
qm=118;
C = c(1:np);
Cs = c(np+1:end);
Cs(1) = 0;
dCdt = zeros(np,1);
dCsdt = zeros(np,1);
dCsdr = zeros(np,1);
% rhalf = zeros(np-1,1);
% DCsDr = zeros(np,1);
D2CsDr2 = zeros(np,1);
%
% rhalf(1:np-1)=(z(1:np-1)+z(2:np))/2;
for i = 2:np-1
dCsdr(i) = (2*dr)*(Cs(i+1)-Cs(i-1));
D2CsDr2(i) = (1/dr^2)*(Cs(i+1)-2*Cs(i)+Cs(i-1));
% dCsdt(i) = 3*Kf/(ep*ap)*(C(i)-Cs(i)); % Needed to change the sign from -ve to +ve here.
end % i loop
for i = 2:np-1
dCdt(i) = (Dz/dz^2)*(C(i+1)-2*C(i)+C(i-1)) - u/(2*dz)*(C(i+1)-C(i-1)) - (1-e)*3*Kf/(e*ap)*(C(i)-Cs(i));
% dCsdt(i) = 3*Kf/(ep*ap)*(C(i)-Cs(i)); % Needed to change the sign from -ve to +ve here.
dCsdt(i) = 1/(1+rho*((1-ep)/ep)*(qm*b/((1+b*Cs(i))^2)))*Dp*(D2CsDr2(i)+((2/ap)*dCsdr(i)));
end % i loop
% For simplicity, set exit rates to those of the immediately previous nodes
dCdt(np) = dCdt(np-1);
dCsdt(np) = dCsdt(np-1);
% Combine rates for the function to return to the calling routine
dcdt = [dCdt; dCsdt];
end % of rates function
Anthony Ozzello
am 28 Dez. 2021
I'm sorry to follow up on an old thread. I'm a relatively new user to Matlab and I'm trying to model a fixed bed reactor like the folks here, but unlike them, I'd like to be able to visualize the values of C and q throughout the depth of the bed as a function of t, so I need the function to return C and q as a matrix of both t and z, I believe. Is there a simple way to do this? I originally thought I would have to use the pdesolver, but this is an interesting way to handle this. Please let me know if there is a simple way to look at the results as both a function of t and z so that I can do cumtraps and surface plots throughout the depth of the fixed bed. Thanks, Tony
9 Kommentare
Anthony Ozzello
am 2 Jan. 2022
Hey again,
I'm working on fitting some of the parameters in the ODE using the problem-based approach outline in the MATLAB documentation here: https://www.mathworks.com/help/optim/ug/fit-ode-problem-based-least-squares.html. I've got everything working properly except for plotting out the results of the fitted data points against the original data. My data exists in time space at non-evenly spaced intervals from 5 to 380 minutes and to get a smooth function, I'm integrating every 0.1 minute. I modified the RtoODE function to only return the solution points that match the original dataset but using this code, where exp_t contains the timestamps of the original data and tspan(2) is the time interval in the integration. There can be multiple components in the analysis all concatented together into y.
function solpts = RtoODE(r,z,n,nComponents,tspan,y0, exp_t)
[T, Y] = ode15s(@(t,y) Myfun(t,y,z,n,nComponents,r),tspan,y0);
for i = 1:length(exp_t)
time = exp_t(i);
j= time/tspan(2)+1;
for k = 1: nComponents
solpts(i,k)=Y(j,k*n);
end
end
end
But when I got to plot the data that comes back, I get an error, where exp_t and exp_y are my raw data, there are 2 components for y, but I'm only trying to plot the first one right now.
myfcn = fcn2optimexpr(@RtoODE,r,z,n,nComponents,tspan,y0, exp_t);
figure
plot (exp_t,exp_y(:,1),'g')
hold on
plot(exp_t,myfcn(:,1),'b')
hold off
Error using plot
Data must be numeric, datetime, duration or an array convertible to double.
It should work based on the example in the documentation, but it doesn't and I can't find anything in any of the previous question and answers or anything elsewhere online that gives me some ideas on how to plot optimization variables if this isn't it.
Any clues?
Thanks,
Tony
Anthony Ozzello
am 10 Jan. 2022
Hi yet again,
I followed Torsten's response, which was very helpful to setting the time span, which enabled me to compare the results of the calculated integration to my experimental data points very well, but I encountered a problem. It seems that my problem requires significant calculation at a lot of points near the origin and then can space out away from it. When I specify the tspan to ensure I capture my experimental points, I have to use more of a linspace to ensure I have the right timestamps and, even if I go with pretty small time intervals, it's nowhere near tight enough near the origin and I'm getting negative concentrations. So it's clear that I have to let MATLAB chose the timespan properly, I think, by just specifying the end points. But then my question is how do I do some kind of fuzzy matching between the timestamps of my experimental data and the timestamps that come back from ode15s for calculating the least squares fit? Can someone point me in a direction or some examples? Thanks, Tony
32 Kommentare
Federico Urbinati
am 19 Okt. 2022
Hi Anthony, have you completed the code related to the fixed bed adsorption of multicomponents?
i'm really stuck, you would do me a favor if i posted it here.
thanks in advance,
Federico
Puteri Ellyana Mohmad Latfi
am 16 Apr. 2022
I am not really good in Matlab but I have to model a fixed bed adsorption to obtain the breakthrough curve of the adsorption. I tried to incorporate my model in Alan's code for Nashwa but I am not really sure if the initial and boundary condition in the code are correct though. Beside that, I got this error when i ran this code.
Warning: Failure at t=2.006182e+02. Unable to meet integration tolerances without reducing the step size below the
smallest value allowed (4.547474e-13) at time t.
> In ode15s (line 655)
In fyp (line 19)
I would really appreciate any of your help. Thank you very much. I have also attached the model's equations in the file.
Cfeed = 0.0224; % Inlet concentration
tf = 2000; % End time
nz = 200; % Number of gateways
np = nz + 1; % Number of gateposts
% Initial Concentrations
Cb = zeros(np,1); Cb(1) = Cfeed;
Cp = zeros(np,1);
% Prepare data to pass to ode solver
c0 = [Cb; Cp];
dt = tf/100;
tspan = 0:dt:tf;
% Call ode solver
[t, c] = ode15s(@rates, tspan, c0);
% Plot results
plot(t, c(:,np)/Cfeed),grid
xlabel('time'), ylabel('Cb/Cfeed')
title('Figure 1: Exit Cb/Cfeed vs time')
figure
plot(t, c(:,2*np)),grid
xlabel('time'), ylabel('Cp')
title('Figure 2: Exit Cp vs time')
figure
plot(t, c(:,np)-c(:,2*np)), grid
xlabel('time'), ylabel('Cb - Cp')
title('Figure 3: Exit Cb-Cp vs time')
k = 598;
q = k*c(:,2*np); % Freundlcih equation
figure
plot(t,q), grid
xlabel('time'), ylabel('q')
title('Figure 4: Exit q vs time')
% Function to calculate rate of change of concentrations with time
function dcdt = rates(~, c)
a=0.00000000261;
b=0.2102;
d=1705.4719;
e=-38.8576;
L=0.08;
nz = 200; % Number of gateways
np = nz + 1; % Number of gateposts
dz = L/nz;
Cb = c(1:np);
Cp = c(np+1:end);
dCbdt = zeros(np,1);
dCpdt = zeros(np,1);
for i = 2:np-1
dCbdt(i) = (a/dz^2)*(Cb(i+1)-2*Cb(i)+Cb(i-1)) - b/(2*dz)*(Cb(i+1)-Cb(i-1)) - d*(Cb(i)-Cp(i));
dCpdt(i) = e*(Cb(i)-Cp(i)); % Needed to change the sign from -ve to +ve here.
end % i loop
% For simplicity, set exit rates to those of the immediately previous nodes
dCbdt(np) = dCbdt(np-1);
dCpdt(np) = dCpdt(np-1);
% Combine rates for the function to return to the calling routine
dcdt = [dCbdt; dCpdt];
end % of rates function
5 Kommentare
Puteri Ellyana Mohmad Latfi
am 13 Mai 2022
Bearbeitet: Puteri Ellyana Mohmad Latfi
am 13 Mai 2022
Hi @Alan Stevens. I tried using the coding for different set of experimental values, however, i didnt get the same graph as the experimental breakthrough graph, which is shown in the image attached (10000ppm).
Cfeed = 0.0227; % Inlet concentration
tf = 100; % End time
nz = 200; % Number of gateways
np = nz + 1; % Number of gateposts
% Initial Concentrations
Cb = zeros(np,1); Cb(1) = Cfeed;
Cp = zeros(np,1);
% Prepare data to pass to ode solver
c0 = [Cb; Cp];
dt = tf/100;
tspan = 0:dt:tf;
% Call ode solver
[t, c] = ode15s(@rates, tspan, c0);
% Plot results
plot(t, c(:,np)/Cfeed),grid
xlabel('time'), ylabel('Cb/Cfeed')
title('Figure 1: Exit Cb/Cfeed vs time')
% figure
% plot(t, c(:,2*np)),grid
% xlabel('time'), ylabel('Cp')
% title('Figure 2: Exit Cp vs time')
%
% figure
% plot(t, c(:,np)-c(:,2*np)), grid
% xlabel('time'), ylabel('Cb - Cp')
% title('Figure 3: Exit Cb-Cp vs time')
%
% k = 56.5;
% q = k*c(:,2*np); % Freundlcih equation
% figure
% plot(t,q), grid
% xlabel('time'), ylabel('q')
% title('Figure 4: Exit q vs time')
% Function to calculate rate of change of concentrations with time
function dcdt = rates(~, c)
a=0.00000001278;
b=0.160556;
d=244.942875;
e=-1.3964;
L=0.065;
nz = 200; % Number of gateways
np = nz + 1; % Number of gateposts
dz = L/nz;
Cb = c(1:np);
Cp = c(np+1:end);
dCbdt = zeros(np,1);
dCpdt = zeros(np,1);
for i = 2:np-1
dCbdt(i) = (a/dz^2)*(Cb(i+1)-2*Cb(i)+Cb(i-1)) - b/(dz)*(Cb(i+1)-Cb(i-1)) - d*(Cb(i)-Cp(i));
dCpdt(i) = -e*(Cb(i)-Cp(i));
end % i loop
% For simplicity, set exit rates to those of the immediately previous nodes
dCbdt(np) = dCbdt(np-1);
dCpdt(np) = dCpdt(np-1);
% Combine rates for the function to return to the calling routine
dcdt = [dCbdt; dCpdt];
end % of rates function
Alan Stevens
am 15 Mai 2022
All I can suggest is check your data carefully. I don't know enough about your system to connect the data values in your graph.png to the constants in the program.
Siehe auch
Kategorien
Mehr zu Mathematics and Optimization 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!