溴化鋰水溶液汽液相平衡的分子模擬

周玉紅，趙宗昌，張曉冬，唐靜，辛緒亮

(大連理工大學 化工與環(huán)境生命學部，遼寧 大連 116023)

吸收式制冷和吸收式熱泵是回收和利用工業(yè)余熱最直接和有效的技術裝備，是節(jié)能減排，實現(xiàn)低碳經(jīng)濟最有效的技術手段之一.吸收循環(huán)熱力學效率與其工作流體的熱力學性質有關.目前，工業(yè)上的吸收制冷機和熱泵機組主要采用溴化鋰水溶液為工作流體，具有廣闊的市場應用前景.溴化鋰水溶液的汽液相平衡性質以及比熱容是重要的熱力學物性參數(shù)，對于吸收循環(huán)的開發(fā)與設計具有重要意義.

計算機分子模擬通過研究體系的微觀性質而得到體系的宏觀熱力學性質，其中Gibbs系綜蒙特卡洛方法(Gibbs ensemble Monte Carlo，GEMC)在預測混合物的汽液相平衡數(shù)據(jù)方面有獨特的優(yōu)勢.目前該方法的應用主要集中在模擬有機物的一元流體或二元溶液的汽液相平衡組成［1-2］和氣體在溶液中的溶解度方面［3］.關于溴化鋰水溶液汽液相平衡性質的分子模擬尚未見報道.

本文采用 NPT-GEMC［4-5］，用 MCCCS Towhee［6］軟件包進行LiBr/H2O二元組分的汽液相平衡性質的模擬;分析了不同壓力、濃度下，液相的構型能變化規(guī)律和溶液的微觀結構的變化規(guī)律;最后模擬了溴化鋰溶液的恒壓比熱容.

1 勢能模型

在溴化鋰水溶液中存在3種分子間力:1)水分子與水分子間的作用力;2)離子與離子間的作用力;3)離子與水之間的作用力.描述粒子間相互作用力的一套參數(shù)化的經(jīng)驗勢函數(shù)，稱為勢能模型.因為溶液的熱力學性質主要由分子間作用力決定，故本文中的勢能模型函數(shù)Uext包含范德華作用項和靜電作用項，不包含分子內作用項，如下:

其中范德華作用項用Lennard-Jones勢能模型描述，靜電作用項用固定點電荷庫侖勢描述.式中:εij、σij分別為粒子i和j范德華作用的勢阱和相互作用半徑;qi和qj分別為粒子i和j的電荷數(shù);rij為粒子i和j的距離;Uext(rij)為粒子i和j的分子間勢能.

水的勢能模型中常用的是固定點電荷模型［7-8］SPC/E(simple point charge/extended)，該模型在分子模擬溴化鋰水溶液的熱力學性質中應用相對廣泛［9-10］.因此水分子間作用選用該模型描述勢能參數(shù)(見表1).溴離子和鋰離子勢能模型參數(shù)分別取自分子模擬軟件MCCCS Towhee中Lybrand Ghosh McCammon 力場［11］和 Aqvist力場［12］(見表 1).不同粒子之間的Lennard-Jones勢能參數(shù)用Lorentz-Berthelot混合規(guī)則.

表1 離子和水的勢能參數(shù)Table 1 Potential parameters of water and ions

2 模擬實驗

在Linux操作系統(tǒng)下，計算均采用MCCCS Towhee 6.2.11 軟件包［6］，運用 NPT-GEMC 系綜方法進行模擬.由于工業(yè)上吸收制冷和吸收熱泵所用LiBr水溶液的質量濃度范圍為40% ～60%，故選取常用操作溫度為80℃、100℃下，所對應的平衡壓力進行模擬.模擬采用周期性邊界條件，范德華作用采用截斷半徑加尾部矯正方法計算，靜電長程相互作用采用Ewald方法計算.為提高計算效率，采用構型偏倚蒙特卡羅方法［13］以提高粒子插入的接收概率;采用聚集體積偏倚蒙特卡羅方法［14］以提高聚集流體的構型取樣速率.

每次模擬進行30×104次循環(huán)，其中前15×104次循環(huán)對體系進行預平衡，后15×104次循環(huán)用于熱力學性質的統(tǒng)計和計算，每次統(tǒng)計體系的構型數(shù)不低于107.每次循環(huán)包括N次MC移動，其中N為模擬體系中的粒子總數(shù).每次移動會選擇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為分子模擬計算得到的液相溴化鋰質量分數(shù)，xref為 ASHRAE手冊［15］中溴化鋰相平衡P-T-x公式，迭代得到的文獻值.用Δ表示模擬結果和文獻值的相對偏差:

用ARD表示平均相對偏差:

通過表2可知，相對平均偏差分別9.67%、11.31%，表明NPT-GEMC方法可以較好地模擬溴化鋰水溶液平衡組成關系P-T-x.將表2中數(shù)據(jù)作圖比較，見圖1.比較圖1(a)、(b)可知，在溴化鋰濃度接近于55%左右時，模擬值和文獻值非常接近;當濃度在55%～70%時，模擬值會高于文獻值;當濃度在40%～55%時，模擬值低于文獻值.該結果表明所選水和離子的勢能模型還不能在所有濃度范圍內完全準確的模擬平衡時的溴化鋰溶液的相平衡關系P-T-x，但是可以應用于工業(yè)上常用的濃度范圍(40% ～60%).

圖1 不同溫度下，不同壓力下的溴化鋰水溶液平衡組成Fig.1 The equilibrium concentrations for aqueous solution of lithium bromide at different temperatures and pressures

3.2 溴化鋰水溶液平衡密度

溴化鋰水溶液的汽液相平衡時，液相密度的分子模擬結果見表3.表中ρcal為分子模擬計算得到的不同溫度壓力下的液相溴化鋰溶液的密度;根據(jù)文獻［16］提供的溴化鋰溶液密度隨溫度t和濃度x的變化曲線圖及數(shù)據(jù)，運用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中數(shù)據(jù)作圖，見圖2.比較圖2(a)、(b)可知，分子模擬值ρcal在所模擬的濃度范圍內均小于文獻值ρref.其原因是:1)3.1節(jié)中汽液平衡組成模擬的誤差，直接導致了平衡密度的不準確;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)靜電作用能均為負，表現(xiàn)為吸引力;范德華作用能均為正，表現(xiàn)為排斥力;總的分子間作用能均為負;3)同一溫度下隨著壓力的增加，溴化鋰飽和溶液濃度降低，分子間的靜電作用能增大，而范德華相互作用能降低;但分子間的總的構型能逐漸增加.

3.4 溴化鋰水溶液平衡時的微觀結構分析

已有研究［17］主要采用分子動力學方法(MD)研究幾個陰(陽)離子存在于水中的水合作用，并未涉及陰陽離子對同時在水中時的溶液結構，不能反映濃溶液的真實水合作用.本文在模擬體系達到平衡時，統(tǒng)計了溶液中不同粒子的徑向分布函數(shù):g(Li+-O)、g(Li+-H)、g(Br--O)以及 g(Br--H)，如圖3、4 所示.

從圖3、4可以看出不同濃度下各粒子間的徑向分布函數(shù)具有大致相同的趨勢.圖3中Li+-O的徑向分布函數(shù)有2個極大值，即2個明顯的配位圈;Li+-H的徑向分布函數(shù)有2個配位圈;g(Li+-O)的第1個峰位在2.033～2.095?A之間，而g(Li+-H)的第1個峰位在2.675～2.755?A之間.這表明溶液中的水分子以帶負電的氧原子靠近Li+.Li+與水具有較強的水合作用.圖4中Br--H的徑向分布函數(shù)有3個明顯的配位圈，Br--O的徑向分布函數(shù)有2個配位圈，g(Br--H)的第 1個峰位在2.305～2.425?A之間，g(Br--O)的第1個峰位在3.325～3.485?A之間.這表明溶液中的水分子以帶正電的氫原子靠近Br-.Br-與水也具有較強的水合作用.

粒子i周圍r處球體內j粒子的配位數(shù)nij的定義為

式中:ρj為j粒子的數(shù)密度.本文根據(jù)gij(r)求出第一配位圈中陰、陽離子與氧原子、氫原子的配位數(shù)nij(r)，見表5.由于離子與氫原子的徑向分布函數(shù)統(tǒng)計的只是離子周圍1號氫原子的分布，故表5中離子周圍氫原子的配位數(shù)擴大2倍.

圖3 t=100℃時不同濃度下鋰離子與水的徑向分布函數(shù)Fig.3 The radial distribution functions for Li+and water at 100℃and different concentrations

圖4 t=100℃時不同濃度下溴離子與水的徑向分布函數(shù)Fig.4 The radial distribution functions for Br-and water at 100℃and different concentrations

表5 t=100℃時不同濃度下的溴化鋰水溶液的徑向分布函數(shù)的第一峰值、位置和配位數(shù)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中可以看出隨著溶液濃度的降低，陰、陽離子與水的徑向分布函數(shù)的第1峰值均呈現(xiàn)出降低的趨勢，但是第1配位圈內的配位數(shù)均呈現(xiàn)增加的趨勢;而陰離子與陽離子的徑向分布函數(shù)的第1峰值呈現(xiàn)增加的趨勢，但第1配位圈內的配位數(shù)呈現(xiàn)降低的趨勢.這表明溴化鋰水溶液濃度降低，陰陽離子間的作用減弱，離子與水的水合作用增強.這是由于離子的水合作用受到異號離子的影響，溶液越稀影響越小，因而離子的水合作用越強，當達到無限稀釋溶液時，離子周圍的配位數(shù)最大，如Li+周圍水的配位數(shù)最大值約為4.4［17］.從表5中還可以看出溶液越稀，水分子之間的徑向分布函數(shù)的第1峰值和第1配位圈內的配位數(shù)均增加，這表明水溶液中水分子之間的作用受到離子的影響，離子的存在減弱了水分子之間的作用.

3.5 溴化鋰水溶液恒壓比熱容的模擬

式中:H為焓.上述MC方法只能計算出體系的分子間勢能，如果不考慮動能和分子內勢能，可以通過構型焓Hcon計算比熱容，構型焓的定義如下:

Kazuhiko Murakami［10］用構型焓計算了溴化鋰溶液的比熱容，其計算的大部分數(shù)據(jù)誤差在-20%～-40%左右.本文采用該方法進行了計算，計算結果CPcon及相對偏差Δcon見表6.

本文用一種更精確的方法對溴化鋰水溶液的比熱容進行了模擬.將焓值分為理想氣體焓Hid和剩余焓Hres:

式中:Uint表示分子內勢能;K表示動能.故熱容的計算公式為

理想氣體熱容可以通過實驗數(shù)據(jù)關聯(lián)、基團貢獻法［18］和量化計算方法得到［19-20］.H20 的理想氣體熱容采用量化計算的方法，用 Gaussian 03軟件包［21］，在 B3LYP/6-31+G＊＊水平下進行構型優(yōu)化和頻率分析，考慮了水分子的平動能和轉動能以及振動能，得到了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-

恒壓熱容的計算公式為單個離子，只需要考慮平動能，因而沒有必要進行量子化學的計算.溫度不影響單原子離子的理想氣體熱容，也不需要矯正.根據(jù)實驗數(shù)據(jù)［22］知2種離子的理想氣體熱容均為:20.786 J·mol-1·K-1.對于不同濃度的溴化鋰溶液可采用公式計算其理想氣體比熱容［18］:

剩余熱容可以在NPT系綜下用MC方法模擬得到:

由于溴化鋰水溶液的熱容隨溫度的變化不大，將上式簡化為

式中:T1、T2為T附近的2個溫度.取 T1=T+10，T2=T-10，并將計算出的恒壓熱容換算成為恒壓比熱容.模擬計算出的比熱容數(shù)據(jù)CP以及文獻關聯(lián)數(shù)據(jù)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，但可以應用于工業(yè)上常用的濃度范圍(40% ～60%).

2)NPT-GEMC方法可以較好的模擬溴化鋰水溶液平衡密度組成關系P-T-ρ，相對平均偏差分別8.65%、8.33%.分子模擬值 ρcal在所模擬的濃度范圍內均小于文獻值ρref.

3)靜電作用能是構型能的主要作用項，隨著溶液濃度的降低，總構型能和靜電作用均逐漸增加.

4)不同濃度下各粒子間的徑向分布函數(shù)具有大致相同的趨勢;陰陽離子與水均具有較強的水合作用;溴化鋰水溶液越稀，離子周圍第一配位圈內水的配位數(shù)增加，水合作用增強.

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］.北京:機械工業(yè)出版社，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/.