溴化锂水溶液气液界面微观结构的分子模拟

朱蓓蓓，高洪涛

(大连海事大学 轮机工程学院，大连 辽宁 116026)

溴化锂吸收式制冷机凭借其耗电少、噪声低、结构简单、运动部件少、安全可靠，可以利用废热、废气、废水和太阳能等低品位热能驱动，对环境和大气臭氧层无害等特性，近年来得到广泛的重视，具有广阔的市场和发展前景.溴化锂吸收式制冷机以水为制冷剂、溴化锂水溶液为吸收剂.为进一步研究溴化锂水溶液气液界面的性质，需要从微观角度进行分析和探讨.近年来，分子模拟作为揭示微米/纳米尺度现象的有效手段广泛应用于电解质水溶液气液界面的微观结构的研究［1－6］.采用分子动力学模拟软件GROMACS v4.0.2［7］分别对不同温度时，不同质量分数的溴化锂水溶液的气液界面的微观结构进行了研究.

1 分子动力学模拟方法

对水分子采用SPC/E模型，溴化锂水溶液中的Li+、Br－采用点电荷位于中心的 Lennard－Jones球.势能函数由L-J项和库仑项组成.其中L－J项由下式给出:

式中:m、n代表离子或者水分子中的氧原子，r代表m与n的距离.对于水－水、水－离子、离子－离子之间的库仑静电相互作用由下面的式子给出:

式中:i、j代表不同的水分子，k、l代表离子，r代表距离.qBr－=－|e|，qLi+=－qBr－.Li+－Li+之间的LJ参数［8］为

Br－－Br－之间的L－J参数［9］为

不同原子之间的L－J参数遵循几何平均，即σmn=

首先将模拟系统(见表1)放在各个方向定义了周期性边界条件的立方体盒子中，能量最小化后，对每个系统进行1 ns的NPT系综(压强为1 atm，温度分别为303.15，363.15，423.15，483.15 K)的模拟用来确定平衡时立方体盒子的边长.恒压方法采用Parrinello－Rahman法，其中时间常数tau_p取0.2 ps;恒温方法采用Nose－Hoover热浴，其中时间常数tau_ t取0.1 ps.不同的温度下，平衡时盒子的边长L不同.然后在z轴方向将2个与平衡时刻的盒子大小相等的真空盒子加在平衡时刻盒子的两边，形成下一步将要进行模拟的初始构型(体积为L×L×3L).使用NVT系综进行模拟.采用蛙跳法积分运动方程，时间步长为1 fs.对水分子使用SETTLE算法.采用Ewald加和方法计算库仑相互作用，库仑相互作用和范德华相互作用的截断半径都取最短盒子边长的一半.模拟用时3 ns，前1 ns用来为系统进行充分的弛豫，系统模拟1 ns达到平衡状态后，每1 000步保存位置坐标、速度坐标、力以及能量，用作后续的结果分析.

表1 系统所含离子、水分子的数目Table 1 The number of lithium ions，bromide ions and water molecules

2 分析和讨论

分别对不同质量分数的溴化锂水溶液不同温度时气液界面的微观结构、密度分布、界面厚度，溴化锂水溶液界面处和液相处离子与水分子中氢、氧原子的径向分布函数以及离子周围水分子取向角的分布函数进行了计算与分析.

2.1 密度分布

图1(a)、(b)、(c)分别表示质量分数分别为20%、40%、60%的溴化锂水溶液303.15 K时，水分子、Li+、Br－沿着界面法线方向的密度分布，其中对水分子密度的单位取 kg/m3，对 Li+、Br－单位取mol/L.发现:在液膜中心，离子的密度曲线达到一个平稳的值，特别是溴化锂水溶液质量分数较小的时候，见图1(a)和(b)，此时，Li+、Br－的密度分布曲线几乎重合，当所含离子数目较多时，见图1(c)，阴离子和阳离子的密度曲线于液相处在某个值附近呈现小幅振荡;离子倾向于分布在偏离界面处，与采用非极化力场所得到的离子在界面处负吸附的结果一致［2－3，5－6］.离子偏离界面的距离大约小于一个水分子直径(水分子直径约0.4 nm)，说明离子在近界面处发生水合作用.而且，随着溴化锂水溶液质量分数的增加，离子的密度曲线出现一个明显的峰值，即当溶液质量分数较高时(图1(b)和(c))，近界面处，离子密度出现极大值，同时，水密度曲线达到极小值，再一次说明Li+、Br－作为水合离子存在于近界面处;与采用SPC/E水模型、Lennard－Jones势函数和库仑点电荷模型计算得到的NaF、NaCl、NaBr水溶液的密度分布的情形类似［5］;与采用极化力场模拟的NaBr、NaI的密度分布不同［4］，当采用极化模型模拟时，离子的密度曲线不重合，而且Br－、I－等阴离子在界面处发生正吸附.表2说明，温度一定时，随着溴化锂水溶液质量分数的增加，体系2、3、4的液相密度逐渐增加，界面厚度逐渐减小.界面厚度的减小是由于液相中存在的离子对界面区域的水分子有吸引作用，这也是导致气液界面张力增加的原因之一.对于质量分数为60%的溴化锂水溶液，随着温度的升高，液相密度减小，气液界面厚度增加.

图1 303.15K时不同质量分数时溴化锂水溶液中的水、Li+、Br－沿着界面法线方向的密度分布Fig.1 The density of water and Li+，Br－along the normal of the interface for different concentrations at 303.15 K

表2 系统液相的密度、界面区域的中心位置和界面厚度的拟合结果Table 2 Fitted bulk liquid density，interfacial width，and the center of interfacial region from the z－dependent density profiles at different temperatures for water

2.2 离子周围水分子的结构

为研究离子周围水分子的结构以及这种局部结构是否受气液界面的出现、温度的改变以及溴化锂水溶液质量分数的影响，本节计算了不同温度时，不同质量分数的溴化锂水溶液气液界面处、液相处离子与水分子中氢、氧原子的径向分布函数以及离子周围水分子取向角的分布.

选取体系4研究303.15K时，质量分数为60%的溴化锂水溶液中，Li+、Br－周围水分子的结构以及这种局部结构是否受气液界面的影响.图2(a)、(b)分别表示位于近界面处、液相的Li+－O、Li+－H、Br－－O、Br－－H间径向分布函数.对于 Li+，Li+－O径向分布函数第一峰值较高，体现了Li+与氧间存在较强的相互作用;Li+－H径向分布函数的第1峰位比Li+－O径向分布函数的第一峰位大，说明Li+周围的水分子这样排布:氧原子朝向Li+，氢原子朝向液相.文献［1］对NaCl水溶液的结构分析也得到了相似的结果:水分子中的氧原子朝向Na+，氢原子面对液相.图2(b)表明，Br－－O、Br－－H径向分布函数第1峰值较小，体现了Br－与水分子间存在较弱的相互作用;Br－－H间径向分布函数存在第2峰，这是由于水分子中有2个氢原子;Br－－O径向分布函数的第1峰位在Br－－H径向分布函数的第1峰位与第2峰位之间，说明，Br－周围的水分子这样排布:水分子的其中1个氢原子朝向Br－.图2(a)中，近界面处与液相处的Li+－O径向分布函数的第1峰位相同，说明界面的出现并没有影响Li+周围水分子的排列;后者的第1峰值比前者略高，这是因为近界面处的水分子数目比液相处少.界面处与液相处的Li+－H、Br－－O、Br－－H径向分布函数也出现了相同的情况，再次说明，界面的出现不影响离子周围水分子的结构.

计算离子周围水分子取向角的分布函数［10］以进一步研究离子周围水分子的取向.取向角是这样定义的:从氧指向离子的向量与水分子偶极向量的夹角θ.取向角分布函数定义为取向角分布的概率.将与离子的距离小于该离子与氧原子之间径向分布函数的第1峰位的水分子取为离子周围的水分子.图3表示的是，体系4近界面处以及液相处，离子周围水分子的取向分布函数.发现:无论近界面处还是液相处的Li+周围的水分子取向分布函数在160.273°出现极大值，说明对于Li+，水分子是以氧靠近离子，氢原子的取向使得水分子的偶极方向指向O－Li+连线所成向量的反向.图3(b)表明:无论近界面处还是液相处的Br－周围的水分子的取向分布函数在49.751°(约为水分子HOH键角的一半)出现极大值，说明对于Br－，意味着水分子的某一氢原子靠近Br－，而且靠近Br－的水分子的氢氧键位于Br－的径向位置，这样的取向占有主要地位;同时，该取向分布函数在125.776°出现较小的峰值，说明还有这样的取向占次要地位:水分子的某一氢原子靠近Br－，与Br－距离较远的水分子的另一氢与氧构成的氢氧键位于Br－的径向位置.

图2 体系4分别位于近界面处(［6.651 nm，6.851 nm］)及液相处(［5.8 nm，6 nm］)的Li+－O、Li+－H、Br－－O、Br－－H的径向分布函数Fig.2 The Li+－O，Li+－H，Br－O and Br－－H radial distribution functions of system 4 near the interface(［6.651 nm，6.851 nm］)and in the bulk liquid(［5.8 nm，6 nm］)

图3 体系4近界面处(［6.651 nm，6.851 nm］)及液相处(［5.8 nm，6 nm］)的Li+、Br－周围水分子的取向分布函数Fig.3 Orientation distribution functions for Li+and Br－of system 4 near the interface or in the bulk liquid，respectively

为研究温度对离子周围水分子结构有何影响，选取体系6来与前面的计算结果进行比较.图4(a)、(b)表示的是423.15 K时，位于近界面处、液相处的Li+、Br－与水分子中氧、氢之间的径向分布函数.发现与303.15K时的结果相比，径向分布函数的强度变小，这是因为随着温度的升高，分子之间的距离会变大;近界面处与液相处的径向分布函数几乎重合，说明随着温度的升高，近界面处与液相处离子周围水分子的结构极为相似.同样考察423.15 K时，离子周围水分子的取向角分布函数.图5表示体系6离子周围水分子的取向角分布函数，发现无论近界面处还是液相处的Li+周围的水分子取向分布函数在157.723°出现极大值;无论近界面处还是液相处的Br－周围的水分子的取向分布函数在大约49.792°出现极大值，在125.16°出现较小的峰值，与303.15 K时的结果相似，只是函数曲线强度变小且更加平滑，说明随着温度的升高，离子周围水分子取向的有序性不再那么明显.

为研究溴化锂水溶液的质量分数对离子周围水分子局部结构的影响，选取体系3来与体系4来进行比较.

图4 体系6位于近界面处(［6.795 nm，6.995 nm］)及液相处(［5.8 nm，6 nm］)的Li+－O、Li+－H、Br－－O、Br－－H的径向分布函数Fig.4 The Li+－O，Li+－H，Br－O and Br－－H radial distribution functions of system 6 near the interface(［6.795 nm，6.995 nm］)and in the bulk liquid(［5.8 nm，6 nm］)

图5 体系6分别位于近界面处及液相处的Li+、Br－周围水分子的取向分布函数Fig.5 Orientation distribution functions of system 6 for Li+and Br－near the interface and in the bulk liquid

图6 体系3分别位于近界面处(［6.199 nm，6.399 nm］)及液相处(［5.8 nm，6 nm］)的Li+－O、Li+－H、Br－－O、Br－－H的径向分布函数Fig.6 The Li+－O，Li+－H，Br－O and Br－－H radial distribution functions of system 3 near the interface (［6.199 nm，6.399 nm］)and in the bulk liquid (［5.8 nm，6 nm］)

图6表明，与体系4的径向分布函数相比，强度变小;而且随着溴化锂水溶液质量分数的减小，界面处与液相处离子周围水分子的局部结构的区别逐渐变小.图7表示体系3离子周围水分子的取向角分布函数，发现无论近界面处还是液相处的Li+周围的水分子取向分布函数在159.656°出现极大值;无论近界面处还是液相处的Br－周围的水分子的取向分布函数在50.409°出现极大值，在131.532°出现较小的峰值，与303.15 K时的结果相似.说明，随着质量分数的减小，离子周围水分子的取向有序性不明显.

图7 体系 3分别位于近界面处(［6.199 nm，6.399 nm］)及液相处(［5.8 nm，6 nm］)的Li+、Br－周围水分子的取向分布函数Fig.7 Orientation distribution functions for Li+and Brnear the interface(［6.199 nm，6.399 nm］)or in the bulk liquid(［5.8 nm，6 nm］)respectively

3 结论

本文采用分子动力学的方法研究了不同温度时，不同质量分数的溴化锂水溶液气液界面的微观结构.对界面法线方向密度分布的研究结果表明，离子在近界面处发生水合作用，当溴化锂水溶液质量分数较大时(60%)，离子密度曲线出现一个明显的峰值，离子在界面处发生负吸附，这是由于本文采用非极化力场进行模拟;温度一定时，随着溴化锂水溶液质量分数的增加，液相密度逐渐增加，界面厚度逐渐减小;随着温度的升高，液相密度减小，气液界面厚度增加.为研究离子周围水分子的结构以及这种局部结构是否受气液界面的影响，分别计算了界面处、液相处离子与水分子中氢、氧的径向分布函数和离子周围水分子的取向分布函数，结果表明，界面的出现并没有影响离子周围水分子的排列:对于Li+，水分子是以氧靠近离子，氢原子的取向使得水分子的偶极方向指向O－Li+连线所成向量的反向;对于Br－，意味着水分子的某一氢原子靠近Br－，而且靠近Br－的水分子的氢氧键位于Br－的径向位置，这样的取向占有主要地位，还有这样的取向占次要地位:水分子的某一氢原子靠近Br－，与Br－距离较远的水分子的另一氢与氧构成的氢氧键位于Br－的径向位置.随着温度的升高或者溴化锂水溶液质量分数的减小，径向分布函数的强度变小;近界面处与液相处离子周围水分子的结构变得相似;离子周围水分子取向的有序性不再那么明显.

［1］JARDON－VALADEZ E，COSTAS M E.Solvation properties of a polarizable water model in a NaCl solution:Monte Carlo isothermal－isobaric ensamble simulations［J］.J Molecular Structure(Theochem)，2004(677):227－236.

［2］PAUL S，CHANDRA A.Dynamics of water molecules at liquid－vapor interfaces of aqueous ionic solutions:effects of ion concentration［J］.Chemical Physics Letters，2003，373:87－93.

［3］BHATT D，NEWMAN J，RADKE C J.Molecular dynamics simulations of surface tensions of aqueous electrolytic solutions［J］.J Phys Chem B，2004，108:9077－9084.

［4］JUNGWIRTH P，TOBIAS D J.Molecular structure of salt solutions:a new view of the interface with implications for heterogeneous atmospheric chemistry［J］.J Phys Chem B，2001，105:10468－10472.

［5］DAIGUJI H，HIHARA E.Molecular dynamics study of the liquid－vapor interface of lithium bromide aqueous solutions［J］.Heat and Mass Transfer，1998，35(3):213－219.

［6］PAUL S，CHANDRA A.Liquid－vapor interfaces of simple electrolyte solutions:molecular dynamics results for ions in stockmayer fluids［J］.J Phys Chem B，2003，107:12705－12712.

［7］HESS B，KUTZNER C，SPOEL D，et al.GROMACS 4:algorithms for highly efficient，load－balanced，and scalable molecular simulation［J］.J Chem Theory Comput，2008，4(3):435－447.

［8］AQVIST J.Ion－water interaction potentials derived from free energy perturbation simulations［J］.J Phys Chem，1990，94:8021－8024.

［9］LYBRAND T P，GHOSH I，MCCAMMON J A.Hydration of chloride and bromide anions:determination of relative free energy by computer simulation［J］.J Am Chem Soc，1985，107(25):7793－7794.

［10］LEE S H，RASAIAH J C.Molecular dynamics simulation of ionic mobility.I.alkali metal cations in water at 25℃［J］.J Chem Phys，1994，101(8):6964－6974.