葛 魁，王 辉，王明军,*，田文喜，苏光辉，秋穗正

(1.西安交通大学 核科学与技术学院，陕西 西安 710049；2.西安交通大学 动力工程多相流国家重点实验室，陕西 西安 710049；3.中国核电工程有限公司，北京 100840)

安全壳是事故工况下防止反应堆内的放射性物质向外泄漏的最后一道屏障，其完整性对核反应堆安全至关重要。第三代核电技术[1-2]因具有很高的安全性和经济性逐渐受到全世界核电行业的认可。第三代核电技术广泛采取了非能动与能动相结合的安全设计理念[3]，满足了“国际最高核电安全标准”的要求[4-5]，可以靠核电厂自身的安全系统来保证安全性，减少发生事故后人为操作不当对核电站安全性的影响。除普通压水堆[6]外，非能动安全系统还被应用于一体化轻水堆[7]中。

黄政[8]基于均相流模型建立了一维非能动安全壳热量导出系统分析程序，利用牛顿迭代法求解，模拟了稳态运行及在事故工况下安全壳和非能动安全壳热量导出系统的瞬态响应过程，得到自然循环系统的流动换热特性。李军等[9]基于计算流体力学(CFD)软件对华龙一号循环水箱升温过程进行了三维流动传热的数值模拟。Guo等[10]基于两相均匀流模型开发了模拟PCS自然循环的瞬态计算程序，并对PCS的循环流量、压降、温度及传热系数等热工水力参数进行分析研究。白晋华等[11]基于Relap5程序从启动时间、运行工况的稳定性等方面对多种PCS设计方案进行了评价。

综上，目前针对PCS热工水力分析大多基于简单的均相流模型，或基于商用Relap5程序进行分析计算。由于均相流模型本身限制，某些工况下的计算结果不能真实反映系统特性及状态参量的变化，而漂移流模型能更好地反映系统特性[12]。因此，本文针对华龙一号PCS，基于两相漂移流模型自主可控开发PCS一维自然循环瞬态计算程序，并利用该程序对PCS的热工水力特性进行分析研究。

1 数学物理模型

华龙一号换热水箱中温度较低流体由于重力作用沿下降管路流入换热器，流体流经换热器时被加热、膨胀，之后流入上升管路最终流回至换热水箱。随着流入换热水箱的流体温度不断升高，当水箱中的流体温度达到饱和时，将有部分水蒸气排出至外界环境，换热水箱内水位降低。当水位降低到一定程度时，水箱补水功能启动，补水管路开始以设定的流速向水箱中注入低温的水，直至水箱水位恢复至原来高度。程序中假设工质在系统内的流动是一维的，即工质的热工水力参数仅沿轴向变化，且忽略了换热器传热管壁沿轴向的导热。

1.1 控制方程

PCS内流体热工水力参数的变化遵循质量、动量和能量守恒[13]的基本规律，再辅以流体的状态方程及相应的辅助方程可构成封闭的方程组，采用Gear算法求解，就可得到上述的状态参数。为能更为准确地模拟两相过程中的瞬态特性，本研究采用漂移流模型[13]。相较于其他模型，漂移流模型采用代表两相介质横向分布的量C0和代表两相之间局部相对速度的量Vgj来描述两相流动特性，且动量方程中的压降除考虑摩擦阻力压降、重力压降和加速压降外，还要考虑漂移流压降。

质量守恒方程为：

(1)

式中：t为时间，s；α为空泡份额；ρf、ρg分别为液相和气相的密度，kg·m-3；A为流通面积，m2；W为流量，kg·s-1；Z为流动方向长度，m。

动量守恒方程为：

(2)

式中：p为压力，MPa；τ为剪切力，N·m-2；Uh为加热周长，m；ρ为流体密度，kg·m-3；g为重力加速度，m·s-2；G为质量流速，kg·m-2·s-1；x为含气率。

在漂移流模型中，动量方程表示的压降除摩擦阻力压降、重力压降和加速压降外，还要考虑漂移流压降梯度：

(3)

(4)

Vgj和C0具体计算公式如下：

(5)

(6)

(7)

Ed=1-exp(0.138ρf/ρg-0.23α/μg)

(8)

(gDLPL(ρf/ρg-1))0.5

(9)

(10)

(11)

式中：σ为表面张力系数，N·m-1；μg和μf分别为气相、液相的动力黏度，N·s·m-2；De为管道的水力直径，m。

能量守恒方程为：

(12)

式中：hf、hg分别为液相和气相的焓值，kJ·kg-1；Vf为液相的流速，m·s-1；Vg为气相的流速，m·s-1；q为热流密度，W·m-2。

自然循环质量流量方程为：

(13)

式中：B为自然循环驱动压头，MPa；Δpf为回路摩擦压降，MPa；Δploc为回路局部阻力压降，MPa；ΔpDG为漂移流压降，MPa；SLA为回路惯量；W为PCS自然循环质量流量，kg·s-1。

换热器由多根传热管组成，程序中将自然循环流量平均分配到各传热管中，且自然循环流量计算中考虑下降管路至多根传热管之间的形阻压降。

1.2 辅助模型

1) 传热模型

在PCS中，安全壳内流体与换热器传热管外壁面、换热器传热管内的流体与传热管内壁面的换热，均采用牛顿冷却定律[14-17]：

Q=h·A·ΔT

(14)

式中：Q为热流量，W；ΔT为传热管内壁面与传热管内流体的温差，K。

安全壳内流体与换热器传热管外壁面通过含不凝结气体的蒸汽冷凝过程进行传热，冷凝换热系数采用Dehbi关系式[18-19]进行计算：

hDehbi=1.25L0.05((3.7+28.7pb)-

(2 443+458.3pb)lgWa)(T1-Tw)-0.25

(15)

式中：hDehbi为冷凝换热系数，W/(m2·K)；L为传热管管长，m；pb为安全壳压力，atm；Wa为安全壳内混合气体中空气质量分数；Tl为安全壳流体温度，K；Tw为传热管壁面温度，K。公式适用范围为：压力0.15～0.45 MPa，管长0.3～3.5 m，壁面过冷度10～50 ℃。

换热器传热管内的流体与传热管内壁面的换热根据不同情况选取不同的换热系数关系式。对于单相流动换热系数采用D-B公式计算：

Nu=0.023Re0.8Prn

(16)

式中，加热流体时n=0.4，冷却流体时n=0.3。

对于两相流动换热采用詹斯-洛特斯(Jens-Lottes)公式计算：

(17)

式中：p为压力，MPa；q为热流密度，W/m2。

高温流体把热量传给换热器传热管，同时传热管内的冷流体在管内流动，带走一部分能量，使换热器传热管冷却。换热器传热管管壁导热方程为：

α2F2(Tw-T2)

(18)

PCS中冷却剂在换热器吸收安全壳内热量后流进换热水箱，与换热水箱中流体混合，导致换热水箱内水温升高。当换热水箱的水温达到饱和时，将有部分水蒸气排出至外界环境，使换热水箱水位降低。

2) 阻力模型

单相流动的摩擦阻力系数f[10]为：

(19)

两相流动摩擦压降的计算要比单相的复杂得多，直接计算有一定难度。通常采用先计算全液相压降梯度，然后再乘以一个因子[13]：

(20)

两相摩擦压降倍增因子选用McAdams关系式[20]：

(21)

式中，vg、vf分别为气相和液相的比体积，m3/kg。

程序针对管排式换热器进行计算，并假设工质的流动是一维的，自然循环总流量平均分给每根换热器管道。每根换热器管道的重力压降、摩擦压降、漂移流压降的计算方法与其他部件管道类似，同时换热器的进口控制体、出口控制体要考虑形阻压降。

3) 数值方法

(22)

y(t0)=y0

(23)

由于核动力系统的复杂性，非线性常微分方程组通常是刚性方程组。本文在求解这类方程组时选用Gear算法[21-22]，它采用向后差分的隐式方法，是求解刚性问题的通用方法，该方法中配备了Adams方法和Gear刚性方法，可根据情况进行变阶或变步长积分，具有较好的稳定性。

2 程序开发

华龙一号PCS示意图如图1所示，主要结构包括上升管路、下降管路、换热器、换热水箱及水箱补水管路。其中PCS换热器通过自然循环将事故后安全壳内的热量导出至换热水箱，水箱补水管路在换热水箱的水位过低时投入运行，为换热水箱补水。

图1 PCS示意图

图2为程序计算流程，首先读取PCS各管道几何参数、安全壳隔间参数、初始热工水力状态等重要参数，之后对程序中各管道的温度、焓值等进行初始化。在程序计算中，首先计算流量、焓值的动态导数，之后调用Gear算法进行求解，并对参数进行更新。当运行到设定的停止时间时，程序结束。

图2 程序计算流程

图3为程序的结构层次图，主程序调用数据输入模块、初始化模块、动态计算模块和系统变量模块。系统变量模块用于COMMON模块的定义，方便各模块间的变量传递。动态计算模块分为微分方程数值求解模块及在每一步长结束后输出数据的数据输出模块。通过计算动态导数，形成雅各比矩阵，并调用Gear算法来求解微分方程组。动态导数计算模块主要计算各控制体的压降导数、焓值导数两部分。此外，动态导数计算模块还调用辅助模块进行物性、换热系数、摩擦阻力等的计算。

3 程序验证

将一维自然循环瞬态计算程序计算结果与Relap5程序计算结果进行对比，以验证程序计算结果的合理性。Relap5程序节点图如图4所示，整个回路的建模包括了PCS换热器、换热水箱、上升管路水平段、上升管路竖直段、下降管路水平段及下降管路竖直段。将回路中下降管路竖直段用部件125模拟，下降管路水平段用部件130模拟；上升管路竖直段用部件140模拟，上升管路水平段用部件145模拟，各部件中包含了所含弯头的阻力系数。Relap5中部件135为换热器，带有热构件用于提供功率，一维自然循环瞬态计算程序也使用恒定功率进行计算。

程序计算中，管道控制体数量划分会对计算结果产生一定影响。由于上升管路竖直段可能会出现两相流动的情形，因此需划分更多的节点。本文测试了4个算例，分别将上升管路竖直段划分为5、10、15、20控制体管道进行计算。4个算例的换热器流体温度如图5所示。

研究结果表明，在前期的单相阶段，控制体划分影响很小，各算例的计算结果基本重合。各算例的差别主要在后期。5控制体管道算例在后期直接发散，无法继续计算；10控制体管道算例的计算结果与15、20控制体管道算例有较为明显差别；15控制体管道与20控制体管道算例计算结果基本重合。

图3 程序结构层次图

图4 Relap5程序节点图

图5 不同上升管路竖直段控制体数目下的换热器流体温度

在本次验证中，换热器模块沿流动方向被划分为5个热工水力控制体；换热水箱模块沿流动方向被划分为5个水力学控制体；由于上升管路竖直段可能会出现两相流动的情形，因此沿流动方向被划分为20个热工水力控制体，上升管路水平段被划分为5个热工水力控制体；下降管路水平段和竖直段均沿流动方向被划分为5个热工水力控制体。利用Relap5程序和本文开发的程序计算了1个换热器功率为332 kW恒定值的工况，环境压力为0.1 MPa，回路初始水温为50 ℃，计算持续时间为50 000 s。

图6示出PCS换热器流体温度。由图6可见，本文程序与Relap5程序的温度变化趋势相同，且换热器温度在计算后期均基本达到稳定值，本文程序计算得到的换热器底部(控制体1)和顶部(控制体5)的流体温度分别为378 K和397 K，而Relap5程序计算得到的相应值为373 K和395 K。两个程序计算结果差异较小，产生差异的原因可能是由于Relap5中使用的某些辅助模型与本文程序不同，以及两流体模型与漂移流模型的差异所导致。由于Relap5计算中重要的辅助模型，如两相摩擦压降模型、换热模型等都无法自由修改，而本文程序针对PCS自然循环这一物理现象选取了更加合适的模型，因此也会导致计算结果产生了一定差异。

图6 换热器流体温度

图7 自然循环流量

图7示出两个程序计算得到的自然循环流量。在计算后期，随换热器流体温度的上升，系统内会出现两相流，导致流量出现震荡。由于Relap5程序使用两流体模型，而本文程序采用漂移流模型，因此流量变化有所差别。在计算后期，两个程序计算得到的流量均基本达到稳定值，本文程序计算得到稳定流量为2.34 kg/s，Relap5程序计算得到的稳定流量为2.43 kg/s。

验证结果表明，本文开发的一维PCS自然循环回路程序能用于分析核反应堆事故下PCS投入运行过程中重要参数的变化，评估PCS带走安全壳内热量的能力。

4 PCS热工水力特性分析

利用本文开发的程序对不同换热功率下PCS长期运行条件下的热工水力特性进行分析研究。图8示出不同的换热功率下PCS换热水箱出口流体的温度及换热器进出口控制体内的流体温度。由于未考虑换热水箱与环境空气的对流换热，换热水箱内水温持续上升，最终达到饱和状态，所以换热水箱出口的流体温度保持不变，为环境压力下的饱和温度。不同加热功率下，换热器进口温度变化不大，由于流体在换热器内被加热，换热器出口流体的温度大于进口温度，且随换热功率呈近似线性增长的趋势。换热功率为292.5 kW时换热器进出口温差约为7 K，而当换热功率为2 044 kW时换热器进出口温差增长到了约12 K。

图9示出不同换热功率下PCS传热管内流体的对流换热系数及PCS自然循环流量。由图10可见，换热系数和自然循环流量均随换热功率的增大而增加。在292.5～2 044 kW范围内，换热系数在1 000 kW·m-2·K-1左右，而流量变化较大，由292.5 kW下的约9 kg·s-1上升至2 044 kW下的约36 kg·s-1。

图8 不同位置的温度

图10示出不同换热功率下PCS上升管路竖直段的含气率。计算中将上升管路竖直段沿高度方向划分为20个控制体，1号控制体在下方，20号控制体在上方。当换热功率较小时，上升管路竖直段内流体保持为单相水。随换热器换热功率的增大，在竖直上升管道的出口附近会首先出现两相的情况，这是由于随高度的增加，管道内压力降低，管道内的流体温度达到了相应压力下的饱和温度，发生沸腾。当换热功率进一步增大，出现两相流的控制体数目不断增加(即出现两相流的管道长度越长)，且出现两相流的控制体内空泡份额逐渐增大。

图11示出瞬态计算时不同换热功率下换热水箱内的水位。水箱补水管路是本文程序针对华龙一号PCS开发的，它可在水箱水位较低时为水箱注水，保证PCS的正常运行。即在瞬态计算中，水箱温度达到饱和后，由于水的蒸发和蒸汽的排放导致水箱水位不断降低，当打开水箱补水开关后，水位降至某一设定值(2.5 m)时，水箱补水功能启动，补水管路开始以设定的流速向水箱中注入低温的水，直至水箱水位恢复至原来高度。由图11可见，换热功率越大，水箱水位下降越快。水箱补满时，水位在一段时间内维持不变的原因是低温水注入水箱后，水箱内的水并未达到饱和温度，未发生蒸发。

图9 换热系数及自然循环流量

图10 上升管路竖直段的含气率

图11 换热水箱水位

5 结论

本文针对华龙一号PCS的结构和特点，基于两相漂移流模型建立了一套适用于计算PCS热工水力现象的数学物理模型，并辅之以传热模型、阻力模型等，开发了专门适用于华龙一号PCS的瞬态热工水力分析程序。经过与Relap5程序结果的对比，二者误差较小，验证了本文程序的正确性和可靠性。采用该程序对PCS内的关键热工水力参数进行分析计算，得到了自然循环流量、换热功率、温度及含气率等的分布规律。

本文开发的PCS瞬态热工水力分析程序为事故后安全壳的安全分析提供了可靠工具，对PCS的设计和改善及对PCS冷却能力的评估均有重要的现实意义，且为后续开发能模拟带有PCS的安全壳内热工水力行为的程序打下基础。