溴化锂水溶液汽液相平衡的分子模拟

周玉红，赵宗昌，张晓冬，唐静，辛绪亮

(大连理工大学 化工与环境生命学部，辽宁 大连 116023)

吸收式制冷和吸收式热泵是回收和利用工业余热最直接和有效的技术装备，是节能减排，实现低碳经济最有效的技术手段之一.吸收循环热力学效率与其工作流体的热力学性质有关.目前，工业上的吸收制冷机和热泵机组主要采用溴化锂水溶液为工作流体，具有广阔的市场应用前景.溴化锂水溶液的汽液相平衡性质以及比热容是重要的热力学物性参数，对于吸收循环的开发与设计具有重要意义.

计算机分子模拟通过研究体系的微观性质而得到体系的宏观热力学性质，其中Gibbs系综蒙特卡洛方法(Gibbs ensemble Monte Carlo，GEMC)在预测混合物的汽液相平衡数据方面有独特的优势.目前该方法的应用主要集中在模拟有机物的一元流体或二元溶液的汽液相平衡组成［1-2］和气体在溶液中的溶解度方面［3］.关于溴化锂水溶液汽液相平衡性质的分子模拟尚未见报道.

本文采用 NPT-GEMC［4-5］，用 MCCCS Towhee［6］软件包进行LiBr/H2O二元组分的汽液相平衡性质的模拟;分析了不同压力、浓度下，液相的构型能变化规律和溶液的微观结构的变化规律;最后模拟了溴化锂溶液的恒压比热容.

1 势能模型

在溴化锂水溶液中存在3种分子间力:1)水分子与水分子间的作用力;2)离子与离子间的作用力;3)离子与水之间的作用力.描述粒子间相互作用力的一套参数化的经验势函数，称为势能模型.因为溶液的热力学性质主要由分子间作用力决定，故本文中的势能模型函数Uext包含范德华作用项和静电作用项，不包含分子内作用项，如下:

其中范德华作用项用Lennard-Jones势能模型描述，静电作用项用固定点电荷库仑势描述.式中:εij、σij分别为粒子i和j范德华作用的势阱和相互作用半径;qi和qj分别为粒子i和j的电荷数;rij为粒子i和j的距离;Uext(rij)为粒子i和j的分子间势能.

水的势能模型中常用的是固定点电荷模型［7-8］SPC/E(simple point charge/extended)，该模型在分子模拟溴化锂水溶液的热力学性质中应用相对广泛［9-10］.因此水分子间作用选用该模型描述势能参数(见表1).溴离子和锂离子势能模型参数分别取自分子模拟软件MCCCS Towhee中Lybrand Ghosh McCammon 力场［11］和 Aqvist力场［12］(见表 1).不同粒子之间的Lennard-Jones势能参数用Lorentz-Berthelot混合规则.

表1 离子和水的势能参数Table 1 Potential parameters of water and ions

2 模拟实验

在Linux操作系统下，计算均采用MCCCS Towhee 6.2.11 软件包［6］，运用 NPT-GEMC 系综方法进行模拟.由于工业上吸收制冷和吸收热泵所用LiBr水溶液的质量浓度范围为40% ～60%，故选取常用操作温度为80℃、100℃下，所对应的平衡压力进行模拟.模拟采用周期性边界条件，范德华作用采用截断半径加尾部矫正方法计算，静电长程相互作用采用Ewald方法计算.为提高计算效率，采用构型偏倚蒙特卡罗方法［13］以提高粒子插入的接收概率;采用聚集体积偏倚蒙特卡罗方法［14］以提高聚集流体的构型取样速率.

每次模拟进行30×104次循环，其中前15×104次循环对体系进行预平衡，后15×104次循环用于热力学性质的统计和计算，每次统计体系的构型数不低于107.每次循环包括N次MC移动，其中N为模拟体系中的粒子总数.每次移动会选择6种MC移动之一，它们分别为:体积变化、构型偏倚交换粒子、体积聚集偏倚、分子构型偏倚重新生长、粒子质心平动、粒子质心转动.

3 模拟结果

3.1 溴化锂水溶液平衡组成

溴化锂水溶液的液相平衡组成的分子模拟结果见表2.

表2 t=80℃、100℃时溴化锂水溶液平衡组成Table 2 The equilibrium concentrations for aqueous solution of lithium bromide at 80℃ and 100℃

表2中xcal为分子模拟计算得到的液相溴化锂质量分数，xref为 ASHRAE手册［15］中溴化锂相平衡P-T-x公式，迭代得到的文献值.用Δ表示模拟结果和文献值的相对偏差:

用ARD表示平均相对偏差:

通过表2可知，相对平均偏差分别9.67%、11.31%，表明NPT-GEMC方法可以较好地模拟溴化锂水溶液平衡组成关系P-T-x.将表2中数据作图比较，见图1.比较图1(a)、(b)可知，在溴化锂浓度接近于55%左右时，模拟值和文献值非常接近;当浓度在55%～70%时，模拟值会高于文献值;当浓度在40%～55%时，模拟值低于文献值.该结果表明所选水和离子的势能模型还不能在所有浓度范围内完全准确的模拟平衡时的溴化锂溶液的相平衡关系P-T-x，但是可以应用于工业上常用的浓度范围(40% ～60%).

图1 不同温度下，不同压力下的溴化锂水溶液平衡组成Fig.1 The equilibrium concentrations for aqueous solution of lithium bromide at different temperatures and pressures

3.2 溴化锂水溶液平衡密度

溴化锂水溶液的汽液相平衡时，液相密度的分子模拟结果见表3.表中ρcal为分子模拟计算得到的不同温度压力下的液相溴化锂溶液的密度;根据文献［16］提供的溴化锂溶液密度随温度t和浓度x的变化曲线图及数据，运用3.1中xref及对应温度得到ρref.用Δ表示模拟结果和文献值的相对偏差.

表3 t=80℃、100℃时溴化锂水溶液平衡时液相密度Table 3 The equilibrium densities for aqueous solution of lithium bromide at 80℃ and 100℃

通过表3可知，平均相对偏差分别为8.65%、8.33%，表明NPT-GEMC方法可以较好地模拟溴化锂水溶液平衡密度关系P-T-ρ.将表3中数据作图，见图2.比较图2(a)、(b)可知，分子模拟值ρcal在所模拟的浓度范围内均小于文献值ρref.其原因是:1)3.1节中汽液平衡组成模拟的误差，直接导致了平衡密度的不准确;2)水和离子的力场模型不是足够的精准，有待于进一步改进.

图2 t=100℃时不同压力下的溴化锂水溶液平衡密度Fig.2 The equilibrium densities for lithium bromide aqueous solution at different temperature pressures

3.3 溴化锂水溶液平衡时的构型能分析

溴化锂水溶液达到平衡时，分子间总构型能Uext、静电作用势能UC以及范德华作用ULJ见表4.其中由式(1)可知:Uext=UC+ULJ.

表4 t=80℃、100℃时溴化锂水溶液的构型能Table 4 The configurational energies for aqueous solution of lithium bromide at 80℃ and 100℃

从表4可以看出，在汽液平衡的溴化锂溶液中:1)分子间的总的构型能Uext主要由静电相互作用能Uc组成，分子间的范德华作用势能ULJ相对较弱;2)静电作用能均为负，表现为吸引力;范德华作用能均为正，表现为排斥力;总的分子间作用能均为负;3)同一温度下随着压力的增加，溴化锂饱和溶液浓度降低，分子间的静电作用能增大，而范德华相互作用能降低;但分子间的总的构型能逐渐增加.

3.4 溴化锂水溶液平衡时的微观结构分析

已有研究［17］主要采用分子动力学方法(MD)研究几个阴(阳)离子存在于水中的水合作用，并未涉及阴阳离子对同时在水中时的溶液结构，不能反映浓溶液的真实水合作用.本文在模拟体系达到平衡时，统计了溶液中不同粒子的径向分布函数:g(Li+-O)、g(Li+-H)、g(Br--O)以及 g(Br--H)，如图3、4 所示.

从图3、4可以看出不同浓度下各粒子间的径向分布函数具有大致相同的趋势.图3中Li+-O的径向分布函数有2个极大值，即2个明显的配位圈;Li+-H的径向分布函数有2个配位圈;g(Li+-O)的第1个峰位在2.033～2.095˚A之间，而g(Li+-H)的第1个峰位在2.675～2.755˚A之间.这表明溶液中的水分子以带负电的氧原子靠近Li+.Li+与水具有较强的水合作用.图4中Br--H的径向分布函数有3个明显的配位圈，Br--O的径向分布函数有2个配位圈，g(Br--H)的第 1个峰位在2.305～2.425˚A之间，g(Br--O)的第1个峰位在3.325～3.485˚A之间.这表明溶液中的水分子以带正电的氢原子靠近Br-.Br-与水也具有较强的水合作用.

粒子i周围r处球体内j粒子的配位数nij的定义为

式中:ρj为j粒子的数密度.本文根据gij(r)求出第一配位圈中阴、阳离子与氧原子、氢原子的配位数nij(r)，见表5.由于离子与氢原子的径向分布函数统计的只是离子周围1号氢原子的分布，故表5中离子周围氢原子的配位数扩大2倍.

图3 t=100℃时不同浓度下锂离子与水的径向分布函数Fig.3 The radial distribution functions for Li+and water at 100℃and different concentrations

图4 t=100℃时不同浓度下溴离子与水的径向分布函数Fig.4 The radial distribution functions for Br-and water at 100℃and different concentrations

表5 t=100℃时不同浓度下的溴化锂水溶液的径向分布函数的第一峰值、位置和配位数Table 5 The first peak values and positions of radial distribution functions and coordination number of lithium bromide aqueous at 100℃and different concentrations

从表5中可以看出随着溶液浓度的降低，阴、阳离子与水的径向分布函数的第1峰值均呈现出降低的趋势，但是第1配位圈内的配位数均呈现增加的趋势;而阴离子与阳离子的径向分布函数的第1峰值呈现增加的趋势，但第1配位圈内的配位数呈现降低的趋势.这表明溴化锂水溶液浓度降低，阴阳离子间的作用减弱，离子与水的水合作用增强.这是由于离子的水合作用受到异号离子的影响，溶液越稀影响越小，因而离子的水合作用越强，当达到无限稀释溶液时，离子周围的配位数最大，如Li+周围水的配位数最大值约为4.4［17］.从表5中还可以看出溶液越稀，水分子之间的径向分布函数的第1峰值和第1配位圈内的配位数均增加，这表明水溶液中水分子之间的作用受到离子的影响，离子的存在减弱了水分子之间的作用.

3.5 溴化锂水溶液恒压比热容的模拟

式中:H为焓.上述MC方法只能计算出体系的分子间势能，如果不考虑动能和分子内势能，可以通过构型焓Hcon计算比热容，构型焓的定义如下:

Kazuhiko Murakami［10］用构型焓计算了溴化锂溶液的比热容，其计算的大部分数据误差在-20%～-40%左右.本文采用该方法进行了计算，计算结果CPcon及相对偏差Δcon见表6.

本文用一种更精确的方法对溴化锂水溶液的比热容进行了模拟.将焓值分为理想气体焓Hid和剩余焓Hres:

式中:Uint表示分子内势能;K表示动能.故热容的计算公式为

理想气体热容可以通过实验数据关联、基团贡献法［18］和量化计算方法得到［19-20］.H20 的理想气体热容采用量化计算的方法，用 Gaussian 03软件包［21］，在 B3LYP/6-31+G＊＊水平下进行构型优化和频率分析，考虑了水分子的平动能和转动能以及振动能，得到了298.15K，1个标题的大气压下水的理想气体热容.本文计算的是373.15K、353.15K的比热容，因此需要按一定的比例矫正.水的校正因子为 0.964 2［22］，故水的理想气体热容分别为:34.03 J·mol-1·K-1、33.87 mol-1K-1.而对于 Li+、Br-

恒压热容的计算公式为单个离子，只需要考虑平动能，因而没有必要进行量子化学的计算.温度不影响单原子离子的理想气体热容，也不需要矫正.根据实验数据［22］知2种离子的理想气体热容均为:20.786 J·mol-1·K-1.对于不同浓度的溴化锂溶液可采用公式计算其理想气体比热容［18］:

剩余热容可以在NPT系综下用MC方法模拟得到:

由于溴化锂水溶液的热容随温度的变化不大，将上式简化为

式中:T1、T2为T附近的2个温度.取 T1=T+10，T2=T-10，并将计算出的恒压热容换算成为恒压比热容.模拟计算出的比热容数据CP以及文献关联数据CrefP和相对偏差Δ见表6.

100℃和80℃的平均相对偏差分别为:5.65%、11.18%，而构型焓方法的平均相对偏差分别为39.08%、30.05%，可以看出该方法可以较好地模拟溴化锂溶液的恒压比热容.

表6 t=80℃、100℃时不同浓度的溴化锂水溶液的恒压比热容Table 6 The constant pressure heat capacity of aqueous solution of lithium bromide at 100℃ and 80℃

4 结论

1)NPT-GEMC方法可以较好的模拟溴化锂水溶液平衡组成关系 P-T-x，相对平均偏差分别9.67%、11.31%.力场模型不能在所有浓度范围内完全准确的模拟平衡时的溴化锂溶液的相平衡关系P-T-x，但可以应用于工业上常用的浓度范围(40% ～60%).

2)NPT-GEMC方法可以较好的模拟溴化锂水溶液平衡密度组成关系P-T-ρ，相对平均偏差分别8.65%、8.33%.分子模拟值 ρcal在所模拟的浓度范围内均小于文献值ρref.

3)静电作用能是构型能的主要作用项，随着溶液浓度的降低，总构型能和静电作用均逐渐增加.

4)不同浓度下各粒子间的径向分布函数具有大致相同的趋势;阴阳离子与水均具有较强的水合作用;溴化锂水溶液越稀，离子周围第一配位圈内水的配位数增加，水合作用增强.

5)将溴化锂溶液的热容分为两部分，用量子化学的方法计算理想气体的热容，用NPT系综MC方法计算剩余热容.该法可以较好的模拟溴化锂溶液的恒压比热容，100℃和80℃的平均相对偏差分别为:5.65%、11.18%.

［1］李晓锋，赵立峰，孙淮.GEMC和GDI方法计算流体气液相平衡的比较流［J］.物理化学学报，2008，24(10):1824-1830.

LI Xiaofeng，ZHAO Lifeng，SUN Huai.Comparison of GEMC and GDI methods in prediction of the fluid vapor-liquid equilibrium［J］.Acta Phys Chim Sin，2008，24(10):1824-1830.

［2］蒲亮，缪强，张志炳，等.醋酸水溶液汽液非接触相平衡的蒙特卡罗模拟研究［J］.计算机与应用化学，2008，25(2):185-188.

PU Liang，MIAO Qiang，ZHANG Zhibing，et al .Monte Carlo simulation on the vapor-liquid equilibrium of the aqueous acetic acid solution［J］.Computers and Applied Chemistry，2008，25(2):185-188.

［3］VORHOLZ J，HARISMIADIS V I，PANAGIOTOPOULOS A Z，et al.Molecular simulation of the solubility of carbon dioxide in aqueous solutions of sodium chloride［J］.Fluid Phase Equilibria，2004，226:237-250.

［4］PANAGIOTOPOULOS A Z.Direct determination of phase coexistence properties of fluids by Monte Carlo simulation in a new ensemble ［J］.Molecular Physics，1987，61(4):813-826.

［5］PANAGIOTOPOULOS A Z.Phase equilibria by simulation in the Gibbs ensemble alternative derivation generalization and application to mixture and membrane equilibria［J］.Molecular Physics，1988，63(4):527-545.

［6］MARTIN M.MCCCS Towhee 6.2.11［EB/OL］.［2010-08-13］.Http://towhee.sourceforge.net/.

［7］GABRIELE R，SADUS R J.Influence of bond flexibility on the vapor-liquid phase equilibria of water［J］.Journal of Chemical Physics，2007，126:044701.

［8］VORHOLZ J，HARISMIADIS V I，RUMPF B，et al.Vapor+liquid equilibrium of water，carbon dioxide，and the binary system，water+carbon dioxide，from molecular simulation［J］.Fluid Phase Equilibria，2000，170:203-234.

［9］HIROFUMI D.Molecular dynamics study of nal-cohols adsorbed on an aqueous electrolyte solution［J］.Journal of Chemical Physics，2001，115(3):1538-1549.

［10］KAZUHIKO M，SHINYA I，MASAHIRO O.Thermodynamic properties of lithium-bromide aqueous by Monte Carlo simulation［J］.The Japan Society of Mechanical Engineers，2002，68:2862-2869.

［11］LYBRAND T P，INDIRA G，MCCAMMON J A，et al.Hydration of Chloride and Bromide An-ions:determination of relative free energy by computer simulation［J］.journal of the american chemical society，1985，107:7193-7794.

［12］AQVIST J.Ion-water interaction potentials derived from free energy perturbation simulations［J］.The Journal of Physics Chemistry，1990，94:8021-8024.

［13］MARTIN M G，Siepmann J I.Novel configurational-bias Monte Carlo method for branched molecules.Transferable potentials for phase equilibria.2.united atom description of branched alkanes［J］.The Journal of Physics Chemistry B，1999，103:4508-4517.

［14］CHEN B J，SIEPMANN I J.Improving the efficiency of the aggregation volume bias Monte Carlo algorithm［J］.The Journal of Physics Chemistry B，2001，105:11275-11282.

［15］CARRIER C.2009 ASHRAE handbook—fundamentals(SI)［M］.Atlanta:American Society of Heating，Refrigerating，and Air Conditioning Engineers Inc，2009:30，70.

［16］戴永庆.溴化锂吸收式制冷空调技术实用手册［M］.北京:机械工业出版社，1999:25-37.

DAI Yongqing.Lithium bromide absorption refrigeration and air conditioner practical manual［M］.Beijing:China Machine Press，1999:25-37.

［17］ZHOU J，LU X H，WANG Y R.Molecular dynamics study on ionic hydration［J］.Fluid Phase Equilibria，2002，194-197:2573-270.

［18］LAGACHE M，UNGERER P，BOUTIN A，et al.Prediction of thermodynamic derivative properties of fluids by Monte Carlo simulation［J］.Physical Chemistry Chemical Physics，2001，3:4333-4339.

［19］CADENA C，ZHAO Q，SNURR R Q，et al.Molecular modeling and experimental studies of the thermo-dynamic and transport properties of pyridinium-based ionic liquids［J］.The Journal of Physics Chemistry B，2006，110:2821-2832.

［20］KELKAR M S，SHI W，MAGINN E J.Determining the accuracy classical force fields for ionic liquids:atomistic simulation of the thermodynamic and transport properties of 1-ethyl-3-methylimidazolium ethylsulfate (［emim］［EtSO4］)and its mixtures with water［J］.Ind Eng Res，2008，47:9115-9126.

［21］FRISCH，TRUCKS，SCHLEGEL.Gaussian 03［Z］.Pittsburgh，PA:Gaussian Inc，2003.

［22］National Institute for Standards and Technology(NIST).Computational chemistry comparison and Benchmark database［EB/OL］.［2011-03-01］.http://srdata.nist.gov/cccbdb/.