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

3 次查看（过去 30 天）
yb
yb 2020-4-29
评论： yb 2020-5-3
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.
  5 个评论
显示 3更早的评论
Bill Greene
Bill Greene 2020-5-1
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
yb 2020-5-1
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

请先登录，再进行评论。

请先登录，再回答此问题。

回答（1 个）

Bill Greene
Bill Greene 2020-5-1
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
  4 个评论
显示 2更早的评论
darova
darova 2020-5-3
编辑：darova 2020-5-3
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
yb 2020-5-3
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);

请先登录，再进行评论。

类别

在 Help Center 和 File Exchange 中查找有关 Boundary Conditions 的更多信息

标签

Community Treasure Hunt

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

Start Hunting!

Translated by