受限纳米通道内气泡核化的分子动力学模拟

2020-06-09苗瑞灿张石重陈占秀杨历

化工进展 2020年5期

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

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

近年来，微纳米科学技术在不同领域得到广泛应用，由于工质核化相变可以有效提高流体的换热特性，将工质核化相变原理应用到微纳尺度中已经成为广大学者研究的热点问题。流体在相变换热时会在通道内生成气泡并形成两相流。因此，研究微尺度下的流体相变及气泡成核生长过程及其机理是提高流体强化换热的一种重要方法。

固体壁面润湿性对气泡核化生长有重要影响，由经典核化理论可知，液体气泡核化沸腾过程与固液接触面粗糙度、液体过热度、固液间的接触角有一定的关系。固液壁面接触角增加，气泡核化沸腾所需要的过热度越低并且气泡发生沸腾所需时间越短。Bourdon等[1]利用实验对壁面润湿性影响大尺度沸腾初始核化情况进行研究，发现固体与壁面间润湿性越强，液体初始核化温度越高。Xu等[2]通过范德瓦尔斯理论研究沸腾核化在异质壁面形成过程，得到在异质情形下壁面润湿性越弱，气泡核化半径越大，更容易形成膜态沸腾。Jo 等[3]研究发现在壁面浸润性一致的情况下，异质核化相较于均质核化更容易发生沸腾现象。在常规尺度下，以流体为连续介质作为理论依据来研究固液之间传热传质机理，而在微纳尺度下流体以微观粒子为依据来研究固液之间的传热传质过程，固液间界面效应明显，导致微尺度下流体相变过程较常规尺度下不同。在微尺度下，已有学者对核化沸腾问题进行研究。Mao 等[4]利用分子动力学研究气泡在均质壁面核化过程，发现微尺度下的核化速率、极限过热度、核化速率与经典核化理论结果相比有较大差异。Kinjo等[5]研究有限尺寸的纳米通道内气泡核化生长过程，结果发现强亲水壁面形成均质核化过程，中等亲水壁面形成气泡异质核化过程以及疏水壁面形成气膜过程。Nagayama等[6]利用分子动力学研究壁面不同浸润性对气泡核化位置的影响，结果发现固液间作用力较强时，流体在通道内发生均质核化，气泡在通道中间的主流液体内生成，固液间相互作用力较弱时，流体在通道内发生异质核化，气泡在固体表面处形成，固液间作用力为超疏水时，在固体表面会形成一层气膜，而没有气泡生成。Hens等[7]发现超疏水壁面上很难形成气泡，而亲水性壁面由于较强的固液相互作用，在近壁面处容易形成类固体层，为气泡核化提供了有利条件。Yamamoto等[8]认为流体与固体之间相互作用力大小影响界面间能量与热量传递，壁面润湿性较强时，热量更容易由固体壁面传到流体中，更有利于核化现象的产生。She 等[9]利用分子动力学方法模拟壁面带有空腔的通道内气泡核化生长过程，发现与光滑壁面相比，带空腔壁面在近壁面处有较大的密度梯度、排斥力以及势能梯度，从而更有利于气泡核化生长。Chen等[10]利用分子动力学模拟壁面存在纳米凸起对气泡核化生长的影响，得到纳米凸起结构可以强化气泡核化，但存在最佳凸起高度。张龙艳等[11-12]利用分子动力学研究微纳通道内气泡核化生长过程，发现壁面润湿性较弱时，气泡在壁面凹槽内形成，壁面润湿性较强时，气泡在通道液体主流区形成。同时发现壁面润湿性越强越有利于气泡核化。此外，广大学者对微纳尺度下通道内流体传热传质过程及机理也有不同观点。郑文秀等[13]研究壁面不同润湿性对沸腾核化的影响，发现壁面润湿性越强，气泡形成位置离壁面越远，此时界面热阻较大不利于核化发生。Novak 等[14]利用分子动力学研究液体在固体壁面处的气泡核化过程，发现壁面润湿性越弱，气泡在壁面处发生异质核化的时间越短。

综上所述，微纳尺度界面效应对流体气泡核化生长及其相变过程有重要影响，学者们在微通道流体传热传质过程及机理分析中持有不同观点，且较少学者涉及受限通道内流体核化生长研究，而生活中大部分应用工程都是有限通道内流体传热传质相变过程，因此研究微尺度下受限通道内流体传热相变过程就显得尤为重要。本文利用分子动力学方法对受限纳米通道内流体气泡核化过程进行模拟，研究不同程度亲疏水性壁面对气泡核化生长的影响，分析固体壁面浸润性对气泡核化生长的作用机制。本模拟采用开源分子动力学模拟软件LAMMPS 实现，原子位型实现可视化采用VMD软件。

1 物理模型及模拟细节

图1 物理模型图

采用分子动力学研究了受限通道内流体核化生长过程，模拟系统如图1 所示，模拟体系尺寸为Lx×Ly×Lz=17.04nm×10.224nm×3.408nm，模拟体系在x、z方向上采用周期性边界条件，y方向上采用固定边界条件。固体壁面原子初始状态按照面心立方晶格（FCC）结构排列，固体壁面共有13714 个原子。固体壁厚h为1.704nm，流体高度H为6.816nm。壁面最外两层粒子固定不动，以维持模拟体系稳定，其余壁面粒子在初始位置做简谐振动。流体原子初始状态也按照照面心立方晶格（FCC）结构排列，系统初始温度设为85K，其流体饱和密度ρsat为1.401×103kg/m3。设置流体初始密度为0.7ρsat，通道内共有5760个流体原子。

流体原子与流体原子间相互作用力采用Lennard-Jones势能来描述，见式(1)。

式中，ε=0.01042eV；σ=0.3405nm；r为相互作用的两个原子之间的距离；rc为截断半径，文中取rc=0.85nm；流体原子质量M=39.948g/mol。固体与流体之间采用修正的L-J势能模型[15]，见式(2)。

固体原子间相互作用参数εs=0.4096eV，σs=0.23377nm，Ms=64g/mol，根据Lorentz-Berthelot公式计算能量参数εls，尺寸参数σls，计算如式(3)、(4)。

固体与固体原子之间相互作用力采用更加精确的EAM 嵌入原子势来进行模拟，文中使用的铜EAM 势能由Mishin 等[16]在2001 年提出，计算如式(5)、(6)。为了更好地观察计算统计结果，在模拟过程中将通道高度平均分成m层，第t层（1≤t≤m）在第(JF-JE)时间段内流体温度计算见式(7)。

式中，Fi为粒子i的嵌入能；φij(rij)为粒子i和粒子j间的对势能；ρi为第i原子处的电子云密度；rij为粒子i和粒子j间的距离；kB为玻尔兹曼常数，数值为1.380649×10-23J/K；m为单个流体质量；i、j为第i、j个流体粒子；a为通道内流体x、y、z三个方向，在系统达到平衡状态后uy、uz的值为0。

流体势能通过将通道高度平均分成n份，通过计算每层的平均势能，得到流体势能在通道内的整体分布，具体计算如式(8)。

式中，rij为粒子i、j之间的距离；uij为两粒子间的势能。

在模拟过程中对系统首先采用正则系综（NVT），维持系统温度为85K，运行50万步系统达到平衡，之后设置上壁面温度为80K，下壁面温度为140K，流体原子采用微正则系综（NVE），系统模拟运行150万步，观察通道内流体分布情况。模拟过程中系统采用Velocity-Verlet 算法求解运动方程，时间步长为1fs。

2 模拟结果分析

2.1 气泡生长核化过程

通过调节修正的L-J 势能函数中的α、β值来控制壁面润湿性，α、β的值越大，流体与固体间相互作用力越强，壁面润湿性越好。图2(a)～(e)为不同壁面润湿性下气泡核化生长过程。由图可见壁面润湿性不同，加热情况下流体在通道内有不同的表现形式。当固液势能参数α=1、β=1 时，流体与固体间作用力最强，流体被强作用力吸附在近壁面处，受热初期小气泡在液体薄层上生成，通道内发生均质核化，随受热时间增加，小气泡在液体薄层上方逐渐融合形成一层气膜，气膜推动通道内液体向上运动。由于通道上壁面温度为80K，且通道高度有限，系统平衡后气膜会在近热壁面处呈片条状，如图2(a)的2000ps快照图所示。当固液势能参数α=0.14、β=0.14 时，如图2(e)所示，近壁面处流体与固体间排斥力最强，此时壁面呈超疏水性，流体从受热初期到终了都不曾生成气泡。当固液势能参数α=0.6、β=1，α=0.14、β=1，加热初期流体在壁面液膜层生成较小的气泡，随时间增加小气泡逐渐融合形成气膜，推动流体向上移动，随着α的减小，气膜尺寸逐渐减小近壁面处液膜厚度也在减薄。当固液势能参数α=0.14、β=0.6 时，流体与固体间有较强排斥力，初期流体在固体壁面处生成小气泡，流体在通道内发生异质核化，随时间增加气泡尺寸逐渐增大，最终气泡大小基本保持不变。

图2 不同壁面润湿性下气泡核化生长过程

图3(a)、(b)中分别给出了势能参数为α=1、β=1，α=0.14、β=0.14 初始快照图，在近壁面处已将“类固体”层局部标出。“类固体”层是通道内流体由于壁面势能的相互作用力在近壁面固体层处流体粒子形成和固体壁面粒子类似的排列结构，这种由流体粒子形成的具有和壁面固体类似的结构层，成为“类固体”层。对比两个不同势能参数下的“类固体”层快照图，可以看出势能参数为α=1、β=1时，近壁面处流体排列方式更加整齐规则，结构近似趋于壁面粒子排列方式。势能参数为α=0.14、β=0.14时，在近壁面相同位置处流体排列方式与图2 相比，整齐度和规则度不高，“类固体”层形成程度较小。这种在近壁面处流体形成的规则结构程度对流体核化生长有一定的促进作用。

图3 不同壁面润湿性下气泡核化生长过程及初始快照图

图4是系统平衡后不同势能参数条件下发现流体质量密度沿通道高度Ly方向分布图，由图可见在近壁面处流体与固体壁面存在一层液体层，随固液势能参数值减小液体层的密度逐渐降低，近热壁面处液膜厚度也在逐渐减薄。在通道高度Ly方向上2.5～4nm 位置，除固液势能参数α=0.14、β=0.14外，流体质量密度均处于最低值，此处为气泡生成位置，由密度值可看出随壁面润湿性减弱气泡尺寸在逐渐减小。当固液势能参数为α=0.14、β=0.14时，流体与固体间排斥力强，近壁面处流体原子受较强排斥力作用，会在近壁面处形成一层流体原子极稀少的气膜层，因此会出现近壁面流体原子质量密度几乎为零的情况，由于固液间强排斥力作用流体原子均匀分布在通道3～7nm 位置处，很少存在于近壁面处。

图4 系统平衡后流体质量密度沿Ly方向分布

图5 是模拟时间为600ps 时流体质量密度沿Ly方向的分布图。图中给出了不同壁面润湿性下同一时刻的流体质量密度分布，由图可知系统模拟相同时间后固液势能参数α=1、β=1 时，在气泡核化位置处流体密度最小，固液势能参数α=0.14、β=0.6时，在气泡核化位置处流体密度最大，随着固液势能参数α、β增加，通道内气泡核化尺寸逐渐增大。固液间较大势能参数使流体在近壁面处有较强的流固作用力，壁面上吸附的流体粒子数越多质量密度越大，形成“类固体”层就越厚，热量更容易通过“类固体”层传给通道内流体，气泡越容易核化。因此，表面润湿性强的固体壁面更有利于流体发生气泡核化。

图5 时间为600ps时流体质量密度沿Ly方向的分布

图6是固液势能参数α=0.6、β=1不同时刻流体质量密度沿Ly方向的分布图。由图可知初始时刻通道内流体密度维持在1.0kg/m3左右，近壁面处流体密度相对较大，流体在近壁面处形成液膜薄层，随着壁面受热时间增加，液膜厚度逐渐减薄，在2.5～4nm 位置流体密度显著下降，流体在此位置发生均质核化生成气泡如图3(a)所示。此外，流体发生均质核化过程主要在受热初期阶段，500～730ps 是流体在通道内核化的主要阶段，气泡在此时间段内生成，之后随受热时间增加气泡尺寸不在发生变化，小气泡会逐渐融合形成气膜，气膜推动流体向上运动。

图6 势能参数α=0.6、β=1流体质量密度分布

2.2 气泡核化机理

图7 不同壁面润湿性下流体温度分布

通过分析不同势能参数下流体温度变化。近一步研究壁面润湿性对流体传热过程的影响。将系统模型沿Ly方向平均分成50 层，计算每层流体粒子的平均温度，只统计流体区域温度变化。图7给出了系统平衡后不同壁面润湿性下流体温度分布情况。由图可见，固液势能参数值越大近壁面处流体温度越高，即壁面润湿性越强，固液传热效率越好，因此，流体在润湿性强的壁面处获得的能量多，更容易达到核化条件形成气泡。势能参数α=1、β=1时，在Ly为2.8nm位置有温度突变，由于系统稳定后气泡融合形成气膜[如图2(a)2000ps时的快照图]，气膜层处含有极少数流体原子，温度从固体壁面传到气膜层再经气膜层传给通道上层流体，由于气膜层传热效率低，极少热量通过气膜层传给上层流体，因此温度会在有气膜层的地方出现突变。

图8 给出了势能参数为α=0.6、β=1，α=0.14、β=1不同时刻流体温度沿Ly方向变化。由图可知壁面温度经过固液接触面传给通道内流体，靠近热壁面处的流体温度率先升高，随通道高度Ly方向增加流体温度逐渐降低。此外，流体温升主要发生在1000ps 之前，之后随时间增加温度几乎保持不变，说明在受热初期流体在近热壁面处发生核化后，气泡核化尺寸逐渐趋于平稳，随加热时间增长温度也达到平衡。由图8(a)、(b)对比发现，在800ps 时，势能参数为α=0.6、β=1，近壁面温度已经接近140K，随时间增加近壁面处温度保持不变。势能参数为α=0.14、β=1的情况下，近壁面处温度随时间增加逐渐升高，最终近壁面处温度接近132K，说明在相同时间内，势能参数值越大，热量越容易传递给通道内流体，即润湿性强的壁面界面热阻小，固液间传热效率高，气泡核化时间短。

图8 不同势能参数、不同时刻流体温度沿Ly方向分布

图9为势能参数α=0.14、β=1时x、y、z方向平均受力分布图，由图可知，在x、z方向上流体粒子受力几乎为零，说明流体粒子在x、z方向没有优先运动，因为x、z方向设置为周期性边界条件。在近热壁面处观察到y方向上力先从0.0058下降到零，之后又从零增加到0.003，最后从0.003一直减小到零左右，力的数值为负代表力的作用方向为y轴负方向，在通道中心力的大小维持在零左右，近冷壁面处力从0 增加到0.003 左右。观察发现在近壁面处有力的梯度出现，在通道中心处几乎没有力的梯度，并且热壁面处力的梯度要大于冷壁面处力的梯度。近壁面处流体所受势能大小相同，但力的梯度却不同，说明高温导致流体在近壁面处表现出较大的力度梯度，从而流体更容易出现核化现象。图9 是不同壁面润湿性下y方向受力分布图，根据图9可知流体粒子在x、z方向上受力几乎为零，所以图10只给出不同壁面润湿性情况下y方向上的受力情况。由图可知，不同壁面润湿性下，流体粒子在近热壁面处力的梯度都要大于冷壁面附近力的梯度，在近热壁面处，不同润湿性壁面力的梯度虽然不同，但相差很小。固液势能参数越大，流体在近壁面处越容易出现“类固体”现象，更容易束缚在近壁面不易进行自由运动，因此近壁面处会出现随势能参数值下降，流体粒子在y方向受力增加的现象。

图9 势能参数为α=0.14、β=1时x、y、z方向受力分布

图10 不同壁面润湿性下y方向受力情况

图11 不同壁面润湿性下流体势能分布

图11 是不同壁面润湿性下流体在通道内势能分布情况，由图可以看出，流体在通道中心势能几乎都保持在0.05左右，在靠近壁面处，流体势能出现较大差异，当势能参数α=1、β=1 时，近壁面处流体势能最大，约为0.25。当势能参数α=0.14、β=0.6 时，近壁面处流体势能最小，几乎接近于零。近壁面处势能大小随α、β值的增加而增加，流体在近壁面处表现的势能强度导致壁面不同程度的润湿性，近壁面处流体势能强度越大，壁面原子对流体原子吸附力越强，壁面润湿性越强。由于壁面势能参数不同，流体在近壁面处表现出不同的现象，当势能参数最强时，流体会在壁面处形成一层液膜，流体排列方式出现“类固体”形式，壁面通过“类固体”将热量传递给通道内流体，传热效果较好。当势能参数较小时，在近壁面处由于壁面粒子与流体粒子间有较强排斥力，导致在近壁面处出现一层气膜，壁面热量透过气膜将热量传递给流体，传热效果相对较差。因此，当势能参数较小时，热量很难传递到流体中，通道内流体将不会成核形成气泡。