模擬盒參數對汽液體系模擬過程的影響

何 川，李少華，時 雯

（重慶大學動力工程學院教育部低品位能源利用技術及系統重點實驗室，重慶400044）

1 引 言

關于汽液界面的研究一直是人們關注的熱門課題，不僅在實際工程領域有著廣泛應用，而且也是基礎理論研究的重要方面.Adamason［1］指出，界面層區域是一種介于氣相區與液相區的過渡態區域.界面層很薄，通常只有幾到幾十個分子的厚度，隨時間和空間變化而不斷漲落起伏，現有的實驗技術在汽液界面的研究上還存在一些困難.

分子動力學模擬是近年來發展起來的數值計算方法，一方面可以從微觀角度模擬分子行為，另一方面也可以通過大量統計計算獲得宏觀熱力學參數，因此變成了汽液界面研究中的有效方法.Chapela［2］等人首次提出了將汽液平衡體系分為氣相區、液相區和氣-液界面區三部分的模擬體系.基于這種方法，國內外學者做了很多相關研究：Holcomb［3］等人的研究表明，汽液界面的密度、張量等熱力學參數是連續變化的，可依據密度變化定義體系中的汽相區、過渡區和液相區；Mecke［4］等人研究了截斷半徑對界面特性的影響，提出LRC方法；劉朝、曾丹苓［5］進一步研究了截斷半徑和分子數對界面張力的影響情況.

通過分析文獻［3］～［5］以及筆者的模擬結果，發現這些傳統的模擬方法與實驗測量相比存在較為明顯的偏差，過去對于分子模擬方法的改進研究大多針對勢函數截斷半徑及模擬分子數，往往忽略了模擬盒參數的影響.通過增加分子數量，延長截斷半徑的確可以提高計算結果的準確性，但是隨之而來的則是計算時間呈幾倍乃至十幾倍的增加，計算效率大大降低，優化模擬的代價很高昂.為了提高模擬結果的準確性，同時也確保不增加模擬時間及計算機壓力，本文從模擬盒參數的角度入手，對飽和汽液共存體系的分子行為進行模擬，研究系統切片數、汽相空間尺寸、液膜厚度以及時間步長等模擬盒因素對結果的影響，為進一步改進和優化分子動力學模擬方法提供依據.

2 模擬方法

模擬對象為常用的Lennard-Jones流體氬，其勢能函數形式為：

其中，σ=3.405×10-10m，ε=1.67×10-21J.

由于受到計算能力限制，通常選取較少的分子數目進行模擬，又為避免因分子太少引起尺度效應，參考文獻［5］，選取1000作為合適的模擬分子數.模擬系綜為NVT 正則系綜，即系統的分子數n、容積V 和溫度T 保持不變.采用直角坐標系，系統在x、y、z方向的尺寸分別為Lx、Ly、Lz.

初始時刻的分子為飽和液相分子，速率呈Maxwell分布，以正方體FCC 結構均勻密集分布在系統中央，形成一塊扁平的薄液膜層，液膜表面垂直于z方向，如圖1所示.液膜質心在原點（0，0，0）處，為防止系統出現宏觀漂移，模擬中需要不斷調整液膜質心為系統的中心位置.勢能截斷半徑為3.5σ，為提高計算效率，采用Verlet-List搜尋法，搜尋半徑為4σ.模擬體系在xyz 方向均采用周期性邊界條件，不考慮外勢場（如重力場）的作用.為進一步提高效率和節約內存，本文采用動態存儲方法，合理的分配和釋放包含分子信息的內存空間.

采用leapfrog差分算法［6］求解粒子的運動方程，在計算過程中，需要不斷對粒子速度進行修正，以確保系統溫度維持在設定值.

圖1 模擬系統示意圖Fig.1 Simulation system diagram

由于從純飽和液態變化到汽液共存態需要較長時間，為達到汽液平衡狀態，每一次模擬長度規定為1×105分子動力學步.當系統處于平衡時，不僅系統的能量和氣相分子數保持不變，而且密度和壓力張量等分布隨時間也不產生變化.平衡后的統計長度為1×105分子動力學步，每隔100步統計一次.將系統沿z方向平行于x-y 平面平均分為NS個切片進行統計，溫度、壓力張量以及表面張力的統計公式如下［7，8，9］：

其中，ρ（k）為第k 切片的分子數密度，T（k）為第k切片的溫度，pN（k）、pT（k）為第k切片的法向應力和切向應力，nk為第k 切片內的分子數，vi為分子i的速度，γ（k）為第k 切片的局部張力，γ 為液膜的表面張力，xij、yij、zij、rij分別為分子i 和j 之間在x、y、z方向及空間r 的距離，Vsl=LxLyLz／NS為一個切片的體積，U＇為勢函數對r 的導數，＜＞為正則系綜統計平均值.統計時，如果i和j 兩個分子都在第k 切片內，則將其值全部計入，如只有其中一個分子位于第k切片內，則計入值的一半.

在模擬過程中，為確保體系各物理量之間的數量級較為接近，減小模擬過程中出現的截斷誤差，將所用到的各物理量進行無量綱化.無量綱化處理后的無因次參數如下：

長度L＊=L／σ；時間t＊=t／（σ（m／ε）1／2）；溫度T＊=KBT／ε；表面張力γ＊=γσ2／ε；密度ρ＊=ρσ3／m；壓力張量P＊=Pσ3／ε；作用勢U＊=U／ε

3 結果與分析

2.1 傳統模擬結果

針對溫度為84K～140K、具有相同分子數的汽液共存系統進行了傳統方法的分子動力學模擬.其他參數設置如下：液膜厚度為16分子數，汽相空間為50σ，系統切片數為100，時間步長為5×10-15s，具體參見表1.

表1 分子動力學模擬的參數Table 1 Parameters of molecular dynamics simulation

圖2給出了不同溫度下的系統的密度分布情況，圖中無量綱尺寸z＊=z／σ.可以看出，隨著溫度的升高，汽液過渡區厚度增加，汽相密度相應增大，液相密度隨之降低.當溫度達到140K 時，汽相與液相密度竟然已經大體相等，很難嚴格區分氣態分子與液態分子，這是由于溫度已經接近氬的臨界溫度（Tc=150K），處于該溫度下的飽和液和飽和氣具有相同的臨界狀態［10］.

圖3給出了不同溫度下的系統的表面張力變化情況，隨著溫度的提高，表面張力逐漸降低.通過與文獻［11］實驗測量值的比較，存在較為明顯的偏差，模擬結果不太理想.因此有必要對其進行優化改進，既要提高結果的準確性和可靠性，也要確保模擬的高效與時間的節約.

圖2 不同溫度下系統的密度分布，（z＊=z／σ）Fig.2 Density distribution of the system with temperature T，（z＊=z／σ）

圖3 不同溫度下系統的表面張力變化Fig.3 Surface tension changes with temperature T

3.2 模擬盒參數影響情況

在傳統分子動力學模擬的基礎上，筆者從模擬盒參數的角度入手，對飽和汽液體系進行模擬，研究當系統溫度為100K 時，系統切片數、汽相空間尺寸、液膜厚度以及時間步長等模擬盒因素對模擬計算結果的影響.

圖4給出的是不同切片數下表面張力的變化情況.可以看出，切片數NS從50增加到200，系統的表面張力基本不變化，這是由于式（2）～（7）均為全遍歷性的統計方式，其計算結果與切片數無關.切片數量僅僅會影響局部參數（如汽液界面密度、局部張力等）的精確性.

圖4 表面張力γ＊隨切片數Ns的變化情況Fig.4 Surface tensionγ＊ changes with the number of slices Ns

圖5、圖6給出的是汽液界面密度以及系統表面張力隨體系z軸尺寸的變化情況，圖中的無量綱尺寸Lz＊=Lz／σ.Lz 較小時，模擬體系的汽相空間較小，無法體現周期性無限空間的特點，分子運動將受到空間大小的限制；Lz 過大則會導致汽相密度增加，用于統計界面張量的液相分子數相應減少，造成統計誤差.圖6的結果表明了，在本文模擬的工況中，當Lz＊為60時，表面張力的計算結果為0.647，與文獻［11］中提供的實驗值0.654十分接近，誤差僅為1.1%.

圖5 不同Lz＊時界面密度的變化情況，（Lz＊=Lz／σ）Fig.5 Interface density changes with Lz＊，（Lz＊=Lz／σ）

圖6 表面張力γ＊隨Lz＊的變化情況，（Lz＊=Lz／σ）Fig.6 Surface tensionγ＊changes with Lz＊，（Lz＊=Lz／σ）

圖7和圖8給出的是相似分子數下，初始時刻不同的液膜厚度對平衡后表面張力及液相密度的影響情況.實際情況中的液膜比較厚，兩側界面的分子不會互相干擾，而模擬計算中，為節約計算時間，液膜通常設置的比較薄，只有幾到十幾個分子厚度，這就會導致液膜兩側界面分子相互作用，產生較大的誤差.從圖中可以看出，表面張力會隨著液膜厚度的增加而接近實驗值，規定無量綱厚度d＊=d／σ，當d＊達到13以后，液相密度趨于穩定，兩側界面分子將不會相互影響，與實際情況相符，滿足模擬要求.

圖9、圖10給出了當系統處于平衡狀態后改變分子動力學步長Δt會對液相密度及張力造成的影響.從圖中可以看出，改變步長不會影響到平衡后的體系密度，但會對表面張力的計算影響較明顯.當步長選擇為11fs（10-15s）時，結果較為接近實驗值.

圖7 表面張力γ＊隨d＊的變化情況，（d＊=d／σ）Fig.7 Surface tensionγ＊changes with d＊，（d＊=d／σ）

圖8 液相密度隨d＊的變化情況，（d＊=d／σ）Fig.8 Liquid density changes with d＊，（d＊=d／σ）

圖9 液相密度隨步長的變化情況Fig.9 Liquid density changes with time step

圖10 表面張力γ＊隨步長的變化情況Fig.10 Surface tensionγ＊changes with time step

4 結 論

采用分子動力學方法模擬飽和汽液共存體系中汽液界面現象，得到了溫度為100K 時，系統切片數、汽相空間尺寸、液膜厚度以及時間步長等模擬盒因素對結果的影響，結果表明：

（1）傳統的分子動力學模擬方法忽略了模擬盒參數的設置，計算結果與實驗測量值相比存在較為明顯的差距.

（2）系統切片數不會影響表面張力以及兩相密度的模擬結果，僅僅會對局部參數的精確性產生一定程度的影響.

（3）汽相空間尺寸會對模擬過程產生較大的影響.無量綱尺寸Lz＊為60時，表面張力的模擬結果與實驗值吻合較好，誤差僅為1.1%.

（4）液膜厚度關系到兩側界面分子是否會相互干擾.當厚度達到13，液相密度趨于穩定，兩側界面分子將不會相互影響.

（5）分子動力學步長不會影響平衡后的體系密度，但對表面張力的計算影響較明顯，當步長為11 fs，計算結果理想.

.

［1］ Adamson A W.Physical chemistry of surface［M］.New York：John Wiley ＆Sons，1982.

［2］ Chapela G A，Saville G，Thompson S M，et al.Computer simulation of a gas-liquid surface ［J］.Journal of the Chemical Society，Faraday Transactions，1977，73（2）：1133.

［3］ Holcomb C D，Clancy P，Zollweg J A.A critical study of the simulation of the liquid-vapor interface of a Lennard-Jones fluid［J］.Mole.Phys.，1993，78（2）：437.

［4］ Mecke M，Winkelmann J，Fischer J.Molecular dynamics simulation of the liquid-vapor interface：The Lennard-Jones fluid［J］.J.Chemi.Phys.，1997，107（21）：9264.

［5］ Liu C，Zeng D L.The effect factor of surface tension of liquid／vapor in molecular dynamics simulation［J］.Journal of Chongqing University：Nat.Sci.ED.，2002，25（4）：76（in Chinese）［劉朝，曾丹苓.汽液界面表面張力模擬中的影響因素［J］.重慶大學學報：自然科學版，2002，25（4）：76.

［6］ Rapaport D C.The art of molecular dynamics simulation［M］.Cambridge：Cambridge University Press，1995.

［7］ Nijmeijer M J P，Bakker A F，Bruin C，et al.A molecular dynamics simulation of the Lennard-Jones liquid–vapor interface［J］.J.Chem.Phys.，1988，89（6）：3789.

［8］ Enders S，Kahl H，Mecke M，et al.Molecular dynamics simulation of the liquid–vapor interface：I.The orientational profile of 2-center Lennard–Jones and of Stockmayer fluid molecules［J］.Journal of Molecular Liquids，2004，115：29.

［9］ Zhdanov E R，Fakhretdinov I A.Molecular dynamics simulation of the liquid–vapor interface of binary mixtures［J］.Journal of Molecular Liquids，2005，120：51.

［10］ Zeng D L，Ao Y，Zhang X M，et al.Engineering thermodynamics［M］.Beijing：Higher Education Press，2002（in Chinese）［曾丹苓，敖越，張新銘，等.工程熱力學［M］.北京：高等教育出版社，2002］

［11］ Ma Q F，Fang R S.Practical thermal physical properties manual［M］.Beijing：China Agricultural Mechanical Press，1986（in Chinese）［馬慶芳，方榮生.實用熱物理性質手冊［M］.北京：中國農業機械出版社，1986］