钠冷快堆瞬态热工水力及安全分析程序开发

2019-10-30秋穗正张大林王式保王心安刘雅鹏

原子能科学技术 2019年10期

秋穗正，张大林，宋苹，王式保，梁禹，王心安，周磊，刘雅鹏

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

钠冷快堆是第4代核反应堆中技术最成熟、运行经验最丰富的堆型，具有闭式燃料循环和嬗变长寿命放射性废料的优势，因而在核能可持续发展战略中占有重要地位。钠冷快堆瞬态热工水力及安全特性的研究在钠冷快堆设计研发和安全评审中具有重要的作用，是国际上从事钠冷快堆研发国家的研究重点和热点。

然而，由于钠冷快堆采用液态金属为工质，其高导热性和低Pr的特征以及钠冷快堆特有的热工水力现象，在压水堆设计和安全分析中广泛适用的瞬态热工水力和安全分析程序不再适用。因此，世界各国针对钠冷快堆瞬态热工水力和安全分析，开发了自主知识产权的一维或三维程序，如美国的SAS4A[1-2]、法国的OASIS[3]和TRIO-U[4]、俄罗斯的RUBIN和GRIF[5]等，并开展了一系列堆内或堆外实验，对开发的程序进行验证，如美国的EBR-Ⅱ失流事故实验[6]、法国的凤凰堆自然循环实验[7]、日本的JOYO堆自然循环[8]和PLANDTL-DHX实验[9]等。由于国外快堆相关软件的核心源码不向中国开放，为适应中国示范快堆研发的需要，中国原子能科学研究院针对示范快堆的结构特点，自主开发了钠冷快堆系统分析程序FASYS[10]和FR-Sdaso[11]，并与华北电力大学合作开发了SAC-CFR程序[12]。哈尔滨工程大学针对中国实验快堆一回路设计，开发了稳态分析程序[13]。西安交通大学在RELAP5 MOD3.2的基础上进行了二次开发[14]，将其拓展到钠冷快堆的计算分析。这些研究对于中国快堆软件的自主可控做出了有益的探索。

西安交通大学在系统梳理国内外最新研究经验的基础上，进一步考虑钠冷快堆组件间耦合传热效应、流量再分配效应及局部三维效应等，开发了一套具有自主知识产权的钠冷快堆热工水力瞬态分析程序THACS(transient thermal-hydraulic analysis code for sodium-cooled fast reactors)，程序经过了丰富的V&V测试验证，具有更强的通用性和扩展性，模型涵盖了钠冷快堆一回路、二回路、余热排出系统中间回路及各回路内所有的重要部件，可用于不同类型钠冷快堆的瞬态热工水力和安全分析。本文主要介绍THACS程序的数学物理模型、程序结构及验证。

1 数学物理模型

为保证开发的钠冷快堆热工水力瞬态分析程序的通用性，以钠冷快堆部件为基础模块，对钠冷快堆一、二回路系统的关键部件进行详细建模。

1.1 堆芯功率模型

快中子反应堆的中子物理方程与热中子反应堆相同，目前大多数的快中子反应堆安全分析程序使用点堆模型[15]，假设堆芯的空间分布不变，时间和空间变量分离求解。为计算瞬态工况下堆芯功率的变化，建立了堆芯裂变功率模型[16]、衰变热模型[17]以及反应性反馈模型。

1) 裂变功率模型

(1)

(2)

式中：n为堆芯裂变功率，W；ρ为堆芯总反应性；β为缓发中子总份额；Λ为瞬发中子寿命，s；λi为第i组缓发中子衰变常量，s-1；Ci为第i组缓发中子裂变功率，W；j为缓发中子组数，一般取6。

2) 衰变功率模型

(3)

式中：hi为第i组衰变热功率，W；βhi为第i组衰变热功率份额；λhi为第i组衰变热衰变常量。

3) 反应性反馈模型

考虑钠冷快堆中主要的反馈包括燃料多普勒反馈、冷却剂密度反馈、燃料轴向膨胀反馈、堆芯径向膨胀反馈和控制棒、安全棒的密度反馈。

燃料的多普勒效应：

(4)

式中：ρD为多普勒效应引入的反应性；k为反应堆增殖因数；KD为多普勒常数；T为堆芯燃料的温度，K。

冷却剂密度反馈：

(5)

式中：ρvoid为冷却剂密度反馈反应性，冷却剂沸腾时也称失钠反应性；WNa为冷却剂密度反馈系数；ρNa为钠密度。

燃料轴向膨胀反馈：

(6)

式中：ρax为燃料的轴向膨胀反馈反应性；Rax为燃料的轴向膨胀反应性系数；Hfuel0为稳态时燃料的轴向高度，m；Hfuel为当前时刻燃料的轴向高度，m。

堆芯径向膨胀反馈：

(7)

式中：ρrad为堆芯的径向膨胀反馈反应性；Rrad为堆芯的径向膨胀反应性系数；R0为稳态时堆芯支撑结构半径，m；R为当前时刻堆芯支撑结构半径，m。

控制棒、安全棒密度反馈：

(8)

式中：ρabs为控制棒、安全棒密度反馈反应性；Wabs为控制棒、安全棒密度反应性系数；ρabsm为控制棒、安全棒密度，kg/m3。

1.2 堆芯热工水力模型

钠冷快堆采用带有六边形组件盒的封闭式组件设计，组件按照六边形排布布置在池式堆芯内，组件与组件之间通过其间的钠进行能量交换，而无质量和动量交换。鉴于钠冷快堆封闭组件特点，开发了一种并联多通道的堆芯热工水力模型，将堆芯按照组件类型设置为多个通道，每个通道代表1种类型的组件，为进一步准确分析堆芯内的关键参数，对燃料元件进行精细建模(图1)。

燃料棒(芯块和包壳)统一导热模型：

(9)

式中：ρ为燃料或包壳的密度，kg/m3；c为燃料或包壳的比热容，J/(kg·K)；T为燃料或包壳的温度，K；r为径向位置，m；λ为燃料或包壳的导热系数，W/(m·K)；QV为体积释热率，W/m3。

冷却剂的质量、动量和能量守恒方程：

(10)

(11)

(12)

式中：w为通道内质量流量，kg/s；A为流通面积，m2；De为当量直径，m；z为冷却剂轴向位置，m；ql为传入冷却剂的热流密度，W/m2；Sc为包壳的湿周，m。

图1 堆芯组件热工水力模型Fig.1 Thermal-hydraulic model of core assembly

1.3 钠池模型

钠池结构是钠冷快堆与压水堆的一个重要区别，钠冷快堆中通常设热钠池和冷钠池两部分。在正常运行工况下，堆芯出口和中间热交换器出口的冷却剂流量很大，冷、热钠池内的分布比较均匀，在自然循环工况下，钠池的冷却剂流量很小，钠池内混合不充分会出现明显的分层现象。为精确模拟钠池，开发了接管式多网格钠池模型。国际上研究发现钠池的分层多发生在冷却剂进出口位置附近[18-19]，因此钠池的分层位置如图2所示。

图2 钠池分层位置示意图Fig.2 Schematic diagram of layer division of sodium pool

借鉴RELAP5大空间流动换热模型，建立了钠池多网格模型，既可描述池内的热分层现象，亦可兼得程序较快的计算速度。多网格模型中，采用控制体的划分方式如图3所示，轴向不同控制体之间通过接管连接，计算流动和换热，同时考虑由于密度差产生的浮升力对钠池内流动的影响，主控制体考虑能量守恒和动量守恒，控制方程如下。

轴向动量方程：

(13)

径向质量守恒方程：

(14)

能量守恒方程：

(15)

图3 热钠池多网格控制体Fig.3 Control volume of multi-grip hot pool

1.4 泵模型

钠泵置于冷池中，汲取冷池内的液体钠，加压后大部分流入高压管，进而进入栅板联箱。程序采用离心泵进行计算，可直接输入扬程、流量、转速或是用四象限类比曲线求解。离心泵的主要参数有泵的扬程H、转矩Thy、体积流量Q和角速度w。典型的四象限类比曲线如图4所示。曲线以表格形式输入，因变量作为自变量的函数由表格查找或线性内插获得。已知主泵的体积流量和转速，利用泵的四象限特性类比曲线即可得到主泵的扬程。

图4 泵的四象限类比曲线Fig.4 Four-quadrant curve of pump

主泵的运行方式通常选用泵转速受控方式，一般有3种方式：1) 常转速，当w=0可模拟泵轴卡死；2) 转速按用户给出的时间表变化；3) 转速由控制系统控制。

1.5 热交换器模型

钠冷快堆系统中存在多种热交换器：一回路与中间回路之间、余热排出系统与钠池之间的钠-钠热交换器，中间回路与水回路之间的钠-水热交换器，余热排出系统与空气之间钠-空气热交换器。尽管各类型热交换器中换热介质及相态不同，但钠冷快堆系统中的热交换器多为管壳式结构，因此可通过建立热交换器一、二次侧不同介质的流动换热模型和管壁的导热模型来构建不同热交换器的数学模型。热交换器一、二次侧传热模型示意图如图5所示。

图5 热交换器一、二次侧传热模型示意图Fig.5 Schematic diagram of heat transfer model of heat exchangers’ primary and secondary sides

一、二次侧流体的流动换热模型遵循流体的基本质量、动量和能量守恒方程，如式(10)～(12)所示。由于传热管壁的厚度很薄，轴向导热可忽略，所以传热管壁的导热方程可简化为如下形式：

h1A1(T1-Tw1)-h2A2(Tw2-T2)

(16)

式中：ρw为管壁密度；Aw为传热面积；cp为管壁比定压热容；Tw为管壁温度；h1为一次侧传热系数；A1为一次侧传热面积；T1为一次侧流体温度；Tw1为靠近一次侧壁面温度；h2为二次侧传热系数；A2为二次侧传热面积；T2为二次侧流体温度；Tw2为靠近二次侧壁面温度。

1.6 辅助模型

辅助模型主要包括各种材料和冷却剂的物性参数模型、换热系数模型和流动阻力模型等。在钠冷快堆系统中包含的冷却剂工质有3种：液态钠、水以及空气，目前钠、水和空气的物性已非常成熟，THACS程序使用的钠物性来自Frin和Leibowitz的钠物性手册[20]，空气物性来自PROPATH[21]12.1版本，水和水蒸气物性关系式来源于IAPWS-IF97国际标准[22]。

钠冷快堆系统中液态金属钠在圆管和棒束通道内流动状态下的换热和阻力模型列于表1、2[23]。

表1 液钠在圆管和棒束内的换热模型Table 1 Heat models of liquid sodium in pipes and bundles

表2 液钠在圆管和棒束内的阻力模型Table 2 Friction models of liquid sodium in pipes and bundles

2 数值算法与程序设计

2.1 数值算法

通过钠冷快堆系统的数学物理模型的建立，构建了一套封闭的求解钠冷快堆系统热工水力特性的方程组，方程组的基本形式为各热工水力参数对时间的导数，通过方程的离散可转化为常微分方程组初值问题的求解：

(17)

本研究所建立的模型包括点堆中子动力学模型和系统热工水力模型，这些模型中瞬发中子的动态过程最快，在10-7s量级，缓发中子次之，为10-4s量级，通道入出口压力、流量、温度的变化相对很慢，不同瞬态过程时间比达108甚至更高量级，该类方程组为病态或刚性方程组，Gear方法[24]采用向后隐式差分方法，设计了一种病态稳定策略，可做到步长与特征值乘积大时是精确的，从而很好地跟踪解的快变部分；而对两者乘积小时又是稳定的，即当特征值很小时也不会失真。Gear方法采用牛顿迭代法进行隐式求解，并相应地利用矩阵系数结构的特点用直接法解线性方程，每前进1个步长解隐式方程组所需的工作量较小，这就加快了计算速度。此外，Gear方法能自启动，易实现变阶和变步长。在Gear方法中还配备了阿当姆斯(Adams)方法，可根据方程组的刚性强弱自动选取合适的算法，既保证了求解精度又提高了计算速度。

2.2 程序架构

图6 THACS程序结构Fig.6 Structure of THACS code

程序采用面向对象的模块化建模方法，为增强程序的可读性，提高程序的可移植性并赋予程序较好的通用性，对不同部件进行独立性较强的模块化设计。各部件的模块均有独立的输入模块、初始化模块和时间导数计算模块，各模块既可独立运行，又可在耦合模块中由主程序调用共同求解，从而方便程序的维护和拓展。THACS程序总体模块化设计如图6所示，其中耦合模块、物性模块、辅助模块以及数值计算模块属于各部件公用的模块，而每个部件模块内部包含每个部件的输入、初始化、导数计算以及边界传入传出部分，删除或添加1个模块不会影响其他模块。

3 程序初步验证

为更全面了解钠冷快堆在瞬态事故过程中的热工水力和安全特性，提高THACS程序的计算可靠性和准确性，使其进一步接近大型商业程序的计算水平，采用国际上公开发表或内部的实验数据，对程序进行了大量的验证和确认(V&V)。本文主要介绍采用美国实验快堆EBR-Ⅱ的有保护失流事故SHRT-17[6]，对THACS程序瞬态计算能力进行的验证。

EBR-Ⅱ电站是美国阿贡国家实验室建造的采用金属燃料的钠冷快堆电站，该电站1964年建成并于1994年退役。在其寿命后期，在EBR-Ⅱ电站上开展了一系列用于验证钠冷快堆固有安全性的试验。本文选用EBR-Ⅱ基准题中的SHRT-17实验对程序进行验证，该实验是有保护失流事故实验，实验目的是验证反应堆依靠自然循环带出热量的能力。EBR-Ⅱ基准题由IAEA启动，全世界10个国家18个研究机构参加，持续4年(2012.5—2016.4)，西安交通大学参与了此项基准题。

EBR-Ⅱ是单池式的钠冷快堆，其额定热功率为62.5 MW，电功率为20 MW，包含3个回路系统，其主回路系统结构如图7所示。图8示出了实验工况下安装了装有热电偶的XX09和XX10测量组件，其中，XX09 为有燃料组件，XX10 为无燃料的反射层组件。

图7 EBR-Ⅱ电站的主回路系统Fig.7 Primary loop system of EBR-Ⅱ

图8 EBR-Ⅱ堆芯及测量组件的布置示意图Fig.8 Configuration of EBR-Ⅱ core and testing assembly

根据提供的反应堆结构和边界条件，使用程序对EBR-Ⅱ反应堆进行模拟，其中主要计算的部件如图9所示。事故设定为：初始状态反应堆在额定功率下运行，事故发生后反应堆两台主泵失去动力惰转，与此同时，反应堆停堆系统启动，反应堆紧急停堆。

图10为高压腔室进口温度，从图中可看出，实验数据抖动剧烈，但稳定值在625 ℃附近，本程序(XJTU)的计算结果和实验值符合较好。

图11示出了Z型管进口温度和IHX二次侧出口温度的计算值与实验值随时间的变化，事故发生后Z型管进口温度迅速上升，达到峰值温度后开始下降，最后达到一个相对稳定值。由于堆芯功率下降，通过IHX传递给中间回路的流量减小，因此IHX二次侧出口温度持续下降。从图11可看出，实验值和计算值随时间的变化趋势一致，其峰值相对误差在10%以内，吻合较好。根据图中堆芯出口温度和Z型管进口温度对比分析发现，Z型管进口温度计算误差是由堆芯出口腔室的温度分层造成的。堆芯上腔室结构如图12所示，在自然循环条件下堆芯上腔室会出现轴向温度分层现象。图11显示Z型管进口温度实验值的最高点出现时间与堆芯出口平均温度最高点出现时间相同，说明从堆芯流出的热流体未与出口腔室内的冷流体充分混合即流向Z型管入口，而计算模型采用集总参数法，堆芯出口的热流体首先与出口腔室内的冷流体充分混合，因此计算结果最高点的出现时间晚且最高峰值温度偏低。