# How to use ODE function for modelling tanks in series

14 visualizzazioni (ultimi 30 giorni)
Andrew Rutter il 28 Dic 2022
Commentato: Jaime Diaz il 24 Mag 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 CommentiMostra -2 commenti meno recentiNascondi -2 commenti meno recenti

Accedi per commentare.

### Risposta accettata

Torsten il 28 Dic 2022
Modificato: Torsten il 28 Dic 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');
xa0 = ones(n,1);
t = 0;
xa_ss = fsolve(@(x)cstr1(t,x,n,xa_in),xa0)
Equation solved. fsolve completed because the vector of function values is near zero as measured by the value of the function tolerance, and the problem appears regular as measured by the gradient.
xa_ss = 5×1
0.3919 0.4811 0.5605 0.6265 0.6786
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 CommentiMostra 13 commenti meno recentiNascondi 13 commenti meno recenti
Jaime Diaz il 14 Apr 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 il 14 Apr 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.

Accedi per commentare.

### Più risposte (1)

Jaime Diaz il 23 Mag 2024
Modificato: Torsten il 23 Mag 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 CommentiMostra 10 commenti meno recentiNascondi 10 commenti meno recenti
Torsten il 24 Mag 2024
Modificato: Torsten il 24 Mag 2024
Then you should rename
C:\Users\Jaime Díaz\Documents\MATLAB\pdepe.m
because MATLAB tries to call this function instead of the MATLAB integrator.
What's the content of the file in your working directory ?
Jaime Diaz il 24 Mag 2024
Great!! It was solved!....many thanks

Accedi per commentare.

### Categorie

Scopri di più su Ordinary Differential Equations in Help Center e File Exchange

R2020a

### Community Treasure Hunt

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

Start Hunting!

Translated by