Adsorption Modelling - Solving PDE - Axial Dispersion Model but STEADY STATE
3 vues (au cours des 30 derniers jours)
Afficher commentaires plus anciens
Hi,
I am modelling a moving bed adsorber (with activated carbon as the adsorbent) for water treatment as a sand filter. Instead of sand in the filter, it wil be filled with activated carbon beads and the water will be moving upwards while the activated carbon moves downwards (countercurrent flow).
To model this, I have made the Axial Dispersion Model steady state. Therefore, the equations become:
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/526079/image.png)
I assumed that the dq/dt isn't cancelled out.
I'm unsure how to write matlab code to solve for c and q and to plot c against z.
I have the code for a fixed bed adsorber but I am not sure how to change it to fit my model.
function test
%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
b = 2.5*10^-5; % Langmuir parameter
qs = 2*10^4; % Saturation capacity
cFeed = 10; % Feed concentration
L = 1; % Column length
t0 = 0; % Initial Time
tf = 1000; % Final time
dt = 0.5; % Time step
z = [0:0.01: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)
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
b = 2.5*10^-5; % Langmuir parameter
qs = 20000; % Saturation capacity
% 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));
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*( ((qs*b*c(1))/(1 + b*c(1))) - q(1) );
DcDt(1) = 0.0;
% Set time derivatives in remaining points
for i=2:n
%Equation: dq/dt = k(q*-q) where q* = qs * (b*c)/(1+b*c)
DqDt(i) = k*( ((qs*b*c(i))/(1 + b*c(i))) - 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
0 commentaires
Réponses (1)
Shadaab Siddiqie
le 23 Fév 2021
From my understanding you want to solve an multi order differential equations. You can go through symbolic variables and expressions and also dsolve for more information.
0 commentaires
Voir également
Catégories
En savoir plus sur Particle & Nuclear Physics dans Help Center et File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!