压水堆燃料棒束通道内过冷沸腾分析

张 蕊，干富军，左巧林，田文喜，苏光辉，秋穗正

（1.西安交通大学动力工程多相流国家重点实验室，陕西西安 710049；2.上海核工程研究设计院，上海 200233）

压水堆燃料棒束通道内过冷沸腾分析

张 蕊1，干富军2，左巧林2，田文喜1，苏光辉1，秋穗正1

（1.西安交通大学动力工程多相流国家重点实验室，陕西西安 710049；2.上海核工程研究设计院，上海 200233）

使用Fluent14.5两流体模型中的RPI（Rensselaer Polytechnic Institute）壁面沸腾模型，对堆芯燃料棒束通道内过冷沸腾现象进行数值模拟，得到了通道内的流场、温度场以及空泡份额的分布，分析了定位格架和搅混翼的存在对热工水力特性的影响。数值结果表明，格架的存在会造成很大的压降，而搅混翼会对流场、温度场和空泡份额分布产生显著影响；RPI壁面沸腾模型的模拟结果与Bartolemei试验数据符合很好。

过冷沸腾；燃料棒束；定位格架

过冷沸腾指壁面温度高于饱和温度，而流道内的液体仍处于过冷状态的流动。核电厂的工程应用设备（如反应堆堆芯、蒸汽发生器）中也有可能出现过冷沸腾。过冷沸腾可极大地提高换热表面的换热系数，但由于换热表面热流密度较高，易发生传热恶化，当过冷沸腾发生偏离核态沸腾（DNB）现象时，由于加热表面被气膜覆盖，加热表面温度瞬间飞升并发生烧毁。因此，进行堆芯流道内过冷沸腾分析很有必要。目前常规的单通道和子通道分析，不能详细地描述燃料组件内的局部流场细节特征。随着计算机技术的飞速发展，大型计算流体动力学商用软件Fluent在诸多工程领域获得了广泛应用。因此，采用CFD方法对反应堆燃料组件过冷沸腾模拟显得十分可行而且必要，同时也具有重要的工程价值和学术意义。

目前燃料棒束通道内CFD模拟主要集中在单相，Navarro等［1］用商用软件CFX11.0模拟了5×5棒束带定位格架结构内的单相流动。Liu等［2］采用商用软件Fluent12.0分析了不同湍流模型（RNG k-ε模型、realizable k-ε模型、低Re k-ε模型、标准k-ε模型、SST k-w模型和雷诺应力模型）对5×5棒束通道计算结果的影响。Li等［3］对含搅混翼燃料组件特征通道内的单相流动进行了模拟研究。

而对燃料棒束通道内两相流动换热现象的研究还十分少见。因此，本文应用Fluent14.5计算软件，采用Kurul等［4］所提出的RPI壁面沸腾模型，对燃料棒束单通道带搅混翼和定位格架过冷沸腾现象进行模拟分析。

1 数学物理模型

过冷沸腾的数值分析基于欧拉两流体模型进行，对气相和液相分别求解质量、动量、能量守恒方程，同时求解相间质量、动量、能量传递模型。

1.1 欧拉两流体模型

质量守恒方程：

动量守恒方程：

能量守恒方程：

式中：αi、ρi、vi、Si、pi、=τi、hi和qi分别为i相的体积份额、密度、速度、源项、压力、应力张量、比焓和热流量；˙mij和Qij分别为i相至j相的质量和能量传递。

1.2 相间动量传递模型

气相、液相间作用力包括相间拽力、升力、壁面润滑力、湍流耗散力，拽力FD为单位体积内气泡施加在液相上的力：

其中：CD为升力系数，由Ishii关系式［5］计算；μf为液相黏度；Aif为相界面密度：

Re为相间相对雷诺数：

式中：dg为气泡直径；vg和vf分别为气相、液相速度；αg为气相体积份额；ρf为液相密度。

相间升力为液相速度梯度施加在气相上的力，由下式计算：

其中，CL为升力系数，由Moraga关系式［6］计算。

壁面润滑力作用在气泡上，推动气泡向主流方向运动，由下式计算：

其中：Cwl为壁面润滑力系数，由Antal模型计算［7］；nw为壁面法向矢量。

湍流耗散力用来描述两相间的湍流相互作用，由下式计算：

其中：CTD＝1；σfg＝0.9。

1.3 相间能量传递

气相脱离壁面后，进入过冷流体时向液体释热，单位体积内的传热量为：

其中：kf为液相导热系数；Nu采用Ranz-Marshall公式［8］计算。

1.4 相间质量传递-RPI壁面沸腾模型

采用RPI壁面沸腾模型计算壁面上的沸腾现象及气泡的冷凝等传质过程。RPI壁面沸腾模型包括基本热流密度分配模型和辅助模型。

1）壁面热流密度分配模型

RPI模型将壁面传递给液体的热流密度qw分为3个部分，即单相液体对流带走的热流密度qc、蒸发热流密度qe和淬灭热流密度qq。

式中：hc为单相对流换热系数；Af为壁面上液相润湿面积；Tw为壁面温度；Tf为液相温度；Vd为基于气泡脱离直径的气泡体积；Nw为气泡核化密度；hfg为汽化潜热；f为气泡脱离频率。

2）辅助模型

RPI模型的辅助模型包括液相润湿面积Af、气泡脱离频率模型f、气泡核化密度Nw、气泡脱离直径Dw等模型，计算关系式分别为：

式中，Ja为过冷度的Jacob数。

2 模型验证

Bartolemei等［9－10］对不同直径的竖直圆管内向上流动的过冷沸腾进行试验研究，试验获得不同压力、热流密度、进口质量流密度、进口过冷度下的空泡份额沿轴向分布以及壁面温度沿轴向分布。选取其中1个工况（d＝15.4 mm，G＝900 kg／（m2·s），q＝570 k W／m2，p＝4.5 MPa）用于模型的验证。图1a为根据试验段参数建立的几何模型及计算工况，图1b为该试验工况下Fluent模拟的沿轴向平均空泡份额、壁面平均温度和液相平均温度结果与试验测量数据的对比。从图中可看出，计算结果和试验数据符合良好，因此RPI模型可用于过冷沸腾的模拟。

图1 模型验证Fig.1 Model validation

3 几何模型及网格划分

本文选取的计算区域为带有格架和搅混翼、由4根1／4棒束围成的冷却剂通道。计算区域的几何模型及网格划分如图2所示。冷却剂从几何模型的底部流入，在燃料棒表面热流的加热作用下，从单相对流换热到发生过冷沸腾。几何模型整体结构简单，在定位格架和搅混翼局部位置几何复杂，因此针对该几何模型采用基于STAR-CCM＋的trimmer网格，整个几何域内为六面体网格，在格架和搅混翼几何较为复杂的区域进行网格劈分和加密。冷却剂通道的4个边界设置为两对周期性边界，允许流体在该边界上流入流出，来考虑临近流道对该流道的影响。

图2 简化几何模型和网格划分Fig.2 Simplified computational geometry and mesh

计算几何的进口边界设为速度进口边界，出口边界为压力出口边界，定位格架和搅混翼设置为无滑移壁面边界条件，燃料棒表面设为给定热流密度边界。具体边界条件设置列于表1。

表1 边界条件设置Table 1 Calculation condition

为分析网格数量对计算结果的影响，分别用5套不同数量的网格对计算工况进行分析计算，网格参数列于表2，采用轴向平均空泡份额作为网格无关性的考量依据。图3为沿轴向平均空泡份额的分布。可看出，第4套网格和第5套网格的结果基本一致，可满足网格无关性要求，因此以下的计算中均采用第4套网格。

表2 网格无关性分析Table 2 Grid independence analysis

图3 沿轴向平均空泡份额的分布Fig.3 Void fraction profile along flow direction

4 结果分析

计算工况中系统的压力选取一回路压力15.8 MPa，燃料棒束表面热流密度根据堆芯功率计算得到，冷却剂的物性查表得到。

图4为流道内的气液两相流线图。从图中可明显看出，在搅混翼上游的流道，流线排列整齐近乎平行分布，经过搅混翼的导向之后，明显观察到通道内出现了涡旋。而气相在搅混翼的下游会向流道中心集中，这是由于搅混翼的存在会在下游形成漩涡，由漩涡产生的离心力会使液相发生旋转，同时使气相积聚在漩涡中心，如图4b所示。

图4 气液两相流线图Fig.4 Streamline presentation of gas-liquid two-phase flow

图5为不同高度横截面上平均压力沿轴向变化曲线。从图中可看出，格架位置压降剧增，这是因为格架的存在会使该位置处的流通面积减小，从而导致很大的压降，而搅混翼的存在对压降影响不大。

以格架条带上表面搅混翼根部所在平面为基准面，图6示出搅混翼下游5、10、20、40 mm高度截面处的液相横向速度矢量图。从图中可看出，流体经过搅混翼的导流作用会在下游形成涡旋，燃料棒束通道的中心会形成涡旋，在搅混翼两侧也会形成与中心涡旋方向相反的两个小漩涡。这是由于搅混翼的存在，部分流体在搅混翼的导向下由周期性边界流出；另一部分流体在搅混翼的导向下流动接触燃料棒流向另一侧边界（图6a中箭头所示）。1对周期性边界上的参数相同，从而形成如图6a所示的涡旋。随着轴向高度的增大，涡旋的强度逐渐减小、涡旋影响面积逐渐增大。在搅混翼下游40 mm处，中心涡旋和两个方向相反的小涡旋仍存在，但强度明显减小。

图5 压力沿轴向变化曲线Fig.5 Pressure change along flow direction

图6 液相横向速度矢量图Fig.6 Liquid lateral velocity vector field distribution

图7为不同高度截面处液相温度的分布。5 mm高度截面上液相温度随着到燃料棒表面的距离增大而减小，距离燃料棒表面越近温度越高。受搅混翼的影响，该截面上的流体按图示方向流进或流出边界。如x正y负象限内的搅混翼会将流体导向x＋和y－方向，致使冷流体流向边界，将热流体推挤出边界。因此，会在4个周期性边界上形成如图所示的分布。在高度方向，随着高度的增加，流体不断被加热，截面上的平均温度不断上升。由图中明显看出，经过搅混翼后流体在流道中会发生旋转，使得温度场的分布也逐渐逆时针方向旋转。

图8为燃料棒束表面的壁面温度和热流密度分布。图中曲线1表示壁面温度随高度的变化，从图中可看出，壁面温度整体呈上升趋势，但在格架和搅混翼位置会出现明显的降低，这是由于格架和搅混翼的存在会对流场起到搅混作用，尤其是搅混翼会将冷流体导向燃料棒，使燃料棒表面温度出现明显降低。曲线2～5分别表示过冷沸腾RPI模型中壁面总热流、蒸发热流、单相对流热流以及淬火热流。过冷流体进入流道，随着过冷度的降低在某一高度处发生过冷沸腾，单相对流换热逐渐减小，蒸发热流逐渐增加。从图中可看出，蒸发热流曲线和单相对流曲线在格架和搅混翼处会分别出现降低和升高，这是由于流场被搅混后会使冷流体流向或直接冲刷壁面，该位置由过冷沸腾转变为单相强迫对流换热，壁面总热流不变的情况下，单相对流换热量上升，从而导致蒸发热流下降。

图7 液相温度的分布Fig.7 Distribution of liquid temperature

图8 壁面温度和热流密度分布Fig.8 Distributions of wall temperature and heat flux

图9为轴向5、10、20、40 mm高度截面处空泡份额分布。空泡份额随着与燃料棒表面距离增大而减小。由于搅混翼的存在，会将流道中心部分冷流体导向燃料棒表面，引起该处空泡份额降低（图9a中1位置）；搅混翼会将冷流体导向边界，而周期性边界的存在，致使该边界上有冷流体的流入（图9a中2位置）。空泡份额的分布随高度变化发生不同程度的旋转，搅混翼会拨动过冷流体流向燃料棒表面。20、40 mm截面上空泡份额的最大值出现在y方向的两个边界中心处，这正是两个小涡旋出现的位置，这是由于流道中该边界处的流通面积最小，该处的流体温度和空泡份额较高，而气泡在涡旋的作用下向涡旋中心聚集所致。

5 结论

本文采用经过试验验证的RPI过冷沸腾模型，对含定位格架和搅混翼的燃料棒束单通道内过冷沸腾现象进行了数值模拟及分析，得到以下结论：

1）通过与过冷沸腾基准试验Bartolemei试验数据进行对比，RPI模型可准确地预测过冷沸腾现象；

2）准确地模拟了含定位格架和搅混翼的单通道内的过冷沸腾现象，得到了流道内的流场、温度场以及空泡份额的分布；

3）定位格架的存在会造成很大压降，而搅混翼的存在会使流道形成涡旋，导致温度和空泡份额的分布沿流道发生旋转。

图9 空泡份额分布Fig.9 Distribution of vapor volume fraction

［1］ NAVARRO M A，SANTOS A A C.Evaluation of a numeric procedure for flow simulation of a 5×5 PWR rod bundle with a mixing vane spacer［J］.Progress in Nuclear Energy，2011，53（8）：1 190-1 196.

［2］ LIU C C，FERNG Y M，SHIH C K.CFD evaluation of turbulence models for flow simulation of the fuel rod bundle with a spacer assembly［J］.Applied Thermal Engineering，2012，40：389-396.

［3］ LI X C，GAO Y.Methods of simulating largescale rod bundle and application to a 17×17 fuel assembly with mixing vane spacer grid［J］.Nuclear Engineering and Design，2014，267：10-22.

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

［5］ ISHII M.Two-fluid model for two-phase flow［J］.Multiphase Science and Technology，1990，5（1-4）：1-63.

［6］ MORAGA F，BONETTO F，LAHEY R.Lateral forces on spheres in turbulent uniform shear flow［J］.International Journal of Multiphase Flow，1999，25（6-7）：1 321-1 372.

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

［8］ RANZ W，MARSHALL W.Evaporation from drops［J］.Chem Eng Prog，1952，48（3）：141-146.

［9］ BARTOLEMEI G G，BRANTOV V G，MOLOCHNIKOV Y S，et al.An experimental investigation of true volumetric vapor content with subcooled boiling in tubes［J］.Thermal Engineering，1982，29（3）：132-135.

［10］BARTOLEMEI 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.

Analysis of Subcooled Boiling in PWR Rod Bundle Channel

ZHANG Rui1，GAN Fu-jun2，ZUO Qiao-lin2，TIAN Wen-xi1，SU Guang-hui1，QIU Sui-zheng1
（1.State Key Laboratory of Multiphase Flow in Power Engineering，Xi’an Jiaotong University，Xi’an 710049，China；
2.Shanghai Nuclear Engineering Research ＆Design Institute，Shanghai 200233，China）

The RPI（Rensselaer Polytechnic Institute）wall boiling model，embedded in the context of two-fluid model of Fluent14.5，was applied to simulate subcooled boiling in core fuel bundle channel with mixing vanes and space grid.The subcooled boiling key parameters including flow field，temperature field and void fraction distribution were obtained.The analysis of the effect of mixing vanes and space grid on flow channel thermal-hydraulic characteristics was conducted.It is shown that the existence of space grid can cause a large pressure drop，and mixing vanes have a great influence on the key parameters of flow field and temperature field.The results of RPI model show a good agreement with Bartolemei experiment data.

subcooled boiling；rod bundle；space grid

TL334

：A

1000-6931（2015）09-1579-07

10.7538／yzk.2015.49.09.1579

2014-05-04；

2014-12-08

国家杰出青年科学基金资助项目（11125522）

张 蕊（1990—），女，陕西咸阳人，博士研究生，核能科学与工程专业