L_J流體自擴散系數(shù)及其與溫度關(guān)系的分子動力學模擬

1、第 5卷第 2期2006年 6月熱 科 學 與 技 術(shù)Journal of Thermal Science and T echnologyV ol. 5N o. 2Jun. 2006文章編號 :167128097(2006 0220101205收稿日期 :2006203202; 修回日期 :2006205208.基金項目 :國家自然科學基金委資助項目 (50475100 .作者簡介 :李維仲 (19562 , 男 , 教授 , 博士 , 博士生導師 , 主要研究方向為計算流體力學 .L 2J 流體自擴散系數(shù)及其與溫度關(guān)系的分子動力學模擬李 維 仲 1, 陳 聰 1, 楊 健2(1. 大連理工大

2、學 能源與動力學院 , 遼寧 大連 116024; 2. 浙江大學 材料與化學工程學院 , 浙江 杭州 310027摘要 :擴散系數(shù)在化工設(shè)計和研究中是不可缺少的傳遞特性 。但其數(shù)據(jù)卻相對缺乏 , 因此需要尋找一種方法 來預(yù)測這個特性就顯得十分重要 。 利用分子動力學方法模擬了簡單流體的自擴散系數(shù) 。 模擬分別采用G reen 2K ubo 法 (VACF :velocity autocorrelation function 和 Einstein 法 (MSD :mean square displacement 。 模擬結(jié)果 與實驗數(shù)據(jù)吻合較好 , 誤差在 10%左右 。 兩種方法的平均值與實

3、驗結(jié)果誤差在 7%左右 。 同時還模擬了流體自 擴散系數(shù)隨溫度的變化關(guān)系 。 結(jié)果表明 , 自擴散系數(shù)與溫度滿足 Arrenhius 關(guān)系 , 數(shù)據(jù)相關(guān)性在 0. 99以上 , 計算 得到的自擴散激活能分別為 1258J m ol (VACF 、 1272J m ol (MSD 和平均值 1265J m ol 。關(guān)鍵詞 :L 2J 流體 ; 分子動力學 ; 自擴散系數(shù) ; ; 中圖分類號 :Q2文獻標識碼 :A0引 擴散系數(shù)是化學工程中設(shè)計和研究不可缺少 的數(shù)據(jù) 。 分子動力學方法最近幾十年來獲得了極 大的發(fā)展 , 理論模型已經(jīng)比較成熟 , 已經(jīng)成為在分子水平上進行數(shù)值模擬的重要手段 1。 利

4、用平衡 分子動力學方法模擬擴散系數(shù) , 主要通過兩個關(guān)系式 :Einstein 關(guān) 系 式 和 G reen 2K ubo 關(guān) 系 式 2。 Meier 等3模擬了 Lennard 2Jones 流體的自擴散系數(shù) , 獲得了較好的結(jié)果 , 但只是采用了 Einstein 關(guān)系 ; 孫煒等 4用兩種方法模擬了簡單流體的自擴 散系數(shù) , 兩種方法得到的結(jié)果相同 , 但其誤差偏 大 , 最大誤差達到 11. 8%。 本文分別利用兩種方 法模擬自擴散系數(shù) , 并考察了自擴散系數(shù)隨溫度 變化情況 。1模擬方法1. 1基本理論經(jīng)典分子動力學方法可以分為平衡分子動力 學方法和非平衡分子動力學方法 , 本文采

5、用平衡 分子動力學方法 。 1. 1. 1 Einstein 關(guān)系D =limt 6N m tNmj =1r j (t -r j (0 2(1式中 :N m 為模擬粒子數(shù)目 ; t 為模擬時間 ; 為 表示系統(tǒng)平均 ; r j (t 為 t 時刻第 j 個粒子的真實 位移 , 由于 MD 模擬中采用周期性邊界條件 , 必須 變換得到粒子的真實位移。 假設(shè)每個時間步長內(nèi) 粒子的位移遠遠小于系統(tǒng)尺寸 (實際模擬中往往可以滿足 , 則模擬粒子的真實坐標可以表示為 5r jx (t =r jx (t +nint (r jx (t -t -r jx (t L x L x其中 :r jx (t 為粒子真實

6、位移 , r jx (t 為 MD 程序中粒子位移 ,nint (x 為最接近于 x 的整數(shù) , r j (0 =r j (0 , L x 為系統(tǒng)方向上的尺寸 , t 為時間步長。 1. 1. 2 G reen 2K ubo 關(guān)系D =3N m Nmj =1v j(t v j(0 d t(2式中 :v j (t 為 t 時刻第 j 個粒子的速度 , 其他參數(shù) 與 Einstein 關(guān)系式相同 。擴散系數(shù)的 Einstein 關(guān)系式和 G reen 2K ubo 關(guān) 系式 在 理 論 上 是 等 價 。 Einstein 關(guān) 系 式 中 的 Nmj =1r j (t -r j (0 2項被稱為均

7、方根位移 (meansquare displacement ,MSD ,G reen 2K ubo 關(guān)系式中的 Nmj =1v j(t v j (0 項被稱為 速 度 自 適 應(yīng) 函 數(shù)(velocity autocorrelation function ,VACF , 都有 很重要的研究價值 6。1. 2模擬體系模擬對象為三維空間內(nèi)的 500個氬原子 , 采 用正則系統(tǒng) (N ,V ,T 不變 , 三個方向均采用周期 性邊界條件 , 利用 Verlet 算法求解粒子運動方程 , 采用鄰近表法計算粒子間作用力 , 取截斷半徑 2. 5, 鄰近表半徑 0. 3, 每 5個步長更新一次鄰近表 。

8、初始構(gòu)型采用 FCC 結(jié)構(gòu) , 時間步長為 210-15s , 模擬總步數(shù) 110000, 前 10000步采用溫度 速 度調(diào)節(jié)使系統(tǒng)達到平衡 , 后 100000步用于統(tǒng)計計 算流體各種性質(zhì)。為了計算方便 , 模擬過程中采用無量綱化參 數(shù) , 即T =kT , V =V 3, p p 3, 3Ar (Ar , 其中 :k B =5nm , 氬的相對原子質(zhì)量 Ar (Ar =39. 948, 時間采用無量綱單位t =210-12, 本文中未指明時 , 均為無量綱單位。1. 3模擬程序LAMMPS (large 2scale atomic m olecularmassively parallel

9、 simulator 是新近出現(xiàn)的一種分 子動力學模擬軟件 , 不僅適用于大尺度并行計算 , 也非常適用于單個處理器計算。 本文數(shù)據(jù)全部在 奔騰 4PC 上利用 LAMMPS 計算得到 。2模擬結(jié)果2. 1模擬參數(shù)選擇 2. 1. 1模擬粒子數(shù)目為了比較不同粒子數(shù)目對模擬的影響 , 考察了密度為 1160kg m 3, 溫度 120K, 模擬粒子從 108到 864, 結(jié)果 1所示 , , 粒子總能 , , 擴散系數(shù)模擬 , 但模擬花費時間大幅度增加。 綜合考慮 , 本文取模擬粒子數(shù)目 500。表 1不同粒子數(shù)目時系統(tǒng)位能、 自擴散系數(shù)和時間消耗比較Tab 11 C omparison of

10、system potential energy , self 2diffusion coefficient and time cost under different atom numbersND VACF (10-9m 2 s -1D MS D (10-9m 2 s -1D AV (10-9m 2 s -1UUt mVACFMS DAV1085. 595. 395. 491. 488590. 0063631-7. 7611. 06-9. 412565. 865. 645. 751. 496630. 003191. 8-3. 36. 93-5. 125006. 516. 136. 321. 4

11、99500. 0023145. 67. 431. 164. 298646. 3986. 0596. 2351. 500770. 0015318. 45. 58-0. 022. 78注 :D VACF 、 D MS D 分別表示由 G reen 2K ubo 關(guān)系和 Einstein 關(guān)系求得的擴散系數(shù) ; D AV 為 D VACF 和 D MS D 的算術(shù)平均值。 U 為每個模擬粒子的 總能量 ; U 為模擬粒子總能量的波動值 ; t m 為模擬花費時間 , 表中數(shù)據(jù)為相對于粒子數(shù)為 108時的相對值 ; U 、 U 、 t m 均為無量綱數(shù) ; VACF 、 MS D 、 AV 為三個擴散

12、系數(shù)對應(yīng)的百分比誤差。2. 1. 2松弛時間和截斷半徑約化的速度自適應(yīng)函數(shù)為 : Nmj =1v j(t v j(0 Nmj =1v j(0 2圖 1為密度 1374kg m 3, 溫度 90K 時的速度自適應(yīng)函數(shù) 。 由圖 1可見 , 隨著時間演變 , 速度自 適應(yīng)函數(shù)值迅速下降 , 在 0. 4ps (圖 1中 t =0. 2 時 , 已變得很小 , 表明此時速度與前面速度的相關(guān) 性已很小 , 在 2ps (t =1 時基本上已經(jīng)達到零 , 隨 后在零附近波動 , 這段時間也稱為松弛時間。 計算 擴散系數(shù)時 , 每個種子取樣最長時間必須大于松 弛時間 。 為了保證結(jié)果的精確性 , 本文中最

13、長取樣 時間為 24ps 。圖 1速度自適應(yīng)函數(shù) (X VACF 隨時間變化關(guān)系Fig 11 Variation of velocity autocorrelation functionwith time L 2J 勢能曲線如圖 2所示 , 可見 , 粒子間距離 r201 熱 科 學 與 技 術(shù) 第 5卷 圖 2 L 2J 勢能 (u 隨原子間距變化關(guān)系Fig 12 Variation of L 2J potential energy with distance=2. 0時 , 勢能已接近零 , 本文取截斷半徑為 2. 5。2. 2自擴散系數(shù)模擬模擬了流體氬在若干狀態(tài)點的自擴散系數(shù) ,并與文獻

14、 7的實驗值和文獻 4的模擬值進行 了比較 。 模擬狀態(tài)點及模擬結(jié)果見表 2。由表 2可見 , 氬的自擴散系數(shù)模擬值與文獻 7的實驗值基本吻合 , 兩種方法模擬結(jié)果誤差 均在 10%左右 。 而對兩種方法獲得的結(jié)果求平 均 , 則可以有效的減小誤差 , 均在 7%左右 。表 2不同狀態(tài)點流體氬自擴散系數(shù)的模擬值與實驗值比較Tab 12 C omparison between simulation and experimental results of self 2diffusion coefficient of argon on different statesT Kp MPa (kg m -

15、3 D EXP D VACF D MS D D AV D M DVACFMS DAVM D900. 13413742. 432. 202. 412. 3052. 209. 470. 5. 14-9. 421000. 32513093. 543. 393. 163. 275257. 49-8. 191100. 66712384. 804. 574. 364. 519. 6. 98-6. 041201. 21311606. 066. 516. 13-1. 16-4. 291. 001302. 02510657. 457. -10. 47-4. 16 -7. 3211. 80注 :D EXP 為擴散

16、系數(shù)實驗值 VACF MS D 、 D AV 、 D M D 的單位均為 10-9m 2s 。2. 3 自擴散系數(shù)與溫度存在 Arrhenius 關(guān)系8:D =D 0expRT(3式中 :D 0為溫度 T 時的極限擴散系數(shù) ; E 為 擴散激活能 ; R 為氣體常數(shù)。 可見 , 擴散系數(shù)的自 然對數(shù)與溫度的倒數(shù)成線性關(guān)系。模擬了密度 1374kg m 3, 溫度從 83. 86K (無量綱值 ,0. 7 到 143. 76K (無量綱值 ,1. 2 時的自 擴散系數(shù) 。 模擬結(jié)果如圖 35所示 。 圖 3、 4分別 圖 3擴散系數(shù)與溫度關(guān)系 (VACF Fig 13 Relation ofse

17、lf 2diffusioncoefficientwithtemperature (VACF 圖 4擴散系數(shù)與溫度關(guān)系 (MSD Fig 14 Relation ofself 2diffusioncoefficientwithtemperature (MSD 圖 5擴散系數(shù)與溫度關(guān)系 (平均值 Fig 15 Relation ofself 2diffusioncoefficientwithtemperature (average 為利用 G reen 2K ubo 關(guān)系和 Einstein 關(guān)系獲得的結(jié)301 第 2期 李維仲等 :L 2J 流體自擴散系數(shù)及其與溫度關(guān)系的分子動力學模擬果 , 圖

18、 5則是通過對兩種方法模擬結(jié)果取算術(shù)平 均值 得 到 。 三 種 方 法 的 數(shù) 據(jù) 相 關(guān) 系 數(shù) 分 別 為 0. 997、 0. 999和 0. 998, 計算得到的自擴散激活能 大小分別為 1258、 1272和 1265J m ol 。3結(jié) 論本文分別用兩種方法模擬了模型流體在若干 狀態(tài)點下的自擴散系數(shù) , 同時還模擬了擴散系數(shù) 隨溫度的變化關(guān)系 , 并與前人的實驗結(jié)果和模擬 結(jié)果進行了比較 。1 利用兩種方法均能獲得滿意的結(jié)果 , 兩種 方法模擬結(jié)果誤差均在 10%左右 , 而對兩種方法 獲得的結(jié)果求平均 , 則可以有效的將誤差減少到 7%左右 。 與文獻 4相比 , 在模擬總時長

19、基本一 致的情況下 , 本文平均誤差更小 。 通過對比發(fā)現(xiàn) , 本文中采用的步長為文獻 4的 0. 2為其 12倍 。,零 , 。2 利用 VACF 法 、 MSD 法以及兩種方法的平 均值求得的流體自擴散激活能大小分別為 1258、 1272和 1265, 與文獻 9中的 1263相比 , 誤差分 別為 0. 4%、 0. 7%和 0. 16%。 而文獻 10的模擬 結(jié)果為 1294, 誤差為 2. 45%。 顯然 , 本文的模擬結(jié) 果誤差較小 。3 在實驗數(shù)據(jù)缺乏的情況下 , 可以利用分子 動力學方法來模擬流體的自擴散系數(shù)和自擴散激 活能 。參考文獻 (R eferences :1楊小震

20、. 分子模擬與高分子材料 M.北京 :科學出 版社 , 2003. 2McQUARRIE D A. Statistical mechanics M.New Y ork :Harper and Row , 1976.3MEIER K, LAESECKE A , K ABE LAC S. A m olecular dynamics simulation study of the self 2duffusion coefficient and viscosity of the lennard 2jones fluid J.Int J of Thermophys , 2001, 22(1 :161.4

21、孫 煒 , 陳 中 , 黃素逸 , 等 . L 2J 流體擴散系數(shù)的分 子動力學模擬 J.華中科技大學學報 :自然科學版 , 2004, 32(5 :84.SUN Wei , CHEN Zhong , HUANG Su 2yi , et al . M olecular dynamics simulation of Lennard 2Jones fluid diffusion coefficients J.J of Huazhong Univ of Sci and Tech :Natural Sci , 2004, 32(5 :84.(in Chinese 5RAPAPORT D C. The

22、art m olecular dynamics M.:Univ. Press ,D J. C omputer simulation of M.Ox ford :Clarendon Press , 1987.7NAGHIZ ADEH J , RICE S A. K inetic theory of dense fluids. X. measurement and interpretation of self 2diffusion in liquid Ar , K r , Xe , and CH4J.The J of Chem Phys , 1962, 36(10 :2710.8方俊鑫 , 陸 棟

23、 . 固體物理學 M.上海 :科學技術(shù) 出版社 , 1980.9ROW LEY R L , PAINTER M M. Diffusion and viscosity equations of state for a Lennard 2Jones fluid obtained from m olecular dynamics simulations J.Int J of Thermophys , 1997, 18(5 :1109.10孫 煒 , 陳 中 , 黃素逸 . 模型流體擴散系數(shù)與溫 度關(guān)系的分子動力學模擬 J.武漢化工學院學報 , 2005, 27(4 :124.SUN Wei , CH

24、EN Zhong , HUANG Su 2yi. M olecular dynamics simulation of Lennar 2Jones fluid :calculation of diffusion coefficients J.J of Wuhan Inst of Chem Tech , 2005, 27(4 :124. (in Chinese 401 熱 科 學 與 技 術(shù) 第 5卷 Molecular dynamics simulation of self 2diffusion coefficient and its relation with temperature usin

25、g simple Lennard 2Jones potentialLI Wei 2zhong 1, CHEN Cong 1, Y AN G Jian2(1. School of Energy and P ower Eng. , Dalian Univ. of T ech. , Dalian 116024, China ; 2. C ollege of Materials Sci. and Chem. Eng. , Zhejiang Univ. , Hangzhou 310027, China Abstract :Diffusion coefficient is indispensable to chemical engineering design and research. H owever , the