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

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

2)(中國(guó)科學(xué)院力學(xué)研究所,非線性力學(xué)國(guó)家重點(diǎn)實(shí)驗(yàn)室,北京 100190)

3)(中國(guó)科學(xué)院力學(xué)研究所,微重力國(guó)家實(shí)驗(yàn)室,北京 100190)

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

1 引 言

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

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

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

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

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

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

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

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

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

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

2 非均勻系統(tǒng)局部平均壓力張量?jī)煞N形式的簡(jiǎn)潔推導(dǎo)

2.1 第一種平均形式的推導(dǎo)

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

在(5)式中

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

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

2.2 第二種平均形式的推導(dǎo)

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

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

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

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

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

3.1 理論分析

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

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

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

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

在長(zhǎng)方體邊長(zhǎng)足夠大,遠(yuǎn)超過分子間有效作用程長(zhǎng)時(shí),M1,M2,M3分別與長(zhǎng)方體的體積、表面積及菱長(zhǎng)近似成正比.為方便計(jì)算,無(wú)論長(zhǎng)方體邊長(zhǎng)取多少,這些粒子對(duì)對(duì)位形壓力張量的貢獻(xiàn)都分別簡(jiǎn)稱為體貢獻(xiàn)、面貢獻(xiàn)和線貢獻(xiàn)

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

3.2 均勻流體系統(tǒng)位形壓力三項(xiàng)貢獻(xiàn)的分子動(dòng)力學(xué)模擬

3.2.1 分子動(dòng)力學(xué)模擬的方案

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

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

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

3.2.2 模擬結(jié)果

在180 K下的氬系統(tǒng)弛豫平衡后的100萬(wàn)時(shí)步內(nèi),在體系中取邊長(zhǎng)L?不同的立方體,對(duì)相關(guān)物理量分別進(jìn)行統(tǒng)計(jì)平均,得到該溫度下各個(gè)尺寸的列于表1中.L?3與立方體含有的粒子平均數(shù)成正比.的關(guān)系如圖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 ? 的關(guān)系Fig.4.Relation between ,,, and L ?.

3.3 均勻流體位形壓力三部分貢獻(xiàn)的尺度分析

3.3.1 大L?分析

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

面貢獻(xiàn)為

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

線貢獻(xiàn)為

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

根據(jù)力學(xué)平衡條件,總平均壓力與所取的L?無(wú)關(guān),由(11)式可知,平均動(dòng)壓力也與L?無(wú)關(guān),所以總平均位形壓力與所取的L?無(wú)關(guān).平均位形壓力的體貢獻(xiàn)可以表示為

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

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

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

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

3.3.2 較小尺度分析

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

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

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

3.4 溫度對(duì)位形壓力的影響

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

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

4 結(jié) 論

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