基于RELAP5和子通道程序的熔盐冷却快堆多尺度热工流体耦合程序开发及应用

宋诗阳 程懋松 林 铭 戴志敏

1（中国科学院上海应用物理研究所 上海 201800）

2（中国科学院大学 北京 100049）

目前，广泛使用的系统安全分析程序大多采用一维模型，不能有效地模拟反应堆堆芯的流场现象，即使最新的RELAP5-3D，三维仿真功能也具有较大的局限性［1］。相比之下，子通道程序能够有效建模堆芯内部各种通道、棒束、绕丝等对象，能够计算堆芯或燃料组件内部冷却剂质量流、横向流、压力以及温度分布等关键热工流体参数。近年来，多尺度热工流体分析程序的开发成为反应堆系统设计和安全分析的重要内容，随着计算机技术的高速发展，国际上已有众多学者进行系统分析程序和三维热工流体程序的耦合研究。如韩国的Jeong 等［2］耦合COBRA/RELAP5 （Computational Brewing Application/ Reactor Excursion and Leak Analysis Program）程序，使用LOFT L2-3 破口事故（Loss of Coolant Accident，LOCA）基准题数据对其验证，计算结果与单独使用RELAP 程序计算结果和实验值三者相符；德国卡尔斯鲁厄理工学院开发了KITATHLET（Karlsruher Institute für Technologie Analyse der Thermohydraulik von lecks und Transienten）与OpenFOAM（Open Field Operation And Manipulation）多尺度耦合程序；法国原子能委员会CEA（The French Atomic Energy Commission）和中国核动力设计研究院开发了CATHARE（Code for Analysis of THermalhydraulics during an Accident of Reactor and safety Evaluation）与TRIO_U 耦合程序，两者共同与钠冷快堆PHENIX 的自然循环实验进行对比验证［3-5］；中国原子能科学研究院的王军等［6］开发了RELAP5/MOD3 和THAS-PC4（Thermal-Hydraulic Analysis code for Personal Computer）的耦合程序，并且通过秦山核电厂失流事故进行验证；清华大学刘余等［7］采用并行虚拟机技术（Process Virtual Machine，PVM），耦合了RELAP5 与CFX 和COBRA程序，并验证了耦合程序的正确性。但目前开发的多尺度热工流体耦合程序主要是针对传统压水堆、钠冷堆、铅冷堆等堆型［8-9］，缺乏适用于熔盐冷却堆的多尺度热工流体耦合程序。

RELAP5系统分析程序适用于瞬态或事故工况的快速预测与分析，中国科学院上海应用物理研究所开发了适用于熔盐堆计算分析的RELAP5-TMSR版本，可应用于熔盐堆等核能系统的瞬态和事故分析。但全自然循环熔盐堆堆芯建模的精细度有着更高的要求，而RELAP5-TMSR具有一维局限性，不能对堆芯进行三维精细建模，无法准确地反映熔盐冷却堆堆芯的三维热工流体现象［10-11］。通过耦合，已经开发且经过验证的适用于氯盐快堆的子通道分析程序ThorSUBTH可解决该问题。ThorSUBTH子通道分析程序选用了合适的高温熔盐工质换热模型、压降模型和湍流交混模型，可以对堆芯内部三维热工流体现象进行更为精确的分析［12-13］。通过耦合的方式，综合利用两者的优点，有效提高反应堆系统分析程序的准确性，相对于单独的RELAP5-TMSR 系统分析程序，能够得到堆芯的三维流场分布、温度分布和热点位置等，更准确地获得具体三维物理场信息和揭示熔盐堆堆芯的三维物理现象。

1 多尺度耦合原理及方法

1.1 耦合方法

常见的耦合方法主要有显式、隐式和半隐式，目前国际上基本采用显式或半隐式的方式，具体实现方式有：1）利用DLL（Dynamic Link Library）技术和三维热工流体程序进行数据交换，此方法涉及C 语言和Fortran 的混合编程，实现方式相对较为复杂［14］，该方式采用显式或半隐式耦合；2）使用PVM并行虚拟机编写接口控制程序来实现数据交换，适用于显式或半隐式耦合方式［7］；3）直接对两程序进行二次开发，添加数据交互模块，此方式实现成本较低，适用于显式耦合［15］。

本文对系统分析程序RELAP5-TMSR和子通道程序ThorSUBTH 进行外部显式耦合。耦合程序的流程如图1所示。在RELAP5-TMSR计算时间步长完成后将耦合所需传递的数据保存，并在ThorSUBTH每个时间步长计算开始前，将耦合参数读入，作为该时间步长的边界条件参数。等待该时间步长计算完成，将需要传递的参数保存传递至RELAP5-TMSR，作为下一个时间步长的边界条件，以此方式重复迭代。

图1 耦合程序流程图Fig.1 Flowchart of coupling program

为实现两程序数据交互和控制，对计算域划分，目前有域分解和域重叠两种方式［16］。图2为系统程序和子通道程序模拟区域划分的两个方式，其中（a）为仅使用系统程序不进行划分，（b）为域分解方式，（c）为域重叠方式，本文采用（b）域分解的方式。

图2 计算域划分方式Fig.2 Division of computational domain

RELAP5-TMSR 和ThorSUBTH 数 据 传 递 参 数如图3 所示，其中：Ti为RELAP5-TMSR 系统出口流体温度，Wi为RELAP5-TMSR系统出口质量流量，Pi为RELAP5-TMSR 系统入口流体压力，Po为ThorSUBTH 入口流体压力，To为ThorSUBTH 出口流体温度，Wo为ThorSUBTH出口质量流量。

图3 耦合程序数据交换示意图Fig.3 Parameters exchange of coupled program

由于ThorSUBTH是三维程序，耦合界面处数据作为一维RELAP5-TMSR 程序的入口边界条件，需要进行处理，其中质量流量为各通道的总和如式（1），温度和压力差采用式（2）和（3）加权的方式。

1.2 耦合时间步长控制

耦合程序使用显示耦合迭代的方式，相较于隐式迭代在每个计算时间步内满足一定的收敛条件后再进行下一时间步计算的方式，时间上具有一阶精度，计算速度相对更快，会牺牲一定的稳定性。但通过合理控制时间步长，显式迭代法也可以得到较好的计算结果［17］。

显式耦合的数据传递只在每一个时间步结束后发生，两程序顺序进行计算，但由于RELAP5-TMSR和ThorSUBTH采用了不同的时间步长控制策略，需要对两程序的时间耦合进行协调控制。RELAP5-TMSR 程序采用自适应时间步长，计算时间步长根据具体Courant数进行确定，程序在本时间步检测到计算未达到收敛要求时，会将时间步长减半，重复进行该时间步的计算，直到满足要求后进行下一个时间步的计算。而ThorSUBTH 程序采用固定时间步长，每个时间步结束后都会进行下一个时间步的计算。

子通道计算时间步长通常大于系统分析程序，两个程序不同的时间步长控制策略会导致耦合时间发生不同步的现象，因此本文中将子通道时间步长作为耦合程序数据传递的时间步长，即子通道程序固定时间步长，将RELAP5-TMSR 的数据交换方式从每个时间步都交换数据变更为时间步积累到指定固定时间长度后再进行交换，该方法可保证两程序的计算时间同步，如图4所示。

图4 耦合程序时间步长迭代方法Fig.4 Time step iteration method for coupled program

2 耦合程序验证

为了确保耦合方法的正确性，需要对程序进行验证。通常验证有两种方式：实验数据验证和程序之间计算结果的对比验证（CODE-TO-CODE验证），本文采用CODE-TO-CODE 验证方式。选择中国核动力研究设计院刘余等［18］采用的开式水平管道流动问题和一个自定义的闭式循环回路流动问题开展耦合程序验证。由于RELAP5-TMSR 和ThorSUBTH在程序的耦合开发过程中没有涉及原程序的计算方法和模型的修改，只针对数据的输入和输出方式进行修改，因此只需要验证两个程序数据传输的正确性即可。

2.1 开式水平圆管瞬态验证

2.1.1验证模型

将本文耦合程序RELAP5-TMSR/ThorSUBTH计算结果与文献［18］中计算结果进行对比验证。模型为水平方向上长度为16 m圆管，管内流通面积为0.4 m2，管内流动工质为300 K 液态水，壁面为绝热条件，出口压力为0.15 MPa。管入口侧流速V和时间t关系式如下：

单独模型和分段模型如图5 所示，每一部分划分20节块，分别采用不同程序进行建模和运行。

图5 单独和耦合计算模型Fig.5 Single and coupled calculation model

2.1.2结果分析

首先分别使用RELAP5-TMSR 和ThorSUBTH对该水平圆管模型进行单独建模计算，管道划分为40 个节块。得到管前部节块1 和中间节块20 处的压力随时间变化曲线，与单独使用RELAP5 和CFX计算的节块1 和节块20 处的压强进行对比，结果如图6所示，图6括号内为分段编号，可以看出，数据吻合较好，以此确保对比验证模型的一致性。

图6 单独程序计算结果Fig.6 Comparison of calculation results of separate single code

然后将圆管进行分段模拟，其中使用ThorSUBTH 模拟Part1，RELAP5-TMSR 模拟Part2。节块1和节块10、节块21和节块30处压力随时间变化曲线如图7（a）、（b）所示，并且将结果和单独使用RELAP5和CFX/RELAP5耦合程序的计算结果进行对比。流量在1 s前和11 s后都不发生变化，因此压力为恒定值，但在1~11 s中，由于管入口处流速V随时间t正弦变化，在1 s时刻流速变化率最大，对应管道头部和尾部压力差最大，在约2.2 s 时刻，流速变化率最小，对应管道头部和尾部压力差最小。结果显示，耦合程序和文献［18］单独使用RELAP5 程序以及CFX与RELAP5耦合程序的结果基本一致。

图7 耦合程序计算结果Fig.7 Calculation result of coupling code

相反地，使用RELAP5-TMSR 模拟Part1，ThorSUBTH 模拟Part2，在节块1 和节块10、节块21和节块30 处的压力随时间变化曲线如图7（c）、（d）所示，结果显示，耦合程序与文献［18］单独使用RELAP5程序以及CFX/RELAP5耦合程序的结果也基本一致。通过该算例初步验证了本耦合程序数据传输的正确性。

2.2 简单闭式回路瞬态验证

2.2.1验证模型

本节采用简单闭式氯盐回路瞬态工况验证多尺度耦合程序。图8 为简单闭式氯盐回路示意图，其中Pipe1~Pipe7 长度分别为5 m、2 m、6 m、2 m、5 m、5 m和5 m，管道截面面积均为0.4 m2。

图8 简单闭式氯盐回路示意图Fig.8 Schematic diagram of simple chloride salt closed loop

单独RELAP5-TMSR和耦合程序建模节点图如图9所示，耦合程序中将Pipe3区域更换为子通道模型。在Pipe5末端设置压力控制器TDV，在Pipe6和Pipe7中间设有一个熔盐泵作为冷却剂驱动装置，冷却剂顺时针流动。闭合回路的初始温度设置为500 ℃，稳压器压力为一个大气压101 kPa。

图9 简单闭式氯盐回路节点示意图Fig.9 Nodal diagram of simple chloride salt closed loop

2.2.2结果分析

回路中流速由泵进行控制，流速随时间变化，在0 s 时刻回路初始流速为2 m·s-1，10 s 时刻开始至40 s 时刻线性增加流速至4 m·s-1，稳定后在60 s 时开始到80 s时刻流速线性降低至1 m·s-1。RELAP5-TMSR 模型中Pipe3 入口处以及耦合模型中子通道入口处的流速变化曲线如图10所示，两程序计算结果基本一致，说明耦合接口处的数据传递正确。

图10 闭合回路流速变化曲线Fig.10 Change curve of closed loop velocity

图11表示Pipe2 出口处压力随时间变化曲线，在10~40 s 时压力曲线随回路的流速增加而升高，40 s 后趋于稳定值；出口压力随着回路在60 s 时流速快速降低而迅速下降，在80 s 时开始稳定。结果表明：RELAP5-TMSR 程序与RELAP5-TMSR/ThorSUBTH耦合程序计算结果吻合较好，说明耦合程序能够较好地计算闭合回路的压力变化。

图11 Pipe2出口压力变化曲线Fig.11 Change curve of Pipe2 outlet pressure

上述开式和闭式回路的瞬态工况验证结果显示 ，RELAP5-TMSR 程 序 与 RELAP5-TMSR/ThorSUBTH耦合程序计算结果吻合较好，初步验证耦合程序具有足够的准确性以及可信度，可以用于开展氯盐冷却快堆的正常工况以及瞬态事故工况的多尺度热工流体特性研究。

3 多尺度热工流体耦合程序应用

为验证该耦合程序的适用性，以小型自然循环氯盐冷却快堆（Small Natural Circulation Chloride Cooled Fast Reactor，SN3CFR）为对象，使用耦合程序开展正常工况和瞬态事故工况的计算与分析。一方面对该反应堆热工设计和安全性进行计算与分析，另一方面进一步说明了耦合程序的有效性和适用性。

3.1 45 MW自然循环氯盐冷却快堆

采用自然循环非能动技术的SN3CFR 具有高温、常压、较硬的中子能谱等特点，在灵活性、经济性和固有安全性上具有极大的优势和潜力。堆芯设计借鉴MCCFR［19-22］，SN3CFR 的热功率为45 MW，燃料包壳采用SiC，结构材料采用Hastelloy-N 合金。堆芯为无组件固态燃料棒，分为5个燃料区域，其中1、2区具有较低的富集度，3、4、5区具有较高的富集度。反应堆的活性区直径和高度分别为130 cm 和125 cm，设计详细参数如表1 所示。堆芯部分使用子通道程序进行建模，在轴向分为共18 层，1/12 堆芯如图12所示。

表1 堆芯设计参数Table 1 Parameters of the core

图12 1/12堆芯示意图Fig.12 Configuration diagram of 1/12 core

图13给出了SN3CFR 系统程序的模型节点图，其中100TDV和110TDJ组成系统程序的入口边界，320TDV 和310TDJ 组成系统程序的出口边界，160和260分别为上升和下降通道，180和280分别为一回路的上腔室和下腔室，200为稳压器，220和430分别为换热器的一次侧和二次侧，410TDV 和420TDJ为二回路入口边界，450TDV 为二回路出口边界。其中一回路采用全自然循环冷却，二回路采用强迫循环。

图13 SN3CFR RELAP5-TMSR模型节块Fig.13 RELAP5-TMSR nodalization of SN3CFR

3.2 额定功率工况分析

采用45 MW SN3CFR 额定满功率工况对耦合程序进行适用性验证，并且与单独使用RELAP5-TMSR程序计算结果对比。为了保证耦合程序计算精度和提高计算效率，也已开展了耦合程序数据交互时间步长的敏感性验证分析，选取0.01 s、0.02 s、0.05 s、0.1 s、0.2 s、0.5 s 一组时间步长，对圆管算例进行了模拟，结果表明：除0.5 s数据有很小差距外，其他时间步长结果均保持了很好的一致性。综合考虑计算精度和计算效率，选取0.1 s作为后续瞬态计算（包括弹棒事故）的耦合时间步长。

表2给定了反应堆设计限制，满功率设计参数与耦合程序计算值如表3 所列，正常稳态运行时1/12 堆芯三维温度分布如图14 所示。在完全自然循环的情况下，堆芯进出口温差达到103.88 ℃。燃料棒的归一化功率分布如图15 所示，1 和2 区富集度较低，3、4 和5 区富集度较高，相对功率最高棒束功率因子为1.327 8，位于第4 区。燃料棒最热层为第10 层，其温度分布如图16 所示，燃料棒包壳最高温度位置位于270号燃料棒，热点温度为780.82 ℃，远小于900 ℃的设计限值，冷却剂出口温度子通道对应的编号为289，温度为705.33 ℃。

图14 1/12堆芯温度分布Fig.14 Temperature distribution of fuel rods in 1/12 reactor core

图15 燃料棒的径向功率分布Fig.15 Radial power distribution of fuel rods

图16 最热层温度分布Fig.16 Temperature distribution of the hottest layer of fuel rods

表2 反应堆设计限值[23]Table 2 Design limits of reactor[23]

表3 稳态运行参数Table 3 Operation parameters of steady state

3.3 瞬态工况分析

根据SN3CFR采用熔盐作为冷却剂和使用自然循环非能动技术的特点，可确定其设计基准事故主要有两类：1）反应性引入事故；2）丧失热阱事故。本文在满功率稳态运行的基础上，主要针对反应性引入事故进行瞬态分析。事故具体序列为：从正常运行工况开始，在第100 s时，2 s内线性引入一个最大控制棒价值（1.38×10-3）反应性，并且未进行紧急停堆。

反应性和反应堆相对功率随时间变化曲线如图17 所示，可以看出，由于反应性引入，在第100 s 时，反应性急剧上升，导致功率在前期瞬态最高值达到68.23 MW，由于负反馈效应，功率将缓慢回落，之后建立新的稳定自然循环，稳定后功率为56.55 MW。由于功率的上升，导致堆芯冷却剂出口温度随之上升，从603.88 ℃上升至630.37 ℃，入口温度从500.00 ℃上升至509.94 ℃，变化曲线如图18 所示，冷却剂需要经过一回路循环再次回到堆芯，因此堆芯入口冷却剂温度的上升相对会有一定的时间延迟。

图17 反应性及反应堆相对功率变化曲线Fig.17 Change curve of reactivity and reactor relative power

瞬态过程中热点出现在第270 号燃料棒，对应最高燃料包壳和燃料中心线温度随时间变化如图19 所示，在反应性引入时急剧上升，由于此时堆芯内冷却剂流量上升速度与温度上升速度相比较慢，因此包壳和中心线的热量持续累积，温度仍会有一定程度上升。燃料包壳最高温度从780.98 ℃上升至瞬态最高温度831.47 ℃，小于瞬态包壳温度设计限值1 100 ℃。燃料中心线最高温度从875.58 ℃上升至瞬态最高温度953.69℃，远小于燃料设计最高温度3 035 ℃，瞬态最高温度时刻1/12 堆芯温度分布如图20所示。但随着反应堆堆芯温度升高，自然循环能力加强，堆芯流量从441.66 kg·s-1上升至477.98 kg·s-1，如图18所示。在无人干预情况下，进入新的自然循环稳态，在1 000 s左右时温度将会逐渐下降至稳定值，燃料包壳最高温度稳定在828.63 ℃，燃料中心线最高温度稳定在947.60 ℃。在瞬态过程中，耦合程序与单独RELAP5-TMSR 程序计算结果相比较，堆芯出口温度约低1 ℃，分析是由于子通道程序计算堆芯横流交混后，流速和温度分布的差异所导致，该差异在合理范围内，吻合较好，验证了该耦合程序的适用性，另外冷却剂、燃料包壳、燃料中心线温度全程都在设计限度内，并且具有较大的裕度，因此本堆具有良好的稳定性和安全性。

图18 耦合程序与单独RELAP5-TMSR进出口温度及质量流量变化曲线Fig.18 Change curves of inlet outlet temperature and mass flow rate in couple program and RELAP5-TMSR

图19 最热燃料棒包壳和中心线温度变化曲线Fig.19 Temperature change curves of hottest rod clad and central

图20 1/12堆芯温度分布Fig.20 Temperature distribution of 1/12 reactor core at transient maximum moment

4 结语

熔盐冷却固态燃料快堆以其独特的安全性和经济性越来越受到人们关注。为了更为准确有效地进行瞬态分析和安全评估，国际上广泛开展多尺度热工流体耦合程序的研究。本文为提高RELAP5-TMSR程序在熔盐冷却快堆瞬态分析和安全评估中的适用性和准确性，以系统分析程序RELAP5-TMSR和子通道程序ThorSUBTH为基础，使用外部显式方法对两者进行耦合，开发了RELAP5-TMSR/ThorSUBTH多尺度热工流体耦合程序，并采用水平圆管模型和简单闭合回路模型验证了耦合程序的正确性。将多尺度热工流体耦合程序应用于小型自然循环氯盐冷却快堆SN3CFR 稳态及瞬态分析，并与单独RELAP5-TMSR 计算结果进行对比，额定功率工况下计算值与设计值符合良好，在事故工况下满足设计限值，进一步体现了多尺度热工流体耦合程序良好的适用性。相较于单独一维系统分析程序，多尺度热工流体耦合程序可以充分发挥不同尺度程序特点，即保证了堆芯重要热工流体参数的计算精度，又确保系统计算分析的速度。多尺度热工流体耦合程序的开发及验证能够为熔盐冷却固态燃料熔盐堆的系统设计、优化和安全分析提供有力的支撑工具，具有重要的工程意义。

作者贡献声明宋诗阳：设计研究方案、负责研究方案具体实施、进行模拟计算、数据分析及论文的撰写；程懋松：提出研究思路、稿件的审阅与修订以及研究进度的监督；林铭：负责提供技术支持与指导以及模型的协助；戴志敏：负责研究项目管理、研究资金获取。