How to use ODE function for modelling tanks in series
15 次查看(过去 30 天)
显示 更早的评论
Hi,
I'm learning how to use functions and ODEs and am trying to recreate a tanks in series model that i made in a programme called Berkeley Madona many years ago. It's purpose is to model a series of cascaded stirred tank reactors where the output of one is passed to the next. Hence i need to use a series of ODEs and pass the output of one to the next. The cascade of CSTRs is used to approximate a PFR.
The basic is N tanks in series of identical volume V. A reacting to B and B reacting back to A.
The overall model is based on Octave Levinspiel's Chemcial Reaction Engineering Text, tanks in series model, and some good starting point for a single tank on youtube.
looking forward to your hints and guidance for a matlab novice...many thanks.
For the first tank....the function cstr1.m
function xa_dot = cstr1(t,xa,xa_in)
xb = 1 - xa; % a + b = 1, no other components.
F= 0.1; %flow in m3/s
V = 0.1; %tank vol m3
k1 = 0.45; % rate constant 1
k2 = 0.1; % rate constant 2
xa_dot = (F/V * (xa_in-xa)) + k1*xa*xb - k2 * xa^2;
%function to model tanks in series 1..N
% Molar Balance
%Tank 1 d/dt(Ca[1] = (F/V)*(Ca_in - Ca[1]) - Ra[1]
%Tank 2..N d/dt(Ca[N] = (F/V*(Ca[N-1]-Ca[N]) - Ra[N]
% Ca = XaC, xa + xb = 1, If C = 1 then Ca = xa (X = molar fracctional conversion)
%hence dxa/dt = (F/V * (xa[N-1] - xa[N]) - Ra[N]
% ra = sum of making a (k1*xa*xb) - destroying a by reverse reaction (k2*xa^2).
% xa_dot = d(xa)/dt.
main.m to model 1 tank....
clear all; close all; clc
xa_in = 0.3; % inital concentration of A flowing in.
xa = 1.0; % inital concentration of A in the reactor 1.
run = 5; %runtime seconds
[t1,xa1] = ode15s(@(t,x)cstr1(t,x,xa_in),[0 run],xa); % (Odefun, timespan, y0)
%xa1 represents the conversion in the first tank....
The first tank works and I can calculate xa1, I now need to find xa in the tanks 2..N, how do I do this...?
Each tank N uses the input from the previous (N-1) tank for the concentration of A and B.
Do I create a second function, or nest it in the first? and I presume I need to initalise the values of xa(2..N)
I'm not sure how to set up a function for a matrix of outputs 2..N, and in the main function how to relate the CSTR1 with CSTR2..N
0 个评论
采纳的回答
Torsten
2022-12-28
编辑:Torsten
2022-12-28
xa_in = 0.3; % inital concentration of A flowing in.
xa = 1.0; % inital concentration of A in the reactor 1.
tend = 15; %runtime seconds
n = 5; % number of tanks
[t,xa] = ode15s(@(t,x)cstr1(t,x,n,xa_in),[0 tend],xa*ones(n,1)); % (Odefun, timespan, y0)
plot(t,xa(:,1),'r',t,xa(:,2),'g',t,xa(:,3),'b',t,xa(:,4),'c',t,xa(:,5),'k');
% Solve for steady-state
xa0 = ones(n,1);
t = 0;
xa_ss = fsolve(@(x)cstr1(t,x,n,xa_in),xa0)
hold on
plot(tend,xa_ss(1),'o','Color','r');
plot(tend,xa_ss(2),'o','Color','g');
plot(tend,xa_ss(3),'o','Color','b');
plot(tend,xa_ss(4),'o','Color','c');
plot(tend,xa_ss(5),'o','Color','k');
function xa_dot = cstr1(t,xa,n,xa_in)
xb = 1 - xa; % a + b = 1, no other components.
F = 0.1; %flow in m3/s
V = 0.1; %tank vol m3
k1 = 0.45; % rate constant 1
k2 = 0.1; % rate constant 2
xa_dot = zeros(n,1);
xa_dot(1) = F/V * (xa_in-xa(1)) + k1*xa(1)*xb(1) - k2 * xa(1)^2;
for i = 2:n
xa_dot(i) = F/V * (xa(i-1)-xa(i)) + k1*xa(i)*xb(i) - k2*xa(i)^2;
end
end
15 个评论
Jaime Diaz
2024-4-14
Dear Torsten
why have i to use the traspose dy = [dCA,dCB].'; here..could i use a zero vector column and if so how will be the instruction?
Torsten
2024-4-14
why have i to use the traspose dy = [dCA,dCB].'; here
The vector of time derivatives must be a column vector for MATLAB's integrators.
could i use a zero vector column and if so how will be the instruction?
zeros(n,1) gives you a column vector of zeros of length n.
更多回答(1 个)
Jaime Diaz
2024-5-23
编辑:Torsten
2024-5-23
Dear Torsten
I have a problem usin the pdepe of Matlab...it says toomany input arguments using pdepe
m=0;
x=linspace(0,3,25);
tend=0.7;
t=linspace(0,tend,25);
sol=pdepe(m,@ecuacion1pde,@ecuacion1ic,@ecuacion1bc,x,t);
C=sol(:,:,1);
plot (x,C(end,:));
xlabel('x(m)');
ylabel('C(g/m^3)');
hold on
function [g,f,s]=ecuacion1pde(x,t,C,DcDx)
E=3;
v=10;
kd=10;
g=1;
f=E*DcDx;
s=-v*DcDx-kd*C;
end
function C0=ecuacion1ic(x)
C0=0;
end
function [pl,ql,pr,qr]=ecuacion1bc(xl,cl,xr,cr,t)
Ce=4;
pl=cl-Ce;
ql=0;
E=3;
pr=0;
qr=1/E;
end
Could you give some indication about the problem?
12 个评论
另请参阅
类别
在 Help Center 和 File Exchange 中查找有关 Ordinary Differential Equations 的更多信息
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!