1. 引言
核态沸腾(Nucleate boiling heat transfer)是一种公认的高效换热方式。近些年,许多研究者采用各种方法提高沸腾换热效率,例如改变流体性质 [1],改变表面润湿度 [2],或者对表面进行改造 [3] 等。其中,表面改造方法较为简单,成本较低,因此得到越来越多的关注。
通常来说,表面改造是将原本平整光滑的表面,通过一系列加工方式,在表面形成空腔或者障碍物,这些空腔或者障碍物可以作为沸腾中气泡的成核位点 [4]。Kandlikar [5] 的实验结果表明,和平直表面相比,空腔中的传热系数提高了8倍,临界热流密度(Critical heat flux, CHF)提高了2.5倍。不平整的表面促进了表面液体和蒸汽的流动,更有利于加热表面被液体润湿,因此改造后表面有更高的换热性能。尽管表面改造对核态沸腾换热有明显增强,但是此方面研究尚有许多不足。由于表面结构复杂,难以定量研究表面特征对核态沸腾的影响。因此,许多研究者对表面进行规则化处理。Dong等人 [6] 在具有规则排列的圆形纳米空腔表面进行了沸腾实验,结果表明,减小空腔直径能使成核位点增多,但提高传热效果并不显著。类似地,Ma和Cheng [7] 对比了规则排列的微腔表面和微柱表面,发现微柱表面更容易产生毛细润湿现象,即加热表面更容易被液体润湿,因此具有更高的CHF。Shen等人 [8] 研究了规则排列的矩形柱表面的核态沸腾换热,发现柱高度对换热性能的影响大于柱宽度的影响。通过对已有文献的回顾,可以看出在研究规则排列表面的沸腾时,大多数文献考虑了矩形和圆柱形。然而,其他形状障碍物对核态沸腾换热也有明显影响,例如锥形障碍物。Liter和Kaviany [9] 在沸腾表面添加了锥形结构的障碍物,通过实验和理论分析,发现和平直表面相比,锥形结构障碍物表面液体流动阻力更小,即气泡容易从表面脱离并且壁面更容易被液体润湿,所以锥形结构表面有更好的换热性能和更高的CHF。然而,由于实验的局限性,他们并没有揭示气泡在锥形结构表面的成核,生长和脱离过程,一些锥形结构的参数例如结构间距,结构尺寸对核态沸腾的影响尚未得知。因此,本文对锥形结构表面的沸腾研究做进一步完善。
在实际的沸腾实验中,由于实验条件的限制,气泡的形成、生长和脱离很难被观察到,难以准确地分析浮力、壁面等因素对气泡和壁面传热的影响。因此,许多研究者使用数值模拟的方法来研究核态沸腾。近些年,格子玻尔兹曼方法(lattice Boltzmann method, LBM)被广泛用于沸腾研究 [10] [11] [12]。在LBM中,每个计算格点的状态(气态或液态)由所给的状态方程决定,因此不需要额外追踪相界面。本文选择的相变LBM模型,主要优点包括 [11]:1) 可以模拟整个沸腾过程,包括气泡的成核,2) 边界条件易于实现,3) 模拟中无需假设沸腾等待时间,4) 计算量较小。
在本文中,采用最近提出相变LB模型,在具有锥形结构的加热表面模拟了核态沸腾。具体讨论了包括锥形结构表面气泡的成核、生长及脱离过程,锥形结构间距对气泡生长和CHF的影响以及锥形结构尺寸对气泡生长和CHF的影响。
2. 数值方法与模型验证
2.1. 数值方法
本文采用了Gong和Cheng提出的相变LB模型 [11],在此模型中,密度分布函数的演化方程为:
(1)
其中
是格子时间,
是密度分布函数,
是离散速度,i是离散速度方向,
是外力项,
是松弛时间,
是平衡态分布函数,其表达式为:
(2)
在公式(2)中,ρ为流体密度,u为平衡态速度,cs为格子声速,其值为
,
是权系数,其取值与选取的离散模型有关,本文选取D2Q9 (二维九速)模型,其对应权系数为:当i = 0时,
;当i = 1~4时,
;当i = 5~8时,
。ei为离散速度,D2Q9的离散速度模型为:
(3)
其中c为格子速度取c = 1,在公式(1)中,外力项
表达式为:
(4)
其中
,
为合力,其由三部分组成:
(5)
其中
是重力或者浮力,
为液-液相互作用力。
表达式为:
(6)
其中G为相互作用力强度,
为“有效密度”,
,p为压力。在本文中,选取Redlich-Kwong (R-K)状态方程因为其有良好的精确性和简便性,其表达式为:
(7)
重力(浮力)
的表达式为:
(8)
其中,
是重力加速度,
是整个计算区域每个计算步的流体平均密度。
另一方面,温度分布函数的演化方程可以由以下公式给出:
(9)
其中,
是温度分布函数,
是温度松弛时间,
是温度平衡态分布函数,其表达式为:
(10)
其中T为温度,U为流体速度。在公式(10)中,源项
与相变有关,其表达式为:
(11)
其中cv为比热。宏观密度,速度和温度由以下公式给出:
(12)
2.2. 模型验证
根据前人的工作 [13],在水平光滑表面,气泡从表面脱离时,其脱离时的直径与重力大小有关,并满足
。接下来,将利用上述相变LB模型,模拟水平光滑表面单个气泡的生长脱离过程,对该经验关系进行验证。
初始时,整个计算区域充满饱和液体,其密度和温度分别为ρl和Tsat。计算区域底部为无滑移壁面,在底部中心有一个长度为5格子的热源,其温度恒定为Tw,底部其他区域为无滑移绝热边界。上表面是等温边界并且采用对流边界使流体可以自由流出 [14]。左右两侧为周期边界。在本文中,无滑移边界和等温(绝热)边界均采用半反弹处理格式 [15] [16]。计算区域大小为150 × 450,其他参数设置为:Tsat = 0.9Tc,a = 2/49,b = 2/21,R = 1.0,τ = 1.0,τT = 1.0,cv = 6.0,ρl = 5.426,Tw = 0.98Tc。
气泡脱离直径(气泡从壁面脱离时的当量直径)和重力关系如图1所示。图中各点代表了利用上述LB模型在不同重力下得到的脱离直径,红色曲线为拟合曲线。脱离直径与重力满足
,与实验关系吻合良好 [13]。
Figure 1. The relationship between bubble departure diameter and gravity
图1. 气泡脱离直径和重力关系
3. 初始条件及边界条件
如图2所示,计算区域大小为Lx × Ly = 490 × 1000,其中Lx为横向长度,Ly为纵向长度。初始时,计算区域内充满静止的饱和液体,温度为Tsat,密度为ρl。底部有一个表面被加工为锥形结构的加热器,表面锥形结构等间距排列,加热器底部采用恒温加热,温度为Tw,两端为绝热,加热器上表面及两侧均为无滑移边界。加热器总长为L,绝热边高为H,表面的锥形结构数量为N,宽为D,高为0.5 D,锥形结构的间距为P。计算区域底部无加热器部分为绝热边界,顶部为自由出口采用对流边界条件 [14],左右两边为周期边界。本文中,无滑移边界和等温(绝热)边界均采用半反弹处理格式 [15] [16]。若无其他声明,本文剩余部分参数设置均为:Tsat = 0.9Tc,ρl = 5.426,a = 2/49,b = 2/21,R = 1.0,τ = 1.0,τT = 1.0,cv = 6.0,L = 390,H = 10。
4. 结果与讨论
4.1. 锥形结构表面气泡的生长脱离
图3给出了气泡在锥形表面的形成,生长及脱离过程,其他参数为:DTw = 0.095Tc,D = 30,P = 54,其中DTw = Tw − Tsat。如图所示,在t = 6000时刻,气泡成核出现在锥形结构的两侧角点处,因为固体的热扩散率大于液体的热扩散率,这导致角点处侧壁的温度大于凹槽内表面液体的温度。因此,角点处的液体先发生了相变。t = 13,000时刻,可以看到最两侧的气泡先发生了脱离而且有向中间流动的趋势,因为加热器两侧的液体向加热器中心流动,从而对气泡产生向中心拖拽力 [17]。另一方面,此时凹槽内未脱离的气泡被浮力拉伸,凹槽内没有发生气泡合并。在t = 17,000时刻,中心的气泡从壁面脱离并且发生合并,两侧的气泡直接脱离没有发生合并。这与常见的矩形结构表面不同 [18],矩形结构的凹槽中气泡发生合并,合并后的气泡不会完全从壁面脱离,会堆积在凹槽中。和矩形结构相比,可以看出锥形结构表面的气泡更容易脱离,因此更有利与液体对加热表面的润湿,锥形结构表面具有更佳的传热效果,这和实验得到的结果相同 [19] [20]。
(a) t = 6000 (b) t = 13,000 (c) t = 17,000
Figure 3. Bubble growth and departure process on conical structure surface
图3. 锥形结构表面气泡生长脱离过程
为了进一步分析锥形结构表面气泡生长过程,图4给出了气泡在角点处生长时的受力情况。如图所示,水平加热壁面和倾斜加热壁面分别对气泡产生蒸发动量力,F1和F2,方向均为平行于壁面且向气泡外部。当气泡生长时,水平壁面温度高于倾斜壁面温度,即F1 > F2,所以,F2在水平方向上的分力小于F1,导致气泡在水平方向上受力不平衡,从而气泡在壁面发生移动。在矩形结构表面,F2在水平方向上的分力为0,气泡在F1的作用下更容易向凹槽内生长,因此气泡在矩形结构表面的凹槽内更容易发生合并 [18]。
Figure 4. Bubble evaporation momentum force Model domain
图4. 气泡蒸发动量力示意图
4.2. 锥形结构间距的影响
图5给出了锥形结构间距P = 110时的气泡生长脱离过程,其他参数为:Tw = 0.095 Tc,D = 30。在t = 4100时刻,障碍物角点和凹槽内液体都发生了相变,并且角点处的气泡较大,这是因为角点处温度较高。需要注意的是,由于障碍物间距较大,凹槽内的壁面能加热较多的液体,因此凹槽内也产生了相变。凹槽内的蒸汽将两边角点处的气泡连通,形成一个大气泡,如t = 6000时刻所示。随着气泡的生长,t = 9000时,气泡从壁面脱离。与此同时,在脱离的气泡两侧产生了新的气泡。与P = 54的情况相比,气泡脱离时刻发生了提前,因为表面障碍物数量减少,促进了表面液体的流动。可以看出,气泡总是在凹槽内(包括角点)内产生,锥形障碍物表面不会发生气泡成核,因此,凹槽内的换热效率大于障碍物表面。
(a) t = 41,000 (b) t = 6000 (c) t = 9000
Figure 5. Bubble growth and departure process at P = 110
图5. 锥形障碍物间距110时的气泡生长脱离过程
为了进一步研究锥形结构表面的沸腾现象,下面我们得到了锥形结构表面的沸腾曲线。其中为时间–空间平均热流密度
的表达式为:
(13)
其中t1和t2表示沸腾开始后的很长一段时间,
为空间平均热流密度,其表达式为:
(14)
其中角标s-l表示固液界面,l为固体的导热率。
图6给出了不同锥形结构间距的沸腾曲线,图中曲线的峰值为临界热流密度CHF,CHF之前为核态沸腾阶段。如图所示,在核态沸腾阶段,壁面的传热性能随着间距的增加而增加。从前文分析可以得知,障碍物间距增加,表面凹槽的总长度增加,有效换热区域增加,因此障碍物间距较长的表面有较好的换热性能。图中每条曲线的峰值即对应的CHF,可以看出,CHF随着间距的增加而增加,这是因为障碍物间距较大时,气泡更容易从表面脱离,液体更容易润湿表面,因此障碍物间距较大的表面有着较高的CHF。
Figure 6. Boiling curves for different distances
图6. 不同间距下的沸腾曲线
4.3. 锥形结构尺寸的影响
图7给出了锥形结构尺寸为D = 20时的气泡生长过程,其他参数为:Tw = 0.095 Tc,P = 62。t = 4100时,障碍物角点和凹槽内液体都发生了相变,随着沸腾过程的进行,角点处的气泡与凹槽内的气泡发生相互作用,合并成一个大气泡,如t = 6000时刻所示。最终,在t = 8000时刻气泡发生了脱离。与D = 30的情况相比,气泡脱离发生了提前,因为障碍物尺寸减小,液体更容易润湿表面,加快了气泡的脱离。
图8给出了不同锥形结构尺寸下的沸腾曲线,图中曲线的峰值为临界热流密度CHF,CHF之前为核态沸腾阶段。如图所示,在核态沸腾阶段,壁面的传热性能随着尺寸的增加而降低。从前文分析可以得知,障碍物尺寸增加,导致凹槽间距减小,相应的表面凹槽的总长度减小,有效换热区域减小,因此障碍物尺寸较大的表面换热性能较差。图中每条曲线的峰值即对应结构尺寸的CHF,可以看出,CHF随着结构尺寸的增加而减小,这是因为障碍物尺寸较小时,凹槽间距变大,气泡更容易从表面脱离,液体更容易润湿表面,因此障碍物尺寸较小的表面有着较高的CHF。
(a) t = 41,000 (b) t = 6000 (c) t = 9000
Figure 7. Bubble growth and departure process at D = 20
图7. 锥形障碍物尺寸20时的气泡生长脱离过程
Figure 8. Boiling curves for different sizes
图8. 不同尺寸下的沸腾曲线
5. 结果与讨论
本文采用相变LBM模型,在锥形结构表面模拟了气泡的成核、生长和脱离,并且得到了锥形结构表面的沸腾曲线,具体结论如下:
1) 在锥形结构表面,气泡成核发生在表面障碍物角点处和凹槽内。与矩形结构表面相比,锥形结构表面气泡脱离更容易,液体更容易润湿表面。
2) 随着锥形结构间距的增加,气泡脱离时间提前,核态沸腾阶段的传热效率增加,临界热流密度也增加。
3) 随着锥形结构尺寸的增加,气泡脱离时间变长,核态沸腾阶段的传热效率变差,临界热流密度减小。
基金项目
本研究由国家自然科学基金(51976128)及上海市自然科学委员会项目(19ZR1435700)资助。
NOTES
*通讯作者。