CN107066745B - 获取快中子堆堆芯瞬态过程三维中子通量密度分布的方法 - Google Patents

Abstract

获取快中子堆堆芯瞬态过程三维中子通量密度分布的方法，对快堆堆芯各向异性的中子通量密度采用多边形棱柱网格进行全三维输运空间离散，并以六十度区域的半球面为单位进行交替扫描，减弱角度并行过程中的迭代格式退化；考虑快堆堆芯的中子通量密度局部效应较弱以及全局空间耦合的特点，对瞬态过程中子通量密度变化采用预估校正准静态策略进行时间离散，分解中子通量密度中随时间不同变化率的成分，并在不同时间尺度上分离求解，同时避免上述成分间的非线性迭代，提高计算效率；本发明中获取快中子堆堆芯瞬态过程三维中子通量密度分布的方法计算准确度高，计算量合理。

Description

获取快中子堆堆芯瞬态过程三维中子通量密度分布的方法

技术领域

本发明涉及快中子堆堆芯运行和安全技术领域，具体涉及获取快中子堆堆芯瞬态过程三维中子通量密度分布的方法。

背景技术

快中子堆是第四代核能系统中的主要堆型，其堆芯平均中子能量比传统热中子裂变占主导的压水堆高百万倍以上，同时，快中子引发裂变反应产生的中子数更多。富余的中子既可以增殖核燃料，提高铀资源的利用率；也可以用于嬗变核废料中高放射性长寿命的次锕系核素(镎、镅、锔等的同位素)，能够极大缩短核废料地质储藏所需的年限。因此，发展快中子堆是解决核能可持续发展瓶颈问题的重要途径。

快中子堆和压水堆具有不同的特性，从堆芯中子学和瞬态计算的角度出发，快中子堆最主要的特点是采用原子数更高的材料做冷却剂和载热剂，导致堆芯中子能谱相比于轻水冷却的压水堆偏硬许多。能谱的差异导致堆芯的设计理念以及瞬态安全相关的特征区别明显，比如：快中子堆堆芯进出口的温升更高，使得堆内材料膨胀引入的反应性更加重要；功率反应性系数相比轻水堆偏小；冷却剂的空泡效应可能为正。因此，快中子堆的瞬态中子通量密度的计算对其瞬态安全以及运行的经济性至关重要。

快中子堆的瞬态中子通量密度计算存在以下特征：堆芯中的反应性控制主要由控制棒提供，且控制棒占据整个组件范围，导致该组件附近的中子通量密度各向异性较强；较高的中子能量会导致散射反应的各向异性增加；对于时间相关特征，燃料装载形式以及硬中子能谱一方面使得堆芯的缓发中子份额减少，同时也使得中子代时间从数十微秒降低到约0.1到1微秒量级，瞬发中子和缓发中子的时间行为差距增大，系统的刚性增强。这些物理特征都对现有的快中子堆瞬态计算过程提出挑战，所以需要发明一种高精度的、计算效率高的获取快中子堆堆芯瞬态过程三维中子通量密度分布的方法。

发明内容

为了克服上述物理现象的困难、解决现有技术存在的问题，本发明提供了获取快中子堆堆芯瞬态过程三维中子通量密度分布的方法，通过分解中子通量密度中随时间不同变化率的成分，并在不同的时间尺度上分离求解，减少计算量。同时将上述过程应用于更严格的中子输运计算，从而保证计算精度。该方法能够用于快中子堆的瞬态安全分析，获取更加准确的设计裕量。

为了实现上述目的，本发明采用了以下技术方案予以实施：

获取快中子堆堆芯瞬态过程三维中子通量密度分布的方法，其中，所述快中子堆堆芯为典型的六角形几何排布或异形非结构几何排布，所述中子通量密度分布计算基于全三维输运计算，该方法包括以下步骤：

步骤1：读取快中子堆的瞬态计算所需参数，包括反应堆堆芯几何尺寸、堆芯材料布置、预制的均匀化截面参数、预制的均匀化动力学参数、堆芯入口温度和堆芯流量分布的边界条件；同时对所涉及的反应堆堆芯进行计算预处理，包括进行几何建模、划分计算区域、生成计算网格、离散角度空间，其中，针对不同几何排布的快堆堆芯采用相适应的三棱柱网格进行几何建模；

步骤2：指定各计算区域的材料，获取区域的宏观截面参数，设置反应堆堆芯计算的边界条件、计算迭代过程的初始值、计算迭代过程的收敛限值，计算稳态条件下各区域相关的中子通量密度分布，并且按照堆芯的输入功率水平进行归一化处理；

步骤3：进行快中子堆芯的热工水力计算，包括燃料棒的导热计算以及冷却剂通道与燃料棒的对流换热计算；以此得到均匀化截面以及动力学参数的反馈参数，包括燃料多普勒温度T_f.e以及冷却剂密度D_c；

步骤4：对步骤2和步骤3进行固定点迭代，直至反馈参数T_f,e和D_c收敛并满足设定的收敛限值；

步骤5：利用收敛后的截面参数进行稳态的特征值共轭方程计算，以得到稳态的共轭中子通量密度分布；

步骤6：将瞬态输运方程的中子通量密度以及缓发中子先驱浓度时间导数分别做准静态近似，即将其分解为各自的幅度函数和形状函数，其中，中子通量密度形状函数的空间导数保留角度相关性，时间导数采用各向同性近似简化角度相关性；利用收敛后的截面参数、动力学参数以及共轭中子通量密度获取上述幅度函数以及形状函数的初始值；

步骤7：根据瞬态过程的设定获取堆芯的扰动状态，如冷却剂温度变化、冷却剂流量变化以及控制棒位置变化，并基于当前的堆芯状态获得对应的截面参数和动力学参数；

步骤8：进行预估校正准静态计算，获取新一时刻t_n+1的中子通量密度以及缓发中子先驱核浓度分布；

1)根据上一时刻t_n的中子通量密度分布采用较大时间步长预估t_n+1时的中子通量密度；

2)由归一化条件生成t_n+1时刻的形状函数；

3)根据t_n和t_n+1时刻的形状函数在更小的时间步长上的t′时刻采用线性插值得到当前的形状函数，取消大时间步与小时间步长上计算的非线性迭代，并以此计算幅度函数方程组的系数矩阵；

4)利用生成的系数矩阵计算t′时刻的幅度函数；

5)校正通量密度：利用步骤4)中得到的幅度函数和步骤3)中得到的形状函数得到校正的中子通量密度分布；

6)利用校正的中子通量密度更新缓发中子先驱核浓度；

步骤9：获取新一时刻t_n+1的真实功率分布，进行格式化输出；

步骤10：重复步骤7至步骤9，直至瞬态过程结束。

与现有技术相比，本发明有如下突出优点：

1.借助准静态近似处理强刚性特征，取消幅度和形状函数间的非线性迭代，同时，采用中子通量密度形状函数的时间导数各向同性简化，提高计算效率；

2.利用勒让德-等权重求积组处理角度的六十度对称关系，同时基于线程级的角度并行技术提高计算效率；

3.采用正三棱柱和任意三棱柱网格进行几何建模，可以考虑各类快堆堆芯设计的几何差异，适应性强；

4.采用全三维中子输运计算得出瞬态过程中的三维空间的中子通量密度分布，可以考虑角度相关的中子通量密度随时间的变化率，计算精度高；

5.能量采用多群近似，显式区分瞬发和缓发中子的能谱差异，不必采用有效值形式的缓发中子份额。

附图说明

图1是快中子堆瞬态计算流程图。

图2是某快中子堆径向六分之一堆芯排布图。

图中标记：1.低富集度燃料组件；2.中富集度燃料组件；3.高富集度燃料组件；4.功率调节控制组件；5.停堆控制组件；6.第一排的屏蔽组件；7.其他的屏蔽组件；8.径向发射层组件。

图3是输运计算的角度离散图。其中图3-(a)是上半球的离散方法分布，图3-(b)是下半球的离散方向分布。

图中标记：1-48.离散方向编号。

图4是预估校正准静态方法时间步进的步骤图。

图5是本方法与传统方法(点堆近似模型)的效果比较图。

具体实施方式

本发明通过分解中子通量密度中随时间不同变化率的成分，并在不同的时间尺度上分离求解，极大降低计算量。同时将上述过程应用于更严格的中子输运计算，从而保证计算精度。方法流程如图1所示，该发明包括以下方面：

步骤1：读取快中子堆的瞬态计算所需参数，包括反应堆堆芯几何尺寸、堆芯材料布置、预制的均匀化截面参数、预制的均匀化动力学参数、堆芯入口温度和堆芯流量分布的边界条件。同时对所涉及的反应堆堆芯进行计算预处理，包括进行正三棱柱或任意三棱柱的几何建模、划分计算区域、生成计算网格、离散角度空间；

步骤2：指定各计算区域的材料，获取区域的宏观截面参数，设置反应堆堆芯计算的边界条件、计算迭代过程的初始值、计算迭代过程的收敛限值，计算稳态条件下各区域相关的中子通量密度分布，并且按照堆芯的输入功率水平进行归一化处理；

步骤3：进行快中子堆芯的热工水力计算，包括燃料棒的导热计算以及冷却剂通道与燃料棒的对流换热计算。以此得到均匀化截面以及动力学参数的反馈参数，包括燃料多普勒温度T_f.e以及冷却剂密度D_c；

步骤4：对步骤2和步骤3进行固定点迭代，直至反馈参数T_f,e和D_c收敛并满足设定的收敛限值；

步骤5：利用收敛后的截面参数进行稳态的特征值共轭方程计算，以得到稳态的共轭中子通量密度分布；

步骤6：将中子通量密度以及缓发中子先驱浓度的时间导数分别做准静态近似，即将其分解为各自的幅度函数和形状函数，形状函数的时间导数取各向同性。利用收敛后的截面参数、动力学参数以及共轭中子通量密度获取上述幅度函数以及形状函数的初始值；

步骤7：根据瞬态过程的设定获取堆芯的扰动状态，如冷却剂温度变化、冷却剂流量变化以及控制棒位置变化，并基于当前的堆芯状态获得对应的截面参数和动力学参数；

步骤8：进行预估校正准静态计算，获取新一时刻t_n+1的中子通量密度以及缓发中子先驱核浓度分布；

1)根据上一时刻t_n的中子通量密度分布采用较大时间步长预估t_n+1时的中子通量密度；

2)由归一化条件生成t_n+1时刻的形状函数；

3)根据t_n和t_n+1时刻的形状函数在更小的时间步长上的t′时刻采用线性插值取代非线性迭代得到当前的形状函数，并以此计算幅度函数方程组的系数矩阵；

4)利用生成的系数矩阵计算t′时刻的幅度函数；

5)校正通量密度：利用步骤4)中得到的幅度函数和步骤3)中得到的形状函数得到校正的中子通量密度分布；

6)利用校正的中子通量密度更新缓发中子先驱核浓度；

步骤9：获取新一时刻t_n+1的真实功率分布，进行格式化输出；

步骤10：重复步骤7至步骤9，直至瞬态过程结束。

步骤1中的典型快中子堆堆芯几何如图2所示，其计算的基本单元为六边形的快中子堆组件，空间计算可以采用任意形式的网格。

步骤2中求解的是稳态形式的多群中子输运方程，其来源如下。

瞬态过程中不仅需要考虑中子通量密度分布，同时需要考虑缓发中子浓度的分布，两者满足如下的关系式：

式中：r,Ω,t——空间变量、角度变量以及时间变量；g,h——能群标识；G——能群数目；m——缓发中子分组标识；M——缓发中子分组数；ψ_g(r,Ω,t)——中子角通量密度；Σ_t,g(r,t)，Σ_s,h→g(r,Ω′·Ω,t)，νΣ_f,h(r,t)——宏观总截面、宏观散射截面以及宏观中子产生截面；C_m(r,t)——缓发中子先驱核浓度；β_m(r,t)，λ_m(r,t)——缓发中子份额、缓发中子先驱核衰变常数；

——瞬发中子裂变谱、缓发中子裂变谱；ν_g(r,t)——中子速度；

——有效增殖因子；

考虑初始时刻，中子角通量密度以及缓发中子先驱核浓度随时间不变化，有：

于是，得到(初始时刻省略时间变量t)：

以及：

方程(4)为特征值方程，系统的功率由最低阶的特征向量以及归一化常数C决定。计算得到的中子通量密度为相对值，需要利用堆芯实际功率水平得到绝对中子通量密度

式中：κΣ_f,h(r)——宏观能量产生截面；P——堆芯设定的功率水平。

上述计算过程中采用离散纵标方法处理角度变量Ω。快中子堆的堆芯通常呈60度、120度、180度对称或者旋转对称，经典的层对称求积组的离散方向无法满足相应的对称关系，为了考虑这些情况，此处采用了勒让德-等权重求积组(P_N-EW,Nth-orderPolynomial Equal Weight)，其在极角方向仍为高斯-勒让德求积点，而幅角方向选用等间距等权重求积点。每一层的总权重为：

式中：P_N(x)——N阶勒让德多项式；ξ_i——N阶勒让德多项式P_N(x)的N/2个由小到大排列的正根。对于第i个极角层的第j的方向，其幅角为：

其中，j的最大值取12。对于P_N-EW求积组，S₆角度离散下各方向的具体分布如图3所示，从图中可以较为方便的获得60度区域的角度对称关系。

对z轴方向的每一层，P_N-EW求积组选取12个幅角点，即总的离散方向数目为12×N。考虑到各方向的泄漏项仅在边界条件处进行关联，即各方向上大部分的计算完全独立，故基于OpenMP进行了角度离散方向的并行化处理。角度扫描以60区域的半球面为单位进行，采用(+ξ)半球面、(-ξ)半球面、再(+ξ)半球面、再(-ξ)半球面的顺序扫描，可以保证始终存在两个60区域的角度是完全独立的。故对于S₆离散，每个扫描单位包含6个离散角度，可以保证小于等于12个线程同时计算时，迭代收敛过程和串行计算完全相同。

步骤3中的反馈参数选取燃料多普勒温度T_f,e以及冷却剂密度D_c，利用上述参数进行均匀化的截面参数以及动力学参数进行插值。其中，燃料多普勒温度的计算选取如下关系式：

式中：T_f(r)——导热计算得到的燃料芯块内的径向温度分布。

步骤4中的收敛准则如下：

式中：ε₁,ε₂——设置的燃料多普勒温度以及冷却剂密度的收敛限值；

——第n次迭代时的燃料多普勒温度；

——第n次迭代时的冷却剂密度。

步骤5中共轭中子通量密度满足的方程如下：

式中：

——共轭中子通量密度；

——共轭中子通量密度方程的特征值。

其中：

步骤5的求解可以采用步骤2中相似的过程。

步骤6中分别将中子通量密度和先驱核浓度分解成为幅度函数以及形状函数，并计算初始值。采用如下记号：

将先驱核浓度方程以及中子通量密度方程进行全隐式(FI,Fully Implicit)时间离散，得到：

同时，分别将中子通量密度和先驱核浓度分解成为幅度函数以及形状函数的形式：

式中：n(t)——中子通量密度的幅度函数；

——中子通量密度的形状函数；c_m(t)——缓发中子先驱核的幅度函数；——缓发中子先驱核的形状函数。

为了保证上式的唯一性，认为其满足如下条件：

将式(19)带入式(1)中并在两端乘以

将式(19)带入式(2)中并在两端乘以

在所有相空间积分，并利用共轭通量方程以及归一化条件。经过一定的变换，得到幅度函数满足的方程：

式中：ρ(t)——动态反应性；

——有效缓发中子份额；Λ(t)——有效中子代时间；

——缓发中子衰变常数。

计算的初始值包括中子通量密度和缓发中子先驱核浓度的真实分布、幅度函数以及形状函数：

1)真实分布：中子通量密度的初值通过求解稳态中子输运方程得到，同时令式(1)左端的时间导数项为零，得到缓发中子先驱核浓度的初值：

2)幅度函数：中子通量密度的幅度函数的初值取为n＝1，缓发中子先驱核的幅度函数的初值由式(21)时间导数为零得到：

3)形状函数：由真实分布和幅度函数直接得到，通量密度的形状函数为：

先驱核浓度的初始形状函数为：

步骤8中的预估校正准静态方法，需要计算的主要变量包括形状函数、幅度函数以及幅度函数方程组中的积分动力学参数。其中，形状函数随时间变化较为缓慢，幅度函数与时间紧密相关，而幅度函数方程组中的参数对时间的依赖程度介于两者之间，因此对应的，该方法中采用了三种尺度的时间步长类型，具体的时间离散如图4所示。图中，Δt和δt分别为形状函数和幅度函数方程组参数的计算时间步，依次称为大时间步(MATS,Macro TimeStep)、中间时间步(METS,Medium Time Step)。同时，幅度函数方程求解时通常将δt进一步的细分为小时间步(MITS,Micro Time Step)，但MITS通常隐含于对应的常微分方程组求解过程中，对于预估校正准静态的时间步进流程而言通常是不可见的。

图4中给出了具体的时间步进流程，如下：

1)计算通量密度：根据上一时刻中子通量密度分布ψ_g(r,Ω,t_n)预估t_n+1＝t_n+Δt时刻中子通量密度；

2)生成形状函数：由归一化条件生成t_n+1时刻的形状函数

3)插值形状函数：根据t_n和t_n+1时刻形状函数在METS上的t′时刻采用线性插值得到

并以此计算幅度函数方程组的系数矩阵；

4)计算幅度函数：利用生成的系数矩阵计算t′时刻的幅度函数；

5)校正通量密度：利用步骤4)中得到的幅度函数和步骤3)中得到的形状函数得到校正的通量密度函数，如下：

6)计算先驱核：利用校正的通量密度更新先驱核浓度；

重复步骤3)到步骤6)直至t′＝t_n+1时即完成了一个METS的步进过程。在下一个Δt的MATS上重复上面整个过程。利用预估的中子通量密度直接产生形状函数，再使用幅度函数校正预估得到的中子通量密度，这样的步进过程能够避免传统改进准静态方法中幅度函数和形状函数的非线性迭代过程，简化计算流程，提高计算效率。

步骤9中利用下式得到真实功率分布：

式中：C——稳态计算得到的功率归一系数。

为验证本发明的有效性，采用如图2所示的典型快中子堆控制棒弹出瞬态过程算例。假设堆芯中的一组控制棒在2秒的时间内提出12厘米，采用本发明中的获取瞬态过程三维中子通量密度分布方法以及传统点堆近似方法分别计算，结果如图5所示。两者的具体数值存在一定差异，本发明中的方法考虑了更加精确的空间效应以及角度的输运效应，能够得到更加准确的控制棒引入的反应性以及瞬态过程规律。

Claims (1)

1.获取快中子堆堆芯瞬态过程三维中子通量密度分布的方法，其中，所述快中子堆堆芯为典型的六角形几何排布或异形非结构几何排布，所述中子通量密度分布计算基于全三维输运计算，其特征在于，包含以下步骤：

步骤1：读取所述快中子堆的瞬态计算所需参数，包括反应堆堆芯几何尺寸、堆芯材料布置、预制的均匀化截面参数、预制的均匀化动力学参数、堆芯入口温度和堆芯流量分布的边界条件；同时对所涉及的反应堆堆芯进行计算预处理，包括进行几何建模、划分计算区域、生成计算网格、离散角度空间，其中，针对不同几何排布的快堆堆芯采用相适应的三棱柱网格进行几何建模；

步骤2：指定各计算区域的材料，获取区域的宏观截面参数，设置反应堆堆芯计算的边界条件、计算迭代过程的初始值、计算迭代过程的收敛限值，计算稳态条件下各区域相关的中子通量密度分布，并且按照堆芯的输入功率水平进行归一化处理；

步骤3：进行快中子堆芯的热工水力计算，包括燃料棒的导热计算以及冷却剂通道与燃料棒的对流换热计算；以此得到均匀化截面以及动力学参数的反馈参数，包括燃料多普勒温度T_f,e以及冷却剂密度D_c；

步骤4：对步骤2和步骤3进行固定点迭代，直至反馈参数收敛并满足设定的收敛限值；

步骤5：利用收敛后的截面参数进行稳态的特征值共轭方程计算，得到稳态的共轭中子通量密度分布；

步骤6：将瞬态输运方程的中子通量密度以及缓发中子先驱浓度的时间导数分别做准静态近似，即将其分解为各自的幅度函数和形状函数，其中，中子通量密度形状函数的空间导数保留角度相关性，时间导数采用各向同性近似简化角度相关性；利用收敛后的截面参数、动力学参数以及共轭中子通量密度获取上述幅度函数以及形状函数的初始值；

步骤7：根据瞬态过程的设定获取堆芯的扰动状态，并基于当前的堆芯状态获得对应的截面参数和动力学参数；

步骤8：进行预估校正准静态计算，获取新一时刻的中子通量密度以及缓发中子先驱核浓度分布；

1)根据上一时刻t_n的中子通量密度分布采用较大时间步长预估t_n+1时的中子通量密度；

2)由归一化条件生成t_n+1时刻的形状函数；

3)根据t_n和t_n+1时刻的形状函数在更小的时间步长上的t′时刻采用线性插值得到当前的形状函数，取消大时间步与小时间步长上计算的非线性迭代，并以此计算幅度函数方程组的系数矩阵；

4)利用生成的系数矩阵计算t′时刻的幅度函数；

5)校正通量密度：利用步骤4)中得到的幅度函数和步骤3)中得到的形状函数得到校正的中子通量密度分布；

6)利用校正的中子通量密度更新缓发中子先驱核浓度；

步骤9：获取新一时刻的真实功率分布，进行格式化输出；

步骤10：重复步骤7至步骤9，直至瞬态过程结束。

Images