How to use ODE function for modelling tanks in series

Question

Andrew Rutter 2022-12-28

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1885392-how-to-use-ode-function-for-modelling-tanks-in-series

评论： Jaime Diaz 2024-5-24

IMG_3836.jpg

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 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

Torsten 2022-12-28

1
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1885392-how-to-use-ode-function-for-modelling-tanks-in-series#answer_1137547

编辑：Torsten 2022-12-28

在 MATLAB Online 中打开

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

15 个评论
显示 13更早的评论隐藏 13更早的评论

Torsten 2024-4-5

在 MATLAB Online 中打开

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

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

Torsten 2024-4-6

编辑：Torsten 2024-4-7

在 MATLAB Online 中打开

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 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.

请先登录，再进行评论。

Answer 2

Jaime Diaz 2024-5-23

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1885392-how-to-use-ode-function-for-modelling-tanks-in-series#answer_1462476

编辑：Torsten 2024-5-23

在 MATLAB Online 中打开

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 个评论
显示 10更早的评论隐藏 10更早的评论

Torsten 2024-5-24

编辑：Torsten 2024-5-24

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 2024-5-24

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

请先登录，再进行评论。

How to use ODE function for modelling tanks in series

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

采纳的回答

15 个评论
显示 13更早的评论隐藏 13更早的评论

更多回答（1 个）

12 个评论
显示 10更早的评论隐藏 10更早的评论

另请参阅

类别

标签

产品

版本

Community Treasure Hunt

How to use ODE function for modelling tanks in series

0 个评论 显示 -2更早的评论隐藏 -2更早的评论

采纳的回答

15 个评论 显示 13更早的评论隐藏 13更早的评论

更多回答（1 个）

12 个评论 显示 10更早的评论隐藏 10更早的评论

另请参阅

类别

标签

产品

版本

Community Treasure Hunt

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

15 个评论
显示 13更早的评论隐藏 13更早的评论

12 个评论
显示 10更早的评论隐藏 10更早的评论