Vol.28 No.12 刘育松等：纳米尺度孔隙内气体导热系数的分子动力学摸拟 。1183。 子数为608个.在分子动力学模拟的过程中，壁 0.1MPa条件下进行的，采用Vedet算法的速度 面内的氮气分子均保持固定不动并且和孔隙内 形式来求解分子的牛顿运动方程，L山作用势的截 的氮气分子相互之间具有作用势.由于所研究的 断半径取为10o.模拟首先在NVT系统中进行， 对象是孔隙内的氮气分子，所以在X,Y和Z三 即模拟系统中的分子数(N)、体积(V)和温度 个方向上采用周期性边界条件， (T)要保持不变用以确定各部分能量及热流. 采用Green-Kubo关系式来计算孔隙内氮气 各个分子初速度设定为零，然后通过一调温系数 的导热系数，即： 来自动调节温度，时间步长△t由△tNe/m/6= 1 k=3WBTJ。<J0J》d (2) 00001计算.当NVT系统达到平衡，即在一定 的模拟时间内，温度T保持恒定，则启动NVE系 其中，V是模拟系统的体积T为系统的温度，符 统下的平衡MD模拟.在NVE系统下，模拟系统 号()表示的是在整个模拟时间内系统的平均值， 中的分子数(N)、体积(V)和能量(E)都保持不 J(t)是热流 变.图2为NVT系统达到平衡后，各个分子在不 =空+片是，5》 (3) 同的分速度和速率区间分子数的分布.'m为分 子的平均速率.从图中可以看出，通过平衡MD 其中，：是分子i的速度矢量，F是分子i由于 方法模拟得到的分子分速度和速率的分布与统计 L山作用势的作用受到邻近分子j的作用力.e 力学得到的Maxwell速度和速率分布曲线基本一 可以通过下式给定： 致.从模拟结果计算得出，此时分子的平均速率 (4) 为460ms1.通过对分子在不同时间的位置和 其中，m是氮气分子质量. 速度的模拟计算可知，分子的平均自由程受立方 体的边长严格限制，远小于自由空间分子的平均 2 模拟结果与讨论 自由程. 图3为NVE系统的能量变化曲线.可以看 平衡MD模拟计算是在温度300K和压力 35[( 30 30l (b) ■模拟结果 25 ■模拟结果 -Maxwe分布 -Maxwell分布 20 ■ 15 15 ◆◆ ■ 10 10 ◆ 5 ■ 要■ -0.6 -0.2020.61.0 86 020.40.60.81.01.21.4 (瑞) 图2氨气分子的分速度和速率的分布.()氨气分子沿x方向速度分布：(b)氮气分子的速率分布 Fig.2 Distributions of the velocity component and speed of nitrogen molecues:(a velocity component distribution in x direction:(b) speed distribution 出，在经过1.4ps后，系统各部分能量达到平衡. Zeng等9依据理想气体模型和动力学理论 J(t)以及时间相关函数的变化曲线分别如图4 推导出了多孔介质中气相导热系数： 和5所示.在此基础上，采用式(2)计算得出孔隙 kg= 内氮气的导热系数为6.8×103W·(mK)1. Q225r-1.25)0.461mg (P/knT)(8kBT/m)2 MD模拟结果仅为同样条件下自由空间的氮气导 0.255.pr (P/kg T d 热系数(2.4X102W(mK)-1)的1/3左右，其 (5) 原因主要有：一是由于孔隙壁的存在，使得分子平 其中，S,是固相的比表面积Ppr为密度，中为气 均自由程减小限制了分子输运能量的能力：二是 孔率，P为压力，mg为分子质量，dg为分子直径 由于分子之间以及分子与孔隙壁分子之间的相互 =cp/cm,cp和cr分别为气体比定压热容和比 作用势，使分子在孔隙壁上的漫反射得到加强. 定容热容.为了比较MD模拟结果、式(5)的计算子数为 608 个 .在分子动力学模拟的过程中, 壁 面内的氮气分子均保持固定不动, 并且和孔隙内 的氮气分子相互之间具有作用势 .由于所研究的 对象是孔隙内的氮气分子, 所以在 X , Y 和Z 三 个方向上采用周期性边界条件 . 采用Green-Kubo 关系式来计算孔隙内氮气 的导热系数[ 8] ,即 : k = 1 3 Vk B T 2∫ ∞ 0 <J(0)J(t)>dt (2) 其中 , V 是模拟系统的体积, T 为系统的温度 ,符 号〈〉表示的是在整个模拟时间内系统的平均值, J(t)是热流 . J(t)= ∑i viεi +1 2 i ,∑ j , i ≠j rij(Fi jvi) (3) 其中, vi 是分子 i 的速度矢量, Fij 是分子 i 由于 LJ 作用势的作用受到邻近分子 j 的作用力 .εi 可以通过下式给定: εi = 1 2 m vi 2 + 1 2 ∑i ≠j u(rij) (4) 其中 , m 是氮气分子质量. 2 模拟结果与讨论 平衡 MD 模拟计算是在温度 300 K 和压力 0.1 MPa 条件下进行的 , 采用 Verlet 算法的速度 形式来求解分子的牛顿运动方程, LJ 作用势的截 断半径取为 10σ.模拟首先在 NV T 系统中进行, 即模拟系统中的分子数(N )、体积(V)和温度 (T)要保持不变, 用以确定各部分能量及热流. 各个分子初速度设定为零 ,然后通过一调温系数 来自动调节温度,时间步长 Δt 由 Δt ε/ m/ σ= 0.000 1计算.当 NV T 系统达到平衡 , 即在一定 的模拟时间内 ,温度 T 保持恒定 ,则启动 NVE 系 统下的平衡M D 模拟 .在 NVE 系统下 ,模拟系统 中的分子数(N)、体积(V)和能量(E)都保持不 变.图 2 为 NVT 系统达到平衡后 ,各个分子在不 同的分速度和速率区间分子数的分布 .v m 为分 子的平均速率.从图中可以看出 , 通过平衡 M D 方法模拟得到的分子分速度和速率的分布与统计 力学得到的 Maxwell 速度和速率分布曲线基本一 致.从模拟结果计算得出, 此时分子的平均速率 为460 m·s -1 .通过对分子在不同时间的位置和 速度的模拟计算可知 , 分子的平均自由程受立方 体的边长严格限制 , 远小于自由空间分子的平均 自由程. 图3为NVE系统的能量变化曲线 .可以看 图 2 氮气分子的分速度和速率的分布.(a)氮气分子沿 x 方向速度分布;(b)氮气分子的速率分布 Fig.2 Distributions of the velocity component and speed of nitrogen molecul es:(a)velocity component distribution in x direction;(b) speed distribution 出,在经过 1.4 ps 后 ,系统各部分能量达到平衡. J (t)以及时间相关函数的变化曲线分别如图 4 和5 所示.在此基础上, 采用式(2)计算得出孔隙 内氮气的导热系数为 6.8 ×10 -3 W·(m·K)-1 . M D 模拟结果仅为同样条件下自由空间的氮气导 热系数(2.4 ×10 -2 W·(m·K)-1)的 1/3 左右 ,其 原因主要有:一是由于孔隙壁的存在,使得分子平 均自由程减小, 限制了分子输运能量的能力;二是 由于分子之间以及分子与孔隙壁分子之间的相互 作用势,使分子在孔隙壁上的漫反射得到加强 . Zeng 等[ 9] 依据理想气体模型和动力学理论 推导出了多孔介质中气相导热系数 : k g = (2.25 r -1.25)0.461 mg(P/ k B T)(8 k B T/πm g) 1/2 0.25S s ρpo r -1 + 2(P/ k B T)πd 2 g c V (5) 其中 , S s 是固相的比表面积, ρpor为密度, 为气 孔率 , P 为压力, mg 为分子质量, dg 为分子直径, r =cp/ cV , cp 和cV 分别为气体比定压热容和比 定容热容 .为了比较 MD 模拟结果 、式(5)的计算 Vol.28 No.12 刘育松等:纳米尺度孔隙内气体导热系数的分子动力学模拟 · 1183 ·