压水堆蒸汽发生器三维两相热工水力分析程序开发与验证

(上海核工程研究设计院有限公司，上海 200233)

0 引言

蒸汽发生器作为核电厂一、二回路的枢纽，是核岛最重要的设备之一。其是否安全运行直接关系到整个核电厂的安全。蒸汽发生器二次侧流场数据是蒸汽发生器设计和传热管流致振动和磨损分析的重要输入。由于蒸汽发生器结构复杂、体积庞大，采用试验手段获得二次侧流场数据代价太大，难以实现；同样由于结构复杂，直接针对管束建模分析很难实现。因此本文将二次侧管束区作为多孔介质，引入分布阻力概念，对蒸汽发生器二次侧三维两相热工水力分析程序进行开发，并开展小规模管束的流动换热试验，以完成对程序的试验验证。

最早将多孔介质模型应用于换热器数值模拟分析的是Patankar等[1]，他们于1974年将管壳式换热器的壳侧作为多孔介质，完成了对管壳式换热器壳侧流场的数值模拟。在此基础上，Sha 等[2]模拟了SG和反应堆堆芯中冷却剂的流动；Prithiviraj等[3-4]模拟了三维换热器中的流动。Ferng等[5-7]采用CFX 软件多孔介质模型对SG管束区的二次侧流场进行了分析。1977年美国电力研究院(EPRI)开始开发蒸汽发生器热工水力计算程序，后来通过试验验证和实堆反馈数据的不断修正，于1984年正式推出了适用于U形管蒸汽发生器(UTSG)和直流蒸汽发生器(OTSG)二次侧三维流场分析的CFD程序ATHOS[8]。Chan等[9]于1986年开发了UTSG稳态、瞬态三维流场分析程序PORTHOS。法国电力公司(EDF)为研究PWR和LMFBR堆芯和组件内的单相、两相流场而开发了THYC程序[10]，之后推出适用于传热器流场分析的版本THYC-EXCHANGER[11]来分析蒸汽发生器二次侧管束区两相流场。法国原子能委员会(CEA)开发并验证了UTSG三维热工水力分析程序GENEPI[12]。可以看出，国际上已有多个蒸汽发生器三维热工水力分析程序，这些软件均基于多孔介质模型，只是其两相流控制方程、管内外传热、相变及流动阻力计算方法有所不同。国内也有学者[13-16]采用CFD软件对蒸汽发生器二次侧进行三维热工水力分析，大多属于方法性的探索研究，未形成专用分析程序。因此，本文基于已有的理论研究基础，弥补现有分析程序的缺陷，开发可用于压水堆蒸汽发生器三维两相热工水力分析程序。

1 程序开发

与常规流体的两相流动类似，宏观尺度的多孔介质内的两相流动研究方法也分为三种：均相流模型、漂移流模型和两流体模型。在早期的研究工作中，为简化计算，多采用均相流模型，如CALIPSOS[17-18]、THEDA-1[19]、THEDA-2[20-21]等程序以及Ferng等[22-23]的一系列工作。均相流模型将两相流体等效为单相流体，不能描述两相间的相对速度，因此在计算竖直通道内两相流动的时候，会出现空泡份额偏大的问题。从理论上来说，两流体模型可以准确地描述两相间的热力学和水力学不平衡特性，但由于管束外流动的相间作用力及相间传热等模型目前尚不完备，采用该模型需引入较多假设，导致两流体模型的计算精度不能保证。漂移流模型同样将两相流体视为混合相，但是为描述相间的速度差，在均相方程的基础上引入相间滑移速度。虽然漂移流模型从理论上相对于两流体模型精度较差，但由于该模型需引入的假设条件较少，计算精度可满足SG内模拟的要求[8]。因此本程序开发中采用漂移流模型。多孔介质内的漂移流模型控制方程如下。

质量守恒方程：

(1)

式中t——时间，s；

β——孔隙率；

ρm——混合物密度，kg/m3,ρm=αgρg+αlρl；

动量守恒方程：

(2)

式中p——压力，Pa；

μm,eff——混合物有效黏度，Pa·s,

μm=αgμg+αlμl;

αg——汽相体积份额；

ρg,ρl——汽相、液相密度，kg/m3；

αl——液相体积份额；

能量守恒方程：

(3)

式中Hg——汽相焓，J/kg；

Hl——液相焓，J/kg；

km——混合物导热系数，W/(m·K),km=αgkg+αlkl;

SE——能量源项,W/m3。

空泡份额输运方程：

(4)

式中Sg——汽相质量源项,kg/(m3·s)。

程序动量方程源项中考虑的阻力主要由以下部件引入，如传热管束、汽水分离器、下降通道、支承板等，其中管束、汽水分离器和下降通道的阻力以分布阻力形式添加到这些部件所在网格中，支承板阻力以集中阻力形式添加到支承板所在位置的网格界面上。根据流动方向不同，将流动方向沿管束轴向和横向分解，在轴向和横向分别引入动量源项。对于顺流和横掠阻力分别采用管束外顺流和横掠[24-25]阻力经验关系式计算。支承板、下降通道和汽水分离器的阻力系数，根据蒸汽发生器设计结构直接给出。

能量方程的源项，即一次侧向二次侧释热，被简化为一维分布计算。程序采取一维-三维耦合传热计算方式。该方法一、二次侧网格示意如图1所示，计算思路如下：将一次侧沿流体流动方向划分为J个控制体，局部网格编号用j表示；将二次侧划分为K个控制体，局部网格编号为k；通过DEFINE_INIT宏中的几何处理模块、根据网格位置将一、二次侧网格进行匹配。

一次侧为单相液体对流换热，采用Dittus-Boelter 公式[26]计算；二次侧换热包括单相对流换热、过冷沸腾换热及饱和沸腾换热，对于单相对流换热，由于目前的经验关系式多为管外顺流和横掠关系式，无倾斜冲刷管束的换热关系式，因此将流速分解为顺流和横掠流速，分别计算顺流换热系数[24-25]和横掠换热系数[27]，并将两者的加权和[8]作为实际的斜掠换热系数。在计算中，给定一次侧入口的质量流速和温度，便可计算一、二次侧间的流动换热。

图1 程序的一、二次侧网格示意

考虑到多孔介质模型对局部参数的平均效应，以及漂移流模型对相界面密度的依赖性较低，因此，程序中对沸腾模型进行简化，仅根据能量平衡计算两相流温度及蒸发率，不计算相界面输运过程。引入热量密度二分方法，将总热流密度分为蒸发热流密度和对流热流密度，分别用于蒸发液相和加热混合相。

空泡份额方程中的源项表示气相质量源项，其值等于蒸发率与冷凝率之差，对于稳态计算忽略SG内的非稳态冷凝现象，气相质量源项等于液相蒸发率。

对于湍流模型，已公开发表的文献多是针对纯流体区域的精细网格提出[28]，Teruel等[29-31]从宏观尺度研究多孔介质内的湍流流动，他们在宏观网格内对k-e方程进行积分，得到宏观k-e二方程湍流模型，并且分别采用微观湍流模型和宏观湍流模型计算几何规则的多孔介质内的流动，从而验证这些宏观k-e湍流模型的正确性。但是这些模型均针对多孔介质内的单相流体进行开发，并不适用于多孔介质内的两相流动。因此，本程序开发中采用了商用蒸汽发生器热工水力分析程序ATHOS[32]和GENEPI[33]中采用的零方程湍流模型计算两相混合物的有效黏度。

针对压水堆蒸汽发生器结构特点，以CFD软件多孔介质模型为基础，采用UDF将以上传热、阻力模型添加到求解器中，并通过UDF求解一次侧流场及一、二次侧耦合传热特性，完成蒸汽发生器二次侧三维两相热工水力分析程序的开发。

2 程序验证

在本文开发的程序用于工程之前，需对其进行试验验证。由于蒸汽发生器结构庞大，对原型蒸汽发生器进行试验研究需要极大的蒸汽流量和加热功率，成本很高。因此，基于模化分析，采用小规模管束进行流动传热试验，模拟蒸汽发生器一、二次侧流场，获得试验件一、二次侧流场信息，同时采用该程序对试验件建模分析，预测试验工况下的试验件的一、二次流场数据，通过与试验测量结果的比较，达到对程序验证的目的。

2.1 蒸汽发生器二次侧流动传热试验

试验中采用40根传热管模拟二次侧管束的流动换热状态，管束横截面如图2所示。为准确体现蒸汽发生器原型的热工水力参数，试验件管束高度上与原型保持相等，试验模拟体中采用的传热管材料、管径、管间距、布置方式、管束长度和支承板结构与原型相同。试验中一、二次侧进口温度、系统运行压力和循环倍率等参数与原型相同，分别为：传热管规格为∅17.48 mm×1.01 mm，直段高度9.2 m，相邻传热管的中心距24.89 mm，呈三角形排列，一次侧压力15.5 MPa，进口温度322 ℃；二次侧压力5.6 MPa，给水温度226.7 ℃。试验模拟体如图3所示，包括上下筒体、上下封头、管板、管束、下降管、管束支承板、汽水分离装置等部件。为了避免测量探头同时穿过筒体和套筒两层壁面，采用两根圆形下降管代替蒸汽发生器环形下降通道，同时在下降管上安装了调节球阀，调节下降管阻力，以便更好地控制试验模拟体运动的循环倍率，下降管与管束区入口的连接如图4所示。试验体运行及设计参数如表1,2所示。

图2 管束横截面示意

图3 试验模拟体示意

图4 下降管与管束区入口连接示意

表1 试验体额定运行参数

表2 试验体设计参数

试验中测量的参数包括流量、温度、压力、差压、空泡份额、流体速度等。一次侧进、出口温度通过在一次侧进、出口管道上安装铠装热电偶测量；一次侧流量通过在一次侧进口管道上安装孔板流量计测量。二次侧出口蒸汽压力在蒸发器出口管上利用压力传感器测量，采用在蒸汽出口管道上安装的孔板流量计测量蒸汽出口流量，给水流量测量与蒸汽流量测量方法一致，在给水管道上安装铠装热电偶测量给水温度。采用光纤探针和伽马射线法测量空泡份额。通过测量得到如下宏观参数：热功率、一回路压力、一回路进、出口温度、给水流量、给水温度、给水压力、蒸汽流量、蒸汽温度、蒸汽压力、二次侧水位等。下降通道上的流量采用超声波流量计进行测量。在下降段底部布置2个热电偶用于测量进入管束区的流体温度。

对于一次侧沿程流体温度，采用K型铠装热电偶，热电偶采用管接头方式从传热管内部引出，再采用管接头将其引出筒体。安装方式如图5所示。管接头作为铠装热电偶穿过压力边界的密封元件，在与传热管焊接安装前进行焊接工艺试验，验证了密封可靠性。

图5 一次侧测温热电偶管接头安装方式

空泡份额的测量分别采用光纤探针法、伽马射线法和压差法获得局部点、截面和体积含汽率。

光纤探针法是利用光学的全反射原理，传输光在光纤探针敏感头界面与媒质界面处符合全反射条件时，这一部分光将返回敏感头内部，如果敏感头接触的媒质折射率产生变化时，由于可能将不满足全反射条件，此时不能产生全反射现象。通过检测经敏感头反射回接收端的光强的强弱变化而产生的高低不同的电信号来分辨光纤探针敏感头测量到的是气相还是液相，从而根据产生的连续变化的电压信号测量局部截面持气率，光纤探针敏感头测量原理图如图6所示。试验中测点位置见图7。具体测点位置和探针数量见表3。

射线法是一种非接触式的测量方法，其原理是根据汽、液两种物质对射线的吸收程度不同推算出射线照射区域的含汽率。该方法不会破坏管道中的流场和温度场的自然分布，适用于高温高压条件下测量，原理如图8所示。试验测点位置如图9所示。

图6 光纤探针法测量含汽率的原理示意

图7 空泡份额测点布置

表3 光纤探针空泡份额测点位置

图8 射线法原理示意

图9 射线法测点位置示意

压差法则是通过测量两个取压截面的压力差，计算两个截面之间的流体体积含汽率。

为了获得较多的试验数据，试验工况围绕额度工况，分别改变一次侧进出口温度、给水温度、加热功率和循环倍率等，完成了11个试验工况。

试验测量的一、二次侧温度分布如图10所示。三种方法空泡份额测量结果见表4。

2.2 试验件模拟分析

采用本文开发的程序对试验件进行数值模拟分析，预测二次流场数据。计算中选取试验件管板以上至分离器入口为分析对象。对计算区域进行建模，管束区采用多孔介质模型进行模拟，为了方便设置支承板阻力系数，支承板阻力采用多孔阶跃边界描述，在边界上设定支承板的阻力系数参数和厚度。为了保证数值计算稳定性与准确性，采用一体化网格，对整个计算区域进行建模，然后划分网格，所划分网格全为六面体网格，网格总数为196 107，计算模型网格如图11所示。

(a)一次侧温度分布

(b)二次侧温度分布

表4 空泡份额测量结果对比

根据试验件和试验工况，计算中入口为下部开口，出口为分离器入口，入口边界设为速度入口，不同工况给出不同入口流速，出口边界设为压力出口。计算中采用COUPLE算法求解控制方程组。本分析采用UDF输入自定义模块，如附加源项、各附加阻力项，因此计算过程中，除求解控制方程，还需调用UDF，在User-Defined Profile定义计算区域孔隙率；在User-Defined Init中定义热源分布、网格润湿面积，以及其他需要在计算进行之前计算的变量；在User-Defined Adjust中计算求解阻力源项所需的各个变量，如各控制体的阻力系数。

计算结果如图12所示。为了图片更加清晰，图片在高度方向缩短了15倍。

图11 计算模型网格

图12 工况1的二次侧温度云图及气相体积份额云图

2.3 程序对比验证

为了验证程序分析的准确性，将相同工况下程序预测结果与试验测量结果进行比较。二次侧温度进行了归一化处理，用测量点的温度除以该工况下的饱和温度，热侧和冷侧归一化比较结果如图13所示。空泡份额比较如图14所示。

(a)热侧

(b)冷侧

(a)射线法测量的空泡份额

(b)探针法测量的二次侧空泡份额

3 结论

本文针对压水堆核电厂蒸汽发生器结构特点，以CFD软件中多孔介质模型为基础，采用UDF对CFD软件进行了二次开发，完成了蒸汽发生器二次侧三维两相热工水力分析程序的开发。并基于模化分析，采用小规模管束进行流动传热试验，模拟蒸汽发生器一、二次侧流场，获得试验体试验工况下的一、二次侧流场数据，采用该程序对试验件建模分析。结果表明：程序能够实现对蒸汽发生器试验件的模拟，预测结果与试验测量结果符合良好，证明了该程序计算结果的可靠性。

未来蒸汽发生器二次侧三维两相热工水力分析程序还需开展进一步工作。经后续不断改进完善，程序将实现在工程上的应用，该程序可成为我国核电厂蒸汽发生器设计研发的重要工具。