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

Po vsádkové destilaci ^[1] 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 ^[2] 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 ^[3]:

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

Určení DAE indexu soustavy

Matematický model vsádkového reaktoru pro 4 složky

: $\begin{aligned}\displaystyle &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 \quad i \in (\text{A, B, C, D})\\ &\sum_{i \in (\text{A, B, C, D})} \frac{c_i M_i}{\rho_i} = 1\end{aligned}$

můžeme napsat v maticovém tvaru

: $\displaystyle\mathbf{M}(t,\mathbf{y}) \frac{\mathrm{d}\mathbf{y}}{\mathrm{d}t} = \mathbf{f}(t,\mathbf{y})$

takto

: $\begin{bmatrix} V & 0 & 0 & 0 & c_\mathrm{A} \\ 0 & V & 0 & 0 & c_\mathrm{B} \\ 0 & 0 & V & 0 & c_\mathrm{C} \\ 0 & 0 & 0 & V & c_\mathrm{D} \\ 0 & 0 & 0 & 0 & 0 \end{bmatrix} \begin{bmatrix} \mathrm{d}c_\mathrm{A}/\mathrm{d}t \\ \mathrm{d}c_\mathrm{B}/\mathrm{d}t \\ \mathrm{d}c_\mathrm{C}/\mathrm{d}t \\ \mathrm{d}c_\mathrm{D}/\mathrm{d}t \\ \mathrm{d}V/\mathrm{d}t \end{bmatrix} = \begin{bmatrix} c_\mathrm{A}^\mathrm{in} \dot{V}^\mathrm{in} + r_\mathrm{A}V \\ c_\mathrm{B}^\mathrm{in} \dot{V}^\mathrm{in} + r_\mathrm{B}V \\ c_\mathrm{C}^\mathrm{in} \dot{V}^\mathrm{in} + r_\mathrm{C}V \\ c_\mathrm{D}^\mathrm{in} \dot{V}^\mathrm{in} + r_\mathrm{D}V \\ 1-\sum_{i\in(\text{A, B, C, D})}c_i M_i /\rho_i \end{bmatrix}$ .

Vidíme, že matice M je kvůli poslední rovnici singulární, jedná se tedy o DAE. Index je definován jako počet derivací algebraických rovnic, které je nutné provést, abychom DAE převedli na ODE. Derivace poslední, algebraické, rovnice je

: $\displaystyle \frac{\mathrm{d}c_\mathrm{A}}{\mathrm{d}t} \frac{M_\mathrm{A}}{\rho_\mathrm{A}} + \frac{\mathrm{d}c_\mathrm{B}}{\mathrm{d}t} \frac{M_\mathrm{B}}{\rho_\mathrm{B}} + \frac{\mathrm{d}c_\mathrm{C}}{\mathrm{d}t} \frac{M_\mathrm{C}}{\rho_\mathrm{C}} + \frac{\mathrm{d}c_\mathrm{D}}{\mathrm{d}t} \frac{M_\mathrm{D}}{\rho_\mathrm{D}} = 0$ .

Nahradíme-li v modelu algebraickou rovnici její derivací, dostaneme

: $\begin{bmatrix} V & 0 & 0 & 0 & c_\mathrm{A} \\ 0 & V & 0 & 0 & c_\mathrm{B} \\ 0 & 0 & V & 0 & c_\mathrm{C} \\ 0 & 0 & 0 & V & c_\mathrm{D} \\ M_\mathrm{A}/\rho_\mathrm{A} & M_\mathrm{B}/\rho_\mathrm{B} & M_\mathrm{C}/\rho_\mathrm{C} & M_\mathrm{D}/\rho_\mathrm{D} & 0 \end{bmatrix} \begin{bmatrix} \mathrm{d}c_\mathrm{A}/\mathrm{d}t \\ \mathrm{d}c_\mathrm{B}/\mathrm{d}t \\ \mathrm{d}c_\mathrm{C}/\mathrm{d}t \\ \mathrm{d}c_\mathrm{D}/\mathrm{d}t \\ \mathrm{d}V/\mathrm{d}t \end{bmatrix} = \begin{bmatrix} c_\mathrm{A}^\mathrm{in} \dot{V}^\mathrm{in} + r_\mathrm{A}V \\ c_\mathrm{B}^\mathrm{in} \dot{V}^\mathrm{in} + r_\mathrm{B}V \\ c_\mathrm{C}^\mathrm{in} \dot{V}^\mathrm{in} + r_\mathrm{C}V \\ c_\mathrm{D}^\mathrm{in} \dot{V}^\mathrm{in} + r_\mathrm{D}V \\ 0 \end{bmatrix}$ .

Pomocí Gaussovy eliminace se podívejme, zda je matice M stále singulární. Přičteme-li k poslednímu řádku postupně -M_i/(ρ_iV) násobky všech předchozích řádků, a výsledek násobíme V, získáme

: $\begin{bmatrix} V & 0 & 0 & 0 & c_\mathrm{A} \\ 0 & V & 0 & 0 & c_\mathrm{B} \\ 0 & 0 & V & 0 & c_\mathrm{C} \\ 0 & 0 & 0 & V & c_\mathrm{D} \\ 0 & 0 & 0 & 0 & \sum_{i \in (\text{A,B,C,D})}\frac{c_i M_i}{\rho_i} \end{bmatrix} \begin{bmatrix} \mathrm{d}c_\mathrm{A}/\mathrm{d}t \\ \mathrm{d}c_\mathrm{B}/\mathrm{d}t \\ \mathrm{d}c_\mathrm{C}/\mathrm{d}t \\ \mathrm{d}c_\mathrm{D}/\mathrm{d}t \\ \mathrm{d}V/\mathrm{d}t \end{bmatrix} = \begin{bmatrix} c_\mathrm{A}^\mathrm{in} \dot{V}^\mathrm{in} + r_\mathrm{A}V \\ c_\mathrm{B}^\mathrm{in} \dot{V}^\mathrm{in} + r_\mathrm{B}V \\ c_\mathrm{C}^\mathrm{in} \dot{V}^\mathrm{in} + r_\mathrm{C}V \\ c_\mathrm{D}^\mathrm{in} \dot{V}^\mathrm{in} + r_\mathrm{D}V \\ \sum_{i \in (\text{A, B, C, D})}\left[(c_\mathrm{i}^\mathrm{in} \dot{V}^\mathrm{in} + r_\mathrm{i}V)\frac{M_i}{\rho_i}\right] \end{bmatrix}$ .

Rovnici múžeme zjednodušit, pokud si všimneme, že výraz v posledním řádku matice je roven 1 (viz stavová rovnice)

: $\begin{bmatrix} V & 0 & 0 & 0 & c_\mathrm{A} \\ 0 & V & 0 & 0 & c_\mathrm{B} \\ 0 & 0 & V & 0 & c_\mathrm{C} \\ 0 & 0 & 0 & V & c_\mathrm{D} \\ 0 & 0 & 0 & 0 & 1 \end{bmatrix} \begin{bmatrix} \mathrm{d}c_\mathrm{A}/\mathrm{d}t \\ \mathrm{d}c_\mathrm{B}/\mathrm{d}t \\ \mathrm{d}c_\mathrm{C}/\mathrm{d}t \\ \mathrm{d}c_\mathrm{D}/\mathrm{d}t \\ \mathrm{d}V/\mathrm{d}t \end{bmatrix} = \begin{bmatrix} c_\mathrm{A}^\mathrm{in} \dot{V}^\mathrm{in} + r_\mathrm{A}V \\ c_\mathrm{B}^\mathrm{in} \dot{V}^\mathrm{in} + r_\mathrm{B}V \\ c_\mathrm{C}^\mathrm{in} \dot{V}^\mathrm{in} + r_\mathrm{C}V \\ c_\mathrm{D}^\mathrm{in} \dot{V}^\mathrm{in} + r_\mathrm{D}V \\ \sum_{i \in (\text{A, B, C, D})}\left[(c_\mathrm{i}^\mathrm{in} \dot{V}^\mathrm{in} + r_\mathrm{i}V)\frac{M_i}{\rho_i}\right] \end{bmatrix}$ .

Vidíme, že matice M je regulární a že pomocí zpětné substituce již snadno odvodíme soustavu obyčejných diferenciálních rovnic. Jelikož jsme k transformaci DAE na ODE potřebovali jednou derivovat, lze tento příspěvek uzavřít s tím, že DAE index našeho modelu je roven 1.