72 views (last 30 days)

Show older comments

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:

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

Shadaab Siddiqie
on 23 Feb 2021

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

Start Hunting!