Fixed Bed Adsorption Column Using Dimensionless Equations

My group and I are trying to write a code for fixed bed adsorption with dimensionless equations and Damkohler's number. I attached an image of the equations we're trying to use and variables. We have a general idea and some code written down of what the variables are meaning but would like some help on where to go with the ODE/PDE code as we're not familiar with them. Any help, guidance, or advice would be greatly appreciated. Thank you very much in advance.

Answer 1

[Torsten]

The usual procedure is to define a grid

0 = X(1) < X(2) < ... < X(n) = 1

approximate the spatial derivative dA/dX in grid point i by

dA(i)/dx ~ (A(i)-A(i-1))/(X(i)-X(i-1)), 

keep the time derivatives in their continuous form,

write the PDE system as a system of 2*n ordinary differential equations

A(1) = A_in
dA(i)/dt = - (A(i)-A(i-1))/(X(i)-X(i-1)) - q0/A_in * s(q(i),A(i))        2<=i<=n
dq(i)/dt = s(q(i),A(i))                                                  1<=i<=n

and use ode15s to solve the system.

If you need further details, look up "method-of-lines".

[tori lansing]

Torsten, Thank you for the help! We just talked with our professor more and he advised us to not use ODE code since these are PDE's. He stated that using finite difference and midpoint slopes would be better. We're trying to do some code with that now, but do you have any other advice or guidance? Thank you again.

[Torsten]

The usual way to solve PDEs is to transform them to a system of ODEs in the spatial discretization points. So you apply finite differencing to the spatial derivatives, leave the temporal derivatives unchanged and solve the resulting system of ordinary differential equations in time by an ODE integrator.

I don't know exactly which method your professor is referring to, but maybe he means to discretize simultaneously in time and space and advance the solution by a fixed time step dt.

See Explicit Method, Implicit Method or Crank-Nicolson Method here:

https://en.wikipedia.org/wiki/Finite_difference_method

[tori lansing]

Torsten, that last sentence "discretize siimultaneously in time and space and advance the solition by a fixed time step dt" is exactly what we need to do. We also need to define paramaters Da, K, & q0/Ain to change in the equation. Right now we are looking at using gradients and numerical differentiation techniques out of the book and we have little but of code done and working. Are you saying that eventually we will need to use and ODE solver? If so can you give more guidance as to what you mean/where we should go? Also thank you for the link to understanding finite difference method. Thanks again for all your help!

[Torsten]

Right now we are looking at using gradients and numerical differentiation techniques out of the book and we have little but of code done and working.

I already did this for you. For an explicit method see the Wiki article under "heat equation, explicit method".

(A_j^(n+1)-A_j^(n))/dt = -(A_j^(n)-A_(j-1)^(n))/dx - q0/Ain *s(q_j^(n),A_j^(n))
(q_j^(n+1)-q_j^(n))/dt = s(q_j^(n),A_j^(n))

Solve for A_j^(n+1) and q_j^(n+1) , and you can start immediately.

[tori lansing]

Torsten, First I want to thank you or all your help and quick responses, it's helping us a lot. Secondly, when we're running what youve given us with our parameters we're running into some isues. Can you take a look and see if maybe we're defining something wrong? We are getting an error that the statement is incomplete. Again, Thank you so much! We appreciate all your help.

tf=3000;
dt=3000/200;
tspan=0:dt:tf; % creates array for length of time
cf=0.015; %initial concentration
nz=200;
np=nz+1;
Cb=zeros(np,1); Cb(1)=cf;
Cp=zeros(np,1);
C0=[Cb;Cp]; % this may need some adjusting, im not sure if this comes out right
% I beleive this establishes the things for the ODE solver
% following the ode solver, the plots are created
% then a function is created for rates that is then recalled in ode solver
% then a for loop to establish the rest of the equations
A_in=input("Input a value for A_in:")
if (A_in>0)&&(A_in<1)
    disp(A_in)
else
    disp('Invalid value for inlet concentration')
end
K=input('Enter a value for the rate constant:')
q_o=input('Enter a value for the total exchange capacity of the system:')
k=0.001
L=10; % meters
K=0.001;
A=0:0.01:A_in;
Ac=A/A_in;
V=0.001;
t=L/V;
delta=tspan./t;
Da=k*A_in*t;
dqdt=zeros(np,1);
s=q_o/A_in;
x=0:0.01:10;
X=x./L;
0 == X(1) < X(2) < ... < X(n) = 1

This statement is incomplete.

[Torsten]

To get a foundation in MATLAB programming, I suggest the MATLAB Onramp course to learn about MATLAB basics in an online course free of costs:

https://matlabacademy.mathworks.com/details/matlab-onramp/gettingstarted

[Torsten]

Here is a short code with the explicit method for the problem:

dy/dt + a*dy/dx = 0
a > 0
0 <= x <= 1
y(x,t=0) = 0 for all x
y(x=0,t) = 1 for all t

It's already very similar to yours.

a = 0.1;
tstart = 0.0;
tend = 2.0;
nt = 1000;
xstart = 0.0;
xend = 1.0;
nx = 1000;
x = linspace(xstart,xend,nx).';
dx = x(2)-x(1);
t = linspace(tstart,tend,nt);
dt = t(2)-t(1);
y = zeros(nx,nt);
y(:,1) = 0.0;
y(1,:) = 1.0;
for it = 2:nt
  y(2:nx,it) = y(2:nx,it-1) - dt * a * (y(2:nx,it-1)-y(1:nx-1,it-1))./(x(2:nx)-x(1:nx-1));
end
plot(x,[y(:,1),y(:,500),y(:,1000)])