吳　東，祁影霞，楊　喜

（上海理工大學(xué)，上海 200093）

?

HFO-1234yf/HFC-32氣-液相平衡特性研究

吳東*，祁影霞，楊喜

（上海理工大學(xué)，上海 200093）

[摘 要]本文簡單介紹了吉布斯蒙特卡洛（GEMC）模擬法原理并采用GEMC模擬了混合制冷劑HFO-1234yf/HFC-32質(zhì)量分?jǐn)?shù)分別為50％/50％和80％/20％，在溫度(273～333) K的氣液相平衡特性。本文計(jì)算結(jié)果與PR狀態(tài)方程計(jì)算值的比較表明，飽和氣液相密度相對(duì)偏差分別為-3％～+2％、-4.5％～+3.5％。結(jié)果表明，GEMC模擬方法誤差小，為進(jìn)一步研究混合制冷劑HFO-1234yf/HFC-32氣液相平衡特性奠定了基礎(chǔ)。

[關(guān)鍵字]模擬；HFO-1234yf/HFC-32；氣液相平衡；狀態(tài)方程

*吳東（1988-），男，碩士在讀，主要從事制冷劑熱物性研究。通信作者：祁影霞，女，副教授，地址：上海市楊浦區(qū)軍工路516號(hào)上海理工大學(xué)。E-mail：qipeggy@126.com。

0 引言

隨著CFCs和HCFCs類制冷劑逐步的禁止使用，國際上形成了兩種不同的制冷劑替代研究路線[1]：一條以美日為代表，主要采用合成HFCs類制冷劑，如R134a、R410A和R407C等；另外一條以德國和北歐一些國家為代表，主張使用天然制冷劑，如CO2、H2O和HCs等。當(dāng)前，汽車空調(diào)依然廣泛采用R134a作為制冷劑，主要是基于其ODP為0，安全性能高；根據(jù)歐盟F-gas法規(guī)，自2011年起禁止使用GWP大于150的制冷劑，而其GWP為1,430[2-4]。陸岷山等[5]通過實(shí)驗(yàn)的方法對(duì)混合制冷劑R134a/R600a的PVT性質(zhì)進(jìn)行了實(shí)驗(yàn)研究。

從上可以看出，由于理想制冷劑應(yīng)符合熱力學(xué)性質(zhì)、遷移性質(zhì)、物理化學(xué)性質(zhì)和其他性質(zhì)等諸多方面的要求，因此新型混合制冷劑的研究將有助于找到更為理想的替代制冷劑[6-8]。近年，由于計(jì)算機(jī)技術(shù)的巨大進(jìn)步，建立在量子力學(xué)和統(tǒng)計(jì)熱力學(xué)基礎(chǔ)上的分子模擬方法已成為預(yù)測新型混合制冷劑的一項(xiàng)突破性技術(shù)[9]。當(dāng)前分子模擬方法主要有：蒙特卡洛（Monte Carlo，MC）、分子動(dòng)力學(xué)（Molecular Dynamics，MD）、分子力學(xué)（Molecular Mechanics，MM）和量子力學(xué)（Quantum Mechanics，QM）[10]。其中，蒙特卡洛法是隨機(jī)性的分子模擬方法，執(zhí)行了Gibbs的統(tǒng)計(jì)力學(xué)，還將分子動(dòng)力學(xué)方法建立在經(jīng)典力學(xué)算符理論的基礎(chǔ)上，使分子模擬有了一個(gè)嚴(yán)格的理論依據(jù)。重慶大學(xué)蔣國柱等[11]應(yīng)用該方法對(duì)L-J流體的氣液相平衡開展了模擬研究。

混合制冷劑氣液相平衡研究的方法主要有模擬仿真、狀態(tài)方程和實(shí)驗(yàn)。目前應(yīng)用最廣泛的氣液相平衡模擬方法[12]就是吉布斯系統(tǒng)蒙特卡洛方法（GibbsEnsemble Monte Carlo，GEMC），通過分別模擬氣相和液相避免了對(duì)相界面的直接處理，對(duì)純工質(zhì)及混合工質(zhì)的氣液相平衡模擬結(jié)果與實(shí)驗(yàn)具有較高的重合度。2009年，RAABE等[13]開發(fā)了一種適用于預(yù)測R1234yf熱力學(xué)性質(zhì)的新型力場，并用GEMC方法模擬了其氣液相平衡特性。目前，由于R32和R1234yf混合氣液相實(shí)驗(yàn)數(shù)據(jù)較少，本文分別采用GEMC法和PR狀態(tài)方程法對(duì)混合制冷劑R32/R1234yf進(jìn)行氣液相平衡模擬，并將MC模擬和狀態(tài)方程計(jì)算的結(jié)果進(jìn)行比較分析。

1 GEMC模擬

1.1 Monte Carlo基本原理

MC方法又稱為隨機(jī)取樣法，它以概率論中的大數(shù)定理和中心極限定理為基礎(chǔ)，通過不斷產(chǎn)生隨機(jī)數(shù)序列來模擬物體的運(yùn)動(dòng)過程。MC模擬方法通過粒子的隨機(jī)移動(dòng)使構(gòu)型發(fā)生改變[14]，進(jìn)行系綜上的平均，當(dāng)給初始構(gòu)型賦予體系后，對(duì)構(gòu)型內(nèi)每個(gè)粒子進(jìn)行隨機(jī)移動(dòng)，比較移動(dòng)前后體系狀態(tài)的變化，按照一定的接受準(zhǔn)則對(duì)移動(dòng)進(jìn)行判斷，若移動(dòng)滿足這一準(zhǔn)則就接受，若不滿足就拒絕，重復(fù)移動(dòng)多次，就能得到體系的真實(shí)構(gòu)型，通過對(duì)系綜上的平均求出體系的宏觀性質(zhì)。圖1為MC模擬流程圖。

圖1　MC模擬流程圖

1.1.1 系綜

系綜是許多同樣結(jié)構(gòu)并處在給定宏觀條件之下的系統(tǒng)集合[15]，其中每個(gè)系統(tǒng)各處于某一個(gè)微觀狀態(tài)，并且彼此獨(dú)立；在相空間內(nèi)微觀運(yùn)動(dòng)狀態(tài)組成一個(gè)連續(xù)的區(qū)域，宏觀量與微觀量相對(duì)應(yīng)，在一定的宏觀條件下是所有可能的運(yùn)動(dòng)狀態(tài)的平均值。在傳統(tǒng)MC模擬方法的應(yīng)用中，存在著兩種系統(tǒng)狀態(tài)的模擬：平衡態(tài)的模擬和非平衡態(tài)的模擬。在氣液相平衡的計(jì)算中，主要采用平衡態(tài)的模擬。平衡態(tài)系綜主要有以下幾個(gè)類型：微正則系綜（N、V、E）、正則系綜（N、V、T）、等溫等壓系綜（N、P、T）、巨正則系綜和吉布斯（Gibbs）系綜等。

1.1.2 力場

力場方法可以看成是勢能函數(shù)的經(jīng)驗(yàn)表達(dá)式，是MC模擬內(nèi)容的基礎(chǔ)。分子力場大致上可分為兩類：全原子力場和聯(lián)合原子力場。此外，也可根據(jù)選擇的函數(shù)和力場參數(shù)分為以下3類：傳統(tǒng)力場，如AMBER、CHARMM、CVFF和MM等系列力場；第二代力場，如CFF91、PCFF、CFF95和MMF94等系列力場；通用力場，如ESFF、UFF和DREDING等系列力場。

1.2 GEMC方法計(jì)算原理

GEMC方法能夠同時(shí)在氣相和液相兩個(gè)模擬盒子中進(jìn)行MC模擬，并且兩者保持相對(duì)獨(dú)立，但是滿足相平衡條件（溫度T、壓力P和化學(xué)勢μ相等），總粒子數(shù)N、總體積V和溫度保持不變；為了能夠滿足相平衡條件，在模擬過程中需要進(jìn)行3種不同的MC移動(dòng)，以不同的接受概率進(jìn)行判斷。

1）氣相和液相盒子中的粒子分別在盒子內(nèi)自由移動(dòng)，如轉(zhuǎn)動(dòng)、平動(dòng)等，使盒子達(dá)到平衡狀態(tài)，粒子移動(dòng)被接受的概率為：

2）在總體積不變的情況下，氣相和液相的體積發(fā)生變化，使兩個(gè)模擬盒子的壓力平衡，粒子移動(dòng)被接受的概率為：

3）在粒子總數(shù)N保持不變的情況下，在氣相和液相兩個(gè)模擬盒子中進(jìn)行交換，使兩個(gè)盒子的化學(xué)勢達(dá)到平衡狀態(tài)，粒子移動(dòng)被接受的概率為：

式中：

E——盒子的能量；

V——體積；

N——粒子數(shù)；

k——玻爾茲曼常數(shù)。

1.3 模擬方法

本文采用NPT+GEMC方法對(duì)R32/R1234yf進(jìn)行氣液相平衡的模擬，模擬建立氣相和液相兩個(gè)盒子，兩個(gè)盒子的初始結(jié)構(gòu)為面心立方（FCC）結(jié)構(gòu)。根據(jù)混合制冷劑的飽和蒸氣壓數(shù)據(jù)，可以確定在質(zhì)量比20％/80％時(shí)，液相分子數(shù)分別為20/44，氣相分子數(shù)分為為83/182；在質(zhì)量比25/75時(shí)，液相分子數(shù)分別為23/40，氣相分子數(shù)分別為82/150。所有的Monte Carlo計(jì)算均采用MAPS（Materials and Process Simulation）的Towhee模塊進(jìn)行。GEMC模擬采用通用UFF力場，電荷計(jì)算采用Gasteiger方法，系統(tǒng)總的電荷量為0，圖2為模擬所建立的氣相與液相盒子。

圖2　氣相與液相模擬盒子

模擬溫度范圍255 K～320 K，溫度間隔為3 K～4 K，初始設(shè)定的壓力P為溫度T對(duì)應(yīng)的飽和蒸汽壓，R32/R1234yf的飽和蒸氣壓數(shù)據(jù)采用KAMIAKA[16]的實(shí)驗(yàn)值。在模擬開始前需要對(duì)兩個(gè)盒子進(jìn)行優(yōu)化，優(yōu)化方法為最速下降法，運(yùn)行步數(shù)為1,000步，精度為0.01，van de Waals截?cái)嘁蜃訛?.5，花鍵移位因子為1.8。模擬循環(huán)的次數(shù)為10萬次，前5萬次循環(huán)用于使模擬體系達(dá)到氣液相平衡，后5萬次循環(huán)用于統(tǒng)計(jì)氣液相平衡數(shù)據(jù)并取平均值。構(gòu)型偏倚公式采用Martin和Frischknecht（2006）形式，構(gòu)型偏倚設(shè)置按照Martin和Frischknecht方法，加入靜電作用項(xiàng)，電介質(zhì)常數(shù)為1 Ang，靜電作用項(xiàng)的計(jì)算采用Ewald求和法。Van de Waals截?cái)喟霃饺橐合嗪虚L的一半，截?cái)嗪蟛捎梦膊啃拚椒ā?/p>

模擬中Monte Carlo的移動(dòng)形式為：分子質(zhì)量中心的平移、質(zhì)量中心的轉(zhuǎn)動(dòng)、分子內(nèi)單原子的平動(dòng)、兩個(gè)盒子分子的構(gòu)型偏倚平動(dòng)、兩個(gè)盒子分子的構(gòu)型偏倚轉(zhuǎn)動(dòng)、各向異性體積移動(dòng)、各向同性體積移動(dòng)。移動(dòng)概率見表1。

表1　R32/R1234yf的MC移動(dòng)概率

1.4 模擬結(jié)果

不同配比的混合制冷劑模擬氣液相平衡密度如表2和表3所示。

由于缺乏R32/R1234yf的氣液相平衡密度的直接實(shí)驗(yàn)數(shù)據(jù)，下節(jié)將與狀態(tài)方程在相同狀態(tài)下的計(jì)算結(jié)果進(jìn)行比較驗(yàn)證。

表2　R32/R1234yf(20/80,wt％)模擬氣液相平衡密度

表3　R32/R1234yf(50/50,wt％)模擬氣液相平衡密度

2 氣液相平衡特性計(jì)算及結(jié)果分析

目前，在NIST（National Institute of Standards and Technology）所發(fā)布的軟件Refprop9.0中還沒有R32/R1234yf的實(shí)驗(yàn)數(shù)據(jù)，并且由于混合工質(zhì)的種類和配比繁多，實(shí)驗(yàn)測定工作量大，實(shí)驗(yàn)時(shí)間長，不可能一一加以實(shí)驗(yàn)研究。另外，根據(jù)KIM等[17]的研究可知，PR狀態(tài)方程計(jì)算結(jié)果與實(shí)驗(yàn)數(shù)據(jù)有很好的一致性。因此本文將用參數(shù)較少的PR方程對(duì)R1234yf/R32的氣液相平衡性質(zhì)進(jìn)行研究。

2.1 PR狀態(tài)方程介紹

PR方程形式如下所示[18]：

式(4)可寫成式(5)的形式：

式(5)根據(jù)系統(tǒng)中的相數(shù)會(huì)產(chǎn)生一個(gè)或兩個(gè)根，在兩相區(qū)，最大的根是氣體的壓縮因子，最小的正根是液體的壓縮因子。式(4)應(yīng)用到臨界點(diǎn)，有：

在非臨界點(diǎn)有：

式中：

Tr——對(duì)比溫度，T/Tc，

ω——偏心因子。

α(Tr,ω)的形式如下：

2.2 混合法則

一般來說，狀態(tài)方程只適用于純工質(zhì)計(jì)算，混合工質(zhì)則必須采用一定的混合法則以保證計(jì)算的精度。常用的混合法則有Van De Waals（VDW）、Huron-Vidal（HV）、Stryjek-Vera（SV）和MHV1等。其中，SV[19]混合法則是VDW的改進(jìn)形式，比VDW具有較高精確度。但是當(dāng)二元工質(zhì)混合時(shí)，SV法則即可簡化為VDW法則。因此，本文采用VDW混合法則計(jì)算混合工質(zhì)R32/R1234yf的氣液相平衡特性。

VDW混合法則的形式如下：

式中kij為二元相互作用系數(shù)，xi為組分i的摩爾分?jǐn)?shù)，ai和bi分別為組分i狀態(tài)方程參數(shù)。

在混合法則中，二元相互作用系數(shù)kij有著重要的影響。kij表示混合物的兩種分子間的非理想作用的特性和相互作用性質(zhì)[20]，kij的準(zhǔn)確度對(duì)混合物氣液相平衡的計(jì)算至關(guān)重要，kij一般是由氣液相平衡數(shù)據(jù)優(yōu)化回歸得到。這里，R32/R1234yf的二元相互作用系數(shù)采用文獻(xiàn)[16]的實(shí)驗(yàn)數(shù)據(jù)擬合值0.037。

2.3 計(jì)算結(jié)果及分析

不同配比的混合二元制冷劑的飽和氣液相密度如表4和表5。由圖3和4可知，用GEMC模擬與PR狀態(tài)方程計(jì)算的R32/R1234yf飽和氣相密度相對(duì)偏差在-3.0％～+2.0％范圍內(nèi)，相對(duì)偏差絕對(duì)值的平均偏差1.44％。飽和液相密度相對(duì)偏差在-4.5％～+3.5％范圍內(nèi)，相對(duì)偏差絕對(duì)值的平均偏差2.25％。

表4　R32/R1234yf(20/80,wt％)飽和氣液相密度

表5　R32/R1234yf(50/50,wt％)飽和氣液相密度

圖3　R32/R1234yf的MC飽和氣相密度與PR數(shù)據(jù)的偏差

圖4　R32/R1234yf的MC飽和液相密度與PR數(shù)據(jù)的偏差

3 結(jié)論

GEMC模擬和PR狀態(tài)方程對(duì)于氣液相平衡性質(zhì)的計(jì)算都具有相當(dāng)高的精度，本文選擇了適當(dāng)?shù)膮?shù)對(duì)二元混合制冷劑R32/R1234yf分別進(jìn)行了氣液相平衡計(jì)算。結(jié)果表明，GEMC模擬方法誤差小，為進(jìn)一步研究混合制冷劑R32/R1234yf氣液相平衡特性奠定了基礎(chǔ)。

參考文獻(xiàn)：

[1]陳斌,陳光明.R407C、R410A制冷系統(tǒng)相關(guān)特性研究進(jìn)展[J].制冷與空調(diào),2003,31(2):42-45.

[2]汪軍,李美玲,周文濤.核電站空調(diào)系統(tǒng)HCFC22冷水機(jī)組替代工質(zhì)性質(zhì)的研究[J].上海理工大學(xué)學(xué)報(bào),2001,23(1):8-13.

[3]葛志松,吳獻(xiàn)忠,崔曉鈺,等.三元非共沸混合工質(zhì)變組成容量控制空調(diào)系統(tǒng)可行性分析[J].上海理工大學(xué)學(xué)報(bào),2005,27(4):336-340.

[4]LONGO G A,ZILIO C,RIGHETTI G.Condensation of the low GWP refrigerant HFC152a inside a brazed plate heat exchanger[J].Experimental Thermal and Fluid Science,2015,68:509-515.

[5]陸岷山,祁影霞,楊喜,等.混合工質(zhì)R134a/R600a 氣相PVTx性質(zhì)的實(shí)驗(yàn)研究[J].制冷技術(shù),2014,34(6):27-31,36.

[6]鄭賢德.制冷原理與裝置[M].北京:機(jī)械工業(yè)出版社,2010.

[7]劉圣春,饒志明,楊旭凱,等.新型制冷劑R1234yf的性能分析[J].制冷技術(shù),2013,33(1):56-59.

[8]肖學(xué)智,周曉芳,徐浩陽,等.低GWP制冷劑研究現(xiàn)狀綜述[J].制冷技術(shù),2014,34(6):37-42.

[9]GLASS C W,REISER S,RUTKAI G,et al.ms2:A molecular simulation tool for thermodynamic properties,new version release[J].Computer Physics Communications,2014,185(12):3302-3306.

[10]莊昌青,岳紅,張慧軍.分子模擬方法及模擬軟件Material Studio在高分子材料中的應(yīng)用[J].塑料,2010,39(4):81-84.

[11]蔣國柱,薛榮書,魏順安.吉布斯系綜Monte Carlo技術(shù)模擬L-J流體汽液相平衡[J].重慶大學(xué)學(xué)報(bào),2003,26(8):64-67.

[12]李曉峰,趙立峰,孫淮.GEMC和GDI方法計(jì)算流體氣液相平衡的比較[J].物理化學(xué)學(xué)報(bào),2008,24(10):1824-1830.

[13]RAABE G,MAGINN E J.Molecular modeling of the vapor-liquid equilibrium properties of the alternative refrigerant 2,3,3,3-tetrafluoro-1-propene (HFO-1234yf)[J].Journal of Physical Chemistry Letters,2010,1(1):93-96.

[14]胡慧娟.吉布斯系綜模擬氮?dú)?水和氧氣-水的氣液相平衡[D].南京:南京工業(yè)大學(xué),2005.

[15]陳舜麟.計(jì)算材料科學(xué)[M].北京:化學(xué)工業(yè)出版社,2005.

[16]KAMIAKA T,DANG C,HIHARA E.Vapor-liquid equilibrium measurements for binary mixtures of R1234yf with R32,R125,and R134a[J].International Journal of Refrigeration,2013,36(3):965-971.

[17]JU H K,MIN S K,KIM Y.Vapor-liquid equilibria for pentafluoroethane + propane and difluoromethane + propane systems over a temperature range from 253.15 to 323.15 K[J].Fluid Phase Equilibria,2003,211(2):273-287.

[18]PENG D Y,ROBINSON D B.A new two-constant equation of state[J].Industrial and Engineering Chemistry Fundamentals,1976,15(1):59-64.

[19]STRYJEK R,VERA J H.PRSV2:A cubic equation of state for accurate vapor-liquid equilibria calculations[J].The Canadian Journal of Chemical Engineering,1986,64:820-826.

[20]胡芃,高靜軒,陳龍祥,等.PR方程結(jié)合vdw混合法則推算二元及三元HFC/HC混合工質(zhì)氣液相平衡性質(zhì)[J].化工學(xué)報(bào),2012,63(2):350-355.

Investigation on Vapor-liquid Equilibrium Characteristics of HFO-1234yf/HFC-32

WU Dong*,QI Ying-xia,YANG Xi
(University of Shanghai for Science and Technology,Shanghai,200093)

[Abstract]Gibbs Ensemble Monte Carlo (GEMC) simulation is briefly introduced and adopted to model the vapor-liquid equilibrium characteristics of binary refrigerant HFO-1234yf/HFC-32 with mass fractions of 50％/50％ and 80％/20％ in the range of (273～333) K.The comparison between the simulated values and those calculated by PR state equations shows that,the relative deviations of the saturated vapor-liquid densites are –3％～+2％ and –4.5％～+3.5％,respectively.It shows that,the deviation of the GEMC simulation method is small,which provides the basis for further researching the liquid-vapor equilibrium characteristics for the mixed refrigerant of HFO-1234yf/HFC-32.

[Keywords]Simulation; HFO-1234yf/HFC-32; Vapor-liquid equilibrium; Equation of state

基金項(xiàng)目：上海市教育委員會(huì)科研創(chuàng)新項(xiàng)目（No.11YZ119）。

doi：10.3969/j.issn.2095-4468.2016.01.106