How to use ODE function for modelling tanks in series
15 visualizaciones (últimos 30 días)
Mostrar comentarios más antiguos
Andrew Rutter
el 28 de Dic. de 2022
Comentada: Jaime Diaz
el 24 de Mayo de 2024
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 comentarios
Respuesta aceptada
Torsten
el 28 de Dic. de 2022
Editada: Torsten
el 28 de Dic. de 2022
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 comentarios
Jaime Diaz
el 14 de Abr. de 2024
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
el 14 de Abr. de 2024
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.
Más respuestas (1)
Jaime Diaz
el 23 de Mayo de 2024
Editada: Torsten
el 23 de Mayo de 2024
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 comentarios
Ver también
Categorías
Más información sobre Ordinary Differential Equations en Help Center y File Exchange.
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!