Modelování 4: Mžiková destilace

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

  1. Vyzkoušejte měnit teplotu nástřiku a tlak v zařízení a sledujte podíl odpařeného nástřiku nV/nF. 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.

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

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

  4. 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 nebo v tomto dokumentu. V našem případě by se počet rovnic v soustavě snížil na dvě: Rachford-Riceho a entalpickou bilanci.
Print This Post Print This Post
Pages: Page 1, Page 2

Published by

Zdenek Grof

I am administrator of this site.

One thought on “Modelování 4: Mžiková destilace”

  1. Pingback: Modelování 6: Vsádková destilace (DAE) – Scorpion Blog

Leave a Reply

Your email address will not be published. Required fields are marked *