内螺纹管临界热流密度预测与数值计算

临界热流密度(CHF)在核反应堆等换热系统的设计和安全运行中是最重要的热工限值之一，对于研究事故工况下核反应堆堆芯的安全具有实际意义。超临界水冷堆启停堆过程和事故工况下，工质必通过近临界压力区，因此必须避免发生偏离泡核沸腾(DNB)。针对光管的CHF现象，国内外学者开展了大量的实验和理论研究，得到了较多经验关联式和半机理模型[1]，但由于适用范围受限难以推广。内螺纹管相较于光管具有良好的传热特性，已得到广泛应用，但其强化传热作用在近临界压力区会受到抑制，以致DNB型CHF也会在内螺纹管中发生[2]。而目前对于内螺纹管的CHF现象尤其是DNB型CHF现象的研究甚少。

自20世纪70年代开始，大量学者对CHF进行了深入的基础研究，并尝试对其发生机理进行理论分析，基于不同假设提出了数十种机理或半机理模型[3]，其中，应用最广的是壁面汽泡壅塞模型和汽泡底部微液层蒸干模型。1983年，Weisman和Pei[4]最先提出壁面汽泡壅塞模型，他们认为当近壁面处形成的汽泡层区内汽泡足够密集时，会阻碍冷流体冷却壁面，使CHF发生。Chang和Lee[5]基于CHF发生的物理判据提出了一个适用于低干度的CHF预测模型，耦合求解两区的动量守恒方程得到极限传递质量流速。之后，Kwon和Chang[6]进一步拓展了该模型，根据大量实验数据拟合得到汽泡层区临界截面含气率的关系式，并通过水和氟利昂实验数据验证了模型。Saha和Zuber[7]针对过冷沸腾预测了净蒸汽产生点和截面含气率，认为净蒸汽产生点取决于局部条件即热力条件和流体动力条件。潘杰等[8]针对低质量流速下两相沸腾换热，建立了适用于低含气率的CHF预测模型。喻娜等[9]基于矩形通道改进了Weisman提出的汽泡壅塞模型。

目前研究发生在过冷沸腾和低干度区域的DNB型CHF的文献主要集中于光管或矩形通道等，而对内螺纹管DNB现象的研究较少，相应的CHF数据也非常少见，因此扩展内螺纹管CHF数据库，对内螺纹管在近临界压力下的DNB型CHF现象进行理论研究是很有必要的。

本文针对垂直上升内螺纹管的DNB型CHF现象，对壁面汽泡壅塞模型进行改进和修正，建立适用于近临界压力下低干度区域的CHF数值计算模型。使用Fortran语言对本文模型进行编译，研究压力、质量流速、热平衡干度和进口欠焓对CHF的影响，用CHF查表值评估本文模型，并将实验得到的内螺纹管CHF数据与Bowring模型、Katto模型、Shah模型和本文模型的CHF计算值进行比较。

1 数学模型

1.1 模型基本假设

为能建立合理简化的CHF理论模型，作如下假设：1) 流体为稳态流动；2) 过冷沸腾和低干度时，CHF是局部现象；3) 壁面汽泡层区和主流区的流动都是均匀的；4) 当汽泡层区截面含气率达到某一临界值时，径向热输运速率受到界面流动的限制，CHF发生；5) 发生CHF时，汽泡在近壁面紧密均匀排列，以单一汽泡直径作为汽泡层区有效厚度；6) 与系统压力相比，轴向压力变化可忽略不计。

1.2 控制方程

汽泡壅塞模型假定流场可分为中心主流区和壁面汽泡层区，湍流脉动使它们在两者交界面存在质量、动量和能量交换。内螺纹管的汽泡壅塞模型示意图如图1所示。

图1 内螺纹管的汽泡壅塞模型示意图
Fig.1 Schematic diagram for bubble coalescence model in rifled tube

考虑两相流模型和本构方程的影响来分析内螺纹管中的一维稳态过冷和低干度流动沸腾。汽泡层区与主流区的横向质量流速传递如图2、3所示。

图2 壁面汽泡层区控制体
Fig.2 Control volume for wall bubbly layer

将质量守恒、能量守恒和动量守恒方程应用于如图2所示的汽泡层区控制体。汽泡层区的总质量平衡方程和液体质量平衡方程如下：

(1)

[Gbc(1-xb)-Gcb(1-xc)]

(2)

式中：Gcb和Gbc分别为通过壁面汽泡层区和主流区界面流入和流出汽泡层区控制体的质量流速，kg/(m2·s)；ξi为汽泡层区和主流区界面周长，ξi=π(din-2δb)，m，δb为汽泡层区厚度，m，本文模型假设δb=Db，Db为汽泡脱离直径，din为内螺纹管有效内径，m；ξw为受热管内壁周长，ξw=πdin，m；qeav为汽泡层区中用于蒸发汽化的热流密度，kW/m2；x为干度，下标b代表汽泡层区，下标c代表主流区；hfg为汽化潜热，kJ/kg；A为横截面积，m2。

图3 主流区控制体
Fig.3 Control volume for core

CHF发生时，由于汽泡层区停滞，界面的质量传递速率是一定的，即：

Gcb=Gbc=G*

(3)

式中，G*为界面横向极限质量传递速率，kg/(m2·s)。

联立式(1)和式(2)可得：

qeva=ξihfgG*(xb-xc)/ξw

(4)

汽泡层区控制体的能量平衡方程如下：

qwξw=dGbAbhb/dz+ξi(Gbchb-Gcbhc)

(5)

式中：qw为内壁热流密度，kW/m2；h为流体焓，kJ/kg。

由式(1)和式(5)可得：

qw=Gcb(hb-hc)ξi/ξw

(6)

基于分相流模型的一维动量方程[10]，汽泡层区和主流区控制体的动量方程为：

(7)

(8)

式中：p为压力，MPa；ρ为密度，kg/m3；u为流速，m/s；τi为主流区流体对汽泡层区的黏性剪切力，N/m2；τw,v为壁面和汽泡层间的黏性剪切力，N/m2；Fd为单个汽泡受到的形态曳力，N；Nb为汽泡数量；β′为汽泡间阻塞系数。

联立式(7)和式(8)并对动量方程进行简化，忽略剪切力和加速度项的影响，可得壁面汽泡层区和主流区界面的极限传递质量流速为：

(9)

其中，ηc为主流区的体积分数：

(10)

对于给定的流道结构和进口流体条件，根据式(6)和式(9)可得CHF计算公式[7]如下：

(11)

式中，qCHF为临界热流密度，kW/m2。

分析式(11)中CHF与各参数间的关系发现，局部汽泡形态曳力、局部流体焓、质量流速、管道几何结构和系统压力是影响CHF的主要参数，因此水力和热力因素共同决定了CHF的发生[11]。

1.3 本构关系

1) 汽泡脱离点和汽泡脱离直径模型

本文模型采用Levy等[12]提出的模型来计算汽泡脱离点位置的流体欠焓和热平衡含气率。

当时，有：

(12)

当时，有：

(13)

当时，有：

(14)

式中：Δhsub,d为汽泡脱离点处的流体欠焓，kJ/kg；cpf为流体比定压热容，kJ/(kg·K)；hs为单相换热系数，kW/(m2·K)；Pr为普朗特数；下标f代表流体；G为质量流速，kg/(m2·s)；f为表面摩擦系数。

其中，壁面到汽泡顶端的无量纲距离为：

(15)

由于内螺纹管的表面摩擦系数较光管的大，不可采用光管的摩擦系数经验公式。在本文模型中，采用了工程应用中使用广泛的内螺纹管表面摩擦系数关联式[13]：

(16)

式中，Ref为雷诺数。

汽泡脱离点的热平衡含气率为：

xd,OBD=-Δhsub,d/hfg

(17)

本文假设以汽泡层区的汽泡直径为汽泡脱离直径，采用了Staub模型[14]中提出的考虑汽泡接触角影响的汽泡脱离直径计算模型，具体如下所示：

(18)

式中：σ为表面张力系数，N/m；f(θ)为汽泡接触角θ的函数，Staub建议取0.02～0.03，本文模型中取为0.025。

2) 壁面汽泡曳力和表面摩擦系数模型

壁面汽泡曳力根据如下公式[15]计算：

(19)

式中：为主流区平均流速，m/s；λ为表面摩擦系数。

考虑到内螺纹管粗糙度的影响，本文模型采用Chang和Lee[5]推荐的结合了两相雷诺数的湍流表面摩擦系数模型，即：

(20)

其中，两相雷诺数Re2φ=Gdin/μ2φ。主流区流体黏度μ2φ根据Beattie-Whalley公式计算：

μ2φ=μl(1-αtr)(1+2.5αtr)+μgαtr

(21)

式中：μ为动力黏度，N·s/m2；下标l表示饱和液体，下标g表示饱和蒸汽；αtr为真实截面含气率。

3) 真实干度和截面含气率分布模型

在汽泡脱离点下游，内螺纹管中流体强制对流的真实干度分布可通过Saha-Zuber模型[7]计算：

(22)

式中：xtr为真实干度；xeq为热平衡干度；xd,OBD为汽泡脱离点的热平衡干度。

对应的管道汽泡脱离点下游截面含气率分布则采用Zuber和Findlay[16]提出的基于两相流漂移通量模型的计算公式：

(23)

式中：C0为分布系数；Vgj为漂移速度，m/s。采用如下公式[10]计算：

C0=β[1+(1/β-1)b]

(24)

(25)

式中：b=(ρg/ρl)0.1；β=x/[x+(1-x)ρg/ρl]。

4) 汽泡层区与主流区平均速度模型

由于壁面汽泡层区厚度非常薄，可假设其速度呈线性分布。因此汽泡层区的平均速度是主流区速度即在两区界面处速度的1/2。

(26)

本文模型假设卡门单相湍流对数速度分布也适用于主流区，有：

(27)

式中：U+为无量纲速度，为汽泡层区与主流区界面的流速，m/s；Y+为无量纲距离，为表观壁面剪切力，

主流区平均流速可根据式(26)采用下式计算：

(28)

1.4 补充方程

1) 两相流动物性

基于主流区和壁面汽泡层区均为均匀流动的假设，汽泡层区和主流区流体物性计算公式如下：

(29)

ρb=ρl(1-αb)+ρgαb

(30)

ρc=ρf(1-αc)+ρgαc

(31)

hb=hl(1-xb)+hgxb

(32)

hc=hf(1-xc)+hgxc

(33)

xb=ρgαb/[ρl-(ρl-ρg)αb]

(34)

xc=ρgαc/[ρf-(ρf-ρg)αc]

(35)

式中：αb为汽泡层区临界截面含气率；αc为主流区截面含气率。

2) 汽泡层区临界截面含气率

汽泡层区临界截面含气率αb是CHF发生的判据，可能会受到加热表面条件和流动参数的影响，是汽泡壅塞模型中关键的中间参量，因此对其的计算异常重要。不同的研究者选取的αb不同，Weisman[4]假定汽泡呈椭球形紧密排列，长短轴之比为3∶1，αb为定值0.82。Chang和Lee[5]根据实验得到的水的CHF数据假定αb为定值0.75。Kwon和Chang[6]则通过分析CHF实验数据提出了一个αb的关联式，关联式中αb是热平衡干度的增函数。Kwon模型中的αb的关联式如下：

αb=0.83-0.29exp(-4.71xeq-1.89)

(36)

在本文CHF数值模型中，在Kwon提出的αb关联式的基础上还考虑了质量流速的影响，假设αb是关于真实热平衡干度和质量流速的函数，根据如图4所示Becker等[17]、Thompson和Macbeth[18]、Kwon和Chang[6]以及Zenkevich[19]所报道的大量CHF实验数据拟合得到的αb关联式如下：

αb=0.85-0.410 7×

exp(-8.172 3xtr-0.358 6G/1 000)

(37)

图4 临界热流密度与汽泡层区临界截面含气率实验值
Fig.4 Critical heat flux and critical void fraction experimental value of bubbly layer

2 数值计算方法

2.1 迭代计算过程

本文以Fortran语言为基础编写了CHF计算程序。已知压力p、质量流速G、内螺纹管径din和热平衡干度xeq，首先假设1个热流密度初值q，然后分别计算汽泡分离点位置的热平衡干度、CHF发生位置的真实干度和平均截面含气率、汽泡分离直径、表面摩擦系数、两区的平均流速、汽泡曳力及两相流物性参数等，最后由式(11)计算得到1个新的热流密度值qCHF，并与假设的初值比较，直至二者误差在规定范围内，迭代计算结束。

2.2 模型验证

CHF的Look-up Table查表法近年来被广泛接受，常被作为CHF理论研究及工程设计的参考。本文根据Groeneveld等[20]提出的CHF查询表中的压力、质量流量、干度和流道几何结构条件，采用本文的CHF模型计算相应的qCHF，并与CHF查表值进行对比以验证模型，具体对比情况如图5所示。本文模型计算值相对于查表值的平均绝对误差为23.6%。由于Groeneveld等提出的CHF查询表是针对光管的实验数据，本文模型是改进后的适用于内螺纹管的CHF模型，因此计算值较查询表中的光管CHF偏大。

图5 CHF计算值与查表值的比较
Fig.5 Comparison of calculated qCHF and Look-up Table qCHF

对φ35 mm×5.67 mm的六头内螺纹管在近临界压力范围内进行了CHF实验，实验参数范围压力p=18～21 MPa、质量流速G=500～1 000 kg/(m2·s)、进口过冷度ΔTin=3～5 ℃、内壁热负荷q=40～960 kW/m2，关于实验的详细描述见文献[21]。图6示出了采用3个经典的CHF预测模型和本文模型计算得到的qCHF与内螺纹管实验值的比较。

图6 CHF计算值与实验值的比较
Fig.6 Comparison of calculated qCHF and experimental qCHF

其中，Bowring[22]模型的平均绝对误差为16.8%，Katto[23]模型的平均绝对误差为16.3%，Shah[24]模型的平均绝对误差为12.4%，本文模型的平均绝对误差为7.0%。因此本文模型对内螺纹管的CHF预测结果误差最小，表明CHF计算值和实验值吻合较好。

3 计算结果与分析

1) 压力对CHF的影响

基于上述模型对φ35 mm×5.67 mm内螺纹管在近临界区的临界热流密度进行了预测。内螺纹管内径取实验使用的内螺纹管水力直径，din=0.021 66 m，内螺纹管长2 m，计算参数范围为压力p=18～21 MPa、质量流速G=500～1 000 kg/(m2·s)、热平衡干度xeq=-0.2～0.5。

图7示出了当质量流速为600 kg/(m2·s)时不同压力条件下的临界热流密度随真实干度的变化。

由图7可见，在低干度区域，qCHF随压力的增加而显著减小，但随真实质量含气率的增大，压力对qCHF的影响逐渐减弱。分析表明，越靠近临界压力，CHF越易发生。在近临界压力区内螺纹管的强化传热作用会被抑制，并且越靠近临界压力抑制作用越明显。

图7 G=600 kg/(m2·s)时压力对临界热流密度的影响
Fig.7 Effect of pressure on CHF at G=600 kg/(m2·s)

2) 质量流速对CHF的影响

图8示出了当热平衡干度为0.15时不同压力条件下的临界热流密度随质量流速的变化。由图8可见，临界热流密度随质量流速的增大而增大。

图8 xeq=0.15时质量流速对临界热流密度的影响
Fig.8 Effect of mass flux on CHF at xeq=0.15

3) 热平衡干度对CHF的影响

图9示出了压力为18 MPa时不同质量流速条件下的临界热流密度随热平衡干度的变化。

由图9可见，qCHF随热平衡干度的增大而逐渐减小。从图中还可发现，随热平衡干度的增大，质量流速对qCHF的影响呈减弱趋势。

4) 进口欠焓对CHF的影响

图10示出了当压力为19 MPa时不同质量流速条件下的临界热流密度随进口欠焓的变化。由图10可见，在其他影响参数不变的前提下，随进口欠焓的增大，qCHF逐渐增加。

图9 p=18 MPa时热平衡干度对临界热流密度的影响
Fig.9 Effect of equilibrium quality on CHF at p=18 MPa

图10 进口欠焓对临界热流密度的影响
Fig.10 Effect of inlet subcooling on CHF

4 结论

本文基于壁面汽泡壅塞理论，对近临界压力区垂直上升内螺纹管中以水为工质的流动沸腾DNB型CHF现象进行数值计算研究，建立了适用于低干度情况下的内螺纹管CHF预测模型，并使用Fortran语言进行编译，主要结论如下。

1) 结合内螺纹管的结构特性，对基于汽泡壅塞理论建立的适用于光管的CHF预测模型进行改进，通过CHF查表值和3个经典CHF预测模型评估本文改进后的模型，发现改进后的模型对内螺纹管CHF值的预测精度较高。

2) 近临界压力下，在低含气率区域，qCHF随压力的增大而显著减小，越靠近临界压力qCHF越小。随干度的增大，压力对qCHF的影响逐渐减弱。在相同压力和质量流速条件下，qCHF随热平衡干度的增大而减小。

3) 质量流速与进口欠焓均与qCHF呈正比关系。qCHF随质量流速和进口欠焓的增加而显著增长。随质量含气率的增大，质量流速对qCHF的影响呈减弱趋势。

4) 本文模型根据大量CHF实验数据拟合得到新的αb关联式，汽泡脱离直径的计算考虑了汽泡接触角的影响，计算表面摩擦系数采用了内螺纹管中使用广泛的经典关联式。后续将通过更多的内螺纹管实验数据得到适用于本研究所用管型的表面摩擦系数。