何 帆，蔡翔舟，郭 威，何 龙，崔 蕾，赵 恒

(1.中国科学院 上海应用物理研究所，上海 201800；2.中国科学院大学，北京 100049； 3.中国科学院 先进核能创新研究院，上海 201800)

RELAP5是美国爱德华国家工程实验室为美国核管会开发的轻水堆冷却系统事故工况的瞬态行为最佳估算程序[1-3]，已被广泛应用于反应堆的瞬态事故分析和安全评审等方面，并被核工业接受成为安全分析工具。在RELAP5/MOD4.0版本中，添加了铅铋合金、熔盐等多种计算流体，其可应用于熔盐堆等第4代核能系统的瞬态事故分析[4-5]。同时，随着流体动力学方法(如商业化程序FLUENT)和计算机性能的高速发展，有关反应堆的三维数值模拟越来越多[6-8]，堆芯上下腔室、冷却剂通道等局部构件计算流体动力学(CFD)模拟已被报道。以FLUENT为代表的CFD程序将会成为开发第4代核反应堆的一个强大工具。为综合利用系统代码和CFD程序的优点，国际上有学者已开始系统程序与CFD程序间的耦合研究，开发了一些系统/CFD耦合程序，并进行了初步验证，如Aumiller等[9]首先利用并行虚拟机技术(PVM)实现了RELAP5-3D和CFD代码的耦合，从原理上证明了RELAP5-3D/CFD显式耦合是可行的；Papukchiev等[10]将ATHLET程序和CFX软件耦合，模拟了流体在管道流动、传热等物理过程；Anderson等[11]针对超高温气冷堆冷却剂流出堆芯进入出口腔室的热混合问题，利用PVM耦合了RELAP5和CFD代码，在出口腔室建立了三维流动效应模型。

本文讨论显式耦合RELAP5/MOD4.0和FLUENT方法，通过RELAP5/MOD4.0源代码的二次开发和FLUENT用户自定义函数(UDF)功能，RELAP5和FLUENT会分别在每个时间步开始时读入耦合计算所需要的变量进行计算，并在该时间步结束时输出耦合计算所需要的变量，然后再进行下一时间步的计算。本文为测试耦合程序的正确性，分别使用RELAP5/FLUENT耦合程序和独立的RELAP5或FLUENT对水平圆管进行模拟分析。最后，利用耦合程序对2 MW熔盐堆进行稳态工况模拟和功率突变的瞬态分析。

1 耦合方法

程序间的数据信息交换是耦合程序的本质[12]。在已有的一维系统程序和三维CFD程序耦合的研究成果中，主要是基于PVM来实现。在PVM中，额外编写接口控制程序，利用PVM内部所提供的Send和Recv函数在每个时间步的计算中调用系统程序和CFD程序的边界参数，并将这些参数值分别赋予相应的程序作为边界条件以进行下一时间步的计算。但对RELAP5/FLUENT显式耦合而言，这种方法相对较为复杂，且需较高的编程能力[13-14]。在本文中直接对RELAP5/MOD4.0源代码进行二次开发输入输出接口模型，同时FLUENT基于UDF功能直接读入和输出边界参数，这种方法相对简单，也较易实现。

图1为RELAP5和FLUENT在耦合计算时需实现交换的相关参数。FLUENT模拟区域的耦合出口边界的输出参数主要包括流体压力p、空泡份额α、流体温度T和质量流量W，将作为下游RELAP5计算区域的入口边界参数；同时RELAP5模拟区域也需反馈流体压力p、流体回流温度T及空泡份额α来作为FLUENT压力出口边界条件。对FLUENT的速度入口边界而言，需输入上游RELAP5流体压力p、空泡份额α、流体温度T和流体流速v等参数，同时反馈给上游RELAP5包括流体压力p、空泡份额α、流体温度T等参数来作为RELAP5计算的边界条件。

图1 耦合程序交换参数示意图Fig.1 Exchanged parameter of coupled program

(1)

(2)

其中：Ai为FLUENT边界单元格的面积；n为时间步数。输出的平均温度采用质量加权的方式：

(3)

RELAP5整体结构如图2所示，其主要分为3部分：输入部分(INPUT)，瞬态/稳态计算部分(TRNCTL)，输出部分(OUTPUT)。瞬态/稳态计算部分负责将输入数据按照变量所属类型进行相应的热工水力、控制逻辑等方面的计算，并将计算结果保存在输出文件中，它由TRNSET、TRAN和TRNFIN等子程序组成。其中，TRAN控制瞬态求解的步进，进行矩阵运算，是程序最重要的部分，也是程序计算过程中耗费时间最长的部分。在RELAP5瞬态计算过程中，TRAN子程序下的Dtstep子程序模块控制时间步长的大小与计算结果输出、绘图编辑的频率。在瞬态/稳态计算部分中，Dtstep子程序模块中的Majout子程序控制大编辑输出(图3)，其特点如下：1) 该模块全程参与瞬态/稳态部分的计算；2) 按照输入卡设定将计算过程中相关变量的数据输出到输出文件中；3) 程序中含有大量的变量信息，并可按照一定格式将其输出到输出文件中；4) 该子程序模块具备一定的扩展性，可添加变量，并进行赋值、逻辑判断等操作。

图2 RELAP5主程序3大部分Fig.2 Three parts of RELAP5 main program

图3 瞬态/稳态计算部分的模块结构Fig.3 Modular structure of transient/steady state calculation part

基于RELAP5的上述特点，本文将程序所需要的输入输出模块放置在该子程序模块中，即从外部(如FLUENT)读取数据，并赋值给程序变量；或将程序变量的值输出到外部。在热工参数输出功能上，RELAP5提供了小编辑输出功能，通过设置输出变量，可将数据输出到输出文件中，本文耦合程序动态输出功能则借鉴了该特点。通过修改小编辑卡相关模块，禁止小编辑卡输出，确保小编辑卡不会对耦合程序的输出产生干扰。同时，将小编辑卡所含变量的数据传递给输出变量；从而完成小编辑卡的输出变量向动态输出变量的转换。

如原RELAP5输入卡301 mflowj 110010000命令，将会输出编号为110的管道的质量流到RELAP5输出文件中，程序修改后的该命令将直接把质量流输出到外部。为了与原输入卡小编辑输出区分，同时扩展动态输出数据的数量，本文进一步修改RELAP5小编辑模块，在rmiedt子程序模块中将输入卡编号范围从301～399修改为5 001～5 999；在Majout子程序中添加小编辑输出代码。程序修改后可通过输入卡命令5001 mflowj 110010000直接将结果输出到外部。

在RELAP5中，边界条件主要是由时间相关变量来设置，包括时间相关控制体(TDV)、时间相关接管(TDJ)和随时间变化的功率等。这些边界条件参数主要包括液体/气体的流速或流量，液体/气体的温度、压力以及流体的空泡率等，其在输入卡中以表格的形式进行输入，表值在RELAP5读入输入文件后保存在对应的变量中，在计算过程中通过插值函数来获得不同时刻的变量值。因此，可通过修改表值来实现修改边界条件的目的。以RELAP5输入卡110TDJ质量流为例：1100201 0.0 0.0，1100201 1.0 2.0，表示在0.0～1.0 s期间，110TDJ线性引入了2.0 kg/s的质量流量。若在计算中动态地将该表值0.0和2.0修改为3.0，则表示当前1.0 s时，110TDJ的质量流量为3.0 kg/s。其他边界条件参数如流体的温度、压力等参数也可以通过这种方式输入。通过这种动态修改表值的方式可实现将FLUENT边界条件传递给RELAP5。

耦合程序计算流程示意图如图4所示。在耦合程序中，RELAP5和FLUENT分别读入各自的输入文件，并进行初始化。鉴于RELAP5是系统级程序而FLUENT多应用于局部构件的分析，RELAP5首先进行第1个时间步的计算更为适宜。RELAP5在第1个时间步计算结束后，将FLUENT所需要的边界参数传递给FLUENT；FLUENT在进行该时间步的计算后将边界参数传递给RELAP5，以便RELAP5进行下一时间步的计算。通过这种方式在每个时间步不停地更新边界参数，直到最后1个时间步完成整个瞬态计算过程。

2 验证及应用

2.1 flibe熔盐在水平圆形管道流动问题

1) 圆管流动问题描述

首先利用flibe熔盐在水平圆管流动问题来测试耦合程序。假设一水平圆管长1.0 m，流通面积为3.14 159×10-4m2，初始压力为0.101 325 MPa，初始流速为1.0 kg/s，初始流体温度为873.15 K，管壁粗糙度为10-6m且绝热。整个瞬态模拟时间为10 s，水平圆管入口处flibe熔盐的流速初始时刻为1.0 kg/s，在2～4 s期间线性增加至1.5 kg/s，6～8 s期间线性减少至1.0 kg/s并保持不变直到10 s结束。

图4 耦合程序计算流程示意图Fig.4 Schematic of coupled program calculation process

水平圆管入口处flibe熔盐的温度初始时刻为873.15 K，4～6 s期间温度线性增加到893.15 K并保持不变直到10 s结束。为了保持一致性，FLUENT采用RELAP5中的flibe熔盐物性。图5为flibe熔盐圆管流动的模型示意图，分别采用RELAP5、RELAP5/FLUENT耦合程序及FLUENT对圆管进行模拟。

2) 结果与讨论

RELAP5、RELAP5/FLUENT耦合程序及FLUENT在圆管中耦合界面1和耦合界面2处的质量流量和温度变化分别如图6、7所示。在整个瞬态模拟分析中，RELAP5/FLUENT耦合程序质量流量先保持不变，然后线性增加到1.5 kg/s后保持约2 s不变，最后线性减少到1.0 kg/s并保持不变。其变化趋势符合预期，也与RELAP5和FLUENT单独计算的结果保持了较好的一致性。同时，耦合程序的流体温度先保持873.15 K不变；当圆管入口处温度线性

图5 flibe熔盐圆管流动的模型示意图Fig.5 Model schematic of flibe salt flow in round tube

图6 耦合界面1(a)和耦合界面2(b)的质量流量Fig.6 Mass flow rate at coupled-boundary 1 (a) and coupled-boundary 2 (b)

图7 耦合界面1(a)和耦合界面2(b)的流体温度Fig.7 Fluid temperature at coupled-boundary 1 (a) and coupled-boundary 2 (b)

增加后，由于温度热传导和流体流动需一定时间，流体温度在耦合界面1和耦合界面2位置处会在t=4.1 s时刻温度开始线性增加到893.15 K，并保持893.15 K不变。图8为在耦合界面1和耦合界面2处RELAP5/FLUENT耦合程序与RELAP5、FLUNET的温度差值随时间变化的情况，其中，C-R表示RELAP5/FLUENT耦合程序和RELAP5的温度差值；C-F表示RELAP5/FLUENT耦合程序和FLUENT的温度差值。RELAP5/FLUENT耦合程序与RELAP5的温度差值均小于0.01 K，与FLUENT的温度差值也均小于0.30 K，相对误差较小。其变化趋势符合预期，与RELAP5和FLUENT单独计算的结果保持了很好的一致性。

图8 耦合界面处耦合程序与RELAP5、 FLUNET的温度差值Fig.8 Temperature difference between coupled program and RELAP5, FLUNET at coupled-boundary 1 and coupled-boundary 2

RELAP5、RELAP5/FLUENT耦合程序及FLUENT在耦合界面1和耦合界面2位置处的管道中流体压力变化分别如图9所示。在t=2～4 s过程中，由于flibe熔盐流速增加，管道中流体压降变大，在出口压力不变的情况下，管道中流体的压力会逐渐变大；在t=4～6 s过程中flibe熔盐温度增加，黏性降低导致压降变小，在耦合界面1和耦合界面2处中流体的压力会降低；在t=6～8 s过程中，由于flibe熔盐流速降低，管道中流体压降变小，管道中流体的压力会逐渐变小。RELAP5、RELAP5/FLUENT耦合程序和FLUENT在不同时刻耦合界面1和耦合界面2之间的压降列于表1。RELAP5/FLUENT耦合程序与FLUENT计算的压降比较相符，在t=4 s时最大相对误差为2.0%，在可接受的范围内。RELAP5/FLUENT耦合程序和FLUENT计算的压降，与RELAP5计算的压降差距较大，这是由于RELAP5基于流体充分发展情况下计算流体压降而RELAP5/FLUENT耦合程序和FLUENT采用低Re湍流模型，在进出口效应的影响下，耦合界面1处的流体尚未充分发展，导致压降计算相对偏大。总体而言，RELAP5/FLUENT耦合程序与RELAP5、FLUENT计算的压力变化比较一致。

2.2 2 MW熔盐堆稳态模拟

1) 2 MW熔盐堆简介

液态燃料熔盐堆基于在线干法处理技术,可实现钍铀增殖,是一种重要的研究堆型[15]。本文以一种2 MW液态燃料熔盐堆的初步概念设计为研究对象。该液态熔盐堆设计热功率为2 MW，堆芯采用石墨作为慢化剂，热工水力设计堆本体熔盐的进口温度为873.15 K，出口温度为893.15 K。一回路采用含有高富集度7Li的LiF-BeF2-ThF4-UF4熔盐作为燃料盐，二回路采用FLiNaK作为载热剂，一回路质量流量为55.4 kg/s。在利用RELAP5/FLUENT耦合程序模拟分析时，堆芯部分选取了1/4模型进行CFD建模，堆芯外的系统回路采用RELAP5建模，如图10a所示。下腔室设计采用喇叭状结构而流量分配装置处于下腔室正中心，其厚度为18 mm，其结构如图10b所示。堆芯熔盐通道1/4模型由图11a所示的堆芯石墨组件构成192个熔盐通道，熔盐通道的分布如图11b所示。由于实验堆堆芯几何结构较为复杂，对实验堆几何模型进行了简化处理，CFD建模计算时省略了流量分配装置支撑连接结构。熔盐堆部分材料的物性参数列于表2。在CFD模型中，对几何结构较为规则的堆芯活性区，采用六面体结构化网格进行划分；对于堆芯上下熔盐腔室为椭球形结构，采用非结构化四面体网格进行划分。网格质量最差为0.45，约85%的网格质量控制在0.8以上。网格无关性分析表明，当计算模型网格数量到达3 000万及以上时，可满足CFD精确计算要求。

图9 耦合界面1(a)和耦合界面2(b)的流体压力Fig.9 Fluid pressure at coupled-boundary 1 (a) and coupled-boundary 2 (b)

表1 熔盐在耦合界面1和耦合界面2之间的压降Table 1 Salt pressure drop between coupled-boundary 1 and coupled-boundary 2

2) 结果分析

利用RELAP5/FLUENT耦合程序对2 MW熔盐堆系统进行分析，相对单独的RELAP5系统分析结果而言，耦合程序可通过FLUENT模型计算获得堆芯更加详细的温度分布和流场分布。堆芯的温度分布如图12所示，熔盐在流经堆芯的过程中，温度逐渐上升并在堆芯活化区出口处达到最大值约905 K，随后熔盐在上腔室进行混合，在出口处温度约893 K。

图10 熔盐堆模型示意图(a)与流量分配装置结构(b)Fig.10 Schematic of molten salt reactor (a) and structure of flow distribution device (b)

图11 石墨组件(a)和熔盐通道分布示意图(b)Fig.11 Schematic of graphite assembly (a) and salt channel distribution (b)

表2 物性参数Table 2 Physical parameter

图12 堆芯温度分布示意图Fig.12 Temperature distribution of reactor core

熔盐堆下腔室流线示意图如图13所示，熔盐由较窄的进口管进入较宽阔的下腔室，熔盐经历一突然扩展的流态，流线形状近似于喇叭状结构，有利于熔盐的流动，避免了熔盐在下腔室的涡流现象。下腔室设置流量分配板，有效抑制了熔盐在下腔室中心区域的流速，使得熔盐流体不是直接由进口冲击中心位置及其附近熔盐通道，而是流经分配板和其孔道向四周扩散流动。因而下腔室内由涡流产生的流动死区基本得到消除，在堆芯径向流量分配相对均匀。

在堆芯活性区不同高度(z=0.0，0.4，0.8，1.1 m)下堆芯截面温度分布如图14所示。由于下腔室喇叭状结构和流量分配板的作用，熔盐在进入堆芯活性区(z=0.0 m处)各通道的流速和温度相对分布均匀。熔盐在流经堆芯活性区过程中，熔盐温度逐渐上升。在同一高度下，堆芯中心区域熔盐温度略高于堆芯外围熔盐温度。堆芯石墨受到熔盐加热，温度会略比通道中的熔盐温度高。在z=0.0 m处堆芯石墨温度约在878～882 K，随着堆芯熔盐加热和石墨自身功率加热，堆芯石墨温度逐渐上升，在z=1.1 m熔盐出口处，石墨温度达到894～906 K。

a——0.0 m；b——0.4 m；c——0.8 m，d——1.1 m图14 在不同高度下堆芯活性区截面温度分布Fig.14 Temperature distribution of active core at different heights

2.3 2 MW熔盐堆功率突变瞬态分析

功率突变瞬态分析是在上述稳态分析的基础上进行的。假定在t=0 s时堆芯功率突然增加10%并保持不变，熔盐堆系统无其他安全设施引入。考虑到熔盐堆系统的复杂性及计算机计算性能的限制，整个瞬态分析过程中时间步长设置为0.01 s，模拟时间为200 s。突然引入的功率变化会导致堆芯和系统的温度分布发生变化。耦合程序模拟分析的系统质量流量和堆芯进出口温度随时间的变化如图15所示。在瞬态计算中，一回路系统的质量流量为55.4 kg/s保持不变；由于熔盐具有较大的比热容及在堆系统中的滞留效应，堆芯进口温度在前期基本不变，出口温度逐渐上升。在t=20 s时刻，堆芯出口的平均温度为896.05 K，堆芯进出口温差为22.9 K。由于二回路换热能力有限，流出堆芯的熔盐温度逐渐上升后，在流经换热器冷却后温度会高于873.15 K，且使得堆芯入口温度在约t=30 s后开始缓慢上升，在t=200 s时堆芯进口温度达到875.0 K。

图16为瞬态计算过程中，在t=5、10、15、20 s时刻堆芯的温度分布。在初始时刻稳态情况下，堆芯内部最高温度约906 K。由于堆芯功率的增加，堆芯内部的最高温度也发生变化。在t=5、10、15、20 s时刻，堆芯内部最高温度分别达到909.0、910.4、910.7和911.0 K。初始时刻，堆芯下腔室熔盐平均温度在874 K，局部区域平均温度会达到875 K；在t=20 s时刻，由于堆芯功率增加，下腔室熔盐平均温度增加，局部区域温度可达879 K；堆芯内石墨温度也相对升高，靠近堆芯上腔室附近的石墨平均温度可达903 K。

图15 系统质量流量和堆芯进出口温度的变化Fig.15 Mass flow rate of system and temperature of reactor core inlet and outlet

a——5 s；b——10 s；c——15 s；d——20 s图16 不同时刻堆芯截面温度分布Fig.16 Temperature distribution of reactor core at different time

3 结论

以RELAP5与FLUENT程序为基础，利用对RELAP5源代码的二次开发和FLUENT的用户自定义函数进行编程，采用显式耦合的方式，开发了RELAP5/FLUENT耦合程序，实现两个代码之间数据交换。RELAP5在每个时间步进行瞬时计算后，返回FLUENT所需的边界条件参数，FLUENT在该时间步结束后返回RELAP5下一步计算所需的边界条件参数，重复此过程，直到达到预定问题的计算时间为止。利用RELAP5/FLUENT耦合程序对2 MW熔盐堆系统进行数值模拟分析，在进行系统分析的同时可获得堆芯详细的三维热工水力效应。该耦合程序可适用于包括flibe熔盐等多种流体的管道流动情况下的热工水力分析；在进行熔盐堆系统分析中，堆芯及上下腔室等存在三维复杂结构，局部区域存在明显的三维流动或需要获得三维温度分布，也可基于RELAP5/FLUENT耦合程序进行熔盐堆的热工水力分析。