L_J流体自扩散系数及其与温度关系的分子动力学模拟

1、第 5卷第 2期2006年 6月热 科 学 与 技 术Journal of Thermal Science and T echnologyV ol. 5N o. 2Jun. 2006文章编号 :167128097(2006 0220101205收稿日期 :2006203202; 修回日期 :2006205208.基金项目 :国家自然科学基金委资助项目 (50475100 .作者简介 :李维仲 (19562 , 男 , 教授 , 博士 , 博士生导师 , 主要研究方向为计算流体力学 .L 2J 流体自扩散系数及其与温度关系的分子动力学模拟李 维 仲 1, 陈 聪 1, 杨 健2(1. 大连理工大

2、学 能源与动力学院 , 辽宁 大连 116024; 2. 浙江大学 材料与化学工程学院 , 浙江 杭州 310027摘要 :扩散系数在化工设计和研究中是不可缺少的传递特性 。但其数据却相对缺乏 , 因此需要寻找一种方法 来预测这个特性就显得十分重要 。 利用分子动力学方法模拟了简单流体的自扩散系数 。 模拟分别采用G reen 2K ubo 法 (VACF :velocity autocorrelation function 和 Einstein 法 (MSD :mean square displacement 。 模拟结果 与实验数据吻合较好 , 误差在 10%左右 。 两种方法的平均值与实

3、验结果误差在 7%左右 。 同时还模拟了流体自 扩散系数随温度的变化关系 。 结果表明 , 自扩散系数与温度满足 Arrenhius 关系 , 数据相关性在 0. 99以上 , 计算 得到的自扩散激活能分别为 1258J m ol (VACF 、 1272J m ol (MSD 和平均值 1265J m ol 。关键词 :L 2J 流体 ; 分子动力学 ; 自扩散系数 ; ; 中图分类号 :Q2文献标识码 :A0引 扩散系数是化学工程中设计和研究不可缺少 的数据 。 分子动力学方法最近几十年来获得了极 大的发展 , 理论模型已经比较成熟 , 已经成为在分子水平上进行数值模拟的重要手段 1。 利

4、用平衡 分子动力学方法模拟扩散系数 , 主要通过两个关系式 :Einstein 关 系 式 和 G reen 2K ubo 关 系 式 2。 Meier 等3模拟了 Lennard 2Jones 流体的自扩散系数 , 获得了较好的结果 , 但只是采用了 Einstein 关系 ; 孙炜等 4用两种方法模拟了简单流体的自扩 散系数 , 两种方法得到的结果相同 , 但其误差偏 大 , 最大误差达到 11. 8%。 本文分别利用两种方 法模拟自扩散系数 , 并考察了自扩散系数随温度 变化情况 。1模拟方法1. 1基本理论经典分子动力学方法可以分为平衡分子动力 学方法和非平衡分子动力学方法 , 本文采

5、用平衡 分子动力学方法 。 1. 1. 1 Einstein 关系D =limt 6N m tNmj =1r j (t -r j (0 2(1式中 :N m 为模拟粒子数目 ; t 为模拟时间 ; 为 表示系统平均 ; r j (t 为 t 时刻第 j 个粒子的真实 位移 , 由于 MD 模拟中采用周期性边界条件 , 必须 变换得到粒子的真实位移。 假设每个时间步长内 粒子的位移远远小于系统尺寸 (实际模拟中往往可以满足 , 则模拟粒子的真实坐标可以表示为 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 的整数 , r j (0 =r j (0 , L x 为系统方向上的尺寸 , t 为时间步长。 1. 1. 2 G reen 2K ubo 关系D =3N m Nmj =1v j(t v j(0 d t(2式中 :v j (t 为 t 时刻第 j 个粒子的速度 , 其他参数 与 Einstein 关系式相同 。扩散系数的 Einstein 关系式和 G reen 2K ubo 关 系式 在 理 论 上 是 等 价 。 Einstein 关 系 式 中 的 Nmj =1r j (t -r j (0 2项被称为均

7、方根位移 (meansquare displacement ,MSD ,G reen 2K ubo 关系式中的 Nmj =1v j(t v j (0 项被称为 速 度 自 适 应 函 数(velocity autocorrelation function ,VACF , 都有 很重要的研究价值 6。1. 2模拟体系模拟对象为三维空间内的 500个氩原子 , 采 用正则系统 (N ,V ,T 不变 , 三个方向均采用周期 性边界条件 , 利用 Verlet 算法求解粒子运动方程 , 采用邻近表法计算粒子间作用力 , 取截断半径 2. 5, 邻近表半径 0. 3, 每 5个步长更新一次邻近表 。

8、初始构型采用 FCC 结构 , 时间步长为 210-15s , 模拟总步数 110000, 前 10000步采用温度 速 度调节使系统达到平衡 , 后 100000步用于统计计 算流体各种性质。为了计算方便 , 模拟过程中采用无量纲化参 数 , 即T =kT , V =V 3, p p 3, 3Ar (Ar , 其中 :k B =5nm , 氩的相对原子质量 Ar (Ar =39. 948, 时间采用无量纲单位t =210-12, 本文中未指明时 , 均为无量纲单位。1. 3模拟程序LAMMPS (large 2scale atomic m olecularmassively parallel

9、 simulator 是新近出现的一种分 子动力学模拟软件 , 不仅适用于大尺度并行计算 , 也非常适用于单个处理器计算。 本文数据全部在 奔腾 4PC 上利用 LAMMPS 计算得到 。2模拟结果2. 1模拟参数选择 2. 1. 1模拟粒子数目为了比较不同粒子数目对模拟的影响 , 考察了密度为 1160kg m 3, 温度 120K, 模拟粒子从 108到 864, 结果 1所示 , , 粒子总能 , , 扩散系数模拟 , 但模拟花费时间大幅度增加。 综合考虑 , 本文取模拟粒子数目 500。表 1不同粒子数目时系统位能、 自扩散系数和时间消耗比较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 关系和 Einstein 关系求得的扩散系数 ; D AV 为 D VACF 和 D MS D 的算术平均值。 U 为每个模拟粒子的 总能量 ; U 为模拟粒子总能量的波动值 ; t m 为模拟花费时间 , 表中数据为相对于粒子数为 108时的相对值 ; U 、 U 、 t m 均为无量纲数 ; VACF 、 MS D 、 AV 为三个扩散

12、系数对应的百分比误差。2. 1. 2松弛时间和截断半径约化的速度自适应函数为 : Nmj =1v j(t v j(0 Nmj =1v j(0 2图 1为密度 1374kg m 3, 温度 90K 时的速度自适应函数 。 由图 1可见 , 随着时间演变 , 速度自 适应函数值迅速下降 , 在 0. 4ps (图 1中 t =0. 2 时 , 已变得很小 , 表明此时速度与前面速度的相关 性已很小 , 在 2ps (t =1 时基本上已经达到零 , 随 后在零附近波动 , 这段时间也称为松弛时间。 计算 扩散系数时 , 每个种子取样最长时间必须大于松 弛时间 。 为了保证结果的精确性 , 本文中最

13、长取样 时间为 24ps 。图 1速度自适应函数 (X VACF 随时间变化关系Fig 11 Variation of velocity autocorrelation functionwith time L 2J 势能曲线如图 2所示 , 可见 , 粒子间距离 r201 热 科 学 与 技 术 第 5卷 图 2 L 2J 势能 (u 随原子间距变化关系Fig 12 Variation of L 2J potential energy with distance=2. 0时 , 势能已接近零 , 本文取截断半径为 2. 5。2. 2自扩散系数模拟模拟了流体氩在若干状态点的自扩散系数 ,并与文献

14、 7的实验值和文献 4的模拟值进行 了比较 。 模拟状态点及模拟结果见表 2。由表 2可见 , 氩的自扩散系数模拟值与文献 7的实验值基本吻合 , 两种方法模拟结果误差 均在 10%左右 。 而对两种方法获得的结果求平 均 , 则可以有效的减小误差 , 均在 7%左右 。表 2不同状态点流体氩自扩散系数的模拟值与实验值比较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、系数实验值 VACF MS D 、 D AV 、 D M D 的单位均为 10-9m 2s 。2. 3 自扩散系数与温度存在 Arrhenius 关系8:D =D 0expRT(3式中 :D 0为温度 T 时的极限扩散系数 ; E 为 扩散激活能 ; R 为气体常数。 可见 , 扩散系数的自 然对数与温度的倒数成线性关系。模拟了密度 1374kg m 3, 温度从 83. 86K (无量纲值 ,0. 7 到 143. 76K (无量纲值 ,1. 2 时的自 扩散系数 。 模拟结果如图 35所示 。 图 3、 4分别 图 3扩散系数与温度关系 (VACF Fig 13 Relation ofse

17、lf 2diffusioncoefficientwithtemperature (VACF 图 4扩散系数与温度关系 (MSD Fig 14 Relation ofself 2diffusioncoefficientwithtemperature (MSD 图 5扩散系数与温度关系 (平均值 Fig 15 Relation ofself 2diffusioncoefficientwithtemperature (average 为利用 G reen 2K ubo 关系和 Einstein 关系获得的结301 第 2期 李维仲等 :L 2J 流体自扩散系数及其与温度关系的分子动力学模拟果 , 图

18、 5则是通过对两种方法模拟结果取算术平 均值 得 到 。 三 种 方 法 的 数 据 相 关 系 数 分 别 为 0. 997、 0. 999和 0. 998, 计算得到的自扩散激活能 大小分别为 1258、 1272和 1265J m ol 。3结 论本文分别用两种方法模拟了模型流体在若干 状态点下的自扩散系数 , 同时还模拟了扩散系数 随温度的变化关系 , 并与前人的实验结果和模拟 结果进行了比较 。1 利用两种方法均能获得满意的结果 , 两种 方法模拟结果误差均在 10%左右 , 而对两种方法 获得的结果求平均 , 则可以有效的将误差减少到 7%左右 。 与文献 4相比 , 在模拟总时长

19、基本一 致的情况下 , 本文平均误差更小 。 通过对比发现 , 本文中采用的步长为文献 4的 0. 2为其 12倍 。,零 , 。2 利用 VACF 法 、 MSD 法以及两种方法的平 均值求得的流体自扩散激活能大小分别为 1258、 1272和 1265, 与文献 9中的 1263相比 , 误差分 别为 0. 4%、 0. 7%和 0. 16%。 而文献 10的模拟 结果为 1294, 误差为 2. 45%。 显然 , 本文的模拟结 果误差较小 。3 在实验数据缺乏的情况下 , 可以利用分子 动力学方法来模拟流体的自扩散系数和自扩散激 活能 。参考文献 (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 流体扩散系数的分 子动力学模拟 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.上海 :科学技术 出版社 , 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孙 炜 , 陈 中 , 黄素逸 . 模型流体扩散系数与温 度关系的分子动力学模拟 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 热 科 学 与 技 术 第 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