纳米结构影响沸腾换热特性的分子动力学模拟

张子剑， 陈占秀， 杨 历， 张石重

(河北工业大学 能源与环境工程学院， 天津 300401)

1 引 言

快速沸腾传热技术能够迅速降低表面温度，并带走热量，广泛应用于微电子设备散热领域[1-3]. 表面改性是提高沸腾换热性能的主要手段，其中纳米结构表面作为表面改性的一种重要方式具有增强沸腾换热的巨大潜力，合理设计表面纳米结构几何参数和排布方式能够进一步提高表面沸腾换热特性[4]. 表面纳米结构对固体-流体界面热传递的影响成为目前研究热点.

近年来在微尺度蒸发和沸腾研究方面取得了一些进展，但是由于纳米尺度上物理现象的复杂性，以实验方式研究纳米结构对沸腾换热过程的热量传递机理仍然是一个巨大的挑战. 许多研究已经表明，纳米结构的形貌、尺寸以及排列方式是快速沸腾过程中重要影响的因素. 实验结果均表明表面纳米结构可以显著增强沸腾传热性能[5-8]. 随着计算机科学的发展，分子动力学模拟成为研究沸腾现象的有力工具. Morshed等人[9]利用分子动力学模拟研究了圆柱形纳米结构高度对沸腾换热的影响，发现固液之间的传热速率随着纳米结构高度的增加而增加. Seyf等人[10]研究了锥形纳米结构高度对沸腾换热的影响，得到相同的结论. Chen等人[11]研究了三角形纳米结构高度对液体沸腾过程的影响，发现热流密度随纳米结构的增加而增加，但存在最佳纳米结构高度. Yin等人[12, 13]研究了沉积在固体壁面的纳米颗粒对沸腾换热的影响，认为沉积纳米颗粒能够增强固液之间的作用力，减小固液界面热阻. 唐元政等人[14]研究了有无碳纳米管对液体沸腾换热的影响，发现加入碳纳米管能够促进沸腾换热. Chen等人[15]研究了圆柱体、长方体和锥体纳米结构对沸腾换热的影响，发现圆柱纳米结构是提高核沸腾效率的最佳选择. Zhang等人[16]研究了凸形和凹形纳米结构对沸腾换热的影响，结果发现凹形纳米结构的换热效果要强与凸形纳米结构，凹形纳米结构能够有效的在凹槽内积累热量. 梅涛等人[17]研究了表面粗糙度对流体流动行为的影响，发现粗糙度的增强增加了凹槽内液体原子. Wang等人[18]比较长方体和指状纳米结构表面上水的过冷沸腾性能，在纳米结构内的空腔中观察到了纳米级垂直对流，该对流增强了沸腾传热性能. Qanbaria等人[19]研究了锥形纳米结构阵列对沸腾换热的影响，发现当纳米结构表面与液体接触面积一定时，纳米结构的尺寸越小，沸腾换热效果越好. Liu等人[20]发现通过增加随机粗糙表面的粗糙度可以提高沸腾过程的临界热通量.

综上所述，纳米结构的形貌、尺寸以及排列方式都对液体沸腾换热过程有重要的影响. 对于纳米结构表面粗糙度以及纳米结构排列方式的研究有较少专家涉及[18-20]，学者们仅研究了粗糙纳米结构表面对沸腾换热的影响，尚未有学者研究具有相同粗糙度的纳米结构表面其不同纳米结构排列方式对沸腾换热的影响. 本文通过研究纳米结构表面粗糙度以及栏栅形和棋盘形两种纳米结构排列方式对沸腾换热过程的影响，分析液体沸腾换热过程中的沸腾起始时间以及界面热阻，为微尺度纳米结构沸腾换热提供理论支持.

2 物理模型

物理模型如图1所示，模型尺寸为12.5 nm(X)×12.5 nm(Y)×58.8 nm(Z)，模拟体系由固体区域、液体区域和气体区域三部分构成. 固体区域原子按照面心立方(FCC)晶格排列，晶格常数为0.392 nm. 如图2所示，固体区域分为4层，最底层为固定层，由4096个铂(Pt)原子组成，底层固定以防止底板受热变形. 加热层位于固定层上方，由6144个铂原子组成. 加热层上方为导热层，由14336个铂原子组成，用来与液体进行热交换. 导热层上方为纳米结构层，纳米结构层通过周期性地以单位晶格增量去除原子而形成栏栅形和棋盘形纳米结构，在纳米结构层中共有一半的原子被提取出来. 在所有铂原子上施加一个弹簧力，以确保铂原子在模拟过程中围绕它们自己的初始位置振动，在每个时间步，每个原子上的力的大小是-Kr，其中r是原子从其当前位置到其初始位置的位移. K是弹簧的系数因子，幅度为0.786 eV/nm. 液体层厚度为15.5 nm，液体原子采用氩原子(Ar)，初始氩原子按照面心立方(FCC)晶格排列，晶格常数为0.572 nm，共有52272个液体原子. 最上层为气体区域，气体原子按照面心立方(FCC)晶格排列，晶格常数为3.672 nm，共有768个气氩原子. X和Z方向为周期性边界条件，Y方向是简单的固定性边界条件. 盒子上边界设有反射墙，当有原子撞击到上壁面时，会以反方向原速弹回，在这个过程中没有能量交换和能量损失.

氩原子之间的相互作用由Lennard-Jones(12 -6)势能模型，表达式为：

(1)

对于铂和氩原子的相互作用，采用Lorentz-Berthlot 混合规则：

(2)

(3)

其中ε为能量参数，σ为尺寸参数，两个参数取决于原子的类型. 本文通过控制εPt-Ar的值来改变固体表面的润湿性，润湿性的强弱可以用接触角来衡量，接触角越大润湿性越弱，接触角越小润湿性越强，接触角公式如下：

(4)

本文所采用的原子相互作用参数及对应的接触角如表1所示.

表1 原子相互作用参数及对应的接触角

图3为不同纳米结构表面示意图，黄色原子为纳米结构层原子，蓝色原子为导热层原子. 在建立模型过程中，保证纳米结构层中原子数不变，只改变纳米结构表面粗糙度以及纳米结构排列方式，图中r表示为纳米结构表面的粗糙度，粗糙度定义为纳米结构表面实际面积与纳米结构表面在XY平面投影面积之比[21]，粗糙度计算公式如下：

(5)

(6)

纳米结构表面尺寸参数及表面粗糙度如表2所示.

模拟过程分为平衡阶段，固体表面加热阶段和沸腾换热阶段. 平衡温度设置为90 K，时间步长和截断半径分别为5 fs和3.5σ. 对每个时间步长积分运动方程采用Velocity Verlet算法. 为加速平衡过程，采用共轭梯度法将系统能量最小化， 平衡阶段对所有液体原子采用NVT系综运行2 ns，

表2 纳米结构表面尺寸参数及粗糙度

对加热层施加langevin温度控制器，设置温度为90 K，运行2 ns. 固体表面加热阶段将液体原子上的NVT系综撤掉，加热层采用langevin温度控制器控制壁面温度，设置为300 K，运行1 ns. 最后保持加热层温度不变，整个系统采用NVE系综，再运行10 ns进行数据采集. 本模拟采用开源分子动力学模拟软件LAMMPS[22]，原子可视化采用OVITO软件[23].

3 模拟结果讨论

3.1 不同纳米结构表面液体的快速沸腾过程

本文以栏栅形和棋盘形排列纳米结构表面为研究对象，观察不同纳米结构表面粗糙度以及排列方式对液体快速沸腾过程及其换热特性的影响. 图4为不同纳米结构表面液体快速沸腾过程的快照图. 图4(a)为栏栅形排列纳米结构表面液体快速沸腾过程，图中显示 Case1在0 ps-900 ps液体层逐渐膨胀，在此过程中主要是液体原子蒸发过程，在900 ps时液体在纳米结构的凹槽内出现气核，这与张龙艳等人[24]研究结果气核更容易出现在壁面凹槽内一致，在1000 ps时气核生长形成蒸气层，1000 ps后蒸气层推动液体层向上移动. 比较栏栅形排列纳米结构表面Case1、Case2和Case3三种情况，其表面粗糙度分别为1.25、1.50和2.00，发现蒸气层出现的时间分别为1000 ps、800 ps和600 ps. 图4(b)为棋盘形排列纳米结构表面液体沸腾过程，其表面粗糙度也分别为1.25、1.50和2.00. 从图4(b)中可以看出，随着棋盘形表面粗糙度增加，蒸气层分别出现在900 ps、700 ps和500 ps. 由图4可知，液体快速沸腾过程可以分为蒸发阶段和沸腾阶段. 在蒸发过程中液体从表面吸收热量进行蒸发. 当液体吸收足够多的热量，纳米结构表面出现蒸气层，液体发生沸腾. 随着栏栅形和棋盘形纳米结构表面粗糙度的增加，液体沸腾时间越来越早. 对于粗糙度相同的栏栅形和棋盘形纳米结构表面，棋盘形纳米结构表面液体发生沸腾的时间较早.

图3 纳米结构表面. (a) 栏栅形, (b) 棋盘形Fig. 3 Nanostructure surfaces. (a) fence-shape, (b) checkerboard-shape

图4 不同纳米结构表面液体快速沸腾过程. (a)栏栅形, (b)棋盘形Fig. 4 Rapid boilings of liquid on surfaces of different nanostructures. (a) fence-shape, (b) checkerboard-shape

图5 不同纳米结构表面氩原子动能和势能分布. (a)栏栅形, (b)棋盘形Fig. 5 The distributions of kinetic energy and potential energy of argon atoms on the surfaces of different nanostructures. (a) fence-shape, (b) checkerboard-shape

将模拟系统在XZ平面划分成二维小网格，每个网格尺寸为σAr-Ar×σAr-Ar，每隔10ps统计一次网格内氩原子动能(Ke)和势能(Pe)分布. 图5为栏栅形和棋盘形纳米结构表面氩原子的动能和势能分布云图. 从动能图中可以看出，凹槽上方氩原子动能较大. 从势能分布可以看出，势能值为负值，其绝对值表示为氩原子之间相互作用的大小，势能的绝对值越大，氩原子之间作用力越强，越不容易形成气泡核[24]. 当沸腾开始时，固体表面附近氩原子势能绝对值仍然较大，这是由于亲水表面与氩原子之间的作用力较大，在固体表面附近出现一层类固体层. 由图可知凹槽上方氩原子具有较大动能和较低势能绝对值，为气泡核的形成提供有利条件. 对比Case1和Case4中氩原子动能和势能分布可以发现，在表面粗糙度为1. 25时，棋盘形纳米结构表面氩原子动能更高，出现更多的气泡核. 随着表面粗糙度的增加，增加了固液接触面积，增加了更多的凹槽，增加了沸腾过程的气化核心. 对于具有相同粗糙度的纳米结构表面，棋盘形排列纳米结构提高了纳米结构表面附近液体的动能，使液体能够较早的发生沸腾.

液体在固体表面上的快速沸腾，显著的特点是液体层离开固体表面飞升起来，本文通过统计不同纳米结构表面液体质心位置随时间的变化情况反映液体层离开固体表面的快慢程度，进一步说明不同纳米结构表面对液体快速沸腾过程的影响，液体质心位置截止到高度为30 nm. 图6为不同纳米结构表面液体质心位置(COM)随时间的变化情况. 由图可知，在蒸发阶段，液体受热膨胀导致质心位置缓慢上升. 当发生沸腾时，蒸气层推动液体层迅速上升离开固体表面. 根据质心位置数据采集计算液体沸腾起始时间tOB.每隔20 ps取一次液体质心位置，取t2=t1+Δt，Δt=20 ps，上升距离为Δh=h2-h1.当时间周期为tn-1-tn时，hn-hn-1<0.2 nm，同时在时间周期为tn-tn+1时，hn+1-hn>0.2 nm,则tOB=tn，结合图3液体快速沸腾快照图得到液体沸腾起始时间tOB，如表3所示. 从表3中得出，栏栅形和棋盘形排列纳米结构表面液体沸腾起始时间均随着粗糙度的增加而减小. 对于具有相同粗糙度的不同纳米结构排列的表面，棋盘形纳米结构表面液体质心上升速率更快，液体层离开固体表面的时间更早.

3.2 不同纳米结构表面对液体沸腾换热特性的影响

首先统计了不同纳米结构表面液体沸腾过程中沿液体高度方向的温度分布，如图 7所示，其为模拟时间为100 ps时不同纳米结构表面液体层在z方向的温度分布情况. 液体层温度梯度越大，热量越容易在壁面附近积累，越有利于气泡核的形成[25]. 由图可知，热壁面附近的液体温度率先升高，液体温度沿液体层高度方向逐渐降低. 随着表面粗糙度的增强，栏栅形和棋盘形纳米结构表面附近液体温度均增大，液体层高度方向温度梯度增加. 当表面粗糙度相同时，棋盘形纳米结构表面附近液体温度比栏栅形表面液体温度高，棋盘形纳米结构表面液体温度梯度更大.

表3 不同纳米结构表面液体沸腾起始时间

图6 不同纳米结构表面液体质心位置随时间的变化Fig. 6 Positions of the liquid center of mass on the surfaces of different nanostructures over time

不同纳米结构表面液体温度随时间的变化情况如图8所示. 从图中可以看出，对于不同纳米结构表面上的液体温度均在蒸发阶段迅速增加，而后逐渐趋于平缓. 这是由于蒸发阶段表面液体能够与固体表面进行较为强烈的热量交换；而随着温度的升高，液体发生膨胀，氩原子之间的距离增大，并在固体表面形成气泡核心，之后固体表面附近形成蒸气层，液体层逐渐脱离固体表面，在固体表面与液体层之间形成的蒸发层阻碍了热量向液体的传递，液体温度趋于平缓. 在沸腾起始时间之前，液体温度在相同时间上随着粗糙度的增加而增加. 当栏栅形和棋盘形纳米结构表面粗糙度相同时，在沸腾开始前棋盘形纳米结构表面液体温度均大于栏栅形纳米结构表面液体温度，更高的温度使液体质心位置上升速率加快，使液体层较早的离开固体表面.

图7 不同纳米结构表面液体在Z方向的温度分布Fig. 7 Temperature distributions of liquid on the surfaces of different nanostructures in the Z direction

图8 不同纳米结构表面液体温度随时间的变化Fig. 8 Temperature of liquid on the surfaces of different nanostructures over time

进一步计算了不同纳米结构表面液体层的热通量，热通量采用单位时间内通过液体原子的总能量增量除以液体层的横截面获得. 图9为不同纳米结构表面液体热通量随时间的变化情况，从图中可以看出，热通量整体呈现下降的趋势. 对于栏栅形和棋盘形纳米结构表面，液体初始时刻的热通量均随着粗糙度的增加而增加. 当表面粗糙度相同时，棋盘形纳米结构初始时刻热通量较大，之后热通量下降较快；栏栅形纳米结构表面初始时刻液体热通量值较小，之后热通量下降较慢，出现图示的热通量变化过程中交叉现象；这是由于棋盘形纳米结构表面的液体层较早出现气化核心，蒸气层形成的时间较早，从而影响表面热量向液体层的传递.

图9 不同纳米结构表面液体热通量随时间的变化Fig. 9 Heat flux of liquid on the surfaces of different nanostructures over time

对于微纳尺度的换热，界面热阻对换热率影响很大[26]. 界面热阻r=(TW-TL)/J， 其中TW是固体表面温度，TL是液体温度，J是热通量[14, 25]. 界面热阻是由两种相邻材料在其界面处的振动态密度不匹配引起的[27]. 本文以固体表面上液体沸腾传热过程计算的热通量和固液温差角度反映界面热阻的变化，图10为不同纳米结构表面液体沸腾过程界面热阻. 从图中可以看出，在蒸发阶段界面热阻基本保持恒定值，当液体开始形成气泡，蒸气层逐渐形成，界面热阻逐渐增加. 图11为模拟时间在200 ps内不同纳米结构表面平均界面热阻. 为了与现有文献研究结果进行验证，本文模拟了平板上液体沸腾过程，计算了平板表面的平均界面热阻. 由图11可以看出，本文得出的平板界面热阻和Yin等人研究结果[12]基本一致. 栏栅形和棋盘形纳米结构表面平均界面热阻均随着表面粗糙度的增大而减小. 对于具有相同粗糙度的不同纳米结构表面，棋盘形纳米结构表面的界面热阻比栏栅形纳米结构表面界面热阻小.

图10 不同纳米结构表面液体沸腾过程界面热阻随时间的变化Fig. 10 Interface thermal resistances of liquid on the surfaces of different nanostructures over time

图11 模拟时间在200 ps内不同纳米结构表面液体的平均界面热阻Fig. 11 Average interface thermal resistances of liquid on the surfaces of different nanostructures at 200 ps

4 结 论

本文采用分子动力学方法模拟了液体在不同纳米结构表面的快速沸腾行为，分析了不同纳米结构对液体沸腾换热特性的影响，得到如下结论：

(1)液体在快速沸腾过程分为蒸发过程和沸腾过程，在蒸发过程中，液体与纳米结构表面之间的界面热阻基本不变，热量迅速向液体传递导致液体温度迅速升高. 当液体发生沸腾时，在表面附近出现一层密度稀薄的蒸气层，使界面热阻迅速增大，阻碍表面热量向液体的传递.

(2)纳米结构表面粗糙度对液体沸腾过程有重要影响. 随着表面粗糙度的增加，液体与固体表面接触面积增大，提高了液体初始时刻的热通量，降低了固液平均界面热阻，缩短了液体沸腾起始时间.

(3)纳米结构的排列方式能够影响液体沸腾换热特性，对于具有相同粗糙度的栏栅形和棋盘形两种不同纳米结构表面，棋盘形纳米结构表面在沸腾开始前界面热阻较小，有利于热量向液体层的传递，沸腾起始时间缩短. 表面改性可以通过调整纳米结构排列方式及其粗糙度来改变沸腾起始时间，从而更有利于调整不同固体表面液体快速沸腾换热过程.