反应堆多物理场耦合计算中的不确定性研究概述

Review of Uncertainty Research in Reactor Multi-Physics Coupling Calculation

摘要: 核反应堆是一个由粒子场、温度场、流体场、应力场等多个物理过程相互紧密耦合的系统，先进核反应堆数值模拟需要解决多尺度、大规模、多物理场耦合问题。近年来，随着计算机技术和现代数值仿真技术的进步，使得通过多物理场、多尺度等耦合框架实现对反应堆综合性能的大规模数值模拟成为可能，数值反应堆乃至反应堆数字孪生系统成为数值模拟的总体目标。随着最佳估算加不确定性分析方法的推广，不确定性分析成为核反应堆设计和安全分析的重要组成部分。相比起针对单一程序或单一物理场的不确定性分析，对于多物理场耦合计算的不确定性分析则更为复杂，除了各个物理场模拟本身的不确定性外，还必须考虑程序耦合过程中引入的不确定性，该方面的研究尚不够深入。本文从单物理场计算程序的不确定性源及其分析方法出发，对目前已开展的耦合程序计算不确定性分析相关工作进行了总结和阐述，并对多物理场耦合计算中的混合不确定性研究进行了初步分析。

Abstract: Nuclear reactor is a system which is closely coupled with many physical processes such as particle field, temperature field, flow field, and stress field. The numerical simulation of advanced nuclear reactor needs to solve the problems of multi-scale, large-scale and multi-physical field coupling. In recent years, with the progress of computer technology and modern numerical simulation technology, it is possible to realize large-scale numerical simulation of comprehensive reactor performance through multi-physical field and multi-scale coupling framework, numerical reactor and even reactor digital twin system become the overall target of numerical simulation. With the popularization of the method of best estimate plus uncertainty analysis, uncertainty analysis has become an important part of nuclear reactor design and safety analysis. Compared with the uncertainty analysis for a single code or a single physical field, the uncertainty analysis for the multi-physical field coupling calculation is more complex. In addition to the uncertainty of each physical field simulation itself, the uncertainty caused by the code coupling process must be considered, which has not been studied deeply enough. Based on the uncertainty sources and analysis methods of single physical field computing code, this paper summarizes and expounds the existing work on the uncertainty analysis of coupling codes, and makes a preliminary analysis on the mixed uncertainty research of multi-physics coupling calculation.

1. 引言

核反应堆是一个涉及到中子物理、热工流体、力学等多学科的复杂系统，随着核反应堆工程和技术的不断创新发展，其对安全性、可靠性、经济性和可持续性提出了更高的要求。尽可能准确地模拟核电厂设计基准事故及超设计基准事故是核安全分析的核心目标之一。目前，许多数值模拟程序能够对不同类型的核反应堆进行广泛的事故模拟和安全分析，例如RELAP5 (Reactor Excursion and Leak Analysis Program) [1] 、TRACE (TRAC/RELAP Advanced Computational Engine) [2] 、PARCS (Purdue Advanced Reactor Core Simulator) [3] 、COBRA (Coolant-Boiling in Rod Arrays) [4] 等。

近年来，随着计算机技术和现代数值仿真方法的进步，通过多物理场、多尺度等耦合框架实现对反应堆综合性能的高保真数值模拟成为可能。数值反应堆诞生在这样的背景下，它结合了核反应堆物理、热工水力、材料、力学、软件工程、算法等多个学科，以多物理耦合计算、高性能并行计算为基础，应用多种先进的计算程序，对反应堆内真实的物理过程进行全面精细化建模、分析和计算。20世纪以来，美国和欧盟先后启动了CONVERGENCE、NURESIM系列 [5] [6] 、CORTEX (CORe Monitoring Techniques and Experimental Validation and Demonstration) [7] 、HPMC (High Performance Monte Carlo Methods for Core Analysis) [7] 、McSAFE (High Performance Monte Carlo Methods for Safety Analysis) [7] 、CASL (Consortium for Advanced Simulation of LWRs) [8] 、NEAMS (Nuclear Energy Advanced Modeling and Simulation) [9] 、GEMINA (Generating Electricity Managed by Intelligent Nuclear Assets) [10] 等多个数值反应堆研发项目，在“数值反应堆原型系统开发及示范应用”项目支持下，我国也开展了虚拟反应堆CVR1.0 [11] 的开发，如图1所示。除此之外，中广核研究院、中国核动力研究设计院等相关单位也开展了较多数值耦合程序的相关研究，例如作为数值反应堆中必不可少的物理和热工部分中广核研究院有限公司开发了三维物理热工耦合分析软件 [12] ，中国核动力研究设计院开发了数字化反应堆高保真中子学程序SHARK [13] 。

Figure 1. Overview of numerical reactor R & D projects [14]

图1. 数值反应堆研发项目概况 [14]

由于认知的局限性和真实系统的复杂性，任何建模与模拟计算都存在不确定性，对其进行量化是所有数值模拟技术必须解决的重要问题，其中数值反应堆耦合系统的不确定性量化研究已经逐渐成为国际上的研究热点 [15] 。对于耦合程序而言，由于每个物理场都有各自不同的不确定性来源，且耦合过程，例如程序之间参数的传递、耦合接口的应用、数值迭代的匹配等也会引入不确定性。不确定性源的识别与量化及其在核反应堆多物理场中的传递机理是一个涉及多学科交叉的复杂问题。因此对数值反应堆多物理耦合计算进行不确定性的量化分析是保障耦合计算结果可靠性的重要基础。

本文从单物理场程序计算的不确定性源及不确定性分析方法出发，对目前国内外已开展的多物理场耦合程序不确定性分析相关工作进行了介绍。此外，结合已开展的不确定性分析工作，本文对未来多物理场耦合计算的不确定性研究进行了初步分析和展望，并提出了一种针对混合不确定性的分析方法。

2. 单物理场程序不确定性研究

目前国内外针对单个计算程序进行不确定性分析的研究较多，不确定性分析包括“识别、描述相关输入参数(即不确定性源项分析)”与“量化仿真计算结果的不确定度”两部分。本节将主要介绍单物理场程序计算时的不确定性来源及其不确定性分析方法。

2.1. 不确定性来源

对于热工水力计算而言，其不确定性来源主要总结为以下几项 [16] ：程序与模型的不确定性；电厂数据的不确定性；比例的不确定性；用户效应的不确定性。

而对于反应堆物理计算而言，其不确定性来源也主要总结为以下几项 [17] ：建模的不确定性；数值计算方法的不确定性；电厂数据的不确定性。

此外，国际上亦存在一种常用的定性分类方式，可分为随机不确定性(Aleatory Uncertainty)、认知不确定性(Epistemic Uncertainty)、数值不确定性(Numerical Uncertainty)，部分研究中将数值不确定性归类于认知不确定性，如图2所示。在很多实际情况下，通常随机不确定性与认知不确定性是同时存在的，此时被称为混合不确定性。

Figure 2. Classification of uncertainty

图2. 不确定性分类

2.2. 不确定性分析方法

不确定性分析方法中最重要的一步就是对上述不确定性来源进行选择和量化。选择指的是通过专家判断，如现象识别及排序表(Phenomena Identification and Ranking Table, PIRT)或敏感性分析确认；量化是需要给出选择的输入不确定性参数的概率分布或区间分布。

不确定性量化分析是指评估输入不确定性对模型输出的影响。如图3所示，不确定性分析方法若依据采用的数学理论的不同，可以分为统计性分析方法以及确定性分析方法 [18] ，两者之间最大的区别在于：统计性分析方法通常是先进行不确定性计算，再进行敏感性计算。而确定性分析方法通常先进行敏感性计算以获取每个输入参数的敏感性因子，然后进行不确定性计算。根据不确定性传播途径的不同，统计性不确定性分析方法又可以分为：基于输入参数不确定性传播方法和基于输出参数不确定性外推方法。目前国际上主要采用的方法是基于输入参数不确定性传播方法。

Figure 3. Classification of uncertainty analysis methods

图3. 不确定性分析方法分类

近年来，确定性分析方法多用于反应堆物理方面的不确定性分析研究工作，主要包括微扰理论和直接数值扰动方法等 [19] 。微扰理论方法首先需要求得积分参量对于核截面数据的灵敏度函数，然后结合协方差矩阵求得积分参数的不确定性，分为经典微扰理论(CPT)和广义微扰理论(GPT)。直接数值扰动的计算思路与微扰理论大致相同，首先进行敏感性分析，每次对一个输入参数进行微小扰动，观察输出结果的变化，然后该输入参数返回其初始值，接着对另一个输入参数进行微小扰动观察输出结果的变化，直到所研究的所有输入参数都经过微小扰动且观察其输出结果，此方法适用于输入参数较少的情况。基于此，国内外研发了大量的反应堆物理不确定性分析计算程序，具体如表1所示。

目前，基于实验数据的逆不确定性量化(Inverse Uncertainty Quantification, IUQ)分析方法也得到广泛的关注。对于输入参数的实验数据无法测量得到的情况，可以通过比较程序输出的计算值与对应的实验测量值，反向量化得到输入参数的不确定性分布。近几年来一些研究机构基于反向求解结合统计学算法提出了本构模型的IUQ分析方法，如法国AREVA、EDF、IRSN和CEA开发的用于CATHARE程序本构模型的CIRCE [20] 方法，韩国KAERI开发的MCDA [21] 方法，意大利比萨大学研究机构GRNSPG-UNIPI开发的FFTBM [22] 方法，法国IRSN基于CATHARE程序开发的DIPE [23] 方法，比利时Tractebel公司基于DAKOTA程序开发的Tractebel IUQ [24] 方法等。

Table 1. Reactor physics calculation uncertainty analysis code

表1. 反应堆物理计算不确定性分析程序

3. 多物理场耦合程序不确定性研究

在反应堆系统中，中子物理–热工水力耦合(核热耦合)是最受关注的多物理场耦合形式。广义上的热工水力学计算程序按照其计算尺度可以分为三种 [41] ：热工水力系统程序(系统级)、热工水力子通道程序(部件级)、计算流体力学程序(Computational Fluid Dynamics, CFD) (局部级)。本节将主要介绍目前国内外已开展的耦合程序不确定性研究，并总结了程序使用的不确定性分析方法，以及选取的不确定源的异同。

3.1. 国内外研究现状

国外对采用热工水力系统程序与中子物理程序进行耦合来开展不确定性的研究相比于国内发展的较早且较多，如RELAP5、TRACE等与PARCS、BARS等进行耦合。国外学者已经利用TRACE-PARCS、RELAP5-PARCS等耦合程序进行了多次不确定性研究，国内相对而言研究较少。除此之外，国外一些学者还采用RELAP5-BARS、KIKO3D-TRABCO、ATHLET-PARCS、ATHLET-DYN3D、APOLLO2-FLICA4-CRONOS2、APOLLO3-FLICA4-ALCYONE、FLICA4-CRONOS等耦合程序进行了不确定性研究。具体的耦合程序、不确定性分析方法以及所考虑的不确定性源如表2所示。

Table 2. Thermal hydraulic system code-neutron physics code coupling

表2. 热工水力系统程序–中子物理程序耦合

除了采用热工水力系统程序与中子物理程序进行耦合，国内学者利用CFD与中子物理耦合程序也开展了不确定性研究。具体的耦合程序、不确定性分析方法以及不确定性源如表3所示。

对于热工水力子通道程序的耦合不确定性研究而言，目前国内外开展的研究还并不多。2022年美国的Delipei [63] 等人采用CTF-PARCS耦合程序对UAM-LWR基准题的阶段三的三个练习进行计算，包括稳态计算、损耗计算和瞬态计算，开发的耦合程序嵌入在不确定性量化框架中执行不确定性分析。在该框架中，Sampler用于对多群微观截面、裂变产率和制造参数进行扰动，DAKOTA程序用于对CTF输入参数和边界条件进行采样，开发Python脚本用于自动化和模块化预处理和后处理，该框架允许通过中子学和热工水力学建模对输入进行一致的扰动。

从以上研究中可以看出，目前国际上通常采用热工水力程序与中子物理学程序进行耦合计算并开展不确定性研究。所研究的事故大多都是反应性引入事故(Reactivity Initiated Accident, RIA)，例如弹棒事故、控制棒失控抽出事故、蒸汽管线断裂事故和硼稀释事故等。这是由于在发生反应性引入事故的特殊情况下，

Table 3. CFD-neutron physics code coupling

表3. CFD-中子物理程序耦合

各种物理场之间存在较强烈的多物理耦合效应，因此需要进行多物理耦合分析来捕捉物理场之间的相互作用效应。所采用的不确定性分析方法大部分为CSAU、GRS等方法，这与对单物理场程序进行不确定性分析时所采用的方法并没有很大的不同之处。所考虑的不确定性来源仅局限在单物理场程序的输入参数、边界条件及几何尺寸参数等，并没有考虑由于程序耦合过程中引入的不确定性源。

3.2. 耦合过程引起的不确定性

单物理场的不确定性分析已有较可行的分析方法，而由于多物理场耦合过程所引起的不确定性分析则还缺乏深入的研究。以核热耦合为例，目前所考虑的不确定性来源大多局限在单物理场程序的输入参数、边界条件及几何尺寸参数等随机不确定性，但是程序在耦合过程中会引入一些额外的认知不确定性，如网格映射、程序时间步长调整等，且此类不确定性源在单一物理场计算中并不存在，因此较难识别与量化。程序耦合过程主要包含耦合方式、空间网格映射、耦合时间步长算法及耦合数值算法四个关键环节，这其中的每一个环节都有可能产生额外不确定性。本节将对这四个方面目前所用的技术以及可能产生的不确定性进行初步分析。

3.2.1. 耦合方式

如图4所示，根据数值求解的方法不同可将耦合方式分为强耦合(或称紧耦合)和弱耦合(或称松耦合)。强耦合是指将堆芯三维中子学方程与热工水力学基本方程进行联立求解，涉及程序的重新开发，而弱耦合是指将现有程序作为独立的求解单元，并通过相互传递数据实现耦合。其中弱耦合又可具体分为内部耦合和外部耦合，内部耦合是指将一个程序作为子程序嵌入到另一个程序中，不同程序之间通过内存访问进行数据的交换；外部耦合则是指在现有程序的基础上独立开发外部耦合接口进行耦合参数的传递。

根据数据传递的特点可将耦合方式分为串行耦合和并行耦合。串行耦合指的是中子物理程序或热工水力学程序在单个或多个CPU上依次运行计算，而并行耦合指的是程序在不同的CPU上独立运行，同步进行中子学与热工水力学计算，并在计算过程中完成耦合参数的传递与交换，目前通常采用并行虚拟机(PVM)或消息传递接口(MPI)来进行数据的交换。

不同的程序耦合方式会带来不确定性，例如在采用强耦合方式进行程序耦合时，需将不同程序的线性方程进行联立求解，目前国际上通常采用牛顿迭代法和Krylov子空间方法求解关于雅可比矩阵的线性方程，在利用不同的求解方法时会产生不确定性。

Figure 4. Different coupling methods bring uncertainty

图4. 不同耦合方式带来不确定性

3.2.2. 空间网格映射

在耦合程序中，中子物理程序和热工水力程序各有一套网格。这两套网格通过网格映射相互传递参数，中子物理程序计算出的堆芯功率分布映射到热工水力程序网格中用于计算热工参数，热工水力学计算则为中子物理学程序提供燃料温度场、密度场、冷却剂温度场等截面反馈参数，如图5所示。

在进行空间网格划分时，中子物理学程序与热工水力学程序往往采用不同的网格划分方式，因此需要将两个程序的空间网格进行一一对应。理论上来说网格划分越精细计算所得到的结果越精确，但是网格划分过细在计算过程中需要消耗大量的计算时间及计算资源。所以在实现耦合程序的空间网格对应时，需要依据精度需求以及计算能力做出综合考虑。不同的网格映射方案会引入不同程度的不确定性。

一般情况下，中子物理计算的空间网格要比堆芯热工计算(例如子通道)的网格要更为精细，因此在进行网格映射的参数传递时，需要进行一定程度的平均化处理。一般来说这种精细度差异的程度越大，所引起的不确定性越高。

3.2.3. 耦合时间步长算法

程序时间步长的选取主要包括中子物理与热工水力程序各自计算的时间步长选取与耦合参数传递的时间点的选取。由于中子物理程序与热工水力学程序都有自己的时间步长选取算法，通常与它们各自的计算精度有关，中子物理学计算通常采用很小的时间步长，一般为毫秒级，但热工水力学计算过程较慢，通常采用较大的时间步长，因此时间步长上的耦合及同步是十分重要的。当时间步长过大时，细微的瞬态现象就不能被捕捉，使得计算过程不容易收敛，而时间步长过小时，则会消耗大量的计算时间。因此，与精细化程度不同的空间网格类似，不同的时间步长选取算法在耦合程序计算过程中也会带来不同程度的不确定性，如图6所示。

在耦合过程中，由于时间步长的差异，热工水力学计算所接收到的中子计算结果通常是时域平均的，而中子物理计算从热工水力学计算收到的反馈参数则要在一定时域内进行归一化处理，这些过程可能引入更多的不确定性。

Figure 5. Different grid mapping schemes bring uncertainty

图5. 不同网格映射方案带来不确定性

Figure 6. Different time steps bring uncertainty

图6. 不同时间步长带来不确定性

3.2.4. 耦合数值算法

耦合计算通常分为显式和隐式两种数值算法，显式耦合在每个时间步结束后不需要进行迭代，无需判断其收敛性而直接进入下一时间步计算，这种方法计算速度快但稳定性极差；隐式耦合需在每个时间步结束后通过反复迭代直到满足功率或温度的收敛准则再开始下一时间步的计算，如图7所示。

显式耦合数值算法计算速度快但稳定性极差，而隐式耦合数值算法计算得到的结果更加稳定但需要耗费大量的时间，因此不同的耦合数值算法也会带来不确定性。例如，假设中子物理计算采用隐式算法，而与之耦合的热工水力学计算采取显示算法，则在耦合算法中如何处理，也会影响到结果的不确定性。

Figure 7. Different numerical algorithms bring uncertainty

图7. 不同数值算法带来不确定性

4. 混合不确定性分析方法

4.1. 不确定性的量化

在对多物理场耦合计算进行不确定性分析时，通常随机不确定性与认知不确定性是同时存在的，此时被称为混合不确定性。随机不确定性被认为是系统的固有随机性，常用概率论中的概率模型对其进行不确定性量化，通常分为前向和反向不确定性量化方法。前向不确定性量化方法的原理相对简单，只要存在足够的实验数据即可得到十分理想的结果，而在实际情况中，可能存在部分模型输入或输出参数无法直接测量得到的情况，此时难以使用前向量化方法，反向不确定性量化方法则可适用于该情况，对于输入参数实验数据无法测量得到的情况，可通过比较模型输出的计算值与对应的实验测量值，反向量化得到输入参数的不确定性分布，贝叶斯校准则是典型的反向不确定性量化方法。

而认知不确定性是由于知识的不完备以及数据的缺乏造成的，此时基于概率的方法将不再适用，可采用一些非概率模型来度量，如凸集模型 [64] 、模糊集模型 [65] 、证据模型 [66] 等。

针对随机和认知同时存在的混合不确定性，学术界提出了概率盒(Probability Box, P-Box)模型 [67] 来对其进行量化，使用上下两条累积分布函数(Cumulative Distribution Function, CDF)作为边界来对混合不确定性进行量化，并进一步结合不确定性传播方法进行不确定性传播分析。

4.2. 双层抽样方法

对于混合不确定性的传播可采用双层抽样方法 [68] 来将输入参数的不确定性传递至输出响应，外层为参数层采样，用于获得区间分布参数的样本；内层为概率层采样，用于计算输出响应的概率或者统计矩。其大致思路如下：将认知不确定性归为外层，将随机不确定性归为内层，先进行一次的外层抽样再进行N次的内层抽样，此时获得输出响应的一条累积分布函数，如此重复直到进行了M次外层抽样，最终可以获得输出响应的M条累积分布函数，即输出响应的一个P-Box，也就将输入参数的不确定性传递至了输出。

4.3. 混合不确定性分析方法流程

对多物理场耦合系统进行不确定性分析，在执行具体的不确定性传播时，可采用双层嵌套抽样方法将输入参数的不确定性传递至输出响应，将认知不确定性归为外层，随机不确定性归为内层，如图8所示，具体流程如下：

1) 第一进行不确定性源的识别，通过敏感性分析及PIRT表、专家经验等识别出重要的输入参数；

2) 第二进行不确定性的量化，随机不确定性采用贝叶斯校准方法进行反向量化，用不确定性区间进行表征，认知不确定性采用证据理论进行量化，用概率盒进行表征；

3) 第三进行不确定性的传播计算，将随机不确定性与认知不确定性分为内外两层，采用双层抽样方法将输入不确定性传递至输出响应。

Figure 8. Mixed uncertainty analysis methodbased on double layer sampling

图8. 基于双层抽样的混合不确定性分析方法

5. 总结

本文首先对单物理场程序计算的不确定性源及不确定性分析方法进行了简要概述，并介绍了目前国际上对多物理场耦合程序不确定性分析的研究现状，最后对未来耦合程序不确定性研究进行了初步分析，得出以下结论：

1) 当前对耦合程序不确定性的研究所采用的耦合程序大多是核热耦合计算程序，如热工水力程序(RELAP5、TRACE、CFD、CTF等)与中子物理学程序(PARCS、BARS等)进行耦合，所研究的事故大部分为反应性引入事故。

2) 在对耦合程序进行不确定性分析时，目前所做研究中所采用的不确定性分析方法大部分为CSAU、GRS等方法，该类方法所考虑的不确定性来源仅局限在单物理场程序计算时所引入的随机不确定性，但是程序耦合过程中所引入的认知不确定性，例如网格映射、程序时间步长调整等，目前尚未有充分的研究，未来应该对此类不确定性源进一步识别与量化。

3) 在对耦合程序进行不确定性分析时，通常随机不确定性与认知不确定性是同时存在的，对于这种混合不确定性，可以采用基于双层抽样的混合不确定性分析方法：首先通过PIRT表或敏感性分析方法进行不确定性源的识别；对于随机不确定性采用贝叶斯校准等反向方法进行量化，对于认知不确定性采用证据理论等方法进行量化；然后通过双层嵌套抽样方法将输入参数的不确定性传播至输出响应，将认知不确定性归为外层，随机不确定性归为内层。

致谢

感谢王诗琦老师对文章英文摘要写作的指导。

基金项目

系统级耦合程序不确定性评估方法及验证研究(IDR2-05-FW-HT-20220006)。

NOTES

^*通讯作者。