趙勝喜 祁影霞 周生平 張 華

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

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

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

1 分子動力學的基本原理

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

由牛頓運動定律可得原子的加速度為：

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

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

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

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

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

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

2 勢能模型

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

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

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

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

3 模擬方法

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

4 模擬結果

圖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)熱物理性質Tab.2 Saturated liquid thermodynamic properties of refrigerant ammonia

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

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

5 結論

采用site-site勢能模型對制冷劑氨的飽和液態(tài)密度和比焓進行了分子動力學模擬，模擬結果與NIST數(shù)據(jù)庫的最大相對偏差分別在1.5%(密度)以內和3.2%(比焓)以內。表明采用合理的勢能模型和參數(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]陳正隆，徐為人，湯立達.分子模擬的理論與實踐[M].北京：化學工業(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性質的分子動力學模擬[J].工程熱物理學報,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.