Adsorption Modelling - Solving PDE - Axial Dispersion Model but STEADY STATE

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

Answer 1

[Shadaab Siddiqie]

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.