导热系数的分子动力学模拟研究及相关问题的探讨.pdf

中国工程热物理学会 传热传质学学术会议 2 0 0 0 3 1 1 4 导热系数的分子动力学模拟研究及相关问题的探讨 冯晓利t 戒 润华大学工理力攀系北京1 李志信过场元 仪 盯 8 a T a W i L x 1 加心 抽勿m m a i l R i n 沙u e e d u s u 摘 要 对于 采用分子动 力学 方法 研究 导热系效的 背最 研究现状及存在的问趁进行了 林述和分析 总 1 1 1攀 亡 lr牙产 结了 通过分子动 力学 模拟 方法求祥 导热 系数的物 理模型和墓本算法 究中引入分子动力学棋拟方法铭要考察的影响因素和几个生要问翅 关 词 徽尺度传热 导热系数 分子动力学 讨论了目 前在徽尺 度导热问 趁的 研 1 在热科学中引入分子动力学方法 气卜 JI 近来 具有嫩小特征空间 尺度和时间尺度的传热问 题倍受研究者的关注 在徽尺度 情况下 经典传热理论的适用性需要拷问 也需要尝试新的 研究方法 以 直接探求微尺 度传热的特殊规律 徽电子 光电子技术和徽电子机械系统 M E M S 的迅速进步 还在 不断地推进器件尺寸的徽小化 在空间尺度为纳米jk级和时间尺度为飞秒 1 0 1 S l级的 条件下 都有与传热相关的重要课题 许多问题甚至播要从分子I 原子的层面来进行研究 P 1 分 子动 力 学 M o l e c u l a r D y n a m ic s 记为M D 方 法 就 是 将系 统中 的 分 子 原 子 视为 经典 粒子 从统计力学出发研究其性质的计算机棋拟方淤l A M D方法的基本原理 是 通过求解系统中受特定分 f I 原子间作用势和系统外加约束 作用的分子I 原子的 N e w t o n运动方程 模拟系统随时间演进的徽观过程并统计计算系统 的 平衡今数和 钧运性质 由 此 可以 从系统的徽观细节求得工程中关心的宏观参数 M D 摸拟的迅速发展和应用得益于计算机运算能力的提高 目 前 采用徽机和普通工作站计 算的 粒子数目 可达 1 0 0 0 1 0 0 0 0 1级 而最近报道的 在超级计算机上的 井行 模拟实验 中的 原 子数目 为1 0 0 0 0 0 0 13 1 在计算 物理 化学 材料和 化工 等领域 M D方法已 经得到 广泛的 应用 是比 较成熟的 研究手段 在热科学中 采用 M D方法的研究正逐步开展起来 在徽尺度系统的物运过程 相变与界面现象 热脉冲加热 薄膜沉积和徽加工等问题的 研究中 MD摸拟都是有效和有前景的方法 导热系数是控制徽尺度系统尤其是固体材料导热的重耍参数 在微尺度实验费用高 昂 甚或难以 进行的 情况下 M D棋拟则可能是预报导热系数的 有重要价值的 方法 本文将 集中讨论采用M D棋拟方法研究导热系数的一 些相关问 题 2 平衡分子动力学 E MD 与非平衡分子动力学 N E MD 方法 E M D方法意味着通过模拟平衡系统求得导热系数 K u b 的线性响应理论 l i n e a r r e s p o n s e t h e o r y 为 将导热系 数 表 述为 可直接模拟的 信息 提供了 必要的 理论 框架 在线 性响 应范围内 G r e e n K u b o 关系式给出导热系数可由 系统徽观热流的时间相关得到1 2 1 一 1T 2 J J Q t J 23k1 T 2 e 0 咖 1 1Y L erv 告 E T f oi j 告 m v 告 u yj 这里k a 为B o lt z m a n n 叹 d日月 粉 u 为 双体作用 势能 以 下 标i 和J 标识 不同的 分 子 E M D模拟通常 都是 在徽正则 或 正则系综中进行 并且上述热流表达式也仅适用于双体作用势能 早期的工作中所用的 势能模型主要是硬球势能模型 后来广泛使用的是较为 接近真实系统的 L e n n a r d J o n e s双 本课 为国家自 招科学签金 N a 5 9 fl6 0 1 3 资 助项目 丁一 一一 气一 Y i 潇 有 甲 体势能摸型 目 前 有针对性的 更为复杂的势能模型也 逐渐应用于 对各种材料的模拟 现实中更多需要考察导热系数的情况是系统中 存在温度梯度 井有不可逆热流发生 在这种惰况下 E M D方法不再适用 须采用N E M D方法模拟不可逆非平衡过程 N E M D 的核心算法是对系统施加扰动 建立非平衡的抽运过程 根据系统在空间上的均匀性和 边 界 条件 N E M D方法 可 分为 各向同 性的 h o m o g e n e o u s 和 各向 异性的 in h o m o g e n e o u s 所谓各向同性的 N E M D 是指非平衡的系统在空间上具有某种均匀性 通常是指与周期 性边界条件相容 E v a n s 发展了 一种与周期性边界条件相容的各向 同 性 N E M D算法 在表征平衡系统 的正则 H a m i l t o n i a n描述中加入依赖于时间的扰动项 同时保证系统的 动t恒定 由 此荃 于推广的线性响应理论 可在外加大扰动下考察系统的 非平衡晌应 计算出系统的导热 系数 这种方案的特点在于全面修正运动方程 而系统中并未实现非平衡的热流或沮度 分布 另一类各向同性N E MD方案中 具有周期性边界条件的元系统中有非平衡热流发生 在计算获得沮度分布和热流后 即 可由F o u ri e r 定律计算系统的导热系数 在这种方案中 既可以给定系统中冷端与热端的温度 再计算热流 也可以 给定热流后 再计算沮度梯 度 由于计算热流更难收敛 目 前大多数算法中一般都先设定热流值 再统计沮度分布 I k e a h o j i 和H a f a r j o ld ill 以 及M l l l le r P la t h e b 1等先 后独立 发 展了 这种各向 同 性 E r方 案 图1 为M u ll e r P l a t h e 的 各向 同 性N E M D 算法的 示意图 系 统在三个方向 均 采用 周期 性边界 将模拟单元 元系统 在x 方向 分为N 层 其中标记为0 的层为冷端 标记为N 1 2 的 层为热端 从N 1 2 层到N 1 层与前一半成镜像 对称 在模拟中 每隔w个时间步长 选取0 层中 速率最大的原子和N 1 2 层中 速率最小的原 子 交换其速度矢t 由 于原子速率分布的 温 度范围总是大于这两个层的沮度之差 所以这 样的一对原子总能找到 即冷层中速率最大的 原子的速率大于热层中速率最小的原子的速 率 系统达到稳态的时候 该速率交换过程产 翻叻二0t几 v at 八 c 2F 411枯 1 圈1一 种典型的H o m o O W G N E b n 方案 生的热流将与0 N 1 2 层间的反向的导热热流相平衡 于是 又二一 m 22 乙 万价 一 V c 少 J 川曰晒 2 tL 州a T l a x 2 由 于原子间交换的是其速度矢且 所以该算法保证了 总动量 总动能和总能t守恒 交换 原子速率的同时 势能未受影响 该算法中 热流可以 根据原子交换造成的动能转移直 接确定 并且由w表征热 流的 大小 而 温度分布的计算也较为容易 H o o v e r 等171 最早提出 各向异性N E M D方案 来模拟边界非周期性的 在空间上各向异性的 系统中的稳态导热 如图 z所示 系统在x方 向 位于高温与低温恒温瑞之间 在Y和 z 方向 采用周期性边界条件 在 x方向分层并计算各 层中的局域温度 该方案的核心在于实现系统 边界与恒温墙相接触的定温算法 在与恒握墙 了 一 卜 圈2 典 型的加 加m 叩 加 N 曰 M D方 案 3 5 声 相接的 壁面层中 为粒子的运动方程引 入附加项 3 十 11 1二1卜 仁卜沪 卜1 ltl 井满足 4 其中 二 令 一bm v m v 一 2 军 M V 一 2 k8T 一 朴 吞 f vt b T v 5 C i e c o t t i 等I i 发展了 各向 异性N E M D算法的变种 物理摸型与图2 相同 但在摸拟系 统与恒沮绪的接触时采用随机方案 系统中粒子的 运动方程为 JZ 一N m 粉 Y V U 栩 f r 二J 州 6 其中f为随 机作用力 强制较子反射回系统中 仅当 靠近系统边界的 较子试图穿越恒沮绪时 施加该随机作用力 反射粒子的 速率由 恒温绮的沮度确定 在各向 异性 N E MD算法 中 值沮瑞间 通过系统的热流与沮度分布都可以 直接获得 从而可以由 F o u r i e r 定律求得 导热系致 由于系统中存在沮 度梯 度 各向 异性 N E MD得到的是对应于系统平均沮度的 导热系数 K o t a k e 等m 所发展的各向 异性N E M D算法中 系统在热流方向采用镜面边界条 件 3 几个盆要问 瓜的讨论 在徽尺度导热问 翅的M D模拟研究中 须探讨并解决一些荃本问 题 1 所有上述的E M 旧算 法 和N E W 算法 都以与系统的某种统计力学系综描述相容 为目 标 也都以 计算宏观系统的导热系数为目 标 在 E M D和均匀性 N E MD算法的发展 中 都认为与周期性边界条件相容意味着可以模拟宏观系统的性质 井以 此作为评价算 法的 准则 之一 然而周期性的边界条 件惫味 着系统有限 无界 周期性系统的计算结果是 一 1 育 否代表宏观性质不能一概而论 2 经典的 E和 N E M D算法都针对祠密流体的导热发展起来 徽尺度传A中 更 为多见是固体材料的导热 在固体中经典的 B o l t z n r a n n统计通常不再适用 因此播要研 究固体中晶格导热的t子描述 并由 此对M 刀模拟结果进行修正 3 即 使对单原 子分 子组 成的 模型 系统 不同的 研究 者给出 的 热流计算式仍有差 别 须从理论上解决热流计算问题并在M D棋拟中 详尽 考察对热流有贡献的几项的相对大小 4 上述算 法都以 单原子系 统的 双体势能模型为 荃础 对于更为复杂和常见的物 质 需要根据其可兼的 势能模型重新研究热流计算公式及导热系数的计算式 参 考 文 献 C h o u E 氏 L a k e s J R L i s o g X G e t a l A n n ua l R e v ie w o f H e a t J 3 m 叮 e r 1 0 1 9 9 9 1 4 1 A l l e n MY T i l d e s le y D I C o m p u te r S im u l a t i o n 矿钩a id s O x f o r d C l e r e d o n P r e s s 1 9 8 7 h 印刃 访 w w c s m a a l g o v G C 吮 四 胡e h IP I C S I n n ml E v a a s D 工P 切 L e f t A 9 1 1 9 8 2 4 5 7 压 咧1 毛 出自 匆 匕B M o l P 切 8 1 1 9 9 4 2 5 1 M Me r F l a t h e F J C L P h y s 1 9 6 1 9 9 7 6 0 8 2 H o o v e r W G A m m