Kinetika fázového přechodu pára–kapalina a její numerická simulace

Hjshj Rwlhh^{†}, Hjshj Rwlhh^{‡}

† Ústav termomechaniky AVČR ‡ Ústav termomechaniky AVČR

23. ledna 2001

Seznam tabulek

1 Teplotní rozsahy experimentálních dat z [11] použité při tvorbě funkcí fitujících druhý resp. třetí viriální koeficient. . . . . . . . . . . . . . . . 4

2 Teplotní rozsahy experimentálních dat z [11] použité při tvorbě funkcí fitujících druhý resp. třetí viriální koeficient. . . . . . . . . . . . . . . . 5

Abstrakt

Je vytvořen jednoduchý model povrchového napětí sféricky zakřiveného povrchu a je ukázáno, že jeho aplikace v Lotheho–Poundově modelu homogenní nukleace značně zmenší dřívější nesoulad tohoto modelu s experimentálními výsledky pro převážně polární látky. Dále je sestaven systém rovnic popisující proudění s homogenní kondenzací přičemž se respektuje skutečnost, že v objemovém elementu tekutiny se mohou nalézat kapky různých velikostí. Dále je vytvořen kód pro simulaci jednorozměrného proudění s homogenní kondenzací při použití modelu růstu kapky v přiblížení proudu volných molekul a při použití Hillovy aproximace. Výsledky simulací jsou porovnány s výsledky experimentů pro studium kondenzace na rázové trubici provedené v ÚT AVČR. Rozdíly mezi výsledky simulací a experimentálním výsledkem jsou diskutovány.

Přehled současného stavu

Homogenní kondenzace se skládá ze dvou charakteristických stádií – homogenní nukleace a růstu kapky. Model růstu kapky (alespoň v případě, kdy kondenzuje jen jedna látka) je poměrně dobře propracován [1] a není bezprostředním cílem našeho zkoumání.

Teoreticky zajímavějším je stádium nukleace. Je to tvoření velmi malých kapek – zárodků – nuklejí. Tyto zárodky mají velkou pravděpodobnost, že dále porostou a jejich počet a velikosti jsou klíčové pro celý proces kondenzace. V jednotce objemu se za jednotku času vytvoří jisté množství těchto zárodků, které se nazývá nukleační rychlost. Jedná se o stádium, kdy hrajou roli individuální vlastnosti molekul, kontinuální popis ţlátky�, z níž se zárodek sestává, přestává být adekvátní a objevují se různé, v mnoha otázkách navzájem se rozcházející teorie^{1}. Nukleační rychlost závisí na podmínkách (tlaku a teplotě) a na látkových vlastnostech. Typickým znakem současných nukleačních teorií je snaha vystačit s makroskopicky dobře měřitelnými vlastnostmi – např. makroskopickou hustotou, křivkou nasycených par a povrchovým napětím. Historicky nejstarší je tzv. klasická teorie nukleace, jejíž vznik je spojen se jmény Becker, Döring, Zeldovič a Frenkel. Tato teorie obstála při kvantitativním popisu nukleační rychlosti vody. V případě nepolárních látek je však souhlas s experimentem zcela neuspokojivý^{2}. Později se vyskytly vážné námitky proti její konzistentnosti se statistickou termodynamikou a vznikly teorie které jsou se statistickou termodynamikou konzistentnější^{300}. Jedná se např. o Lotheho–Poundovu teorii nukleace [2] nebo teorii založenou na Fisherově modelu volné enerie klastrů [3]. Nedávno se objevil pokus o ţnápravu� mnoha rozporů klasické teorie nukleace její mírnou modifikací [4]. V poslední době se objevily práce, které se opírají o Fisherův výraz pro volnou klastru s povrchovým napětím závislým na velikosti klastru [5, 6]. Avšak ani tyto konzistentnější teorie nepopisují uspokojivě nukleaci současně polárních i nepolárních látek. Tyto důvody neustále iniciují další experimentální i teoretický výzkum homogenní nukleace. Podrobný výklad nejdůležitějších nukleačních teorií bude uveden ve druhé kapitole s jistým důrazem na historické souvislosti.

1první poznámka

2druhá poznámka

I. Cíl práce

Cílem této práce je:

1. prohloubit znalosti o podstatě fázového přechodu plyn–kapalina na mikroskopické úrovni

2. sestavení systému rovnic, který dostatečně přesně popisuje chování systému s prouděním a kondenzací pro účely studia kondenzace pomocí rázové trubice. Tento systém musí být přitom stále řešitelný současnými výpočetními prostředky (alespoň pro případy jednodimenzionálního a dvoudimenzionálního proudění)

3. vytvoření spolehlivého počítačového programu pro případ jednodimenzionálního proudění, který by měl sloužit pro lepší řízení experimentů na rázové trubici prováděné v současné době v ÚT AVČR a srovnáním naměřených a vypočtených hodnot umožnit verifikaci různých fyzikálních modelů kondenzace, především pak jeho důležité fáze – nukleace.

2 Fisherův polofenomenologický model klastrů

Fisher vychází ve své studii Helmholtzovy volné energie klastrů z mřížkového modelu plynu (Lattice Gas modelu). Mřížkový model plynu [7] stručně řečeno modeluje reálný plyn jako by molekuly mohly být pouze v mřížkových polohách jakési pomyslné mřížky, jejíž mřížková konstanta odpovídá hustotě kapalné fáze. Molekula buď může být v uzlovém bodu mřížky přítomna nebo ne. Předpokládá se, že spolu interagují jen sousední molekuly (pokud jsou v příslušných mřížkových bodech přítomny) a předpokládá se, že tato interakce je popsána jen konstantní interakční energií. Molekuly, které jsou od sebe dále spolu neinteragují. Numericky se dá spočítat stavová suma – sumuje přes všechny možné konfigurace molekul a příspěvky ke stavové sumě odpovídající různým stavům v impulsovém prostoru se při nekvantovém (přesněji řečeno kvaziklasickém) popisu dají převést na jednoduchý integrál, který se dá z celé stavové sumy vytknout a snadno spočítat. I když je tento model reálného plynu v podstatě dost primitivní, přesto respektuje nejen rotační a translační stupně volnosti, ale i replacement free energy a vlastně i všechny možné tvary klastrů!

Fisher [8] na základě zkušenností z výpočtu stavové sumy metodou Monte–Carlo navrhuje následující tvar pro Helmholtzovu energii klastru

F_{i} = if_{∞} + σs_{1}i^{η} + τkT lni − kT ln(q_{0}V)  (1)

Využívá dále představy, že neideální plyn je vlastně směsí ideálních ţplynů klastrů� různých velikostí (Mayer conjecture [9]) a pro celkový počet molekul v jednotce objemu c a tlak p_{v} platí

c = \frac{1}{V} ∑_{i=1}^{∞}ie^{− \frac{ΔG(i)}{kT}} ,  (2)

p_{v} = \frac{kT}{V} ∑_{i=1}^{∞}e^{− \frac{ΔG(i)}{kT}} .  (3)

Předpokládá, že tato představa je dobře použitelná až do okolí kritického bodu. Provádí dále analýzu chování spojení rovnic (x.33), (1), (2) a (3) a pomocí teorie kritických exponentů ukazuje, že

δ = \frac{1}{τ − 2},  (4)

kde δ je koeficient z kritické izotermy

p_{c} − p ∼ (ρ_{c} − ρ)^{δ} .  (5)

300 třístá poznámka

ρ je hustota a intex c znamená hodnotu v kritickém bodě. Z numerických výsledků mu vychází

η ∼ 0.615 ∼ 0.640,  2\frac{5}{26} < τ < 2\frac{1}{5}. (6)

Van der Waalsova rovnice podle teorie kritických exponentů dává [8]

η = \frac{2}{3},  τ = 2\frac{1}{3}.  (7)

V kritickém bodě podle (2) (s využitím toho, že povrchové napětí je v kritickém bodě nulové a tlak je roven tlaku nasycených par) platí

ρ_{c} = q_{0} ∑_{i=1}^{∞}i^{1−τ}.  (8)

p_{c} = q_{0}kT ∑_{i=1}^{∞}i^{−τ} .  (9)

Řada (8) je vlastně Riemannova funkce proměnné 1 − τ a konverguje pro τ > 2.

Kiang [10] z experimentálních hodnot kritického tlaku a hustoty (dvě hodnoty) určuje obě konstanty τ, q_{0} pomocí spojení rovnic (8) a (9) . Ukazuje, že koeficient δ z rovnice (4) souhlasí velmi dobře s experimentální hodnotou pro nepolární látky. Kiang a další [3] pak navrhují použít tohoto postupu určení τ a q_{0} pro výpočet Gibbsovy energie tvorby klastru a pro výpočet nukleační rychlosti (podle kinetického schématu a aproximací, jako v klasické teorii) jak pro polární, tak pro nepolární látky.

3 Dillmannův–Meierův model nukleace

Dillmann a Meier [5] navrhli vylepšení vztahu (1) tak, že povrchové napětí se uvažuje závislé na počtu molekul v klastru a navrhují

σ_{i} = (1 + α_{1}i^{− \frac{1}{3}} + α_{2}i^{− \frac{2}{3}} ) σ_{∞},  (10)

kde σ_{∞} je povrchové napětí rovného povrchu a α_{1}, α_{2} jsou funkce teploty. Podle (x.9), (x.34), (1) a (2) lze pro malé tlaky pro kvazirovnovážné rozdělení psát

c = q_{0} ∑_{i=1}^{∞}(\frac{p_{v}}{p_{s∞}} )^{i} e^{− \frac{σ_{i}s_{i}}{kT} −(τ−1) ln i}.  (11)

Z této rovnice se dá odvodit vztah mezi korekčním součinitelem pro povrchové napětí klastru sestávajícího se z i molekul podle definice

κ_{i} ≡ \frac{σ_{i}}{σ_{∞}}  (12)

a derivací tlaku podle hustoty při stálé teplotě a nulovém tlaku

κ_{i} = −\frac{1}{ϑi^{2/3}} ln[\frac{p_{s∞}^{i}}{q_{0}}i^{τ−1} \frac{1}{i!} lim_{p→0} (\frac{∂^{i}c}{∂p^{i}} )_{T} ] , (13)

kde

ϑ = \frac{σ_{∞}s_{1}}{kT}  (14)

Derivace se dají určit z empirické viriálové stavové rovnice

\frac{p_{v}}{c} = kT + Bp_{v} + Cp_{v}^{2} + ....  (15)

Po dosazení těchto derivací do (13) se dostane

κ_{1} = −\frac{1}{ϑ} ln \frac{p_{s∞}}{q_{0}kT}  (16)

a

κ_{2} = −\frac{1}{2^{\frac{2}{3}} ϑ} ln[− (\frac{p_{s∞}}{q_{0}kT} )^{2} 2^{τ−1}q_{0}B]. (17)

Podmínkou (16) je zaručen správný počet monomerů při nízkých tlacích. Koeficienty α_{1} a α_{2} se určí z podmínky, aby vztah (10) dal správné hodnoty κ_{1} a κ_{2} [5]:

α_{1} = \frac{(κ_{2} − 1) − (κ_{1} − 1)2^{− \frac{2}{3}}}{2^{− \frac{1}{3}} − 2^{− \frac{2}{3}}}  (18)

a

α_{2} = −\frac{(κ_{2} − 1) − (κ_{1} − 1)2^{− \frac{1}{3}}}{2^{− \frac{1}{3}} − 2^{− \frac{2}{3}}}  (19)

Tím jsou oba koeficienty ve výrazu pro korekci povrchového napětí (10) určeny z podmínky správného počtu monomerů při nízkých tlacích a druhého viriálního koeficientu B^{3}. Analogickým postupem, jako v Kiangově metodě (jediný rozdíl je vlastně jen v aplikaci korekce povrchového napětí)

3S myšlenkou modifikovat Kiangovu metodu změnou povrchového napětí určenou z druhého viriálního koeficientu (avšak nepříliš dokonalou) přišel již dříve Hamill a další [12]

látka  rozsah pro B’ [K] rozsah pro C’ [K]
Voda  373–473  473–598
Methanol  298–423  423–498
Čpavek  273–396  ——
Chloroform  316–398  ——
Tetrachlormethan  318–350  353–419

Tabulka 1: Teplotní rozsahy experimentálních dat z [11] použité při tvorbě funkcí fitujících druhý resp. třetí viriální koeficient.

se pak použitím klasických aproximací kinetiky, které vedly ke vztahu (x.20), dostane (viz. [5])

J = \frac{1}{3}\sqrt{\frac{ϑ}{π} (1 + α_{1}i^{∗}^{− \frac{1}{3}} + \frac{9τ}{2ϑ}i^{∗}^{− \frac{2}{3}} )} ×

×\frac{s_{1}p_{v}}{\sqrt{2πm_{1}kT}} e^{− \frac{Δ^{∗}}{kT}} ,  (20)

přičemž počet molekul v kritickém klastru i^{∗} je určen reálným kořenem kubické rovnice

τx^{3} + \frac{1}{3}α_{1}ϑx^{2} + \frac{2}{3}ϑx − ln \frac{p_{v}}{p_{s∞}} = 0,  (21)

kde

x = i^{∗}^{− \frac{1}{3}} .  (22)

Veličina Δ^{∗} je pomocí Gibbsovy energie tvorby kritického klastru ΔG^{∗} ((2.34) pro i = i^{∗} a F_{i} je dáno vztahem (1)) definována následujícím vztahem

c(i^{∗}) = \frac{1}{V} e^{− \frac{ΔG^{∗}}{kT}} = e^{− \frac{Δ^{∗}}{kT}} ,  (23)

nebo–li

Δ^{∗} ≡ kT lnV + ΔG^{∗}  (24)

Tento model dává při vyšších přesyceních pro mnoho látek realističtější hodnoty než klasický model (x.23) [5]. V poslední době se objevila kritika Dilmann–Meier modelu – viz. např. [13]. Autoři ukazují, že předpoklad (x.9) založený na představě ideálního plynu, vede k velké chybě při fitování κ_{2} na druhý viriální koeficient a odvozují korekci povrchového napětí bez tohoto předpokladu. Finální vzorce pak vedou k podstatě horšímu souladu s experimentem, než vzorce Dillmanna a Meiera.

Je užitečné si všimnout, že Dillmannova– Meierova metoda určení korekčních součinitelů α_{1} a α_{2} vyžaduje, aby poměrně mírná závislost (10) vystihovala korekci povrchového napětí (a tím Gibbsovu energii tvorby klastru) až do extrémně malých klastrů (o jedné molekule), kde zřejmě přestává mít pojem povrchového napětí dobrý smysl. To jestli je závislost (10) použitelná pro takto malé klastry je podle mého názoru velmi diskutabilní. TOTO je pokusný odstavec demonstrující, jak

Obrázek

lze pomocí TeXovského příkazu \parshape, příkazu \marginpar a vhodným nastavením délkových registrů marginparwidth a marginparsep

dosáhnout ţvkousnutí� objektu do textu ze strany. Je ovšem potřeba také trochu manipulovat s posazením objektu pomocí příkazu \raisebox. Pokud jste soubor art.tex upravovali, je dost pravděpodobné, že poloha ţobrázku� není v pořádku. Pak je potřeba vhodně nastavit marginparsep a možná i parametr zvýšení v příkazu \raisebox. Pro pohodlnější práci je pro takovéto účely určen specializovaný balík picinpar.

Následuje další odstavec. Ze vzhledu je vidět, že marginparsep a marginpar jsme nastavili na hod-

pozn. na okraji

noty platné před předcházejícím odstavcem.

A Materiálové vlastnosti použité k výpočtům poprvé

\frac{B′p_{c}}{RT_{c}} = f_{0}(T_{r}) + ωf_{1}(T_{r})  (A.1),

kde

f_{0}(T_{r}) = 0.1445 − \frac{0.330}{T_{r}} − \frac{0.1385}{T_{r}^{2}} − \frac{0.0121}{T_{r}^{3}}

−\frac{0.000607}{T_{r}^{8}},

f_{1}(T_{r}) = 0.0637 − \frac{0.331}{T_{r}^{2}} − \frac{0.423}{T_{r}^{3}} − \frac{0.008}{T_{r}^{8}}.

\frac{B′p_{c}}{RT_{c}} = \frac{9}{128} (1 − \frac{6}{T_{r}^{2}} ) − aT_{r}e^{\frac{b}{T}}  (25)

látka  rozsah pro B’ [K] rozsah pro C’ [K]
Voda  373–473  473–598
Methanol  298–423  423–498
Čpavek  273–396  ——
Chloroform  316–398  ——
Tetrachlormethan  318–350  353–419

Tabulka 2: Teplotní rozsahy experimentálních dat z [11] použité při tvorbě funkcí fitujících druhý resp. třetí viriální koeficient.

TOTO JE POKUSNÝ ODKAZ: první řádek hned za rovnicí 1 je v odstavci 2 na straně 2.

B Materiálové vlastnosti použité k výpočtům podruhé

\frac{B′p_{c}}{RT_{c}} = \frac{9}{128} (1 − \frac{6}{T_{r}^{2}} ) − aT_{r}e^{\frac{b}{T}}  (26)

Materiálové vlastnosti použité k výpočtům potřetí

\frac{B′p_{c}}{RT_{c}} = \frac{9}{128} (1 − \frac{6}{T_{r}^{2}} ) − aT_{r}e^{\frac{b}{T}}  (27)

Reference

[1] G.Gyarmathy. The spherical droplet in gaseous carrier streams. In G.F.Hewitt, J.M.Delhaye, and N. Zuber, editors, Multiphase Science and Technology. Hemisphere Publishing Corporation, McGraw-Hill C., 1982.

[2] J.Lothe and G.M.Pound. Reconsiderations of nucleation theory. Journal of Chemical Physics, 36:2080–2085, 1962.

[3] C.S.Kiang,  D.Stauffer,  G.H.Walker,

O.P.Puri, J.D.Wise Jr., and E.M.Patterson. J.Atmospheric Sci., 28:1222, 1971.

[4] S.L.Girshick and Chia-Pin-Chiu. Kinetic nucleation theory: A new expression for the rate of homogeneous nucleation from an ideal supersaturated vapor. Journal of Chemical Physics, 93:1273–1277, 1990.

[5] A.Dillmann and G.E.Meier. A refined droplet approach to the problem of homogeneous nucleation from the vapor phase. Journal of Chemical Physics, 94:3872–3884, 1991.

[6] V.I.Kalikmanov and M.E.H. van Dongen. Semiphenological theory of homogeneous vapor– liquid nucleation. Journal of Chemical Physics, 103:1–6, 1995.

[7] R.J.Baxter. Exactly Solved Models in Statistical Mechanics. Academic Press, London, New York, 1982.

[8] M.E.Fisher. The theory of condensation and the critical point. Physics, 3:255–283, 1967.

[9] J.E.Mayer and M.G.Mayer. Statistical Mechanics. John Wiley, New York, 1940.

[10] C.S.Kiang. Use of liquid–droplet model in calculations of the critical exponent δ. Physical Review Letters, 24:47–50, 1970.

[11] J.H.Dymond and E.B.Smith. The Virial Coefficients of Pure Gases and Mixtures. Clarendon Press, Oxford, 1980.

[12] P.Hamill, D.Stauffer, and C.S.Kiang. Nucleation theory: Fisher’s droplet picture and microscopic surface tension. Chemical Physics Letters, 28:209–212, 1974.

[13] A.Laaksonen, I.J.Ford, and M.Kulmala. Revised parametrization of the Dillmann–Meier theory of homogeneous nucleation. Physical Review E, 49:5517–5524, 1994.