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

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.

Thanks in advance.

[darova]

I see that and depend only on t (2 ode)

And C is function of 2 variables (t and z)

I have a question: how two solutions ( and C for example ) are connected?

Since C can be surface and is a curve

[yb]

Hi darova

you're right there is 2 odes and only one pde , thanks a lot

the 2 variables c and q are connected since the overall balance for a solute in the bed obeys

[darova]

The question still is not clear: if C changes in z direction and doesn't. How they can be connected?

C can have several different values for only one value of

[Bill Greene]

Technically, all three equations are PDE. does change as a function of z because C is a function of z. This is a common misunderstanding of the definition of PDE.

[yb]

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

thanks in advance

Answer 1

[Bill Greene]

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.

function matlabAnswers_5_1_2020
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]

how can I break it into sub-expressions

[yb]

Hi Bill

Please, can you tell me how can I break my second equation into sub-expressions?

Thanks in advance

[darova]

Im pretty sure that Bill means something like that

% 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
for i = 2:n
    c1 = kc*ap/(1-Epsilon)/Ro;
    c2 = c(i) - qs(i)/b/(qm-qs(i));
    dqbdt(i)=c1*c2;
    c1 = ((l^2)*kc*ap)/(3*Ds*Ro*(1-Epsilon));
    c2 = -Us/Epsilon*dcdz(i);
    c3 = 3*Ds/l^2;
    c4 = kc*ap/Epsilon;
    c5 = qs(i)/b/(qm-qs(i));
    c6 = 1+qs(i)/(qm-qs(i));
    c7 = 3*Ds*Ro*(1-Epsilon)*b*(qm-qs(i));
    c8 = l^2*kc*ap/c7;
    dqsdt(i)=c1*(c2+((c3-c4)*(c(i)-c5)))/(1+c8*c6);
    c1 = -Us/Epsilon*dcdz(i);
    c2 = c(i)-qs(i)/b/(qm-qs(i));
    dcdt(i) = c1 - kc*ap/Epsilon * c2;
end

[yb]

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);