反应堆系统瞬态仿真方法和装置、系统、存储介质

本发明属于反应堆系统仿真，具体涉及一种反应堆系统瞬态仿真方法和装置、系统、存储介质。

背景技术：

1、反应堆系统仿真主要涉及到冷却剂流动、堆芯释热、燃料传热、压力边界、换热构件及安全系统动作等多个过程。目前的系统分析程序大多采用控制体、连接件、热构件和边界条件对反应堆一回路、二回路及相关设备进行建模，然后通过质量守恒、动量守恒和能量守恒方程进一步描述系统瞬态响应。这种建模方式已广泛用于工程分析中，但在局部沸腾、相出现、相消失、流型转换和强热反馈等工况下，方程组离散后会表现出明显的非线性、强耦合和病态特征。如果继续采用显式推进、半隐式处理或分步迭代方式等方法，会导致计算结果受到时间步长、初始场分布和变量更新顺序的影响，并且出现收敛速度下降或迭代稳定性不足等问题。

2、以反应堆热工水力两相流计算和中子动力学耦合计算为例，计算过程中存在较多的强耦合变量。在热工水力两相流模型中，压力、空泡率、液相速度、气相速度、液相焓、气相焓、密度、相变质量源项、壁面换热、相间换热和界面阻力等变量相互影响以及堆芯中子动力学计算中，中子通量、缓发中子先驱核密度和裂变释热功率又与燃料温度、冷却剂温度、冷却剂密度及气相空泡率存在相互作用的关系。而且材料温度和冷却剂状态的变化也会引起中子扩散系数及宏观反应截面的变化，导致裂变释热功率进一步作为燃料传热方程和热工水力能量方程的源项。上述变量经过空间离散和时间离散后，与闭式关系、源项模型、物性模型及边界条件共同组成大规模非线性离散方程组。目前，在处理上述非线性离散方程组较常见的方法包括 picard 迭代、算符分裂、变量顺序更新、松耦合交替迭代和局部线性化等。这些方法实现较为简单，但其收敛过程容易受到物理场更新顺序、耦合容差、时间步长以及闭式关系切换的影响。在气泡产生、相消失、流型转换、界面阻力突变、换热关系切换、快速功率变化或强反馈瞬态条件下，残差函数可能出现明显非光滑变化，进而导致迭代残差停滞、振荡、耦合场不一致或非物理解等问题。对于需要保持较大时间步长、强耦合一致性和瞬态计算稳定性的反应堆系统仿真，上述问题会限制计算效率和结果可靠性。然而全隐式求解方法能够在同一非线性方程组中同时更新守恒变量、闭式关系相关变量以及多物理场耦合项，具有较好的稳定性，适合处理强耦合、多尺度和强非线性的反应堆系统瞬态仿真问题。因此，在全隐式框架下，采用显式 jacobian 矩阵构造方法求解同一时间步内的大规模非线性离散方程组，会带来较高的存储开销和计算开销，反之，采用无 jacobian 的newton-krylov 方法，则内层 krylov 迭代对预处理矩阵质量较为敏感。对于两相流、传热传质和中子动力学耦合模型，还需要同时考虑方程稀疏结构、变量物理尺度差异、闭式关系非光滑性以及相出现或相消失导致的病态特征对预处理效果的影响。因此，需要一种适用于反应堆系统仿真的全隐式非线性离散方程求解方法，在保持全隐式耦合求解稳定性的同时，降低预处理构造成本，并提高复杂工况下非线性迭代和线性内迭代的收敛可靠性。

技术实现思路

1、为解决现有技术存在的问题，本发明提供一种反应堆系统瞬态仿真方法和装置、系统、存储介质，在保持全隐式耦合求解稳定性的同时，降低预处理构造成本，并提高复杂工况下非线性迭代和线性内迭代的收敛可靠性。

2、为实现上述目的，本发明提供了如下方案：

3、一种反应堆系统瞬态仿真方法，包括：

4、步骤s1、确定反应堆仿真中的待求解对象，并接入与所述待求解对象对应的待求解反应堆物理模型；通过所述待求解反应堆物理模型的模型评价接口，将其状态变量与非线性残差响应表示为统一残差方程组 f(u)=0，其中 u 为当前仿真时间层待求解的状态变量向量。

5、步骤s2、可根据当前仿真时间层的初始条件与边界条件、上一仿真时间层的收敛状态以及稳态计算结果进行初始化统一残差方程组的状态变量和待求解反应堆物理模型的模型参数，并通过所述模型评价接口获得初始非线性残差、模型分支状态和物理约束状态。

6、步骤s3、分析所述统一残差方程组中各未知量与各残差分量之间的稀疏依赖关系，并根据所述稀疏依赖关系建立变量冲突图；其中，所述变量冲突图用于表示两个未知量是否共同影响至少一个残差分量，若两个未知量共同影响至少一个残差分量，则两个未知量不得划分至同一颜色组中同时扰动；

7、步骤s4、对所述变量冲突图进行图着色处理，得到多个颜色组，其中，同一颜色组内的未知量互不冲突，不同颜色组之间的未知量按照颜色顺序分别进行有限差分扰动；

8、步骤s5、利用图着色预处理的jacobian-free newton-krylov方法求解所述非线性残差方程组，具体包括：

9、步骤s51、基于newton-krylov方法进行外循环迭代计算，通过当前的迭代计算将所述非线性残差方程组进行线性化处理，得到用于求解newton修正方向的线性方程组。其中：所述线性化过程是通过状态变量扰动和残差响应所获得的jv乘积，即：jacobian 矩阵j 乘以一个向量v；

10、步骤s52、基于所述颜色组分别生成有限差分扰动向量，对同一颜色组中的未知量同时进行扰动后重新计算残差得到扰动前后的残差变化，根据未知量与残差分量之间的对应关系，将所述残差响应重新组成图着色预处理矩阵；

11、步骤s53、求解所述线性方程组的内循环迭代则用krylov子空间迭代法，并采用所述图着色预处理矩阵对内循环进行预处理，从而获得newton修正方向；

12、步骤s6、在newton迭代过程中对强非线性工况进行一种或多种的稳定化处理，所述稳定化处理包括：闭式关系延迟更新、物理约束、主动集判定、消失相速度限制和线搜索步长控制；

13、步骤s7、据所述newton修正方向和线搜索步长形成修正状态，对所述候选状态施加物理可行域约束、相状态约束和线搜索接受判断；修正状态被接受后，更新待求解反应堆物理模型的状态变量和模型参数；

14、步骤s8、当更新后的非线性残差满足收敛条件时，获得当前仿真时间层的反应堆物理场参数；根据所述反应堆物理场参数确定反应堆系统的冷却剂流动状态、传热状态或者中子功率状态，并将所述反应堆物理场参数作为下一仿真时间层的初始状态；输出反应堆系统在各仿真时间层的瞬态物理场分布，直至达到预设仿真终止条件。

15、作为优选，选取的待求解对象为反应堆热工水力系统中的两相流模型，则采用全隐式两流体方程离散后所得到非线性离散的反应堆物理模型，其中：所述状态变量向量u包括压力、气相空泡率、液相比焓、气相比焓、液相质量流量和气相质量流量；所述非线性残差方程组包括气相质量方程、液相质量方程、气相动量方程、液相动量方程、气相能量方程和液相能量方程。

16、所述反应堆热工水力系统中的两相流模型如下式所示：

17、

18、其中：为相指标，g表示汽相，l表示液相；r为空间方向指标。上标n和n+1分别表示上一时间层和当前隐式时间层；表示时间步长，单位s，下标p表示控制体中心节点。表示控制体体积，单位为；f表示控制体边界面，表示控制体p的所有边界面集合；表示控制面面积,单位为。表示第s相密度,单位为；表示第s相体积分数。表示第s相在r方向上的速度分量,单位为m/s；表示第s相在控制面f上的法向速度,单位为m/s。p表示压力,单位为；表示控制体中心处沿r 方向的压力梯度,单位为pa/m；表示重力加速度在r方向上的分量,单位为。表示第s相比焓，单位为j/kg；表示第s相的相变质量源项，单位为。表示第s相在r方向上的动量源项，单位为。表示第s相能量源项，单位为。式中所带n+1上标的变量均在当前时间层求解，因此该方程组构成全隐式非线性离散方程组。

19、如果其待求解对象为三维中子扩散方程组，其特征在于，状态变量向量包括各计算网格单元的中子通量、缓发中子先驱核密度和裂变释热功率，所述非线性残差方程组至少包括中子通量平衡残差、缓发中子先驱核守恒残差和裂变释热关系残差。其中，中子通量平衡残差依赖当前网格单元及相邻网格单元的中子通量、缓发中子先驱核密度和中子反应参数，裂变释热功率根据中子通量、宏观裂变截面和计算网格体积确定。

20、所述的三维中子扩散方程组如下式所示：

21、

22、其中，分别表示网格单元处的中子通量方程残差、缓发中子先驱核方程残差和裂变释热方程残差；上标表示当前时刻的变量值，上标表示上一时刻的变量值；下标分别表示三维计算网格在方向上的网格编号；表示时间步长，单位为。表示中子通量，单位为表示网格单元的中子净流密度，单位为表示缓发中子先驱核密度，单位为表示中子运动速度，单位为表示缓发中子先驱核衰变常数，单位为表示与中子吸收和裂变反应有关的离散系数；表示与三维空间中子泄漏有关的离散系数；缓发中子先驱核衰变源表示先驱核产生项；表示外加中子源以及其他已知项形成的离散源项。上述离散系数根据中子反应截面、计算网格尺寸、材料温度和材料密度确定。表示网格单元的裂变释热功率，单位为表示单次裂变释放的有效能量，单位为表示宏观裂变截面，单位为表示网格单元体积，单位为。

23、作为优选，所述变量冲突图的建立和图着色预处理矩阵的构造具体包括：假如两个未知量共同影响至少一个残差分量，则在所述两个未知量对应的节点之间建立冲突边。反之，两个未知量不存在共同影响的残差分量，则允许所述两个未知量划分至同一颜色组。针对每一颜色组，生成对应的有限差分扰动向量，该扰动向量仅在这一颜色组对应的未知量位置取非零值，其余位置取零。然后将所述扰动向量叠加到当前状态变量后，重新计算非线性残差，得到该颜色组对应的残差差分。再根据变量冲突图和残差依赖关系，并确定残差差分与各未知量之间的对应关系后，将相应差分结果写入预处理矩阵的非零位置，形成稀疏图着色预处理矩阵；其中，所述有限差分扰动量根据未知量的物理尺度确定，并通过有界缩放方式限制扰动幅值。

24、作为优选，所述强非线性工况包括热工水力闭式关系、中子场反馈参数以及相出现或相消失工况；所述稳定化处理包括：在图着色预处理矩阵构造所进行的分组有限差分残差评价时，保持当前newton状态下的闭式关系参数和主动集状态不变。当 newton 修正步满足线搜索接受条件后，根据更新后的状态变量重新计算所述热工水力闭式关系参数和中子场反馈参数，newton 线性方程右端项、线搜索接受准则和收敛判断采用修正状态重新评价后的原始残差。对修正状态变量施加物理可行域约束，使空泡率、压力、焓、相速度、中子通量、缓发中子先驱核密度和裂变释热功率保持在预设物理范围内。

25、本发明还提供一种反应堆系统瞬态仿真装置，包括：

26、第一处理模块，用于确定反应堆仿真中的待求解对象，并接入与所述待求解对象对应的待求解反应堆物理模型；通过所述待求解反应堆物理模型的模型评价接口，将其状态变量与非线性残差响应表示为统一残差方程组；

27、第二处理模块，用于分析所述统一残差方程组中各未知量与各残差分量之间的稀疏依赖关系，并根据稀疏依赖关系建立变量冲突图；

28、第三处理模块，用于对所述变量冲突图进行图着色处理，得到多个颜色组，其中，同一颜色组内的未知量互不冲突，不同颜色组之间的未知量按照颜色顺序分别进行有限差分扰动；

29、第四处理模块，用于利用图着色预处理的jacobian-free newton-krylov方法求解所述非线性残差方程组；

30、第五处理模块，用于在newton迭代过程中对强非线性工况进行一种或多种的稳定化处理，所述稳定化处理包括：闭式关系延迟更新、物理约束、主动集判定、消失相速度限制和线搜索步长控制；

31、第六处理模块，用于根据所述newton修正方向和线搜索步长形成修正状态，对所述候选状态施加物理可行域约束、相状态约束和线搜索接受判断；修正状态被接受后，更新待求解反应堆物理模型的状态变量和模型参数；

32、第七处理模块，用于当更新后的非线性残差满足收敛条件时，获得当前仿真时间层的反应堆物理场参数；根据所述反应堆物理场参数确定反应堆系统的冷却剂流动状态、传热状态或者中子功率状态，并将所述反应堆物理场参数作为下一仿真时间层的初始状态；输出反应堆系统在各仿真时间层的瞬态物理场分布，直至达到预设仿真终止条件。

33、本发明还提供一种反应堆系统瞬态仿真系统，包括：存储器和处理器，所述存储器上存储有由所述处理器运行的计算机程序，所述计算机程序在被所述处理器运行时执行反应堆系统瞬态仿真方法。

34、本发明还提供一种存储介质，所述存储介质上存储有计算机程序，所述计算机程序在运行时执行反应堆系统瞬态仿真方法。

35、与现有技术相比，本发明的有益效果为：

36、（1）本发明将待求解反应堆物理模型中的守恒方程、闭式关系、源项、边界条件及多物理场耦合关系统一构建全隐式非线性残差方程组，并采用jacobian-free newton-krylov方法在同一时间步内联立求解对应的状态变量。相较于显式推进、算符分裂或变量顺序更新方法，可减少由分步处理、物理场更新顺序及耦合容差额外引入的耦合误差，提高强耦合瞬态计算中各状态变量及物理场之间的一致性。

37、（2）本发明根据残差分量与状态变量之间的稀疏依赖关系建立变量冲突图，通过图着色将不存在共同残差依赖的状态变量划分至同一颜色组，并按颜色组执行有限差分残差评价，以构造与所述稀疏依赖关系相对应的预处理矩阵。在未知变量总数为、颜色组数量为，且的情况下，每次预处理矩阵构造所需的残差评价次数可由逐变量扰动时约次减少至约次，从而降低预处理矩阵构造的计算开销。结合jacobian-free矩阵—向量乘积，可避免显式组装完整jacobian矩阵，并有利于改善krylov内迭代的预处理效果和收敛性能。

38、（3）针对相出现、相消失、近单相、流型转换及闭式关系非连续切换等强非线性工况，本发明在一次图着色预处理矩阵构造过程中保持当前newton状态下的闭式关系参数和主动集状态不变，在newton修正步满足接受条件后再根据接受状态更新相关参数。该处理可降低有限差分扰动触发闭式关系非连续切换对预处理矩阵近似精度的影响。进一步结合物理可行域约束、空泡率主动集约束、消失相速度限制和线搜索步长控制，可降低残差振荡、相速度发散及非物理解出现的风险，提高复杂工况下非线性迭代的稳定性和求解结果的物理可接受性。

39、（4）本发明以非线性残差方程、残差与状态变量之间的稀疏依赖关系、变量物理尺度、模型参数更新规则及物理约束作为模型接入信息。在接入不同待求解反应堆物理模型，并配置其残差评价、变量依赖关系、变量尺度、模型参数更新规则及物理约束信息后，图着色预处理和jacobian-free newton-krylov求解流程能够保持不变。因此，该求解框架可用于反应堆热工水力、两相流、传热传质、中子动力学及其多物理场耦合模型的全隐式非线性求解，并支持不同反应堆物理模型在统一求解框架下的应用。