Modelování 6: Vsádková destilace (DAE)

V další části se budeme věnovat řešení systému diferenciálně algebraických rovnic (DAE) v Matlabu. Jako příklad DAE systému může posloužit model vsádkové destilace, jehož řešení podrobně rozebereme.

Diferenciálně algebraické rovnice

DAE = ODE + AE

Na DAE systémy je možné pohlížet jako na kombinaci soustavy obyčejných diferenciálních rovnic (ODE) a algebraických rovnic (AE). V chemickém inženýrství se s DAE problémy setkáme v případech, kdy v modelu máme

vztahy vyjadřujících zákony zachování hmoty/energie (tj. bilance), v neustálených systémech jim odpovídají diferenciální rovnice, a zároveň
vztahy vyjadřující nějakou vaznou podmínku (například rovnováhu, stavovou rovnici, atd.), které jsou ve formě algebraických rovnic.

Obecný systém DAE je formálně možné vyjádřit implicitním zápisem

: $\displaystyle\mathbf{f}\left(t,\mathbf{y},\frac{\mathrm{d}\mathbf{y}}{\mathrm{d}t}\right)=\mathbf{0}.$

Jeho podmnožinou jsou problémy vyjádřitelné lineárně-implicitním zápisem

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

kde matice M se nazývá “mass matrix”. Pokud existuje inverzní matice M^-1 (tj. M není singulární), tak se vlastně jedná pouze o ODE

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

“Opravdový” DAE systém tedy vyžaduje singuární mass matrix. Některé DAE systémy je možné vyjádřit tzv. semi-implicitním
zápisem, kde závislé stavové proměnné y lze rozdělit na diferenciální y_d a algebraické y_a stavové proměnné

: $\begin{aligned}\displaystyle\frac{\mathrm{d}\mathbf{y}_d}{\mathrm{d}t} &= \mathbf{f}(t,\mathbf{y}_d,\mathbf{y}_a)\\ \mathbf{g}(t,\mathbf{y}_d,\mathbf{y}_a) &= \mathbf{0} \end{aligned}$

DAE-index

Obyčejné diferencíální rovnice (ODE) múžeme považovat za DAE s indexem 0. Index systému DAE určuje, jak obtížné bude řešení soustavy z hlediska numerického algoritmu. Například implicitní Euler je algoritmus, který má pro obyčejné diferenciální rovnice přesnost druhého řádu O(h₂). Jinými slovy, zkrácení délky kroku na polovinu vede ke řádově čtyřikrát menší odchylce od správné hodnoty. Každé zvýšení indexu vede u tohoto algoritmu ke ztrátě řádu přesnosti: DAE s indexem 1,2 a 3 mají O(h), 0(1) a O(1/h). To znamená, že DAE s řádem 2 a vyšším budou tímto algoritmem již neřešitelné, protože zkrácení délky kroku už nesníží velikost chyby a algorimtus nebude konvergovat k řešení.

Je důležité upozornit, že ani pokročilejší algoritmy používané v Matlabu ve funkci ode15i neumějí (dle dokumentace) řešit DAE s indexem 2 a vyšším.

Pro semi-implicitní zápis, je index je definován jako číslo, které udává kolikrát je nutné derivovat algebraické rovnice g(t,y_d,y_a), abychom DAE převedli na ODE

: $\begin{aligned}\displaystyle\frac{\mathrm{d}\mathbf{y}_d}{\mathrm{d}t} &= \mathbf{f}(t,\mathbf{y}_d,\mathbf{y}_a)\\ \frac{\mathrm{d}\mathbf{y}_a}{\mathrm{d}t} &= \mathbf{r}(t,\mathbf{y}_d,\mathbf{y}_a) \end{aligned}$

K příkladu postupu určení indexu pro konkrétní soustavy se vrátíme později.

DAE systémy s vyšším indexem jsou velmi obtížně řešitelné. Tyto překážky se dají překonat přeformulováním soustavy, které sníží její index. Tato operace ale může být pracná, hlavné kvůli nutnosti analytického derivování a přeodvozování. Poznamenejme, že Symbolic Math Toolbox verze Matlabu 2014 má k dispozici řadu funkcí pro analýzu a snížení indexu DAE systémů.

Model vsádkové destilace

Předpokládejme vsádkovou destilaci N-složkové směsi (N=3) probíhající za konstantního tlaku. Bilancovaným systémem bude kapalina ve varné nádobě. Potřebujeme tyto stavové proměnné: Látkové množství kapaliny $n_\mathrm{L}(t)$ a její složení vyjádřené pomocí látkových množství prvních N-1 složek $x_i(t)$ , molární tok $\dot{n}_\mathrm{V}(t)$ a složení odcházejících par $y_i(t)$ a teplotu kapaliny $T(t)$ , která je stejná jako teplota par. Celkem máme tedy 2N+1 stavových proměnných.

Parametry modelu, které musíme určit budou tlak v nádobě $p$ a tepelný tok dodávaný k udržení varu kapaliny $\dot Q$ . Uvažovaný děj se bude řídit molárními a entalpickou bilancí a dále rovnicemi popisujícími rovnováhu kapalina/pára. Pro jednoduchost opět využijeme Raoultův zákon, který je použitelný například pro směs hexan/heptan/oktan. Budeme mít tedy 7 rovnic pro 7 stavových proměnných

: $\begin{aligned} \frac{\mathrm{d}n_\mathrm{L}}{\mathrm{d}t} &= -\dot{n}_\mathrm{V}\\ \frac{\mathrm{d}}{\mathrm{d}t}\left(x_\mathrm{A} n_\mathrm{L}\right)&= -y_\mathrm{A}\dot{n}_\mathrm{V}\\ \frac{\mathrm{d}}{\mathrm{d}t}\left(x_\mathrm{B} n_\mathrm{L}\right)&=-y_\mathrm{B}\dot{n}_\mathrm{V}\\ \frac{\mathrm{d}}{\mathrm{d}t}\left(h_\mathrm{L} n_\mathrm{L}\right) &= \dot{Q} - h_\mathrm{V}\dot{n}_\mathrm{V} \\ 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 (1-y_\mathrm{A}-y_\mathrm{B}) &= p^\circ_\mathrm{C}(T) (1-x_\mathrm{A}-x_\mathrm{B}) \end{aligned}$

Bilance destilátu

Kromě kapaliny ve varné baňce (tj. zbytku) nás bude zajímat i složení předestilované kapaliny (tj. destilátu). Mohli bychom tedy do modelu přidat tři další stavové proměnné $n_\mathrm{D}(t)$ , $x_\mathrm{AD}(t)$ a $x_\mathrm{BD}(t)$ a model doplnit o molární bilance destilátu.

Když si ale uvědomíme, že předestilovaná kapalina už nijak nemůže ovlivnit děje ve varné baňce, a že složení a obsah destilátu lze snadno dopočítat, pokud známe množství suroviny na počátku

: $\begin{aligned} n_\mathrm{D}(t) &= n_\mathrm{L}(0) - n_\mathrm{L}(t) \\ x_\mathrm{AD}(t)n_\mathrm{D}(t) &= x_\mathrm{A}(0)n_\mathrm{L}(0) - x_\mathrm{A}(t)n_\mathrm{L}(t) \\ x_\mathrm{BD}(t)n_\mathrm{D}(t) &= x_\mathrm{B}(0)n_\mathrm{L}(0) - x_\mathrm{B}(t)n_\mathrm{L}(t) \end{aligned}$

můžeme, pro zjednodušení, bilance destilátu z modelu vynechat a dopočítat až po dokončení integrace v rámci tzv. post-processingu výsledných dat.

Entalpická bilance

Před implementací modelu ve formě kódu v MATLABU je potřeba upravit entalpickou bilanci do použitelnější formy. Rozepsání obou členů u derivace součinu funkcí dává

: $\displaystyle h_\mathrm{L}\frac{\mathrm{d}n_\mathrm{L}}{\mathrm{d}t}+ n_\mathrm{L}\frac{\mathrm{d}h_\mathrm{L}}{\mathrm{d}t}=\dot{Q}-h_\mathrm{V}\dot{n}_\mathrm{V}$

a protože z celkové bilance $\frac{\mathrm{d}n_\mathrm{L}}{\mathrm{d}t}=-\dot{n}_\mathrm{V}$ je možné psát

: $\displaystyle n_\mathrm{L}\frac{\mathrm{d}h_\mathrm{L}}{\mathrm{d}t}=\dot{Q}-(h_\mathrm{V}-h_\mathrm{L})\dot{n}_\mathrm{V}$

Teď je potřeba vyjádřit molární entalpie kapaliny a páry, která bude záviset na jejich složení, teplotě a skupenství

: $\begin{aligned} h_\mathrm{L} &= h_\mathrm{l}(\mathbf{x},T)\\ h_\mathrm{V} &= h_\mathrm{g}(\mathbf{y},T) = h_\mathrm{l}(\mathbf{y},T) + \Delta_{lg} h(\mathbf{y},T) \end{aligned}$

Derivace molární entalpie lze zapsat pomocí derivace složené funkce

: $\displaystyle \frac{\mathrm{d}h_\mathrm{L}}{\mathrm{d}t} = \frac{\partial h}{\partial \mathbf{x}}\frac{\mathrm{d}\mathbf{x}}{\mathrm{d}t} + \frac{\partial h}{\partial T}\frac{\mathrm{d}T}{\mathrm{d}t} = 0 + C_p(\mathbf{x},T) \frac{\mathrm{d}T}{\mathrm{d}t}$

kde jsme předpokládali nulové směšovací entalpie (nulový první člen) a použili definici tepelné kapacity (druhý člen).

Pokud zvolíme jako referenční teplotu $T_\mathrm{ref}=T$ , referenčni skupenství kapalné, předpokládáme nulové směšovací entalpie, potom

: $\displaystyle h_\mathrm{l}(\mathbf{x},T) = h_\mathrm{l}(\mathbf{y},T) = 0$

Po dosazení dostaneme výslednou podobu entalpické bilance

: $\displaystyle n_\mathrm{L}C_p(\mathbf{x},T)\frac{\mathrm{d}T}{\mathrm{d}t} = \dot{Q} - \Delta_{lg} h(\mathbf{y},T) \dot{n}_\mathrm{V}$

kde fyzikální vlastnosti směsi je možné vyjádřit váženým průměrem fyzikálních vlastností jednotlivých složek

: $\displaystyle n_\mathrm{L} \left(\sum_{i \in \mathrm{A,B,C}} x_i C_{p,i}(T)\right)\frac{\mathrm{d}T}{\mathrm{d}t} = \dot{Q} - \left( \sum_{i \in \mathrm{A,B,C}} y_i \Delta_{lg} h_i(T)\right) \dot{n}_\mathrm{V}$

Výpočet fyzikálních vlastností

Předtím, než začneme psát vlastní program pro model vsádkové destlilace je nutné si připravit funkce pro výpočet fyzikálních vlastností.

Pro výpočet tenze par jsme již dříve vytvořili funkce ptension a ptension_coef.

Také pro výpočet výparné entalpie již máme funkce dhevap a dhevap_coef z minulé části.

Analogickým způsobem si připravíme i funkce pro výpočet tepelné kapacity cpliq a cpliq2

function cp = cpliq(t,C)
%CPLIQ Liquid heat capacities (VII-3)
%   Calculate the heat capacity of liquid phase mixture
%   Using "VII-3" equation from "Holecek: Chemicko-inzenyrske tabulky"
%
%   CP = CPLIQ(T,C)
%
%   T   - temperature in Celsius
%   C   - [TK, A1, A2, A3, A4] constants from thermodynamic table
%   CP  - heat capacity in kJ/mol.K
 
ta   = t+273.15;   % temperature 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
 
cp = A1.*log(1-ta./TK) + A2./(1-ta./TK) + A3 + A4.*ta./TK;
cp = cp * Rgas / 1000; % kJ/(mol.K)
end

function cp = cpliq2(t,C)
%CPLIQ2 Liquid heat capacities (VII-2)
%   Calculate the heat capacity of liquid phase mixture
%   Using "VII-2" equation from "Holecek: Chemicko-inzenyrske tabulky"
%
%   CP = CPLIQ2(T,C)
%
%   T  - temperature in Celsius
%   C  - [TK, A1, A2, A3, A4] constants from thermodynamic table
%   CP - heat capacity in (kJ/mol.K)
 
ta   = t+273.15;   % temperature 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
 
cp = A1 + A2.*(ta./100) + A3.*(ta./100).^2 + A4.*(ta./100).^3;
cp = cp * Rgas / 1000; % kJ/(mol.K)
end

Pro načtení koeficientů empirických rovnic lze využít již vytvořenou funkci dhheat_coef.

Implementace DAE modelu

Pro řešení DAE systému použijeme Matlabovskou funkci ode15i. Začneme vytvořením funkce pro výpočet residuálu $\mathbf{F}(t,\mathbf{y},\mathbf{y}')$ , kterou pojmenujeme batdes_res. Vstupní proměnné této funkce jsou: nezávislá proměnná (čas), vektor stavových proměnných a vektor derivací stavových proměnných. Funkce by měla vrátit sloupcový vektor residuálů. Náš model je autonomní, tj. tvar rovnic se nemění v čase a proto první vstupní proměnná není ve funkci použitá. Zároveň jsme do funkce přidali proměnnou pars, ve které předáváme parametry modelu.

function res = res_batdes(~,y,yp,pars)
%RES_BATDES Calculate residual for batch distillation model
%   Batch distillation model
%      RES = RES_BATDES(~,Y,YP,PARS)
%   Y, YP - vectors of state variables and their derivatives
%      [nL, xA, xB, nVdot, yA, yB, t]
%   PARS  - model parameters
%      [p, Qdot]
 
% unpack parameters
p    = pars(1); % pressure [kPa]
Qdot = pars(2); % heat input [kW]
 
% unpack state variables (differential)
nL = y(1);    % liquid amount [mol]
xA = y(2);
xB = y(3);
xC = 1-xA-xB; % liquid mole fractions [mol/mol]
t  = y(7);    % temperature [C]
 
dnLdt = yp(1);
dxAdt = yp(2);
dxBdt = yp(3);
dtdt  = yp(7); % derivatives of differential state variables
 
% unpack state variables (algebraic)
nVdot = y(4); % vapour flux [mol/s]
yA    = y(5);
yB    = y(6);
yC = 1-yA-yB; % vapour mole fractions [mol/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]
 
% molar heat capacity of liquid phase
coef = dhheat_coef();
cpA = cpliq (t,coef(1,:));
cpB = cpliq2(t,coef(2,:));
cpC = cpliq (t,coef(3,:));
cpABC = xA*cpA + xB*cpB + xC*cpC; % kJ/(mol.K)
%cpABC = xA*0.221+xB*0.257+xC*0.287;
 
% molar evaporation enthalpy
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
%evapABC = yA*26.618+yB*31.654+yC*36.353;
 
% calculate residuals
% mole balances
res(1) = dnLdt + nVdot;
res(2) = xA*dnLdt + dxAdt*nL + yA*nVdot;
res(3) = xB*dnLdt + dxBdt*nL + yB*nVdot;
% 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) = nL*cpABC*dtdt - Qdot + evapABC*nVdot; % kJ/s
%res(7) = nVdot - 0.03;
 
res = res'; % make residual column vector
end

Kromě nové proměnné (vektor derivací stavových proměnných) se struktura kódu a postup řešení DAE zatím podobá řešení ODE nebo AE. Teď se ale dostáváme k problému, se kterým jsme u AE nebo ODE nesetkali, a tím je určení konzistentních počátečních podmínek.

Konzistentní počáteční podmínky

U soustavy obyčejných diferenciálních rovnic (ODE) můžeme (a musíme) zvolit počáteční hodnotu pro všechny stavové proměnné v soustavě. U DAE systému je situace složitější: počet stupňů volnosti je menší než počet stavových proměnných. To znamená, že volbou některých stavových proměnných jsou již určeny hodnoty ostatních stavových proměnných. Jinými slovy, počáteční stavové proměnné nelze volit libovolně, ale jejich hodnoty musí být konzistentní. U DAE s indexem 1 stačí, aby byla splněna podmínka, že

: $\mathbf{F}(t_0,\mathbf{y}_0,\mathbf{y}'_0) = \mathbf{0}$

Musíme najít takové vektory y0 a yp0, abychom po jejich použití ve funkci residuálů dostali nulový vektor. U DAE vyších řádů než 1 se mohou vyskytnout další omezení, takže výše uvedená podmínka nemusí být postačující.

K nalezení počátečních podmínek pomáha Matlabovská funkce decic. Na vstupu zadáme odhad stavových proměnných a odhad jejich derivací. Funkce se pokusí upravit zadané hodnoty tak, aby na výstupu již byla splněna výše uvedená podmínka. Některé stavové proměnné a některé derivace lze označit jako zafixované, tj. funkce decic je ponechá beze změny. Protože máme 2N proměnných (N stavových proměnných a N derivací stavových proměnných) a funkce residuálů obsahuje N rovnic, můžeme zafixovat maximálně N proměnných. Některé proměnné není možné fixovat a je nutné je nechat volné i v případech, když známe jejich hodnotu. Jak je to se stavovými proměnnými v příkladu vsádkové destilace?

Složení a množství kapalné fáze (n_L, x_A, x_B): Tyto proměnné je možné zvolit libovolně, protože množství a složení vsádky na počátku musí být zadáno. Tyto proměnné je tedy nutné zafixovat. Derivace těchto proměnných necháme volné, aby se mohly spočítat.
Složení parní fáze (y_A, y_B): Složení parní fáze je určeno již zadáním složení kapalné fáze. Tyto proměnné tedy nelze zafixovat, můžeme nastřelit nějaké hodnoty a nechat je dopočítat. Derivace těchto proměnných se ve funkci residuálů neobjevují a není je možné nijak vypočítat. Můžeme je proto nastavit nulové a zafixovat. To samé platí pro derivaci toku páry.
Teplota a tok parní fáze: Teplota je podobně jako složení páry určena již složením kapaliny a musí odpovídat teplotě varu kapaliny. Teplotu tedy nelze zafixovat. Zbývá vyřešit, co s derivací teploty a tokem parní fáze. Zatím jsme zafixovali 6 hodnot ze 7 maximálně zafixovatelných. Takže není možné určit obě zbývající. Zkusme tedy zafixovat derivaci teploty na nule. Také je možné nechat obě proměnné nefixované

Tady je část kódu, ve které hledáme počáteční podmínky. Jedničky ve vektorech fix znamenají, že daná proměnná je fixována.

>> y = [100, 0.33, 0.33, 0.1, 0.5, 0.4, 50]; >> yfix = [ 1, 1, 1, 0, 0, 0, 0];

Nastavili jsme počáteční složení kapaliny (33%-hexan / 33%-heptan / 34%-oktan) a 100 mol vsádky. Odhadli jsme tok parní fáze 0.1 mol/s, složení páry a teplotu 50°C.

>> yp = [0, 0, 0, 0, 0, 0, 0]; >> ypfix = [0, 0, 0, 1, 1, 1, 1];

Zde jsme zafixovali derivace stavových proměnných, jejichž derivace se v modelu nevyskytují a (spekulativně) i derivaci teploty. Zde je vzor volání funkce decic.

>> [y0,yp0] = decic(@fce, t0, y, yfix, yp, ypfix, options) >> fce(t0,y0,yp0)

Výsledkem jsou vypočítané počáteční podmínky a ověření, že residuál je na počátku malé číslo.

Integrace

Vlastní integrace se provádí pomocí funkce ode15i. Postup volání je stejný jako u ostatních funkcí z ODE-knihovny pouze s rozdílem, že se navíc zadává i vektor derivací stavových proměnných.

Po dokončení integrace zbývá dopočítat složení destilátu a dokreslit příslušné grafy. Bez dalšího komentáře uvedeme hlavní skript.

% batch distillation (main script)
%
% model parameters
p    = 101.325; % kPa
Qdot = 1;       % kW
pars = [p, Qdot];
 
% initial conditions (known and being fixed)
time0 = 0.; timef = 3300; % s
nL0 = 100;  % mol
xA0 = 0.33; xB0 = 0.33; xC0 = 1-xA0-xB0;
% initial conditions (guessed and being calculated)
yAg = 0.5; yBg = 0.4;
nVdotg = 0.1; % mol/s
tg = 50; % C
 
% calculate consistent initial conditions
options = odeset('RelTol',1e-8);
yg     = [nL0, xA0, xB0, nVdotg, yAg, yBg, tg];
ygfix  = [1,   1,   1,   0,      0,   0,   0 ];
 
ypg    = [0.0, 0.0, 0.0, 0.0,    0.0, 0.0, 0.];
ypgfix = [0,   0,   0,   1,      1,   1,   1 ];
% note1: Derivatives of algebraic states should be fixed at ZERO value.
% note2: It is recommended to fix the initial temperature derivation
%        to zero. Otherwise "decic" wont find good initial condition.
 
[y0,yp0] = decic(@(t,y,yp) res_batdes(t,y,yp,pars), ...
                 time0, yg, ygfix, ypg, ypgfix, options)
res_batdes(time0,y0,yp0,pars); % check calculated initial conditions
 
% integrate DAE
[tt,yy] = ode15i(@(t,y,yp) res_batdes(t,y,yp,pars), ...
                 [time0, timef], y0, yp0, options)
 
% post-processing (calculate destilate composition)
nL = yy(:,1);
xA = yy(:,2); xB = yy(:,3); xC = 1-xA-xB;
yA = yy(:,5); yB = yy(:,6); yC = 1-yA-yB;
nVdot = yy(:,4);
t = yy(:,7);
 
nD = nL0 - nL;                  % destilate amount
yAD = (xA0.*nL0 - xA.*nL)./nD;
yBD = (xB0.*nL0 - xB.*nL)./nD;
yCD = 1-yAD-yBD;                % mole fractions in the destilate
 
% plot graphs
subplot(2,2,1) % composition
plot(tt, xA, 'r-',  tt, xB,  'g-',  tt, xC,  'b-', ...
     tt, yA, 'ro',  tt, yB,  'go',  tt, yC,  'bo', ...
     tt, yAD,'r--', tt, yBD, 'g--', tt, yCD, 'b--')
 xlabel('time / s')
 ylabel('mol / mol')
 
 subplot(2,2,2) % liquid amount
 plot(tt, nL, 'k-d')
 xlabel('time / s')
 ylabel('nL / mol')
 
 subplot(2,2,3) % temperature
 plot(tt, t, 'ro')
 xlabel('time / s')
 ylabel('temperature / C')
 
 subplot(2,2,4) % vapour flux
 plot(tt, nVdot, 'b+')
 xlabel('time / s')
 ylabel('nVdot / mol/s')

% batch distillation (main script) % % model parameters p = 101.325; % kPa Qdot = 1; % kW pars = [p, Qdot]; % initial conditions (known and being fixed) time0 = 0.; timef = 3300; % s nL0 = 100; % mol xA0 = 0.33; xB0 = 0.33; xC0 = 1-xA0-xB0; % initial conditions (guessed and being calculated) yAg = 0.5; yBg = 0.4; nVdotg = 0.1; % mol/s tg = 50; % C % calculate consistent initial conditions options = odeset('RelTol',1e-8); yg = [nL0, xA0, xB0, nVdotg, yAg, yBg, tg]; ygfix = [1, 1, 1, 0, 0, 0, 0 ]; ypg = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.]; ypgfix = [0, 0, 0, 1, 1, 1, 1 ]; % note1: Derivatives of algebraic states should be fixed at ZERO value. % note2: It is recommended to fix the initial temperature derivation % to zero. Otherwise "decic" wont find good initial condition. [y0,yp0] = decic(@(t,y,yp) res_batdes(t,y,yp,pars), ... time0, yg, ygfix, ypg, ypgfix, options) res_batdes(time0,y0,yp0,pars); % check calculated initial conditions % integrate DAE [tt,yy] = ode15i(@(t,y,yp) res_batdes(t,y,yp,pars), ... [time0, timef], y0, yp0, options) % post-processing (calculate destilate composition) nL = yy(:,1); xA = yy(:,2); xB = yy(:,3); xC = 1-xA-xB; yA = yy(:,5); yB = yy(:,6); yC = 1-yA-yB; nVdot = yy(:,4); t = yy(:,7); nD = nL0 - nL; % destilate amount yAD = (xA0.*nL0 - xA.*nL)./nD; yBD = (xB0.*nL0 - xB.*nL)./nD; yCD = 1-yAD-yBD; % mole fractions in the destilate % plot graphs subplot(2,2,1) % composition plot(tt, xA, 'r-', tt, xB, 'g-', tt, xC, 'b-', ... tt, yA, 'ro', tt, yB, 'go', tt, yC, 'bo', ... tt, yAD,'r--', tt, yBD, 'g--', tt, yCD, 'b--') xlabel('time / s') ylabel('mol / mol') subplot(2,2,2) % liquid amount plot(tt, nL, 'k-d') xlabel('time / s') ylabel('nL / mol') subplot(2,2,3) % temperature plot(tt, t, 'ro') xlabel('time / s') ylabel('temperature / C') subplot(2,2,4) % vapour flux plot(tt, nVdot, 'b+') xlabel('time / s') ylabel('nVdot / mol/s')

Pages: 12