高温气冷堆全耦合系统直接联立求解的方法研究和程序开发

张 汉，郭 炯，邬颖杰，王毅箴，刘保坤，崔梦蕾，孔勃然，朱凯杰，刘礼勋，江 卓，窦沁榕，唐焕燃，李 富

(清华大学 核能与新能源技术研究院，先进反应堆工程与安全教育部重点实验室，北京 100084)

核电站系统的耦合计算是对现有反应堆各领域数值技术的融合、集成和进一步提升，是核能领域的前沿方向之一。完整的核电站系统同时具有多种特性不同的耦合方式，在堆芯存在由反应堆物理、热工水力构成的基于共用计算域的多物理耦合，堆芯、蒸汽发生器、主泵和连接管线构成了基于拓扑网络的一回路多部件耦合，一、二回路通过蒸汽发生器构成了基于换热界面的多回路耦合。多物理、多部件、多回路是所有类型反应堆共有的耦合特征。对于高温气冷堆核电站，在上述特征基础上，由于采用“多堆带一机”的运行模式而形成了更大规模、更复杂的多模块耦合；堆芯的宏观球床温度、介观燃料球内部温度与微观颗粒内部温度构成了多尺度耦合；同时由于采用响应速度快的螺旋管式直流蒸汽发生器而使一、二回路耦合更加紧密。因此，高温气冷堆核电站更复杂，是一个具有多物理、多尺度、多部件、多回路、多模块特征的非线性强耦合系统，要精确模拟核电站的动态特性，需同时求解以上多种耦合关系，更需要稳定、准确、高效的耦合求解。

算符分解和Picard迭代等传统耦合方法至今仍然被核工业界广泛应用，其基本思路是将耦合系统解耦、并通过传递耦合参数的方式考虑耦合效应，这样能充分利用已有计算模块，给出一些工程急需的结果。但对于更复杂、更大规模的完整核电站耦合系统，传统耦合方法受限于稳定性条件、线性收敛速度等本质缺陷，无法胜任。以Jacobian-free Newton-Krylov(JFNK)/Newton-Krylov(NK)为代表的计算性能更优越的直接联立耦合方法是最具潜力的计算框架。在美国能源部CASL/NEAMS等计划支持下，对直接联立计算方法进行了一系列研究，有力地证明了直接联立耦合方法的有效性[1-5]，但至今尚未实现在统一框架下的完整核电站全耦合系统直接联立求解。特别是MOOSE耦合平台，其内核采用了先进的JFNK耦合方法，但受限于核心算法和程序架构，被证明本质上只能用于求解同区域、单层级、单一方式的多物理耦合问题，对于堆芯级-部件级-系统级的多层级耦合结构问题则还需要在算法研究和程序架构设计上进一步突破[6-8]。

本文在综述国内外反应堆耦合计算研究进展的基础上，介绍清华大学核能与新能源技术研究院在高温气冷堆直接联立耦合方法及程序开发方面的研究工作，提出针对多尺度问题的非线性消去直接联立方法等关键技术，尝试给出具备同时解决不同耦合机制的统一算法框架，并进一步抽象为可描述多层级耦合结构的统一程序框架。

1 国内外研究现状

1.1 国外研究现状

反应堆耦合计算对于反应堆设计和安全分析至关重要，在美国能源部CASL/NEAMS等计划[1-2,9-10]和欧盟NURSIM/NURISP/NURESAFE系列计划[11-12]支持下，相关领域得到了快速发展。耦合计算方法逐渐从传统的算符分解、Picard迭代思路向更先进的所有物理场在统一耦合算法下直接联立求解的方向发展；同时模拟范围也尝试从单一堆芯多物理耦合问题向统一软件框架下的耦合部件更多、耦合关系更复杂的完整核电站全耦合系统的方向发展。美国已涌现出一系列基于MOOSE耦合平台[3-4]、采用先进直接联立方法的新一代先进反应堆耦合计算软件，如中子输运程序RATTLESNAKE[13]、系统热工水力程序RELAP7[14]、三维多孔介质热工水力程序PRONGHORN[15]、燃料性能分析程序BISON[16]等。这些基于MOOSE平台开发的耦合程序采用通用的建模方式，同时采用有限元方法完成耦合场自动数值离散，具有强大的自动化、标准化建模求解能力，同时采用先进的高效直接联立JFNK耦合方法以得到更高的计算性能。Park等[17]在MOOSE框架下开发了高温气冷堆堆芯物理-热工程序PRONGHORN，该工作根据对物理-热工现象的理解，在预处理矩阵构建中保留耦合的关键信息，忽略弱耦合的、难于显式构建的元素，在保证预处理效果的情况下减少构建代价，JFNK的计算效率是传统Picard迭代的6～8倍左右；在大时间步长下，Picard迭代可能会发散，但JFNK仍能快速收敛，具有更好的收敛性。Pope等[18]研究了压水堆瞬态多物理耦合计算问题，JFNK方法的自适应时间步长不受CFL数限制，仅根据计算精度进行确定，可达到数十、乃至数百倍CFL数，计算效率可更高。可见，以JFNK为代表的联立耦合方法在计算稳定性、计算精度、计算效率等方面展现出优越性。

但这些基于MOOSE开发的耦合计算模块本质上只能求解同区域、单层级、单一耦合方式的多物理问题，并不能解决多尺度、多部件和多回路耦合问题。2021年，加州大学伯克利分校联合爱达荷国家实验室尝试基于PRONGHORN程序解决熔盐球床堆的颗粒-燃料球-球床多尺度耦合问题，虽然每个子场都具备求解瞬态问题的能力，但受限于MOOSE平台对多尺度耦合问题描述能力不足，仅能通过解耦热源进行简化，计算较简单的稳态问题。即使采用较简单的Picard迭代方法，都没有实现需要考虑完整耦合关系的瞬态多尺度问题，更没有实现直接联立求解方法[7]。2015年爱达荷国家实验室尝试利用RELAP7、BISON、RATTLESNAKE模拟AP1000的多部件耦合问题[6]，2021年尝试利用PRONGHORN和RELAP7模拟高温气冷堆多部件问题[8]，同样受到MOOSE平台的固有限制，仅采用基于MultiApp程序模块的传统Picard迭代耦合方法，并未实现多部件的直接联立求解。对于多回路耦合问题，蒸汽发生器是连接一、二回路的关键设备，由于没有有效解决基于有限元方法的多相流相变计算问题，RELAP7至今未开发出真实可用的蒸汽发生器模型，同样没有实现多回路耦合直接联立计算[14]。

1.2 国内研究现状

国内在反应堆耦合计算，特别是直接联立求解方面的研究工作起步较晚。2015年，张汉等[19-20]研究了直接联立求解多物理、多尺度耦合问题，实现了直接联立求解球床-燃料球两级多尺度问题，进一步提出了基于非线性消去JFNK方法的更优方案，并应用于高温气冷堆瞬态分析程序TINTE，证明了直接联立耦合方法的优势。邬颖杰等[21-22]进一步提出了多层级非线性消去JFNK方法，实现了球床-燃料球-燃料颗粒3级多尺度问题，此外，对蒸汽发生器等关键部件的直接联立求解进行了探索[23]。2017年，Hu等[24]采用JFNK方法求解沸腾两相流动问题，利用半隐式求解器作为基于物理预处理，结果表明JFNK方法具有更好的收敛性且计算效率优于Newton方法。2018年，西安交通大学吴宏春、曹良志课题组对解析构建Jacobian矩阵的NK方法进行了探索，并利用物理-热工耦合程序对OECD NEACRP控制棒弹出事故基准题进行了分析，NK方法的收敛速度明显快于Picard迭代；2020年，将JFNK方法进一步扩展到堆芯物理-热工-力学耦合问题，并利用非线性消去方法将力学模型进行了局部消去，数值计算结果表明JFNK方法的计算效率优于Picard迭代[25-27]。2020年，上海交通大学刘晓晶课题组利用具有“黑箱耦合”特点的近似块Newton方法耦合中子程序SKETCH-N和热工程序COBRA，用于铅基反应堆弹棒事故分析，数值计算结果表明在相同计算精度下，近似块Newton方法的计算效率优于传统的算符分解方法和Picard迭代方法[28-29]。2021年，西安交通大学苏光辉、秋穗正课题组利用MOOSE耦合平台，分别开发了五方程两相流分析程序和棒状燃料元件性能瞬态分析程序，并进行了程序验证[30-32]。2021年，清华大学工程物理系王侃课题组利用MOOSE耦合平台分别开发了中子计算程序和热工CFD计算程序，并分析了气冷微堆的甩负荷ATWS事故[33-34]。

针对完整核电站的高效耦合计算问题，清华大学核能与新能源技术研究院以完整高温气冷堆核电站全耦合系统为研究对象，探索了多种算法，包括常规的算符分解、Picard、JFNK、NK、非线性消去JFNK/NK方法等，或不同的算法组合，并探索不同的实施方案，包括修改扩展已有的计算程序、基于通用平台重新定义物理问题、自主开发针对特定问题的具有可扩展性的耦合平台等，现已开发了多个中间版本的耦合程序。

2 耦合计算方法

耦合计算方法是反应堆耦合系统数值模拟的核心。简便起见，以2个子物理场的耦合系统为例介绍，可假设2个子物理场分别是反应堆物理和热工水力，其中x1为中子通量，f1(x1,x2)=0为中子平衡方程；x2为温度、压力、速度等热工参数，f2(x1,x2)=0为热工控制方程组。

2.1 算符分解耦合方法

算符分解是最简单的耦合计算方法，其基本思路是将整个耦合系统分解为若干个子系统，然后依次求解各子系统[19]。这种耦合方法的优点是可通过传递耦合参数的形式，方便地利用已有子系统程序进行耦合计算。缺点是存在部分物理量信息滞后的问题，如在求解中子场方程f1(x1,x2)=0时，热工参数x2来自于上一个时间步，热工反馈信息更新不及时，可证明这将导致算符分解耦合方法存在计算精度低(一阶时间精度)、稳定性差(要求很小的时间步长)、计算效率低(时间步长小，时间步循环次数多)等问题[19]，难以推广到完整高温气冷堆核电站的耦合计算。

2.2 Picard迭代耦合方法

Picard迭代耦合方法是算符分解方法的改进，其基本思路是在原方法基础上添加一个外迭代过程，即Picard迭代，通过外迭代使所有物理场全部收敛[19]。在计算中子场方程f1(x1,x2)=0时，可获得同一时间步(不是同一个迭代步)的热工参数x2信息，因此可实现高阶精度，解决了计算精度低的问题[19]。但当Picard的迭代矩阵谱半径大于1时，Picard迭代会数值发散，无法收敛，即使是可以收敛的情况，Picard迭代也只有线性收敛速度，往往导致迭代次数多、计算效率低[35]。Picard迭代的进一步发展是Anderson加速方法，Anderson加速方法利用多个Picard迭代步的信息构建新的搜索方向，达到加速收敛的目的，但目前来看加速效果可能并不理想[2]。同时，可证明当Picard迭代不收敛时，Anderson加速方法同样不收敛，没有解决收敛性的问题[36]。

2.3 直接联立耦合方法

算符分解方法、Picard方法主要采用系统解耦、分别求解的计算思路，各物理量异步更新。以NK、JFNK为代表的耦合方法打破传统耦合计算思路，将所有子系统联立为一个大规模的非线性方程组，整体求解、同步更新，具有更高的收敛速度。

1)NK方法

NK方法是一种典型的直接联立耦合方法，具有外-内迭代的2层迭代结构。NK方法的外迭代是基于Newton方法的高效非线性求解方法，求解整个系统形成的大规模非线性方程组；内迭代是针对线性问题的Krylov子空间方法，高效求解Newton线性化形成的大规模线性方程。NK方法需要显式建立、存储Jacobian矩阵，高效计算Jacobian矩阵是关键技术。由于Jacobian矩阵元素的解析表达式依赖于物性参数等的具体表达式，解析法受限于通用性，一般较少使用。基于有限差分方法构建Jacobian矩阵是一种广泛使用的方法[37]。

2)JFNK方法

JFNK方法是在NK方法的基础上进一步的发展，同样具有Newton迭代-Krylov子空间迭代构成的2层迭代结构。不同于NK方法显式构建Jacobian矩阵，JFNK方法利用有限差分直接近似Jacobian矩阵-向量积，避免了显式构建、存储Jacobian矩阵[19，38]。

3)非线性消去JFNK/NK方法

非线性消去JFNK/NK方法是JFNK或NK方法的一种变体[19,20-22]，其思路是：(1)将耦合系统中的一部分，如子物理系统x1作为中间变量，而不是最终解向量；(2)通过求解方程f1(x1,x2)=0，对子物理系统x1进行非线性消去，得到x1=h(x2)；(3)再利用JFNK或NK方法求解原问题的等效方程f2(h(x2),x2)=0得到x2[19]。这种耦合方法的特点是：一方面，可缩小Newton迭代、Krylov子空间迭代求解问题的规模，减少迭代次数，这对计算效率是有利的；另一方面，由于有附加的非线性消去过程x1=h(x2)，增加了额外的计算代价。因此，总的计算效率是两个因素的综合，这与问题特点相关。对于高温气冷堆多尺度耦合问题，带非线性消去过程的JFNK/NK方法计算效率较一般的JFNK/NK更优。近似块Newton方法可看作非线性消去JFNK/NK方法的一种近似处理和变体[29，39]，两者区别主要体现在对f1(x1,x2)=0的消去计算上。非线性消去JFNK/NK方法要求准确求解非线性方程f1(x1,x2)=0。近似块Newton迭代只求解f1(x1,x2)=0的Newton线性化方程(即∂f1/∂x1·δx1=-f1(x1,x2))，优势是减少了消去过程的附加计算量，缺点是没有完全求解f1(x1,x2)=0，Newton迭代次数会增加，总计算效率是由两个因素综合决定的。

3 高温气冷堆联立求解的探索

3.1 高温气冷堆核电站的耦合特点与需求

球床模块式高温气冷堆核电站采用氦气作为冷却剂、石墨作为慢化剂、TRISO包覆颗粒作为球形燃料元件，具有固有安全性好、堆芯出口温度高、可提供多种工艺热的特点。以中国20万千瓦级模块式高温气冷堆示范电站(HTR-PM)为例，其采用的是两个核蒸汽供应模块的“两堆带一机”运行模式，完整核电站是多种耦合方式的组合效应，具有多层级的耦合结构，如图1所示。

图1 高温气冷堆耦合系统

在堆芯方面，由于采用球形燃料元件随机堆砌形成的球床堆芯，导致氦气流动呈现显著三维特征，各方向流动耦合更加紧密且出入口温度变化范围大，与反应堆物理构成紧密的多物理耦合；考虑进水进气事故中高温石墨与水蒸气或氧气发生化学反应，产生CO、CO2、H2并与冷却剂氦气相互混合，构成了更复杂的反应堆物理-热工水力-化学反应多物理耦合效应。球床堆芯由数十万个球形燃料元件随机堆积而成，每个燃料元件内随机弥散着1万多个TRISO包覆颗粒，而每个包覆颗粒具有燃料核芯、多个包覆层等内部结构。在反应堆物理-热工计算中，通常采用多孔介质球床模型描述宏观球床导热及与氦气的对流换热，而核截面温度反馈则需要更加精细的介观燃料球和微观燃料颗粒内部温度信息，因此需全局多孔介质球床温度-介观燃料球温度-微观燃料颗粒温度的耦合计算，这构成了高温气冷堆多尺度耦合。

在部件方面，反应堆堆芯是一个具有多物理、多尺度内部耦合特征的重要部件。此外，螺旋管式直流蒸汽发生器部件具有快速响应，紧密耦合的一、二次侧构成的内部结构。同时，由于一次通过的运行模式使二次侧给水直接被加热至过热状态，存在复杂相变和非线性更强的相间耦合。在高温气冷堆一回路，堆芯、蒸汽发生器、主氦风机和热氦导管等具有内部耦合特征的部件，通过依次流经各部件的氦气紧密的耦合，形成更高层级的多部件耦合。

在系统层面，由多个部件组成的一回路系统和二回路系统通过直流蒸汽发生器相耦合，构成多回路耦合。此外，舱室冷却系统与一回路系统之间通过压力容器表面的热辐射实现了另一种多回路耦合关系。由于“多堆带一机”的模式，不同核蒸汽供应模块系统(每个核蒸汽供应模块由1个堆芯部件和1个蒸汽发生器部件组成)产生的蒸汽在共用的蒸汽母管中汇合，通过共用的二回路系统为不同核蒸汽供应模块供应给水，这样虽然每个核蒸汽供应系统模块的一回路是独立的，但给水温度、蒸汽压力通过共用二回路相耦合，构成超大规模、更复杂的多模块耦合。

高温气冷堆是一个具有多物理、多尺度、多部件、多回路、多模块特征的复杂、非线性、强耦合系统，呈现出堆芯级-部件级-系统级的多层级嵌套耦合结构，如图2所示。传统耦合方法可实现堆芯级，甚至部分实现部件级问题的耦合计算，但受限于数学本质缺陷，很难进一步扩展到更高层次、更复杂的系统级问题，更为先进的直接联立耦合方法是具有潜力的解决方案[40-41]。直接联立耦合方法的基本思路是，根据各层级问题的自身特点，在非线性求解、线性求解、预处理等方面进行分层级处理，在整体上利用Newton方法进行全局求解，实现稳定、高效、准确求解。清华大学核能与新能源技术研究院以高温气冷堆全耦合系统为研究对象，以直接联立方法为框架，对关键算法和耦合平台开展了长期研究。

图2 高温气冷堆多层级耦合结构

3.2 完整核电站的联立耦合算法探索

直接联立耦合方法将整个耦合系统作为一个整体进行求解，需要求解一个大规模的非线性方程组，还需考虑高温气冷堆核电站具有多种耦合方式、多层级耦合结构的特点。因此，需要开发更高效的联立耦合算法。

1)预处理方法

直接联立耦合方法是采用Krylov子空间迭代方法求解Newton线性化后得到的大规模线性方程组。预处理方法是保证Krylov子空间迭代计算效率的关键，是反应堆耦合计算的共有问题，包括预处理方式的选择及预处理矩阵的构建两个核心问题。对于预处理方式，主要包括线性预处理和非线性预处理[42-45]。线性预处理是对线性方程组的预处理，左预处理如式(1)所示，其中P是预处理矩阵，详细讨论参见文献[42-45]。而非线性预处理则是通过对非线性方程F(x)=0本身作等价变换，使Jacobian矩阵具有良好性质，起到加速效果，如式(2)所示。由于非线性预处理可写成不动点迭代格式，因此具有黑箱耦合能力，易通过改造现有代码实现。对于预处理矩阵P，在耦合问题中，一般采用基于物理预处理方法，仅保留强耦合项，如最简单的块Jacobian预处理。代表子物理场本身的对角子矩阵块一般也是近似求解的，如采用不完全LU分解(ILU(k))等方法。

P-1·J·δx=-P-1·F(x)

(1)

0⟺x-G(x)=0

(2)

文献[44-47]针对堆芯物理-热工耦合问题，研究了不同预处理矩阵下线性/非线性预处理的计算效率[44](表1)。线性预处理的计算效率一致高于相同预处理矩阵下的非线性预处理方法，这是因为非线性预处理需在每个Krylov子空间步计算一次矩阵求逆，而线性预处理只需在每个Newton步进行矩阵求逆[44]。随填充水平k的增加，Krylov子空间迭代次数减少，而每个迭代步的计算代价增多。非线性预处理对矩阵求逆代价更敏感，不适用代价高而迭代次数少的预处理方式。

表1 线性/非线性预处理的计算性能[44]

2)尺度变换技术

不同物理场的物理量数值差异非常大，如热群中子通量约为1018m-2·s-1，而燃料代表温度约为103℃。对于所有物理场整体求解的直接联立方法，这就导致表征整体耦合关系的Jacobian矩阵各元素之间的数值差异非常大，其条件数可达到约1032，极其病态，是反应堆耦合计算的共性问题。通过左、右尺度变换技术可有效地改善Jacobian矩阵条件数(式(3))[48]。右尺度变换SR的含义是使所有物理量的数值归一化到相同量级；左尺度变换SL的含义是使所有方程的残差归一化到相同量级，问题收敛时各方程达到相同的收敛水平。尺度变换因子的选取是关键，一种经验性尺度变换因子选取方式如式(4)所示[48]，右尺度变换因子采用各物理场的典型数量级的倒数，左尺度变换因子可采用与右尺度变换相同的因子。但这种方法不能保证得到的是最优尺度变换因子。另一种确定因子的方法是将其转换为最优化问题，优化变量是左、右变换尺度因子，优化目标是Jacobian矩阵条件数最小，如式(5)所示。为了提高计算效率，一般采用Ruiz算法等启发式方法求解这个最优化问题，得到最优尺度变换因子[49]。针对堆芯物理-热工耦合问题，经验性尺度变换因子可使Jacobian条件数降低到约1010，最优尺度变换因子可进一步下降到约103，同时Krylov子空间的迭代次数也会进一步减少[49]。

(3)

(4)

(5)

3)多尺度耦合方法

高温气冷堆堆芯同时具有多物理、多尺度2种耦合方式，但对其他堆型有参考。高温气冷堆需要求解宏观球床温度-不同批次的介观燃料球温度-微观颗粒温度，具有2个不同于简单多物理耦合问题的特点：(1)介观和微观变量个数远多于宏观变量个数，若考虑15个不同批次的燃料球，每个燃料球需15个网格，而每个燃料球网格内包含1个需要5个网格描述的颗粒模型，则微观/介观待求量约是宏观待求量的170倍；(2)介观和微观变量相对于堆芯问题更易求解，这是因为燃料球和颗粒都是一维模型，同时燃料球和颗粒的非线性仅是由于导热系数是温度的函数导致的，非线性较弱，这2个特点使得问题更易求解。基于以上2个特点，提出了非线性消去方法消去介观燃料温度和微观颗粒温度，这可大幅减少问题规模，进而减少Krylov迭代次数；同时由于介观和微观方程都是弱非线性的一维问题，消去过程的额外计算代价很小，因此总的计算效率可显著提高。针对宏观球床温度-介观燃料球温度二级多尺度问题，计算效率相比于一般的JFNK直接联立方法可提高3～4倍[20]。对于更复杂的宏观-介观-微观三级多尺度问题，介观-微观的非线性消去过程本身仍是多尺度问题，其系数矩阵如图3所示，线性迭代次数可减少90%甚至更多，计算效率可进一步提高到4～7倍[22]。

图3 介观-微观系数矩阵的非零元素结构[22]

4)蒸汽发生器部件

高温气冷堆采用螺旋管式直流蒸汽发生器，其一次侧氦气的热惯性很小、二次侧水装量少，比自然循环蒸汽发生器响应更快、耦合更紧密，二次侧流动换热现象更复杂，入口过冷水被直接加热到过热蒸汽，存在单相液到两相流、再到单相过热蒸汽的2次复杂相变的强非线性过程。构建高效预处理矩阵是蒸汽发生器联立求解的关键，不同于基于共同计算域的多物理耦合，蒸汽发生器是通过共用的换热界面相耦合，根据换热过程可构建高效的基于物理预处理矩阵[23,50]。在直接联立耦合方法中，相产生(空泡份额α=0)/相消失(α=1)都会导致缺失相所对应的Jacobian子矩阵块奇异，缺失相的方程具有无穷多解。动态残差法是一种有效的处理方法[23]，使耦合系统具有符合物理意义的唯一解。表2为高温气冷堆HTR-10蒸汽发生器的设计值和直接联立计算值。图4为蒸汽发生器不同耦合算法对比，可看出，预处理JFNK方法效率显著高于Picard迭代方法。

表2 高温气冷堆HTR-10蒸汽发生器模拟结果[23]

图4 蒸汽发生器不同耦合算法对比[23]

5)其他关键技术

(1)Jacobian矩阵的快速计算：高效计算Jacobian矩阵是NK方法中的关键技术。基于有限差分方法构建Jacobian矩阵是一种有效方法，利用矩阵稀疏性，采用图着色方法可有效减少有限差分计算次数，大幅提高计算效率[51-53]。图着色NK方法的效率有潜力优于JFNK方法[52-53]。

(2)统一的高精度离散格式：耦合系统中各物理场都是各自发展，有必要研究统一的高阶离散格式。清华大学核能与新能源技术研究院推广了节块法，开发了针对对流方程和瞬态扩散方程的高精度节块法离散格式，并实现了基于统一高精度节块法的反应堆物理-热工联立耦合计算[54-59]。

(3)缓发中子方程、中子毒物方程的处理：缓发中子和中子毒物由线性常微分方程描述，因而反应堆物理-热工水力耦合系统是偏微分-常微分混合方程组。可将消去方法推广到线性方程，将缓发中子先驱核和中子毒物浓度作为中间变量，能有效减小问题规模，提高计算效率[21，60]。

3.3 完整核电站的耦合软件框架探索

1)基于已有程序的改造

TINTE程序是高温气冷堆安全分析程序，具有精细的物理模型，但只采用传统的算符分解方法。清华大学核能与新能源技术研究院基于TINTE程序的原各物理场的计算模块，实现了JFNK直接联立求解物理、传热、流动这种多物理耦合，同时考虑了球床堆芯宏观分布与燃料球内的介观尺度分布的2级多尺度耦合问题。为了对原计算程序改动尽量小，采用了非线性预处理方式。为了解决空间多尺度耦合问题，发展了带有非线性消去过程的JFNK方法。研究结果表明，JFNK联立耦合方法具有更高的精度、更强的稳定性，可采用更大的时间步长达到更高的计算效率[48,61]，如图5所示。

a——热中子通量；b——计算精度；c——计算稳定性

针对多回路、多模块耦合方式，为了简化问题、突出重点，对堆芯模型和二回路部件模型进行了简化，同时采用已有的蒸汽发生器程序BLAST，开发了轻量级的耦合计算程序。同样地，为了减少对原计算程序的修改，采用非线性预处理方式。作为探索性尝试，采用Picard和JFNK的混合求解方案，初步说明了用JFNK实现多回路、多模块耦合系统的技术可行性[62-63]。经过以上探索，基于已有程序的扩展和改造可实现JFNK直接联立耦合计算，并验证了JFNK方法在计算稳定性、数值精度、计算效率等方面的优势。但受限于原有程序计算流程，只能采用具有黑箱耦合特性的非线性预处理，而不是效率更高的线性预处理，不能充分发挥直接联立方法的优势。同时，原有程序框架对更大范围、更复杂耦合问题的扩展也存在诸多限制，难以扩展到完整高温气冷堆核电站耦合系统的直接联立求解，需要开发更具可扩展性的耦合平台型程序。

2)基于通用耦合程序的开发

相比于修改原有代码，基于以MOOSE为代表开源耦合平台进行多物理耦合计算具有标准化、模块化、可扩展、可并行的特点。清华大学核能与新能源技术研究院基于MOOSE通用耦合框架开发了具有高温气冷堆堆芯多物理特征的中子程序模块和热工水力程序模块[64]，如图6所示。新开发的多物理耦合程序采用了更高效的线性预处理方法，并对预处理矩阵的构建方式、预处理矩阵的分解方式、稳态问题初值设置等问题进行了系统研究和探索[64]。基于MOOSE平台开发了介观燃料球-微观燃料颗粒的2级多尺度耦合程序，但没有实现完整的3级多尺度直接联立计算[64]。

图6 基于MOOSE的物理-热工计算模块[64]

MOOSE采用基于面向对象描述方式的耦合平台技术路线，用户只需要提供各个物理场的偏微分方程算符表达式，就可方便快速地构造出一个复杂耦合系统，极大地缩短了求解耦合问题的程序开发周期。但在框架设计阶段，MOOSE的定位是针对不同领域的通用耦合平台，对反应堆耦合问题的特殊性考虑不足。现阶段高温气冷堆的某些复杂精细模型，如氦气滞止区、控制棒、球床-反射层换热等，在MOOSE中实现难度很大。此外，MOOSE耦合平台的本质是对同一区域、单层级、单一耦合方式多物理耦合问题的直接联立计算。但高温气冷堆核电站全耦合模型涉及多回路、多模块耦合，同时叠加多部件、多尺度、多物理耦合，MOOSE不具有模拟能力。

3)自主耦合平台的开发

经过已有程序改造、通用平台开发等尝试，清华大学核能与新能源技术研究院最终确定了针对高温气冷堆核电站的特性，自主开发可扩展耦合计算平台的技术路线。自主耦合计算平台的开发是在前期工作的基础上开展的，原有计算程序的深度理解和解析提供了核心计算模块的物理计算模型，通用耦合平台MOOSE的工作对自主平台设计提供了重要参考。自主耦合计算平台针对反应堆耦合特性而设计，从直接联立耦合方法的机制出发，将整个耦合问题自上而下逐层次分解，在不同层次的耦合机制中凝练出基本的、共同的问题描述，并进行抽象化、规范化，然后就可形成一个针对高温气冷堆核电站多层级耦合结构的、具有扩展性的平台化耦合计算程序框架，如图7所示。

图7 高温气冷堆耦合计算平台

现阶段已初步实现针对高温气冷堆一回路多物理、多尺度、多部件耦合系统的计算框架和核心计算模块的开发，自上而下包括物理模型层、离散格式层及更底层的数值计算层等。物理模型层现已开发2个层次，第1层是部件层次，主要包括反应堆堆芯、蒸汽发生器、主氦风机和相应的连接管道。下一个层次是物理场层次，主要包括堆芯部件中的堆芯中子场、堆芯热工场、介观燃料球、微观燃料颗粒和中子毒物等。数值离散层的功能是将物理场控制方程离散为代数方程组，由1个基础类和1个抽象类完成，主要包括时间项离散、空间项离散、不同方程间耦合项、边界条件离散。数值计算层的功能是离散方程的数值求解，包括非线性方程组求解、线性方程组求解、高效线性预处理计算、稀疏矩阵运算等。在此基础上，未来将进一步完善一回路多部件耦合模型，开发多回路、多模块等更高层次的物理模型层的计算模块，实现完整高温气冷堆核电站全耦合系统的直接联立计算。在开发新的物理模块时，由于方程的离散和求解可由已有的数值离散层和数值计算层完成，因此只需关注物理层的程序开发，可大幅缩短程序开发周期，这是平台化程序的优点。

通过上述定义的模块及基础类，用户可快速地组装出需要模拟的球床式高温气冷堆耦合问题，完成中子物理、热工流体、多尺度燃料及蒸汽发生器等特殊部件的直接联立耦合计算，目前可实现的功能包括堆芯内部的瞬态/稳态计算、球床-燃料球-颗粒空间温度分布的计算、堆芯功率调节瞬态工况等，主要的模拟结果如图8所示。未来将利用HTR-PM实测数据进行程序验证。

a——堆芯固体温度；b——介观-微观燃料元件温度；c——蒸汽发生器温度

4 结论

本文在综述国内外反应堆耦合计算研究的基础上，以高温气冷堆核电站为研究对象，介绍了清华大学核能与新能源技术研究院在直接联立求解方法及程序开发方面的研究工作。介绍了非线性消去直接联立方法、高效预处理方法、多相流联立求解等关键技术，讨论了修改扩展已有的计算程序、基于通用平台重新定义物理问题、自主开发针对特定问题的具有可扩展性的耦合平台等不同的实施方案和技术路线。结果表明，直接联立计算方法潜力很大，但还有很多关键技术需要研究、攻关，才能最终解决此问题。

感谢王登营、邓志红、周夏峰、卢佳楠和牛进林等为高温气冷堆耦合计算做出的重要贡献。