Modelování 6A: Polovsádkový reaktor (DAE)

Po vsádkové destilaci se vracíme k dalšímu příkladu řešení diferenciálně-algebraických rovnic (DAE). Odvodíme model polovsádkového reaktoru, ukážeme řešení této soustavy v Matlabu a zkusíme zjistit DAE-index systému.

Model polovsádkového reaktoru

Mějme polovsádkový reaktor na počátku obsahující pouze složku A, do kterého přidáváme konstantní rychlostí složku B. V reaktoru probíhá vratná reakce A + B <--> C + D. Důležitý detail je, že pro určení hustoty směsi budeme předpokládat ideální chování, tj. platnost aditivity objemů jednotlivých složek.

Základem modelu budou látkové bilance všech složek, které postupně upravíme

: $\begin{aligned}\displaystyle \frac{\mathrm{d}n_i}{\mathrm{d}t} &= \dot{n}_i^\mathrm{in} + r_i V \\ \frac{\mathrm{d}}{\mathrm{d}t}(c_i V) &= c_i^\mathrm{in}\dot{V}^\mathrm{in} + r_i V \\ V\frac{\mathrm{d}c_i}{\mathrm{d}t}+c_i\frac{\mathrm{d}V}{\mathrm{d}t} &= c_i^\mathrm{in}\dot{V}^\mathrm{in} + r_i V\end{aligned}$

Pro N složek máme N+1 stavových proměnných: objem vsádky v reaktoru V a N koncentrací složek c_i. Dalšími veličinami v rovnicích jsou objemový průtok přítoku do reaktoru, V̇ⁱⁿ, koncentrace složek v tomto proudu, c_iⁱⁿ a rychlosti vzniku/zániku složky součtem dopředné a zpětné reakce, r_i.

Pro N+1 stavových proměnných potřebujeme N+1 rovnic. Zatím máme pouze N bilancí. Chybí nám jedna rovnice a tou je vztah mezi objemem směsi a koncentracemi jednotlivých složek (stavová rovnice). Tu odvodíme z předpokladu platnosti aditivity objemů

: $\displaystyle V = \sum_{i=1}^N V_i$ .

Za předpokladu platnosti následujících jednoduchých vztahů

: $\displaystyle V_i = \frac{m_i}{\rho_i}\quad m_i = n_i M_i = c_i V M_i$ ,

dostaneme postupně

: $\displaystyle V = \sum_{i=1}^N V_i = \sum_{i=1}^N \frac{m_i}{\rho_i} = V\sum_{i=1}^N \frac{c_i M_i}{\rho_i}$ ,

což po vydělení objemem dá hledanou rovnici

: $\displaystyle 1 = \sum_{i=1}^N \frac{c_i M_i}{\rho_i}$ .

Jelikož platí

: $\displaystyle c_i M_i = \frac{m w_i}{V}$ ,

je snadné stavovou rovnici přeformulovat do alternativního tvaru

: $\displaystyle\frac{1}{\rho}= \sum_{i=1}^N \frac{w_i}{\rho_i}$

což jinými slovy znamená, že měrný objem (tj. převrácená hodnota hustoty) směsi je hmotnostními zlomky vážený průměr měrných objemů jednotlivých složek. Jelikož pracujeme s koncentracemi složek a nikoli s hmotnostními zlomky, je ale původní tvar stavové rovnice pro naše účely výhodnější.

Dodejme, že kdybychom potřebovali z jakéhokoliv důvodu přidat hustotu směsi jako další stavovou proměnnou, museli bychom matematický model doplnit o bilanci objemu, kterou odvodíme z hmotnostní bilance

: $\begin{aligned}\displaystyle \frac{\mathrm{d}m}{\mathrm{d}t} &= \dot{m}^\mathrm{in}\\ \frac{\mathrm{d}}{\mathrm{d}t}(\rho V)&=\rho^\mathrm{in}\dot{V}^\mathrm{in}\\ V\frac{\mathrm{d}\rho}{\mathrm{d}t} + \rho\frac{\mathrm{d}V}{\mathrm{d}t} &= \rho^\mathrm{in}\dot{V}^\mathrm{in}\end{aligned}$

Model bude fungovat stejně dobře s nebo bez této bilance objemu a hustoty směsi jako stavové proměnné. Mnohem jednodušší je ale pro hustotu směsi využít tento vztah

: $\displaystyle\rho = \frac{m}{V} = \sum_{i=1}^N \frac{m_i}{V} = \sum_{i=1}^N \frac{n_i M_i}{V} = \sum_{i=1}^N c_i M_i$

Implementace modelu v Matlabu

Nejdříve je potřeba napsat funkci, která vypočítá residuál F(t,y,y‘). Tady je příklad této funkce pro model polovsádkového reaktoru.

function res = res_semibatch(~,y,yp,pars)
%RES_SEMIBATCH Calculate residual for semi-batch reactor model
%   Semi-batch reactor model
%      RES = RES_SEMIBATCH(~,Y,YP,PARS)
%   Y, YP - vectors of state variables and their derivatives
%      [cA, cB, cC, cD, V]
%   PARS - model parameters
%      [MA, MB, MC, MD, RHOA, RHOB, RHOC, RHOD, CAIN, CBIN, CCIN, CDIN, ...
%       VIN, K1, K2]
 
% unpack parameters
M   = pars(1:4);      % molar weight [kg-i/mol-i]
rho = pars(5:8);      % density      [kg-i/m3-i]
cin = pars(9:12);     % concentration in inflow stream [mol-i/m3]
Vin = pars(13);       % inflow stream flow rate [m3/s]
k1  = pars(14);       % forward reaction rate coefficient [mol/m3.s]
k2  = pars(15);       % backward reaction rate coefficient [mol/m3.s]
 
% unpack state variables ...
c = y(1:4);        % concentration [mol-i/m3]
V = y(5);          % liquid in reactor volume [m3]
 
% ... and their derivatives
dcdt = yp(1:4);
dVdt = yp(5);
 
% calculate reaction rate
r1 = k1*c(1)*c(2); % forward reaction
r2 = k2*c(3)*c(4); % backward reaction
 
% calculate residuals
% mole balances
res(1) = Vin*cin(1) + (-r1+r2)*V - c(1)*dVdt - V*dcdt(1);
res(2) = Vin*cin(2) + (-r1+r2)*V - c(2)*dVdt - V*dcdt(2);
res(3) = Vin*cin(3) + (+r1-r2)*V - c(3)*dVdt - V*dcdt(3);
res(4) = Vin*cin(4) + (+r1-r2)*V - c(4)*dVdt - V*dcdt(4);
% equation of state
res(5) = c(1)*M(1)/rho(1) + c(2)*M(2)/rho(2) + c(3)*M(3)/rho(3) ...
       + c(4)*M(4)/rho(4) - 1;
 
res = res'; % make residual column vector
end

Dalším krokem bude hlavní skript, který uvedeme po částech. Začneme zadáním konkrétních hodnot parametrů, tj. fyzikální vlastnosti jednotlivých složek a složení a množství přítoku do reaktoru.

% semi-batch reactor (main script)
%
% species: (1)-CH3COOH, (2)-C2H5OH, (3)-CH3COOC2H5, (4)-H2O
 
% model parameters
M   = [60.05 46.07 88.11 18.02] * 1e-3; % kg/mol
rho = [1044.6 789.3 900.3 1000];        % kg/m3
cin = [0, rho(2)/M(2), 0, 0];           % mol/m3
Vin = 0.4e-3/60;                        % m3/s
k1  = 8.00e-9;                          % m3/mol.s
k2  = 2.55e-9;
 
pars = [M, rho, cin, Vin, k1, k2];      % 15 parameters

Nyní je potřeba určit počáteční podmínky. U DAE systémů je to komplikovanější, protože kromě hodnot stavových proměnných potřebujeme určit i jejich derivace. Hodnoty stavových proměnných nejsou navíc všechny nezávislé, ale musí splňovat vazné podmínky dané modelem. V našem případě například koncetrace musí splňovat stavovou rovnici, tj. u 4 složkové směsi můžeme zafixovat koncentrace pouze 3 složek a koncentraci poslední složky musíme jen odhadnout a nechat dopočítat. Protože víme, že v reaktoru je na počátku pouze složka A, fixujeme koncetrace ostatních složek na nule a c_A(0) necháme dopočítat.

% initial conditions (known and being fixed)
cB0 = 0;
cC0 = 0;
cD0 = 0;
V0 = 500e-3; % m3
% note: cA0 canot be fixed, because SUM(ci0*Mi/rhoi) = 1 must be held
 
% initial conditions (guessed and being calculated)
time0 = 0.;
timef = 48*3600; % s
cAg = rho(1)/M(1); % mol/m3
 
% calculate consistent initial conditions
options = odeset('RelTol',1e-8);
yg     = [cAg, cB0, cC0, cD0, V0];
ygfix  = [0,   1,   1,   1,   1 ];
 
ypg    = [0.0, 0.0, 0.0, 0.0, 0.0];
ypgfix = [0,   0,   0,   0,   0  ];
 
[y0,yp0] = decic(@(t,y,yp) res_semibatch(t,y,yp,pars), ...
                 time0, yg, ygfix, ypg, ypgfix, options)
res_semibatch(0,y0,yp0,pars) % check calculated initial conditions

K výpočtu počátečních podmínek jsme použili funkci decic a na konci zkontrolovali, že nalezené podmínky dávají nulový residuál, tj. jsou konzistentní. Následuje výpočetně nejnáročnější část: vlastní integrace soustavy pomocí funkce ode15i:

% integrate DAE
[tt,yy] = ode15i(@(t,y,yp) res_semibatch(t,y,yp,pars), ...
                 [time0 timef], y0, yp0, options);

Poslední část kódu se zabývá vizualizací získaných výsledků. Pokud chceme znát, jak se měnila hustota směsi, můžeme ji dopočítat ze znalosti stavových proměnných dle výše uvedeného vztahu.

% post-processing (calculate density of reaction mixture)
cA = yy(:,1);
cB = yy(:,2);
cC = yy(:,3);
cD = yy(:,4);
V = yy(:,5);
rho = cA*M(1) + cB*M(2) + cC*M(3) + cD*M(4); % kg/m3
 
% plot graphs
tt = tt/3600; % h
subplot(2,2,1) % concentration
plot(tt,cA,tt,cB,tt,cC,tt,cD)
xlabel('time / h')
ylabel('concentration / mol/m3')
legend('CH3COOH','C2H5OH','ACETATE','H2O')
 
subplot(2,2,3) % density
plot(tt,rho)
xlabel('time / h')
ylabel('density / kg/m3')
 
subplot(2,2,4) % volume
plot(tt,V)
xlabel('time / h')
ylabel('volume / m3')
% end of main script

Na další stánce se pokusíme určit DAE-index systému.

Pages: 12

Modelování 6A: Polovsádkový reaktor (DAE)

Model polovsádkového reaktoru

Implementace modelu v Matlabu

Published by

Zdenek Grof

Leave a Reply Cancel reply