吴　东，祁影霞，杨　喜

（上海理工大学，上海 200093）



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

吴东*，祁影霞，杨喜

（上海理工大学，上海 200093）

[摘 要]本文简单介绍了吉布斯蒙特卡洛（GEMC）模拟法原理并采用GEMC模拟了混合制冷剂HFO-1234yf/HFC-32质量分数分别为50％/50％和80％/20％，在温度(273～333) K的气液相平衡特性。本文计算结果与PR状态方程计算值的比较表明，饱和气液相密度相对偏差分别为-3％～+2％、-4.5％～+3.5％。结果表明，GEMC模拟方法误差小，为进一步研究混合制冷剂HFO-1234yf/HFC-32气液相平衡特性奠定了基础。

[关键字]模拟；HFO-1234yf/HFC-32；气液相平衡；状态方程

*吴东（1988-），男，硕士在读，主要从事制冷剂热物性研究。通信作者：祁影霞，女，副教授，地址：上海市杨浦区军工路516号上海理工大学。E-mail：qipeggy@126.com。

0 引言

随着CFCs和HCFCs类制冷剂逐步的禁止使用，国际上形成了两种不同的制冷剂替代研究路线[1]：一条以美日为代表，主要采用合成HFCs类制冷剂，如R134a、R410A和R407C等；另外一条以德国和北欧一些国家为代表，主张使用天然制冷剂，如CO2、H2O和HCs等。当前，汽车空调依然广泛采用R134a作为制冷剂，主要是基于其ODP为0，安全性能高；根据欧盟F-gas法规，自2011年起禁止使用GWP大于150的制冷剂，而其GWP为1,430[2-4]。陆岷山等[5]通过实验的方法对混合制冷剂R134a/R600a的PVT性质进行了实验研究。

从上可以看出，由于理想制冷剂应符合热力学性质、迁移性质、物理化学性质和其他性质等诸多方面的要求，因此新型混合制冷剂的研究将有助于找到更为理想的替代制冷剂[6-8]。近年，由于计算机技术的巨大进步，建立在量子力学和统计热力学基础上的分子模拟方法已成为预测新型混合制冷剂的一项突破性技术[9]。当前分子模拟方法主要有：蒙特卡洛（Monte Carlo，MC）、分子动力学（Molecular Dynamics，MD）、分子力学（Molecular Mechanics，MM）和量子力学（Quantum Mechanics，QM）[10]。其中，蒙特卡洛法是随机性的分子模拟方法，执行了Gibbs的统计力学，还将分子动力学方法建立在经典力学算符理论的基础上，使分子模拟有了一个严格的理论依据。重庆大学蒋国柱等[11]应用该方法对L-J流体的气液相平衡开展了模拟研究。

混合制冷剂气液相平衡研究的方法主要有模拟仿真、状态方程和实验。目前应用最广泛的气液相平衡模拟方法[12]就是吉布斯系统蒙特卡洛方法（GibbsEnsemble Monte Carlo，GEMC），通过分别模拟气相和液相避免了对相界面的直接处理，对纯工质及混合工质的气液相平衡模拟结果与实验具有较高的重合度。2009年，RAABE等[13]开发了一种适用于预测R1234yf热力学性质的新型力场，并用GEMC方法模拟了其气液相平衡特性。目前，由于R32和R1234yf混合气液相实验数据较少，本文分别采用GEMC法和PR状态方程法对混合制冷剂R32/R1234yf进行气液相平衡模拟，并将MC模拟和状态方程计算的结果进行比较分析。

1 GEMC模拟

1.1 Monte Carlo基本原理

MC方法又称为随机取样法，它以概率论中的大数定理和中心极限定理为基础，通过不断产生随机数序列来模拟物体的运动过程。MC模拟方法通过粒子的随机移动使构型发生改变[14]，进行系综上的平均，当给初始构型赋予体系后，对构型内每个粒子进行随机移动，比较移动前后体系状态的变化，按照一定的接受准则对移动进行判断，若移动满足这一准则就接受，若不满足就拒绝，重复移动多次，就能得到体系的真实构型，通过对系综上的平均求出体系的宏观性质。图1为MC模拟流程图。

图1　MC模拟流程图

1.1.1 系综

系综是许多同样结构并处在给定宏观条件之下的系统集合[15]，其中每个系统各处于某一个微观状态，并且彼此独立；在相空间内微观运动状态组成一个连续的区域，宏观量与微观量相对应，在一定的宏观条件下是所有可能的运动状态的平均值。在传统MC模拟方法的应用中，存在着两种系统状态的模拟：平衡态的模拟和非平衡态的模拟。在气液相平衡的计算中，主要采用平衡态的模拟。平衡态系综主要有以下几个类型：微正则系综（N、V、E）、正则系综（N、V、T）、等温等压系综（N、P、T）、巨正则系综和吉布斯（Gibbs）系综等。

1.1.2 力场

力场方法可以看成是势能函数的经验表达式，是MC模拟内容的基础。分子力场大致上可分为两类：全原子力场和联合原子力场。此外，也可根据选择的函数和力场参数分为以下3类：传统力场，如AMBER、CHARMM、CVFF和MM等系列力场；第二代力场，如CFF91、PCFF、CFF95和MMF94等系列力场；通用力场，如ESFF、UFF和DREDING等系列力场。

1.2 GEMC方法计算原理

GEMC方法能够同时在气相和液相两个模拟盒子中进行MC模拟，并且两者保持相对独立，但是满足相平衡条件（温度T、压力P和化学势µ相等），总粒子数N、总体积V和温度保持不变；为了能够满足相平衡条件，在模拟过程中需要进行3种不同的MC移动，以不同的接受概率进行判断。

1）气相和液相盒子中的粒子分别在盒子内自由移动，如转动、平动等，使盒子达到平衡状态，粒子移动被接受的概率为：

2）在总体积不变的情况下，气相和液相的体积发生变化，使两个模拟盒子的压力平衡，粒子移动被接受的概率为：

3）在粒子总数N保持不变的情况下，在气相和液相两个模拟盒子中进行交换，使两个盒子的化学势达到平衡状态，粒子移动被接受的概率为：

式中：

E——盒子的能量；

V——体积；

N——粒子数；

k——玻尔兹曼常数。

1.3 模拟方法

本文采用NPT+GEMC方法对R32/R1234yf进行气液相平衡的模拟，模拟建立气相和液相两个盒子，两个盒子的初始结构为面心立方（FCC）结构。根据混合制冷剂的饱和蒸气压数据，可以确定在质量比20％/80％时，液相分子数分别为20/44，气相分子数分为为83/182；在质量比25/75时，液相分子数分别为23/40，气相分子数分别为82/150。所有的Monte Carlo计算均采用MAPS（Materials and Process Simulation）的Towhee模块进行。GEMC模拟采用通用UFF力场，电荷计算采用Gasteiger方法，系统总的电荷量为0，图2为模拟所建立的气相与液相盒子。

图2　气相与液相模拟盒子

模拟温度范围255 K～320 K，温度间隔为3 K～4 K，初始设定的压力P为温度T对应的饱和蒸汽压，R32/R1234yf的饱和蒸气压数据采用KAMIAKA[16]的实验值。在模拟开始前需要对两个盒子进行优化，优化方法为最速下降法，运行步数为1,000步，精度为0.01，van de Waals截断因子为2.5，花键移位因子为1.8。模拟循环的次数为10万次，前5万次循环用于使模拟体系达到气液相平衡，后5万次循环用于统计气液相平衡数据并取平均值。构型偏倚公式采用Martin和Frischknecht（2006）形式，构型偏倚设置按照Martin和Frischknecht方法，加入静电作用项，电介质常数为1 Ang，静电作用项的计算采用Ewald求和法。Van de Waals截断半径取为液相盒长的一半，截断后采用尾部修正方法。

模拟中Monte Carlo的移动形式为：分子质量中心的平移、质量中心的转动、分子内单原子的平动、两个盒子分子的构型偏倚平动、两个盒子分子的构型偏倚转动、各向异性体积移动、各向同性体积移动。移动概率见表1。

表1　R32/R1234yf的MC移动概率

1.4 模拟结果

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

由于缺乏R32/R1234yf的气液相平衡密度的直接实验数据，下节将与状态方程在相同状态下的计算结果进行比较验证。

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

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

2 气液相平衡特性计算及结果分析

目前，在NIST（National Institute of Standards and Technology）所发布的软件Refprop9.0中还没有R32/R1234yf的实验数据，并且由于混合工质的种类和配比繁多，实验测定工作量大，实验时间长，不可能一一加以实验研究。另外，根据KIM等[17]的研究可知，PR状态方程计算结果与实验数据有很好的一致性。因此本文将用参数较少的PR方程对R1234yf/R32的气液相平衡性质进行研究。

2.1 PR状态方程介绍

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

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

式(5)根据系统中的相数会产生一个或两个根，在两相区，最大的根是气体的压缩因子，最小的正根是液体的压缩因子。式(4)应用到临界点，有：

在非临界点有：

式中：

Tr——对比温度，T/Tc，

ω——偏心因子。

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

2.2 混合法则

一般来说，状态方程只适用于纯工质计算，混合工质则必须采用一定的混合法则以保证计算的精度。常用的混合法则有Van De Waals（VDW）、Huron-Vidal（HV）、Stryjek-Vera（SV）和MHV1等。其中，SV[19]混合法则是VDW的改进形式，比VDW具有较高精确度。但是当二元工质混合时，SV法则即可简化为VDW法则。因此，本文采用VDW混合法则计算混合工质R32/R1234yf的气液相平衡特性。

VDW混合法则的形式如下：

式中kij为二元相互作用系数，xi为组分i的摩尔分数，ai和bi分别为组分i状态方程参数。

在混合法则中，二元相互作用系数kij有着重要的影响。kij表示混合物的两种分子间的非理想作用的特性和相互作用性质[20]，kij的准确度对混合物气液相平衡的计算至关重要，kij一般是由气液相平衡数据优化回归得到。这里，R32/R1234yf的二元相互作用系数采用文献[16]的实验数据拟合值0.037。

2.3 计算结果及分析

不同配比的混合二元制冷剂的饱和气液相密度如表4和表5。由图3和4可知，用GEMC模拟与PR状态方程计算的R32/R1234yf饱和气相密度相对偏差在-3.0％～+2.0％范围内，相对偏差绝对值的平均偏差1.44％。饱和液相密度相对偏差在-4.5％～+3.5％范围内，相对偏差绝对值的平均偏差2.25％。

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

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

图3　R32/R1234yf的MC饱和气相密度与PR数据的偏差

图4　R32/R1234yf的MC饱和液相密度与PR数据的偏差

3 结论

GEMC模拟和PR状态方程对于气液相平衡性质的计算都具有相当高的精度，本文选择了适当的参数对二元混合制冷剂R32/R1234yf分别进行了气液相平衡计算。结果表明，GEMC模拟方法误差小，为进一步研究混合制冷剂R32/R1234yf气液相平衡特性奠定了基础。

参考文献：

[1]陈斌,陈光明.R407C、R410A制冷系统相关特性研究进展[J].制冷与空调,2003,31(2):42-45.

[2]汪军,李美玲,周文涛.核电站空调系统HCFC22冷水机组替代工质性质的研究[J].上海理工大学学报,2001,23(1):8-13.

[3]葛志松,吴献忠,崔晓钰,等.三元非共沸混合工质变组成容量控制空调系统可行性分析[J].上海理工大学学报,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]陆岷山,祁影霞,杨喜,等.混合工质R134a/R600a 气相PVTx性质的实验研究[J].制冷技术,2014,34(6):27-31,36.

[6]郑贤德.制冷原理与装置[M].北京:机械工业出版社,2010.

[7]刘圣春,饶志明,杨旭凯,等.新型制冷剂R1234yf的性能分析[J].制冷技术,2013,33(1):56-59.

[8]肖学智,周晓芳,徐浩阳,等.低GWP制冷剂研究现状综述[J].制冷技术,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在高分子材料中的应用[J].塑料,2010,39(4):81-84.

[11]蒋国柱,薛荣书,魏顺安.吉布斯系综Monte Carlo技术模拟L-J流体汽液相平衡[J].重庆大学学报,2003,26(8):64-67.

[12]李晓峰,赵立峰,孙淮.GEMC和GDI方法计算流体气液相平衡的比较[J].物理化学学报,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].南京:南京工业大学,2005.

[15]陈舜麟.计算材料科学[M].北京:化学工业出版社,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方程结合vdw混合法则推算二元及三元HFC/HC混合工质气液相平衡性质[J].化工学报,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

基金项目：上海市教育委员会科研创新项目（No.11YZ119）。

doi：10.3969/j.issn.2095-4468.2016.01.106