自由分子区内纳米颗粒的热泳力计算*

崔杰 苏俊杰 王军† 夏国栋 李志刚

1) (北京工业大学能源与动力工程学院, 传热强化与过程节能教育部重点实验室, 传热与能源利用北京市重点实验室, 北京 100124)

2) (香港科技大学机械及航空航天工程学系, 香港)

基于非平衡态分子动力学模拟方法, 研究了自由分子区内纳米颗粒的热泳特性.理论研究表明, 纳米颗粒与周围气体分子之间的非刚体碰撞效应会明显地改变其热泳特性, 经典的Waldmann 热泳理论并不适用,但尚未有定量的直接验证.模拟计算结果表明: 对于纳米颗粒而言, 当气−固相互作用势能较弱或气体温度较高时, 气体分子与纳米颗粒之间的非刚体碰撞效应可以忽略, Waldmann 热泳理论与分子动力学模拟结果吻合较好; 当气−固相互作用势能较强或气体温度较低时, 非刚体碰撞效应较为明显, Waldmann 热泳理论与模拟结果存在较大误差.基于分子动力学模拟结果, 对纳米颗粒的等效粒径进行了修正, 并考虑了气体分子与纳米颗粒之间的非刚体碰撞效应, 理论计算结果与分子动力学模拟结果吻合较好.

1 引 言

微细颗粒在温度不均匀的气体中运动时, 会受到周围气体分子的差异性碰撞, 从而受到一个沿温度梯度相反方向上力的作用, 即为热泳力.颗粒在热泳力作用下的运动, 称为热泳[1−4].通常情况下,来自高温区域的气体分子与颗粒之间的碰撞更为剧烈, 这便导致颗粒受到的热泳力方向与温度梯度方向相反, 指向冷端.在某些情况下, 也可能出现热泳力与温度梯度方向相同的情况, 即反向热泳[2].对热泳现象的研究在薄膜制备、微纳米制造、环境科学、气溶胶等领域中都有重要的意义[5−10].近年来, 针对颗粒输运现象的理论与实验研究越来越多, 但是, 相关现象的内在机理还有待深入研究和分析[11−18].

颗粒的输运特性与其Knudsen 数Kn (Kn =λ/ L )密切相关.其中, λ 为气体分子的平均自由程, L 表示颗粒的特征长度.对于球形颗粒来说,特征长度为颗粒半径R, 并使用符号KnR表示.根据Knudsen 数的大小, 可以将颗粒的动力学行为粗略地划分为连续介质区(continuum regime)、过渡区(transition regime)和自由分子区(free mole−cule regime).在连续介质区(KnR≪ 1), 流体分子的自由程很小, 分子间相互作用的频率很高, 此时, 可以把流体看作连续介质处理, 使用经典的Navier−Stokes 方程来对颗粒的输运特性进行分析计算.Epstein[19]基于此方法研究了连续介质区的热泳现象, 但求解过程中采用的边界条件并不合适.Brock[20]采用完备的滑移边界条件得到了连续介质区热泳力的一阶近似解.其他学者们也尝试开展了高阶近似解的计算[21].然而, 此方法也受到了一定的质疑[22−24].在自由分子区(KnR≫ 1), 由于气体分子的平均自由程远大于颗粒的特征长度, 被颗粒反射后的分子需飞行很长的距离才会与其他分子相碰, 因此可以忽略入射气体分子与被颗粒反射气体分子间的相互作用.此时, 可假设气体分子的速度分布函数不会受到悬浮颗粒的影响.1956 年,Waldmann[25]基于气体分子与颗粒之间的刚体碰撞模型, 得到了自由分子区中颗粒在单原子气体中所受热泳力的计算公式:

式中, mg为气体分子的质量, kB为Boltzmann 常数, T 为温度, κ 为气体热导率, R 为纳米颗粒半径.∇ T 为气体环境的温度梯度.越来越多的理论与实验研究表明, 随着KnR增大, 颗粒所受的热泳力值逐渐接近Waldmann 公式的计算结果[2], 因此, (1)式也被称为热泳力的自由分子极限.在连续介质区(KnR～ 1), 气体分子的离散特性刚开始体现, 气体分子的速度分布受颗粒运动的影响较大, 但又不够达到可以看作连续介质的程度, 此时,求解Boltzmann 输运方程是十分困难的.有学者曾尝试将自由分子区的理论扩展到过渡区, 但结果并不理想[26].Talbot 等[27]通过改变几个Brock 公式中的参数, 使Brock 的连续介质区理论公式在KnR≫ 1 时收敛于(1)式.

通常情况下, 气体分子的平均自由程为100 nm量级.对于粒径较小的纳米颗粒而言, 其特征尺寸远小于气体分子的平均自由程, 因此纳米颗粒在气体中的动力学行为一般处于自由分子区.Waldmann公式是基于气体分子与颗粒之间的刚体碰撞模型得到的, 但是, 随着颗粒尺寸减小至纳米尺度, 颗粒与气体分子之间的非刚体碰撞效应将成为影响颗粒与气体分子之间动量传递的重要影响因素之一[16].文献[17]中, Li 和Wang 基于非刚体碰撞模型, 得到了自由分子区内纳米颗粒所受热泳力的理论计算公式:

式中, mp为颗粒的质量; mr= mgmp/(mg+ mp),为气体分子与颗粒的约化质量; Ω(1,1)*与Ω(1,2)*为无量纲碰撞积分[11].对于刚体碰撞来说, Ω(1,1)*=Ω(1,2)*= 1, (2)式退化为(1)式.换言之, (1)式为(2)式在刚体碰撞假设下的特殊情形.但是, (2)式的准确性尚未在实验或者模拟中得到定量的验证.一方面受限于纳米尺度下测量技术的不足, 另一方面颗粒的热泳运动一般较弱, 很容易被其布朗运动所掩盖, 因此关于纳米颗粒热泳现象的实验研究往往比较困难[28].分子动力学模拟方法可以直观地研究并分析颗粒的输运及受力特性[28−33], 在模拟中, 可以在纳米颗粒上人为地施加简谐势来对纳米颗粒进行约束, 从而明显降低了颗粒布朗运动对其热泳力计算的影响.

对于实际纳米颗粒而言, 其形状多为非球形(例如碳纳米管等).非球形颗粒在流场中的取向会对其受力产生重要的影响, 颗粒的瞬时热泳力的大小和方向都与其形状和取向有关, 较为复杂[34].在自由分子区, 气体分子的平均自由程远大于颗粒粒径, 因此, 气体分子与颗粒的随机碰撞会使得颗粒高速旋转.颗粒的宏观受力体现为气体分子与颗粒之间微观动量传递量在长时间内的统计平均结果.如果颗粒没有受到较强外势场的作用, 颗粒取向分布较为均匀, 可以采用球形颗粒模型近似研究非球形颗粒的热泳特性.

本文基于非平衡态分子动力学模拟方法, 研究了自由分子区内球形纳米颗粒的热泳现象, 计算得到了纳米颗粒所受热泳力的大小和方向.计算结果表明: 当气−固相互作用势能较强或气体温度较低时, 纳米颗粒热泳力的计算需要考虑气体分子与纳米颗粒之间的非刚体碰撞效应; 进行颗粒粒径修正后, 基于Li−Wang 公式(2)的计算结果与分子动力模拟结果吻合较好.

2 分子动力学建模

2.1 分子动力学模拟系统

本文采用分子动力学模拟方法计算纳米颗粒在自由分子区内所受热泳力.分子动力学模拟系统如图1 所示, 直径和质量分别为DP和mP的纳米颗粒被浸入到气体分子数为N 的模拟区域中, 选择固体分子与气体分子质量相同, 均为m.在系统y 和z 方向采用周期性边界条件.在x 方向的边界上采用镜面反射边界条件, 并在邻近边界的两端分别建立温度控制区域(高低温热浴), 其对应的温度分别为Th和Tl.温度控制区域在x 方向上的尺寸为lT= 0.1lx, 此处, lx为x 方向的系统尺寸.

图1 分子动力学模型图Fig.1.Snapshot for the MD simulation system.

采用经典的Lennard−Jones(L−J)势函数来模拟气−气、气−固、固−固分子间的相互作用:

式中, 参数ε 和σ 分别为L−J 势参数, rij为原子i 和原子j 之间的距离.系统中的原子(固体原子、气体原子)参数均采用Ar 原子的参数[35]: σ =3.405 Å, ε/kB= 114 K.截断半径rc= 2.5σ.研究气−固结合强度对颗粒所受热泳的影响发现, 气−固结合强度εGP可以在一定范围内变化,

其中kij为可在一定范围内变化的参数.此外, 在固体颗粒内部相邻原子之间额外引入finite exten−sible nonlinear elastic(FENE)势函数以保持纳米颗粒为准球形, 其势能表达式为[30]

式中, 势能参数设定为R0= 1.5σ; kFENE= 30ε/σ2,为FENE 势能的弹性系数.kFENE的大小决定了纳米颗粒的导热特性, 文献[30]对热泳力随kFENE变化情况进行了研究分析, 结果表明纳米颗粒的热导率对其所受热泳力影响不大.球形纳米颗粒的直径使用下式进行计算[30]:

式中, xc为纳米颗粒的质心, NS为组成纳米颗粒的固体原子数量, 研究过程中选择纳米颗粒的固体原子数量至少为40, 纳米颗粒的质量mP= NSm.Galliero 和Volz[30]的研究表明, 当mP/m ＞ 10 时,颗粒的热泳特性将不再受到颗粒质量的影响.相比于气体分子, 纳米颗粒的质量与尺寸已足够大, 因此可忽略颗粒质量对热泳的影响[36,37].

由于纳米颗粒的布朗运动, 颗粒会在系统中随机行走, 并受到一个瞬时的曳力(一般情况下曳力的大小比热泳力高2—3 个数量级), 为热泳力的计算带来较大的困难.因此, 本文在颗粒质心与系统中心引入额外简谐回复势能Uhar:

式中, rc和ro分别为颗粒质心与系统中心坐标,khar为简谐回复势能的弹性系数.作用在纳米颗粒上的简谐回复力将均匀施加在纳米颗粒的每一个原子上.简谐回复势在系统三个方向上给颗粒所施加的回复力分别为F = (Fx, Fy, Fz).当系统达到稳定后, 可得颗粒所受到的热泳力为

在热泳力的计算中引入简谐势, 可以一定程度限制颗粒的随机运动, 消除颗粒运动过程中所受曳力对其热泳力的影响.文献[32]研究了热泳力随弹性系数khar的变化情况, 结果表明, 只要弹性系数khar大于0.1ε/σ2, 其大小对统计结果影响就可以忽略不计.本文取khar= 30ε/σ2.

本文中, 除特别说明, 所有参数将采用无量纲形式(用上标“*”表示), 采用下列各物理量进行无量纲化: 长度x0= σ, 温度T0= ε/kB, 时间t0=(mσ2/ε)0.5, 力F0= ε/σ.例如, 温度的无量纲形式为T*= T/T0.

根据气体分子运动论, 气体分子的平均自由程为

式中, n 为气体的分子数密度, d 为气体分子的有效硬球直径.根据参考文献[38], 温度T*= 2.63时氩原子气体分子的有效硬球直径约为0.94σ.为了验证模拟系统中气体的状态, 从模拟系统中选择一组气体的状态参数与气体的状态方程进行对比,

式中, P*为NVT 系综的压力, Bi(T*)是维里系数.取ΔT*=0.92.P*随n 的变化情况如图2 所示.维里系数的取值分别为= 0.3734,= 0.0421.

图2 压力P*随数密度n 变化图Fig.2.The pressure P* versus the gas molecular number density n.

分子动力学模拟中, 使用Velocity−Verlet 算法对颗粒的运动方程进行积分, 为了保证系统能量的稳定, 取时间步长Δt*= 10—3(Δt = 2 fs).初始时刻, 气体分子随机分布于模拟系统中, 初始速度为温度T*下的麦克斯韦分布.系统首先在NVT 系综下进行弛豫, 在高低温热源温度的算术平均温度下运行直到t*== 103, 达到稳定.然后, 采用速度校正法[39]将高低温热浴的温度分别控制为和, 直到t*== 1.5 × 104, 从而得到一个恒定的温度差ΔT*, 系统的采样时间将在t*=到时间段内进行(= 8 × 105).本文所给出的任何物理量都是在此时间范围进行的长时间统计平均得到的.为了提高结果的准确性, 在相同的热力学宏观状态下, 取不同初始速度分布, 至少8 次独立模拟结果的平均值作为最终结果.

3 有限空间效应

Waldmann 热泳力公式的一个前提为与颗粒相互作用的气体分子是不受空间尺寸束缚的(空间体积无穷大).在这种情况下, 气体热物性不会受到空间边界的影响.然而, 在模拟研究中, 气体分子的活动范围总是被具有一定空间尺寸的边界所束缚, 在自由分子区中, 气体分子的平均自由程可能会与空间尺寸的量级相当, 此时必须考虑有限的体积效应对气体热物性的影响[40].取气体环境所处的空间尺寸L 作为特征长度, 此时努森数可用KnL表示.在图1 所示的分子动力学模型中, L 为x 方向高低温热浴的间距, 即L = 0.8lx.1972 年,Phillips[41]借助Chapman−Enskog 分布函数对Bo−ltzmann 输运方程进行了求解, 得到了有限空间内颗粒的热泳力表达式:

式中, αmt和αmn分别为颗粒表面切向和法向的动量协调系数, α1和α2分别为高低温壁面的温度协调系数.在本文的分子动力学模拟系统中, 可以认为气体分子是完全协调的, 即= α2≈ 1 以及αmt= αmn≈ 1.

分 别 采 用 Phillips (11)式 和 Waldmann(1)式进行热泳力理论计算, 并与分子动力学模拟计算结果进行对比分析.Phillips 公式的计算结果表示为FT,Phillips, 分子动力学模拟结果为FT,MD.Waldmann 公式将分别使用气体的宏观热导率κ 和有效热导率κ'计算, 计算结果分别用FT,B和FT,E表示.图3(a)与图3(b)分别表示在kij= 5.26 (强气−固相互作用下)和kij= 1.0 (弱气−固相互作用下)情况下FT/FT,B随l 的变化图.由图3 可以看出, 模拟系统中“壁面”间距是影响颗粒所受热泳力的重要因素之一, 这说明有限空间效应确实存在.FT,E与FT,Phillips吻合较好, 表明Phillips 对Wald−mann 热泳力计算式的修正是可以通过对气体热导率的修正来实现的.

表1 分子动力学模拟系统的几何特征参数Table 1.Geometric and characteristic parameters of the simulation systems.

图3 FT, MD/FT, B, FT, E/FT, B 和FT, Phillips/FT, B 随模拟空间尺寸的变化 (a) kij = 5.26; (b) kij =1Fig.3.Influence of the size of simulation box on FT, MD/FT, B,FT, E/FT, B and FT, Phillips/FT, B: (a) kij =5.26, (b) kij =1.

对于强气−固相互作用(图3(a)), 分子动力学结果FT,MD与理论结果FT,E和FT,Phillips存在较大误差, 这种偏差是由于气体分子与纳米颗粒之间的非刚体碰撞所引起的.对于弱气−固相互作用(图3(b)), 气体分子与纳米颗粒之间的非刚体碰撞效应较弱, 此时, 基于刚体碰撞模型假设的Wald−mann 理论仍然适用于纳米颗粒, FT,E和FT,Phillips都与分子动力学结果FT,MD吻合较好.图4 为热泳力随气体有效热导率κ' 的变化曲线.与前文结果类似, 对于弱气固耦合的情况(kij= 1.0), FT,E和与分子动力学结果FT,MD吻合较好; 对于强气固耦合(kij= 5.26), 分子动力学模拟得到的热泳力则明显高于基于有效热导率的理论计算结果.

图4 热泳力FT,MD 和FT,E 随气体有效热导率κ' 的变化Fig.4.Influence of the effective thermal conductivity of the media gas on the thermophoretic forces.

4 结果与讨论

4.1 系统温度(气-固结合强度)影响

由前文所述可知, 将气体的有效热导率引入到Waldmann 公式中, 可得到有限空间中气体的Waldmann 自由分子极限.但是, 对于纳米颗粒来说, 特别是当气−固相互作用较强时, Waldmann 热泳力计算式的误差较大.基于分子动力学模拟结果,热泳力随系统平均温度气−固结合强度kij和温度梯度 ∇ T 的变化情况如图5 所示(图5(a)中的插图为在相同系统参数下本文分子动力学结果与文献[28]中结果的对比).由图5 可以看出, 对于纳米颗粒来说, Waldmann 热泳力计算式的适用性受系统平均温度与气−固结合强度kij的影响, 当系统温度较低或气−固结合强度较大时, Waldmann公式给出的热泳力计算结果误差较大; 随着系统温度的升高或气−固结合强度的减小, 纳米颗粒所受热泳力的分子动力学计算结果逐渐收敛于Wald−mann 公式的预测值.这是由气体分子与纳米颗粒间非刚体碰撞效应的强弱受温度或气−固结合强度影响所导致的.当结合强度较大或温度较低时, 势能在气体分子与纳米颗粒进行动量交换的运动轨迹中占主导地位; 而当结合强度较小或温度较高时, 气−固相互作用势能的影响十分微弱, 刚体碰撞模型近似成立.

图5 不同参数下热泳力的分子动力学结果FT,MD 与Waldmann 公式结果FT,E 比较图 (a) 环境温度T*; (b) 气−固结合强度kij; (c) 温度梯度∇T∗Fig.5.The variation of thermophoretic force FT between present MD simulation result FT, MD and Waldmann equa−tion result FT, E under different parameters: (a) The environ−ment temperature ; (b) the intensity of gas−particle in−teraction kij; (c) temperature gradient ∇ T∗.

4.2 颗粒尺寸影响

颗粒尺寸的大小不仅影响着颗粒在高低温两侧的动量交换, 还影响着颗粒与气体分子之间的碰撞模型, 是影响颗粒所受热泳力的重要因素之一.为了便于比较, 定义无量纲热泳力:

式中, FT为热泳力的理论计算结果, FT= FT,E为基于有效热导率的Waldmann 理论计算结果, FT=FT,LW为考虑非刚体碰撞效应之后(2)式的计算结果, 对应的无量纲热泳力分别表示为和图6 展示无量纲热泳力和随颗粒直径的变化图.模拟过程中系统的参数为:= 55,ΔT*= 0.92, n*=0.0054.颗粒粒径变化范围为 DP∗≈ 2.5 到≈9.98, 且模拟系统处于自由分子区.由图6(a)可以看出, 随着颗粒尺寸逐渐接近分子尺寸(= 1),无量纲热泳力单调增加, 而对于颗粒尺寸较大的纳米颗粒, 颗粒所受热泳力的分子动力学结果则逐渐收敛于Waldmann 公式的预测值(由于计算能力的限制, 并没有给出更大尺寸颗粒的计算结果).对于纳米颗粒而言, 由于纳米颗粒与气体分子之间的非刚体碰撞效应, 热泳力结果并不能收敛于经典的Waldmann 理论值.(2)式对Waldmann理论进行了修正, 修正后的无量纲热泳力随颗粒尺寸变化结果如图6(b)所示.相对于Wald−mann 理论的计算结果, 无量纲热泳力计算误差明星降低, 这验证了(2)式的准确性.但是, 不同气−固结合强度下的计算结果还存在一定的误差,这是由于当气−固结合强度较大时, 气体分子的动能不足以克服气−固结合势能的束缚, 气体分子会被吸附于纳米颗粒表面, 从而引起纳米颗粒尺寸的增加, 造成理论与模拟结果的偏差.

图6 无量纲热泳力 随颗粒直径 变化图 (a)刚体碰撞模型; (b) 非刚体碰撞模型FT = FT, LW1; (c) 非刚体碰撞模型FT = FT, LW2(考虑气体分子吸附引起的颗粒尺寸变化)Fig.6.The dimensionless thermophoretic force as a function of for different gas−particle interaction intensity:(a) The rigid body collision model FT = FT, E; (b) the non−rigid body collision model FT = FT, LW1; (c) FT = FT, LW2(wherein the vary of particle size is considered).

图7 对颗粒表面气−固相互作用势能以及气体分子动能分布进行了统计, 同时计算了颗粒表面气体分子密度分布函数g.图7 中, 黑色与蓝色虚线分别为颗粒表面第一层与第二层原子的位置, r'=.由图7 可以看出, 当气−固结合强度较小时(图7(a)), 势井深度较浅, 气体分子所携带的动能足以克服势能的束缚.此时, 纳米颗粒表面的气体分子数密度与环境中的气体分子数密度相当, 并没有发生吸附现象.当气−固结合强度较大时(图7(b)), 颗粒表面第一层内的势井深度较深,而气体分子的动能相对太低, 使得该层气体分子无法摆脱势能的束缚, 从而被吸附在纳米颗粒表面,这就造成了纳米颗粒表面的第一层内的气体分子数密度远远高于环境中的气体分子数密度.通过对纳米颗粒的粒径进行修正, 可以得到修正后的等效粒径, 将等效粒径代入 (2)式, 无量纲后的结果表示为如图6(c)为无量纲热泳力随颗粒尺寸的变化结果.结果表明, 计算结果与分子动力学模拟结果吻合较好.

图7 纳米颗粒表面气体原子的势能Ep, 动能Ek 及密度函数g 分布图Fig.7.The potential energy Ep, kinetic energy Ek and dens−ity function g on the surface of nanoparticles.

5 结 论

本文研究了自由分子区内纳米颗粒的热泳特性, 基于非平衡态分子动力学模拟方法计算了纳米颗粒所受热泳力.考虑有限空间效应对气体热物性的影响, 引入气体的有效热导率计算得到了较为准确的Waldmann 自由分子极限, 并与分子动力学模拟结果进行了对比.结果表明, 对于纳米颗粒来说, 当气−固相互作用势能较弱或气体温度较高时,Waldmann 热泳力公式与分子动力学模拟结果吻合较好; 当气−固相互作用势能较强或气体温度较低时, 由于气体分子与纳米颗粒之间的非刚体碰撞效应较为明显, 基于刚体碰撞模型假设的Wald−mann 计算式与模拟结果存在较大误差.基于Li−Wang 计算式, 并考虑由于纳米颗粒对气体分子的吸附效应引起颗粒等效尺寸的变化, 计算结果与分子动力学模拟结果较为吻合.