Cvièení { 2D Clausiova-Clapeyronova rovnice
1/12 pch03
Evropský sociální fond þPraha & EU: Investujeme do va¹í budoucnostiÿ Inovace pøedmìtu Poèítaèová chemie je podporována projektem CHEMnote (Inovace bakaláøského studijního programu Chemie { moderní vzdìlávání podpoøené pou¾itím notebookù { CZ.2.17/3.1.00/33248) v rámci Operaèního programu PRAHA { ADAPTABILITA.
Clausiova{Clapeyronova rovnice ve 2D
Úkol:
2/12 pch03
Ovìøit platnost Clausiovy{Clapeyronovy rovnice na 2D modelu kapaliny a plynu
Model:
Potenciál typu 8-4 (≈ Lennard Jones ve 2D): 1 1 u(r) = 8 − 4 r r
Pevné pøita¾livé èi odpudivé zdi (potenciál typu 7{3, resp. jen jeho odpudivá èást ∝ 1/r7). Redukované jednotky: kB = R/NA = 1, energie a teplota mají stejné jednotky Velièiny budeme udávat na 1 atom, ne na 1 mol (index at)
Postup
3/12 pch03
V systému s dvìma fázemi oddìlenými rovinným rozhraním stanovte tlak nasycených par v závislosti na teplotì pro dvì teploty. Pou¾ijte MD s libovolným termostatem. Stanovte støední teplotu a tlak, simulací páry v MC pak stanovte kompresibilitní faktor. Z Clausiovy{Clapeyronovy rovnice (s opravou na neidealitu páry) stanovte výparnou entalpii vèetnì odhadu chyby. Stanovte výparnou entalpii ze støední potenciální energie kapaliny v periodických okrajových podmínkách. Porovnejte obì hodnoty.
4/12 pch03
Simulaèní metody
Simulace startuje automaticky z náhodné kon gurace pomocí MC (odstranìní pøekryvù molekul), pak se automaticky pøepne na MD. MD simulace (leap-frog) + Berendsenùv termostat (lze pou¾ít i jiný). Lze volit i Monte Carlo (Metropolis) { pro nìkteré výpoèty vhodnìj¹í. Dvì metody mìøení (tensoru) tlaku (ve smìru y): Z prùmìrné síly pùsobící na pevnou stìnu: pstìna =
Z viriálu sil:
fstìna L
,
pyy = ρkBT + ρ = N/V , V = L2, L
L
1 DV
= délka stìny
DX
ryfy
E
= délka stìny, D = dimenze, sèítá se pøes v¹echny párové síly (èástice{èástice, stìna{èástice) h·i = prùmìrování okam¾itých hodnot v prùbìhu simulace
5/12
Výparná entalpie z Clausiovy{Clapeyronovy rovnice pch03
Jak jistì víte, Clausiova{Clapeyronova rovnice ln(p1/p2) ∆vapHm = − 1/RT1 − 1/RT2 se odvodí z Clapeyronovy rovnice s pou¾itím následujících zjednodu¹ení: Výparná entalpie nezávisí na teplotì. Zanedbáváme objem kapaliny proti objemu plynu. Platí stavová rovnice ideálního plynu. Za podmínek simulace jsou chyby prvních dvou pøedpokladù malé (2 %), ale poslední pøedpoklad je velmi nepøesný (a¾ 15 %). Pøesnìj¹í je pøedpoklad: Kompresibilitní faktor plynu Z = p/ρkBT za tlaku nasycených par nezávisí na teplotì. Korigovaná Clausiova{Clapeyronova rovnice: ln(p1/p2) ∆vapHat ≈ −ZkB 1/T1 − 1/T2 kde Z aproximujeme hodnotou v T = (T1 + T2)/2
Výparná entalpie ze støední potenciální energie
6/12 pch03
Stanovíme ji ze vzorce Epot ∆vapHat = ∆vapUat + pVat = − N
+ ZkBT
SIMOLANT { instalace (Windows)
http://www.vscht.cz/fch/software/simolant Stáhnìte simolant-win32.zip Rozbalte do vhodné slo¾ky (nespou¹tìjte pøímo z Spus»te simolant.exe
simolant-win32.zip)
7/12
Tlak nasycených par { simulace pro T1 pch03 Na pomalej¹ím poèítaèi sni¾te poèet atomù (slider \N"), ale ne pod 150 Menu: Prepare system → Vapor-liquid equilibrium
Menu: Measure,show,record → Extended (Energy, Temperature, Pressures, γ) Slider \simulation speed" (vpravo dole) dejte na maximum (zobrazuje a zpracovává se pouze ka¾dá 15. kon gurace) Slider \measurement block" dejte na max. (blok = prùmìr z 100 bodù) Nastavte teplotu (slider \T" { ne \τ") na T1 ∈ (0.15, 0.16) { èím ni¾¹í teplota, tím rychlej¹í poèítaè potøebujete { Tip: jemný pohyb slideru: ctrl-↑ a ctrl-↓ { Tip: také lze zadat pøíkaz do okénka cmd: T=0.155 + Enter Nechte dùkladnì ustálit Stisknìte record v panelu \Measure". Nemìòte parametry simulace v prùbìhu mìøení! Po chvíli stisknìte record znovu. Zobrazí se výsledky. Vhodný poèet blokù (n= vpravo nahoøe) je aspoò 50, lépe pøes 100; relativní chyba velièiny P(top wall) pøíp. Pyy by mìla být men¹í ne¾ cca 10 %. Nevyhovuje-li, stisknìte continue . Jste-li spokojeni, ulo¾te pomocí save (overwrite... .
8/12
Tlak nasycených par { simulace pro T2 pch03 Opakujte pro teplotu T2 ∈ (0.19, 0.20) { staèí cca polovina blokù, proto¾e tlak za vy¹¹í teploty je vy¹¹í a statistická chyba men¹í (roste ale neidealita) Zapi¹te pomocí record → append to... . Analýza dat I
Výsledky najdete v souboru simolant.txt. Pokud jste vícekrát stiskli append to... , najdete v souboru více sad výsledkù. Musíte se vyznat. . . Najdìte hodnotu p1 pro teplotu T1: P (top wall) pøíp. Pyy z 1. tabulky (výsledky by mìly být ≈ stejné) Stejnì najdìte p2 pro teplotu T2 z 2. tabulky Vypoètìte støední teplotu a tlak takto: T + T2 T= 1 ,
p=
√ p1 p2
2 Vypoètìte þstøední èíselnou hustotuÿ páry podle stavové rovnice ideálního plynu: ρ = p/T
Výpoèet kompresibilitního faktoru (Monte Carlo)
9/12 pch03
provedeme v periodických okrajových podmínkách Menu: Boundary conditions → Periodic Menu: Measure,show,record
→
Extended (Energy, Temperature, Pressures, γ)
Menu: Simulation method → Monte Carlo (Metropolis) Nastavte teplotu na T = (T1 + T2)/2 (nejlépe jako \T=èíslo " v okénku \cmd:") Nastavte hustotu na ρ: (nejlépe jako \rho=èíslo " v okénku \cmd:") Nechte ustálit Stisknìte record a simulujte aspoò 10 blokù V pøíslu¹ném bloku dat v simolant.txt najdìte hodnotu Z Chcete-li být pøesnìj¹í, zkuste zmìnit hustotu tak, aby se výsledný tlak rovnal p, napø. zvolte novou hustotu podle ρ = ρ/Z a opakujte Standardní je pou¾ít simulaci v NPT souboru se zadaným tlakem p = p, tento soubor není v¹ak z rùzných dùvodù v SIMOLANTu implementován
10/12 pch03
Pokud se budete nudit. . . MC nebo MD?
Opakujte pøedchozí výpoèet s Menu: Simulation method → MC→MD (Berendsen) Porovnejte chyby obou metod. Proè je zde MD ménì pøesná? Pokud se budete nudit. . . Alternativní postup získání Z
ρ(y)
Místo simulace páry v periodických okrajových podmínkách je mo¾né stanovit hustoty páry rovnou v obou simulacích rovnováhy: Menu: Measure,show,record → extended + Vertical density pro le 1 V souboru simolant.txt najdete hustotní pro l, který zobrazíte kapalina (napø. po naètení do Excelu) a z èásti odpovídající plynu získáte 0.5 jeho hustotu ρ1, pak Z1 = p1/ρ1 a stejnì pro druhou teplotu. Z obou hodnot Z pak udìláte plyn prùmìr. 0 0
10
20 y
30
40
Výpoèet výparné entalpie z tlakù nasycených par
11/12 pch03
Ze získaných hodnot vypoètìte (pou¾íváme kB = 1) ln(p1/p2) ∆vapHat = −Z 1/T1 − 1/T2 Nezapomeòte spoèítat chybu výsledku! Ve výsledcích jsou uvedeny odhady standardních chyb∗ stanovené z blokù dat. Staèí uva¾ovat chyby v p1 a p2, proto¾e chyby v teplotách a Z jsou relativnì malé: q δ(∆vapHat) =
∗
δrel(p1)2 + δrel(p2)2 |1/T1 − 1/T2|
standardní chyba = smìrodatná odchylka prùmìru
12/12 pch03
Výparná entalpie z vniøní energie
Menu: Boundary conditions
→
Periodic
Menu: Measure,show,record → Extended (Energy, T., P., γ) to u¾ máte nastaveno Nastavte teplotu na T = (T1 + T2)/2 Mù¾ete pou¾ít MC i MD Pro vìt¹í rychlost mù¾ete trochu zmen¹it slider \measurement block" Posunujte (nejprve rychle, pak pomalu) slider \ρ" (hustota), a¾ tlak (P=) kolísá okolo nuly (pøesnì má být p = √p1p2). Kon gurace musí být homogenní tekutina bez dìr!
Zaznamenejte hodnotu Epot (chyba je oproti výpoètu z tlakù zanedbatelná) Spoèítejte
∆vapHat = −
Epot N
+ ZT
(pou¾íváme
kB = 1)
Srovnejte s hodnotou podle Clausiovy{Clapeyronovy rovnice
