不同润湿表面液体沸腾的分子动力学模拟

张石重，陈占秀，杨历，苗瑞灿，张子剑

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

电子器件集成度的日趋增加，使散热问题已经对电子器件发展形成制约。相变过程被认为是解决高热流散热的重要途径[1]。相变表面的润湿性对相变过程有着重要影响，随着微纳米表面强化沸腾技术相关研究的深入进行，表面亲疏水性对沸腾影响的研究逐步展开[2]。郑晓欢等[3]对常压下光表面、超亲水和超疏水表面的池沸腾进行实验研究，认为超亲水表面亲水疏气，产生气泡体积小但气泡脱离频率较快；而超疏水表面亲气疏水，使得激活活化点更为容易。马强等[4]发现换热系数与表面过热度有关，当过热度较高，亲水和超亲水表面换热系数大于疏水和超疏水表面；但当过热度较低时，疏水和超疏水表面换热系数比亲水和超亲水表面大，这是由于疏水和超疏水表面有着较低的沸腾起点。在微纳尺度条件下，由于尺寸限制实验难以开展，利用基于经典力学原理的分子动力学（molecular dynamics，MD）模拟可以深入了解原子的微观动力学行为，将有助于理解在纳米尺度上相变的机制。已有学者对润湿性的影响进行了分子动力学模拟研究。Nagayama等[5]利用分子动力学研究了受限纳米通道中的气泡行为，发现亲水表面为均相成核，疏水表面却为异相成核。郑文秀等[6]认为润湿性越弱时，靠近壁面的分子向其他分子的能量传递越容易。相反的，张龙艳等[7]通过研究壁面浸润性对气泡初始核化和生长的影响，认为壁面浸润性越强，能量传递越快，气泡在固壁处越容易核化，气泡生长也更为迅速。Peng等[8]认为界面强度在形成吸收层中具有重要作用，蒸发液膜的热传递强烈依赖于吸收层。诸多学者对于纳米尺寸下的相变过程进行了研究。Chen等[9]通过模拟在具有非均匀润湿性纳米结构表面上的气泡成核行为，发现较大面积的强亲水区有利于气泡核的形成和生长。Hens等[10]认为亲水表面为气泡成核和蒸汽薄膜的形成提供了有利条件。Mao 等[11]研究了在有限空间内由热铜板加热的液态水膜的快速沸腾，发现靠近板的液态水分子立即过热并快速相变。Seyf等[12]研究了壁面材料（铝和银）的种类对超薄液态氩膜在纳米结构上爆炸沸腾的影响，发现由于固液相互作用的降低，固体银壁附近的原子的速度增加更多，导致液滴的上升速度增加，液滴与表面的分离发生得更快。与之不同，Hasan 等[13]在亲水表面观察到爆炸沸腾，但是对于疏水表面，观察到扩散蒸发，且亲水表面传热高于疏水表面的传热。

综上所述，壁面润湿性对固壁表面上液膜的沸腾过程有着重要影响。尽管对于纳米尺度下相变的研究已经取得了诸多成果，但到目前为止，不同学者对壁面润湿性影响传热过程机理并未达成一致，对于由于壁面润湿性不同带来的运动状态不同的研究涉及较少。因此，需要进一步研究纳米级相变传热和沸腾现象的行为和机理。本文利用分子动力学对不同润湿表面液膜沸腾过程进行模拟，研究不同程度亲疏水性壁面对沸腾的影响，重点分析不同润湿壁面沸腾过程中能量转换和液体运动中的差异。分析润湿性对其影响，有助于更好地理解纳米尺度下沸腾与核化的机理，以期望能够填补实际应用理论支撑的缺陷。本模拟采用开源分子动力学模拟软件LAMMPS 实现，原子位型实现可视化采用OVITO软件。

1 物理模型及模拟细节

建立图1所示模拟体系，对固体壁面及置于系统内氩原子进行分子动力学模拟。该模拟体系为Lx×Ly×Lz=34.57σ×227.15σ×6.91σ 的长方体盒子（σ为氩原子尺寸参数），由固体区、液体膜和蒸气区域组成。系统底部为固体几何模型，4680 个铂原子置于模拟系统的底部，以面心立方结构（FCC）排列方式构建金属固体，密度为ρσ3=2.613。固体几何模型从下至上分别设置为3部分：将最底面的一层原子固定，防止模拟过程基板发生变形；中间两层固体原子设置为“虚拟原子”，与外界恒温加热器相连作为热源；上部金属原子是真正的铂原子，与氩原子进行热量交换。对除固定原子外的固体铂原子施加弹簧力，将其约束在初始位置，弹簧系数为3249.1εσ-2。选择氩原子作为液体原子，7380 个液态氩原子被放置于固体铂表面，密度为ρσ3=0.78，厚度为34.5σ。模拟系统其余部分充满蒸气氩原子。x 和z 方向上应用周期性边界条件，在y方向应用固定边界条件。模拟框y+方向顶部设置反射墙，当氩原子试图穿过反射墙时，它将被反弹，该过程不会产生能量交换和损失。

图1 物理模型图

原子对之间的相互作用，本文采用Lennard-Jones（LJ）势函数(1)描述。

式中，下角标i、j可取为l或s，l和s分别代表流体原子和固体原子；ε为能量参数，取值分别为εll=1（对应的实际能量参数为1.6×10-21J）和εss=50εll（对应的实际能量参数为8.36×10-20J）；σ是尺寸参数， σll=1 （对应的实际尺寸参数为0.3405nm），σss=0.73σll（对应的实际尺寸参数为0.2475nm）。上述参数均来自文献[7]的报告。r是两个原子之间的距离，取截断半径为3.5σAr-Ar，超过截断半径时忽略其原子间作用力。

根据Lorentz-Berthelot 公式计算能量参数εls、尺寸参数σls，计算公式如式(2)、式(3)。

为获得不同润湿性的壁面，使用参数α对固液间势能进行调节，见式(4)。

表1 列出了模拟过程中不同粒子对之间的L-J参数和对应的接触角，本文中Case ε7.07表示固液作用为7.07的表面，其余两种类推。为了确定不同表面的润湿能力，如图1(b)和图1(c)所示，模拟了6nm 直径氩液滴在不同的润湿性基底上的扩散过程，以观察接触角。模拟过程在NVT系综下进行，温度控制为90K，模拟时间为2ns，以确保液滴充分在表面上蔓延。

采用velocity-Verlet 算法来求解整个粒子运动方程，时间步长为0.001τ，其中τ=(mσ2/ε)1/2为LJ约化时间常数。整个模拟过程分为3步：首先，使用共轭梯度（CG）算法实现能量最小化，以确保模型的初始状态是稳定的。采用正则系综（NVT）控制系统在T=0.744εκB-1（相对应实际温度为90K）下运行50 万步，在模拟过程中监测整个系统的温度、压力和能量变化，确保系统在NVT 过程结束时达到平衡状态。然后撤掉整体热浴，仅对基板进行控温，将基板的温度提高到T=1.57εκB-1（相对应实际温度为190K），此温度能够保证液体沸腾的同时达不到固体熔点。当将固体壁面加热至所研究的状态后，将T=0.744εκB-1的氩与T=1.57εκB-1高温壁面在微正则系综（NVE）进行集合并采集所需数据，在此时从0时刻开始重新计时。

表1 不同润湿性表面参数

2 模拟结果分析

2.1 液膜沸腾过程密度变化

疏水性和亲水性通常与液体为水时的润湿性现象有关，本文中虽然液体不是水，但遵循相同的惯例，仍用疏水性和亲水性表示，它们实际上意味着疏液性和亲液性。为明确显示液体薄膜在加热过程中的动态行为，图2给出了不同润湿表面上模拟系统运动过程快照。与前人的工作结果[10-13]一致：在最初阶段，液态氩在固体铂板表面被加热；随着加热的进行，氩原子开始蒸发；在加热进行一定程度时，近固体表面处的液体开始成核出现气泡并逐渐扩展成密度较小的气体层［图2(a)］，而距壁面较远的液氩分子仍为液体状态，近壁面处气体层推动未蒸发的液氩上升，出现固液分离现象。在液体区域呈现出3个部分：近壁面气体层，其上的液体区域以及系统最上面的蒸气区域。ε=7.07的表面，在150τ 左右时出现了分离现象，液体区域向上运动至顶部与反射墙碰撞被弹回，由于近固体表面蒸气产生的压力导致液体薄膜重新向上运动，在运动过程中液体层几乎保持不变，而ε=0.50 表面液体层在上升过程出现了明显缩小过程。随着润湿性的减弱，固液分离的发生时间延后，ε=1.00 表面与ε=0.50表面分别在900τ～1000τ和4800τ～4900τ之间出现了固液分离现象。基于本文与文献[10]在模型以及加热过程中出现的现象上的相似性，故认为本文结果是可靠的。

图2 不同润湿性壁面上沸腾过程

为了更进一步了解氩在y方向上运动情况和密度分布情况，将模拟体系沿着y 方向平均分成200个切片，计算在每个切片的原子数以获得每个切片的平均数密度。不同时刻y方向的数密度变化显示在图3中，纵坐标为原子的量纲为1数密度ρσ3，横坐标为y方向上的层数。随着加热时间的增加，近壁面处形成了低密度区，距离基板较远处氩密度增加。系统中液体氩的密度峰值逐渐减小，体积也逐渐缩小，缩小速率随着润湿性减弱而增大。Case ε7.07 液膜体积缩小速率极为缓慢，液体厚度和密度峰值基本不变，Case ε0.50 峰值密度降低速率更明显，意味着液体蒸发更快。Case ε7.07 在1300τ时数密度大于900τ，是因为液体运动到顶部受到挤压使得密度增大。同时应当指出，初始阶段壁面附近液氩数密度略高于液体整体密度，这是由于固液间相对作用力较强时，壁面附近会吸附氩原子形成“类固体”[14-15]。

图3 不同模型沿y方向氩数密度分布

图4 不同壁面润湿性下近壁面氩数密度分布

整个过程中壁面附近3.5σAr-Ar内氩原子数密度如图4所示，可以发现壁面附近氩原子数密度随着加热的进行呈显出先降低后稳定的趋势，润湿性强的壁面附近氩原子数密度下降较快，达到近似稳定状态也较早；随着润湿性减弱，下降速率也随之降低，达到近似稳定需要时间更长。从图2(a)可以看出，壁面附近始终存在氩原子，根据润湿性的不同，3 种表面上氩原子数密度分别稳定在0.453、0.136和0.102附近。随着亲水性的减弱，固液之间作用力降低，壁面附近氩原子数密度逐渐降低。

2.2 沸腾过程中能量变化

图5 不同壁面润湿性下流体温度变化

图5(a)给出了不同壁面润湿性下流体温度随时间变化规律，3种润湿表面氩温度均呈现先迅速增加后缓慢增加趋势，Case ε7.07 润湿性最强，温度快速上升最快，升温速率出现拐点时间最早，最终温度最低，随着润湿性减弱，升温速率随之减弱，温度上升变缓的时刻也随之延后，但最终的平均温度更高。这说明润湿性较强Case ε7.07在模拟开始传热更好，温度上升更快，发生固液分离较早；发生固液分离后气体层的存在阻碍了热量传递，使得温度上升速率减缓。随着润湿性减弱，初始阶段传热效果变差，温度上升减缓，发生固液分离时间延后。

图5(b)给出了加热开始120τ时不同壁面润湿性下流体温度沿厚度方向变化规律，此时气体层尚未出现。可以发现，靠近壁面处温度并不一致，润湿性最强的Case ε7.07近壁面温度也最高，润湿性降低，温度随之降低；液氩层远离固体壁一侧的温度近似一致。通过图5(b)可以看出，润湿性越强，壁面与流体间的温降越小，传热效果越好；润湿性强，液体中的温度梯度更大，近壁面气体更容易达到沸腾条件产生气体层，出现固液分离现象。

图6 不同壁面润湿性下流体能量变化

温度变化体现出氩原子吸收能量的差异。图6(a)给出了不同壁面润湿性下流体能量随时间变化规律，能量增量指每个时间步长上从固体基板流入流体的能量，通过统计所有氩原子的每原子势能和动能来获得能量增量。可以看出，能量增量与润湿性的关系和温度类似，先迅速增加，然后由于固液分离，能量交换变慢。此外，在模拟后期流体能量以较低的增长率保持增长，表明基板和液态氩之间仍有能量交换，这解释了为什么图5中温度能够继续上升。Case ε7.07 润湿性最强，能量上升最快，上升变缓时间最早，最终能量最低，随着润湿性减弱，能量增加速率随之减弱，但最终能量增量更大。结合图6(b)，可以得知吸收的能量主要转换为了势能。

不同润湿性壁面相变过程在时间上有很大差异，为探究差异的原因，计算了热通量的变化。通过氩原子总能量的变化除以液态氩薄膜的横截面积来获得热通量，即q=(1/A)×(∂EAr/∂T)，式中，EAr和A分别表示整个系统氩原子能量增量和系统y方向的截面面积。Cao 等[16]同样采用了这种方法计算了润湿性对传热的影响。不同壁面的热通量变化展示在图7中，不同表面的平均热通量在0τ达到最大值，然后随时间逐渐减小。这是因为初始时刻固液之间温差是最大的，随着传热过程的进行，固液温差减小，使得热通量减小，在150τ 左右的时候，ε=7.07表面热通量开始小于ε=1.00表面热通量，在约1000τ 的时候，ε=1.00 表面热通量开始小于ε=0.50表面热通量。初始阶段亲水表面的平均热通量明显高于疏水表面，这说明固液间相互作用力越强，传热速率越快。这种传热速率上的差异导致了润湿性强的表面发生爆炸沸腾的时间要早于润湿性弱的表面。

图7 不同壁面润湿性下热通量随时间变化规律

2.3 沸腾过程中液膜运动

氩原子平均质心位置随时间变化趋势展示在图8中，以便更为全面地了解整个过程运动情况。在ε=7.07 表面，在0～150τ 时间内质心在y 方向以一个较缓慢速度上升，在150τ 左右上升速度突然增加，直至上升到215σ，此时液氩层已到达模拟体系顶部。这与图3 中液体氩原子先缓慢受热膨胀、然后突然爆炸沸腾的现象具有一致性。由于顶部反射墙的存在，液氩原子运动方向在碰撞后发生了改变，开始向下运动，在向底部运动一定距离后，再次向上运动，做往返运动。对于ε=1.00 表面，氩原子质心上升到1580σ开始第1次下降，此时液氩层尚未达到顶部，而ε=0.50 表面上，在计算时间内并未发现质心有明显下降。此外，对于不同润湿性ε=7.07 与ε=0.50，质心能够到达的位置分别为215σ 和115σ。液体质心随时间变化的斜率可反映其上升速度，润湿性强的ε=7.07 表面上升速度要高于润湿性较弱的表面。

图8 不同壁面润湿性下氩质心变化情况

图9给出了整个过程上部气体层和下部气体层压力随时间的变化情况，P 和PL 分别用来指代液体层上方和下方的气体的压力。由于气体层出现时间不同，下层压力起始点并非0τ。可以看出，对于ε=7.07 表面，其下部压力远大于且一直大于上部气体压力，这可以解释液氩的迅速上升行为以及碰撞后液体层能够再次上升的原因；随着润湿性减小，在初始上升阶段的压差也随之减小。

ε=7.07 和ε=1.00 两种表面从固液分离到液氩首次改变运动方向时间段内压力变化如图10所示，ε=7.07表面下部压力要远大于上部压力，随着液体向上移动，下部压力在一定时间内呈下降趋势，上部压力增大；ε=1.00 表面上压力出现了交叉，在1200τ～1600τ 时间区间内，下部气体压力要低于上部气体，推测是由于蒸发使得上部压力升高而初始阶段下部压力较小。这种不同可以用于解释ε=7.07表面上的液氩有着很高的上升速度且到达顶部发生碰撞才改变方向，而ε=1.00 表面首次改变方向时并未到达顶部。

图9 不同壁面润湿性下压力变化情况

图10 不同壁面润湿性下流体压力

3 结论

通过分子动力学模拟研究固体基板上氩的相变过程，分析了壁面润湿性对沸腾过程中能量传递和液体运动情况的影响，得出了以下结论。

（1）在近壁面附近存在未蒸发液体原子，其数密度随着润湿性增加而增加。3种表面均发生了固液分离现象，出现固液分离的时间随润湿性增加而减小。固液分离后液体氩的密度、体积逐渐减小，减小速率随润湿性增加而降低。

（2）润湿性影响基板与液体间传热过程。润湿性强，能量交换速率快，热通量大，升温迅速且液体中温度梯度大。气膜的产生会阻碍传热，模拟结束时壁面润湿性越强，系统内流体温度越低。氩原子吸收的能量主要转换为势能。

（3）出现固液分离现象后，下部气体压力大于上部气体压力，推动液体层上升。润湿性越强，上下压差越大。在ε=1.00 表面，会在一段时间发现上部气体压力较大的情况。液体质心能够达到高度随着润湿性增强而增加。