CN106202866B - 一种稳定精确的反应堆物理热工耦合计算方法 - Google Patents

Abstract

一种稳定精确的反应堆物理热工耦合计算方法，是一种在求解过程中同时得到中子学场，流体场和传热场状态的方法，由于所有变量在迭代过程中同时更新，因此本方法较传统耦合方法具有更佳的数值稳定性，本方法通过联立时空中子扩散方程和热工方程的离散方程，得到一个新的更大规模的非线性方程组，对于这个方程组，应用牛顿迭代法进行数值求解，由于牛顿迭代法需要在每一个迭代步计算方程组的雅阁比矩阵，对于物理热工耦合计算较为困难，本方法通过建立扩散方程和热工方程间的数值联系，显式计算出雅阁比矩阵，避免了由于近似计算雅阁矩阵导致的问题，同时，本方法继承了牛顿法的特性，有望具备更快的收敛速度。

Description

一种稳定精确的反应堆物理热工耦合计算方法

技术领域

本发明涉及反应堆堆芯设计和安全分析技术领域，具体涉及一种稳定精确的反应堆物理热工耦合计算方法，称之为全耦合法。

背景技术

反应堆是一个多物理场耦合的复杂系统，由于其复杂性，针对各个物理场形成了相对独立的学科，其中核反应堆物理分析以及核反应堆热工分析是反应堆设计和安全分析的两大重要领域。

所谓物理热工耦合，是指堆内的物理参量间互相存在反馈效应，这些反馈效应在反应堆发生瞬态安全性事故时会更加显著，因此准确地考虑这些反馈效应则需要物理热工的耦合分析计算。

目前广泛采用的物理热工耦合计算方法是一种简单耦合，国际上亦称之为“代码耦合”或“松耦合”，物理方程和热工方程的求解彼此是解耦进行的，其实现方式是通过外部脚本或代码合并，将一个专门求解反应堆中子学问题的程序和一个专门求解反应堆热工学问题的程序联接，物理程序将计算得到的功率分布传递给热工程序，热工程序将计算得到的热工参量反馈给物理程序，这一过程中仅存在数据交互，而核心求解过程彼此独立，互为“黑匣子”，在此基础上，在瞬态分析计算中又分为显示耦合方法和隐式耦合方法，其区别在于显示耦合在一个时间步内只进行一次物理计算和一次热工计算，而隐式耦合需要进行物理热工间的迭代，直到计算结果在一个时间步内满足一定的收敛判据。

现有的物理热工耦合方法，最大的优势在于易于实现，因为它几乎不需要程序再开发，然而，由于这种“松耦合”的算法本质，其收敛速度有限，若要获得高收敛精度的结果需要进行很大的物理热工迭代次数，因此需要在计算精度和计算效率间平衡，此外现有耦合方法在某些情况下亦会出现计算震荡甚至发散的问题，虽然通过增加物理热工耦合计算间的阻尼因子可以显著缓解该现象，但依旧不能跳出“代码耦合”的思想框架，物理方程和热工方程的依旧是解耦求解的。

毫无疑问，物理方程和热工方程的解耦是导致现有耦合计算方法存在问题的重要源头，只有将物理热工方程联立，一次计算中同时求解中子学方程和热工方程才能从根本上解决上述问题。

发明内容

为了解决上述现有技术存在的问题，本发明的目的在于提供一种稳定精确的反应堆物理热工耦合计算方法，全耦合法，该方法利用牛顿迭代法，下文简称为牛顿法，求解联立的物理热工耦合的非线性方程组，实现了物理热工的耦合计算分析，是一种与传统耦合方法迥异的新的耦合计算方法，理论上具有良好的数值稳定性和数值求解精度，同时牛顿法的超线性收敛速度亦能保证计算效率。

牛顿法求解非线性方程组的一般过程如下：

设F(x)＝0为一非线性方程组，求F(x)＝0的解。

给定初值x＝x⁰

作do循环

J_F(x^k)δx^k＝-F(x^k)

x^k+1＝x^k+δx^k

当满足收敛判据时循环结束

其中，

的雅阁比矩阵，

k--牛顿迭代步数。

为了实现全耦合方法，本发明采用如下技术方案：

一种稳定精确的反应堆物理热工耦合计算方法，包括如下步骤：

步骤1：确定中子学方程，即时空中子扩散方程的离散形式，对时空中子扩散方程采用非线性迭代半解析节块法进行空间上的离散，时间项采用隐式差分，能量取两群近似，得到带节块耦合修正因子项的节块平衡方程，即带修正因子的CMFD方程，其表达式如下：

式中，

G——能群总数；

Δt——时间步长；

——第g群n时刻节块m内的平均中子通量；

——第g群n+1时刻节块m内的平均中子通量；

v_g——第g群中子速度；

χ_pg——瞬发中子裂变谱；

β_p——瞬发中子份额；

χ_dg——缓发中子裂变谱；

ν∑_f——中子产生截面；

∑_g'g——g’群到g群宏观散射截面；

∑_t——宏观总截面；

λ_l——瞬发中子份额；

C_l——第l群缓发中子先驱核浓度；

——由先驱核浓度线性近似推导得到的常数项系数；

——第g群节块m内的中子泄漏项；

——第g群节块m内u方向左右表面的中子静流；

Δu——u方向的节块宽度；

——节块等效扩散系数；

——非线性迭代节块耦合修正因子；

——m节块在u方向左侧面的不连续因子；

确定热工方程的离散形式，对热工方程采用单通道模型，只考虑径向导热，空间上采用有限体积法离散，时间项采用隐式差分，得到有限差分格式的热工离散方程：

式中

A——通道流动面积；

ρ_j——第j网格冷却剂密度；

m_j——第j网格质量流量；

h_j——第j网格冷却剂比焓；

H——传热系数；

P_r——释热周长；

T_w——壁面温度；

T_b——冷却剂温度；

ρ_f——燃料密度；

c_p——燃料比热；

V——传热网格体积；

K——等效导热率；

q'''——体积释热率；

步骤2：由于中子学方程与热工方程在数学表达式中并无显示的联系，因此需要对求解的节块平衡方程与热工离散方程建立数学关系，在节块平衡方程中截面是慢化剂密度和燃料温度的隐式函数∑_x＝g(ρ,T_f)，热工离散方程中体积释热项是关于中子通量和裂变截面的函数从而根据这两个关系，便能构建封闭的中子学方程和热工方程耦合的非线性方程组，形式如下：

式中，

f_φ(x),f_T(x),f_h(x),f_m(x)分别为残差形式的中子方程，传热方程，能量方程和质量方程；

F(x)便是联立的物理热工方程组，x为解向量，分别为中子通量，冷却剂比焓，冷却剂质量流量，燃料温度；

步骤3：建立了物理热工耦合的非线性方程组F(x)，在用牛顿法求解F(x)时，需在每一个牛顿迭代步构造F(x)的雅阁比矩阵，如下所示，

在雅阁比矩阵元素的计算中，需要求解中子残差方程对比焓的偏导数这实际上是要计算截面关于比焓的偏导数由于截面是关于慢化剂密度和燃料温度的隐式函数∑_x＝g(ρ,T_f)，因此使用求导的链式法则计算得到，

对于雅阁比矩阵中剩下的元素的计算，均能够通过各个残差方程得到解析表达式；

通过以上运算便得到了雅阁比矩阵各个元素的数值，为了获得每个牛顿迭代步的更新量δx，需要求解以下代数方程组，

J_Fδx＝-F(x)

J_F是一个大型的稀疏矩阵，可以使用许多成熟的代数方程组迭代方法，如GMRES算法；

以k作为牛顿迭代步记号，当||F(x^k)||≤ε_Newton时认为牛顿迭代收敛， ||F(x^k)||为F(x^k)的范数，ε_Newton为牛顿迭代收敛判据；

步骤4：物理热工耦合的非线性方程组F(x)求解完毕后，需要根据最新解得的节块中子通量更新节块耦合修正因子D^k,nod，因此是在牛顿迭代法外层再涵盖一层用于更新D^{k ,nod}的非线性迭代步；节块耦合修正因子D^k,nod的计算与半解析节块法相同。

与现有技术相比，本发明具有如下几个突出特点：

1.统一求解联立的物理热工方程组，同时解得中子通量，慢化剂比焓，燃料温度。

2.将理论上具备二阶收敛速度的牛顿法用于物理热工耦合计算分析。

3.显示构造求物理热工非线性方程组的雅各比矩阵，避免了近似计算雅阁比矩阵对迭代收敛性的影响。

4.三重迭代格式，最外层为非线性迭代，计算节块耦合修正因子，中间层为牛顿迭代，计算物理热工耦合非线性方程组，内层为雅阁比矩阵求解的线性方程组迭代，非线性迭代可以保证中子方程的解收敛于高阶节块法，牛顿迭代则保证了耦合计算的精度和效率。

附图说明

图1三重迭代格式结构示意图。

具体实施方式

首先确定所需联立的方程形式，并建立封闭的物理热工联立方程组。

进行空间时间能量离散，缓发中子先驱核浓度进行时间上的线性近似后，节块内的2群时空中子扩散方程为

式中：

G——能群总数；

Δt——时间步长；

——第g群n时刻节块m内的平均中子通量；

——第g群n+1时刻节块m内的平均中子通量；

v_g——第g群中子速度；

χ_pg——瞬发中子裂变谱；

β_p——瞬发中子份额；

χ_dg——缓发中子裂变谱；

ν∑_f——中子产生截面；

∑_g'g——g’群到g群宏观散射截面；

∑_t——宏观总截面；

λ_l——瞬发中子份额；

C_l——第l群缓发中子先驱核浓度；

——由先驱核浓度线性近似推导得到的常数项系数；

——第g群节块m内的中子泄漏项；

——第g群节块m内u方向左右表面的中子静流；

Δu——u方向的节块宽度；

——节块等效扩散系数；

——非线性迭代节块耦合修正因子；

——m节块在u方向左侧面的不连续因子；

由缓发先驱核浓度方程(2)式可以看出，新一时刻的缓发先驱核浓度可在中子通量求出后代入得到，因此全耦合方法中无需将缓发先驱核浓度方程放入最终的大方程组中联立求解，所需求解的未知量为节块中子平均通量

此处需特别指出，对于非线性迭代节块耦合修正因子仅需知道它是通过对每个节块求解一个局部两节块问题得到的，其相关理论十分成熟，与非线性迭代半解析节块法完全一致。

基于单通道模型的离散后的热工方程组如下

(3)～(6)式依次为质量守恒方程，动量守恒方程，能量守恒方程，传热方程，式中

j——流体场网格标识；

i——传热场网格标识

A——通道流动面积；

ρ_j——第j网格冷却剂密度；

m_j——第j网格质量流量；

h_j——第j网格冷却剂比焓；

U'——冷却剂流速；

P——压力；

g——重力加速度；

f——摩擦系数；

φ²——两相摩擦倍增因子；

D_h——水力学直径；

ρ_l——冷却剂液相密度；

H——传热系数；

P_r——释热周长；

T_w——壁面温度；

T_b——冷却剂温度；

ρ_f——燃料密度；

c_p——燃料比热；

V——传热网格体积；

K——等效导热率；

q'''——体积释热率；

热工方程中未知量为h,m,P,T，若忽略压力对流体物性的影响，则由动量守恒方程(4)式可以看出，当求得h,m后，便可由(4)式求得压力P，从而全耦合方法最终联立的大方程组只包含(3)，(5)， (6)式，所需求解未知量为h,m,T。

经过以上讨论，将(2)，(3)，(5)，(6)式联立，并将每个方程的所有项移到等式同一端得到各个方程的残差形式，如下

如图1所示，全耦合方法中，最外层非线性迭代更新D^NOD，得到新的D^NOD值后，内层的牛顿迭代实际上求解了如下问题

其雅阁比矩阵为

根据方程性质，雅阁比矩阵实际上是稀疏的，这是由于中子平衡方程与质量流量m并无直接关系，因此雅阁比矩阵中元素值为零，同理，值为零的元素还有因此雅阁比矩阵中实际上需要计算的元素为

对于某一个牛顿迭代步，其计算逻辑为

J_F(x^k)δx^k＝-F(x^k)

x^k+1＝x^k+δ_x ^k

求解上式是一个线性方程组求解问题，当问题规模不大时可采用直接法，如LU分解，但一般说来对于堆芯量级的问题，雅阁比矩阵规模会特别巨大，可采用GMRES，双共轭梯度法等一系列数学上成熟的大型线性方程组数值方法求解。

当满足如下收敛判据时，牛顿迭代结束

||F(x)||<ε_Newton

当满足非线性迭代层收敛判据时，非线性迭代结束，至此全耦合法完成一次完整计算。

Claims (1)

1.一种稳定精确的反应堆物理热工耦合计算方法，其特征在于：包括如下步骤：

步骤1：确定中子学方程，即时空中子扩散方程的离散形式，对时空中子扩散方程采用非线性迭代半解析节块法进行空间上的离散，时间项采用隐式差分，能量取两群近似，得到带节块耦合修正因子项的节块平衡方程，即带修正因子的CMFD方程，其表达式如下：

式中，

G——能群总数；

Δt——时间步长；

——第g群n时刻节块m内的平均中子通量；

——第g群n+1时刻节块m内的平均中子通量；

v_g——第g群中子速度；

χ_pg——瞬发中子裂变谱；

β_p——瞬发中子份额；

χ_dg——缓发中子裂变谱；

ν∑_f——中子产生截面；

∑_g'g——g’群到g群宏观散射截面；

∑_t——宏观总截面；

λ_l——瞬发中子份额；

C_l——第l群缓发中子先驱核浓度；

——由先驱核浓度线性近似推导得到的常数项系数；

——第g群节块m内的中子泄漏项；

——第g群节块m内u方向左右表面的中子静流；

Δu——u方向的节块宽度；

——节块等效扩散系数；

——非线性迭代节块耦合修正因子；

——m节块在u方向左侧面的不连续因子；

确定热工方程的离散形式，对热工方程采用单通道模型，只考虑径向导热，空间上采用有限体积法离散，时间项采用隐式差分，得到有限差分格式的热工离散方程：

式中

A——通道流动面积；

ρ_j——第j网格冷却剂密度；

m_j——第j网格质量流量；

h_j——第j网格冷却剂比焓；

H——传热系数；

P_r——释热周长；

T_w——壁面温度；T_b——冷却剂温度；

ρ_f——燃料密度；

c_p——燃料比热；

V——传热网格体积；

K——等效导热率；

q'''——体积释热率；

步骤2：由于中子学方程与热工方程在数学表达式中并无显示的联系，因此需要对求解的节块平衡方程与热工离散方程建立数学关系，在节块平衡方程中截面是慢化剂密度和燃料温度的隐式函数∑_x＝g(ρ,T_f)，热工离散方程中体积释热项是关于中子通量和裂变截面的函数从而根据这两个关系，便能构建封闭的中子学方程和热工方程耦合的非线性方程组，形式如下：

式中，

f_φ(x),f_T(x),f_h(x),f_m(x)分别为残差形式的中子方程，传热方程，能量方程和质量方程；

F(x)便是联立的物理热工方程组，x为解向量，分别为中子通量，冷却剂比焓，冷却剂质量流量，燃料温度；

步骤3：建立了物理热工耦合的非线性方程组F(x)，在用牛顿法求解F(x)时，需在每一个牛顿迭代步构造F(x)的雅阁比矩阵，如下所示，

在雅阁比矩阵元素的计算中，需要求解中子残差方程对比焓的偏导数这实际上是要计算截面关于比焓的偏导数由于截面是关于慢化剂密度和燃料温度的隐式函数∑_x＝g(ρ,T_f)，因此使用求导的链式法则计算得到，

对于雅阁比矩阵中剩下的元素的计算，均能够通过各个残差方程得到解析表达式；

通过以上运算便得到了雅阁比矩阵各个元素的数值，为了获得每个牛顿迭代步的更新量δx，需要求解以下代数方程组，

J_Fδx＝-F(x)

J_F是一个大型的稀疏矩阵；

以k作为牛顿迭代步记号，当||F(x^k)||≤ε_Newton时认为牛顿迭代收敛，||F(x^k)||为F(x^k)的范数，ε_Newton为牛顿迭代收敛判据；

步骤4：物理热工耦合的非线性方程组F(x)求解完毕后，根据最新解得的节块中子通量更新节块耦合修正因子D^k,nod，因此是在牛顿迭代法外层再涵盖一层用于更新D^k,nod的非线性迭代步；节块耦合修正因子D^k,nod的计算与半解析节块法相同。