传热学 微尺度 ppt.ppt

微米 纳米尺度传热学 第四章微米 纳米尺度传热学中的基本分析方法 第四章微米 纳米尺度传热学中的基本分析方法 西安电子科技大学 一 导言 二 Boltzmann输运理论 三 分子动力学理论 四 计算流体流动问题的直接Monte Carlo模拟方法 五 量子分子动力学方法 一 导言 微米 纳米尺度传热问题本身的微观特点使得传统分析方法受到极大挑战 此时建立在宏观经验上的唯象模型不再十分有效 虽然在某些问题上 对二些传统流体力学 传热学理论及其相应的基本方程和界面条件作适度修正后 也可达到分析某些微系统传热问题的目的 但这种应用的范围受到很大的限制 要认识微米 纳米尺度范围内的传热规律 需要从微观的能量输运本质着手 以便揭示材料微结构中的动量和能量输运机制 按照从连续介质现象到量子现象的特征尺寸 迄今比较适合于分析微传热和流动问题的主要方法有如下几类 Boltzmann方程方法 分子动力学方法 直接Monte Carlo模拟方法及量子分子动力学方法等 其中Boltzmannn方法被公认为是一种 一 导言 极具普适性和有效性的工具 而分子动力学方法则用于揭示那些量子力学效应不明显时的物理现象的分子特征 它们也对分子统计理论 如Boltzmann方法及直接Monte Carlo模拟法 提供分子碰撞动力学方面的知识 直接Monte Carlo模拟则是一种计算微尺度器件内 通常其Knudsen数较大 尤其是稀薄气体流的流动和传热问题的方法 对于具有量子效应的物理过程 如光与物质的相互作用 金属材料中的热传导问题等 应采用量子分子动力学方法 并通过同时求解分子动力学方程及量子力学方程如Schrodinger方程来加以分析 本章内容将简要介绍这几类方法的要旨 它们是开展微米 纳米尺度传热学研究的重要理论基础 二 Boltzmann输运理论 1 概述 众所周知 在动力学理论中 空间和时间内的局域热平衡是一个隐含的固有假没 设体积的特征长度为lr 时间尺度为 r 则当物体的尺寸L lr 或真实时间t r 也或二者兼有时 则动力学理论不再成立 这是因为此时局域平衡假设不再有效 为此需要一个更基本的理沦 Boltzmann输运方程正是这样一种理论 它被认为是现有方法中用来分析微尺度能量输运现象的最具有普遍适用性 最基本和强有力的工具 虽然其最初的主要目的是用作气体研究 但发展至今已被推广用于范围极其广泛的各类介质 二 Boltzmann输运理论 该理论所具备的普适性令人惊讶 这是因为几乎所有宏观输运方程 如Fourier定律 Ohm定律 Fick定律及双曲型热传导方程等 均可由该方程导出 而且一些输运方程 如辐射输运方程及质量 动量及能量守恒方程等 也均可从Boltzmann方程导出 且对于流体 固体 多相系统等均具有良好的适应性 建立和求解Boltzmann方程的主要动机源于两类应用 其一是为了在当材料内能量载子的平均自由程远小于问题的特征尺寸时 能够从微观模型导出介质的宏观行为 所以这些应用是统计力学基本问题的一种特殊情形 而统计力学的任务就是在物 二 Boltzmann输运理论 质原子结构及其宏观连续介质行为之间建立一个桥梁 此方面的典型应用是解释气体的宏观行为 并从分子对相互作用定理计算出黏度及热传导系数 Boltzmann方程的第二类应用是在平均自由程与特征几何尺寸相比不再能忽略时对宏观介质的描述 很明显 在这样的条件下 人们不再能指望介质的 宏观行为 能够轻易地用密度 比热容 热导率等量来描述 虽然这些概念仍然是有意义的 且最后结果仍要借助于物体的可测量如温度来衡量 所以 在这些条件下 Boltanann方程作为一个可涵盖整个传热行为的方程而占有特别重要的地位 2 Boltzmann方程的简单推导 二 Boltzmann输运理论 这里以气体介质为例来加以说明 不过 如下推导对于流体和固体介质也是适用的 在推导过程中 分子之间的碰撞假设仅占其生命周期的非常小的一部分 这意味着只有双分子碰撞是重要的 考虑气体中每一分子受外力ma m为分子质量 a为分子加速度 作用 其大小可以是位置r和时间t但非速度v的函数 在时间t和t dt山之间 不与其它分子发生碰撞的分子的速度v将变为v adt 且其位置矢量r变为了r vdt 则在时刻 二 Boltzmann输运理论 t时落人体积单元r dr及速度范围在v dv内的分子数为f v r t dvdr 其中f为分子的分布函数 经过时间间隙dt后 若分子碰撞的效应可以忽略 则同样的分子而非其他分子的集合将占据体积r vdt dr 且速度在v adt dv 范围 这时的分子集合数为f v adt r vdt t dt dvdr 后期集合中的分子数一般与前期集合不同 因为碰撞分子会使初始集合中的分子过程偏离 它也会导致其他分子偏转从而使之成为末态集合分子 所以 后期集合中的分子净增量与dvdrdt成正比 并可由表示 于是有 二 Boltzmann输运理论 两边同除以dvdrdt并令dt趋于零 则可导出关于f的Boltzmann方程 即 或 其中Df表示式 4 1 左边项 按矢量符号表示为 二 Boltzmann输运理论 其中vx vy vz为速度矢量v的分量 而ax ay az为加速度矢量口的分量 对于混合气体 广义的Bolzmann方程可写作 其中s代表气体的种类 msas为在r t处作用在分子ms上的力 表示速度分布函数fs由于碰撞而发生的改变率 二 Boltzmann输运理论 该方程可修正后适用于更一般的分子模型 在具有球形对称的转动分子情况f下只依赖于v r t 及角速度 且f的方程与式 4 1 在形式上相同 对于更一般的模型 f将包含进一步的变量 以刻画分子的方位及其他性质 与这些变量相关的项必须一般性地出现在Bolmnann方程中 Bolmnann方程的一般形式通常可简洁地表述为 二 Boltzmann输运理论 其中f r p t 为随时间t 位置矢量r及动量p变化的系综粒子的统计分布函数 F为作用在粒子上的力 上式左边项称作漂移项 右边项则称为散射项 Boltzmann方程适用于服从某种统计分布的所有系综粒子 如电子 离子 声子 气体分子等 比如 对于受电场作用的电子 上式左边第三项可写作 其中E为电场矢量 e为电子电荷 k为电子波矢 二 Boltzmann输运理论 式 4 2 的右边项为由碰撞或散射引起的分布函数改变率 其严格表达式十分复杂 因为碰撞会使粒子从一坐标 r p 转变到另一坐标系 r p 于是可写出 其中W p p 为从状态p 到p的改变率 求和中的首项是由p 态到p态引起的 第二项则相反 散射率W通常为p的函数 式 4 6 可转化为 一个积分表达式 这使得Boltznmann方程是一个同时含有积分和微分项的方程 求解起来十分困难 因此常常要对Boltznmann方程进行适当简化 以实现一定程度的 理论分析 其中一种十分有效的途径是以下介绍的碰撞间隙 或松弛时间近似 理论 二 Boltzmann输运理论 3 碰撞间隙理论 不同粒子的散射及碰撞机制通常十分复杂 最常用的简化措施是通过引入松弛时间 来近似实现 即将Boltzmann方程中的散射项写作 二 Boltzmann输运理论 其中 0为平衡态下的分布函数 r p 为松弛时间 它是位置及动量的函数 在金属材料中 输运参数 如电导率 的温度依赖特性来自 随温度的改变量 在半导体中 口随温度的改变则主要由载荷子数目的改变引起 上述方法是一种基于碰撞间隙 而非平均路径的近似处理 该理论的有效性令人赞赏 其基本假设是在时间dt内 在给定的微小体积内有dt 个分子受到碰撞 并将其分布函数 改变为 0 这样的近似使得Boltzmann方程得以线性化 从而大大简 二 Boltzmann输运理论 化了方程的求解 它意味着 若系统偏离平衡态即 0不为零时 则碰撞将使该动力学过程按指数衰减 0 exp t 的形式恢复到平衡态 于是 无场效应 a 0 时的Boltzmann方程即为 其精确解可求得为 这实际上也可对时刻t t 进入容积并在时间到达 二 Boltzmann输运理论 体积r dr的v dv范围内的分子数进行归类而直接写出 对上式中 按t 的幂次展开并忽略t 2及更高项后 可得到另一简化式 碰撞间隙理论是一种非常简化的理论 其适用性是有限的 据Qiu和Tien关于金属介质中电子和声子相互作用的严格分析表明 碰撞项实际上由两部分组成 即 二 Boltzmann输运理论 头一项当电子与声子 晶格 的温度相同时为零 一般而言 它表明电子和声子之间的能量交换是一个非弹性散射过程 第二项当电子处于热平衡时变为零 所以 时间松弛假设很大时则不成立 并且松弛时间取决于晶格温度而非电子温度 实际上仅对部分散射过程有效 当电子和晶格温度差别 4 Boitzmann分布的场效应及碰撞效应 二 Boltzmann输运理论 Boltzmann方程最初用于处理非平衡态经典气体性质 气体状态由其分布函数 x y z px py pz 定义 可由六维空间x y z px py pz决定 因此 dxdydzdpxdpydpz为点x y z处在容积dxdydz 动量分量px py pz的dpxdpydpz范围内的粒子数 所以若要对气体进行完整描述需要求解 的改变受不同的外场影响 通常容易作到的是对 场 如电场或磁场或更一般的温度梯度等等 及碰撞的影响进行区分 它可写作 二 Boltzmann输运理论 稳态下 df dt 0 所以 在一给定时间间隙 由于外场引起的 改变必须由碰撞导致的 改变量来平衡 1 均匀温度场问题 温度梯度为零 处理电子时 人们最常采用的是波矢及其分量而非动量 对于空间均匀条件下 如没有温度梯度 的固体 则f关于x y z的依赖关系可忽略 它只取决于kx ky kz 于是就只需考察k空间 从而可极大地简化讨论 在处理电阻及热电性时 只要是均匀场及均匀样品 则可采用这一简化 围绕点kx ky kz处体积元dkxdkydkz的状态数由下式给定 二 Boltzmann输运理论 众所周知 在温度T下处于平衡的电子气f的表达式由著名的Fermi Dirac分布函数表示 此时其为f0 即 其中 为电子能量 f为电子平衡时的能量 kB为Boltzmann常数 这里应该指出 平衡分布f0的形式并不固定 如对于气体分子其满足Maxwell Boltzmann分布规律 对于电子则满足上述的Fermi Dirac分布 而对于光子和声子则满足Bose Einstein分布 二 Boltzmann输运理论 2 场的影响 若气体在x方向受均匀电场Ex的作用 则f的改变量很容易确定 在时间 t内 所有k空间内的占据态可均匀地由相同量代替 而kx由运动定理可知随时间均匀增加 即 于是 二 Boltzmann输运理论 这样 新的分布函数f与f0形式相同 但改变了 kx 即 由于实际应用中 我们关心的只是k空间内关于平衡分布的非常小的位移 于是可以写出新的分布函数f 二 Boltzmann输运理论 而且 f0只通过能量方程 4 12 而依赖于kx 所以有 若认为电子是一种粒子 则必须将其考虑为由我们所感兴趣的频率w和波数k附近的一定频率及波数组成的波包 于是电子速度v为这些波的群速度 其分量可定义为 由于电子能量 所以 二 Boltzmann输运理论 于是式 4 16 变为 最后 有 二 Boltzmann输运理论 3 稳态分布 利用稳态下的Boltzmann方程 并将式 4 18 及 4 19 代入 则有 于是 该式描述了受稳态电场及任意散射过程 由松弛时间 刻画 共同影响的电子数分布 它可用于计算各种输运参数 二 Boltzmann输运理论 5 Boltzmann输运理论导出的传热和流动守恒方程 由于Boltzmann的普适性 它可用来导出微尺度传热分析中所关心的几乎所有守恒及本构方程 如下给出其中的一些推导过程 用碰撞间隙理论写出的Boltzmann方程的一般表达式为 为研究粒子的能量输运 需要求解Boltzmann方程以获得分布函数f r p t 于是单位面积的能量流率或能流可写作 二 Boltzmann输运理论 其中q r t 为能流矢量 v r t 为速度矢量 p 是作为动量函数的粒子能量 注意 f r p t 的单位是单位体积单位动量内的个数 动量空间内的求和可转化为一个积分 该积分在引入状态密度D 后也可写成能量的积分 于是能流矢量可写作 二 Boltzmann输运理论 Fourier定律虽然求解Boltzmann方程并非易事 但可进行数种简化 若假设t r r为时间范围 则最通常的简化是不讨论式 4 7 中的时间变化量 此外 若假设L l lr 其中L为所考察尺度 lr为长度范围 l为平均自由程 则梯度项可近似为 f f0 沿x方向的一维Boltzmarm方程可求出为 该式称为准平衡假设 其中惟一包含非平衡因素的项是散射项 局域热力学平衡实际上是隐含在 二 Boltzmann输运理论 在近似df dx df0 dx中的 不过 由于局域平衡f0只能在长度范围lr内定义 该近似最后将变为df dx f0 lr 这一近似式及时间尺度内的近似在动力学理论中也得到采用 所以我们可期待得到类似的结果 由于平衡分布是温度的函数 于是有 由此可导出能流 二 Boltzmann输运理论 包含f0的第一项因在所有方向上的积分为零而可以消去 式 4 28 即为Fourier导热定律 其积分部分即为热导率 若假设松弛时间及速度均独立于粒子能量 则积分变为 这恰恰是动力学理论导出的结果 Cvl 3 在质量传递方面 采用类似推导 也可得出Fick扩散定律 二 Boltzmann输运理论 双曲热传导方程对Boltzmann方程两边同乘一个因子vx D d 并关于能量积分 则方程转换为 该式中的加速度项已被消去 考虑到这样的情形即L l lr及t r可作如下假设 i 松弛时间独立于粒子能量 ii 对项采用准平衡假设 则式 4 30 变为 二 Boltzmann输运理论 此即Cattaneo方程 将其与如下能量守恒方程结合 即可导出双曲型热传导方程 其中C为介质热容 式 4 33 的解是一种波的形式 它表明温度场的传播是以波的形式进行的 二 Boltzmann输运理论 值得注意的是在推导双曲型热传导方程时所作的假设 这里 感兴趣的时间尺度在松弛时间量级 而长度则远大于局域热力学平衡下的特征尺寸 在推导双曲型热传导方程及Fourier定律时的惟一差别是前者保留了瞬态项 这使得双曲型热传导方程在时间上是非局域的而在空间上则不是 Boltmmnn方程当然更具有一般性 因而可用于空间和时间内的非局域性和非均匀性的研究 二 Boltzmann输运理论 流体黏度在小量 情况下 方程隐含认为气体状态随时间变化不快 即f f0足够小 因而式 4 9 左边f可由f0代替 则有 直角坐标下 对于密度和温度均匀 流动沿ox轴且质量速度v0仅为z的函数的气体 式 4 34 变换为 二 Boltzmann输运理论 穿过z 常数平面的沿x方向的黏性应力满足 代人的表达式 4 35 则由于是的奇函数 式 4 36 简化为 其中 二 Boltzmann输运理论 该式大括号内的积分由于被积项是vx v0的奇函数而消去 第二个积分式可表示为p 于是 黏度 可写出为 质量 动量及能量守恒方程体力学 传热学及电子输运中所遇到的守恒方程均可通过Boltzmann方程导出 考虑函数为粒子动量的幂 即 其中正数n 0 1 2 其平均值可写作 二 Boltzmann输运理论 其中 为粒子数密度 将Boltzmann方程乘以 p 关于动量积分 则得到动量方程的一般形式 注意 每一粒子动量可分作两部分 即p pd pr其中pd为外场梯度作用下粒子集合运动的平均或漂移动量 pr为由热运动引起的随机动量分量 它代表耗散项 由于在动量空间所有随机动量分量的平均值为零 于是 p pd 则在零阶动量 n 0及 p 为常数 情况下 可得到连续或数守恒方程 即 二 Boltzmann输运理论 其中vd为漂移速度 pd m So是粒子的源或产生率 Si为粒子沉或移走率 动量守恒方程可在一阶动量即 p p mv情况下获得 即 第二项为张量的平均 由于对pr的奇次幂的平均为零 可得 其中 ij为单位张量 二 Boltzmann输运理论 左边第三项即为流体力学中所指的体积力 也许更合适的一种说法是将其看作势梯度项 这是因为热力学力DJ写作任意势梯度的函数 即F U 而势U是诸如重力势G 电化学势 等的总和 式 4 42 的右边为散射项 在松弛时间近似下 右边项可假设满足 其中 m为动量松弛时间 于是 动量守恒方程变为 二 Boltzmann输运理论 左边第三项包含随机粒子的运动动能形式 它代表粒子的压力 于是式 4 44a 可写为如下形式 左边第二项通常指水平对流项 当它可忽略时 式 4 4b 在零加速度下简化为 二 Boltzmann输运理论 在流体输运情况下 有 这与多孔介质中流体流动的情况等价 显然 式 4 42 除包含碰撞项外具有与Navier Stokes方程类似的形式 Navier Stokes方程可利用Boltzmann方程通过Chapman Enskog近似导出 且由式 4 42 右边可推导出耗散项 二 Boltzmann输运理论 由于能量 p2 2m 若取二阶动量 p p2 则由式 4 40 可导出能量守恒方程 即 其中 为能量密度 单位为J m3 J 为能流矢量 单位为W m2 它的一般式可表示为 二 Boltzmann输运理论 vd 能量的水平对流项 它由漂移项引起 q为耗散引起的热流 它由粒子的随机运动引起 于是可导出如下的简化能量方程 其中U为前面讨论过的所有势的总和 式 4 47 中的散射项分作能量源及沉 这里对其作一简短讨论 右边第一项是外力对粒子所作的功 因而必须出现在能量守恒方程中 为获得关于q的关系 需要用到更高阶动量的Boltzmann方程 不过 人们通常假设Fourier定律q T来使问题封闭 二 Boltzmann输运理论 但应记住这样的事实 即Fourier定律实际上是在空间和时间准平衡假设的基础上导出的 它并不总是有效的 考虑时间非局域而空间准平衡时 导出的更高阶关系是由式 4 31 所描述的Cattaneo热流方程 粒子系统的能流密度 来自熵运动及漂移的贡献 写作 二 Boltzmann输运理论 注意到因子3 2仅对具有三个运动自由度的单原子气体及分子有效 每一自由度具有能量kBT 2 对动量守恒式 4 44 乘以vd真并将其从能量守恒方程 4 49 减去 则可得到热能守恒方程为 注意到作功项 vd U由于功只增加机械能而不增加熵或温度因而可以消去 只有当该功通过散射而耗散掉 则系统熵及温度增加 式 4 51 中的散射项可写作如下形式 二 Boltzmann输运理论 其中T0为库温度 为能量松弛时间 右边第一项只是对应于平衡温度T0的能量松弛项 第二项是由于动量和能量松弛过程之间的差别引起的 由于碰撞会改变粒子动量而非能量 因而能量松弛时间与动量松弛时间存在差别 即使这两个时间相同 该项也不会为零 它来自粒子动能对温升的贡献 所以这部分是耗散掉的功 它导致熵增及温升 注意 虽然作功项 v F不出现 式 4 52 中的项代表了将热能加到系统中的耗散功 若考虑外电 二 Boltzmann输运理论 场对电子所作的功 电子 声子的相互作用最终会将该功耗散掉 从而导致能量损失到声子 所以库温度是声子的温度 比值称为流体Prandtl数Pr 这是因为流体扩散率与其对应的松弛时间成反比 光子和声子辐射输运方程光子和声子不像电子及分子那样满足数的守恒 不过 它们遵循能量守恒规律 光子和声子强度可定义如下 二 Boltzmann输运理论 其中Ik为波矢k下的强度 v为波矢k下的速度 s为极坐标 为能量 强度也可用频率w及矢量k方向的极坐标角度 定义 即 若对Boltzmann方程 式 4 2 乘以因子可得如下方程 二 Boltzmann输运理论 其中 这里 每一个W为一个散射率 显然 散射项非常复杂 需要解释一下 二 Boltzmann输运理论 式 4 54 是在频率w 极坐标s及方向 处基于强度的能量守恒方程 现在考虑式 4 55 中的第一个求和 由于发生在频率w w 极坐标s s及方向 下的散射 Iw r w s t 的强度得以增加 第二项则由于朝其他频率 极坐标及方 向散射而引起强度的Iw r w s t 损失 注意到 若考察的是光子 则该项表示光子 光子之间的散射 但它并不代表光子与其他粒子之间的散射 所以 该项表明的是粒子散射的类型 即或者是光子 光子 或者是声子 声子之间的散射 这是人们所熟知的内散射项 在光子辐射输运中 散射常常只考虑沿 方向而不沿频率及极坐标方向发生 二 Boltzmann输运理论 这是因为非弹性的光子散射在工程计算中通常是忽略掉的 除非在Raman散射中才涉及 然而 对于声子散射 非弹性散射如正常及倒易过程则很常见 且必须在该项中予以考虑 此外 这样的声子 声子散射常常是在不同的声子极化坐标LO 光学 纵向 TO 横向光学 LA 纵向声学 及TA 横向声学 进行 它满足碰撞过程中的声子能量及动量守恒 式 4 55 中的第二项加和由于与其他粒子的相互作用而导致强度Iw r w s t 增加或降低 粒子类型由j标记 由动量及方向定义的相空间标记 如金属或半导体导带中的高能电子 二 Boltzmann输运理论 可通过电子 声子之间的相互作用在一定极坐标下发射声子 如LO声子 来降低能量 这里电子由j标记 而声子由 标记 该声子的频率 方向和极坐标是由散射过程的能量及动量守恒决定的 在光子辐射输运中 该项常常称为黑体源项 这对黑体辐射而言是正确的 不过 在半导体激光或光发射二极管这样的器件中 光子并不按黑体光谱发射 而是在由半导体电子带状结构决定的某一黑体光谱带内发射 所以 该项在式 4 55 中是作为一般辐射项存在的 类似地 当声子或光子被另一粒子吸收并从系统中移走时将出现损失项 二 Boltzmann输运理论 很清楚 在由式 4 55 所描述的最一般形式中 散射项增加了求解的困难 所以 人们常常采用松弛时间近似以便于求解 这样 辐射输运方程即简化为 右边第一项是朝外散射的项 s作为其散射松弛时间 第二项为光子 声子吸收 或朝其他粒子如电子 或光子到声子 或声子到声子等的能量传输 其中 s为吸收时间 二 Boltzmann输运理论 第三项为发射项 1 e为其发射率 这里 来自其他粒子的能量贡献给强度Iw 最后一项为来自其他频率及立体角 的内散射项 辐射输运方程能以一种更简单的方式写出 其中内散射项可全部忽略 但假设w w散射以使平衡恢复 即达到 这是声子辐射输运中经常所作的假设 那时内在频率散射并不使声子平衡恢复 二 Boltzmann输运理论 在光广散射中 对于光学上薄的或弹性光子极限 可以得到两个黑体表面间的辐射热流为 对于声子情形 它通常称为Casimir极限 在光学厚的或扩散极限 热流变为 其中为光子热导率 类似的结果也可从声子辐射输运问题中得到 三 分子动力学理论 实现分子模拟的技术步骤十分直观 它采用一组具有指定粒子对作用规律的模型分子对结构空间进行采样 采样可以是随机的 通过Monte Carlo方法 也可以是确定性的 通过分子动力学经典运动方程 本节介绍分子动力学所涉及的重要问题 分子动力学计算是按照分子系统的时间演化进行的 由此可产生相互作用分子的详细轨道图景 它比之Monte Carlo方法的优点在于凝聚相的动力学性质 如自扩散率及黏度 可在平衡或远离平衡 受大的外场作用时 的情况下求出 1 概述 三 分子动力学理论 由于计算机软硬件技术的巨大发展 微尺度模拟方法 如分子动力学方法 Monte Carlo方法 正在成为分析微尺度物理问题的通用工具 其在估计微结构材料的输运性方面起到十分重要的作用 这些经典方法正被广泛地看作除理论和实验之外的第三种认识微尺度现象的重要方法 当电子的动力学行为变得显著时 则尚需求解复杂的瞬态量子力学方程 此即所称为的量子动力学方法 三 分子动力学理论 在分子动力学模拟中 分子的动力学行为通常假设遵循经典运动方程 这在典型分于的山Broglie波长远小于分子平均间距 即式 4 58 的平移运动中确实是一种令人满意的近似 其中h为普朗克常数 m为粒子质量 V为体积 N为颗粒数 三 分子动力学理论 在此限制下 存在一个容易达到的能态的基本连续分布 转动问题也可按经典方式处理 只要保证转动能隙远小于kBT且分子处于其振动基态 这些准则对大多数分子是正确的 但似乎不包括含有H或D的小分子 如HCl及H20 尽管这些小分子也可按分子动力学方法作例行处理且也获得了很大成功 分子动力学方法在时间和空间范围内考察一系列相互作用的被模拟分子 其最常见的形式是对牛顿运动方程a d2r dt2 F m 其中r为分子矢量位置 F为分子上的净作用力 m为分子质量 进行数值积分 微观规范系综为分子动力学方法的本质系综 三 分子动力学理论 假设系统中的电子处于基态 且精确地在分子动力学尺度内跟随核运动 于是只需考察核的轨迹 在某些情况下 更方便的作法是采用其他经典动力学表达式 如发生动力学系综 作用在每一分子上的力由相互作用势 r 其中 r r 为分子i和j的作用中心间距 即r ri rj 给定 在不同的作用势下运动方程可采用不同积分格式 只要确知所考察系统的相互作用势 则分子动力学方法不受问题的相态及热力学条件限制 它可用于处理各种相变问题 虽然计算量很大 分子动力学方法无需对诸如几何对称 输运性或热力学行为作先验假设 三 分子动力学理论 所有计算均基于分子相互作用势在第一原理基础上进行 热力学及输运性质是作为计算结果而非假设条件 分子动力学方法已被用于化学 生物学 物理学及材料科学等领域 目前 它被用于模拟超临界环境下一个微滴的蒸发问题 在热科学领域内 分子动力学方法正被尝试用于从数值上预测固体薄膜的热导率以及单晶硅和硅纳米线的热导率 以及验证一些经典传热学规律的准确性 分子动力学方法也被用于模拟kennard Jones型流体的法向激波的内部结构 三 分子动力学理论 正确选择相互作用势对于分子动力学方法的准确模拟至关重要 现今已找到一些材料的相应作用势 如Stillinger Weber作用势已经被若干作者证实是对硅材料的一种很好的描述而硅材料构筑成了现代电子技术的基石 可以预见 该作用势对微尺度器件中的一些热问题分析将具有特别重要的意义 所以这里特别予以列出 但由于Lennard Jones势是迄今应用最为广泛的作用势 所以在后文讨论中仍以其为主 顺便指出 由于目前能反应真实的物质问相互作用势极为罕见 微米 纳米尺度传热学将来可能有的进展也许也就源于此 这是因为关于各类情况下液 液 液 气 液 固及固 气等情况下的作用势方程的构建 测量及其应用将提供无数的研究机会 三 分子动力学理论 Stillinger Weber作用势包括二体及三体项 其形式分别为 其中rab为原子a b之间的距离 r a时 而r a时h 不存在 cab为矢量rab和rac之间的夹角 三 分子动力学理论 在式 4 59 4 61 中 r以距离 2 0951 为一单位 能量的单位以 3 4723 10 12erg 1erg 10 7J 表示 Stillinger与Weber又寸式 4 59 4 61 中的参数进行研究后 得到A 7 04956 B 0 602225 p 4 q 0 1 a 1 8 21 1 2 对于理想四面体角度cos 1 3以及零温度理想钻石晶格 则不存在三体能量项 分子动力学模拟的一般推进格式按照预报一修正方式进行 可归纳如下 a 利用位置 速度 加速度等的当前值预测其在下一个时刻t t的相应值 b 在新位置处评估力及加速度ai fi mi c 采用新的加速度对所预测的位置 速度 加速度等进行修正 2 硬球分子动力学 三 分子动力学理论 硬球的作用势在硬球直径r 处呈尖锐的不连续性 超出这一范围 则无相互作用 而对于r 范围 则作用势无穷大 该模型分子的特殊性使得其分子动力学计算步骤与渐进微作用势的情况非常不同 硬球分子动力学数值上利用气体有限容积 分子 的动力学理论 它遵循一系列双体碰撞 硬球分子动力学也可能出现三体碰撞情形 即三个球同时碰撞 对于精度有限的计算机而言 这在理论上是可能实现的 三 分子动力学理论 d 在回到步骤 a 之前 计算所需变量 如能量 维里参数及序参数等 一个比较理想的分子动力学模拟格式应具有如下特征 a 应尽可能快速且采用尽可能大的时间间隙 并且所要求的计算机内存应最小 b 应允许使用较长的时间步 c 应尽可能精确地重复经典轨道 d 应满足已知的能量和动量守恒定律 且时间上可逆 e 在形式上尽可能简单和易于操作 分子动力学方法最早被用于理想硬球分子流体 如下对这一类系统作一讨沦 由于二体作用势在模拟中最简单也最常用 所以本书仅就这种情况加以介绍 关于三体相互作用势及其他问题将不作为讨论的内容 三 分子动力学理论 硬球分子动力学考虑系统中N N 1 2的每一对硬球 给定一对硬球 以指标i和j表示 在任意瞬间的坐标分布 计算其下一碰撞时间tij 于是对每一对硬球产生一系列碰撞时间 从而第一个碰撞对可由此确立 然后将所有硬球按第一对硬球碰撞所需时间以常速度在空间内移动 之后可假设碰撞对按经典理想弹性球力学理论经历一次碰撞 这涉及到碰撞后硬球速度的计算 下一个碰撞对可通过测试碰撞对系列来找到 考虑前两个硬球的碰撞伙伴 按时间顺序对每一次碰撞重复多次 三 分子动力学理论 硬球分子动力学的问题之一是需要一个单独的步骤来处理上一结构的密度增大问题 与固定容积内模拟单元内硬球直径增加等价的密度逐渐增长会导致硬球的交叠 这是与硬球势的内在动力学格式不相符合的 如下描述的运动方程将基于一个更近期的基本硬球碰撞方程 它考虑了密度的逐渐增长问题 其中定义了一个压缩因子d dt 并作为输人参数写到计算机程序中 模拟单元体积假设是固定的 但粒子直径 认为随时间以速率 增加 或降低 系统中每一硬球直径在每一时刻是相同的 三 分子动力学理论 给定一系列直径为 坐标为r 速度为v的硬球 则对两硬球i j的表面在经过时间tij后接触时 有 其中 0是碰撞对移动之前的 值 在碰撞间隙即 自由飞行 区间 其半径按同样速率变化 从而当其发生碰撞时 直径将变得相同 将式 4 62 展开并求解二次方程式 可得到碰撞时间tij 三 分子动力学理论 其中 对于负数或复数tij情形 则两个硬球不在碰撞路径上 碰撞时 可通过变换 其中 tij 来影响碰撞 3 连续势分子动力学 三 分子动力学理论 在连续势分子动力学中 经典动力学方程采用有限差分格式来积分 其中分子位置和速度按相等时间间隔 t顺序计算 这是与硬球分子动力学或Monte Carlo方法非常不同的地方 因为二者分别是一次处理一个分子或者按不同的时间间隙进行 在连续作用分子动力学中 需要知道分子速度和加速度 以便在时间山内将所有分子同时移到一系列新位置 这里介绍一种常用的处理运动方程的格式 该格式由Verlet引入分子动力学中 它利用相邻两个时刻分子的坐标来产生下一个时刻的分子位置 考虑三个时刻t t t及t t的分子i ri 其中 t用作模拟的时间步长 则有 三 分子动力学理论 该式利用了ri t 的Taylor展开 O t3 表示正比于 t3的高阶级数 类似地 按逆向时间 即时间间隔为 t 展开时 有 将式 4 64 与式 4 65 相加 有 2020 3 15 75 可编辑 三 分子动力学理论 所以式 4 66 比之式 4 64 更准确 式 4 8 在大的时间步长下能得到稳定格式 这也是它之所以获得广泛采用的原因 对式 4 66 按如下格式重新组合即成为位置Vedet格式 人们感兴趣的量是分子在时间t t的位置ri t t 它是该时刻惟一的未知量 在时刻t 方程右边所有量均为已知 式 4 67 于实际应用很不方便 因其在任意时刻 需要记录时刻t和t t的坐标 而且周期边界输运问题 盒子边长要加到分子坐标中或从中减去 可能发生在时间t和t t之间 则t t时的坐标将不得 三 分子动力学理论 不在 S范围内作调整 这使得位置Verlet格式 式 4 67 的使用非常不方便 一个更方便的不受到周期边界输运影响的形式是采用如下定义的所谓半时间步长速度 即 它称为半时间步长速度 是因为若粒子以一个定常的或平均的速度在时间t t t内运动 则其速度可由式 4 68 定义 所以 我们有ri t t ri t vi t t 2 t对反方向时间步长类似有ri t t ri t vi t t 2 t 于是按半时间步长速度 式 4 67 可重新写出 三 分子动力学理论 或 其中 若将式 4 65 从式 4 64 中减去 可得 三 分子动力学理论 利用式 4 72 可定义所谓的Verlet 时间步 速度vi t 注意 时间步速度的精度比起半时间步长速度低一个山量级 所以对实际应用中需要知道时间步速度的情况要尽量避免采用 不过 需要绝对知道同样时刻 即t n t n为一整数 下的位置和速度的情形 至少在乎衡态 并不多 因此 半时间步长速度在计算中常被采用 式 4 70 及 4 71 一起即为人们所熟知的速度Verlet蛙跳格式 它是一种稳定性和适应性很好的格式 采用这些方程 则分子路径可在一系列连续时间步长下得到 从而由此可利用统计力学表达式来求出物理及热力学性质 一个更精确的时间步速度格式可以表示为 三 分子动力学理论 蛙跳格式事实上为Hamilton方程的时间中心格式 及 一个一直有争议的问题是 速度Verlet和位置Verlet格式会导致不同的相空间轨迹 速度Verlet 三 分子动力学理论 格式在用于某种自由度高度热化的系统的快速及慢速 或 热 及 冷 动力学时非常直观 它可通过引人多时间步长格式来实现 其中自由度更新频率在大时间步长 t下稍微慢些 而热自由度以更小的时间步长 t传播 其中 t n t 其中n 1 分子动力学计算格式还有许多变形 Satoh就各种分子动力学格式的稳定性所作的对比研究表明 最优越的格式应是速度Verlet格式 之后才依次是蛙跳格式 Beeman格式等 渠跳格式比起这些格式均要差一些 因其更复杂且稳定性和能量波动性较差 三 分子动力学理论 加速度可依据牛顿第三定律方程由作用力求出 其中叫mi为分子i的质量 Fi为作用在其上的净力 我们有 其中rxi为笛卡尔直角坐标ri的x方向的分量 ix为x方向的单位矢量 更具体一些 我们考虑加速度和力的x方向分量 三 分子动力学理论 其中rij ri rj rxij是rij在x方向的分量 由于kennard Jones势应用最为广泛 作为一般性的讨论 如下将主要以其为例 LJ表示为 其中 为势作用最大深度 它常称为 阱深 为势穿过r轴即 时的r值 三 分子动力学理论 参数 通常称为分子直径 粒子对作用力 可通过作用势的解析形式得到 即 若定义 则 通常 分子动力学程序90 以上的计算时间都花在评估作用力或其他成对的附加性质上 三 分子动力学理论 在LJ程序中 通常为方便起见 将长度按 表示 能量按 表示 而m表示为粒子质量 所有其他量均可用这些基本单位表示 如压力可表示为 3 而时间为 m 1 2 事实上 LJ势只是如下一系列可能的所谓n m势 其中n为式中排斥项的指数 m为吸引力部分的指数 中的一种形式 对于LJ势 要附加一个限制 并且 若其在r rmin处最小 则 三 分子动力学理论 由这些限制可导出 例如 对n 36及m 18的情形 可有A 4及 这与12 6形式的LJ势清情况下所得到的 不同 模拟时 粒子初始坐标经常选在晶格上 这样可以避免引起过度的分子交叠 从而导致灾难性的巨大作用力 若物质密度太低而不能保持晶格特性时 它会导致出现融化现象 从而形成流体 这对于研究低密度或高密度流体的情况十分有用 高密度 三 分子动力学理论 流体可通过以小步长 如对几千个时间步长增加0 05 逐渐增加初始低密度流体 0 4 的密度来有效地产生 以这种方式 系统能更快地在所希望的高密度流体态达到平衡 并远离初始给定的晶格坐标 图4 1示意了这一步骤 它比之等待晶格在此密度下融化来得到所需高密度的情况更为优越 典型小系统所包含的周期边界条件倾向于阻碍晶格融化 三 分子动力学理论 图4 1利用分子动力学生成高密度流体的融化步骤 图形示出厂初始不稳定的分子 在几千个时间步长融化而形成稀释流体态 晶格 之后分子坐标及元胞边长以均匀尺度持续下降 在每一密度处使之平衡 重复该过程直至达到所需密度 压力在初始态下为较大负压 之后随时间逐渐增加 三 分子动力学理论 初始速度由能产生Msxwell Boltnlmm 速度分布的高斯型分布随机数发生器选择 为避免分子动力学盒穿越空间 这样会导致伪动力学性质 在模拟初始时盒的净动量设为零 之后程序在计算机误差范围内尽可能使之保持零动量 比如将粒子对的作用力同时加到Fxi中而又从Fxj几中减去 因此 只要考虑粒子盒的总动量 任何力的误差均会抵消 该步骤可列写如下 初始时 即t 0 速度为vi t t 2 随机数Ri 即 三 分子动力学理论 所模拟元胞的总动量为 于是对每一分子i 从其初始设定的质量速度中心减去上述动量的N 1 即有 该式对元胞的全部动量均作了零处理 从而质量中心动量不随时间变化 在LJ程序中 所有的粒子质量均相同 且通常以简化单位表示 分子动力学胞中的分子动能必须与所设定的初始温度协调 三 分子动力学理论 模拟中的粒子对相互作用通常在距离r rc 处截断 这里rc L 2且为L元胞最小边长 超出这一距离的所谓 可忽略 相互 作用对于热力学性质的计算非常重要 这在作用范围更大的情况F尤其如此 由于流体具有平均几何对称性 其分子轨迹对这些截断距离以外的相互作用的敏感性较少一些 对相互作用能 有 其中对能量的长程修正为 三 分子动力学理论 该式假定g r rc 于是 在此近似下 ULRC可写为一个易于计算的表达式 类似地 对压力的长程修正为 对LJ势 我们有 三 分子动力学理论 对于压力的长程修正特别大 所以修正必须包括在任何常压模拟中 否则 所设定的系统密度可能会是错误的 三 分子动力学理论 分子动力学方法最重要的应用之一在于物质输运系数的计算 分子水平上的动力学行为在宏观上体现为实验可测的输运性质 如白扩散系数及剪切应力等 分子动力学模拟已被用于理解分子水平上的松弛过程 以及分子结构及化学特性对这些可测性质的影响 其目标之一是对统计性质的模拟 且已获得较大成功 输运性质刻画流体在外加宏观梯度作用下的动力学响应行为 扩散描述的是分子通过热驱动从系统中的一部分移动到另一部分的过程 数学上 4 物质输运系数及某些热物理量的计算 三 分子动力学理论 物质的扩散趋势由白扩散系数D刻画 其单位为 长度 2 时间 1 考虑这样一个系统 即体积V中含有N个分子 且分子i速度为vi 则系统的流定义为 各向同性流体中穿过单位面积的分子输运量与其平面法向的浓度梯度成正比 即 此即正规定义D的Fick第一定律 三 分子动力学理论 系统中的热传导也是固体和流体微观动力学行为的一种体现 系统内热非平衡性的描述可由热导率k定量化 为测定k 常常需要建立一个关于温度梯度的初始态 于是Fourier定律被用于定义热导率 另一个与流体状态强烈相关的输运系数是刻画流动阻力的黏度 牛顿最早将其定义为联系平板单位面积上及其平行方向的速度梯度的比例系数 剪切应力 为单位面积A上的力F 它阻碍与静止平板相距h的速度为v的平板的运动 则有 三 分子动力学理论 其中 s为牛顿黏度 分子动力学模拟最显著的特点之一是其允许 至少在原理上 对一个N体系统输运性质进行计算 其中控制方程采用Green Kubo线性响应公式 如下给出计算白扩散系数 黏度 剪切黏度及体黏度 及热导率等的Green Kubo公式 流体能量和动量的传输可通过三种机制进行 第一种机制涉及分子的实际运动 它在低密度下起决定作用 且可通过动力学理论相当准确地处理 第二种机制在高密度情况下起决定作用 它 三 分子动力学理论 依赖于相互作用力的影响来实现 分子之间互相施加作用力 从而导致动量沿空间传输 第三种机制涉及到这些动力学及相互作用项之间的关联 事实上 所有的输运系数均具有纯动力学及相互作用的成分 幸运的是 控制输运性质的分子过程通常发生在皮秒时间尺度 这正好落人分子动力学可处理的日引司尺度内 任意输运系数 可写作如下两种等价公式中的一种 三 分子动力学理论 其中a为常数 为积分变量 式 4 97 中的第一个表达式称为Einstein表述 第二个表达式称作Green Kubo表述 在后一种表述中 输运系数写为一个自修正函数的积分面积 两种形式在统计上的等价性可通过改变积分变量来证明 自扩散系数对于自扩散系数D 有 D a 1 3及 ri 则其可由均方位移得到 三 分子动力学理论 实用过程中 采用 得到的结果更好 由相应Green Kubo公式得到的等价自扩散系数写为 其中i表示模拟中的任意粒子 表示在一时间区域t0内的平均 即 三 分子动力学理论 于是 自扩散系数既可通过粒子的均方位移得到 通过考察每一分子在不同时刻的位置 也可通过依赖于时间的分子平移速度的修正函数得到 Green Kubo公式的本质在于输运性质是通过一个平衡态的涨落决定的 它不涉及非零的宏观输运 剪切黏度为求得剪切黏度 s 回到式 4 97 可有 s a V KBT及 于是 三 分子动力学理论 其中对于单原子流体 有 以上选择应变张量xy分量是随意的 但也可采用xz及yz方向的分量进行分析 事实上 实用中 三个分量常常被同时采用以提高统计精度 研究表明 在热力学极限内 运动方程是否恒热并不重要 它们均能给出相同的自扩散系数及剪切黏度 牛顿剪切黏度在实验上可通过几种毛细流动装置或晶体共振法测得 剪切黏度常常写作一个频率无穷小的弹性模量G 及特征松弛时间的乘积 它为模化剪切 三 分子动力学理论 应力修正函数下的面积 即 s G s真 G 是一个统计量 体黏度对于体或压缩 膨胀黏度 B 有 B a V KBT及 则有 三 分子动力学理论 其中 虽然已有很多方法用于剪切黏度的测量 但至今尚无直接测定体黏度的方法 体黏度通常由声波吸收部分传统表述中的修正项得到 热导率对于热导率 有 及 其中 三 分子动力学理论 及 由于粒子在离开元胞及从元胞对面返回时的位移等于盒子的边长 在周期边界条件中不能采用均方位移方法 于是 在采用连续势模拟时 也许 三 分子动力学理论 除了自扩散系数外 通常均要采用Green Kubo表述 因为可按上述形式消去而不受由于周期边界条件引起的位移制约 因为它们只涉及到速度或者能被写成粒子对的分离项 硬球流体的黏度 s真是用于计算时间修正函数的重要时间尺度 它是声波穿过边长为L的模拟元胞的时间 这会成为系统误差的一个来源 因而在处理量级为L cs 其中cs为声速 上的时间修正函