Zatím jsme při modelování potřebné fyzikální vlastnosti počítali z empirických vztahů, jako je např. Antoineho rovnice. V mnoha případech jsou ale fyzikální data k dispozici pouze ve formě tabulek. V této části ukážeme, jak je možné i s těmito daty pracovat a začlenit je do modelu.
Představme si, že pro závislost určité veličiny y na jiné veličině x máme tabulku, tj, dvojice hodnot [x, y]. Pokud potřebujeme napsat funkci, která pro zadané x najde odpovídající hodnotu y(x), je nutné proložit tabulované hodnoty nějakou funkcí a pomocí regrese určit příslušné koeficienty. Často používanými funkcemi jsou polynomy vyšších řádů, ale pro určité typy dat mohou být výhodnější i jiné tvary funkcí. Nevýhodou tohoto přístupu je to, že nefunguje univerzálně a správnost proložení je nutné pro každý případ zkontrolovat. Proto může být výhodnější prokládat data pomocí tzv. splinů.
Splines závislosti (křivky) jsou konstruovány tak, že splňují tyto podmínky:
- křivky procházejí všemi prokládanými body,
- křivky i jejich první derivace jsou hladké a spojité.
Těmto podmínkám vyhovují kubické funkce y = ai x3 + bi x2 + ci x + di, které mají pro každý interval mezi prokládanými body různé koeficienty.
Přeskočme nyní vlastní teorii i rovnice, ze kterých se počítají koeficienty splinů a ukažme si použití funkce spline v praxi.
Příklad tabulovaných dat může být rovnováha kapalina-pára pro systém methanol-voda uložená v souboru “methanol-voda.txt”
0.00 0.0000 100.0 0.0000 100.0 0.01 0.0763 98.1 0.0770 98.9 0.03 0.1948 94.8 0.1994 96.8 0.05 0.2822 92.3 0.2922 95.0 0.10 0.4257 87.5 0.4486 91.3 0.15 0.5144 84.3 0.5461 88.3 0.20 0.5767 81.8 0.6134 85.9 0.25 0.6247 79.9 0.6633 83.9 0.30 0.6640 78.3 0.7025 82.1 0.35 0.6979 76.8 0.7348 80.5 0.40 0.7282 75.5 0.7624 79.1 0.45 0.7558 74.4 0.7868 77.8 0.50 0.7817 73.3 0.8089 76.6 0.55 0.8062 72.2 0.8295 75.4 0.60 0.8296 71.3 0.8491 74.2 0.65 0.8523 70.3 0.8680 73.0 0.70 0.8744 69.4 0.8864 71.9 0.75 0.8960 68.5 0.9048 70.7 0.80 0.9173 67.7 0.9231 69.5 0.85 0.9382 66.9 0.9417 68.3 0.90 0.9590 66.1 0.9606 67.1 0.95 0.9796 65.3 0.9800 65.8 0.97 0.9878 65.0 0.9879 65.3 0.99 0.9959 64.7 0.9959 64.8 1.00 1.0000 64.5 1.0000 64.5
Připomeňme, že v prvních dvou sloupcích je molární zlomek metanolu v kapaliné a parní fázi a ve třetím sloupci je odpovídající teplota. Čtvrtý a pátý sloupec, které obsahují rovnovážná data vyjádřená hmotnostními zlomky, potřebovat nebudeme. Pro začátek:
>> dat = dlmread('methanol-voda.txt');
>> xtab = dat(:,1);
>> ytab = dat(:,2);
>> ttab = dat(:,3);
Příkaz dlmread načte čísla ze souboru do matice “dat”, které pro přehlednost zkopírujeme do příslušně pojmenovaných vektorů. Chceme-li nyní najít složení parní fáze yA odpovídající kapalné fázi o složení např. xA=0.50:
>> yA = spline(xtab, ytab, 0.50)
Tato funkce vytvoří spline procházející dvojicemi hodnot ve vektorech “xtab” a “ytab” a vrátí jeho hodnotu v bodě 0.5. Vidíme, že výsledek odpovídá údajům v tabulce.
Mžiková destilace metanol-voda
Rovnováhu zadanou tabulkou je možné využít při modelování mžikové destilace metanol-voda. Pro zjednodušení zpracujeme model bez entalpické bilance, tj. nebudeme počítat závislost mezi teplotou na vstupu a teplotou v zařízení. Zde je seznam všech stavových proměnných:
- složení a molární tok nástřiku (nF, zA)
- složení a molární tok destilátu (nV, yA)
- složení a molární tok zbytku (nL, xA)
- teplota v zařízení (t)
Protože uvažujeme dvousložkovou směs, tak molární zlomky druhé složky lze snadno spočítat a nebudeme s nimi počítat jako se stavovými proměnnými. Tlak v tomto případě musí být atmosférický a odpovídat tabelované rovnováze. Celkem stavových proměnných: 7
Vztahy mezi stavovými proměnnými:
- Molární bilance jednotlivých složek (2 rovnice)
- Vztahy pro rovnováhu kapalina/pára (2 rovnice)
Celkem rovnic: 4
Volných parametrů zůstává tedy 7-4=3. Zadáme tedy např. složení a průtok nástřiku a teplotu v zařízení. Ostatní stavové proměnné budou neznámé. Ukažme příklad funkce pro výpočet vektoru residuálů tohoto modelu:
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, nV, yA] % PARS - model parameters % [nF, zA, t] % unpack parameters nF = pars(1); % feed mole flux [mol/s] zA = pars(2); % [mol/mol] zB = 1-zA; % methanol/water mole fractions in feed t = pars(3); % temperature in the vessel [C] % unpack uknown variables nL = xvec(1); % liquid phase flux [mol/s] xA = xvec(2); xB = 1-xA; % liquid mole fractions [mol/mol] nV = xvec(3); % vapour phase flux [mol/s] yA = xvec(4); yB = 1-yA; % vapour mole fractions [mol/mol] % read vapour/liquid equilibrium from table data = dlmread('methanol-voda.txt'); xdat = data(:,1); % first column - xA ydat = data(:,2); % second column - yA tdat = data(:,3); % third column - t ym = spline(xdat,ydat,xA); tm = spline(xdat,tdat,xA); % calculate residuals % mole balances res(1) = nF*zA - nL*xA - nV*yA; res(2) = nF*zB - nL*xB - nV*yB; % vap/liq equilibrium (actual values should be consistent with tabulated data) res(3) = yA - ym; res(4) = t - tm; res = res'; % make residual column vector end |
Za povšimnutí stojí jakým rovnicemi jsme do modelu zahrnuli rovnováhu:
kde funkce f1 a f2 vyhledají pro konkrétní hodnotu xA v prvním sloupci tabulky příslušné hodnoty ve druhém a třetím sloupci. Má-li platit rovnováha kapalina/pára, pak hodnoty (xA, yA, t) musejí odpovídat hodnotám v tabulce a právě to zajišťují uvedené rovnice.
Úlohu dokončíme ukázkou hlavního programu (skriptu)
% flash distillation (main script) % methanol/water, without enthalpy balance % vapour/liquid equilibrium from tabulated data % % model parameters zA = 0.4; % methanol mole fraction zB = 1-zA; % water mole fraction nF = 10; % feed flux mol/s t = 80; % temperature in the vessel [C] % note: pressure is atmospheric (depends on table used) pars = [nF, zA, t]; % guess values for uknowns nL = 5; nV = 5; % [mol/s] xA = 0.3; yA = 0.6; % [mol/mol] xvec0 = [nL, xA, nV, yA]; % 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 = 1-xA; yA = xsol(4); yB = 1-yA; nL = xsol(1); nV = xsol(3); xyz = [xA,xB; yA,yB; zA,zB] a = nV/nF % evaporated ratio |
A to je všechno. Jako cvičení se můžete pokusit doplnit do modelu entalpickou bilanci. Spliny se dají také použít např. při výpočtech extrakce.
Uložení spline a jeho derivace
Někdo by mohl namítnout, že námi navržené použití splinu v modelu je neefektivní, protože při každém volání pravých stran se musí:
- otevřít a načíst soubor z disku,
- vypočítat koeficienty splinu,
- určit hodnotu splinu v daném bodě.
Přitom by stačilo první dva body provést pouze jednou: koeficienty splinu budou pořád stejné. Bude mít pravdu. Nicméně naše verze je efektivní z hlediska jednoduchosti: není nutné řešit způsob uložení a předávání koeficinetů na příslušné místo v programu. Efektivita z hlediska výpočetního času nás v tomto případě netrápí, protže zdržení při dnešních rychlostech počítače nepoznáme. Pokud byste narazili na problém, kde začíná být výpočetní čas kritický, tak je možné použití splinu rozdělit na dva kroky
>> pp = spline(xdat,ydat)
>> y = ppval(pp,x)
První příkaz vypočítá hodnoty koeficientů a uloží je do struktury
pp. Později, kde je to zapotřebí, se použije příkaz ppval, který už jen vyhodnotí hodnotu splinu v příslušném bodě bez nutnosti znovu počítat koeficienty.
Kvůli referenci přidám návod, jak vytvořit derivaci splinu, jehož koeficienty jsou uloženy ve struktuře pp.
>> [breaks,coefs,l,k,d] = unmkpp(pp); >> pp1 = mkpp(breaks,repmat(k-1:-1:1,d*l,1).*coefs(:,1:k-1),d); >> dydx = ppval(pp1,x) |
Ve struktuře pp1 je spline, který je derivací splinu pp. Analogicky je možné doplnit i druhou derivaci
>> [breaks,coefs,l,k,d] = unmkpp(pp1); >> pp2 = mkpp(breaks,repmat(k-1:-1:1,d*l,1).*coefs(:,1:k-1),d); >> d2ydx2 = ppval(pp2,x) |
Protože původní spline byl složen z kubických polynomů (4 parametry na každý úsek), jeho první derivace se skládá z kvadratických (3 parametry) a druhá derivace lineárních úseků (2 parametry).
Print This Post
