伪势格子Boltzmann方法的汽化热阱效应模拟研究

金志浩，章 港，战洪仁，李 帅，王立鹏

(沈阳化工大学 机械与动力工程学院， 辽宁 沈阳 110142)

沸腾传热具有低温差、高热流密度的特性，被广泛应用于冶金、芯片制造等领域.强化传热有着巨大的经济效益，因此对沸腾传热机理的研究已经成为国际上研究热点.近些年国内外专家先后提出微液膜蒸发模型[1-2]、瞬态导热模型[3]、三相线模型[4]，但是到目前为止，对于沸腾强化传热机理没有统一解释.

1999年，沈自求在分析异气冲击加热壁面[5]时，提出汽化界面热阱效应，即异种气体冲击加热壁面,使气液界面上的液相快速汽化，吸取大量的汽化潜热，导致气液界面温度降低，并且壁面温度也随着汽化热阱效应的形成而降低.该效应很好地解释了沸腾传热增强原理.沸腾时气泡核与加热壁面之间存在液体微层，汽化热阱效应形成，其极大地加强了液体微层的导热，为我们研究沸腾传热机理提供新的思路.通过实验很难观测和测量到汽化热阱现象的微观过程，数值模拟可以有效记录与检测流场的速度和温度变化.

格子Boltzmann技术是一种介观模拟方法，该方法在处理非线性问题、复杂界面问题、相变问题有着先天优势，因此广受科研工作者青睐.1993年，Shan等首次提出伪势模型，即把作用力代入伪速度势中，该方法可以促使两相自发分离，使伪势模型有着极高的稳定性[6].2011年，Zhang用伪势格子Boltzmann模型与热力学方程相耦合，形成热格子Boltzmann模型，首次模拟了相变过程，但是控制流体状态方程选用范德瓦尔方程，导致模拟结果数值稳定性差[7].2018年，Gong提出把P-R状态方程代入能量方程，通过四阶龙格库塔格式求解能量方程，得出演化后的温度场，该方法数值稳定性强，模拟结果的热力学性质更接近真实流体[8].

为了避免多气泡相互作用而产生的影响，本文选用单气泡模型来解释沸腾传热机理.运用伪势格子Boltzmann方法，对单气泡成核、生长、脱离、重新生成，以及周围液体微层的流体细观流动与传热特性进行了模拟研究，首次结合界面汽化热阱效应对单气泡生长及其传热规律进行分析，结合热流密度变化曲线，讨论其强化传热机理，为研究相变、气泡动力学等行为对沸腾传热特性的影响提供了依据.

1 伪势格子Boltzmann模型

单组分伪势格子Boltzmann模型中流体粒子分布函数的演化方程为[9]

fi(x+eiδt,t+δt)-fi(x,t)=

(1)

(2)

(3)

密度ρ与粒子速度u计算公式如下:

(4)

宏观速度U为

(5)

式中:Δu=Fδt/ρ为速度变量;F为流体受到的合力，其由流体粒子间两相作用力Fint、流体粒子与固体壁面作用力Fs、重力Fg组成，合力F=Fint+Fs+Fg.

引入流体间粒子作用力Fint是为了促使气液两相分离.

(6)

式中[10]：β为权重因子，取值1.16，作用是提高数值稳定性；G为气液两相相互作用力常数，通常取值-1；x′为x的相邻格子点坐标；伪势函数ψ(x)的选取与流体的状态方程有关.

(7)

其中：p为流体压力；g为流体作用系数；c0为格子结构确定系数.

(8)

其中：sw用来判断格点是否是壁面的函数，对于固体壁面取1，对于流体壁面取0；Gw是流体与固体壁面的作用常数，调节Gw可以获得不同接触角壁面.

重力Fg表达式如下:

Fg=g(ρ-ρave).

(9)

式中g为重力加速度.

以上为格子Boltzmann单组分伪势模型，模拟相变过程采用以下温度场的能量方程[11-12].

(10)

式中:λ为导热系数；cv为定压比热；PEOS为流体压力;∇为哈密顿算子.

导热系数为

1.3.1 面积 种公羊面积需求为每只4m2以上，妊娠母羊面积需求为每只1m2以上，空怀母羊面积需求为每只0.8m2以上，育成羊面积需求为每只0.6m2以上。

λ=0.028ρ.

(11)

压力PEOS采用控制流体状态的P-R状态方程来得到.

(12)

然后对温度场应用二阶龙格库塔法进行离散.

Tt+δt=Tt+(h1+2h2).

(13)

其中:Tt为当前时间步长t下的温度场;Tt+δt为下一个时间步长的温度场.

2 模拟与结果

2.1 计算域与物理模型

计算区域如图1所示.计算空间内充满饱和温度为Tf=0.86Tc的气液两相流体，饱和液体密度ρl取值为6.5，饱和气体密度ρv取值为0.38，液相与气相导热系数比为17，无量纲重力加速度g取为-0.000 05.经过网格独立性验证，计算域200×250可保证数值的稳定性.设置加热长度为3个格子长度的恒温热源，其位置在下表面中心处，热源温度Ts=1.3Tc，热源附近壁面设为绝热边界条件[13-14]，容器左右侧壁面取周期边界条件[15-16]，上下固体壁面取Zou-He边界条件[17].

图1 物理模型

为验证程序正确性，对不同重力加速度下气泡脱离进行模拟，此前有数值模拟表示气泡脱离直径D与g-0.5成正比.在程序验证中，取无量纲重力加速度分别为-0.000 05、-0.000 04、-0.000 03、-0.000 02、-0.000 01，得到气泡脱离直径与无量纲重力加速度的关系，如图2所示，横坐标g为无量纲重力加速度，纵坐标D为气泡脱离直径.由图2计算可得，气泡脱离直径与无量纲重力加速度的拟合函数关系为D～g-0.5，符合Fritz理论.

图2 无量纲重力加速度g与脱离直径D的关系曲线

模拟中引入了特征长度l0、特征速度u0和特征时间t0.

(14)

(15)

t0=l0/u0.

(16)

其中:l0为毛细长度，是描述气液两相流的典型特征尺度，表示表面张力与浮力之比，根据拉普拉斯毛细方程获得对应饱和温度下的表面张力σ.

无量纲长度和无量纲时间可定义为：

l*=l/l0,

(17)

t*=t/t0.

(18)

式中:l为格子长度；t为格子时间.

2.2 亲水壁面单气泡生长过程

通过改变公式(8)中Gw的数值，得到接触角为60°的亲水壁面，以其模拟单个气泡成核、生长、脱离过程，如图3所示.图3(a)中液体开始相变形成气泡核，但是由于亲水壁面的亲水疏气性,相变的位置并不是紧贴加热面,而是在加热面正上方的饱和液体中进行，因此加热面与气泡核之间存在一层薄薄的液体微层.从图3(a)至图3(c)，气泡不断生长，期间主要依靠液层微层的相变为气泡长大提供气相，至图3(b)时，气泡下方的气液界面开始与壁面接触，至图3(c)时液体微层全部汽化，随后有流体补充至加热面，气泡继续生长[图3(d)].在图3(c)至图3(d)过程中，气泡在壁面张力的束缚下由圆形逐渐被拉伸成椭圆形，随后完全脱离.图3(e)为气泡脱离后的形态图.

图3 亲水壁面单气泡生长、脱离形态

图4模拟了亲水壁面单个气泡成核、生长、脱离过程的温度场.在气泡成核阶段，随着热源不断地加热，热量首先通过导热传递给附近的流体，由图4(a)可以看出，此时热源附近液体微层的温度大于周围饱和液体的温度，处于过热状态.对比图4中的(a)、(b)、(c)这3个阶段可以看出：近壁处的流体温度呈下降趋势，这是气泡的生长阶段，因为这个状态下的气泡生长主要依靠液体微层的相变，而相变的同时会带走大量的潜热，所以导致近壁处的流体温度下降.

图4 亲水壁面单气泡生长、脱离温度场

从图4(c)、图4(d)可以看出：热源附近的流体温度逐渐升高，这是由于气泡脱离加热面，液体回流至加热壁面，由于液相的热阻远小于气相，所以，热量迅速由热源传递至其附近流体中，使得温度升高.

2.3 疏水壁面单气泡的生长过程

通过改变公式(8)中Gw的数值，得到接触角为120°的疏水壁面，以其模拟单个气泡成核、生长、脱离过程，如图5所示.从如图5(a)中可以看出：初始阶段，由于热源对流体不断加热，此时靠近热源的液体开始发生相变，由于疏水壁面的亲气疏水性，相变是在紧贴壁面处进行.在图5(a)至图5(b)过程中，因为疏水壁面表面张力作用，气泡开始横向生长，三相接触线逐渐变长，此时主要依靠气泡周围的流体相变为气泡生长提供气相.在图5(b)至图5(c)这个阶段，随着气泡的生长其体积不断增大，导致其所受的浮力不断增大，气泡开始垂直向上移动，三相接触线开始变短.从图5(c)至图5(d)过程中可以看出：在壁面张力与气泡浮力的拉扯下，气泡在中间形成颈部并开始脱离.从图5(d)至图5(e)，气泡完成脱离，同时在壁面残留一部分气泡，作为下一个气泡生成的气泡胚核. 对比亲水壁面，疏水壁面气泡生长周期远大于亲水壁面气泡生长周期.

图5 疏水壁面气泡生长形态

图6模拟了疏水壁面单个气泡成核、生长、脱离过程的温度场.从图6(a)中可以看出,此时热源两侧近壁处的温度降低.这是由于疏水壁面的相变紧贴壁面进行，在壁面形成三相接触线，而此时气泡的长大主要依靠三相线处的液体相变，并且这时的三相线刚刚形成，距离热源比较近，所以,导致热源两侧近壁处的温度降低.从图6(b)中可以看出,此时热源附近的流体温度逐渐升高.这是由于随着气泡的生长,三相线逐渐远离热源，此时热源附近被气体覆盖,不会发生相变，所以，随着加热过程的进行，热源附近的流体温度逐渐升高.图6(c)至图6(d)，随着气泡的上移，三相线开始缩短逐渐靠近热源，导致热源附近的流体温度再次降低.

图6 疏水壁面气泡温度场

从模拟的结果中发现：加热亲水表面至其相变所需要时间为t*=1.8，同样加热条件下，加热疏水表面至其相变所需要时间为t*=1.2.由此可见，相同加热温度下，相较于亲水壁面，疏水壁面率先汽化.这是由于疏水表面具有亲气疏水特性，导致流体不会在加热面附近形成过热液体微层，而是直接相变；亲水表面由于其亲水疏气特性，热源处的流体紧贴壁面，在加热初期，相比于疏水壁面这部分流体很难发生相变，从而形成过热液体层，随着加热过程的进行，这部分流体积累能量逐渐增大，导致其在内部发生相变.Gao[18]在实验中，用激光干涉法测得疏水壁面相变不存在微液层，亲水壁面相变存在微液层.模拟结果与实验结果一致.

2.4 亲水壁面气泡生长时气泡气液界面对热源附近流体温度的影响

图7为热源与气泡气液界面温差ΔT和时间t*关系图.A—B段为气泡成核阶段，由于过热液体微层的相变带走大量潜热,导致气液界面处的流体温度下降，使得热源与气液界面的温差逐渐增大.B—C段为气泡纵向生长阶段(此时主要依靠液体微层的相变维持气泡生长)，随着气泡的生长,气液界面逐渐与壁面接触,导致其换热面积逐渐减小，所以,这一阶段热源与汽化界面的温差增长斜率小于A—B段.C—D段为气泡在加热面横向生长阶段(此时主要依靠三相线处的液体相变维持气泡生长)，从C点开始,近壁处的气泡气液界面全部消失,此时热源附近的流体全为气体,即C—D段不会发生相变,所以,热源会直接加热其附近的气体,使其温度升高,从而导致热源与附近流体的温差减小.

图7 热源与气液界面温差ΔT与无量纲时间t*的变化关系

2.5 亲水壁面气泡生长时气泡气液界面对热流密度曲线的影响

为了研究气泡气液界面对沸腾传热的影响，需要计算热源与其附近流体的热流密度.气泡气液界面主要依靠热源的导热传递热量，所以应用傅里叶导热定律进行计算，具体公式如下：

q(x)=-λ(∂T/∂y)|y=0.

图8为近壁处热流密度q与时间t*关系图，A、B、C、D与图7中的4个点相对应.由图7可知A—B流体的温度升高，导致传热温差降低，因此，热流密度曲线在A—B段缓慢下降.B—C段的热流密度包括两部分，即导热和热阱效应，热阱效应起主导作用，由图7可知在B—C段温差逐渐升高，从而导致液体转化为气体的量也在增加，而随后由于气泡气液界面的高度降低，液体逐渐减少，导致相变释放的热量也随之减少，所以，热流密度在B—C段出现最大极值，即气泡气液界面紧贴壁面时，热流密度大幅度上升到C点.C—D段进入脱离阶段，气泡逐渐脱离壁面，气泡气液界面高度逐渐升高，源源不断冷流体补充到壁面，相变减小，因此，C—D段热流密度大幅度减小.

图8 热流密度q与无量纲时间t*的关系

3 结 论

使用混合热格子Boltzmann方法，对亲疏水壁面单气泡成核、生长、脱离进行模拟，结合气液界面汽化热阱效应对传热特性的影响进行分析，得到以下结论：

(1) 相同加热条件下，相比于亲水壁面，疏水壁面率先汽化.这是由于疏水壁面具有亲气疏水性，致使其不会形成过热液体微层，直接发生相变；而亲水壁面由于亲水疏气特性，会先形成过热液体微层，而后发生相变.

(2) 疏水壁面气泡平铺壁面生长，热源附近不存在液体微层，所以，疏水壁面不存在汽化热阱现象；亲水壁面气泡在壁面生长，存在液体微层，因此，亲水壁面存在汽化热阱现象.

(3) 亲水壁面气泡生长脱离周期较疏水壁面短.亲水壁面气泡生长不同阶段，热源与气液界面间热流密度不相同：初始阶段，热源与附近流体温差逐渐降低，因此，热流密度逐渐减小；气泡生长阶段，由于液相相变产生汽化热阱效应，带走大量汽化潜热，大幅度提升热流密度，强化传热;脱离阶段，气泡逐渐脱离壁面，热流密度大幅度减小.