何 川，李少華，時(shí) 雯

（重慶大學(xué)動(dòng)力工程學(xué)院教育部低品位能源利用技術(shù)及系統(tǒng)重點(diǎn)實(shí)驗(yàn)室，重慶400044）

1 引 言

關(guān)于汽液界面的研究一直是人們關(guān)注的熱門(mén)課題，不僅在實(shí)際工程領(lǐng)域有著廣泛應(yīng)用，而且也是基礎(chǔ)理論研究的重要方面.Adamason［1］指出，界面層區(qū)域是一種介于氣相區(qū)與液相區(qū)的過(guò)渡態(tài)區(qū)域.界面層很薄，通常只有幾到幾十個(gè)分子的厚度，隨時(shí)間和空間變化而不斷漲落起伏，現(xiàn)有的實(shí)驗(yàn)技術(shù)在汽液界面的研究上還存在一些困難.

分子動(dòng)力學(xué)模擬是近年來(lái)發(fā)展起來(lái)的數(shù)值計(jì)算方法，一方面可以從微觀角度模擬分子行為，另一方面也可以通過(guò)大量統(tǒng)計(jì)計(jì)算獲得宏觀熱力學(xué)參數(shù)，因此變成了汽液界面研究中的有效方法.Chapela［2］等人首次提出了將汽液平衡體系分為氣相區(qū)、液相區(qū)和氣-液界面區(qū)三部分的模擬體系.基于這種方法，國(guó)內(nèi)外學(xué)者做了很多相關(guān)研究：Holcomb［3］等人的研究表明，汽液界面的密度、張量等熱力學(xué)參數(shù)是連續(xù)變化的，可依據(jù)密度變化定義體系中的汽相區(qū)、過(guò)渡區(qū)和液相區(qū)；Mecke［4］等人研究了截?cái)喟霃綄?duì)界面特性的影響，提出LRC方法；劉朝、曾丹苓［5］進(jìn)一步研究了截?cái)喟霃胶头肿訑?shù)對(duì)界面張力的影響情況.

通過(guò)分析文獻(xiàn)［3］～［5］以及筆者的模擬結(jié)果，發(fā)現(xiàn)這些傳統(tǒng)的模擬方法與實(shí)驗(yàn)測(cè)量相比存在較為明顯的偏差，過(guò)去對(duì)于分子模擬方法的改進(jìn)研究大多針對(duì)勢(shì)函數(shù)截?cái)喟霃郊澳M分子數(shù)，往往忽略了模擬盒參數(shù)的影響.通過(guò)增加分子數(shù)量，延長(zhǎng)截?cái)喟霃降拇_可以提高計(jì)算結(jié)果的準(zhǔn)確性，但是隨之而來(lái)的則是計(jì)算時(shí)間呈幾倍乃至十幾倍的增加，計(jì)算效率大大降低，優(yōu)化模擬的代價(jià)很高昂.為了提高模擬結(jié)果的準(zhǔn)確性，同時(shí)也確保不增加模擬時(shí)間及計(jì)算機(jī)壓力，本文從模擬盒參數(shù)的角度入手，對(duì)飽和汽液共存體系的分子行為進(jìn)行模擬，研究系統(tǒng)切片數(shù)、汽相空間尺寸、液膜厚度以及時(shí)間步長(zhǎng)等模擬盒因素對(duì)結(jié)果的影響，為進(jìn)一步改進(jìn)和優(yōu)化分子動(dòng)力學(xué)模擬方法提供依據(jù).

2 模擬方法

模擬對(duì)象為常用的Lennard-Jones流體氬，其勢(shì)能函數(shù)形式為：

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

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

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

采用leapfrog差分算法［6］求解粒子的運(yùn)動(dòng)方程，在計(jì)算過(guò)程中，需要不斷對(duì)粒子速度進(jìn)行修正，以確保系統(tǒng)溫度維持在設(shè)定值.

圖1 模擬系統(tǒng)示意圖Fig.1 Simulation system diagram

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

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

在模擬過(guò)程中，為確保體系各物理量之間的數(shù)量級(jí)較為接近，減小模擬過(guò)程中出現(xiàn)的截?cái)嗾`差，將所用到的各物理量進(jìn)行無(wú)量綱化.無(wú)量綱化處理后的無(wú)因次參數(shù)如下：

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

3 結(jié)果與分析

2.1 傳統(tǒng)模擬結(jié)果

針對(duì)溫度為84K～140K、具有相同分子數(shù)的汽液共存系統(tǒng)進(jìn)行了傳統(tǒng)方法的分子動(dòng)力學(xué)模擬.其他參數(shù)設(shè)置如下：液膜厚度為16分子數(shù)，汽相空間為50σ，系統(tǒng)切片數(shù)為100，時(shí)間步長(zhǎng)為5×10-15s，具體參見(jiàn)表1.

表1 分子動(dòng)力學(xué)模擬的參數(shù)Table 1 Parameters of molecular dynamics simulation

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

圖3給出了不同溫度下的系統(tǒng)的表面張力變化情況，隨著溫度的提高，表面張力逐漸降低.通過(guò)與文獻(xiàn)［11］實(shí)驗(yàn)測(cè)量值的比較，存在較為明顯的偏差，模擬結(jié)果不太理想.因此有必要對(duì)其進(jìn)行優(yōu)化改進(jìn)，既要提高結(jié)果的準(zhǔn)確性和可靠性，也要確保模擬的高效與時(shí)間的節(jié)約.

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

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

3.2 模擬盒參數(shù)影響情況

在傳統(tǒng)分子動(dòng)力學(xué)模擬的基礎(chǔ)上，筆者從模擬盒參數(shù)的角度入手，對(duì)飽和汽液體系進(jìn)行模擬，研究當(dāng)系統(tǒng)溫度為100K 時(shí)，系統(tǒng)切片數(shù)、汽相空間尺寸、液膜厚度以及時(shí)間步長(zhǎng)等模擬盒因素對(duì)模擬計(jì)算結(jié)果的影響.

圖4給出的是不同切片數(shù)下表面張力的變化情況.可以看出，切片數(shù)NS從50增加到200，系統(tǒng)的表面張力基本不變化，這是由于式（2）～（7）均為全遍歷性的統(tǒng)計(jì)方式，其計(jì)算結(jié)果與切片數(shù)無(wú)關(guān).切片數(shù)量?jī)H僅會(huì)影響局部參數(shù)（如汽液界面密度、局部張力等）的精確性.

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

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

圖5 不同Lz＊時(shí)界面密度的變化情況，（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給出的是相似分子數(shù)下，初始時(shí)刻不同的液膜厚度對(duì)平衡后表面張力及液相密度的影響情況.實(shí)際情況中的液膜比較厚，兩側(cè)界面的分子不會(huì)互相干擾，而模擬計(jì)算中，為節(jié)約計(jì)算時(shí)間，液膜通常設(shè)置的比較薄，只有幾到十幾個(gè)分子厚度，這就會(huì)導(dǎo)致液膜兩側(cè)界面分子相互作用，產(chǎn)生較大的誤差.從圖中可以看出，表面張力會(huì)隨著液膜厚度的增加而接近實(shí)驗(yàn)值，規(guī)定無(wú)量綱厚度d＊=d／σ，當(dāng)d＊達(dá)到13以后，液相密度趨于穩(wěn)定，兩側(cè)界面分子將不會(huì)相互影響，與實(shí)際情況相符，滿足模擬要求.

圖9、圖10給出了當(dāng)系統(tǒng)處于平衡狀態(tài)后改變分子動(dòng)力學(xué)步長(zhǎng)Δt會(huì)對(duì)液相密度及張力造成的影響.從圖中可以看出，改變步長(zhǎng)不會(huì)影響到平衡后的體系密度，但會(huì)對(duì)表面張力的計(jì)算影響較明顯.當(dāng)步長(zhǎng)選擇為11fs（10-15s）時(shí)，結(jié)果較為接近實(shí)驗(yàn)值.

圖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 液相密度隨步長(zhǎng)的變化情況Fig.9 Liquid density changes with time step

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

4 結(jié) 論

采用分子動(dòng)力學(xué)方法模擬飽和汽液共存體系中汽液界面現(xiàn)象，得到了溫度為100K 時(shí)，系統(tǒng)切片數(shù)、汽相空間尺寸、液膜厚度以及時(shí)間步長(zhǎng)等模擬盒因素對(duì)結(jié)果的影響，結(jié)果表明：

（1）傳統(tǒng)的分子動(dòng)力學(xué)模擬方法忽略了模擬盒參數(shù)的設(shè)置，計(jì)算結(jié)果與實(shí)驗(yàn)測(cè)量值相比存在較為明顯的差距.

（2）系統(tǒng)切片數(shù)不會(huì)影響表面張力以及兩相密度的模擬結(jié)果，僅僅會(huì)對(duì)局部參數(shù)的精確性產(chǎn)生一定程度的影響.

（3）汽相空間尺寸會(huì)對(duì)模擬過(guò)程產(chǎn)生較大的影響.無(wú)量綱尺寸Lz＊為60時(shí)，表面張力的模擬結(jié)果與實(shí)驗(yàn)值吻合較好，誤差僅為1.1%.

（4）液膜厚度關(guān)系到兩側(cè)界面分子是否會(huì)相互干擾.當(dāng)厚度達(dá)到13，液相密度趨于穩(wěn)定，兩側(cè)界面分子將不會(huì)相互影響.

（5）分子動(dòng)力學(xué)步長(zhǎng)不會(huì)影響平衡后的體系密度，但對(duì)表面張力的計(jì)算影響較明顯，當(dāng)步長(zhǎng)為11 fs，計(jì)算結(jié)果理想.

.

［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，F(xiàn)araday 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，F(xiàn)ischer 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］.重慶大學(xué)學(xué)報(bào)：自然科學(xué)版，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，F(xiàn)akhretdinov 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）［曾丹苓，敖越，張新銘，等.工程熱力學(xué)［M］.北京：高等教育出版社，2002］

［11］ Ma Q F，F(xiàn)ang R S.Practical thermal physical properties manual［M］.Beijing：China Agricultural Mechanical Press，1986（in Chinese）［馬慶芳，方榮生.實(shí)用熱物理性質(zhì)手冊(cè)［M］.北京：中國(guó)農(nóng)業(yè)機(jī)械出版社，1986］