# how can I code a system of 2 pdes and an ode for modelling a packed bed

8 views (last 30 days)
yb on 29 Apr 2020
Commented: yb on 3 May 2020
Hi community
I want to code a system of 2 odes and a pde for modelling a packed bed. I have already the breakthrough curve that should result but I got a different one. I want to know where is the problem in my code?
I descretized the variable c in space with the method of lines and I used the ode45 solver to solve the system of odes.
You will find in the attached file the sytem of equations , my code and the breakthrough that should result.
yb on 1 May 2020
Hi bill
thank you for your clarification
can you help me to find the problem in my code, I don't know why I get a different breakthrough

Bill Greene on 1 May 2020
I did some experimentation with your code. Here are a couple of comments related to the changes I made:
1. Because the equations are first-order in z, you can specify a boundary condition only at the inlet end. You have to be careful not to specify a condition at the outlet.
2. The number of time points you define for ode45 don't affect the accuracy of the solution; they simply indicate where you would like to see the output. So, typically, you need only enough points to generate a reasonably smooth plot.
My modified version of your code is shown below. The output still doesn't agree with your plot but at least the trends are right. I did not try to verify that your finite difference code agrees with your equations. As you have written it, equation two is VERY difficult to verify; I suggest breaking it into sub-expressions.
Us=1.47*10^-3;
Epsilon=0.35;
ap=6683.3;
kc=1.057*10^-5;
b=0.1461;
qm=93.69;
l=0.075*10^-3;
Ds=2*10^-13;
Ro=374;
cFeed = 150; % Feed concentration
L = 0.12; % Column length
t0 = 0; % Initial Time
tf = 80000; % Final time
tf=1.5e5;
dt = 10; % Time step
dt=1000;
z = [0:0.012:L]; % Mesh generation
n=21;
z=linspace(0,L,n);
t = [t0:dt:tf];% Time vector
nt=100;
t=linspace(t0, tf, nt);
n = numel(z); % Size of mesh grid
c0 = zeros(n,1);% t = 0, c = 0
c0(1)=cFeed;
qb0 = zeros(n,1); % t = 0, q = 0 for all z,
qb0(1) = zeros;
qs0 = zeros(n,1);
qs0(1) = zeros;
y0=[c0 ; qb0 ; qs0];
% Appends conditions together
[t, y] = ode45(@(t,y) f(t,y,z,n),[t0 tf],y0);
figure; plot(t,y(:,n));
xlabel('time')
ylabel('exit conc.')
figure; plot(z, y(end,1:n));
title('C at final time');
end
function dydt=f(t,y,z,n)
Us=0.00024;
Epsilon=0.35;
ap=6661.81;
kc=3.96*10^-6;
b=0.1461;
qm=93.69;
l=7.5*10^-5;
Ds=2*10^-13;
Ro=373.68;
dcdt=zeros(n,1);
dqbdt=zeros(n,1);
dqsdt=zeros(n,1);
%dydt=zeros(3*n,1);
dcdz=zeros(n,1);
c=y(1:n);
qb=y(n+1:2*n);
qs=y(2*n+1:3*n);
%for i=2:n-1
for i=2:n
dcdz(i)= (c(i)-c(i-1))/(z(i)-z(i-1));
end
%dcdz(n) = 0;
%dcdt(1) = 0;
for i=2:n
dqbdt(i)=(kc*ap/((1-Epsilon)*Ro))*(c(i)-(qs(i)/(b*(qm-qs(i)))));
dqsdt(i)=((((l^2)*kc*ap)/(3*Ds*Ro*(1-Epsilon)))*(((-Us/Epsilon)*dcdz(i))+((((3*Ds)/(l^2))-((kc*ap)/Epsilon))*(c(i)-(qs(i)/(b*(qm-qs(i))))))))/(1+(((l^2)*kc*ap)/(3*Ds*Ro*(1-Epsilon)*b*(qm-qs(i))))*(1+(qs(i)/(qm-qs(i)))));
dcdt(i)=((-Us/Epsilon)*dcdz(i))-((kc*ap)/Epsilon)*(c(i)-(qs(i)/(b*(qm-qs(i)))));
end
dydt=[dcdt;dqbdt;dqsdt];
end
yb on 3 May 2020
Hi darova
thank you very much, but I think we should n't give the same name for different expressions (for example c1)
c1 = kc*ap/(1-Epsilon)/Ro;
c1 = ((l^2)*kc*ap)/(3*Ds*Ro*(1-Epsilon));
c1 = -Us/Epsilon*dcdz(i);