全文HTML
--> --> -->图 1 物理模型结构示意图
Figure1. The diagram of physical model structure.
Lennard-Jones (12-6) 势能模型应用于流体原子间以及固液原子间, 其表达式为
界面热阻与温度跳跃ΔT密切相关, 即满足关系ΔT = RK·Q, 其中Q为热通量. 界面处的温度跳跃为壁面温度与界面流体温度之差(Twall–Tfluid), 但由于界面处流体温度分布往往会出现很大的波动, 导致统计温度跳跃时会出现很大误差. 因此, 本文利用Navier型边界条件来研究界面热传导[12], 即利用连续热传导方程来预测纳米通道界面处流体的温度. 预测方式如图2所示, 通过拟合延长通道主流区的温度分布直至壁面处, 将预测得到的流体温度作为界面处流体温度.
图 2 固液界面处流体温度的预测
Figure2. Prediction of fluid temperature at the solid-liquid interface.
在分子动力学模拟中, 固液界面处的温度跳跃可表示为[12]
模拟系统中壁面温差会诱导通过固液界面的热通量Qx, 其计算公式如下[18]:
为了获得更加精确的模拟结果, 每个分子的性质应该接近于真实的分子, 即在初始时刻, 需要给每个分子设定随机速度以及位置. 确定系统的系综为正则系综(NVT), 采用Nose-Hoover恒温器使系统中的所有原子在100 K的温度下随机运动, 速度大小遵循高斯随机分布[19], 公式表示为
模拟计算时利用共轭梯度法使系统能量最小化, 即迅速使系统达到相对合理的状态, 减少模拟计算成本, 设置时间步长为2 fs, 计算100万步(2 ns)使系统驰豫稳定. 随后取消NVT系综, 利用微正则系综(NVE)对系统进行模拟. 为诱导系统产生一恒定的热通量, 需要给壁面设置一个恒定的温度, 分别作为热源和冷源, 如图1所示, 红色代表热壁面, 蓝色代表冷壁面. 采用Langevin恒温器将壁面温度恒定为140和90 K. 在不施加外力的条件下, 控制壁面温度后, 进行200万步(4 ns)的平衡计算, 系统平衡后, 再运行100万步(2 ns)将所需要的数据进行统计, 最后将这些数据进行时间平均. 此外, 还研究了驱动力对界面热阻的影响, 以Poiseuille流动为研究对象, 控制壁温后, 对每一个氩分子施加沿y正方向的驱动力, 运行200万步(4 ns)实现动态平衡, 系统平衡后, 运行计算100万步(2 ns)将所需要的数据进行统计, 最后将这些数据进行时间平均. 本文模拟采用LAMMPS程序编写[20].
模拟工况 | 壁面温度 | 固液势能强度εwf/ε | ||||
工况1 | 热壁面(140 K) | 6.0 | 1.5 | 1.0 | 0.75 | 0.5 |
冷壁面(90 K) | 6.0 | 6.0 | 6.0 | 6.0 | 6.0 | |
工况2 | 热壁面(140 K) | 6.0 | 6.0 | 6.0 | 6.0 | 6.0 |
冷壁面(90 K) | 6.0 | 1.5 | 1.0 | 0.75 | 0.5 |
表1模拟工况
Table1.Simulated conditions
壁面润湿性由固液原子间的作用强度决定, 流体在壁面铺展得越充分, 接触角越小, 可认为润湿性越好, 对于光滑壁面而言, 润湿性与固液势能强度之间的关系为
εwf/ε | 6.0 | 1.5 | 1.0 | 0.75 | 0.5 |
θ/(°) | 0 | 0 | 0 | 60 | 90 |
表2不同固液势能强度下对应的接触角
Table2.Corresponding contact angle under different solid-liquid potential energy intensity
2
3.1.静态流体中的界面热阻
图3(a)和图3(b)分别展示了工况1和工况2中液态氩的温度分布. 可以发现流体温度分布是从热壁面侧向冷壁面侧逐渐递减, 并且近壁面处的温度波动和预测得到的壁面附近流体温度显示了流体温度与壁面温度之间的不连续性. 对于工况1, 如图3(a)所示, 当εwf/ε = 6.0时, 近壁面附近流体温度波动最大, 而在x = 2.11 nm到x = 8.16 nm之间, 可以观察到温度分布呈线性分布. 随着固液势能强度的减弱, 热壁面附近流体温度也逐渐趋向于线性. 同样在工况2中, 如图3(b)所示, 通道主流区域温度分布呈线性分布, 而随着冷壁面固液势能强度的减弱, 冷壁面附近流体温度波动逐渐减缓, 并趋向于线性分布. 这是由于壁面润湿性的减弱会降低近壁面处的流体密度波动, 如图4所示, 其中图4(a)和图4(b)分别表示εwf/ε = 6.0和εwf/ε = 0.75时, 工况1和工况2下近热壁面和近冷壁面处的流体数密度分布. 可以发现, 当εwf/ε = 6.0时, 冷、热壁面附近液态氩的密度波动很大, 相邻两层中密度的巨大变化可能导致两个层中相应的温度变化; 而εwf/ε = 0.75时. 冷、热壁面附近流体密度波动较小, 对应的温度分布较为均匀. 此外, 对应于图3(a)和图4(a), 当εwf/ε = 6.0时, 可以发现在x = 2.0 nm位置处流体的数密度接近于零, 说明该位置上存在极少或不存在液体分子, 导致在该位置上统计流体温度时会出现很大的误差, 为了避免该现象, 本节并没有显示该位置上流体的温度分布.图 3 液态氩在不同固液势能强度εwf/ε下的温度分布 (a)工况1; (b)工况2
Figure3. Temperature distribution of liquid argon at different solid-liquid potential energy εwf/ε: (a) Condition 1; (b) condition 2.
图 4 (a)工况1中近热壁面附近流体数密度分布; (b)工况2中近冷壁面附近流体数密度分布
Figure4. (a) Number density distribution near the hot wall surface in condition 1; (b) number density distribution near the cold wall surface in condition 2.
为准确地计算出图3中冷热壁面附近的温度跳跃, 利用连续性热传导方程进行推算, 结果如图5所示. 对于工况1, 随着固液势能强度的增强, 热壁面处的温度跳跃逐渐减小, 而冷壁面处的温度跳跃却略有些增大, 这是因为热壁面润湿性的增强会增大整体流体温度分布的梯度, 由(2)式可知, 冷壁面润湿性不变的情况下, 温度梯度的增大会增大冷壁面处的温度跳跃. 同样, 工况2中随着冷壁面处温度跳跃的降低, 热壁面处温度跳跃略有增大. 同时, 可以发现工况1中热壁面处的温度跳跃比工况2中冷壁面处的温度跳跃高, 且随着壁面润湿性的增强, 这种差异逐渐减小, 说明壁面温度也可以影响界面温度跳跃, 但影响力随着润湿性的增强而减弱.
图 5 固液势能强度对温度跳跃的影响
Figure5. Effect of solid-liquid potential energy on temperature jump.
计算系统的热通量是求解界面热阻的关键, 图6展示了图3中各个条件下计算的热通量. 可以看出, 随着固液势能强度的增大, 两种工况下的热通量均有所增大, 且增大的趋势近似为线性. 这是因为壁面润湿性的增强会强化壁面与流体间的热量传递. 同时, 除 εwf/ε = 6.0和0.5, 工况1中的热通量要略高于工况2中的热通量, 说明壁面温度对系统热通量也有一定的影响, 但效果不明显.
图 6 固液势能强度对系统热通量的影响
Figure6. Effect of solid-liquid potential energy on heat flux of system.
界面热阻为温度跳跃与热通量的比值, 图7显示了固液势能强度对界面热阻的影响. 结果显示, 随着固液势能强度的增强, 工况1中热壁面的界面热阻与工况2中冷壁面的界面热阻均减小, 该趋势与温度跳跃的变化一致. 需要注意的是, 在本文模拟的条件下, 仅当固液势能强度最弱时(εwf/ε = 0.5), 才可以明显观察到工况1中热壁面处的界面热阻要高于工况2中冷壁面处的界面热阻, 在其他条件下, 壁面温度对界面热阻的影响基本可以忽略. 由于工况1中冷壁面和工况2中热壁面润湿性始终不变, 其对应的界面热阻值也基本处于恒定. 由此可以得出结论, 壁面温度对界面热阻的影响受到壁面润湿性的限制.
图 7 固液势能强度对界面热阻的影响
Figure7. Effect of solid-liquid potential energy on interface thermal resistance.
固体壁面与近壁面流体间的传热决定了界面热阻的大小[21,22], 特别是对于近壁区第一层液体原子而言, 其排列方式及数量分布能够很好地反映界面热阻. 图8展示了工况1中热壁面和工况2中冷壁面上吸附流体分子的快照图. 可以明显发现, 固液势能强度的增大能有效吸附更多的流体分子, 且流体分子的排列更加规整, 当εwf/ε ≥ 1.5时, 流体分子的排列形式逐渐接近于固体壁面, 即“类固体”行为. 壁面吸附流体分子越多, 越规整, 有利于固液界面间的传热, 降低界面传热的阻力. 同时, 将图8中吸附流体分子的数量具体化, 以观察壁面温度对界面热阻影响的本质, 结果如图9所示. 当εwf/ε ≤ 1.5时, 工况1中热壁面上吸附的流体分子数低于工况2中冷壁面上吸附的分子数量, 从分子运动角度来分析, 吸附在高温壁面上的氩分子具有更高的动能, 脱离壁面束缚的可能性更大. 壁面吸附流体分子的减少可能是导致热壁面处界面热阻值大于冷壁面处热阻值的原因.
图 8 壁面吸附流体分快照图
Figure8. Snapshot of wall adsorbed fluid molecules.
图 9 固液势能强度对壁面吸附流体分子数的影响
Figure9. Effect of solid-liquid potential energy on the number of fluid molecules adsorbed on the wall surface.
2
3.2.液态氩的热导率验证
为验证静态流体中计算导热的正确性, 以及进行更深一步的研究, 本节主要就图3中每个模拟工况下对应的氩的热导率进行验证. 图10统计了图3中所有条件下液态氩的平均温度(Tave), 可以发现, 随着固液势能强度的增大, 工况1中流体的平均温度逐渐增大, 而工况2中流体的平均温度逐渐减小, 总的变化区间在97.54—130.82 K范围内. 利用傅里叶定律图 10 液态氩的平均温度分布
Figure10. Average temperature distribution of liquid argon.
图 11 液态氩的热导率验证
Figure11. Thermal conductivity verification of liquid argon.
2
3.3.动态流体中的界面热阻
本节主要以工况1为研究对象, 施加外部力形成Poiseuille流动, 并通过设置不同外部力来研究动态流体下固液界面热阻的变化趋势, 并与静态流体下的界面热阻相比, 总结相关规律.图12显示了不同外力(F)作用下纳米通道内的流体速度分布, 包括(a) F = 0.01ε/σ; (b) F = 0.012ε/σ; (c) F = 0.014ε/σ; (d) F = 0.016ε/σ. 结果显示, 在相同外力作用下, 随着热壁面固液势能强度的减弱, 流体流动速度逐渐增大, 且流体分布的对称性逐渐减弱, 流体速度的峰值由通道中心位置逐渐向热壁面侧偏移. 随着外力的不断增大, 流体速度大小也有明显的上升, 但速度分布的趋势基本一致. 图12中显示速度变化主要集中在热壁面侧, 而冷壁面侧由于固液势能强度不变, 速度分布变化很小.
图 12 工况1中液体氩在不同外力作用下的速度分布 (a) F = 0.01ε/σ; (b) F = 0.012ε/σ; (c) F = 0.014ε/σ; (d) F = 0.016ε/σ
Figure12. Velocity distribution of liquid argon under different external forces in condition 1: (a) F = 0.01ε/σ; (b) F = 0.012ε/σ; (c) F = 0.014ε/σ; (d) F = 0.016ε/σ
通过上述对流体速度场的介绍, 可以更好地理解通道内的流动情况. 进一步分析了不同外力作用下液氩的温度场分布, 结果如图13所示. 由图可知, 液态氩的温度分布依然是从高温壁面向低温壁面处递减, 但随着外力的增大, 温度分布形式明显变化. 如图13(a)所示, 当外力较小时(F = 0.01ε/σ), 流体温度分布与静态流体温度分布相似, 通道主流区温度场变化呈线性分布. 随着外力的增大, 主流区温度场逐渐偏离线性, 温度偏离线性分布的部分为靠近热壁面的一端, 这是因为冷壁面的润湿性恒定且较强, 外力的变化对壁面润湿性较强时的温度分布影响较小, 相反地, 随着热壁面润湿性的减弱, 外力的增大会明显提高近热壁面区域的流体温度, 如图13(b)—(d)所示. 同时可以发现, 随着固液势能强度的减弱, 热壁面附近流体温度偏离线性区域的范围也越大, 当固液势能强度较弱时(εwf/ε = 0.5), 该偏离位置距离壁面约1.19 nm. 这是由于壁面润湿性较弱时界面处更容易产生摩擦粘性热, 使得近热壁面处出现温度升高以及局部脱离线性分布的现象. 且从图13中可以看出, 随着外力的增大, 不同固液势能强度下的流体温度分布差距逐渐减小, 特别是图13(d)中, εwf/ε = 0.5时对应的流体温度与εwf/ε = 0.75和εwf/ε = 1.0时对应的温度分布基本重合.
图 13 工况1中液体氩在不同外力作用下的温度分布 (a) F = 0.01ε/σ; (b) F = 0.012ε/σ; (c) F = 0.014ε/σ; (d) F = 0.016ε/σ
Figure13. Temperature distribution of liquid argon under different external forces in condition 1: (a) F = 0.01ε/σ; (b) F = 0.012ε/σ; (c) F = 0.014ε/σ; (d) F = 0.016ε/σ.
图14直观展示了相同固液势能强度下, 不同外力作用对流体温度分布的影响. 包括(a) εwf/ε = 6.0; (b) εwf/ε = 1.5; (c) εwf/ε = 1.0; (d) εwf/ε = 0.75; (e) εwf/ε = 0.5. 结果显示, 无论固液势能强度取值多大, 流体温度随着外力的增大均有着不同程度的上升. 当固液势能强度较强时(εwf/ε = 6.0), 外力的增大使流体温度略有上升, 且分布形式有偏离线性分布的趋势. 随着固液势能强度的减弱, 外力对流体温度分布的影响逐渐明显, 即不同外力间温度分布差距越来越大, 且温度分布偏离线性的现象越来越明显. 特别是当固液势能强度较弱时(εwf/ε = 0.5), 随着外力的增大, 距离壁面1.19 nm范围内的流体温度分布明显脱离线性, 这是由于较弱固液势能强度下, 固液界面处更容易产生摩擦粘性热导致的. 为了定量地描述1.19 nm范围内界面摩擦对局部温度的贡献, 图15显示了εwf/ε = 0.5时, 不同外力作用下1.19 nm范围内局部平均温度分布. 结果显示, 随着外力的增大, 局部区域流体平均温度逐渐上升, 上升幅度分别为9.37%, 10.68%, 15.81%和21.35%.
图 14 工况1中液体氩在不同固液势能强度下的温度分布 (a) εwf/ε = 6.0; (b) εwf/ε = 1.5; (c) εwf/ε = 1.0; (d) εwf/ε = 0.75; (e) εwf/ε = 0.5
Figure14. Temperature distribution of liquid argon under different solid-liquid potential energy intensity: (a) εwf/ε = 6.0; (b) εwf/ε = 1.5; (c) εwf/ε = 1.0; (d) εwf/ε = 0.75; (e) εwf/ε = 0.5.
图 15 不同外力作用下近热壁面局部(1.19 nm)平均温度分布
Figure15. Local (1.19 nm) average temperature distribution near the hot wall under different external forces.
为了说明产生上述温度分布的原因, 本节计算了工况1中热壁面附近液态氩的速度分布, 目的是为了比较不同外力和固液势能强度下产生的滑移现象, 结果如图16所示. 图中仅显示了距离壁面厚度为0.4 nm的流体区域内速度分布, 结果表明, 当固液势能强度εwf/ε < 1.5时, 随着外力的增大, 近热壁面处流体的速度逐渐增大, 即产生了明显的界面滑移, 并导致界面处产生摩擦粘性热. 而在较强的润湿性下, 由于壁面分子与流体分子间较强的作用力, 使得近热壁面处流体的速度极小, 甚至为零, 即在界面处没有速度滑移, 说明此时在界面处没有产生额外的粘性热或产生的粘性热极少. 这与图14中观察到的现象一致, 即壁面润湿性较弱时, 外力的增大会明显增大界面滑移, 产生明显的温度上升现象, 而在润湿性较强时, 外力的增大不会带来界面滑移, 同时界面处温度变化幅度也不再明显.
图 16 近热壁面处液态氩的速度分布
Figure16. Velocity distribution of liquid argon near the hot wall surface.
图17为工况1中, 不同外力与固液势能强度对系统热通量的影响. 如图所示, 随着外力的增大, 每个固液势能强度下对应的热通量均有所上升, 因为在较高外力作用下, 系统产生了较大的额外能量, 流体分子获得了较大的动能, 固液界面处产生的摩擦粘性热使系统温度升高. 此外可以发现, 随着固液势能强度的增大, 外力对热通量的影响逐渐减小, 这是由于热壁面与流体之间较高的势能强度限制了界面处流体分子的运动, 从而降低或消除了流体与壁面间的摩擦粘性热.
图 17 外力作用下固液势能强度比对系统热通量的影响
Figure17. Effect of solid-liquid potential energy on system heat flux under external force
图18为工况1中, 不同外力与固液势能强度对温度跳跃的影响. 图18(a)显示了热壁面处的温度跳跃(ΔThot), 结果表明, 随着外力的增大, 热壁面处的温度跳跃逐渐减小, 这是因为在所选取的外力范围内, 由Navier型温度边界条件预测的界面处的流体温度均未超过热壁面的温度(140 K), 因此, 外力越大, 界面处的流体温度越接近壁面温度, 并导致温度跳跃的减小. 与静态流体中相同的是, 在固液势能强度较强的条件下, 呈现出较低的温度跳跃. 这是因为虽然在较弱润湿性下, 外力的存在更容易产生摩擦粘性热, 但其产生的额外热量始终小于较强润湿性下壁面与流体间导热产生的热量, 说明了纳米尺度下, 固液原子间作用力是决定界面传热的主导因素. 图18(b)显示了工况1中冷壁面侧的温度跳跃(ΔTcold), 结果表明, 随着壁面润湿性的增强, 冷壁面侧的温度跳跃略有增大; 而随着外力的增大, 冷壁面侧的温度跳跃也有所增大, 但幅度很小, 这可能是由于流体分子内部碰撞而产生了少量的热量, 使得近冷壁面流体温度也略有上升, 增大了与冷壁面间的温度差异.
图 18 工况1中不同外力作用下固液势能强度对温度跳跃的影响 (a) 热壁面; (b) 冷壁面
Figure18. Effect of solid-liquid potential energy on temperature jump under different external forces in condition 1: (a) Hot wall surface; (b) cold wall surface.
根据上述对系统热通量及温度跳跃的分析, 得到了不同外力作用及固液势能强度对界面热阻的影响规律, 结果如图19所示. 其中图19(a)表示热壁面侧界面热阻的变化规律, 可以看出, 随着外力的增大, 热壁面处界面热阻在每个对应的固液势能强度下均有所降低; 同时固液势能强度的增大也有利于降低界面热阻. 当固液势能强度较强时, 各个外力下对应的热阻值趋于一致, 说明在较强壁面润湿性条件下, 外力对界面热阻的影响较小. 图19(b)为冷壁面侧界面热阻的变化规律, 从图中可以发现, 冷壁面侧的界面热阻变化幅度很小, 且规律性较差.
图 19 工况1中不同外力作用下固液势能强度对界面热阻的影响 (a) 热壁面; (b) 冷壁面
Figure19. Effect of solid-liquid potential energy on interface thermal resistance under different external forces in condition 1: (a) Hot wall surface; (b) cold wall surface.
为了探究在外力存在条件下, 壁面吸附流体分子数量对界面热阻的影响, 计算了不同外力下热壁面上吸附的氩分子数量, 结果如图20所示. 从图中可以看出, 固液势能强度的增大可以增加壁面吸附流体分子的数量, 从而导致界面热阻减小, 这与之前所得结论一致. 而外力的变化对壁面吸附流体分子数量的影响较小, 这是因为壁面吸附流体分子数量主要是由固液势能强度决定的. 当固液势能强度较强时, 外力作用很难使界面流体分子脱离壁面的束缚, 界面滑移现象不明显; 当固液势能强度较弱时, 外力作用可能使流体分子脱离壁面束缚, 并获得一定的动能, 但其运动位置始终在固液界面处, 导致固液界面处流体分子数量基本不变, 同时产生明显的滑移现象.
图 20 动态流体中热壁面吸附流体分子数
Figure20. The number of adsorbed fluid molecules on the hot wall in flowing fluid.