飽和液態(tài)制冷劑氨的熱物理性質(zhì)的分子動力學(xué)模擬

趙勝喜 祁影霞 周生平 張 華

(上海理工大學(xué)能源與動力工程學(xué)院 上海 200093)

氨是最早用于人工制冷的天然制冷劑之一，它的ODP和GWP均為零，而且氨泄漏后能夠被雨水吸收返回土壤成為農(nóng)田的肥料，是一種對環(huán)境無害的綠色制冷劑，廣泛應(yīng)用于國內(nèi)外大中型冷庫及其它低溫制冷設(shè)備。氨的熱物理性質(zhì)的研究至關(guān)重要。傳統(tǒng)的研究制冷劑熱物性的方法主要有狀態(tài)方程擬合和實驗測定，但是這些方法既費時、費力又復(fù)雜。在此，提出運(yùn)用分子動力學(xué)方法模擬制冷劑氨的熱物理性質(zhì)。

分子模擬是依據(jù)統(tǒng)計力學(xué)基本原理，將一定數(shù)量的分子輸入計算機(jī)內(nèi)進(jìn)行分子微觀結(jié)構(gòu)的測定和宏觀性質(zhì)計算的模擬方法。分子動力學(xué)模擬[1](Molecular Dynamics Simulation )是分子模擬中應(yīng)用于預(yù)測物質(zhì)的熱物理性質(zhì)中最常用的一種方法。1971年，Rahman A等人[2]首次運(yùn)用分子動力學(xué)模擬了液態(tài)水的性質(zhì)。1989年，Vega C[3]運(yùn)用分子動力學(xué)模擬了制冷劑R152a的熱物理性質(zhì)，Lusting R[4]模擬了液態(tài)乙烷的性質(zhì)。1993年，Lisal M[5-7]等人模擬了制冷劑R143a、R152a、R142b的熱力學(xué)性質(zhì)，并優(yōu)化了替代制冷劑R32和R23的勢能模型函數(shù)。1997年，Seiji Higashi[8]等人運(yùn)用Lennard-Jones勢能模型模擬了制冷劑R32液態(tài)的熱物理性質(zhì)。2003年，清華大學(xué)的余大啟[9]等人運(yùn)用二中心LJ分子嵌入沿偶極矩(2CLJD)的分子模型，對制冷劑R134a的PVT性質(zhì)進(jìn)行了模擬?；谶@些研究，這里在前人發(fā)展的site-site勢能模型的基礎(chǔ)上，模擬了制冷劑氨的飽和液態(tài)熱物理性質(zhì)。

1 分子動力學(xué)的基本原理

分子動力學(xué)方法是按體系內(nèi)部的內(nèi)稟動力學(xué)規(guī)律來計算并確定位形的轉(zhuǎn)變。考慮含有n個分子的運(yùn)動體系，系統(tǒng)中的能量包括分子的動能與總勢能的和。其總勢能為分子中各原子位置的勢能函數(shù)之和，一般用 表示。根據(jù)經(jīng)典力學(xué)，系統(tǒng)中任一原子所受之力為勢能的梯度：

由牛頓運(yùn)動定律可得原子的加速度為：

將牛頓運(yùn)動定律方程式對時間積分，可預(yù)測原子經(jīng)過時間t后的速度與位置：

式中，U代表原子所受的勢能。 和mi分別代表原子所受之力和質(zhì)量， 和 分別代表原子的速度，加速度和位置。上標(biāo)“0”代表各物理量的初始值。

由于加速度在時間上是連續(xù)函數(shù)，在計算中，我們通常把時間表達(dá)成離散的形式：

式中，tn是原子到第n步的時間；n是步數(shù)；h是時間步長，表示一段非常短的時間間隔(一般取0.8fs)。假定在時間步長h范圍內(nèi)粒子受到的力是不變的，給定粒子的初始位置和初始速度，則可以對該方程進(jìn)行數(shù)值求解。計算過程如下圖1所示：

圖1 數(shù)值求解示意圖Fig.1 Schematic diagram of numerical solution

通過以上的反復(fù)循環(huán)計算，可得到各時間下系統(tǒng)中各原子運(yùn)動的位置、速度及加速度等數(shù)據(jù)。然后再利用統(tǒng)計計算方法得到系統(tǒng)的靜態(tài)和動態(tài)特性，從而得到系統(tǒng)的宏觀性質(zhì)，比如壓力，密度，內(nèi)能，焓等熱物理性質(zhì)。

2 勢能模型

氨的勢能模型采用site-site勢能模型[10]，分子間的相互作用力包括Lennard-Jones(6-12)項和庫侖作用項，其形式如下：

其中，r是不同分子中的原子間的距離，σ和ε分別是分子中不同原子間的L-J勢能的能量和尺度參數(shù)。q是分子中各原子所帶的部分電荷，ε0是自由介電常數(shù)，e是電荷的單位。下標(biāo)a、b分別代表不同分子中的原子對。

對于氨分子中的不同原子的σ、ε和q可以運(yùn)用量子力學(xué)的從頭計算方法獲得，模擬采用Bernhard Eckl[11]等人的優(yōu)化結(jié)果，如表1所示。

表1 氨分子的勢能參數(shù)Tab.1 Potential function parameters of ammonia

3 模擬方法

分子動力學(xué)計算之前，必須先估計計算的可行性，一般來說系統(tǒng)的原子個數(shù)越多，模擬的結(jié)果越精確，但是以目前的計算機(jī)水平，分子動力學(xué)計算的容量約為數(shù)千個原子的系統(tǒng)。此次模擬選用500個氨分子進(jìn)行面心立方晶格建模，首先采用等溫等容(NVT)系綜使體系達(dá)到平衡，然后在等溫等壓 (NPT)系綜下得到模擬結(jié)果。對建立的面心立方晶格模型采用周期性邊界條件，利用弛豫作用隨機(jī)得到分子的初始位置。針對各分子間的Lennard-Jones項采用硬球截斷法，截斷半徑為立方晶體邊長的一半，長程靜電力則采用Ewald加和法計算。運(yùn)動方程采用基于預(yù)測-校正積分方法的Gear算法，溫度采用速度標(biāo)定法控制。時間步長采用0.8fs, 模擬步程約為50000步，其中前30000步用于體系達(dá)到平衡，后20000步用于得到氨的熱物理性質(zhì)。

4 模擬結(jié)果

圖2 制冷劑氨的飽和液態(tài)T-h圖線Fig.2 Saturated liquid T-h of refrigerant ammonia

圖3 制冷劑氨的飽和液態(tài)T-ρ圖線Fig.3 Saturated liquid T-ρ of refrigerant ammonia

表2 制冷劑氨的飽和液態(tài)熱物理性質(zhì)Tab.2 Saturated liquid thermodynamic properties of refrigerant ammonia

運(yùn)用分子動力學(xué)方法模擬了制冷劑氨的飽和液態(tài)熱物性，并將模擬結(jié)果與美國國家標(biāo)準(zhǔn)研究所(NIST)的Refprop8.0數(shù)據(jù)庫提供的制冷劑氨的飽和液態(tài)熱物性質(zhì)進(jìn)行了比較，其中比焓是以0℃的飽和液態(tài)氨的比焓為基準(zhǔn)得到(h=-762.75kJ/kg)。主要結(jié)果如表2，圖2和圖3所示。

從表2、圖2及圖3可知，模擬結(jié)果同美國國家標(biāo)準(zhǔn)研究所(NIST)的數(shù)據(jù)庫提供的氨的飽和液態(tài)熱物理性質(zhì)有很好的一致性。氨飽和液態(tài)密度的最大相對偏差在1.5%以內(nèi)，比焓的最大相對偏差在3.2%以內(nèi)，低壓區(qū)的模擬結(jié)果比高壓區(qū)的模擬效果好。

5 結(jié)論

采用site-site勢能模型對制冷劑氨的飽和液態(tài)密度和比焓進(jìn)行了分子動力學(xué)模擬，模擬結(jié)果與NIST數(shù)據(jù)庫的最大相對偏差分別在1.5%(密度)以內(nèi)和3.2%(比焓)以內(nèi)。表明采用合理的勢能模型和參數(shù)，運(yùn)用分子動力學(xué)方法來預(yù)測單一組分工質(zhì)的熱物理性質(zhì)是可行的。

本文受上海市重點學(xué)科建設(shè)項目(S30503)，上海市教育委員會科研創(chuàng)新項目(11YZ119)及上海市研究生創(chuàng)新基金項目(JWCXSL1102)資助。(The project was supported by Shanghai Leading Academic Discipline Project (No.S30503),Shanghai Education Commission Scienti fi c Research Innovation Projects (No.11YZ119) and The Innovation Fund Project For Graduate Student of Shanghai(No.JWCXSL1102).)

[1]陳正隆，徐為人，湯立達(dá).分子模擬的理論與實踐[M].北京：化學(xué)工業(yè)出版社，2007.

[2]Rahman A, Stillinger F H. Molecular dynamics study of liquid water [J]. Chem Phys, 1971, 55: 3336-3359.

[3]Vega C, Saager B, Fischer J. Molecular dynamics studies for the new refrigerant R152a with simple model potentials[J]. Molecular Physics, 1989, 68(5): 1079-1093.

[4]Lusting R, Toro-Labbe A, Steele W A. Molecular dynamics study of the thermodynamics of liquid ethane [J].Fluid Phase Equilibria, 1989, 48:1-7.

[5]Lisal M, Budinsky R, Vacek V, et al.Vapor-liquid equilibria of alternative refrigerants by molecular dynamics simulation [J].Int Thermophysics,1993,20(1):163-169.

[6]Lisal M, V Vacek. Effective potentials for liquid simulation of the alternative refrigerants HFC-32:CH2F2and HFC-23: CHF3[J]. Fluid Phase Equilibria, 1996, 118:61-76.

[7]Lisal M, Vacek V. Molecular dynamics simulations of fluorinated ethane [J]. Molecular Physics, 1996,87(1):167-187.

[8]Seiji Higashi,Akira Takada.Molecular dynamics study of liquid CH2F2(HFC-32)[J].Molecular Physics,1997,92:641-650.

[9]余大啟, 李卓毅, 曹炳陽, 等. R134a PVT性質(zhì)的分子動力學(xué)模擬[J].工程熱物理學(xué)報,2003,24(6):24-26.(Yu Daqi, Li Zhuoyi,Cao Bingyang, et al.Molecular dynamics studies about the PVT properties of R134a[J].Journal of Engineering Thermophysics,2003,24(6):24-26.)

[10]T Kristof, J Vorholz, J Liszi, et al. A simple effective pair potential for the molecular simulation of the thermodynamic properties of ammonia [J]. Molecular Physics, 1997, 97: 1129-1137.

[11]Bernhard Eckl, Martin Horsch, Jadran Vrabec, et al.Molecular modeling and simulation of thermophysical properties: application to pure substances and mixtures[G]// High Performance Computing in Science and Engineering'08, 2009,3:119-133.