李宗坤,范广铭,*,曾晓波,严一鸣,马富兵

(1.哈尔滨工程大学 核科学与技术学院,黑龙江 哈尔滨 150001;2.哈尔滨工程大学 黑龙江省核动力装置性能与设备重点实验室,黑龙江 哈尔滨 150001)

由于相变的存在,过冷沸腾作为一种高效的传热方式被广泛应用,因此许多学者对过冷沸腾的传热计算开展了相关研究。随着传热学的发展,过冷沸腾传热计算已逐渐从纯经验模型过渡到半经验模型,再过渡到如今的机理模型,其中汽泡动力学和壁面热流分配模型在过冷沸腾热流密度计算中的应用受到广泛关注与研究。

在采取壁面热流分配机理模型计算壁面热流密度时,汽泡脱离直径、浮升直径等参数对壁面热流密度计算具有很大的影响。Klausner等[1]和Zeng等[2]的汽泡力平衡模型是目前应用较多的求解汽泡相关参数的模型。在最近的研究中,Akand等[3-4]分别对竖直加热壁面和以IVR(implementation of in-vessel retentio)为背景的倾斜加热壁面的汽泡受力分析开展了相关研究,并对汽泡浮升直径和临界热流密度(CHF)进行了计算,计算得到的汽泡浮升直径的平均误差约为9.89%,CHF的平均误差约为7.2%。基于汽泡动力学研究,王畅等[5]建立了窄矩形通道内的过冷流动沸腾传热系数计算关系式,计算误差可控制在±30%以内。Ren等[6-8]分别对汽化核心密度、窄矩形通道汽泡滑移、结合汽泡动力学的汽泡脱离直径和汽泡浮升直径等开展了研究。

对于壁面热流分配机理模型,经典RPI模型是广为接受的机理模型,该模型最初应用于池式沸腾的热流密度计算,后来进一步发展到过冷流动沸腾中。在RPI模型的基础上,Sateesh等[9]考虑了竖直加热壁面的汽泡滑移效应,Yeoh等[10]基于垂直环形通道中向上强迫对流过冷沸腾开发了壁面热流分配机理模型,闫美月[11]针对矩形通道内强迫对流过冷沸腾建立了壁面热流分配机理模型,以达到计算矩形通道内壁面热流密度的目的。

目前,壁面热流分配模型多数以汽泡力平衡模型为基础,并且在目前所公开的文献中,多数研究集中在垂直或水平平板加热工况。然而,实际情况中存在一些特殊的棒束通道,如倾斜安装的换热器、海洋核动力浮动平台或海洋条件下船舶核动力装置中的棒束通道,棒束长时间处于倾斜状态[12-13]。在倾斜条件下,由于汽泡受力方向的改变,倾斜棒束壁面汽泡可能会直接脱离加热壁面而不存在滑移现象,也可能出现汽泡绕加热棒曲面滑移的现象。然而,壁面热流分配机理模型很大程度上依赖于汽泡相关参数,在计算壁面热流密度的过程中,倾斜条件对汽泡行为以及汽泡参数会产生一定影响,进而对壁面热流分配模型的计算产生很大影响。因此,针对倾斜条件下棒束壁面的传热计算,有必要结合汽泡动力学和壁面热流分配模型开展壁面热流密度计算的进一步研究。本文主要进行倾斜条件下汽泡力平衡模型的建立与计算,在力平衡模型基础上建立壁面热流分配模型,并与实验数据进行对比验证。

1 汽泡力平衡模型

根据Klausner经典汽泡受力理论[1],对倾斜条件下棒束的上加热壁面进行汽泡受力分析,如图1所示,其中汽泡在x方向和y方向的受力方程可表示为:

图1 倾斜加热棒表面汽泡受力分析

∑Fx=Fs,x+Fsl+Fh+

Fcp+Fbsinθinc+Fdu,x

∑Fy=Fs,y+Fqs+Fbcosθinc+Fdu,y

(1)

式中:Fs为表面张力,N;Fsl为剪切升力,N;Fh为水动力压力,N;Fcp为接触压力,N;Fb为浮力,N;θinc为倾斜角度,(°);Fdu为非对称生长力,N;Fqs为流动曳力(准稳态曳力),N;下标x代表垂直于流向的方向,y代表平行于流动方向的方向。

金頔等[14]和肖仁杰等[15]也曾把倾斜角度的影响考虑在汽泡受力分析上,但对于汽泡脱离或浮升的条件,均是x方向或y方向的力打破平衡,即当∑Fx>0或∑Fy>0时,汽泡开始脱离或浮升。但倾斜条件的存在会导致近壁面汽泡行为发生改变。一方面,实验现象拍摄结果表明,在倾斜上加热壁面的汽泡,滑移距离非常短,并且多数汽泡会沿浮力方向(竖直向上)直接脱离加热壁面,如图2所示;另一方面,相比于汽泡尺寸,相邻加热棒壁面对近壁面汽泡行为影响较小,即在浮力引起汽泡脱离的方向上无限制汽泡脱离的其余壁面,且根据Zeng等[2]的研究,浮力在汽泡脱离过程中占据重要地位。因此,对于倾斜上加热壁面汽泡脱离判定条件定为浮力方向(竖直向上)受力总和大于0,即∑Fg>0时,汽泡脱离加热壁面,如式(2)所示:

图2 汽泡脱离示意图

∑Fg=Fb+Fqscosθinc+(Fsl+Fh+Fcp)·

Fs,xsinθinc+Fs,ycosθinc>0

(2)

式中,θi为汽泡倾斜角。

根据Zeng等[2]的研究,相比于其余汽泡受力的量级,水动力压力和接触压力的量级远小于其余汽泡受力量级。因此,在处理倾斜加热壁面汽泡受力时,水动力压力和接触压力忽略不计,即Fh、Fcp→0。简化后的方程如式(3)所示:

∑Fg=Fb+Fqscosθinc+

Fs,xsinθinc+Fs,ycosθinc>0

(3)

1) 浮力

浮力是汽泡在重力场中由于汽液密度差引起的力,计算公式如式(4)所示:

(4)

式中:Vb为汽泡体积,m3;rb为汽泡生长半径,m;ρl为液相密度,kg/m3;ρv为汽相密度,kg/m3;g为重力加速度,m/s2。

2) 非对称生长力

非对称生长力是汽泡由于非对称生长而导致的力,本文采用Klausner方法计算,如式(5)所示:

(5)

式中,r′b和r″b分别为汽泡生长速度和汽泡生长加速度。沿汽泡倾斜角(θi)方向分解为x、y方向的对称生长力,其中汽泡倾斜角根据Ren等[8]的研究取15°。

汽泡生长速度和汽泡生长加速度的计算需要汽泡尺寸和时间的相对关系,即汽泡生长控制方程。目前,Zuber[16]汽泡生长模型应用较广泛,因此,本文采用Zuber模型,如式(6)所示:

(6)

式中:b为考虑汽泡非球状的经验修正系数,Zeng的研究表明b的范围为1～1.73,而且b取1时,其实验汽泡浮升直径和脱离直径符合最好;Ja为雅可比数;α为热扩散率;t为汽泡生长时间,s;kl为液相热导率,W/(m·℃);cpl为比定压热容,kJ/(kg·℃);Tw为加热壁面温度,℃;Tsat为饱和水温度,℃;γ为汽化潜热,kJ/kg。根据现有研究以及最终实验结果,本文在计算倾斜工况时b取1。对R(t)分别求一阶导数和二阶导数即可获得汽泡生长速度和汽泡生长加速度。

3) 剪切升力

由于加热壁面法向速度梯度的存在而导致汽泡在加热壁面法向方向上受到剪切升力。本文采用Mei等[17]推导的剪切升力公式,如式(7)所示:

(7)

式中:Gs为剪切速率;Reb为汽泡雷诺数。

(8)

Reb=2ρ1urrb/μl

(9)

式中:ur为汽泡质心处液相速度与汽泡速度之差,在汽泡脱离前,认为汽泡速度为0,ur即为汽泡质心处液相速度;μl为动力黏度,Pa·s。本文采用Reichardt[18]提出的公式计算近壁面液相速度,如式(10)所示:

(10)

式中:u(x)为壁面法向距离x处的液相流速,m/s;x+为无量纲壁面距离;u*为摩擦速度;δ、c和χ为经验常数,分别取值0.4、7.4和11。计算时,取法向距离为汽泡半径时的近壁面液相速度。

x+和u*计算公式如式(11)所示:

(11)

式中:vl为运动黏度,m2/s;τw为壁面剪应力,kg/(m·s2)。

(12)

式中:u为主流速度,m/s;λ为摩擦系数,由Blasius公式计算。

λ=0.316 4Re-0.25

(13)

4) 流动曳力

由于流体与汽泡存在速度差而引起的阻碍汽泡在流体中发生相对运动所受到的力,称为流动曳力。采用Mei推导的公式进行计算,如式(14)所示:

Fqs=6πvlρlurrb·

(14)

5) 表面张力

根据Klausner的研究,表面张力计算公式如式(15)所示:

(15)

式中:dw为汽泡与壁面的接触直径,m;σ为表面张力系数,N/m;α和β分别为汽泡的前、后接触角,rad。根据Yun等[19]和Jia等[20]的研究,dw=0.134rb,α和β分别取45°和36°。

以时间为变量,利用Zuber汽泡生长模型控制汽泡尺寸,计算得到汽泡受力。根据倾斜工况下汽泡受力的不同与相关计算,依据判别标准判断汽泡是否脱离,进而得到汽泡脱离直径(dd)和汽泡生长时间(tg)等汽泡相关参数。

2 壁面热流分配机理模型

在过冷沸腾传热特性研究过程中,已开发出许多过冷沸腾传热计算模型,主要包含3种壁面热流密度的计算形式:纯经验关系式、壁面热流密度划分半经验关系式、壁面热流密度分配机理模型[21]。本文主要针对第3种机理模型给出倾斜加热棒束壁面热流分配模型。

对于倾斜条件下棒束的上加热壁面,只存有汽泡脱离现象,因此总壁面热流可分为蒸发热流密度qme(kW/m2)、淬冷热流密度qtc(kW/m2)和单相对流热流密度qc(kW/m2)。

1) 蒸发热流密度

蒸发热流是由汽泡底部微液层和过热液体层蒸发而产生的,一方面由汽泡以汽化潜热的形式带离加热壁面,另一方面用于抵抗汽泡顶部冷凝效应以维持汽泡形态与生长。根据徐平昊[22]的相关研究,汽泡顶部冷凝效应相对较小并对其进行了忽略处理,Sateesh等[9]对汽泡顶部冷凝效应也进行了忽略处理,因此本文忽略汽泡顶部冷凝效应,壁面蒸发热流密度计算如式(16)所示:

(16)

式中:dd为汽泡受力分析计算中得到的汽泡脱离直径,m;nA为汽化核心密度,m-2;f为汽泡脱离频率,s-1。

(17)

式中:tg为生长时间,即汽泡从产生到脱离的时间,s;tw为等待时间,即同一成核位点汽泡脱离后至下一汽泡出现的时间,s。在上述汽泡受力分析计算时,可得到汽泡生长时间,根据van Stralen等[23]的表示,即tw=3tg,Sateesh等[9]在进行模型建立时,也采用了上述相关性。

2) 淬冷热流密度

由于汽泡脱离后,汽泡周围冷流体填补到汽泡脱离位置,直接与加热壁面接触而发生瞬态导热产生传热。根据Sateesh等[9]的研究,其计算公式如式(18)所示:

(18)

式中:ΔTb=Tw-Tl,Tl为流体温度,℃;Ab为脱离汽泡影响面积,m2;K为经验常数,用来对汽泡投影面积进行修正得到汽泡影响面积,在Yeoh等[24]的研究成果中K的取值在1.8～2之间,本文中K取2.0。

3) 单相对流热流密度

在不受汽泡影响的区域,会发生单相强制对流,计算公式如式(19)所示:

qc=hc(1-Ab)ΔTb

(19)

式中,hc为单相对流传热系数,根据Sateesh和Chuang等[25]的研究,本文采取Churchill-Chu[26]公式计算:

(20)

式中:Gr为格拉晓夫数;Pr为普朗特数;D为当量直径。

因此,最终总热流密度qs为:

qs=qme+qtc+qc

(21)

辅助汽化核心密度计算选取Chuang等[25]考虑带有倾斜角度影响的加热平板汽化核心密度计算公式,即:

nA=-3 400+5 470ΔTsat+3.64θinc

(22)

式中,ΔTsat=Tw-Tsat,℃。

3 计算流程

计算流程如图3所示,给定初始时刻t0,依据Zuber汽泡生长模型,计算相关汽泡尺寸。

图3 计算流程图

针对倾斜工况,计算汽泡受力,依据汽泡脱离判别标准,判断汽泡受力是否达到脱离条件,若达到,受力计算结束,输出dd和tg,否则增加一步时间步长Δt,利用Zuber模型控制汽泡生长,继续计算汽泡受力,直至达到汽泡脱离判别标准,输出dd和tg。

利用汽泡受力计算结果——dd、tg和tw等参数作为壁面热流分配模型输入项,根据壁面热流分配机理模型,计算各热流密度分量,并得到倾斜工况下最终壁面热流密度计算值。

4 结果与分析

本文利用已有实验装置[27],开展倾斜加热棒束过冷流动沸腾实验。验证实验数据范围如下:倾斜角度为10°～20°;热流密度为50.63～86.08 kW/m2;流体过冷度为7.34～18.54 ℃;壁面过热度为7.53～11.80 ℃;实验段连通大气,因此实验压力为当地大气压;质量流速约为37.95～38.86 kg/(m2·s)。

定义评价指标,即平均相对误差(MRE)如下:

(23)

验证结果如图4所示,根据以前的研究结果[28],热流密度最大相对不确定度仅为0.31%,误差较小。可以看出,热流密度实验值变化趋势与热流密度计算值变化趋势基本一致,93.75%数据计算误差在±30%以内,81.25%数据计算误差在±20%以内,平均相对误差为13.07%,该模型对倾斜条件下棒束通道过冷流动沸腾热流密度具备一定的计算能力。

图4 模型计算热流密度与实验热流密度比较

但是计算模型依旧存在些许偏差,可以看出,当流体过冷度较高时,计算偏差相对较大,而当流体过冷度较低时,计算偏差相对较小。可以预见的是,当流体过冷度更大时,模型计算值将会更低,偏差将会更大。相同的结果也出现在Chuang[25]和肖波齐[29]等的实验与模型中,这说明模型本身存在计算能力上的不足与欠缺。

为了探究模型计算值和实验值偏差的原因,选取模型计算值良好的工况——倾斜10°、壁面过热度10.69 ℃、热流密度75.95 kW/m2工况作为主要参考工况,对主要参数壁面温度进行进一步分析。计算过程中,将所有工况壁温设置为参考工况壁温,保持其余热工参数不变,计算结果如图5所示。可见,在参考工况热流密度之前,由于参考工况壁温高于对应实验工况壁温,因此导致模型计算的热流密度高于原热流密度计算值;在参考工况热流密度之后,由于参考工况壁温低于对应实验工况壁温,因此导致模型计算的热流密度低于原热流密度计算值。可见,壁温变化在模型计算中起着重要作用,在实验工况范围内,遵从高壁温对应高热流密度这一物理规律。

图5 壁温影响分析

为进一步研究过冷沸腾区间壁温变化趋势对模型计算的影响,绘制了如图6所示的典型壁温变化趋势图,当热流密度处于35 kW/m2时,壁温变化趋势发生转折,同时结合实验现象观察,加热棒壁面也开始出现汽泡,即认为此处为过冷沸腾起始点(ONB点)。可以发现,在ONB点后,虽然壁温变化趋势开始由单相时的直线发生变化,但壁温并不是立即变平,而是缓慢上升一段区间,然后再逐渐变得十分平缓。而根据图5可知,在模型计算过程中,同一工况壁温的升高会导致模型计算热流密度的升高。可见,在高过冷度过冷沸腾区,壁温依旧在缓慢上升,相对于低过冷度过冷沸腾区变化十分平缓的壁温,整体水平要低,因此导致模型计算值较实验值低。而在低过冷度过冷沸腾区,壁温变化十分平缓,模型计算值与实验值符合良好。模型对于低过冷度过冷沸腾具备良好的计算能力,对于高过冷度过冷沸腾的计算仍值得进一步探究。

图6 典型沸腾曲线

由诸多汽化核心密度公式可知,汽化核心密度与壁面过热度有很大相关性。如图6所示,相比于低过冷度过冷沸腾,高过冷度过冷沸腾中,在变化同等热流密度的情况下,壁面过热度变化程度更大。根据式(22)可知,更大程度的壁面过热度变化导致汽化核心密度变化程度更大,而汽化核心密度对于模型计算热流密度影响很大。选取倾斜10°、壁面过热度10.69 ℃、热流密度75.95 kW/m2工况作为主要参考工况对汽化核心密度影响进行进一步分析。计算过程中,保持参考工况各参数不变,在实验工况范围内,改变汽化核心密度进行相关计算,结果如图7所示。可见,当汽化核心密度低于参考工况汽化核心密度时,计算热流密度偏低。汽化核心密度低于参考工况汽化核心密度越多,计算热流密度低于参考工况热流密度就越多。汽化核心密度高于参考工况汽化核心密度时同理。

图7 汽化核心密度影响分析

相比于低过冷度过冷沸腾区,高过冷度过冷沸腾区壁温变化更为明显,且壁温整体偏低,进而导致汽化核心密度变化更为明显,且汽化核心密度整体偏低。而较低的汽化核心密度会导致较低的模型热流密度计算值,造成高过冷度过冷沸腾区热流密度计算值存在一定程度的偏低。因此,在高过冷度过冷沸腾区,模型计算热流密度偏差相对较大,而在低过冷度过冷沸腾中,模型计算热流密度与实验值符合较好。汽化核心密度对于壁面过热度的依赖性对模型计算很重要,汽化核心密度与壁面过热度之间的相关变化趋势以及计算关系对模型计算影响很大。

5 结论

本文结合汽泡动力学与壁面热流分配机理模型,对倾斜加热棒束过冷沸腾壁面热流密度进行了研究,得到如下结论。

1) 在倾斜条件下棒束上加热壁面的汽泡受力分析过程中,根据实验现象,制定“浮力方向受力打破平衡,汽泡脱离”的判定准则,从而计算得到汽泡脱离直径dd、汽泡生长时间tg和汽泡等待时间tw等在壁面热流分配机理模型中所需要的相关汽泡参数。

2) 对相同工况同时建立了汽泡力平衡模型和壁面热流分配模型,并与实验值进行了对比验证。计算结果表明,以力平衡模型为基准,壁面热流分配模型计算结果在一定范围内具有良好的计算能力,平均相对误差为13.07%,误差线±30%以内涵盖93.75%的数据,±20%以内涵盖81.25%的数据。

3) 耦合汽泡力平衡模型的壁面热流分配模型在低过冷度过冷沸腾区具备十分良好的计算能力;在高过冷度过冷沸腾区域,计算偏差虽然在可接受范围之内,但相比于低过冷度过冷沸腾区偏差更大且计算值相比于实验值偏低,这是由于过冷沸腾不同区域本身特性所导致。相比于低过冷度过冷沸腾,高过冷度过冷沸腾整体偏低且较为明显的壁温变化趋势导致的汽化核心密度的改变是模型计算值偏低的主要原因。