東南大學能源與環境學院  嚴攀  殷勇高

       【摘  要】LiCl吸濕溶液由于其節能、可循環利用等優點在空氣調節領域得到了廣泛應用，同時電場被多次證實可以用于強化氣液傳質，但除濕及電場強化的機理尚不明確。為了探究LiCl溶液吸濕的過程及為研究靜電場強化傳質效果的機理分析打下基礎，本文采用分子動力學探究靜電場對LiCl溶液的水化結構的影響，首先選取了常用的三種的LiCl水溶液的力場，并根據模擬物性的能力選出最適合LiCl的力場即Koneshan模型。因此本文用Koneshan模型研究靜電場對溶液體系的水化結構和影響傳質的表面張力的影響，發現垂直于氣液界面的電場會使水分子傾向于定向排列，從而降低了LiCl的表面張力，但靜電場并不影響離子-水分子的水合結構。

       【關鍵詞】分子動力學模擬 溶液除濕 電場強化 水合結構

       【中圖分類號】TU831.6

       【基金項目】國家自然科學基金自助項目（52076039）、國家重點研發計劃課題資助項目（2018YFC0705306）。

Abstract: LiCl solution has been widely used in the field of air conditioning due to its energy-saving and recyclability, and electric fields have been repeatedly proved to be able to enhance gas-liquid mass transfer, but the mechanism of dehumidification and electric field enhancement is unclear. To explore the moisture absorption process of LiCl solution and lay the foundation for studying the enhancement mechanism of electrostatic field, molecular dynamics simulation is adopted to explore the effect of electrostatic field on the hydration structure of LiCl solution. Three different force fields of LiCl are explored firstly, and the most suitable one is selected according the ability of reproducing physical properties. Therefore, this paper uses model of Koneshan to study the effect of electrostatic field on the Hydration structure and surface tension of the solution system. It is found that the electrostatic field perpendicular to the gas-liquid interface tends to align the water molecules, thus reducing the surface tension of LiCl, but the electric field does not affect the hydration structure of ion-water molecules.
Keywords: molecular dynamics simulation; solution dehumidification; electric field enhancement; hydration structure

       鹵鹽溶液因其吸濕性、節能性、易實現、可循環利用的特點在空氣調節領域得到了廣泛應用^[1,2]，其中LiCl溶液由于其除濕性能優越得到了廣泛的關注。但目前的研究手段主要停留在宏觀實驗層面和經驗關聯式的擬合^[3]，缺少對分子尺度的探究。同時有研究表明，外部的物理場強可以用來強化氣液傳質，其中靜電場可以通過極化效應對物系產生影響，且易于控制和調節，已經成為近年來化工過程強化領域研究和開發的熱點^[4]，靜電場對溶液除濕/再生的影響及機理鮮有報道，了解除濕/再生過程及其強化過程的微觀機理對于拓寬溶液除濕這一技術的應用和理論發展有重要意義。

       分子動力學模擬（MD）正逐漸成為多種學科范圍內的有效研究手段，已經廣泛應用于模擬生物蛋白、蒸發、冷凝、氣液界面等現象^[5,6]。對于模擬最重要的是選取能夠描述粒子間作用力的勢能函數參數即力場，針對電解質溶液的力場主要包括非極化力場^[7-9]和極化力場^[10,11]。LiCl溶液中包括水氫原子、氧原子以及Li+和Cl-，而離子的力場模型多數是在水模型的基礎上發展而來，通過擬合水合自由能、離子水合半徑、徑向分布函數、晶格能等來獲取參數，力場參數是主要是離子的勢能阱ε、軟球直徑σ以及極化率。Jorgensen[7]等報道了金屬離子F-，Cl-，Br-，以及I-和一些堿性金屬離子的Lennard Jones參數，這些參數是結合TIP4P的水分子模型得出的，同時考慮了不同邊界條件的影響，為了擬合出正確的水化能以及第一殼層的徑向分布函數，同時復現出了和從頭算得出的結果相同的一水合物的能量而優化出的結果，這是第一個完整且自洽的用于模擬的離子參數集。Joung[8]等為了復現離子的一系列性質，包括結構、運動、溶解、結晶態、以及和其他離子相互作用方面的參數，采用集中常用的水分子模型以及Ewald加和的算法，首次在優化σ和ε時把平衡陰陽離子的晶格能LE和晶格常數LC參數納入擬合的考慮范圍內。Koneshan等^[9]使用25°C的SPC/E模型和長程作用的反應場，考慮了水化結構以及滯留時間等因素，結果顯示計算結果較好。極化模型應用較多的是基于SWM4-NDP水模型的Drude模型，其中的離子參數由Yu^[12]擬合得出，但目前只在低濃度中適用良好。綜合來看，非極化力場計算量較小，且被多次證實可以準確地預測溶液的多種熱力學性質^[13]。水合結構不僅僅是發展力場的重要考慮因素，也是水和溶液存在諸多性質差異的根本原因。在電解質溶液中，離子附近形成的電場使部分溶劑水分子在離子周圍發生定向排列，離子水合層中水的微觀結構發生改變，這就是溶液中離子對水結構的影響，而電場的加入則會影響溶液本身的內部電場從而影響水化結構及水分子偶極取向^[14]，帶來宏觀物性表面張力的改變，同時表面張力是影響傳質的重要因素，因此對電場下溶液水化結構開展分子動力學模擬有重要意義。

       從上述綜述中可以看出，針對電解質溶液的力場有多種參數，而適用于LiCl溶液的力場參數并沒有最佳選擇，這也阻礙了后續的鹽溶液除濕過程模擬以及強化機理分析的腳步。考慮到模擬過程的巨大計算量和模擬時間，本文優先選擇計算效率較高的非極化力場，通過比較物理性質的實驗值和分子動力學模擬值來選擇適合LiCl溶液的力場，并應用所選擇力場模擬LiCl溶液的水化結構和對傳質影響較大的表面張力，對比研究有無靜電場下溶液體系的變化，為從分子視角研究鹽溶液除濕過程和電場強化鹽溶液-空氣體系熱質交換的機理奠定基礎。

1 方法

       本文的模擬對象為LiCl水溶液，采用開源軟件GROMACS進行300-340K系列溫度條件下的LiCl溶液的物性和結構模擬，模擬體系如圖1所示，共由596個水分子、141個Li⁺和141個Cl^-組成，針對密度模擬、比熱容等模擬時采用模型a(2.9×2.9×2.9nm)，模擬表面張力時采用模型b(2.9×2.9×8nm)。對于分子數目的限制可能會導致小樣本系統的表面效應掩蓋其本體效應，不能完全反應真實系統的性質和行為，因此在本課題研究中，采取周期性邊界條件的方法模擬濃度為35%的LiCl溶液。在模擬過程中，只考慮分子之間非鍵作用力，包括范德華力E_vdw、靜電庫倫力E_elec（包括了氫鍵）。因此，總的勢能計算如下：

       E_total=E_vdw+E_elec            （1）

       范德華力和經典庫侖力的計算式如下：

                 （2）

                  （3）

       LJ勢中的離子的勢能阱ε以及軟球直徑σ以及電荷量參數q列于表1，分別是來自Jorgensen、Joung和Koneshan的擬合數據^[7-9]。粒子之間作用力的計算方式可以分為幾何混合法則和算術混合法則^[15]，GROMACS已實現通過Settle算法對剛性水分子進行約束，為了模擬過程的不崩潰首先需要對整個系統進行能量最小化以達到勢能最小的構型，并在NPT系綜下進行溫度和壓力的預平衡，壓力在過程中始終控制在1atm，控制溫度的算法采用Velocity-rescaling算法，控制壓力的算法采用Parrinello-Rahman，處理長程靜電相互作用力采用PME的方法，為了減少計算量，使用1nm的截斷距離，時間步長取為2fs。所有的過程采用穩態模擬，分為預平衡階段和數據采集階段。預平衡過程使用2ns來確保溫度和壓力達到穩定，并用額外的1ns每隔1ps來采集數據，計算密度、定壓比熱容及擴散系數。

       為了模擬計算表面張力，需要營造一個氣液界面，將模擬盒子的高度拉伸至8nm，模擬體系在三個方向都采用周期邊界條件。模型的高度計得較長以保證足夠的汽相空間。液膜設計得足夠厚以防止液膜兩側的分子相互作用，模擬時間和步長與前述相同。當考慮靜電場后，在GROMACS中設定不同大小的恒定電場以進行控制變量法比較，如圖1(c)所示，在體系的Z軸方向施加0-2 V/nm的電場，如圖6所示。為了防止電場的引入使模擬坍塌，需要進行時間更長的預平衡，其余計算參數不變。

表1 選取的力場參數


圖1 模型示意圖（其中白色為H，紅色為O，粉色為Li，綠色為Cl）

2 力場評價結果和討論

       35%濃度的LiCl水溶液的密度如圖2所示，溫度范圍為300K至340K，總體而言，三種模型都能夠模擬出密度隨溫度變化的趨勢。與文獻中的實驗值相比，所有模型都可以再現LiCl水溶液的密度^[16]。 Jensen等人的模型給出最佳的密度估計，平均相對誤差為4.0%，其次是Koneshan等人的模型，最大相對誤差為4.4%。圖3顯示了不同溫度的LiCl水溶液中水分子的擴散系數，在300K的溫度下，Jensen模型給出的預測值約為1.26×10^-5 cm²/s，根據文獻實驗[17]，在相同工況下水在LiCl溶液中的擴散系數約為1.0*10^-5 cm²/s，總體上Jensen模型更加接近擴散系數的真實值。根據漲落理論算出定壓比熱容Cp，如圖4所示。將模擬的定壓比熱容與文獻中的實驗值^[16]相比，三種模型的表現都有較大誤差，其中復現能力最差的是Joung模型，Jensen和Koneshan模型的準確程度相似。LiCl水溶液在300 K、320K、340 K時的表面張力模擬值如圖5所示，所有的模型都能正確地模擬出表面張力隨溫度的變化趨勢。由于Jensen對表面張力的模擬誤差過大，只在此處比較Joung和Koneshan模型的準確性，與實驗值^[18]相比， Koneshan力場的模擬結果優于Joung的結果，Koneshan力場的平均相對誤差為8%，而Joung力場的平均相對誤差為10.3%。

圖2  LiCl溶液密度模擬值                 圖3　水在溶液中的擴散系數

       根據密度，比熱容，擴散系數和表面張力的模擬結果，Jensen等人提出的力場能夠以最高的精度再現實驗密度和擴散系數，但模擬表面張力表現最差。Joung等人提出的力場在物性的模擬方面沒有復現比較好的值。與此相反的是Koneshan等人提出的力場，不僅能夠產生與實驗數據相比最佳的表面張力，還具有令人滿意的其他熱力學性質復現的能力。這表明了Koneshan等人的模型優于其他兩者。因此，本文將使用Koneshan力場來進行分子動力學模擬。 

圖4 LiCl的定壓比熱容                圖5 LiCl溶液的表面張力

3 靜電場對LiCl溶液水化結構影響 

       解質溶液中陰離子和陽離子與水分子之間存在著強烈的庫侖相互作用，形成了較大的溶劑化殼層。未加電場時的溶劑-離子結構可以從圖6的徑向分布函數中得出。徑向分布函數（radial distribution function, RDF）的物理意義可由圖6顯示，圖中，中心原子為“參考分子”，與其中心的距離由r→r+dr間的分子數目為dN。定義徑向分布函數g(r)為：

       ρg(r)4πr²=dN            (4)

       式中，ρ為系統密度。徑向分布函數（RDF）可以解釋為系統的區域密度與平均密度的比值，參考分子的附近區域密度不同于系統的平均密度，但與參考分子的距離較遠時應與平均密度相同，即當r值大時徑向分布函數函數的值應該接近1。

                          圖6 徑向分布示意                  圖7  Li⁺、Cl-相對于氧原子的徑向分布函數   

       圖7顯示了陰陽離子和氧原子之間的徑向分布函數，Li⁺與氧原子的徑向分布函數在0.196 nm 處有一個強烈的峰值；Cl^-與氧原子的徑向分布函數在0.32nm處有一個峰值，該處為第一水合層，可以看出水分子中的氧原子與Li⁺更近，其原因是氧原子呈電負性，更容易和帶正電的離子相結合，因此與氧原子與Li離子之間的相互作用力更強。從圖6還可以得出在300-340K溫度區間內，溫度對徑向分布函數的影響很小，這與純水的情況不一樣，在純水中，溫度會使O-O，O-H徑向分布函數發生顯著變化，其原因是高溫帶來的氫鍵的破壞^[19]。

圖8 離子周圍水分子取向（其中a為Cl-周圍水分子取向，b為Li+周圍水分子取向）

       圖8給出了不同大小的電場施加下，Cl^-和Li⁺周圍水分子的取向角在[0°,180°]范圍內的分布概率，橫坐標為離子到氧原子的向量與水分子偶極之間的夾角，從取向角分布可以得出更具體的水化結構信息。從圖8(a)中可以看出，在127.58° 時的概率有峰值，Cl^-周圍水分子取向角大于150° 或小于100° 的概率都很小。從圖8(b)得出Li⁺周圍水分子取向角在約8°時的概率有峰值，并且角度越大，Li⁺周圍水分子取向角在該角度分布的概率越小。由上可見其O-H共價鍵指向陰離子，而氧原子與Li⁺較近，Li⁺接近位于水分子角平分線上。結合徑向分布函數的規律可以描述出整體的水分子和離子的分布情況，如圖9所示。隨著電場強度的增加的增大，離子周圍水分子取向角概率峰值位置不變，大小有輕微的改變。這就說明盡管水分子會受到電場的影響導致總體偶極矩的方向有所改變，但是由于水合作用帶動離子隨溶劑一起轉動從而保證了相對固定的水合結構。

圖9 離子周圍水分子的排布傾向                     圖10 電場對表面張力及總偶極矩的影響

       圖10給出了不同大小的電場對表面張力的影響。能夠得出對于平行電場而言，方向垂直于氣液界面的電場使表面張力減小，這可以歸因于垂直的電場可以使水分子這種極性分子趨向于定向排列^[14]，從而減小了粒子之間的作用力，表現為表面張力減小，定向排列的現象也可由圖8中Z軸方向的總偶極矩分量Mz的變化趨勢推出，由于模擬過程中采用的剛性水分子，所以每個水分子的極矩為固定的2.35D，但增大的總的Z軸偶極矩說明了水分子在Z軸方向電場的施加下傾向于Z軸的定向排列。     

4 結論

       為了從微觀角度探究LiCl吸濕溶液及為研究靜電場強化傳質效果的機理分析打下基礎，本文采用分子動力學方法，探索了描述LiCl水溶液中粒子相互作用的合適力場。根據密度、比熱容、擴散系數、表面張力等模擬數據和實驗值的對比情況，在三個最新的非極化力場中選擇了Koneshan模型，進一步選模擬電場對溶液體系的微觀結構、表面張力的影響，并進行了初步分析，主要結論如下： 

       （1）盡管非極化力場存在一定的缺陷，但由于其計算性能高且擁有復現熱力學參數的良好能力，依然有很大的利用價值。Koneshan等人的模型給出了最佳的表面張力模擬，并顯示出適中的預測其他性能的能力，這表明Koneshan等人的模型可以更好地用于包含氣液界面的模擬。

       （2）垂直于氣液界面的電場會使LiCl的表面張力減小，且靜電場會使水分子傾向于Z軸定向排列，但并不影響離子-水分子的水合結構，Li⁺比Cl^-的水合半徑更小，O-H共價鍵指向陰離子，Li⁺在偶極反方向的水分子角平分線上。

符  號  說  明

參考文獻

       [1] Yin Y G, Zhang X S. Comparative study on internally heated and adiabatic regenerators in liquid desiccant air conditioning system[J]. Building and Environment, 2010, 45(8): 1799-1807.
       [2] Zhang F, Yin Y, Zhang X. Performance analysis of a novel liquid desiccant evaporative cooling fresh air conditioning system with solution recirculation[J]. Building & Environment, 2017, 117(MAY): 218-229.
       [3] Asfand F, Stiriba Y, Bourouis M. CFD simulation to investigate heat and mass transfer processes in a membrane-based absorber for water-LiBr absorption cooling systems[J]. Energy, 2015, 91: 517-530.
       [4] 馬空軍, 賈殿贈, 孫文磊, et al. 物理場強化化工過程的研究進展. 現代化工, 2009: 33-37,39.
       [5] Mieda S. Analysis of the Interaction between a Protein and Polymer Membranes Using Steered Molecular Dynamics Simulation to Interpret the Fouling Behavior[J]. Bulletin of the Chemical Society of Japan, 2020, 93(12): 1443-1448.
       [6] Xue L, Keblinski P, Phillpot S R, et al. Effect of liquid layering at the liquid-solid interface on thermal transport[J]. International Journal of Heat and Mass Transfer, 2004, 47(19-20): 4277-4284.
       [7] Jensen K P, Jorgensen W L. Halide, Ammonium, and Alkali Metal Ion Parameters for Modeling Aqueous Solutions[J]. J Chem Theory Comput, 2006, 2(6): 1499-509.
       [8] Joung I S, Cheatham Iii T E. Determination of alkali and halide monovalent ion parameters for use in explicitly solvated biomolecular simulations[J]. The journal of physical chemistry B, 2008, 112(30): 9020-9041.
       [9] Koneshan S, Rasaiah J, Lynden-Bell R, et al. Solvent Structure, Dynamics, and Ion Mobility in Aqueous Solutions at 25 °C[J]. Journal of Physical Chemistry B, 1998, 102: 4193-4204.
       [10] Dang L. Mechanism and Thermodynamics of Ion Selectivity in Aqueous Solutions of 18-Crown-6 Ether: A Molecular Dynamics Study[J]. Journal of The American Chemical Society - J AM CHEM SOC, 1995, 117.
       [11] Dang L. Development of Nonadditive Intermolecular Potentials Using Molecular-Dynamics - Solvation of Li+ and F- Ions in Polarizable Water[J]. The Journal of Chemical Physics, 1992, 96: 6970-6977.
       [12] Yu H, Whitfield T, Harder E, et al. Simulating Monovalent and Divalent Ions in Aqueous Solution Using a Drude Polarizable Force Field[J]. Journal of chemical theory and computation, 2010, 6: 774-786.
       [13] Vega C, Abascal J L F. Simulating water with rigid non-polarizable models: a general perspective[J]. Physical Chemistry Chemical Physics, 2011, 13(44): 19663-19688.
       [14] Wang B-B, Wang X-D, Duan Y-Y, et al. Molecular dynamics simulation on evaporation of water and aqueous droplets in the presence of electric field[J]. International Journal of Heat and Mass Transfer, 2014, 73: 533-541.
       [15] Kong C L. Combining rules for intermolecular potential parameters. II. Rules for the Lennard‐Jones (12–6) potential and the Morse potential[J]. The Journal of Chemical Physics, 1973, 59(5): 2464-2467.
       [16] Conde M R. Properties of aqueous solutions of lithium and calcium chlorides: formulations for use in air conditioning equipment design[J]. International Journal of Thermal Sciences, 2004, 43(4): 367-382.
       [17] Potnis S V, Lenz T G, Dunlop E H. Measurement of water diffusivity in aqueous lithium bromide and lithium chloride solutions[J]. Chemical Engineering Communications, 1995, 139: 41-49.
       [18] Wen Y, Bjurstrom H, Setterwall F. Surface-Tension of Lithium Bromide Solutions with Heat-Transfer Additives[J]. Journal of Chemical and Engineering Data, 1991, 36(1): 96-98.
       [19] Shvab I. Molecular dynamics simulation of water and aqueous solutions[D]. Swinburne University of Technology, 2014.

       備注：本文收錄于《建筑環境與能源》2021年4月刊 總第42期（第二十屆全國暖通空調模擬學術年會論文集）。版權歸論文作者所有，任何形式轉載請聯系作者。