耦合一、二次侧换热的蒸汽发生器二次侧流场分析

丛腾龙，田文喜，秋穗正，苏光辉

(西安交通大学 动力工程多相流国家重点实验室 核科学与技术系，陕西 西安 710049)

蒸汽发生器(SG)作为压水堆的关键设备，同时承担着一、二回路侧换热以及一回路压力边界的作用，在反应堆运行过程中，需保持传热管的完整性；然而，在反应堆运行期间，传热管容易出现震动、疲劳以及由此带来的磨损甚至破裂。传热管的失效主要是由管束震动引起的磨损、管束热应力损坏以及传热管污垢沉积腐蚀造成的，对以上因素进行分析，需以SG二次侧的速度场、温度场作为输入参数。由于SG结构复杂、体积庞大，对SG二次侧流场进行实验和精细网格的CFD计算难以实现，因此本文采用多孔介质模型对SG二次侧流场进行数值分析。

目前已有一些针对SG二次侧流体热工水力特性的CFD研究，但这些研究仍有待改进。在Ferng等[1-3]的研究中，存在如下问题：1) 进出口边界设置不合理；2) 两相流计算关系式过于粗糙；3) 采用不正确的阻力和湍流计算模型；4) 计算忽略下降段、支承板和汽水分离器；5) 二次侧热源由RELAP5提供，而非耦合计算。这一系列简化或不合理假设，导致计算结果与实际运行情况差别很大，计算结果在工程应用上不具备参考价值。

基于以上分析，本文采用多孔介质方法，对压水堆SG二次侧流场进行数值模拟，同时耦合一、二次侧三维换热。在模拟中，综合考虑了管束、下降段、汽水分离器、支承板的阻力影响。对于汽水分离器，只考虑其阻力效应，不考虑汽水分离。

1 数学物理模型

常用的两相流计算模型包含均匀流模型、漂移流模型和两流体模型。均匀流模型将两相流体等效为单相流体，不能描述两相间的滑移速度等量。两流体模型可准确描述两相间的热力和水力学不平衡性，但由于管束外流动的相间作用力以及相间传热等模型目前还不完备，采用该模型需引入较多假设，导致两流体模型的计算精度不能保证。漂移流模型将两相流体视为混合相，但在相间引入漂移速度，考虑相间的水力学不平衡性，该模型虽从理论上相对于两流体模型精度较差，但由于该模型需引入的假设条件较少，计算精度可满足SG内多孔介质模拟的要求[4]。基于以上考虑，本文采用漂移流模型。多孔介质内的漂移流模型控制方程如下。

质量守恒方程：

β·(βρmvm)=0

动量守恒方程：

β·(βρmvmvm)=-β

能量守恒方程：

β·βαkvk(ρkEk+p))=

空泡份额方程：

β·(βαgρgvm)=

式中：下标m和g分别代表混合物和气相；β为孔隙率，即网格中流体体积与网格总体积之比；ρ为密度；vm为混合物速度；μm,eff为混合物有效黏度，等于混合物分子黏度和湍流黏度之和；αk为第k相体积份额；p为压力；S为附加源项，下标v、E、m分别代表动量、能量、质量；keff为有效热传导系数；vdr,k为第k相的漂移速度，漂移速度与相分布计算关系式参见文献[5]。

动量、能量、空泡份额方程中的源项定义如下。

对于动量方程，其源项表示由于控制体内蒸汽发生器构件的存在导致的附加阻力，包括下降段阻力、管束阻力、汽水分离器阻力、支承板阻力。其中，前3个阻力均为分布阻力，以分布阻力的形式添加到控制体网格，支承板阻力为局部阻力，以多孔阶跃模型的方式添加到支承板所在位置的网格界面上。对顺流和横流阻力分别采用管束外顺流[6]和横流[7]阻力经验关系式计算。支承板、下降段和汽水分离器的阻力，根据蒸汽发生器设计参数给定[4]。

对于能量方程，其源项为一次侧向二次侧释热。将一次侧流场划分为与二次侧相同的网格，对于每个网格，每次迭代开始之前，在UDF中分别计算网格内一次侧向二次侧的释热量。U型管内为单相液体对流，采用Dittus-Boelter公式[8]计算；管外换热包括单相对流换热、过冷沸腾换热以及饱和沸腾换热，对于单相对流换热，将流速分解为顺流流速和横流流速，分别计算顺流换热系数[9-10]和横流换热系数[11]，将这两者的加权和[4]作为实际的斜流换热系数。当主流温度达到饱和温度时，发生饱和沸腾；当管外壁温度大于饱和温度且主流温度低于饱和温度时，有可能发生过冷沸腾，本文在该工况下分别按照单相流动和过冷沸腾计算壁面热流密度，若采用过冷沸腾公式计算得到的热流密度大于采用单相流动的计算值，认为过冷沸腾发生。

在计算中，给定一、二次侧入口的质量流速和温度，即可计算一、二次侧之间的换热。

对于空泡份额方程，其源项表示气相质量源项，其值等于蒸发率与冷凝率之差。动量方程中的漂移速度采用Lellouche-Zolotar公式[5]计算。

对于湍流模型，已公开发表的模型多是针对纯流体区域的精细网格提出的[12]，对于多孔介质内的湍流模型，研究者在多孔介质粗网格内对已有湍流模型(如k-ε)进行积分，得到适用于多孔介质的湍流模型[13-15]，但这些模型仅用于单相流动，并无针对两相流动的湍流模型。本文采用文献[16-17]提出的针对管束间两相流动的零方程湍流模型计算湍流黏度及有效黏度，该模型在ATHOS程序中已被验证[4,18]。

以上模型构成完备的控制方程组，将这些方程在ANSYS FLUENT求解器中求解，求解方法为COUPLE算法。

2 网格划分及边界条件

本文所选取的二次侧计算区域为：从下降段入口向下到围板缺口，然后折流进入管束区，最终进入主汽水分离器。计算区域示意图如图1a所示。对于圆筒形的汽水分离器筒体，为简化网格划分，采用方形通道代替，保持筒体流通面积和中心位置不变。为了提高计算的收敛性并模拟汽水分离器的阻力对上游流场的影响，将汽水分离器添加到计算域中，汽水分离器只考虑其阻力，不考虑汽水分离效果。对计算区域划分网格，所有网格均为六面体网格，在支承板位置添加多孔阶跃边界模拟支承板阻力。网格划分如图1b所示，网格质量大于0.5。图1c为计算区域管束布置示意图。为进行网格独立性验证，分别对4组网格进行计算，获得独立解时的网格数量为111 099，如图2所示。

a——计算区域示意图；b——网格划分示意图；c——管束示意图

图2 网格独立性验证

二次侧网格入口边界设置为速度入口，出口边界设置为压力出口，考虑到多孔介质网格尺度较大，壁面无滑移边界不适用，因此本文采用滑移壁面，将壁面摩擦阻力分布添加到壁面网格内。入口温度由疏水和新给水的流量及焓确定，出口压力设置为蒸汽发生器顶部腔室压力。

一次侧流场和一、二次侧流体的换热均采用UDF进行计算，二次侧流场计算每迭代一步，UDF计算1次一次侧流场和一、二次侧换热量，将换热量作为二次侧流场内热源赋给FLUENT求解器。

3 结果分析

给定一、二次侧入口流量和温度，采用FLUENT计算二次侧流场，同时采用UDF耦合一次侧流场及一、二次侧换热，得到三维的一、二次侧流场，管内外换热系数及U型管温度分布。FLUENT计算值与设计值的对比列于表1。

表1 FLUENT计算值与设计值的对比

3.1 一、二次侧温度分析

图3为SG二次侧流体在对称截面和水平切面上的温度云图，图3a中右侧(x>0)为热侧，左侧为冷侧。由图可看出，热侧流体在进入管束区后，很快达到饱和温度，而冷侧流体缓慢达到饱和，冷侧流体流过蒸汽发生器内1/3高度时才完全达到饱和。达到饱和后，流体温度保持不变。由图3b所示的横切面温度分布可看出，由于流动和几何的对称性，流体温度在截面上对称。

a——对称截面上温度分布；b——不同高度水平截面上的温度分布

图4为径向积分后温度沿管束方向的分布。一次侧温度沿着流动的进行逐渐降低。二次侧温度在冷热侧均随高度的增大而增大，并最终维持在饱和温度。热侧流体上升到约1 m高度时，即达到饱和温度；冷侧流体在上升到约4 m高度时才达到饱和温度。由于二次侧换热系数远大于一次侧换热系数，管外壁温度沿着管束方向基本不变，维持在547 K。管内壁温度沿着管束方向自热侧向冷侧逐渐降低，管内外壁温差沿管束方向逐渐减小，管内壁最高温度约为570.6 K。

图4 温度沿管束方向的分布

图5为一、二次侧换热系数及壁面热流密度沿管束的分布。管内为单相强制对流，换热系数随物性的变化有轻微下降趋势；管外分别为单相对流、过冷沸腾和饱和沸腾，其中沸腾换热系数远大于单相换热系数。一次侧平均换热系数为15 856.5 W/(m2·K)，二次侧平均换热系数为63 623.0 W/(m2·K)，二次侧最大换热系数为122 862.9 W/(m2·K)。由于换热系数和管内外温差的降低，管壁面热流密度沿管束方向也逐渐降低，外壁面平均热流密度为149.9 kW/m2，最大、最小壁面热流密度分别为369.6和44.6 kW/m2。

图5 一、二次侧换热系数及壁面热流密度沿管束方向的分布

3.2 二次侧流动含汽率分析

SG二次侧流体的流动含汽率分布如图6所示。可看出，在管束区入口区域，冷热侧流动含汽率差别较大，热侧流动含汽率小于0.01的区域很小，而冷侧该区域占有较大面积。随着流体沿管束上升，流体不断蒸发，气泡在横向扩散，冷热侧流动含汽率逐渐均衡，且热侧靠近中心处流动含汽率逐渐变为最大。

a——对称截面上的流动含汽率分布；b——不同高度水平截面上的流动含汽率分布

图7为流动含汽率沿高度方向的分布，3条曲线分别表示整个二次侧区域水平截面平均值、热侧平均值、冷侧平均值。由图7可看出，随高度的增加，流动含汽率逐渐增大。热侧流动含汽率远大于冷侧流体。冷、热侧出口平均流动含汽率分别为0.16和0.46，总体出口平均流动含汽率为0.27。在流体进入汽水分离器前，流动含汽率出现降低，这是由于在两相流体流出管束区、流入汽水分离器前，流通截面积突变，导致气相流速降低，从而导致流动含汽率降低。在流体流入汽水分离器之后，流动含汽率保持恒定。

图7 流动含汽率沿高度方向的分布

图8 汽水分离器内流动含汽率分布

图8为汽水分离器内流动含汽率的分布。由图8可看出，对于热侧最内层汽水分离器，其流动含汽率很高，达0.75，而对于冷侧外围的汽水分离器，流动含汽率只有0.07，不同的汽水分离器内混合物流动情况相差很大。

3.3 二次侧压力分析

图9为对称截面上的压力分布，压力选取出口截面为参考面，该截面压力设为0 Pa，图中标注压力为相对压力。在下降段，由于重位压降、下降腔室壁面摩擦压降和加速压降的共同作用，压力沿着下降段向下逐渐增大；流体在围板缺口区域，由于流体转向，压力陡增，之后进入管束区域，流体压力随着高度方向逐渐减小。同一高度截面上，热侧压力总体稍大于冷侧压力。图10为压力沿高度方向的分布。压力沿着高度方向出现明显的阶跃变化，这是由于支承板的存在，引入局部阻力，导致流体经过支承板时，压力突降；支承板导致压降的大小随轴向高度的增大而增大，这是由于随高度的增加，两相流体的流动含汽率增大，两相摩擦压降倍增因子增大，从而在相同的阻力系数下，支承板位置越高，引入的压降越大。位于11～12 m之间的压力阶跃是由于流体流入汽水分离器造成的，流体流通截面积突降，流速增大，导致静压降低。

图9 压力在对称截面上的分布

图10 压力沿高度方向的分布

3.4 二次侧速度分析

图11为对称截面上的速度分布。由图11可看出，热侧流体速度明显大于冷侧，最大速度约7.4 m/s。还可看出，除围板缺口折流部分和弯管区域，流体基本沿管束方向流动，横流较小。

图11 对称截面上的速度矢量图

图12 弯管部分横流速度分布及横流能量分布

图12为弯管部分横流速度分布及横流能量分布。由图12可看出，自冷侧到热侧，横流速度整体呈现先增大后减小的趋势，最大横流速度约为4.06 m/s；在弯管部分顶端，由于汽水分离器的影响，导致汽水分离器控制部分的流速高于弯管最顶端的速度，这与Ferng等[2-3]的不考虑汽水分离器的速度分布不同。根据文献[20]，管束流致振动强度与横流能量(ρu2)呈正比，由于冷侧流体密度较大，冷侧横流能量高于热侧，冷热侧最大横流能量分别为1 145和695 J/m3。

4 结语

本文采用多孔介质模型，对压水堆蒸汽发生器二次侧三维两相流场进行数值模拟，同时耦合一、二次侧三维换热，得到如下结论：

1) 蒸汽发生器内冷热侧流体由于受热不均，导致两侧流体流动状态差别较大；

2) 一、二次侧平均换热系数分别为15 856.5和63 623.0 W/(m2·K)，二次侧最大换热系数为122 862.9 W/(m2·K)，U型管外壁面平均热流密度为149.9 kW/m2；

3) 汽水分离器内流动含汽率分布极不均匀，其值介于0.07～0.75之间，此分布可为汽水分离器设计提供入口边界条件；

4) U型管弯管段最大横流速度约为4.06 m/s；冷侧冲刷U型管的横流能量大于热侧，其最大值分别为1 145和695 J/m3。汽水分离器的位置对U型管横流速度以及横流能量的影响很大。

参考文献：

[1] FERNG Y M. Investigating the distribution characteristics of boiling flow and released nuclide in the steam generator secondary side using CFD methodology[J]. Ann Nucl Energy, 2007, 34(9): 724-731.

[2] FERNG Y M, CHANG H J. CFD investigating the impacts of changing operating conditions on the thermal-hydraulic characteristics in a steam generator[J]. Appl Therm Eng, 2008, 28(5): 414-422.

[3] FERNG Y M, YINPANG M, KANG J C. Thermal-hydraulic simulation of localized flow characteristics in a steam generator[J]. Nucl Technol, 2001, 136(2): 186-196.

[4] KEETON L, SINGHAL A, SRIKANTIAH G. ATHOS3: A computer program for thermal-hydraulic analysis of steam generators[M]. Palo Alto, CA, US: Electric Power Research Institute, 1986.

[5] LELLOUCHE G, ZOLOTAR B. A mechanistic model for predicting two-phase void for water in vertical tubes, channels and rod bundles, EPRI-N-2246[R]. Palo Alto, CA, US: Electric Power Research Institute, 1982.

[6] MacADAMS W H. Heat transmission[M]. New York, US: McGraw Hill, 1954.

[7] GRIMISON E. Correlation and utilization of new data on flow resistance and heat transfer for cross flow of gases over tube banks[J]. Trans ASME, 1937, 59(7): 583-594.

[8] DITTUS F W, BOELTER L M K. Heat transfer in automobile radiators of the tubular type[J]. Int Commun Heat Mass Transfer, 1985, 12(1): 3-22.

[9] DINGEE D A, BELL W, CHASTAIN J, et al. Heat transfer from parallel rods in axial flow[R]. Columbus, Ohio: Battelle Memorial Inst., 1955.

[10] DINGEE D A, CHASTAIN J W. Heat transfer from parallel rods in an axial flow[C]∥ASME Reactor Heat Transfer Conference. Oak Ridge: US Atomic Energy Commission, 1968.

[11] DWYER O, SHEEHAN T, WEISMAN J, et al. Cross flow of water through a tube bank at reynolds numbers up to a million[J]. Industrial & Engineering Chemistry, 1956, 48(10): 1 836-1 846.

[12] 陶文铨. 数值传热学[M]. 第2版. 西安：西安交通大学出版社，2001.

[13] ANTOHE B, LAGE J. A general two-equation macroscopic turbulence model for incompressible flow in porous media[J]. Int J Heat Mass Transfer, 1997, 40(13): 3 013-3 024.

[14] CHANDESRIS M, SERRE G, SAGAUT P. A macroscopic turbulence model for flow in porous media suited for channel, pipe and rod bundle flows[J]. Int J Heat Mass Transfer, 2006, 49(15): 2 739-2 750.

[15] De LEMOS M J S. Turbulence in porous media: Modeling and applications[M]. London: Elsevier, 2012.

[16] LAUNDER B E, SPALDING D B. Lectures in mathematical models of turbulence[M]. London: Academic Press, 1972.

[17] SCHLICHTING H. Boundary-layer theory[M]. 7th ed. Germany: Springer, 1979.

[18] HOPKINS G. Verification of the ATHOS3 code against feedring and preheat steam generator test data[M]. Palo Alto, CA, US: Electric Power Research Institute, 1988.

[19] 林诚格，郁祖盛，欧阳予. 非能动安全先进核电厂AP1000[M]. 北京：原子能出版社，2008.

[20] AXISA F, ANTUNES J, VILLARD B. Overview of numerical methods for predicting flow-induced vibration[J]. J Pressure Vessel Technol, 1988, 110(1): 6-14.