尺寸效应对微通道内固液界面温度边界的影响*

张龙艳 徐进良 雷俊鹏

(华北电力大学, 低品位能源多相流与传热北京市重点实验室, 北京 102206)

(2018 年10 月19日收到; 2018 年11 月19日收到修改稿)

采用非平衡分子动力学方法模拟不同浸润性微通道内液体的传热过程, 分析了尺寸效应对固液界面热阻及温度阶跃的影响. 研究结果表明, 界面热阻随微通道尺寸的变化可分为两个阶段, 即小尺寸微通道的单调递增阶段和大尺寸微通道的恒定值阶段. 随着微通道尺寸的增加, 近壁区液体原子受对侧固体原子的约束程度降低, 微通道中央的液体原子自由移动, 固液原子振动态密度近似不变, 使得尺寸效应的影响忽略不计.上述两种阶段的微通道尺寸过渡阈值受固液作用强度与壁面温度的共同作用: 减弱壁面浸润性, 过渡阈值向大尺寸区域迁移; 相较于低温壁面, 高温壁面处的过渡阈值更大. 增加微通道尺寸, 固液界面温度阶跃呈单调递减趋势, 致使壁面温度边界和宏观尺度下逐渐符合. 探讨尺寸效应有助于深刻理解固液界面能量输运及传递机制.

1 引 言

近年来, 随着微纳技术的高速发展, 微纳设备中流体传热传质机理的研究日益受到学术界的广泛关注[1-3]. 当热流流过固液界面时, 界面热阻的存在导致界面处发生温度不连续现象. 界面热阻的影响程度可由热阻长度LK来衡量(即Kapitza长度其中为固液界面处的阶跃温度,为界面处流体侧的温度梯度. 实验测量和分子模拟的研究结果都表明, 在液固相互作用较弱的界面, 热阻长度可达数十纳米[4-7]. 在常规尺度下, 固液界面处的温度边界条件直接采用固壁边界, 不考虑温度阶跃. 但是在微纳结构中,微尺度效应的存在使得能量输运及传递规律明显区别于常规尺度[8,9], 导致各种物理参量与系统尺寸有着紧密联系. 因此, 研究微纳系统中界面热阻及界面温度阶跃随尺寸的变化关系具有重要的物理意义与实际应用价值.

现有研究表明, 界面热阻和温度阶跃受到众多因素的影响, 包括固液相互作用强度[10-13]、固体结构[14]、表面粗糙度[15]、壁面温度[16,17]、流体剪切率[18]等. 在数值模拟方面, 分子动力学方法作为探测微纳尺度现象的有效手段, 被广泛应用于界面热阻的研究. Liang和Tsai[19]探究了液膜厚度小于1.3 nm时的界面传热过程, 发现固液界面热阻随液膜厚度的增加而急剧增大. Wang等[20]模拟体系中狭缝尺寸为0.361—1.62 nm, 得到与上述类似的规律. 而Kim等[21]却发现当固壁间距为3.24—12.96 nm时, 固液界面热阻并未发生明显的变化. 尽管关于尺寸效应对界面热阻的影响已有一定数量的研究报道, 但是受计算资源的限制, 目前的模拟体系尺寸多数限制在10 nm的数量级. 在实验研究方面, 原子力显微镜等微纳测量技术的发展, 促使界面热阻的实验测量更加精细化. 目前比较先进的实验方法包括TTR (transient thermo reflectance)技术和3ω法[22], 其测试材料的最小尺寸在100 nm数量级.

综上所述, 现有数值模拟研究体系与实验测量之间存在较大的尺寸跨度, 该范围内的界面热阻变化情况很少涉及, 关于温度边界的研究仍处于空白. 此外, 对于小于10 nm的微通道内, 已有模拟结果中, 界面热阻随微通道尺寸的变化趋势存有争议, 观点尚未统一. 因此, 针对尺寸效应如何影响固液界面热阻及温度阶跃有待进一步探究. 本文利用非平衡分子动力学方法模拟在给定温差下不同浸润性微通道内液体传热的过程, 分析了尺寸效应对固液界面热阻及温度阶跃现象的影响, 并深入剖析界面热阻作用机制. 研究结果为从物理意义角度揭示温度边界条件的“宏观”与“微观”尺度的内在机理提供了有益启发. 本模拟采用开源分子动力学模拟软件LAMMPS实现.

2 物理模型及模拟细节

图1为物理模型和界面热阻长度定义示意图.系统沿着x和y方向均采用周期性边界条件,z方向为固定边界条件. 模拟体系尺寸为(h为变化量,σ为氩原子的尺寸参数), 研究表明, 在x和y方向进一步增大尺寸对计算结果没有影响. 固体壁面原子Pt按照面心立方(FCC)晶格排列, 初始密度为对应的晶格常数为1.15σ, 其(100)晶面与液体原子接触, 共有4800个固体原子. 上下壁面每一侧的固体原子共有13层, 靠近流体区域的10层固体原子构成热墙, 与液体区域进行热量传输. 固体板之间充满液体氩. 初始液氩原子也按照FCC晶格排列方式布置, 初始密度为对应晶格常数为1.72σ, 共有2366—10829 (随着h发生变化)个液体原子.

液体氩原子之间的相互作用采用Lennard-Jones (LJ)势能模型, 表达式为

式中r为原子间的距离; 液体氩原子之间的尺寸参数为能量参数为原子质量固体铂原子之间的相互作用同样采用LJ势能模型, 尺寸参数为能量参数为其中下标s表示固体. 液体原子与固体原子的势能参数均来自Nagayama等[23]的研究报告.

图1 (a) 物理模型图; (b) 界面热阻长度定义Fig.1. (a) Physical model of system; (b) definition of thermal resistance length at interface.

固液之间的势能作用通过对LJ势能模型进行修正[24]:

采用Velocity-Verlet算法求解运动方程, 时

为了控制固体墙的温度, 在热墙的最外层增加3层虚拟原子, 以模拟其无限大势场. 虚拟原子类似热墙的外延扩展. 最外两层虚拟原子为固定原子, 运行过程中使其固定在晶格位置处保持原子静止不动, 以便维持系统的稳定. 靠近热墙的虚拟原子用来模拟与热墙链接的虚拟热源, 热浴的实现采用Langevin方法[25], 对其额外施加一个随机力和阻尼力, 并使其满足三维郎之万方程. 这两个力共同确保实现固体墙及此层固体原子的温度恒定. 其动量方程表示为

式中pi为第 i个原子的动量矢量;为阻尼常数; f(t)为原子之间的相互作用力; F(t)为随机作用力. 可以从具有零均值和标准偏差的高斯分布中进行简化, 其标准偏差为在本文中, 考虑固体原子Pt的简谐振动模型, Pt原子的相互作用力为简谐力, 认为固体原子均在其晶格固定位置附近做简谐振动. 弹簧系数为k = 46.8 N/m. 在模拟中实际固体原子与虚拟原子之间靠虚拟弹簧相连, 在水平方向弹簧系数为0.5k, 在垂直方向为2k. 虚拟原子与固定原子之间的弹簧系数在水平方向为2k, 垂直方向为3.5k.

3 结果与讨论

3.1 温度和密度分布

微通道内液体分布状态影响固液界面能量传递, 引起固液界面温度边界条件的变化. 为探索微通道尺寸效应对温度边界的作用机制, 首先需要分析微通道内液体的温度及密度分布特征. 本文为了观察统计液体的密度和温度分布, 将液体区域沿着z方向划分为n层(n随微通道尺寸H发生变化), 第 l个切片层 (1≤ l ≤ n)在 JStart—JEnd时间段内的无量纲粒子密度

式中Nk为第l层的液体原子数目为各液体分层的厚度; A为液体计算区域在xy平面的面积,

第 l个切片层 (1 ≤ k ≤ n)在 JStart—JEnd时间段内的温度为

图2 固液势能参数 α =0.14 时, 流体沿着z方向的温度分 布 和 密 度 分 布 (a) H =8.07σ ; (b) H =57.65σ ;(c)H=103.77σFig.2. Temperature and density profile of liquid along z direction when solid-liquid potential energy parameter α=0.14: (a) H =8.07σ ; (b) H =57.65σ ; (c)H=103.77σ.

观察图2液体的密度分布特征, 分别从近壁区和通道中央两个方面进行分析. 一方面, 近壁区液体原子受固体壁面作用较强而呈现类固体排列方式, 液体密度表现出明显的分层现象. 当微通道尺寸H=8.07σ时, 近热壁区内液体无量纲原子密度分布的第一个波峰值约为冷壁区约为当微通道尺寸增大到H=57.65σ时,液体原子密度波峰值明显降低, 热壁区约为冷壁区约为近壁区液体原子振荡区间宽度约为 6.34σ. 对于微通道H=103.77σ, 热壁区液体无量纲密度波峰值2.358, 冷壁区近壁区密度振荡区间宽度约为6.27σ. 比较上述两种大尺寸微通道内流体近壁区分布微观状态, 发现改变微通道尺寸, 对近壁处振荡液体区的范围及密度波动幅度的影响几乎忽略不计. 另一方面, 对于小尺寸微通道H=8.07σ, 通道中央的液体原子排布仍呈现出微弱的振荡现象, 这是由于两侧固壁原子的共同作用引起的; 而对于大尺寸微通道H=57.65σ及 1 03.77σ,固体壁面原子对微通道中央区域液体原子的约束程度微弱. 液体原子在微通道内可以自由移动, 原子均匀排布, 其平均密度均约为0.79σ3.

综上所述, 微通道尺寸的变化引起内部流体温度和密度分布的改变. 为了分析固液界面传热与微通道尺寸的关系是否受其他因素的影响, 进一步探究了壁面浸润性的作用机制. 图3为微通道尺寸为H=46.12σ时, 不同固液势能参数下液体无量纲密度和温度分布. 图中黑色虚线表示固液界面无温度阶跃时液体的理论温度分布. 固液势能参数决定了固体原子与液体原子的相互作用强度, 它不但改变近壁区液体原子的分布状态, 还会影响固液界面处的能量传递. 当固液势能参数为α=1.0 时, 固液原子相互作用最强; 当α=0.14 时, 固液原子相互作用最弱. 由图3(a)的密度分布可以看出, 在近壁区流体排布结构呈现出类固体的层状排列方式, 当势能参数α=0.14 时, 其热壁面处的四层流体数密度峰值依次为1.413 , 1.036 和 0.876 ,根据已有文献资料[26]可知, 当温度为100 K (无量纲温度为 0.826 )时, 晶体氩的密度为1.48 g/cm3(无量纲数密度为ρ∗=0.881 ), 对比发现, 近壁处的第一层氩原子密度远高于结晶密度, 甚至第四层也与晶体密度相近, 因此, 近壁区的氩液体原子表现出类固体的结构特征. 增强壁面浸润性,近壁区氩原子无量纲数密度波动幅度增大, 即其类固体排列特征更加显著. 然而, 密度振荡区域的范围近似不变, 均为 Δz∗ ≈6.34 , 壁面浸润性的影响甚微. 观察图3(b)液体无量纲温度分布可知, 当α=1.0 时, 固液原子传递热量效率最高, 界面处温度阶跃最低, 尤其在热壁面处几乎与固体壁面温度相等.液体内部温度分布斜率也接近无温度阶跃情况下的温度梯度. 当α=0.14 时, 液体温度曲线变平缓,其明显偏离无阶跃温度分布状态. 在热壁面处无量纲温度阶跃约为0.0796, 是α=1.0 时的4.14倍.

图3 微通道尺寸 H=46.12σ时流体沿着z方向的(a) 密度分布和(b)温度分布Fig.3. (a) Density profile and (b) temperature profile of liquid along zdirection when microchannel height H= 46.12σ.

3.2 温度阶跃

为了直观地描述微通道尺寸效应对固液界面温度阶跃的影响, 图4给出了固液势能参数α=0.14, 2.0时, 固液界面温度阶跃与微通道尺寸的变化规律. 由图4可知, 增加微通道尺寸, 界面温度阶跃随之减小, 微观尺度下固液界面处的阶跃温度边界逐渐向宏观尺度的无阶跃温度边界靠拢.尤其, 强浸润性壁面α=0.2 在微通道尺寸H=103.77σ时, 热壁面处无量纲温度阶跃为0.01175, 界面处液体温度与壁面温度的比值两者几乎接近. 此外, 界面温度阶跃还与壁面浸润性和壁面温度有关. 在壁面温度相同的情况下, 液势能参数α=0.2 的相互作用力更强, 其温度阶跃总小于α=0.14 的值. 当壁面浸润性相同时, 热壁面处的温度阶跃总小于冷壁面的值. 这是因为温度升高, 近壁区的液体原子具有更高的能量挣脱固体原子的束缚, 从而与固体壁面碰撞更加频繁, 能量交换能力增强[27]. 因此, 在高温壁面处界面热阻较小, 固液界面温度阶跃相应的较低.

图4 固液界面温度阶跃随微通道尺寸的变化Fig.4. Variation of temperature jump at solid-liquid interface with microchannel height.

根据液体的温度分布, 给出液体实际温度梯度与微通道尺寸H和冷热壁面处界面热阻长度LK,c,LK,h的理论表达式为

此式表明在界面热阻长度LK不随通道尺寸变化的情况下, 微通道尺寸H越大, 实际温度梯度越接近理论温度梯度, 即在无温度阶跃的情况下的液体温度梯度. 图5展示了实际温度与理论温度的接近程度随微通道尺寸的变化关系, 其中表示温度梯度, 图中纵坐标表示液体实际理论温度梯度与理论温度梯度的差值比理论温度梯度值. 由图5可知, 在小尺寸微通道H=5.19σ时, 实际温度梯度比理论温度梯度几乎小1/2. 增大微通道尺寸, 液体实际温度分布与接近理论状态之间的区别逐渐减小. 综上所述, 随着微通道尺寸的增大, 固液界面温度阶跃呈单调递减趋势, 内部流体实际温度梯度与理论梯度无限逼近, 致使壁面温度边界和宏观尺度下逐渐吻合.

图5 液体实际温度与理论温度的接近程度随微通道尺寸的变化关系Fig.5. Degree of closeness of actual and theoretical temperature of liquid versus microchannel height.

3.3 界面热阻

在宏观尺度下, 无阶跃温度边界被认为是理所当然的; 而微尺度下, 由于界面热阻的存在, 固液界面处发生温度不连续现象, 使得微尺度下温度边界条件明显区别于常规尺度. 如果从物理意义的角度区分宏观与微观的温度边界条件界限, 那么就需要考虑两者之间是否具有内在的联系, 并剖析界面热阻对温度边界的作用机理. 探讨尺寸效应对固液边界热阻的影响, 将有助于加对强温度边界条件物理本质的理解.

图6为不同浸润性下, 冷热壁面处界面热阻长度随微通道尺寸的变化. 由图6(a)可知, 对于较强浸润性壁面, 即固液势能参数α＞0.3 , 当微通道尺寸小于10.95σ的尺寸范围时, 增加微通道尺寸,界面热阻将随之略微增大. 对于较弱浸润性壁面,即势能参数α=0.2 , 0.14的情况, 界面热阻随微通道尺寸的变化更加显著, 其在小尺寸范围内迅速上升. 直至微通道尺寸分别达到12.68σ, 23.06σ时,界面热阻的变化与微通道尺寸近似无关, 尺寸效应可以忽略不计. 总结不同浸润性下微通道尺寸效应对固液界面热阻的影响规律, 可将其归纳为两个部分: 小尺寸范围的单调递增趋势与大尺寸范围的恒定值规律. 上述两种规律的尺寸过渡阈值随壁面浸润性的减弱向大尺寸范围迁移. 观察图6(b), 可以发现冷壁面处的界面热阻变化规律与图6(a)中热壁面的一致, 区别在于微通道尺寸过渡阈值对壁面浸润性的依赖程度减弱, 尤其对于弱固液相互作用, 即势能参数α=0.14 的情形, 其微通道尺寸过渡阈值降低为17.29σ, 相比热壁面, 约小于5.77σ,产生明显的差异. 此外, 保持微通道尺寸不变, 界面热阻随壁面浸润性增强而减小; 而且对于相同浸润性及微通道尺寸的情况, 热壁面处界面热阻总是低于冷壁面处界面热阻值, 该结论与文献[28]结果一致.

图6 固液界面热阻随微通道尺寸的变化 (a) 热壁面; (b) 冷壁面Fig.6. Variation of thermal resistance at solid-liquid interface with microchannel height: (a) Hot wall; (b) cold wall.

固液界面热阻主要受固体与近壁区流体之间传热的影响[19-21], 因此对其机理分析需从近壁区液体原子与固体结构特征出发. 由微通道中流体密度分布特征可以看出, 在近壁区液体原子尤其是第一层液体原子按类固体方式排列, 由固体壁面向氩薄膜传递热量的过程可以认为是载热子的传递过程. 在小尺度下, 微通道中的氩薄膜受两侧固体壁面的影响, 呈现出显著的类固体排列结构特征. 如图2(a)所示, 随着通道尺寸的增大, 近壁区氩原子受对侧壁面的影响减小, 导致氩原子排布的无序程度增加, 使得其类固体结构特征减弱. 因此, 在热量传递过程中, 增加纳米通道的高度, 载热子的边界散射增强, 通过弹道输运的载热子相应地减少.此外, 增加微通道尺寸, 引起短程载热子不能弹道输运穿过界面, 从而增强了散射程度, 使得弹道效应进一步减弱[29,30]. 由此推断得到, 界面热阻随微通道尺寸的增加而增大. 当微通道尺寸增加到一定程度, 近壁区氩原子受对侧固体原子的影响忽略不计, 其类固体排列结构几乎保持不变. 该结论可由图2(b)和图2(c)密度分布得到验证, 在大尺寸通道中, 近壁区氩原子密度分布的振荡幅度及振荡范围变化甚微, 原子排列的有序程度几乎保持不变.载热子穿过界面时的散射程度不再随尺寸变化, 界面热阻维持不变. 改变壁面浸润性, 如图3(a)所示, 固液势能作用越强, 则近壁区氩原子排列越有序, 从而削弱了载热子的界面散射, 即增强了弹道效应. 因此, 当微通道尺寸相同时, 界面热阻随壁面浸润性的增强而减小.

综上所述, 界面热阻由载热子的散射程度决定. 目前, 对于界面热阻的微观机理解释主要有声失配理论(AMM)和散射失配理论(DMM)等. 根据上述理论, 界面热阻是由相邻两种材料的原子之间振动耦合, 即振动态密度(VDOS)失配程度决定的. 固液原子之间的声子VDOS失配程度越低,更多的能量将以简谐振动的方式传递, 界面热阻将越小[31]. 声子VDOS可由速度自相关函数经过傅里叶变换得到:

式中v为原子速度;〈〉表示对不同时间点的统计平均;w为振动频率. 因此固体与近壁区液体原子之间的VDOS失配程度可表示为

图7 固液势能参数为 α=0.14 时固液原子VDOS分布Fig.7. VDOS profile of solid atoms and liquid atoms located at interfaces when solid-liquid potential energy parameter α=0.14 .

计算得到不同微通道尺寸下的固体与近壁区液体原子VDOS分布如图7所示. 由图7可知, 固体壁面原子的VDOS的波峰分布在高频区域, 而近壁区液体原子的VDOS向低频区域迁移. 这是因为相较于固体原子而言, 液体原子之间的结合力更弱. 对于小尺寸微通道, 内部液体受两侧固壁原子的约束作用显著, 通道中液体的密度振荡特征明显, 弹道声子传输占主要地位, 界面处以简谐振动的方式传递能量, 固体与液膜之间的VDOS失配程度降低, 导致固液热阻减小. 增大微通道尺寸,处于微通道中央的液体原子几乎不受固体原子影响, 能够在微通道内自由移动, 近壁区液体原子VDOS受对侧壁面固体原子的束缚甚微, 固液原子振动耦合程度保持不变. 当微通道尺寸为H=34.59σ与H=103.77σ时液体原子VDOS分布几乎没有差异, Δ VDOS 未发生变化; 而H=8.07σ时,固体原子VDOS未发生改变, 液体原子的VDOS波峰值明显下降, 使得固液之间 Δ VDOS 减小.因此, 界面热阻长度随微通道尺寸先增大后保持不变.

4 结 论

采用非平衡分子动力学方法模拟不同浸润性微通道内液体的传热过程, 对液体密度和温度分布展开研究, 分析了尺寸效应对固液界面热阻及温度阶跃的作用机制, 得到如下结论.

1)界面热阻随微通道尺寸的变化规律可划分为两个阶段: 小尺寸微通道的单调递增阶段, 大尺寸微通道的恒定值阶段. 对于小尺寸微通道, 液体原子受固体原子影响显著, 导致固液原子VDOS失配程度较低, 可以有效减小固液界面热阻. 对于大尺寸微通道, 只有近壁区液体原子受同侧固体原子的约束明显, 远离固体壁面的液体原子则自由移动, 固液原子VDOS几乎保持不变, 导致固液界面热阻不再随着微通道尺寸的变化而发生改变.

2)虽然界面热阻随微通道尺寸的变化趋势与壁面浸润性及壁面温度无关, 但是上述两种阶段的微通道尺寸过渡阈值却受二者的影响: 减弱壁面浸润性, 过渡阈值向大尺寸区域迁移; 相较于高温壁面, 低温壁面处的过渡阈值更小.

3)固液界面温度阶跃随微通道尺寸的增加呈减小趋势, 使其逐渐接近宏观尺度下的无阶跃温度边界条件. 探讨这种尺寸效应有助于深刻理解微观和宏观尺度下的温度边界条件的作用机制, 为微纳设备的制造提供理论支撑.