- 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.
- 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.
how can I code a system of 2 pdes and an ode for modelling a packed bed
8 views (last 30 days)
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.
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:
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.
cFeed = 150; % Feed concentration
L = 0.12; % Column length
t0 = 0; % Initial Time
tf = 80000; % Final time
dt = 10; % Time step
z = [0:0.012:L]; % Mesh generation
t = [t0:dt:tf];% Time vector
t=linspace(t0, tf, nt);
n = numel(z); % Size of mesh grid
c0 = zeros(n,1);% t = 0, c = 0
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(z, y(end,1:n));
title('C at final time');
%dcdz(n) = 0;
%dcdt(1) = 0;