基于统一耦合框架的堆芯物理热工耦合程序的开发及验证

姜 荣， 冯文培， 陈红丽， 强胜龙， 李治刚， 潘俊杰， 张喜林， 罗 晓

(1. 中国科学技术大学核科学技术学院， 合肥 230027； 2. 中国核动力研究设计院， 成都 610213)

1 引 言

反应堆系统内中子物理、热工水力、燃料力学等多物理过程间存在复杂的耦合关系. 为了实现反应堆多物理多尺度耦合，提高模拟计算的精确度，扩大计算规模，减少重复性工作[1]，基于统一耦合框架的多物理程序耦合成为研究热点. 这些耦合框架或耦合平台包括ICoCo、MOOSE平台、OpenPALM、preCICE等，其中通过ICoCo理念封装实现多物理程序的耦合方法在国外已开展较多研究.

法国CEA基于ICoCo理念实现了中子学程序CRONOS2、两相流分析程序FLICA4、系统分析程序CATHARE程序间的耦合[2]. 德国学者在NURESAFE项目框架下实现了子通道程序CTF和中子程序DYN3D的耦合[3]. 为了更好地模拟反应堆压力容器中的三维热工水力现象，德国KIT使用ICoCo接口实现了系统程序TRACE和CFD程序TrioCFD之间的耦合，该耦合系统已集成于SALOME平台[4]. 中国科学技术大学和德国KIT合作开发了CFD程序TrioCFD和子通道程序SubChanFlow之间基于ICoCo理念和MEDCoupling库的统一耦合框架的耦合[5]. 为了更好地模拟钠冷快堆热工水力现象，CEA实现了在ASTRID项目框架下系统程序CATHARE、子通道程序TrioMC和CFD程序TrioCFD之间的耦合[6]. 国内对于ICoCo封装理念的研究还较少，因此对其开展更多的研究从而为多物理多尺度耦合提供基础是有必要的.

本文基于ICoCo封装理念和MED库数据传递模型，对中子时空动力学程序CORCA-K和子通道分析程序CORTH进行了基于耦合框架的物理-热工耦合程序的开发. 采用NEACRP国际轻水堆弹棒基准[7]进行了程序测试，包括对封装后程序和耦合程序的测试，并与基准结果以及CORCA-K原有耦合接口的计算结果进行对比.

2 理论依据

ICoCo接口是在NURISP项目范围内[8]开发的面向对象的通用代码耦合接口，首先应用于CATHARE和TrioCFD的多尺度热工水力耦合[9]. ICoCo接口基于C++开发，通过定义Problem通用母类来定义函数，这些函数的功能包括初始化、时间步推进、保存和还原以及场数据传递，表1给出了ICoCo定义的重要接口及功能.

表1 重要的ICoCo接口及功能

基于ICoCo理念对程序进行封装，即对程序开发上述接口. 其中，在进行场数据交换时采用功能强大的网格和区域操作库MEDCoupling实现场数据的交互和映射. MEDCoupling场数据结构包含足够多的信息，可以进行各种插值操作，也可以在不同进程间通过并行模式和分布式模式交换数据结构信息. 进行ICoCo封装时需要对每个程序定义一套或多套MED网格，以进行输入场数据和输出场数据的传递和映射.

为了实现基于统一耦合框架的物理-热工程序的耦合，采用C++ 编写的supervisor程序实现对各子程序的调用. 图1给出了C++ supervisor进行程序调用的示意图.

图1 supervisor调用程序示意图

3 程序开发

3.1 CORCA-K的ICoCo封装

CORCA-K程序是中国核动力研究设计院自主研发的三维瞬态中子学计算软件，其核心功能是进行三维瞬态中子扩散方程的数值求解计算，空间离散采用第二类边界条件格林函数节块法，时间离散采用二级二阶对角线隐式龙格库塔格式[10]. CORTH程序是中国核动力研究设计院自主研发的子通道分析程序，可以进行反应堆热工水力设计与安全分析，也可以用来进行新型燃料组件的设计与研发[11].

基于统一耦合框架的堆芯物理程序与热工程序的耦合，需要对CORCA-K和CORTH进行ICoCo封装. 首先，对CORCA-K的源程序进行功能模块划分，如划分为启动模块、初始化模块、时间步计算模块、中止模块等. 然后根据功能模块开发表1所示的ICoCo接口.

在编写场数据交换与映射相关接口时，根据CORCA-K程序的网格特征，建立三维结构化MEDCoupling网格来进行数据交换. 同时，设置CORCA-K的输入物理场为有效燃料温度、慢化剂温度、冷却剂密度和燃料芯块温度，CORCA-K的输出物理场为功率分布.

完成CORCA-K的ICoCo封装后，编译得到CORCA-K的动态链接库，编写supervisor调用动态链接库，即可实现CORCA-K的计算功能. 按照同样的方法对CORTH进行封装并编译得到CORTH的动态链接库.

3.2 耦合实现

在得到封装后CORCA-K和CORTH的动态链接库后，编写supervisor程序，编译得到基于统一耦合框架的CORCA-K和CORTH的耦合程序，图2为该耦合程序进行计算的流程图.

图2 CORCA-K与CORTH耦合计算流程图Fig.2 The coupling calculation flowchart of CORCA-K and CORTH

图3表示的是CORCA-K和CORTH的MED网格. CORCA-K创建的是三维结构化网格，CORTH创建的是三维非结构化网格. 在进行数据传递时，CORCA-K计算得到堆芯三维功率分布，通过MED网格映射传给CORTH，CORTH将新的功率分布作为输入，重新计算得到平均多普勒燃料温度、慢化剂平均温度、最高燃料温度和冷却剂密度，再通过MED网格映射传给CORCA-K，CORCA-K重新计算得到功率分布，这样循环下去，实现二者的耦合.

CORTH非结构化网格

4 数值验证

4.1 测试算例

选择NEACRP轻水堆弹棒基准中的C1和C2为测试算例，对封装后的CORCA-K和基于耦合框架的耦合程序进行测试. 由OECD/NEA发布的轻水堆弹棒基准，即在热态零功率和热态满功率下，堆芯中一组控制棒快速弹出造成正反应性引入. NEACRP基准中堆芯由157组燃料组件和64组反射组件组成，其中49组燃料组件中布置有控制棒，各组件的径向尺寸是21.606 cm×21.606 cm，轴向总长度为427.3 cm，如图4所示. 瞬态分析计算5 s，0.1 s时发生弹棒事故，0～1 s内，时间步长为0.005 s，1～5 s内，时间步长为0.05 s. 表2给出了NEACRP弹棒基准中C1和C2运行数据的描述.

4.2 封装后CORCA-K测试

通过编写supervisor程序检验每个ICoCo接口是否可以正常调用和实现其功能，从而对封装后的CORCA-K进行测试. 以C1和C2为研究对象，分别采用不同的计算模式如临界计算、弹棒计算等进行计算，全面验证封装后CORCA-K的功能. 经过多次验证计算，对于相同的算例和相同的计算模式，封装后CORCA-K与封装前CORCA-K计算得到的输出卡内容完全一致. 由此可见，CORCA-K的封装是正确的. 由于封装前后输出卡完全一致，因此这里不列出具体的数据.

图4 轻水堆基准C1和C2径向堆芯布置图

表2 C1和C2的堆芯运行数据

4.3 耦合程序测试

使用基于耦合框架的CORCA-K和CORTH耦合程序计算C1和C2，并使用CORCA-K已开发的与CORTH直接耦合的接口CORTHV2计算相同的算例作为对比. 由于CORTHV2耦合接口计算的结果与基准参考解已经存在偏差，基于耦合框架的CORCA-K与CORTH的耦合程序是在CORCA-K和CORTH两程序的基础上开发，因此为验证耦合框架的耦合方式的准确性，同时排除CORCA-K和CORTH本身的计算误差，主要对比基于耦合框架的CORCA-K与CORTH的耦合程序与CORTHV2的计算结果. NEACRP基准的PANTHER修正后的计算结果[12]作为参考.

表3、表4给出了使用两种耦合方式计算C1和C2得到的稳态和瞬态计算结果，下文以“ICoCo”代表基于耦合框架的耦合程序，“CORTHV2”代表CORCA-K原来的耦合接口. 从表3中可知，C1使用两种方式进行稳态计算得到的临界硼浓度均为1 128.79 ppm，与基准结果仅差0.49 ppm，C2使用两种方式计算得到的临界硼浓度均为1 154.73 pcm，与基准结果仅差1.87 pcm. 一方面可以验证ICoCo和CORTHV2两种耦合方式进行稳态计算的准确性，另一方面也可以验证封装后CORCA-K进行临界计算的准确性.

C1中，两种耦合方式最大功率出现的时间仅差0.002 5 s，最大功率和5 s时的功率相差0.228 9，最大功率相差0.23. 考虑到CORCA-K与CORTH之间的数据传递需要经过MED网格的交互与映射，因此可能会造成一定误差，之后可以进行网格敏感性研究，进一步确定误差范围. 对于5 s时的功率，两种耦合方式计算的结果相同，均为0.16. 对于5 s时的平均多普勒燃料温度，二者的计算结果也非常接近，最大不超过0.5 ℃，燃料最高温度相差不超过3 ℃. 同时，两种耦合方式的结果与基准结果也基本符合. C2中，两种耦合方式的计算结果更加一致. 二者最大功率到达时间仅差0.002 5 s，最大功率及5 s时的功率完全一致(1.07和1.03)，5 s时的平均多普勒燃料温度和燃料最高温度相差均不超过1 ℃.

表3 不同耦合方式计算得到的C1弹棒基准的稳态和瞬态结果

图5～图10分别给出了使用两种耦合方式计算得到的相对功率变化、平均多普勒燃料温度、最高燃料温度的变化趋势. 由图中可见，最大功率、最大功率出现的时间、5 s时的功率、平均多普勒燃料温度和燃料最高温度的变化趋势二者皆吻合得很好，与PANTHER的结果略有差距，但相差很小.

结合以上图表和分析可知，基于统一耦合框架的堆芯物理-热工耦合程序的计算结果是准确的.

表4 不同耦合方式计算得到的C2弹棒基准的稳态和瞬态结果

图5 ICoCo和CORTHV2耦合模式计算得到的C1功率变化对比

图6 ICoCo和CORTHV2耦合模式计算得到的C1平均多普勒燃料温度对比

图7 ICoCo和CORTHV2耦合模式计算得到的C1最高燃料温度对比

图9 ICoCo和CORTHV2耦合模式计算得到的C2平均多普勒燃料温度对比

图10 ICoCo和CORTHV2耦合模式计算得到的C2最高燃料温度对比

5 结 论

(1) 本文基于ICoCo统一封装理念和MED库统一数据传递模型，分别对中子时空动力学程序CORCA-K和子通道分析程序CORTH进行了ICoCo封装. 然后基于统一耦合框架实现了CORCA-K和CORTH的耦合.

(2) 以NEACRP轻水堆弹棒基准中C1和C2为研究对象进行耦合程序的测试，并与CORCA-K原有的耦合接口CORTHV2计算的结果进行对比，以基准结果作为参照. 经过对比分析可知，基于统一耦合框架的物理热工耦合程序计算的结果与CORTHV2计算的结果符合地很好，各参数的变化趋势完全一致，由此证明基于统一耦合框架的物理-热工耦合程序的计算是正确的.

(3) 下一步可以进行耦合程序的敏感性测试，进一步确定其误差范围.