How to use ODE function for modelling tanks in series

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

Answer 1

[Torsten]

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)

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

[Andrew Rutter]

Thanks Torsten and Star Strider.

Now I understand - you've added n into the function, that allows you to define an equation for the first tank using xa_in and the subsequent ones to use xa(i-1) for i=2:n.

The output looks as I would expect (i think) as the 5th tank takes a while to see feed so starts to drop in concentration later.....

many thanks.

I shall have a play....happy holidays.

[Jaime Diaz]

Thanks Torsten and Star Strider

How will be the equation if we consider the production of a product

Thanks

Jaime

[Torsten]

Since b is produced, you most probably mean

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,1-xa(:,1),'r',t,1-xa(:,2),'g',t,1-xa(:,3),'b',t,1-xa(:,4),'c',t,1-xa(:,5),'k');
% Solve for steady-state
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 = 5x1

0.3919 0.4811 0.5605 0.6265 0.6786

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

hold on
plot(tend,1-xa_ss(1),'o','Color','r');
plot(tend,1-xa_ss(2),'o','Color','g');
plot(tend,1-xa_ss(3),'o','Color','b');
plot(tend,1-xa_ss(4),'o','Color','c');
plot(tend,1-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

[Jaime Diaz]

thanksTorsten

what i would like to know is how can i introduce the ode's for the reactive and the product in the last function

[Jaime Diaz]

for example A->B in the tank series

[Jaime Diaz]

In Berkeley Madonna i have this

METHOD RK4

STARTTIME = 0

STOPTIME=40

DT = 0.02

{N Reactores de Flujo Continuo y Mezcla Completa conectados en serie}

{Constantes}

K=0.2 {Constante de velocidad de primer orden, 1/h. }

Vt=10 {Volumen total de los N reactores, m^3}

V=Vt/N {Volumen por reactores, m^3}

F=1 {Caudal, m^3/h}

N=3 {Número de reactores >2, -}

{Condiciones iniciales}

INIT CA[1] = 0

INIT CA[2..N] = 0

INIT CB[1] = 0

INIT CB[2..N] = 0

{Balances de masa}

d/dt(CA[1]) = (CA0-CA[1])*F/V-K*CA[1]

d/dt(CA[2..N]) = (CA[i-1]-CA[i])*F/V-K*CA[i]

d/dt(CB[1]) = (CB0-CB[1])*F/V+K*CA[1]

d/dt(CB[2..N]) = (CB[i-1]-CB[i])*F/V+K*CA[i]

{Concentraciones de A y B entrando, g/m^3}

CA0 = 2

CB0=0

[Torsten]

tstart = 0;
tend = 40; % h
N = 5; %number of tanks
F = 1; %flow in m^3/h
Vt = 10; %total volume
V = Vt/N; % Volume per tank [m^3]
K = 0.2; % Reaction constant [1/h]
CA0 = 2; 
CB0 = 0; %Inflow concentrations [g/m^3]
[T,C] = ode15s(@(t,c)cstr(t,c,N,F,V,K,CA0,CB0),[tstart tend],zeros(2*N,1));
CA = C(:,1:N);
CB = C(:,N+1:2*N);
CSS = fsolve(@(c)cstr(0,c,N,F,V,K,CA0,CB0),zeros(2*N,1));

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.

CASS = CSS(1:N);
CBSS = CSS(N+1:2*N);
figure(1)
hold on
plot(T,CA(:,1),'r',T,CA(:,2),'g',T,CA(:,3),'b',T,CA(:,4),'c',T,CA(:,5),'k')
plot(T(end),CASS(1),'o','Color','r')
plot(T(end),CASS(2),'o','Color','g')
plot(T(end),CASS(3),'o','Color','b')
plot(T(end),CASS(4),'o','Color','c')
plot(T(end),CASS(5),'o','Color','k')
hold off
grid on

figure(2)
hold on
plot(T,CB(:,1),'r',T,CB(:,2),'g',T,CB(:,3),'b',T,CB(:,4),'c',T,CB(:,5),'k');
plot(T(end),CBSS(1),'o','Color','r')
plot(T(end),CBSS(2),'o','Color','g')
plot(T(end),CBSS(3),'o','Color','b')
plot(T(end),CBSS(4),'o','Color','c')
plot(T(end),CBSS(5),'o','Color','k')
hold off
grid on

function dy = cstr(t,c,N,F,V,K,CA0,CB0)
  CA = c(1:N);
  CB = c(N+1:2*N);
  dCA(1) = (CA0-CA(1))*F/V-K*CA(1);
  dCB(1) = (CB0-CB(1))*F/V+K*CA(1);
  dCA(2:N) = (CA(1:N-1)-CA(2:N))*F/V-K*CA(2:N);
  dCB(2:N) = (CB(1:N-1)-CB(2:N))*F/V+K*CA(2:N);
  dy = [dCA,dCB].';
end

[Jaime Diaz]

Dear Torsten

Great...many thanks!

[Jaime Diaz]

Dear Torsten

An additional question...how can i introduce the steady state solution both parameters with the fsolve?

[Torsten]

An additional question...how can i introduce the steady state solution both parameters with the fsolve?

Exactly analogous to the original problem (see above).

It would have been a good exercise for you to get used to MATLAB, but unfortunately, you missed it.

[Jaime Diaz]

Sorry i am rather new using matlab, could you please give the example for the calculation of A and B in steady state...thanks

[Torsten]

"fsolve" solves dCA(1:N)/dt = 0 and dCB(1:N)/dt = 0 for CA(1:N) and CB(1:N).

[Jaime Diaz]

Torsten

i saw it already...many thanks

[Jaime Diaz]

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]

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.

Answer 2

[Jaime Diaz]

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?

[Torsten]

Works under R2024a (see above).

[Jaime Diaz]

Sorry still the same error

[Torsten]

Restart MATLAB, load the above code in the editor and click on the green arrow to run.

If it still doesn't work: which MATLAB version do you use ?

[Jaime Diaz]

I did it....Matlab 2024a...it is strange...i put the code in Matlab on line and it works

[Steven Lord]

Please show us the full and exact text of any warning and/or error messages you received when you ran the code (all the text displayed in orange and/or red in the Command Window.) The exact text (no paraphrasing, no summarizing, no trimming, just a direct copy and paste) may be useful in determining what's going on and how to avoid the warning and/or error.

[Jaime Diaz]

Error using pdepe

Too many input arguments.

Error in example (line 6)

sol=pdepe(m,@ecuacion1pde,@ecuacion1ic,@ecuacion1bc,x,t);

[Torsten]

Maybe you created your own .m-file with the name "pdepe.m" on your MATLAB path ?

What do you get if you type

which pdepe

[Jaime Diaz]

example

Error using pdepe

Too many input arguments.

Error in example (line 6)

sol=pdepe(m,@ecuacion1pde,@ecuacion1ic,@ecuacion1bc,x,t);

[Torsten]

This is not an answer to my question.

What do you get if you type

which -all pdepe

/MATLAB/toolbox/matlab/funfun/pdepe.m

? More than this path ?

[Jaime Diaz]

which -all pdepe

C:\Users\Jaime Díaz\Documents\MATLAB\pdepe.m

C:\Program Files\MATLAB\R2024a\toolbox\matlab\funfun\pdepe.m % Shadowed

[Torsten]

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]

Great!! It was solved!....many thanks