Scorpion Blog » Modelování 4: Mžiková destilace

V této části podrobně rozebereme modelování mžikové destilace. Matematický model systému, který odvodíme, vede ke soustavě algebraických rovnic. Podrobněji se podíváme na entalpickou bilanci. Pro rovnováhu kapalina/pára opět použijeme Raoultův zákon.

Matematický model

Stavové proměnné

Analýzu problému uděláme nejdříve pro třísložkovou směs, a v závěru ji zobecníme na systém N složek. Jaké proměnné potřebujeme k tomu, abychom u mžikové destilace určili stav systému?

Množství a složení proudu nástřiku (3 proměnné),
Množství a složení proudu zbytku/kapalné fáze (3 proměnné),
Množství a složení proudu destilátu/parní fáze (3 proměnné),
Teplotu nástřiku,
Teplotu a tlak v dělícím bubnu (tj. za redukčním ventilem).

Zanedbáme vliv tlaku na enthalpii a proto nebudeme potřebovat znát tlak nástřiku. Pouze předpokládáme, že tlak bude takový, aby nástřik byl při své teplotě kapalina. Konkrétní tři proměnné určující množství a složení proudu lze vybrat z několika možností. Jednou z nich jsou molární, případně i hmotnostní toky jednotlivých složek. My vybereme celkový molární tok a molární zlomky prvních dvou složek. Molární zlomek třetí složky lze snadno dopočítat a nebudeme s ním počítat jako se stavovou proměnnou. Teplotu obou výstupních proudů předpokládáme stejnou jako je teplota v dělícím bubnu.

Celkem máme tedy 12 stavových proměnných pro třísložkovou směs. Tyto proměnné ale nejsou nezávislé, protože mezi nimi platí různé vztahy

Bilanční rovnice reprezentující zákon zachování hmoty (3 nezávislé rovnice)
Rovnice rovnováhy kapalina/pára reprezentující podmínku, že chemický potenciál každé složky bude stejný v obou fázích (3 rovnice)
Enthalpická bilance reprezentující zákon zachování energie odvozená z 1. věty termodynamické (1 rovnice)

Jestliže je 12 stavových proměnných vázáno 7 rovnicemi, zůstává 5 volných stavových proměnných. Těchto 5 stavových proměnných je potřeba zvolit a stávají se z nich tzv. parametry modelu. Ze zbývajících 7 proměnných se stávají tzv. neznámé stavové proměnné, které určíme řešením soustavy 7 algebraických rovnic.

Než budeme pokračovat dále, zobecníme závěry na N-složkovou směs:

3N+3: Stavových proměnných celkem
2N+1: Neznámých stavových proměnných vázaných stejným počtem rovnic
N+2: Volných stavových proměnných, které je potřeba zvolit jako parametry modelu. Jinými slovy: počet stupňů volnosti.

Výběr parametrů

Pro třísložkovou směs jsme odvodili, že ze seznamu: n_F, z_A, z_B, n_L, x_A, x_B, n_V, y_A, y_B, T_F, T a p, potřebujeme určit 5 parametrů. Můžeme ale vybrat těchto pět proměnných zcela libovolně? Omezení samozřejmě existují. Například v seznamu parametrů nemohou být zároveň n_F, n_L a n_V. Pokud určíme dvě proměnné, pak určení třetí nepřinese novou informaci, protože všechny tři proměnné musí splňovat celkovou molární bilanci. Fyzikální intuice pomůže odůvodnit podobné zakázané kombinace. Z matematického úhlu pohledu je potřeba si uvědomit, že některé rovnice nezávisí na všech potenciálních neznámých. V celkové bilanci se vyskytují pouze molární toky tří proudů a residuál této rovnice je funkcí jen těchto tří proměnných. Pokud ze všech tří uděláme parametry, potom už je celková bilance zbytečná a musíme vzít v úvahu, že se tím sníží počet použitelných rovnic v soustavé.

V inženýrství existují dva typy úloh: simulační a návrhové. V obou případech považujme množství a složení nástřiku za známé (n_F, z_A, z_B) a tak zbýva určit dva parametry. Simulační úloha má odpovědět na otázku “co bude, když?”, tj. zadáme podmínky v zařízení a počítáme, jaké bude složení a teplota výstupních proudů. Parametry, kterými lze řídit mžikovou destilaci bude tlak v zařízení a teplota nástřiku (T_F, p). U návrhového výpočtu hledáme takové provozní podmínky, aby bylo dosaženo předepsaného výsledku. Parametry určující úlohu mohou být například předepsaný podíl odpařené kapaliny, výtežek konkrétní složky, molární zlomek určité složky v destilátu a podobné. Návrhový výpočet je obtížnější než simulační, protože je velice snadné zadat takové požadavky, které nelze splnit (soustava nemá řešení pro konkrétní hodnoty parametrů). Pro začátek proto vyřešíme simulační úlohu, kde budou zadány tyto parametry

množství a složení nástřiku (n_F, z_A, z_B)
teplota nástřiku T_F
tlak v zařízení p

Zbývající proměnné (n_L, x_A, x_B, n_V, y_A, y_B a T) budou neznámé.

Molární bilance a rovnováha

Vyjádření rovnováhy kapalina/pára jsme diskutovali již dříve. Opět budeme předpokládat platnost Raoultova zákona

$\begin{aligned}p\, y_\mathrm{A} &= p^\circ_\mathrm{A}(T) x_\mathrm{A}\\p\, y_\mathrm{B} &= p^\circ_\mathrm{B}(T) x_\mathrm{B}\\p\, y_\mathrm{C} &= p^\circ_\mathrm{C}(T) x_\mathrm{C}\end{aligned}$

U bilančních rovnic použijeme bilance tří složek

$\begin{aligned}\dot n_\mathrm{F}\, z_\mathrm{A} &= \dot n_\mathrm{L}\, x_\mathrm{A}+\dot n_\mathrm{V}\, y_\mathrm{A}\\\dot n_\mathrm{F}\, z_\mathrm{B} &= \dot n_\mathrm{L}\, x_\mathrm{B}+\dot n_\mathrm{V}\, y_\mathrm{B}\\\dot n_\mathrm{F}\, z_\mathrm{C} &= \dot n_\mathrm{L}\, x_\mathrm{C}+\dot n_\mathrm{V}\, y_\mathrm{C}\end{aligned}$

Entalpická bilance

Budeme se podrobněji věnovat odvození entalpické bilance. Předpokládáme nezávislost hodnoty entalpie na tlaku (což je striktně vzato možné pouze u ideálních plynů nebo nestačitelných kapalin) a tak molární entalpie h proudů závisí pouze na složení, teplotě a skupenství

$\dot n_\mathrm{F} h(\mathbf{z}, T_\mathrm{F}, l) = \dot n_\mathrm{L} h(\mathbf{x}, T, l) + \dot n_\mathrm{V} h(\mathbf{y}, T, g)$

Zvolíme-li referenční teplotu T_ref a referenční skupenství kapalné, lze molární entalpie vyjádřit

$\displaystyle\dot n_\mathrm{F} \int_{T_\mathrm{ref}}^{T_\mathrm{F}} C_p(\mathbf{z},T)\mathrm{d}T = \dot n_\mathrm{L} \int_{T_\mathrm{ref}}^{T} C_p(\mathbf{x},T)\mathrm{d}T + \dot n_\mathrm{V} \left[ \int_{T_\mathrm{ref}}^{T} C_p(\mathbf{y},T)\mathrm{d}T + \Delta_\mathrm{lg}h(\mathbf{y},T)\right]$

Referenční teplotu můžeme volit libovolně. Když zvolíme referenční teplotu stejnou jako teplotu v zařízení (T_ref=T), entalpická bilance se zjednoduší

$\displaystyle\dot n_\mathrm{F} \int_{T}^{T_\mathrm{F}} C_p(\mathbf{z},T)\mathrm{d}T = \dot n_\mathrm{V} \Delta_\mathrm{lg}h(\mathbf{y},T)$

Entalpická bilance v tomto tvaru ilustruje, že energie potřebná na odpaření kapaliny se rovná energii získané ochlazením nástřiku z teploty nástřiku na teplotu v zařízení. Tepelnou kapacitu a výparnou entalpii směsi můžeme vyjádřit jako vážený průměr tepelných kapacit a výparných entalpií jednotlivých složek

$\displaystyle\dot n_\mathrm{F} \sum_{i \in\{A,B,C\}} \left( z_i \int_{T}^{T_\mathrm{F}} C_{p,i}(T)\mathrm{d}T \right) = \dot n_\mathrm{V} \sum_{i \in\{A,B,C\}} y_i \Delta_\mathrm{lg}h_i(T)$

Budeme potřebovat funkce počítající rozdíl entalpií příslušné složky při změně teploty $\int_{T}^{T_\mathrm{F}} C_{p,i}(T)\mathrm{d}T$ a změně skupenství $\Delta_\mathrm{lg}h_i(T)$ . Závislosti tepelné kapacity i výparné entalpie na teplotě pro různé složky jsou tabelovány ve formě empirických vztahů, které lze nalézt například v chemicko-inženýrských tabulkách.

Implementace modelu

Model mžikové destilace implementujeme pro třísložkovou směs hexan/heptan/oktan. Připomeňme si, že z minulých částí máme funkce počítající závislost tenze par na teplotě. Kvůli entalpické bilanci potřebujeme vytvořit další funkce pro výpočet tepelných kapacit a výparných entalpií.

Výpočet integrálu tepelných kapacit

V chemicko-inženýrských tabulkách najdeme empirické vztahy pro výpočet tepelné kapacity při různých teplotách

$\displaystyle\frac{C_{p,i}}{R}=A_1 \ln(1-T/T_k) + \frac{A_2}{1-T/T_k}+A_3+A_4 T/T_k$

Integrál tepelné kapacity je možné vyřešit analyticky bez větších potíží

$\begin{aligned}\displaystyle\int_{T_a}^{T_b} \frac{C_{p,i}}{R}\mathrm{d}T =& \left[A_1(T-T_k)(\ln(1-T/T_k)-1)\right]_{T_a}^{T_b}\\ +& \left[-A_2 T_k \ln(T_k-T)\right]_{T_a}^{T_b} \\ +& \left[A_3 T\right]_{T_a}^{T_b} +\left[A_4/(2T_k) T^2\right]_{T_a}^{T_b}\end{aligned}$

Výpočet integrálu zapíšeme přímo do funkce dhheat.m

function dh = dhheat(t1,t2,C)
%DHHEAT Enthalpy difference due to temperature change (VII-3)
%   Calculate the enthalpy difference of liquid phase component as
%   temperature goes from T1 to T2
%   Using "VII-3" equation from "Holecek: Chemicko-inzenyrske tabulky"
%
%   DH = DHHEAT(T1,T2,C)
%
%   T1, T2 - temperatures in Celsius
%   C      - [TK, A1, A2, A3, A4] constants from thermodynamic table
%   DH     - enthalpy difference in kJ/mol

a = t1+273.15;
b = t2+273.15;     % temperature range in K
Rgas = 8.3144621;  % ideal gas constant J/mol.K

TK = C(1); % unpack constants
A1 = C(2);
A2 = C(3);
A3 = C(4);
A4 = C(5);

% Equation VII-3
% Cp/R = A1*ln(1-T/TK) + A2/(1-T/TK) + A3 + A4*T/TK
% Int{A1*ln(1-T/TK) dT} = [ A1*(T-TK)*(ln(1-T/TK)-1)]
% Int{A2/(1-T/TK) dT} = [ -A2*TK*ln(TK-T)]
% Int{A3 dT} = [A3*T]
% Int{A4*T/TK dT} = [A4/(2*TK)*T^2]

dh = A1*(b-TK)*(log(1-b/TK)-1) - ...
     A1*(a-TK)*(log(1-a/TK)-1);

dh = dh + (-A2*TK*log(TK-b)) - ...
        (-A2*TK*log(TK-a));

dh = dh + (A3*(b-a));
dh = dh + (A4/(2*TK)*(b^2-a^2));
dh = dh * Rgas / 1000; % kJ/mol
end

Protože je heptan tabelován pomocí jiné rovnice než hexan a oktan, musíme bohužel vytvořit funkci zvlášť pro heptan, dhheat2.m

function dh = dhheat2(t1,t2,C)
%DHHEAT2 Enthalpy difference due to temperature change (VII-2)
%   Calculate the enthalpy difference of liquid phase component as system
%   goes from temperature T1 to temperature T2
%   Using "VII-2" equation from "Holecek: Chemicko-inzenyrske tabulky"
%
%   DH = DHHEAT2(T1,T2,C)
%
%   T1, T2 - temperatures in Celsius
%   C      - [TK, A1, A2, A3, A4] constants from thermodynamic table
%   DH     - enthalpy difference in kJ/mol

a = t1+273.15;
b = t2+273.15;     % temperature range in K
Rgas = 8.3144621;  % ideal gas constant J/mol.K

TK = C(1); % unpack constants
A1 = C(2);
A2 = C(3);
A3 = C(4);
A4 = C(5);

% Equation VII-2
% Cp/R = A1 + A2*(T/100) + A3*(T/100)^2 + A4*(T/100)^3
% Int{A1 dT} = [A1*T]
% Int(A2*T/100} = [A2/100 * T^2/2]
% Int(A3*(T/100)^2 = [A3/100^2 * T^3/3]
% Int(A4*(T/100)^3 = [A4/100^3 * T^4/4]

dh = A1*(b-a);
dh = dh + (A2/100   * (b^2-a^2)/2);
dh = dh + (A3/100^2 * (b^3-a^3)/3);
dh = dh + (A4/100^3 * (b^4-a^4)/4);
dh = dh * Rgas / 1000; % J/mol
end

Konstanty empirických rovnic je třeba dodat ve vektoru C. Hodnoty konstant budou uložené v “databází” dhheat_coef.m

function [CC] = dhheat_coef()
%DHHEAD_COEF Coefficients for DHHEAT or DHHEAT2 functions
%   Provide coefficients from thermodynamic tables
%   Rows of CC matrix contains data for individual components
%   Columns of CC matrix contain: TK, A1, A2, A3, A4

% literature: "Holecek, Chemickoinzenyrske tabulky, table VII-2"

% hexane (equation VII-3, validity range 180 K -- 443 K)
TK(1) = 507.90;
A1(1) = -3.44495;
A2(1) = 0.467981;
A3(1) = 15.6153;
A4(1) = 6.333982;

% heptane (equation VII-2, validity range 260 K -- 400 K)
TK(2) = 540.15;
A1(2) = 25.4134;
A2(2) = -4.11613;
A3(2) = 1.82098;
A4(2) = -0.0860461;

% octane (equation VII-3, validity range 222 K -- 426 K)
TK(3) = 568.95;
A1(3) = -5.81193;
A2(3) = 0.916902;
A3(3) = 19.5769;
A4(3) = 9.20996;

% pack coefficients into matrix CC
CC = [TK', A1', A2', A3', A4'];
end

Teď je nutné funkce otestovat, zda fungují správně. Tento krok nepodceňujte, protože případná chyba se najde snadněji teď než až při testování modelu mžikové destilace. Vypočítejme potřebné teplo k ohřátí jednoho molu hexanu/heptanu/oktanu z 50 na 150 °C.

>>c = dhheat_coef() >> dhheat(50,150,c(1,:)) >> dhheat2(50,150,c(2,:)) >> dhheat(50,150,c(3,:))

Výsledky jsou vypisovány v kJ/mol, a mělo by vyjít 22.26; 25.78 a 28.77. Připomeňme, že výraz c(1,:)je první řádek matice c.

Výpočet výparné entalpie

Programování výpočtu fyzikálních vlastností může být zdlouhavé, ale výpočet výparné entalpie, který nám zbývá, je naštěstí jednodušší než integrál tepelných kapacit.

$\Delta_{lg} h_i(T)=A (1-T/T_k)^\beta\exp(-\alpha T/T_k)$

Výparnou entalpii zjištěnou pomocí empirické rovnice zapíšeme ve funkci dhevap.m

function dh = dhevap(t,C)
%DHEVAP Enthalpy of evaporation
%   Calculate evaporation enthalpy
%   Using "IX-1" equation from "Holecek: Chemicko-inzenyrske tabulky"
%
%   DH = DHEVAP(T,C)
%
%   T  - temperatures in Celsius
%   C  - [TK, A, alf, bet] constants from thermodynamic table
%   DH - enthalpy difference in kJ/mol

TK  = C(1); % unpack constants
A   = C(2);
alf = C(3);
bet = C(4);
TR  = (t+273.15)/TK; % reduced temperature

% Equation IX-1 
% hLV = A * (1-TR)^bet * exp(-alf*TR)
dh = A .* ((1-TR).^bet) .* exp(-alf*TR) ;
end

a konstanty uložíme ve funkci dhevap_coef.m

function [CC] = dhevap_coef()
%DHEVAP_COEF Coefficients for DHEVAP function
%   Provide coefficients from thermodynamic tables
%   Rows of CC matrix contains data for individual components
%   Columns of CC matrix contain: TK, A, alpha, beta

% literature: "Holecek, Chemickoinzenyrske tabulky, table IX-1"

% hexane (validity range 298 K -- 444 K)
TK(1) = 507.40; %%507.90;
A(1) = 43.85;
alf(1) = -0.0390;
bet(1) = 0.3970;

% heptane (equation VII-2, validity range 260 K -- 400 K) 
TK(2) = 540.2; %%540.15;
A(2) = 53.66;
alf(2) = 0.2831;
bet(2) = 0.2831;

% octane (equation VII-3, validity range 222 K -- 426 K)
TK(3) = 568.8; %%568.95;
A(3) = 58.46;
alf(3) = 0.1834;
bet(3) = 0.3324;

% pack coefficients into matrix CC
CC = [TK', A', alf', bet'];
end

Otestujte, že výparné entalpie hexanu/heptanu/oktanu při 100°C jsou 26.62; 31.65 a 36.35 kJ/mol.

>> d = dhevap_coef() >> dhevap(100,d(1,:)) >> dhevap(100,d(2,:)) >> dhevap(100,d(3,:))

Funkce residuálů: Vlastní model mžikové destilace

Nyní je vše připraveno k tomu napsat matematické rovnice modelu. Neznámé stavové proměnné a parametry modelu dodáme ve dvou vektorech xvec a pars. Funkci residuálů pojmenujeme res_flash.m

function res = res_flash(xvec,pars)
%RES_FLASH Calculate residual for flash distillation model
%   Flash distillation model
%      RES = RES_FLASH(XVEC,PARS)
%   XVEC - uknown state variables
%      [nL, xA, xB, nV, yA, yB, t]
%   PARS - model parameters
%      [nF, zA, zB, tF, p]

% unpack parameters
nF = pars(1); % feed mole flux [mol/s]
zA = pars(2); % [mol/mol]
zB = pars(3);
zC = 1-zA-zB; % hexane/heptane/oktane mole fractions in feed
tF = pars(4); % feed temperature [C]
p  = pars(5); % pressure [kPa]

% unpack uknown variables
nL = xvec(1); % liquid phase flux [mol/s]
xA = xvec(2);
xB = xvec(3);
xC = 1-xA-xB; % liquid mole fractions [mol/mol]
nV = xvec(4); % vapour phase flux [mol/s]
yA = xvec(5);
yB = xvec(6);
yC = 1-yA-yB; % vapour mole fractions [mol/mol]
t  = xvec(7); % temperature [C]

% ethalpy change (cooling down the feed)
coef = dhheat_coef();
hA = dhheat (t,tF,coef(1,:));
hB = dhheat2(t,tF,coef(2,:));
hC = dhheat (t,tF,coef(3,:));
hABC = zA*hA + zB*hB + zC*hC; % kJ/mol

% enthalpy change (evaporation)
coefevap = dhevap_coef();
evapA = dhevap(t,coefevap(1,:));
evapB = dhevap(t,coefevap(2,:));
evapC = dhevap(t,coefevap(3,:));
evapABC = yA*evapA + yB*evapB + yC*evapC; % kJ/mol

% vapour pressures
[A,B,C] = ptension_coef();
pA = ptension(A(1),B(1),C(1),t);
pB = ptension(A(2),B(2),C(2),t);
pC = ptension(A(3),B(3),C(3),t); % kPa

% calculate residuals
% mole balances
res(1) = nF*zA - nL*xA - nV*yA;
res(2) = nF*zB - nL*xB - nV*yB;
res(3) = nF*zC - nL*xC - nV*yC;
% Raoults law equilibrium
res(4) = p*yA - pA*xA;
res(5) = p*yB - pB*xB;
res(6) = p*yC - pC*xC;
% enthalpy balance
res(7) = nF*hABC - nV*evapABC;

res = res'; % make residual column vector
end

Funkci se snažíme psát co nejpřehledněji. Na začátku zkopírujeme hodnoty parametrů a neznámých do proměnných s výstižnějším pojmenováním (“rozbalíme” parametry a neznámé). Potom dopočítáme na těchto stavových proměnných závislé vlastnosti (integrály tepelných kapacit, výparné entalpie, tenze par) a na závěr vypočítáme residuály soustavy rovnic.

Numerické řešení modelu

Teď je možné model vyřešit použitím funkcí pro řešení soustavy algebraických rovnic (fsolve nebo mmfsolve). Příkazy je možné zadat přímo v příkazovém řádku a po ověření, že vše funguje a soustava konverguje k řešení, dokončíme úlohu zapsáním sekvence potřebných příkazů do skriptu (hlavního programu)

% flash distillation (main script)
% 
% model parameters (N+2 values)
zA = 0.3;     % hexane mole fraction
zB = 0.5;     % heptane mole fraction
zC = 1-zA-zB; % octane mole fraction
nF = 10;      % feed flux mol/s
tF = 80;      % feed temperature [C]
p  = 0.2*101.325; % pressure [kPa]
pars = [nF, zA, zB, tF, p];

% guess values for uknowns (2N+1 values)
nL = 0.5; nV = 0.5; % [mol/s]
xA = 0.3; xB = 0.3; yA = 0.3; yB = 0.3; % [mol/mol]
t = 60; % [Celsius]
xvec0 = [nL, xA, xB, nV, yA, yB, t];

% solve set of uknowns
xsol = mmfsolve(@(x) res_flash(x,pars), xvec0);

% check solution converged
res = res_flash(xsol,pars)

% show results
xA = xsol(2); xB = xsol(3); xC = 1-xA-xB;
yA = xsol(5); yB = xsol(6); yC = 1-yA-yB;
nL = xsol(1); nV = xsol(4); t = xsol(7);
temperatures = [t, tF]
xyz = [xA,xB,xC; yA,yB,yC; zA,zB,zC]
a = nV/nF % evaporated ratio
%ABC_evaporated = [yA/zA*a, yB/zB*a, yC/zC*a]*100

Tak a to je celé, model mžikové destilace je dokončen.

Now the fun begins!

Pokud vás modelování baví, určitě neskončíte s řešením jednoho konkrétního zadání, ale budete si chtít vyzkoušet, jak se bude chování systému měnit v jiných podmínkách. Experimentování s fungujícím modelem Vám pomůže lépe pochopit fungování mžikové destilace. Tady je několik tipů a otázek k zamyšlení.

Vyzkoušejte měnit teplotu nástřiku a tlak v zařízení a sledujte podíl odpařeného nástřiku n_V/n_F. Zjistíte že k tomu, aby došlo k rozdělení na dvě fáze, lze teplotu nástřiku použít pouze z omezeného intervalu. Rozsah tohoto intervalu závisí na tlaku. Zkuste vysvětlit fyzikální význam těchto koncových teplot.
Zjistíte, že vytvořený model je dostatečně robustní a zkonverguje k řešení i v případech, kdy teplota nástřiku je mimo povolený interval. Řešení nemá fyzikální smysl, ale aspoň pomůže určit kterým směrem upravit teplotu nástřiku. Vysvětlete, proč model dává nefyzikální řešení a navrhněte jak by bylo možné získat fyzikální řešení i pro tyto případy.
Můžete se pokusit přepsat program na návrhový výpočet tak, aby našel vhodnou teplotu nástřiku pro zadaný podíl odpařeného nástřiku. Jednotlivé úpravy budou:
- Ve vektoru parametrů nahradí zadaný odpařený podíl teplotu nástřiku
```
aset = pars(4); % set required fraction evaporated (a=nV/nF)
```
- Teplota nástřiku bude novou neznámou stavovou proměnnou
```
tf = xvec(8); % feed temperature /C
```
- Protože tím v systému přibyla jedna neznámá, je nutné přidat novou rovnici (požadovaný odpařený podíl)
```
a = nV/nF; % calculate fraction evaporated using actual state
res(8)=a - aset; % fraction evaporated should be same as required
```
  Nejdříve určíme skutečný odpařený podíl použitím aktuálních hodnot stavových proměnných a potom zapíšeme podmínku, aby se tento podíl rovnal zadané hodnotě
Na mžikovou destilaci existuje postup podle Rachford-Rice metody, která umožní spojit řešení všech materiálových bilancí a rovnováh do jedné nelineární rovnice. Postup odvození je vysvětlen například v tomto videu z Youtube ^[1] nebo v tomto dokumentu ^[2]. V našem případě by se počet rovnic v soustavě snížil na dvě: Rachford-Riceho a entalpickou bilanci.