固定网格下考虑变形效应的反应堆瞬态中子通量模拟方法

1.本发明涉及核反应堆瞬态分析计算技术领域，具体涉及一种在固定网格下考虑变形效应的反应堆瞬态中子通量模拟方法。


背景技术：

2.在反应堆瞬态过程中，由于功率、温度的变化，堆芯会不可避免地发生变形现象，如燃料棒伸长、组件弯曲、堆芯径向膨胀等。尤其对快堆来说，堆芯变形效应是瞬态中不可忽略的现象，它对瞬态中功率、温度的变化影响巨大。因此，传统方法仅采用中子物理-热工耦合分析无法准确刻画瞬态现象，需要耦合变形效应。然而，由于变形效应使得核反应堆几何发生变化，导致瞬态计算方法中的几何建模、网格划分、不同物理场映射变得十分困难，成为了当前国际上快堆瞬态计算分析的“卡脖子”难题。
3.目前，在瞬态分析中考虑变形效应有两类方法：三维直接耦合方法和点堆计算方法。三维直接耦合方法就是将力学计算程序与中子学及热工计算程序耦合，但由于变形效应会导致计算问题几何发生改变，因此物理、热工、力学程序都需要采取细网程序进行耦合，以直接模拟变形后的几何结构。例如，美国阿贡国家实验室开发了名为sharp的程序系统，它由物理程序protues，热工分析程序nek5000和力学分析程序diablo组成，程序基于有限元细网，直接进行精细几何建模计算。它的优点是精度高，然而缺点也十分明显：各物理场均为精细网格模型，计算量巨大，程序系统计算效率十分低下；同时变形导致的网格映射关系复杂、程序鲁棒性差。因此sharp程序实际只能在某一稳定状态下进行核-热-力耦合计算，无法真正应用到堆芯的瞬态分析中去。点堆计算方法基于最简化的点堆模型，只考虑整体效应，无法模拟局部功率变化等效应，计算效率高，但计算精度低、适用范围有限，例如美国阿贡国家实验室开发的快堆瞬态计算程序sas4a/sassys便是采用这种方法。因此，国际上尚无能在固定网格下考虑反应堆变形效应的瞬态中子通量方法。


技术实现要素：

4.为了克服上述现有技术存在的问题，本发明的目的在于提供一种固定网格下考虑变形效应的反应堆瞬态中子通量模拟方法，预先计算获得三类反应堆堆芯变形反应性反馈系数，结合某一时刻堆芯的实际变形量，计算获得总的堆芯变形反应性反馈量，将其加至中子幅值函数中反应性变量的计算公式中，从而考虑堆芯变形效应对中子通量分布的影响，提高数值计算模拟精度。
5.为了实现以上目的，本发明采取如下的技术方案予以实施：步骤1：将反应堆堆芯变形分解为燃料轴向伸长、燃料组件弯曲和堆芯径向膨胀三类变形现象，并在稳态情况下，通过直接法或微扰理论分别求解上述三类反应堆堆芯变形反应性反馈系数，即单位伸长量、单位弯曲量和单位径向膨胀量引起的反应性变化量；步骤2：在瞬态计算的某一时刻，计算堆芯的三维变形场，结合步骤1中获得的三类反应堆堆芯变形反应性反馈系数便能获得堆芯总变形反应性反馈量：
ꢀꢀ
(1)式中： —— 由于变形效应导致的堆芯总反应性反馈量； —— 空间网格i处的材料在方向发生单位变形的堆芯反应性反馈系数；—— 空间网格i处的材料实际在方向的变形量；—— 空间网格总数；步骤3：在预估校正准静态时空动力学方法框架下，在计算幅值函数中的反应性变量时，加入步骤2得到的堆芯总变形反应性反馈量，从而在瞬态计算中耦合变形效应，最终获得准确的中子通量；三维时空动力学方程组为：
ꢀꢀꢀꢀꢀ
(2) (3)式中： —— 第g能群中子速度； —— 在位置、角度和t时刻时，第g能群的中子通量； —— 在位置和t时刻时，第g能群中子宏观总截面； —— 在位置和t时刻时，处于角度和第h能群的中子散射至角度和第g能群的中子宏观散射截面； —— 在位置、角度和t时刻时，第h能群的中子通量； —— 在位置和t时刻时，瞬发中子在第g能群的中子能谱；
ꢀ
—— 在位置和t时刻时，缓发中子的总份额；—— 稳态中子增殖系数； —— 在位置和t时刻时，第h能群中子宏观产生截面； —— 在位置和t时刻时，第j组缓发中子在第g能群的中子能谱；——在位置和t时刻时，第j组缓发中子先驱核衰变常数；——在位置和t时刻时，第j组缓发中子先驱核的密度； —— 在位置和t时刻时，第j组缓发中子的份额；j —— 缓发中子先驱核总组数；g—— 中子总能群数；预估校正准静态时空动力学方法中，对中子通量方程(2)和先驱核浓度方程(3)进行全隐式时间差分离散；同时，采用准静态近似下的因子分裂，将三维时空动力学方程的分布函数分为强时间相关的幅值函数部分与弱时间相关的形状函数部分：
ꢀꢀꢀꢀꢀꢀꢀꢀ
(4)式中： —— t时刻中子通量的幅值函数； —— t时刻第j组缓发中子先驱核的幅值函数； —— 在位置、角度和t时刻时，第g群中子通量的形状函数； —— 在位置和t时刻时，第j组缓发中子先驱核的形状函数；最终中子通量方程（2）和先驱核浓度方程（3）变换为求解形状函数的方程和求解幅值函数的方程；其中，求解幅值函数的方程为：
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(5)
式中： —— t时刻堆芯总反应性； —— t时刻堆芯平均缓发中子总份额；—— t时刻堆芯中子代时间； —— t时刻第j组缓发中子先驱核堆芯平均衰变常数； —— t时刻第j组缓发中子的堆芯平均份额；上述参数计算公式为：
ꢀꢀ
(6)
ꢀꢀ
(7)
ꢀꢀꢀꢀꢀꢀꢀꢀ
(8)
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(9) (10)
ꢀꢀꢀ
(11)式中： —— 堆芯活性区； ——在位置处，t时刻与初始时刻之间第g能群宏观总截面之差； —— 在位置、角度处第g能群的稳态共轭中子通量；
ꢀ
—— 在位置处，处于角度和第h能群的中子散射至角度和第g能群的中子宏观散射截面在t时刻与初始时刻之差； —— 在位置、角度和t时刻时，第h能群中子通量的形状函数；—— 在位置和t时刻时，第g能群的总中子能谱；——在位置处，t时刻与初始时刻之间第h能群中子宏观产生截面之差；公式（6）为堆芯总反应性计算公式，表征了堆芯内由于控制棒移动、材料温度变化原因导致的堆芯截面变化引起的反应性变化，该反应性变化直接影响堆芯中子通量变化速率；针对快堆特性，对变形效应进行耦合，即将堆芯反应性计算公式在公式（6）基础上增加一项，写为：
ꢀꢀꢀꢀꢀꢀꢀ
(12)其中，由步骤2中公式（1）计算获得，公式（12）的物理含义即为：瞬态中堆芯总反应性由两部分组成，一方面是由于控制棒移动、材料温度变化原因导致的堆芯截面变化引起的反应性变化；另一方面是由于堆芯变形效应引起的反应性变化；以上两方面将同时影响中子通量幅值，进而影响后续时间步的通量形状，最终获得精确的中子通量。
6.优选的，利用多物理场仿真软件comsol计算堆芯的三维变形场。
7.与现有技术相比，本发明有如下优点：本发明利用三维时空动力学计算方法模拟瞬态中子通量变化，能在固定网格下考虑反应堆的变形效应。相比于三维直接耦合方法，本发明避免了在每个时间步重新划分网格和物理量映射，显著提升了计算效率和工程实用性；相比于点堆计算方法，本发明保持了三维时空的模拟能力，能保证模拟的精度。因此，本发明将能实现考虑变形效应的反应堆瞬态中子通量高精度、高效率模拟，是瞬态计算中耦合变形效应的全新思路。
附图说明
8.图1为 固定网格下考虑变形效应的反应堆瞬态中子通量模拟方法总体流程图。
9.图2a为反应堆堆芯燃料轴向伸长示意图。
10.图2b为反应堆堆芯燃料组件弯曲示意图。
11.图2c为反应堆堆芯径向膨胀示意图。
具体实施方式
12.下面结合附图和具体实施方式对本发明做进一步详细说明：如图1所示，本发明固定网格下考虑变形效应的反应堆瞬态中子通量模拟方法，包括如下步骤：步骤1：将反应堆堆芯变形分解为燃料轴向伸长、燃料组件弯曲和堆芯径向膨胀三类变形现象，分别如图2a、图2b和图2c所示。在稳态情况下，通过直接法或微扰理论分别求解上述三类反应堆堆芯变形反应性反馈系数，即单位伸长量、单位弯曲量和单位径向膨胀量引起的反应性变化量；步骤2：在瞬态计算的某一时刻，用多物理场仿真软件comsol对堆芯进行建模并计算三维变形场，获得相比于初始状态，堆芯中每一个空间网格内材料轴向伸长量、弯曲量及方向、径向膨胀量，例如堆芯中间组件的最上端网格在轴向伸长了2毫米，向堆芯外法线方向弯曲了1毫米，同时径向膨胀了3毫米。基于所有网格的变形量，并结合步骤1中获得的三类反应堆堆芯变形反应性反馈系数便能由下式获得堆芯总变形反应性反馈量：
ꢀꢀꢀꢀ
(1)式中： —— 由于变形效应导致的堆芯总反应性反馈量； —— 空间网格i处的材料在方向发生单位变形的堆芯反应性反馈系数；—— 空间网格i处的材料实际在方向的变形量；—— 空间网格总数；步骤3：在预估校正准静态时空动力学方法框架下，在计算幅值函数中的反应性变量时，加入步骤2得到的堆芯总变形反应性反馈量，从而在瞬态计算中耦合变形效应，最终获得准确的中子通量；例如，稳态反应堆堆芯中的控制棒突然上提，随后反应堆中子通量的幅值和形状都会发生变化，为了定量计算中子通量随时间和三维空间的变化关系，需求解以下三维时空动力学方程组：
ꢀꢀꢀ
(2)
ꢀ
(3)式中： —— 第g能群中子速度； —— 在位置、角度和t时刻时，第g能群的中子通量； —— 在位置和t时刻时，第g能群中子宏观总截面； —— 在位置和t时刻时，处于角度和第h能群的中子散射至角度和第g能群的中子宏观散射截面； —— 在位置、角度和t时刻时，第h能群的中子通量； —— 在位置和t时刻时，瞬发中子在第g能群的中子能谱； —— 在位置和t时刻时，缓发中子的总份额；—— 稳态中子增殖系数； —— 在位置和t时刻时，第h能群中子宏观产生截面； —— 在位置和t时刻时，第j组缓发中子在第g能群的中子能谱；——在位置和t时刻时，第j组缓发中子先驱核衰变常数；——在位置和t时刻时，第j组缓发中子先驱核的密度； —— 在位置和t时刻时，第j组缓发中子的份额；j —— 缓发中子先驱核总组数；g—— 中子总能群数；预估校正准静态时空动力学方法中，对中子通量方程（2）和先驱核浓度方程（3）进行全隐式时间差分离散；同时，采用准静态近似下的因子分裂，将三维时空动力学方程的分布函数分为强时间相关的幅值函数部分与弱时间相关的形状函数部分：
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(4)式中：
ꢀ
—— t时刻中子通量的幅值函数； —— t时刻第j组缓发中子先驱核的幅值函数； —— 在位置、角度和t时刻时，第g群中子通量的形状函数； —— 在位置和t时刻时，第j组缓发中子先驱核的形状函数；最终中子通量方程（2）和先驱核浓度方程（3）变换为求解形状函数的方程和求解幅值函数的方程；其中，求解幅值函数的方程为：
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(5)式中： —— t时刻堆芯总反应性； —— t时刻堆芯平均缓发中子总份额； —— t时刻堆芯中子代时间； —— t时刻第j组缓发中子先驱核堆芯平均衰变常数； —— t时刻第j组缓发中子的堆芯平均份额；上述参数计算公式为： (6)
ꢀꢀ
(7)
ꢀꢀꢀꢀ
(8)
ꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀꢀ
(9)
ꢀꢀ
(10)
ꢀꢀꢀꢀ
(11)式中： —— 堆芯活性区； ——在位置处，t时刻与初始时刻之间第g能群宏观总截面之差； —— 在位置、角度处第g能群的稳态共轭中子通量； —— 在位置处，处于角度和第h能群的中子散射至角度和第g能群的中子宏观散射截面在t时刻与初始时刻之差； —— 在位置、角度和t时刻时，第h能群中子通量的形状函数；—— 在位置和t时刻时，第g能群的总中子能谱；——在位置处，t时刻与初始时刻之间第h能群中子宏观产生截面之差；公式（6）为堆芯总反应性计算公式，表征了堆芯内由于控制棒移动、材料温度变化原因导致的堆芯截面变化引起的反应性变化，该反应性变化将直接影响堆芯中子通量变化速率，但公式（6）并未考虑堆芯变形效应对反应性变化的影响。本发明对变形效应进行耦合，即将堆芯反应性计算公式在公式（6）基础上增加一项，即步骤2中公式（1）计算获得的堆芯总变形反应性反馈量，最终将公式（6）写为：
ꢀꢀꢀ
(12)公式（12）的物理含义即为：瞬态中堆芯总反应性由两部分组成，一方面是由于控制棒移动、材料温度变化原因导致的堆芯截面变化引起的反应性变化；另一方面是由于堆
芯变形效应引起的反应性变化。以上两方面将同时影响公式（5）计算中子通量幅值，进而影响后续时间步的通量形状，最终获得控制棒提升的瞬态过程中，考虑变形效应的、更为准确的中子通量。