N DVACFΠ(10 - 9 m2 ·s- 1)
108
5. 59
DMSDΠ(10 - 9 m2 ·s- 1) 5. 39
DAVΠ(10 - 9 m2 ·s- 1) 5. 49
U 1. 488 59
ΔU 0. 006 363
tm ηVACF ηMSD
ηAV
1 - 7. 76 11. 06 - 9. 41
表 2 不同状态点流体氩自扩散系数的模拟值与实验值比较
Tab12 Comparison between simulation and experimental results of self2diffusion coefficient of argon on different states
TΠK
注 : DVACF 、DMSD 分别表示由 Green2Kubo 关系和 Einstein 关系求得的扩散系数 ; DAV 为 DVACF 和 DMSD 的算术平均值。U 为每个模拟粒子的 总能量 ;ΔU 为模拟粒子总能量的波动值 ; tm 为模拟花费时间 ,表中数据为相对于粒子数为 108 时的相对值 ; U 、ΔU 、tm 均为无量纲数 ; ηVACF 、ηMSD 、ηAV 为三个扩散系数对应的百分比误差。
模拟了流体氩在若干状态点的自扩散系数 , 并与文献[ 7 ] 的实验值和文献 [ 4 ] 的模拟值进行 了比较源自。模拟状态点及模拟结果见表 2 。
由表 2 可见 ,氩的自扩散系数模拟值与文献 [ 7 ] 的实验值基本吻合 , 两种方法模拟结果误差 均在 10 % 左右 。而对两种方法获得的结果求平 均 ,则可以有效的减小误差 ,均在 7 % 左右 。
1 模拟方法
1. 1 基本理论 经典分子动力学方法可以分为平衡分子动力
1. 1. 1 Einstein 关系
∑ D
=
lim
t →∞
6
1〈
Nm t
Nm j=1
[
rj
(
t)
-
rj (0) ]2 〉
(1)
式中 : Nm 为模拟粒子数目 ; t 为模拟时间〈; …〉为
表示系统平均 ; rj ( t) 为 t 时刻第 j 个粒子的真实
收稿日期 : 2006203202 ; 修回日期 : 2006205208.
基金项目 : 国家自然科学基金委资助项目 (50475100) .
作者简介 : 李维仲 (19562) , 男 , 教授 , 博士 , 博士生导师 , 主要研究方向为计算流体力学.
© 1994-2010 China Academic Journal Electronic Publishing House. All rights reserved.
pΠMPa ρΠ(kg ·m- 3)
D EXP
DVACF
DMSD
DAV
DMD
ηVACF
ηMSD
ηAV
ηMD
90
0. 134
1 374
2. 43
2. 20
2. 41
2. 305
2. 20
9. 47 0. 82 5. 14 - 9. 42
100 0. 325
1 309
3. 54
3. 39
3. 16
3. 275
奔腾 4PC 上利用 LAMMPS 计算得到 。
2 模拟结果
2. 1 模拟参数选择 2. 1. 1 模拟粒子数目
为了比较不同粒子数目对模拟的影响 ,考察 了密度为 1 160 kgΠm3 ,温度 120 K ,模拟粒子从 108 到 864 时系统位能和扩散系数的变化情况 ,结果 如表 1 所示 。
位移 ,由于 MD 模拟中采用周期性边界条件 ,必须
变换得到粒子的真实位移 。假设每个时间步长内
粒子的位移远远小于系统尺寸 (实际模拟中往往
可以满足) ,则模拟粒子的真实坐标可以表示为[5]
′
rjx
(
t)
=
rjx ( t)
+
nint ( [
′
rjx
(
t
-
Δt)
-
rjx ( t) ]ΠL x ) L x
0 j=1
·vj (0) 〉d t
(2)
式中 : vj ( t) 为 t 时刻第 j 个粒子的速度 ,其他参数
学方法和非平衡分子动力学方法 ,本文采用平衡 与 Einstein 关系式相同 。
分子动力学方法 。
扩散系数的 Einstein 关系式和 Green2Kubo 关
系式 在 理 论 上 是 等 价 。Einstein 关 系 式 中 的
李 维 仲1 , 陈 聪1 , 杨 健2
( 1. 大连理工大学 能源与动力学院 , 辽宁 大连 116024 ; 2. 浙江大学 材料与化学工程学院 , 浙江 杭州 310027 )
摘要 : 扩散系数在化工设计和研究中是不可缺少的传递特性 。但其数据却相对缺乏 ,因此需要寻找一种方法
来预测这个特性就 显 得 十 分 重 要 。利 用 分 子 动 力学方法模拟了简单流体的自扩散系数 。模拟分别采用 Green2Kubo 法 (VACF : velocity autocorrelation function) 和 Einstein 法 (MSD : mean square displacement) 。模拟结果 与实验数据吻合较好 ,误差在 10 % 左右 。两种方法的平均值与实验结果误差在 7 % 左右 。同时还模拟了流体自 扩散系数随温度的变化关系 。结果表明 ,自扩散系数与温度满足 Arrenhius 关系 ,数据相关性在0. 99 以上 ,计算 得到的自扩散激活能分别为 1 258 JΠmol (VACF) 、1 272 JΠmol (MSD) 和平均值 1 265 JΠmol 。
3. 25
4. 24 10. 73 7. 49 - 8. 19
110 0. 667
1 238
4. 80
4. 57
4. 36
4. 465
4. 51
4. 79 9. 17 6. 98 - 6. 04
120 1. 213
关键词 : L2J 流体 ; 分子动力学 ; 自扩散系数 ; 数值模拟 ; 激活能 中图分类号 : Q2 文献标识码 : A
0 引 言
扩散系数是化学工程中设计和研究不可缺少 的数据 。分子动力学方法最近几十年来获得了极 大的发展 ,理论模型已经比较成熟 ,已经成为在分 子水平上进行数值模拟的重要手段[1] 。利用平衡 分子动力学方法模拟扩散系数 ,主要通过两个关 系式 : Einstein 关 系 式 和 Green2Kubo 关 系 式[2] 。 Meier 等[3] 模拟了 Lennard2Jones 流体的自扩散系 数 ,获得了较好的结果 ,但只是采用了 Einstein 关 系 ;孙炜等[4] 用两种方法模拟了简单流体的自扩 散系数 , 两种方法得到的结果相同 , 但其误差偏 大 ,最大误差达到 11. 8 %。本文分别利用两种方 法模拟自扩散系数 ,并考察了自扩散系数随温度 变化情况 。
为了计算方便 ,模拟过程中采用无量纲化参 数 ,即 ~T = kTΠε, V~ = VΠσ3 , ~p = pσ3Πε, ρ~ = ρσ3ΠAr (Ar) , 其中 :εΠkB = 119. 8 K , σ = 0. 340 5 nm , 氩的相对
原子质量 Ar (Ar) = 39. 948 ,时间采用无量纲单位
第 5 卷第 2 期 2006 年 6 月
热科学与技术 Journal of Thermal Science and Technology
文章编号 : 167128097 (2006) 0220101205
L2J 流体自扩散系数及其与温度 关系的分子动力学模拟
Vol . 5 No. 2 Jun. 2006
Fig11 Variation of velocity autocorrelation function with time
L2J 势能曲线如图 2 所示 ,可见 ,粒子间距离 r
© 1994-2010 China Academic Journal Electronic Publishing House. All rights reserved.
(velocity autocorrelation function ,VACF) , 都 有 很 重要的研究价值[6] 。
1. 2 模拟体系 模拟对象为三维空间内的 500 个氩原子 , 采
用正则系统 (N ,V ,T 不变) ,三个方向均采用周期 性边界条件 ,利用 Verlet 算法求解粒子运动方程 , 采用邻近表法计算粒子间作用力 , 取截断半径 2. 5σ,邻近表半径 0. 3σ,每 5 个步长更新一次邻近 表 。初始构型采用 FCC 结构 ,时间步长为 2 ×10- 15 s ,模拟总步数 110 000 ,前 10 000 步采用温度 Π速 度调节使系统达到平衡 ,后 100 000 步用于统计计 算流体各种性质 。
其中
:
′
rjx
(
t)
为粒子真实位移
, rjx
(
t)
为 MD 程序中
粒子位移 ,nint ( x) 为最接近于 x 的整数 , r′j (0) =
rj (0) ,Lx 为系统方向上的尺寸 ,Δt 为时间步长 。
1. 1. 2 Green2Kubo 关系
∫ ∑ D
=
1 3 Nm
∞ Nm
〈 vj ( t)
2. 1. 2 松弛时间和截断半径
约化的速度自适应函数为 :
Nm
Nm
∑ ∑ vj ( t) ·vj (0) Π [ vj (0) ]2
j=1
j=1
图 1 为密度 1 374 kgΠm3 ,温度 90 K 时的速度