崔樹穩(wěn) 劉偉偉 朱如曾 錢萍

1)(滄州師范學院物理與信息工程學院,滄州 061001)

2)(中國科學院力學研究所,非線性力學國家重點實驗室,北京 100190)

3)(中國科學院力學研究所,微重力國家實驗室,北京 100190)

4)(北京科技大學數(shù)理學院,北京 100083)

1 引 言

對于平衡均勻無外場系統(tǒng),Clausius[1]和Maxwell[2]基于維里定理,給出了壓力公式:對于體積V中包含N個粒子的宏觀體系,平衡后的壓力為

其中i,j表示粒子;mi是粒子i的質(zhì)量;vi是粒子i的速度;rij是i,j之間距離;fij表示i,j之間的作用力,吸引力為正;〈 〉表示系綜平均,也即時間平均.(1)式包含兩部分,第一項是原子運動所貢獻的壓力,稱為動壓力,第二項源于原子之間的相互作用力,稱為位形壓力.

對于非均勻系統(tǒng),例如氣-液共存時的表面過渡層,壓力應是張量,且與位置有關.應力張量分布函數(shù)σα,β(r)需要定義,使其滿足連續(xù)介質(zhì)的動量方程

其中J是動量密度,n是粒子數(shù)密度,?ext是外場勢.Kirkwood和Buff[3]最早就對勢情況提出空間點R處應力張量的構(gòu)成方法.Irving和Kirkwood[4]給予更簡潔的表述:作用在點R處面積元 dA上的力等于所有連線通過該面積元的分子對之間作用力之和,加上由于分子通過面積元 dA的運動在單位時間內(nèi)交換的動量.Harasima[5]最早認識到滿足連續(xù)介質(zhì)動量方程(2)式的應力張量不唯一,并提出了另一種針對液體表面層構(gòu)造應力張量的方法.Schofield和 Henderson[6]證明將 Irving和 Kirkwood定義中的分子對之間的連線改為曲線,其他不變,也是一種應力張量的正確構(gòu)成方法,他們還給出了適合于多體勢情況的應力張量的構(gòu)成方法,使Irving和Kirkwood(IK)方法和Harasima方法成為他們的特例.顯然,對于對勢情況而言,IK 方法是最自然和最簡單的,因而也是后來對勢情況下,計算局部平均應力張量許多方案的基礎[7,8].

對于由N個粒子組成的非均勻系統(tǒng),從Irving和Kirkwood的定義很容易得到作為空間點(R)函數(shù)的壓力張量的公式,

其中i,j表示粒子;α,β表示方向;mi是粒子i的質(zhì)量;分別表示i粒子的動量在α,β方向的分量;表示矢量rj?ri在β方向上的分量;?(rij)表示粒子i與粒子j之間的相互作用勢.

(3)式的理論意義是顯然的.但是由于其中含有d函數(shù),不能直接用于實際測量和分子動力學模擬,所以被稱為微觀壓力張量[9].Cormier等[9]采用傅里葉變換方法得到了(3)式的局部平均形式,即對于由N個粒子組成的非均勻流體中體積為V的局部小區(qū)域,平均壓力張量為

比較(1)和(4)式可知,就壓力而言,它們表示的都是體積V中的平均值,區(qū)別在于前者的V是被剛性邊界隔離著的均勻系統(tǒng)的總體積,而后者的V只是均勻或不均勻系統(tǒng)中的一小塊的體積,所以兩者明顯不是等價的.一些研究者就非均勻系統(tǒng)中如果用前者代替后者將會引起的顯著誤差進行了討論[10?12].

(4)式適用于固體、液體、氣體等各種系統(tǒng),因此應用價值十分廣闊.例如Li等[13]采用巨正則蒙特卡羅方法模擬了超臨界Lennard-Jones流體在納米狹縫中的吸附行為,他們在(4)式中加進了壁-液勢的影響項,這也相當于對(1)式做了修正,獲得了吸附流體的壓力.對于局部壓力在原子量級的計算,一直是人們研究的熱點,例如Torres-Sánchez等[14]以及Chen[15]對此進行了深入的研究.Yu和Jin[16]用硬核雙Yukawa流體混合物模型研究膠體的熱力學和結(jié)構(gòu)特性時,在(4)式的基礎上計算了膠體體系的主體相壓力.

本文將用更為簡潔的方法推導出適用于均勻、非均勻系統(tǒng)的局部平均壓力張量普遍公式(4).鑒于在應用分子動力學計算(4)式時需要考慮計算耗時的最優(yōu)化問題,但尚未見有關報道.作為一個最簡例子,本文將給出均勻流體系統(tǒng)在以原子直徑為長度單位的局部平均尺寸L?較大條件下,局部平均位形壓力中的三部分貢獻項(體貢獻項、面貢獻項和線貢獻項)與平均尺寸的關系式(含有待定參數(shù)),這是計算最優(yōu)化方案的依據(jù);以氬原子氣體系統(tǒng)為例,用分子動力學模擬方法確定待定參數(shù),并給出位形壓力三部分貢獻項在大尺寸和小尺寸L?下各項行為的特點及溫度的影響.這將為壓力張量的分子動力學模擬計算時選項的最優(yōu)化方法提供一個范例.

2 非均勻系統(tǒng)局部平均壓力張量兩種形式的簡潔推導

2.1 第一種平均形式的推導

第一種平均形式就是(4)式.推導如下:將(3)式對局部小體積V求平均

在(5)式中

將(6)式代入(5)式即可以得到(4)式.(4)式包含兩部分:第一部分是粒子的運動對壓力的貢獻,稱為動壓力,用表示;第二部分來源于粒子之間的相互作用,稱為位形壓力,用表示.

在各向同性的平衡條件下,(7)—(9)式可以簡化為

2.2 第二種平均形式的推導

在由N個分子組成的流體系統(tǒng)中,取體積為V的局部小長方體,如圖1所示.設流體系統(tǒng)中分子對的聯(lián)線與長方體V的交集,所形成的非零長度的線段的總數(shù)為所有這些非零長度線段所構(gòu)成的集合記為于是(9)式可以改寫為

將(13)式代入(7)式,可以得到平均壓力張量,

(14)式就是平均形式(4)的第二種表示形式.容易證明,(4)和(14)式中的體積V可以取任何形狀,而不只限于長方體.只是要注意,此時一個分子對的聯(lián)線與區(qū)域V所交的線段可以超過一個,求和時,應遍及全部線段.

圖1 體積為 V 的長方體系統(tǒng)示意圖Fig.1.Schematic figure of a rectangle with volume V.

3 均勻和非均勻系統(tǒng)位形壓力張量局部平均形式的分析

3.1 理論分析

在由N個分子組成的系統(tǒng)(包括均勻系統(tǒng)和非均勻系統(tǒng)情況)中,取體積為V的小長方體,如圖1所示.在(13)式中,對長方體V內(nèi)的平均壓力張量位形部分有貢獻的粒子對中的兩個粒子之間的距離必須小于分子間有效作用程長,其相對位置有三種主要情況:第一種是兩個粒子都在長方體內(nèi),lij(V)=1,如圖1所示;第二種是一個粒子在長方體內(nèi),一個在長方體外,lij(V)<1,如圖2所示;第三種情況是兩個粒子都不在長方體內(nèi),但交長方體的側(cè)面于兩點,lij(V)<1,如圖3所示.三種情況的數(shù)量分別記為M1,M2,M3.

圖2 一個粒子在 V 內(nèi),一個粒子在 V 外,只有一個交點示意圖Fig.2.Geometry for calculating the contribution to the pressure from a pair of molecules i and j with only one intersection.

圖3 兩個粒子都在V外有兩個交點示意圖Fig.3.Geometry for calculating the contribution to the pressure from a pair of molecules i and j with two intersections.

粒子對中有至少一個粒子位于長方體V的邊界面上或棱上的情況,雖然不是完全不可能,但是概率幾乎為零,對計算平均位形壓力的貢獻可以忽略,故不計入.于是(13)式的求和上限M滿足關系式

在長方體邊長足夠大,遠超過分子間有效作用程長時,M1,M2,M3分別與長方體的體積、表面積及菱長近似成正比.為方便計算,無論長方體邊長取多少,這些粒子對對位形壓力張量的貢獻都分別簡稱為體貢獻、面貢獻和線貢獻

于是V中的平均位形壓力張量總計為

3.2 均勻流體系統(tǒng)位形壓力三項貢獻的分子動力學模擬

3.2.1 分子動力學模擬的方案

采用氬原子的平衡系統(tǒng)作為研究對象進行模擬和分析.氬的初始位型采用簡立方的點陣結(jié)構(gòu),原子間距1.2σ.模擬體系尺寸為:x×y×z=18σ×18σ×18σ,x,y,z方向上均采用鏡像邊界條件.氬原子之間采用截斷距離為 3σ的Lennard-Jones勢能函數(shù)來描述

其中ε為勢能參數(shù),σ為原子直徑,r為原子中心之間的距離.對于氬原子σ=0.3405 nm,ε=kB×120 K,其中kB=1.38×10?23J/K.

模擬中采用無量綱化量,分別以σ,ε和氬原子的質(zhì)量m=6.63382×10–26kg 作為長度、能量和質(zhì)量單位.經(jīng)過無量綱化之后,給出了其他物理量的標度:比如溫度的無量綱量T?=kT/ε,180 K相當于 1.5;時間的無量綱體系的演化采用 Velocity-Verlet方法,截斷長度rcutoff=3.0σ,弛豫過程采用溫度為180 K的NVT系綜,時間步長取 δt=5 fs,弛豫 50 萬個時間步.在達到遲豫平衡之后,采用累積平均方法計算物理量的時間平均值

3.2.2 模擬結(jié)果

在180 K下的氬系統(tǒng)弛豫平衡后的100萬時步內(nèi),在體系中取邊長L?不同的立方體,對相關物理量分別進行統(tǒng)計平均,得到該溫度下各個尺寸的列于表1中.L?3與立方體含有的粒子平均數(shù)成正比.的關系如圖4所示.

表1 ,,和的模擬值Table 1. Values of,, and given by simulation.

表1 ,,和的模擬值Table 1. Values of,, and given by simulation.

L?p?1 p?2 p?3 p?c 0.4 0 0.014–0.082–0.069 0.8 0 0.00814–0.0773–0.069 1.2 0.035–0.04–0.067–0.0705 1.6 0.039–0.071–0.038–0.07 2 0.034–0.072–0.031–0.069 2.4 0.024–0.072–0.021–0.069 2.8 0.018–0.07–0.015–0.069 3.2 0.011–0.068–0.012–0.069 3.6 0.0053–0.067–0.0089–0.069 4–3×10–4–0.064–0.007–0.07 5–0.01–0.06–0.0054–0.07 6–0.02–0.05–0.0038–0.073 8–0.035–0.034–0.0021–0.071 10–0.045–0.024–0.0014–0.07 12–0.051–0.02–9×10–4–0.071 14–0.055–0.015–5×10–4–0.07 16–0.06–0.012–2.6×10–4–0.072 17–0.06–0.01–1×10–4–0.07

圖4 ,,, 與 L ? 的關系Fig.4.Relation between ,,, and L ?.

3.3 均勻流體位形壓力三部分貢獻的尺度分析

3.3.1 大L?分析

隨著所取平均體積的大小不同,(14)式中的平均動壓力保持不變,平均位形壓力包含的三項相對大小會有很大不同.弄清楚它們與計算尺度L?的關系,對于我們在實際應用中節(jié)省計算機時非常重要.在L?較大,即L?>8 的情況下,首先分析位形壓力中與L?之間的依賴關系,再分析與L?之間的依賴關系.

面貢獻為

此式右邊括號中的第一項的來源是,在大L?條件下,作為近似,可以先不考慮棱對面貢獻的影響,根據(jù)(17)式,面貢獻與體積的乘積應與第二種粒子對數(shù)M2成正比,故與V的總面積,即與L?的平方成正比.由于此情況下計算的面貢獻沒有考慮棱對面貢獻的影響,計算的面貢獻與實際的面貢獻有差別,因此應該校正這樣計算的面貢獻:先不計棱兩端的端點效應,棱效應與體積的乘積應與棱長,即邊長L?成正比,這就是(23)式中的第二項.依此類推,還應加上頂角的影響,即第三項它與L?無關,是個待定常數(shù).k1,k2,k3是與具體的系統(tǒng)和條件有關的常數(shù).

線貢獻為

此式右邊括號中的第一項的來源是,在大L?條件下,作為近似,可以先不考慮頂角對線貢獻的影響,依據(jù)(18)式,線貢獻與體積的乘積,與M3成正比,故與V的邊長,即與L?成正比.這樣計算的線貢獻與實際的線貢獻有差別,因此要加上頂角端點的影響,即第二項k5.k4,k5是常數(shù),由具體的系統(tǒng)和條件決定.

根據(jù)力學平衡條件,總平均壓力與所取的L?無關,由(11)式可知,平均動壓力也與L?無關,所以總平均位形壓力與所取的L?無關.平均位形壓力的體貢獻可以表示為

圖5 的擬合曲線Fig.5.Fitting curve of.

圖6 的擬合曲線Fig.6.Fitting curve of.

對于本文模擬的平衡態(tài)氬系統(tǒng),取L?較大時對應的分子動力學模擬值,進行擬合,結(jié)果給出:k1=0.24053,k2=?0.01509,k3=0.0345,k4=0.10328,k5=?0.02699,如圖5和圖6所示.

上述三個方程(23),(24),(25)與分子動力學結(jié)果擬合的方差都很小,這也證明了本文的尺度分析和推理是正確的.

3.3.2 較小尺度分析

從氬系統(tǒng)的模擬結(jié)果圖4可以看出,在L?≤4區(qū)域,三種貢獻的行為有些復雜,其中特征及物理根源分析如下.

雖然上述分析是對氬系統(tǒng)而言,但除去具體的數(shù)字k1,k2,k3,k4,k5及大小尺度L?分界線可能不同外,所有定性性質(zhì)對其他系統(tǒng)不會有實質(zhì)改變.

分析結(jié)果對于正確取舍三部分貢獻是重要的.一些文獻[6?8,12]直接將宏觀大體積整體平均定理(1)式用到局部小體積上是不準確的.原因是:1)遺漏了線貢獻和面貢獻;2)即使對(1)式做了修正,沒有遺漏面貢獻,對于面貢獻項,都用 1/2代替了(4)式右邊的lij(V),即用分子對距離的一半在β方向的投影代替了(17)式中的.對于第1)點,只有當L?遠遠大于分子有效作用距離時,這種遺漏才不會帶來明顯的誤差;對于第2)點,僅對分子有效作用距離范圍內(nèi)密度均勻的情況適用,此時lij(V)的統(tǒng)計平均值為1/2.

3.4 溫度對位形壓力的影響

在討論界面特性及氣液固相變過程中,需要研究界面壓力張量及表面張力隨溫度的變化情況[7,8].因此本節(jié)采用分子動力學模擬,進一步研究了溫度對氬系統(tǒng)位形壓力的影響.分子動力學模擬的細節(jié)與 3.2 節(jié)相同,取L?=6,溫度取值范圍T?=1.4—2.3,模擬結(jié)果如圖7所示.

從圖7可以看出,位形壓力隨著溫度的升高而升高.這是由于溫度升高,分子平均動能增加,使得粒子對處于高斥力區(qū)(rij<σ)的概率增加,斥力增大,從而位形壓力增高.

4 結(jié) 論

對于平衡的非均勻系統(tǒng),人們推導了局部平均壓力張量的表達式.本文用更為簡潔的方法推導了這一表達式.此表達式也適用于均勻流體系統(tǒng).本文給出在局部平均尺寸L?>8 的條件下均勻流體系統(tǒng)平均位形壓力中的三部分貢獻項(體貢獻項、面貢獻項和線貢獻項)與L?的理論關系式(含有待定參數(shù)),并以氬原子氣體系統(tǒng)為例,在溫度180 K,原子數(shù)密度 0.8σ?3下,對分子間采用林納德-瓊斯勢進行了分子動力學模擬,給出了條件下三項貢獻及總位形壓力的模擬曲線,確定了L?>8條件下理論關系式中的待定系數(shù),并做了大L?和小L?下各項行為的特點分析和溫度影響的模擬分析,得到L?足夠大,才可以忽略面貢獻項和線貢獻項,在納米尺度下,忽略面貢獻項和線貢獻項,也就是忽略邊界效應會給計算帶來明顯的誤差.這些結(jié)論對于壓力張量的分子動力學模擬計算時選項的最優(yōu)化是有意義的.