汽泡边界层模型在矩形管道下的适用性分析

何 雯,赵陈儒,韩晋玉,薄涵亮

(清华大学 核能与新能源技术研究院,北京 100084)

过冷流动沸腾由于其较高的换热效率广泛应用于石油、化工、核电等领域[1-2]。沸腾过程伴随着复杂的汽泡行为,如汽泡的核化、生长、浮升和滑移,欲探究过冷流动沸腾下的流动和换热机理,需要对这些汽泡行为进行准确的描述。矩形管道是核反应堆板型燃料组件常用的管道类型,长宽比很大,窄缝宽度通常为2～3 mm[3]。受长宽比的影响,矩形管道下的流动特性和汽泡行为可能与圆管不同,因此值得深入研究。

空泡份额是过冷流动沸腾中一个重要的两相参数,该参数不仅会影响流场的流动和换热特性,如压降和换热效率,还是流型划分的重要依据[4]。研究空泡份额对分析过冷流动沸腾下的汽泡行为和流动特性都具有重要意义。除了实验和数值模拟,理论模型/经验关系式也是预测该参数的一种重要方法。目前,预测空泡份额的模型大致可以分为4种:修正的均相流模型、滑速比模型、漂移流模型和其他经验关系式/模型[5]。修正的均相流模型主要基于空泡份额和含气率的定量关系,在均相流模型的基础上进行修正,如Bankoff模型[6]、Chisholm模型[7]。滑速比模型假设汽液两相以单独的速度流动,采用滑速比S反映两相流速之间的关系,进而获得空泡份额与工况和含气率之间的关系,如Thom模型[8]、Winterton模型[9]。漂移流模型将流场分布的不均匀性和两相局部的速度差考虑在内,采用分布参数C和漂移速度ugj反映这两个因素的影响,如Saha&Zuber模型[10]、Ishii模型[11]。第四类模型则主要基于理论分析或基于实验数据拟合而来,如Levy模型[12]。这些模型和关系式的数学形式及适用范围在文献[5]中进行了详细的整理。由此可见,预测空泡份额的模型很多。然而,这些模型几乎都是一维的,即只能用于描述空泡份额沿轴向的变化情况,并且很少能将汽泡行为考虑在内,而流场的流动和换热特性与汽泡行为息息相关[13-14]。

受壁面热量的影响,近壁面区域发生着汽泡的核化、生长、脱离等,汽泡行为更加复杂。描述近壁面汽泡行为的动力学参数主要有4个:生长速率、脱离直径、脱离频率和活化核心密度[15-17]。随着汽泡动力学研究的逐渐成熟,He等[18-19]针对流动沸腾构建了一种新的理论模型,即汽泡边界层模型。区别于一维模型,汽泡边界层模型将流场沿径向划分为两个区域,即主流区域和汽泡边界层区域,汽泡边界层的厚度等于汽泡脱离直径。采用汽泡动力学参数描述边界层内的汽泡行为,然后通过一组准二维基本方程实现主流和汽泡边界层的双向耦合,进而可以同时获得两个区域内多个两相参数的变化情况。He等[18]首先将模型应用于过冷流动沸腾,模型不仅描述了两个区域内空泡份额等参数沿轴向的变化情况,还获得了两个区域沿径向的质量和能量交换情况,该模型已成功应用于低温供热堆燃料元件内的热工分析[20]。He等[19]继续将模型推广至饱和流动沸腾中,包含弹状流、搅浑流和环状流,模型同样获得多个两相参数的变化规律。由此可见,相较于一维模型/经验关系式,汽泡边界层模型不仅能将近壁面的汽泡行为考虑在内,还能获得不同径向区域多个两相参数的变化情况,为过冷流动沸腾的预测提供了一种新的方法。

目前,过冷沸腾段的汽泡边界层模型仅在常规圆管中进行了验证。对于长宽比较大的矩形管道,汽泡在管内的生长可能受到管道的限制和挤压,尤其是矩形窄缝,表现出不同的流动和换热特性。因此,过冷沸腾段汽泡边界层模型在矩形常规管道和窄缝通道下的适用性还有待进一步评估。

1 汽泡边界层模型

1.1 物理模型

过冷流动沸腾起始于单相过冷液体,当第1个汽泡在壁面核化(ONB点)后,流场进入高过冷沸腾段。在此流段,汽泡仅附着在壁面,不脱离也不滑移[20]。随着热量的不断输入,液体过冷度逐渐下降,汽泡开始脱离核化点(OSV点),此时流场进入低过冷沸腾段。此流段中,汽泡脱离核化点后部分沿壁面滑移,部分浮升进入主流,由于主流液体依然处于过冷状态,进入主流的汽泡会被冷凝。当流体温度达到饱和温度后,流场达到Tsat点,过冷沸腾段结束。

由于低过冷沸腾段汽泡行为更复杂,空泡份额变化更明显,因此,本文主要针对低过冷沸腾段,对矩形管道构建汽泡边界层模型。以单面加热矩形管道为例,具体如图1所示。基于汽泡边界层的思想,矩形管道沿径向划分为主流区域和汽泡边界层两个区域,两个区域的流场不仅沿轴向发生变化,径向上两个区域间也随时发生着质量、能量和动量交换。汽泡尺寸达到脱离直径后,开始脱离核化点,部分沿壁面滑移,部分直接浮升进入主流。由于汽泡滑移过程中其尺寸和形状几乎不发生变化[18],因此,取汽泡边界层的厚度为脱离直径。该值受工况影响,但工况一旦确定后认为其沿程不发生变化。脱离直径采用文献[17]构建的关系式计算,即:

(1)

图1 单面加热矩形管道内的汽泡边界层物理模型

其中:ρl和ρg分别为液体和汽体的密度;σ、cp,l和ilg分别为表面张力、液体比定压热容和汽化潜热;Δtw为壁面过热度。该关系式同时适用于窄缝和常规管道,管道类型包括圆管、环管和矩形管道,适用的工况范围为:水力直径d,1～42.2 mm;压力p,0.101～0.15 MPa;壁面过热度Δtw,0.5～50.1 K;液体过冷度Δtsub,2～50.1 K;质量流速G,67～1 927 kg/(m2·s);普朗克数Pr,1.29～7.68;汽泡雷诺数Reb,31～1 170;过热雅可比数Jaw,0.000 8～0.204。

1.2 基本假设

为了简化计算,需要对低过冷沸腾段流场进行相应的简化。首先,假定流场稳定且分布均匀,并且忽略截面轴向压力波动和径向压力变化对物性的影响。其次,由于主流液体处于过冷状态,因此假设汽泡仅在边界层内产生和生长。随着热量的不断输入,边界层平均温度逐渐上升,汽泡核化的数量也在增加。考虑到边界层很薄,因此假设边界层内的温度和空泡份额均呈线性上升。低过冷沸腾段起始于OSV点,结束于Tsat点,由于边界层内液体温度逐渐增加,因此两点之间边界层内的空泡份额呈上升趋势。对于OSV点,本文采用Okawa[20]的判定标准,Okawa 认为该点后流场空泡份额快速上升的原因与流型的转变有关,于是基于汽泡聚并的条件,认为近壁面相邻汽泡间的距离恰好等于2倍汽泡直径,对应的空泡份额为0.3。而对于Tsat点,由于目前尚不清楚饱和沸腾段汽泡在近壁面的空间分布情况,He等[19]在针对饱和沸腾段构建汽泡边界层模型时,选取不同的边界层空泡份额,结果发现,取值0.5时模型的准确度更高。因此,本文取Tsat点处边界层空泡份额为0.5。此外,汽泡在壁面稳定存在时,其受到的蒸发作用和冷凝作用相平衡。在此基础上,He等[18]假定汽泡稳定附着在壁面(即高过冷沸腾段)时,汽泡周围流场的平均温度恰好为饱和温度。由于OSV点是高过冷沸腾段的终点,汽泡尺寸恰好为脱离直径,因此假定该点处边界层内的平均温度为饱和温度。而对于Tsat点,He等[19]对饱和沸腾段构建了汽泡边界层模型,建议取Tsat点处流体过热度为壁面过热度的1/3。

综上,针对低过冷沸腾段,所有的基本假设如下:1) 流场稳定且分布均匀;2) 忽略截面轴向压力波动和径向压力变化对物性的影响;3) 汽泡仅在汽泡边界层内产生和生长;4) 汽泡边界层内流体平均温度沿轴向线性上升,OSV点处流体温度等于饱和温度,Tsat点处流体过热度等于壁面过热度的1/3;5) 汽泡边界层内空泡份额沿轴向线性上升,OSV点处为0.3,Tsat点处为0.5。

1.3 数学模型

同样以单面加热的矩形管道为例,将图1中的物理模型转化为数学形式,具体如图2所示。由于边界层区域的汽液两相流速和温差较小,这个区域的基本方程基于均相流模型。而主流区域的汽液两相流速和温差均较大,因此这个区域的基本方程基于分相流模型。分相流模型的基本思想是将汽液两相人为地划分开,并认为两相间的质量和能量传递都发生在相界面。因此,主流区域的流场也被二次划分为主流液体区域和主流汽体区域[21]。

图2 单面加热矩形管道的汽泡边界层数学模型

分别针对主流汽体区域、主流液体区域和边界层区域构建质量、动量和能量方程。质量方程主要考虑了径向质量交混、冷凝等对各区域质量流量的影响。这3个区域的质量守恒方程分别如下:

Mbcxb-Mcon=M′cg-Mcg

(2)

Mbc(1-xb)+Mcon-Mcb=M′cl-Mcl

(3)

Mcb-Mbc=M′b-Mb

(4)

其中:Mbc和Mcb分别为边界层进入主流以及主流进入边界层的质量流量;xb为边界层含气率;Mcon为汽泡进入主流后发生冷凝的质量流量;Mcg、Mcl和Mb分别为主流汽体区域、主流液体区域和边界层区域的质量流量。

动量方程则考虑了各力对流场的影响,如流体与壁面间的摩擦应力(τw,v)、主流和边界层之间的剪切应力(τi)以及主流汽液两相之间的曳力(F)。3个区域的动量守恒方程分别如下:

-Acg(p′-p)MbcxbUb-MconUcg-

F-ρcggAcgn=M′cgU′cg-McgHcg

(5)

-Acl(p′-p)+Mbc(1-xb)Ub+

MconUcg-McbUcl-τilζin+F-

ρclgAcln=M′clU′cl-MclHcl

(6)

-Ab(p′-p)+McbUcl-MbcUb+τiζin-

ρbgAbn-τw,vζwn=M′bU′b-MbHb

(7)

其中:p为系统压力;ζi为主流和边界层交界面的周长;ζw为壁面周长。

能量方程则考虑了壁面输入的热量(qw)、各区域间的导热以及径向质量交换带来的热量传递等对流场能量的影响。3个区域的能量守恒方程分别如下:

MbcxbHbg-HcgMcon=M′cgHcg-McgHcg

(8)

Mbc(1-xb)Hbl+HcgMcon-

McbHcl+Qc=M′clH′cl-MclHcl

(9)

qwξwn-Qc+McbHcl-MbcHb=

M′bHb-MbHb

(10)

其中:Hcg和Hcl分别为主流区域汽相和液相对应的焓;Hbg和Hbl分别为边界层区域汽相和液相对应的焓;Qc为汽泡边界层和主流区域间的导热量。

针对矩形管道,导热量Qc计算如下:

Qc=λl(Tb-Tc)l

(11)

其中:Hb和Hc分别为汽泡边界层和主流区域的温度;l为两个区域中心点的距离;λl为液体导热系数。

由于Hcg=Hbg,方程(2)和(8)一致,因此基本方程有8个。而未知数也有8个,分别为3个区域的流速(Ucg,Ucl,Ub)、主流和边界层的径向交换量(Mcb,Mbc)、主流液体温度(Tcl)、主流空泡份额(αc)以及压力(p)。方程中汽泡进入主流后的冷凝量(Mcon)计算如下:

(12)

其中,kcon为冷凝导热量,采用Rouhani[22]经验关系式计算:

(13)

其中:Cs为加热周长;α(z)为不同轴向位置处的截面空泡份额。

此时基本方程数和未知数个数相等,模型可直接求解。然而,动量方程中包含了主流汽液两相间的曳力(F),由于主流内汽泡数量和形状一直在发生变化,两相界面也随时发生着变化,导致F很难获得。因此,在模型求解过程中,将主流汽液两相的动量方程合并,补充一个滑速比公式(He等[19])用于反映主流区域两相流速的变化规律,具体如下:

(14)

综上,数学模型由7个基本方程和1个滑速比方程组成,方程数和未知数个数相等,模型可顺利求解。

2 模型验证

空泡份额是流动沸腾中一个重要的两相参数。对于低过冷沸腾段,汽泡边界层模型曾与常规圆管内的空泡份额实验数据进行对比[18],验证了模型在常规圆管下的准确性。而对于矩形管道,本文同样将模型与空泡份额实验数据[23-25]进行比较。首先,针对常规尺寸的矩形管道,Christensen[23]开展了过冷流动沸腾下的实验,管道四面加热,热流密度恒定,管道尺寸为11.1 mm×44.4 mm,工质为水。实验涉及的工况范围为:压力2.76～6.89 MPa、热流密度211～495 kW/m2、质量流速637～915 kg/(m2·s)。由于加热壁面处有汽泡产生,因此4个加热面均存在汽泡边界层,厚度等于汽泡脱离直径,而其他区域则为主流区域,分别对两个区域构建质量、能量和动量基本方程。

区别于现有一维模型/经验关系式,模型可以同时获得两个区域内空泡份额的沿程变化情况。截面平均空泡份额(α)等于两个区域空泡份额的面积平均,即:

(15)

其中:Ab和Ac分别为汽泡边界层和主流区域的截面积;αb和αc分别为两个区域的空泡份额;A为整个管道的截面积。

不同工况下,αb、αc和α的变化示于图3。根据基本假设(5),αb从0.3线性增长到0.5。由于实验压力较大,导致汽泡脱离直径(Dd)数值较小[17]。因此,相较于管道尺寸,汽泡边界层厚度很小,主流区域占比较大,使得主流空泡份额αc和截面平均空泡份额α的变化曲线几乎重合。模型的误差则采用截面平均空泡份额α与实验数据的误差来定量反映,相对误差(eR)的绝对值计算如下:

图3 汽泡边界层模型的空泡份额预测值与Christensen[23]实验值对比

(16)

其中,αexp和αpre分别为空泡份额的实验值和预测值。每组工况对应的相对误差列于表1。

表1 模型空泡份额预测误差与Cai模型的对比

此外,还将模型应用于矩形窄缝通道(通常指宽度小于/等于3 mm的矩形管道)中,验证模型在窄缝通道下的准确性。Kureta等[24]针对工质水开展了其在矩形窄通道下的过冷流动沸腾实验,实验中按0.89 ms的时间间隔对流场进行拍照,进而获得空泡份额的空间分布。管道尺寸为3 mm×20 mm,沿20 mm的长边进行单面加热。基于此,汽泡边界层仅存在于该加热面上,剩余区域均为主流区域。Martin[25]同样针对矩形窄缝通道,开展了工质水在高压下的流动沸腾实验。测量设备主要由1个西门子X射线发生器和2个光电倍增器组成,整个系统连接在一个水平支架上。支架可以在垂直于流动的方向上以可调节的速度(1 mm/min)移动,采用X射线每25 μm对流道进行1次记录,进而获得空泡份额的分布。管道尺寸为2.8 mm×50 mm,沿长边双面加热,此时汽泡边界层同时存在于两个长加热面,剩余流域为主流区域。将模型应用到两组实验下,各区域的空泡份额变化分别示于图4和图5。

图4 汽泡边界层模型的空泡份额预测值与Kureta[24]实验值的对比

图5 汽泡边界层模型空泡份额预测值与Martin[25]实验值对比

从图4、5可见,尽管管道尺寸变小,汽泡边界层模型依然能准确地描述不同工况下空泡份额的变化规律,并且涵盖的工况范围很大,其中压力为0.1～8 MPa,热流密度为575～1 700 kW/m2,质量流速为240～2 200 kg/(m2·s)。不同工况下模型对应的相对误差列于表1。尽管在部分工况下模型的误差较大,如对于Martin-3工况,eA达到70.7%,但从图5c可以发现,模型对该工况下空泡份额的变化规律描述得很准确。误差较大的原因在于该流段下空泡份额的数值很小,较小的差异都会造成相对误差的数值较大。由此可见,汽泡边界层模型受管道类型和尺寸的影响很小,在矩形常规管道和窄缝通道中都具有较高的预测精度。原因主要有以下两点:一方面,汽泡边界层模型基于理论分析而来,数学模型也以基本方程为主,模型本身受管道结构的影响很小;另一方面,模型将流场划分为主流和汽泡边界层两个区域,认为两个区域内的流场和汽泡行为相对独立。尽管管道变成了矩形窄缝通道,但由于过冷沸腾下流场大多处于泡状流,汽泡的尺寸不大,数量也不多,基本不会出现汽泡被管道挤压后充斥整个管道的情况。因此,模型的划分方式依然有效。

为更直观反映汽泡边界层模型的准确性,除了与实验数据对比外,还将模型与现有的Cai模型[5]进行了对比。针对过冷流动沸腾,Cai等将大量现有理论模型/经验关系式的准确度进行了对比,选出一套准确度最高的模型用于预测空泡份额的变化。该模型由多个理论模型/经验关系式组合而成:首先,采用Saha&Zuber关系式[10]确定OSV点的位置;然后,结合流场的热平衡干度,采用Saha&Zuber关系式[10]得到含气率的变化规律;最后,通过Cai自行构建的模型或Thom关系式[8]得到空泡份额的变化规律。Cai模型应用于各组工况下的误差也列于表1。从对比结果可见,对于矩形管道下的过冷流动沸腾,汽泡边界层模型不仅较Cai模型的准确度更高,相对误差分别为32.7%和52.6%,还将空泡份额的描述从一维变成准二维,能同时描述主流和汽泡边界层区域的空泡份额变化情况,对流场的描述更加细致。

3 其他两相参数

除了空泡份额,汽泡边界层模型还能获得多个两相参数的沿程变化规律。流速的变化不仅可以为流型划分提供重要的参考依据,还能反映主流汽液两相间曳力以及壁面摩擦力对流场流动特性的影响。以矩形窄缝通道Martin-1工况为例,其低过冷沸腾段两相参数变化示于图6。图6a为主流液体、主流汽体和汽泡边界层平均流速的变化。从图6a可见,随着热量的不断输入,流场汽相含量逐渐增加,3个区域的流速均增加。受浮力的影响,主流汽相的流速最大,并与主流液相形成速度差,由滑速比S体现(图6b)。由于过冷沸腾段流场大多处于泡状流,汽相含量不大,因此滑速比的增长范围也不大。图6c为主流和汽泡边界层区域温度的变化。由于该工况下壁面过热度不高,汽泡边界层的温度在这个流段内变化不大。相较之下,主流区域的温度上升明显,流体从过冷状态一直被加热到饱和温度。图6d为压力的变化,综合反映了重力、摩擦力等对流场的影响。通常情况下,实验很难准确获得近壁面的流场情况,尤其是对于矩形窄缝通道,而数值模拟方法的计算量较大。相较之下,汽泡边界层模型可以提供过冷沸腾段详细的流场信息,为过冷流动沸腾的计算提供了一种新的方法。此外,过冷沸腾段沸腾起始点(ONB点)、汽泡脱离点(OSV点)等都与近壁面的流场息息相关,因此,模型获得的近壁面流场信息还能为这两个点的分析提供基础数据。

图6 矩形窄缝通道内低过冷沸腾段两相参数的变化

4 结论

汽泡边界层模型是预测过冷流动沸腾的一种新理论模型,为验证模型的准确性,本文将模型同时应用于矩形常规管道和窄缝通道内空泡份额的预测中。结果显示,汽泡边界层模型在矩形管道下的准确度依然较高,与空泡份额实验数据相比相对误差为32.7%,而现有空泡份额预测模型的相对误差仅52.6%。验证工况范围为:压力0.1～8 MPa、热流密度211～1 700 kW/m2、质量流速240～2 200 kg/(m2·s)、液体雷诺数4.42×103～3.72×105、普朗克数0.83～1.76。由此可见,相较于现有预测模型,汽泡边界层模型不仅准确度更高,受管道结构和尺寸的影响较小,还能同时获得主流和汽泡边界层两个区域内多个两相参数的变化规律,能对流场进行更细致的描述,为未来进一步应用于核反应堆矩形燃料组件内的热工分析和计算奠定了基础。