欠热沸腾模型参数敏感性分析

郜 冶，李小畅，刘平安

（哈尔滨工程大学 航天与建筑工程学院，黑龙江 哈尔滨 150001）

近十年来，基于欧拉-欧拉两流体六方程模型的RPI（Rensselaer Polytechnic Institute）欠热沸腾数值模型［1］得到了较大的发展，越来越多地应用于反应堆工程热工水力研究［2-3］。然而，由于沸腾机理尚不十分明确，现有的基于RPI模型的欠热沸腾数值计算方法主要建立在众多经验或半经验封闭子模型的基础上，这些子模型适用范围有限，且模型参数无通用确定方法，因此，对工况的依赖性较大。封闭子模型中的气泡脱离壁面直径、气泡脱离频率及成核面密度模型关系到气泡份额和壁面温度等重要热工参数计算结果的准确性。许多学者针对不同实验工况及流动工质提出了不同的封闭子模型，但均缺乏普适性［4］。另外，气泡脱离壁面后的运动主要受气液两相间相互作用力影响，气泡的运动反过来又会影响近壁面气泡份额和换热效率，因此两相间相互作用力模型同样十分重要。可见基于现有模型研究各封闭子模型及相间作用力模型的具体影响方式对正确使用RPI模型进行数值计算意义重大。

本文基于Nylund等［5］的FT-6a实验，研究3个重要封闭子模型（气泡脱离壁面直径、气泡脱离频率及气泡成核面密度）及两个重要相间非曳力模型（升力和湍流耗散力）对气泡轴向及径向分布和壁面过热度等计算结果的影响方式。

1 数学模型

RPI欠热沸腾数值计算方法基于欧拉-欧拉两流体六方程模型、壁面热流量分区模型及相间传输模型。本文主要给出壁面热流量分区模型及相间传输模型，欧拉-欧拉两流体六方程模型作为最基本的守恒方程本文不再给出。

1.1 壁面热流量分区模型

Kurul等［1］提出的RPI模型的核心机理是基于Judd等［6］的实验和理论研究：通道内发生壁面沸腾时可将壁面热流量分为单相对流传热、蒸发传热及气液置换传热3部分。壁面热流量分区模型基本方程如下：式中：qW、qC、qQ和qE分别为壁面总热流量、单相对流传热量、气液置换传热量和蒸发传热量；hC为单相对流换热系数，通过壁面函数求解；Tw、Tl分别为壁面温度与近壁面液相温度；kl、ρl、cpl分别为液相导热系数、密度与比定压热容；tw、f 分 别 为 气 泡 脱 离 周 期 与 频 率；ρv 为 气相密度；F 为气泡对壁面的影响因子，取值为2；Nw为壁面成核面密度；Dw为气泡脱离壁面直径；Av、Al分别为壁面上气相与液相的单位覆盖面积，Av=min（1，πF2DwNw／4），Al=1-Av；hfv为汽化潜热。由于该模型在计算中将气相温度固定为系统压力下的饱和温度，因此仅适用于气泡体积分数较低（αmax＜0.8）的情况。

式（1）中气泡脱离壁面直径、气泡脱离频率及壁面成核面密度未知，需要相应子模型对其进行封闭，本文选用如下使用相对较为广泛的模型［7］：

式中：Tsub、Tsup分别为液相欠热度和壁面过热度，Tsup=Tw-Tsat，Tsub=Tsat-Tl；Tref为参考温度；Dref为气泡脱离壁面参考直径；Nref为壁面成核参考面密度；Cf为气泡脱离阻力因子；Dmax为气泡脱离壁面时的最大直径。式（2）中多处涉及近壁面液相温度，为减小对近壁面网格的依赖，本文取y+=250时第1层网格节点上的液相温度。

1.2 相间传输模型

1）质量传输

气液两相间的质量传输包括过热壁面上的汽化和主流区的蒸发或冷凝。过热壁面上液相的汽化速率可由式（1）中的蒸发传热量得到：

气泡脱离壁面后在主流区的蒸发和冷凝速率如下：

式中：Гlv、Гvl分别为蒸发和冷凝速率；Ai为相界面积，根据气泡直径Db求解；hi为相界面换热系数，根据Ranz-Marshall关系式获得［8］。

2）能量传输

气液两相间的能量传输以相界面作为参考，两相间的热量传输为：

式中，qlv与qvl分别为液相向气相传热量及气相向液相传热量，二者大小相同，正负相反。

3）动量传输

两相间动量传输用相互作用力表示，流动沸腾中总的相互作用力FT包括曳力FD和非曳力FND两类，其中非曳力又包含升力FL、湍流耗散力FTP及壁面润滑力FWL：

各相互作用力的求解模型为：

式中：CD为曳力系数，由Ishii-Zuber关系式确定［9］；v为速度矢 量；Δ为哈密顿算子；Db为气泡直径；CL为升力系数；yw为壁面法向距离；Cw1、Cw2为 无 量 纲 参 数，分 别 取-0.01 和0.05［10］；FTD，Favre、FTD，Lopze分别为湍流耗散力 的Favre平 均 模 型［11］和Lopze模 型［12］，其 中CD与曳力模型相同，CTD为湍流耗散系数；νt、σt分别为液相湍流黏性系数及湍流Schmidt数；kl为液相湍流动能。

4）气泡直径及相界面积

本文采用Kurul等［1］推荐的基于当地液相欠热度气泡直径模型：

基于该气泡直径模型，相界面积则可采用如下对称模型获得：

2 数值计算方法

数值计算基于ANSYS CFX 14.5软件包，该程序提供了欧拉-欧拉两流体六方程模型框架，在加热壁面边界条件的处理中嵌入了RPI模型，并提供了部分相间传输模型，未提供的模型可利用User Fortran以源项的形式加入到守恒方程中。

本文以Nylund等［5］的FT-6a实验为研究对象，实验模型及尺寸参数如图1a所示。实验以水为工质，系统压力为5.0MPa，入口欠热度为4.5K，入口流量为1 163kg／（m2·K），均匀壁面热流量为522kW／m2。取实验模型周向的1／10和轴向加热段前2m 作为数值计算对象，如图1a中阴影部分所示。气相物性参数取系统压力下的饱和蒸汽参数；液相则取进出口的平均值。液相湍流模拟采用k-ω SST 模型［13］，气相则采用零方程模型。计算中建立了3套网格：286×50、391×100、605×200，各网格计算的平均气泡份额差值小于3%，选取第2套网格作为计算对象，如图1b所示。

图1 FT-6a几何模型及计算网格Fig.1 Sketch of FT-6ageometry and computational mesh

3 计算结果分析

分析RPI模型中3 个子模型（Dw、Nw及f）及两个相间非曳力模型（FL及FTD）对计算结果的影响。从式（1）及（2）可看出，RPI模型主要由Dw、Nw及f 确定，对于一具体的流动沸腾，这些封闭子模型主要由式（2）中的Dref、Cf及Nref决定。从式（7）则可看出，升力主要受升力系数CL影响，Favre湍流耗散力模型参数CD已由曳力模型确定，Lopze湍流耗散力模型则主要受耗散系数CTD影响。因此，针对上述所选定的RPI封闭子模型及相间非曳力模型的研究可通过分析如下5个参数对数值计算结果的 敏 感 性 来 实 现：Dref、Cf、Nref、CL及CTD。根据文献［7］，5个参数推荐值列于表1。

表1 模型参数参考值Table 1 Reference values of model parameters

FT-6a实验测量了气泡沿轴向和径向的分布，未给出壁面过热度，本文采用Jens 关联式［15］估算充分发展的欠热流动沸腾壁面过热度：

式中：q 为 壁 面 热 流 量，MW／m2；p 为 系 统 压力，Pa。基于式（10）所计算的FT-6a实验壁面过热度为9.49K。另外，Bartolomej等［16］基于圆管的类似工况（4.5 MPa，0.57 MW／m2）下的欠热沸腾实验结果约为10K，因此可基本确定Jens关联式计算结果的准确性。

3.1 子模型参数敏感性分析

计算时升力及湍流耗散力分别采用Tomiyama和Favre模型，其他子模型参数取表1中的值。图2为Dref对轴向平均气泡份额及壁面过热度的影响。从图2可看出，随参考直径的增大，平均气泡份额逐渐增大，壁面过热度则随之降低。这是由于随壁面气泡直径的增大，蒸发传热量及气液置换传热量随之增大，由蒸发传热产生的气泡份额增加，此时传热效率升高，壁面温度降低。参考直径达1.3mm 时，平均气泡份额与实验结果符合较好，但壁面过热度大幅低于Jens关联式计算值，参考直径取0.9mm 时壁面过热度与关联式符合较好。该图表明Dw对平均气泡份额和壁面温度的计算结果影响均较大。

图3示出Nref对平均气泡份额及壁面过热度的影响。图3表明，随Nref的增大，气泡份额逐渐增大，壁面过热度则随之减小。显然这与物理事实相符，壁面成核点越多，气泡产生量越多，传热效率提高，壁面温度降低。当Nref＞8×105m-2时，计算结果不再变化，这是由于式（1）中的qE受max函数限制。在Nref的有效影响范围内，成核面密度对气泡份额及壁面温度影响均较明显。

图4为Cf对平均气泡份额及壁面过热度的影响。从图4 可看出，当减小阻力因子时（＜1），该参数对气泡份额的计算结果几乎无影响，对壁面过热度的影响也较小。当增大阻力因子时（＞1），气泡份额有所下降，壁面过热度则大幅升高。

图2 Dref对轴向平均气泡份额及壁面过热度的影响Fig.2 Influence of bubble departure diameter on axial average void fraction and wall superheat

图3 Nref对平均气泡份额及壁面过热度的影响Fig.3 Influence of wall nucleation site density on axial average void fraction and wall superheat

图4 Cf 对平均气泡份额及壁面过热度的影响Fig.4 Influence of bubble detachment frequency on axial average void fraction and wall superheat

图5 RPI子模型对线平均气泡份额径向分布的影响Fig.5 Influence of RPI sub-models on radial distribution of line-averaged void fraction

图5为RPI各子模型计算的线平均气泡份额在1 598mm 高处径向分布与实验结果的对比。从图5可看出，合适的模型参数可得到与实验数据基本一致的计算结果。图中各曲线基本平行说明各封闭子模型对气泡的径向分布无任何影响，只影响气泡份额大小。

综合上述，脱离壁面气泡直径及成核密度对气泡份额及壁面过热度计算结果影响相对较大，气泡脱离频率对二者的影响相对较小。不能单独对比某一参数的实验值来确定计算结果的可靠性。对于FT-6a 实验，推荐参考直径0.9mm，参考成核面密度8×105m-2，气泡脱离频率阻力因子1.0。

3.2 非曳力模型参数敏感性分析

通过气泡径向分布计算结果来分析非曳力模型参数敏感性。图6为采用不同升力系数计算得到的气泡径向分布与实验结果的对比。从图6可看出，不同升力系数主要影响棒束间隙（0.01m）及绝热壁面附近（0.035m）气泡份额大小。升力系数越大，棒束间隙处气泡份额越小，绝热壁面附近气泡份额越大。Tomiyama关系式与实验符合较好，CL取0.5 时与Tomiyama关系式的计算结果接近。

图6 升力系数敏感性分析Fig.6 Sensitivity analysis of lift coefficient

图7为不同湍流耗散力模型计算结果与实验结果的对比。从图7可看出，与升力不同，湍流耗散力主要影响加热壁面附近气泡份额大小。从Lopez模型可看出，湍流耗散系数越大，加热壁面附近气泡份额越小。这意味着湍流耗散力有促进气泡远离加热壁面的作用。Favre模型与CTD取0.5时的Lopez模型计算结果接近，且与实验数据符合较好。

图7 湍流耗散力模型敏感性分析Fig.7 Sensitivity analysis of turbulent dispersion force

为进一步观察相间非曳力模型对气泡径向分布的影响，图8示出了1 598 mm 截面处升力及湍流耗散力对气相分布的影响。从图8a可看出，CL较小（0.01）时，气泡更多地向主流区运动；CL较大（1.0）时，气泡被束缚在壁面附近，Tomiyama模型计算的结果则处于二者之间。图8b表明，CTD较小（0.1）时，气泡聚集在壁面附近；CTD较大（0.5）时，气泡则向主流区扩散。分析结果与图7基本一致。

图8 升力（a）及湍流耗散力（b）对气相分布的影响Fig.8 Influence of lift force（a）and turbulent dispersion force（b）on void fraction contours

4 结论

1）气泡脱离壁面直径模型对气泡份额及壁面温度的计算影响均较大，随壁面气泡直径的增大，气泡份额增加，壁面温度降低。

2）成核面密度在其有效影响范围内对壁面过热度的影响较大，对气泡份额相对较小。成核面密度越大，气泡份额越大，壁面过热度越低。

3）对气泡脱离频率而言，当气泡脱离阻力因子在小于1的范围内变化时，对气泡份额的计算结果几乎无影响，对壁面过热度的影响也较小；当阻力因子大于1时，则对二者的影响均较大。

4）各RPI封闭子模型对气泡径向分布无任何影响，只影响气泡份额大小。

5）不能通过对比单个参数的实验测量值来验证计算的可靠性，应综合对比多个实验值，以确定最佳模型参数。

6）相间非曳力模型对气泡径向分布的计算影响较大，升力有阻止气泡离开加热壁面的作用，湍流耗散力有促进气泡向主流区运动的作用。

本文可为今后改善欠热沸腾数值计算的准确性在模型及参数选取方面提供参考。

［1］ KURUL N，PODOWSKI M Z.On the modeling of multidimensional effects in boiling channels［C］∥27th National Heat Transfer Conference.USA：ANS，1991：28-31.

［2］ LO S，OSMAN J.CFD modeling of boiling flow in PSBT 5×5bundle［J］.Science and Technology of Nuclear Installations，2012（1）：1-8.

［3］ 李小畅，郜冶.压水堆子通道欠热沸腾数值验证及交混翼研究［J］.原子能科学技术，2013，47（12）：2 208-2 215.LI Xiaochang，GAO Ye.Numerical validation of subcooled boiling in subchannel of PWR and Research on mixing vane［J］.At Energy Sci Technol，2013，47（12）：2 208-2 215（in Chinese）.

［4］ CHEUNG S C P，VAHAJI S，YEOH G H，et al.Modeling subcooled flow boiling in vertical channels at low pressures，Part 1：Assessment of empirical correlations［J］.International Journal of Heat and Mass Transfer，2014，75（1）：736-753.

［5］ NYLUND K M，BECKET R，EKLUND O，et al.Measurements of hydrodynamics characteris-tics，instability thresholds，and burnout limits for 6-rod clusters in natural and forced circulation［R］.Sweden：FRtGG-I，1967.

［6］ JUDD R L，HWANG K S.A comprehensive model for nucleate pool boiling heat transfer including microlayer evaporation［J］.ASME J Heat Transfer，1976，94（4）：623-629.

［7］ KREPPER E，KONCAR B，EGOROV Y.CFD modeling of subcooled boiling concept，validation and application to fuel assembly design［J］.Nucl Eng Des，2006，237（7）：716-731.

［8］ RANZ W E，MARSHALL W R.Evaporation from drops［J］.Chem Eng Prog，1952，48（4）：173-180.

［9］ ISHII M，ZUBER N.Drag coefficient and relative velocity in bubbly，droplet or particulate flows［J］.AIChE Journal，1979，25（5）：843-855.

［10］ANTAL S P，LAHEY R T，FLAHERTY J E.Analysis of phase distribution in fully developed laminar bubbly two-phase flow［J］.Multiphase Flow，1991，17（5）：635-652.

［11］BURNS A D B，FRANK，T，HAMILL I，et al.The Favre averaged drag model for turbulent dispersion in Eulerian［C］∥5th International Conference on Multiphase Flow.Yokohama，Japan：ICMF，2004：1-17.

［12］LOPEZ D B.Turbulent bubbly flow in a triangular duct［D］.New York：Rensselaer Polytechnic Institute，1991.

［13］KONCAR B，KREPPER E，EGOROV Y.CFD modeling of subcooled flow boiling for nuclear engineering applications［C］∥International Conference of Nuclear Energy for New Europe.Bled Slovenia：NSS，2005：1-14.

［14］TOMIYAMA A，TAMAI H，ZUN I，et al.Transverse migration of single bubbles in simple shear flows［J］.Chem Eng Sci，2002，57（11）：1 849-1 858.

［15］JENS W H，LOTTES P A.Analysis of heat transfer，burnout，pressure drop and density data for high pressure water，ANL-4627［R］.Michigan：ANL，1951.

［16］BARTOLOMEJ G G，CHANTURIYA V M.Experimental study of true void fraction when boiling subcooled water in vertical tubes［J］.Thermal Engineering，1967，14（2）：123-128.