李 群 ，王寶和

(大連理工大學化工學院，遼寧大連 116024)

擴散系數(shù)是表征物質(zhì)運輸過程的重要參數(shù)，但采用常規(guī)實驗手段很難準確測量得到。對于受限空間的流體，擴散系數(shù)的研究更加困難。隨著計算機和分子動力學模擬技術(shù)的發(fā)展，從分子水平研究流體的擴散規(guī)律已經(jīng)引起了國內(nèi)外許多學者的極大關(guān)注［1－3］。Meier等采用平衡分子動力學方法，模擬得到了Lennard－Jones(L－J)流體的自擴散系數(shù)和黏度［2］。葛宋和陳民采用平衡分子動力學方法，通過均方位移計算得到了超臨界條件下，L－J流體的自擴散系數(shù)隨溫度的變化規(guī)律;同時，利用Green－Kubo法計算了超臨界L－J流體混合物的擴散性質(zhì)［4］。本文擬采用L－J模型，探究非受限空間中溫度、截斷半徑、模擬粒子數(shù)和受限空間中能量系數(shù) 、狹縫寬度及溫度對氬流體自擴散系數(shù)的影響。

1 模擬方法

目前，采用分子動力學模擬技術(shù)，計算自擴散系數(shù)的方法主要有兩種，即Green－Kubo法和Einstein法，分別如式(1)和(2)所示［1］。本文采用Einstein法計算氬流體的自擴散系數(shù)。

式中，Dself為自擴散系數(shù)，m2/s;N為粒子總數(shù);ui(t)和ui(0)為第i個粒子在t和0時刻的速度，m/s;ri(t)和ri(0)為粒子i在t和0時刻的位置，m;t為時間，s;〈*〉表示系綜平均。

1.1 非受限空間模擬體系的建立

非受限空間的模擬盒子如圖1所示。在模擬的初始時刻，氬流體粒子以面心立方晶格(FCC)均勻分布在邊長為L的正方體模擬盒中。模擬過程中，系統(tǒng)采用正則系綜(NVT)，并不斷調(diào)整體系質(zhì)心使之處于坐標原點;控溫方式采用Nose－Hoover控溫法;在x、y、z方向上均采用周期性邊界條件。

1.2 受限空間模擬體系建立

對于受限空間的分子動力學模擬(圖1a)，壁面有著極其重要的作用。壁面大致分為虛擬熱壁面和結(jié)構(gòu)壁面。在虛擬熱壁面體系中，不考慮壁面的微觀粒子結(jié)構(gòu)和壁面與流體之間的相互作用，入射到壁面的流體粒子被彈回到流體系統(tǒng)，流體粒子彈回速度的大小和方向由壁面的溫度和粗糙度決定。結(jié)構(gòu)壁面模擬真實壁面的原子晶格結(jié)構(gòu)，除了要考慮流體粒子自身兩兩相互作用外，還要考慮流體粒子與壁面粒子的相互作用。結(jié)構(gòu)壁面又分為剛性壁面，Phantom壁面和Einstein固體壁面［5］。

本文受限空中的壁面是由固定不動的虛擬LJ原子組成的剛性壁面，模擬盒子如圖1b所示(上面和下面的三層原子為壁面原子，中間夾層的六層粒子為氬原子)。與非受限空間一樣，在受限空間中也采用正則系綜(NVT)。x和z方向為周期性邊界條件，y方向為固體壁面邊界條件。

圖1 空間模擬體系的示意圖

1.3 勢能模型的選取

分子動力學模擬結(jié)果的準確性，主要取決于描述分子間作用力的的勢能函數(shù);因此，選擇合適的勢函數(shù)至關(guān)重要。勢函數(shù)可分為對勢與多體勢。本文采用L－J(12－6)對勢函數(shù)，如式(3)所示。

式中:uij為對勢勢能，J;ε能量參數(shù)，J;σ長度參數(shù)，m;rij為粒子i與j之間的距離，m。在本文中，流體與固體壁面之間的L－J勢能函數(shù)如式(4)所示。

式中:εls、σls分別表示流體粒子和固體壁面粒子之間的能量參數(shù)和距離參數(shù)，分別由式(5)和(6)確定。

式中:下標l為氬流體，s為固體壁面，ls為固體壁面與氬流體的相互作用勢能參數(shù)，α為壁面粒子與流體粒子相互作用能量系數(shù)(本文簡稱為能量系數(shù))，β為壁面粒子與流體粒子相互作用長度系數(shù)，本文β取1.03來研究壁面粒子能量系數(shù)α對氬流體在受限空間中自擴散系數(shù)的影響［6］。氬流體的σAr－Ar=0.3405 nm，εAr－Ar/kB=119.8 K。

1.4 模擬體系的運算

在分子動力學模擬中，需要對各物理量進行無量綱化，如表1所示。表中，上標*代表無量綱的物理量;kB為波爾茲曼常數(shù)，1.3810－23J/K;m 為質(zhì)量，kg;T為溫度，K;ρ為密度，kg/m3;H為狹縫寬度，m。初始速度由速度發(fā)生器隨機給定，控溫方法為Nose－Hoover控溫法。時間步長為2fs，總模擬時間400 ps，總步數(shù)為200000步，每1000步輸出一個數(shù)據(jù)。本文模擬數(shù)據(jù)均采用LAMMPS(Largescale Atomic/Molecular Massively Parallel Simulator)軟件計算得到。

表1 物理量的無量綱化

2 模擬結(jié)果與討論

2.1 非受限空間的氬流體自擴散系數(shù)

2.1.1 截斷半徑對自擴散系數(shù)的影響

在粒子數(shù)Num=864，溫度T=90 K條件下，對截斷半徑 rc=0.85、1.02、1.19、1.36、1.53、1.70 nm進行分子動力學模擬，得到粒子的均方位移(MSD)，再由均方位移計算得到氬流體的自擴散系數(shù)(每個條件進行三次模擬計算，并取其平均值)，結(jié)果如表2所示。從表2可以看出，所有模擬值與實驗值的相對誤差都小于7%。其中rc=1.36 nm時，模擬值與實驗值最接近，相對誤差為0.41%，因此，下文均取截斷半徑rc=1.36 nm進行模擬計算。

表2 截斷半徑對氬流體自擴散系數(shù)的影響

2.1.2 粒子數(shù)對自擴散系數(shù)的影響

在T=100 K，截斷半徑rc=1.36 nm的條件下，粒子數(shù) Num=256、500、864、1372 和2048 時的氬流體自擴散系數(shù)如表3所示。由表3可知，模擬平均值與實驗值的相對誤差均小于10%，且隨著粒子數(shù)增大，相對誤差的絕對值逐漸減小，由此可知，隨粒子數(shù)增加模擬結(jié)果越接近于實驗值。但模擬耗時隨著粒子數(shù)增多急劇增大，綜合考慮計算時間、計算機內(nèi)存和相對誤差，下文中的模擬粒子數(shù)均采用864。

表3 粒子數(shù)對氬流體自擴散系數(shù)的影響

2.1.3 溫度對自擴散系數(shù)的影響

在粒子數(shù)Num=864，截斷半徑rc=1.36 nm，溫度 T=90、100、110、120、130 K 時，氬流體的自擴散系數(shù)模擬結(jié)果如表4所示。由表4可見，隨著溫度升高，自擴散系數(shù)逐漸增大。各溫度的自擴散系數(shù)模擬值和實驗值基本吻合，除130 K外相對誤差都在10%以內(nèi)，但越接近于臨界狀態(tài)，自擴散系數(shù)的模擬誤差越大。

表4 不同溫度下氬自擴散系數(shù)模擬結(jié)果

可見，氬流體的自擴散系數(shù)隨溫度的升高而增大。假定氬流體的自擴散系數(shù)與溫度的關(guān)系服從Arrhenius方程(8)［7］。

式中，D0為極限自擴散系數(shù)，m2/s;E為擴散活化能，J/mol;R氣體常數(shù)。

對式(8)取對數(shù)得到式(9)。可見，lnD與1/T呈線性關(guān)系，其斜率為－E/R，截距為1nD0。如圖2所示，直線斜率為－E/R=－3.695，計算得到自擴散活化能 E 為3072.19 J/mol;截距 1nD0=4.945，則其極限自擴散系數(shù)D0為139.77×10－9m/s2。因此，可以認為氬流體的自擴散系數(shù)與溫度的關(guān)系符合Arrhenius方程。

圖2 lnD與1/T關(guān)系曲線

2.1.4 徑向分布函數(shù)

徑向分布函數(shù)g(r)可以直接反映物質(zhì)分子的聚集特性。在粒子數(shù)Num=864，rc=1.36 nm條件下，對 T=90、100、110、120、130 K 的氬流體體系進行模擬，得到粒子的徑向分布函數(shù)圖像如圖3所示。由圖3可以看出，當r＞4.0σ時，各圖中徑向分布函數(shù)的值都趨近于1，r＞0.9σ時徑向分布函數(shù)的值均為零，這表示該流體系統(tǒng)中分子與分子之間的最近距離不可能小于0.9σ，而且當分子間的距離超過4.0σ時，與均勻流體性質(zhì)相同。此外，徑向分布函數(shù)的最大峰值均出現(xiàn)在r約為1.2σ處，這表明離分子1.2σ處附近出現(xiàn)其他分子的概率最大，此處的區(qū)域密度最大，但峰值卻隨著溫度升高而減小。

圖3 不同溫度下氬流體的徑向分布函數(shù)

2.2 受限空間的氬流體自擴散系

2.2.1 能量系數(shù)對自擴散系數(shù)的影響

當狹縫寬度 H=2.99nm，rc=1.36 nm，能量系數(shù) α =0.1、0.6、1.0、3.0、5.0、7.0 時，模擬得到的氬流體自擴散系數(shù)如圖 4所示(T=90、100、110、120 K)。由圖4可見，隨著溫度增大，氬流體自擴散系數(shù)逐漸增大，但隨著能量系數(shù)的增大，自擴散系數(shù)逐漸減小。

圖4 能量系數(shù)α對自擴散系數(shù)的影響

2.2.2 狹縫寬度對自擴散系數(shù)的影響

在溫度T=100 K，rc=1.36 nm，狹縫寬度 H=2.99、4.11、5.28、6.45、8.21 nm 的條件下，模擬得到的氬流體自擴散系數(shù)如圖5所示。從圖中可以看出，隨著狹縫寬度的增大，氬流體自擴散系數(shù)逐漸增大，且增大趨勢漸趨于平緩，最終接近于非受限空間中的自擴散系數(shù)。

圖5 狹縫寬度對氬流體自擴散系數(shù)的影響

2.2.3 溫度對自擴散系數(shù)的影響

在狹縫寬度 H=2.99 nm，rc=1.36 nm 的受限空間中，當溫度 T=90、100、110、120、130 K 時，氬流體自擴散系數(shù)的模擬結(jié)果如圖6所示，并與非受限空間中自擴散系數(shù)的模擬值及實驗值進行了比較。從圖6可見，在受限條件下，氬流體自擴散系數(shù)的模擬值與實驗值都隨著溫度升高而增大，氬流體自擴散系數(shù)模擬值隨溫度的變化也滿足Arrhenius關(guān)系，計算得到的自擴散活化能E=2362.85 J/mol，極限自擴散系數(shù)D0=51.21×10－9m2/s。但在相同溫度的條件下，受限空間氬流體的自擴散系數(shù)比非受限空間的要低。

圖6 溫度對受限空間氬流體自擴散系數(shù)的影響

3 結(jié)論

本文采用分子動力學模擬方法，考察了非受限空間中截斷半徑、粒子數(shù)、溫度等和受限空間中能量系數(shù)、狹縫寬度及溫度等對氬流體自擴散系數(shù)的影響，并得到了不同溫度下氬流體的徑向分布函數(shù)。通過對模擬結(jié)果的分析，得到以下結(jié)論:

在非受限空間中，隨溫度T升高，自擴散系數(shù)增大，且符合Arrhenius方程，自擴散活化能為 3072.19 J/mol。

在受限空間中，隨著能量系數(shù)α的增大，氬流體的自擴散系數(shù)逐漸減小;隨著溫度的升高，氬流體的自擴散系數(shù)增大，且滿足Arrhenius方程，擴散活化能為2362.85 J/mol;隨著狹縫寬度的增加，氬流體的自擴散系數(shù)也逐漸增大，最后趨于實驗值;在相同溫度的條件下，受限空間氬流體的自擴散系數(shù)比非受限空間的要低。

［1］Allen M P，Tildesley D J.Computer simulation of liquids［M］.Oxford:Oxford University Press，1987.

［2］Meier K，Laesecke A，Kabela S C.A molecular dynamics simulation study of the self－diffusion coefficient and viscosity of the Lennard － Jones fluid［J］.International Journal of Thermophysics，2001，22(1):161 －173.

［3］嚴六明，嚴琪良，劉洪來.Ar－Kr溶液擴散系數(shù)的分子動力學模擬及其與溫度的關(guān)系［J］.高等學?；瘜W學報，1997，18:2026 －2029.

［4］葛 宋，陳 民.超臨界L－J流體混合物互擴散性質(zhì)的分子動力學模擬［J］.高等化學學報，2011，32:2593－2597.

［5］曹炳陽.速度滑移及其對微納尺度流動影響的分子動力學研究［D］.北京:清華大學航天航空學院，2005.

［6］秦 星，張秉堅，張 暉，等.微孔中簡單流體混合物自擴散系數(shù)的分子動力學模擬與關(guān)聯(lián)［J］.化學物理學報，2003，16(6):467 －471.

［7］方俊鑫，陸 棟.固體物理學［M］.上海:科學技術(shù)出版社，1981.