乏燃料溶液系统临界安全分析计算方法研究

银华北 王永平 苟军利 刘国明 祖铁军 尹 文 郑友琦 杜夏楠

1（西安交通大学 西安 710049）

2（中国核电工程有限公司 北京 100840）

为实现我国闭式循环目标，自主发展后处理技术、修建后处理厂具有重要意义。由于乏燃料后处理过程中核燃料会发生固-液相变、富集等复杂过程，涉及众多设备和环节，存在发生系统达临界事故的风险。因此，在乏燃料后处理厂中，对核临界事故的预防和防护，是保证核安全的关键内容。

在乏燃料溶液系统的临界事故中，临界瞬间功率快速上升，伴随剧烈的辐解、相变传质过程，针对固体燃料的核反应堆堆芯瞬态分析程序无法模拟上述过程。传统的核临界安全事故分析方法采用经验公式估计裂变次数，精度及适用范围有限。

针对溶液系统的上述特征，国内外机构开展了瞬态分析研究。英国原子能委员会首先发布了CRITEX程序，将溶液罐视作圆管，沿轴向分层进行一维热工水力-辐解气体计算，与点堆动力学耦合进行瞬态分析。法国与日本分别利用其溶液系统临界实验装置SILENE［1］与TRACY［2］进行实验研究溶液系统临界事故下的功率释放与溶液变化，并开发了AGNES2［3］、TRACE［4］、FECTH［5］等 程 序。其 中，TRACE与AGNES2程序由日本原子能研究所先后开发，其中子学均采用点堆动力学计算，而热工-水力采用R-Z二维轴对称模型，将溶液对流近似等效为导热，二者区别在辐解气体模型，TRACE程序基于CRITEX模型，将气泡速度近似为仅功率变化有关的函数，而AGNES2程序由气泡在溶液中的受力平衡导出气泡速度关系式；AGNES2程序在实验验证中精度优于TRACE程序。FETCH程序由法国开发，采用球谐函数展开下的有限元方法进行中子学计算，并直接使用两相流计算流体力学（Computational Fluid Dynamics，CFD）软件进行热工水力-辐解气体建模计算，且考虑了溶液流动对缓发中子先驱核分布的影响，其验证结果与实验值符合较好，然而由于采用了CFD模型，导致计算效率低，工程实用性差。

中国核电工程公司开发了CACCS程序［6］，动力学计算采用点堆模型，热工水力与辐解气体模型与AGNES2程序相似。清华大学开发了TCCHAR程序［7］，其中子学计算采用蒙特卡罗方法，热工模型采用集总模型，辐解气体模型在考虑产生、流入与流出的基础上还考虑了气泡大小在轴向的变化，其验证结果与实验值误差在15%以内。然而，在瞬态计算中，基于蒙特卡罗的中子学计算显然存在计算效率低、瞬态中的小反应性变化易被统计方差淹没的缺点。

可见，国内外现有程序中大部分采用点堆动力学进行中子学计算，使用余弦函数与贝塞尔函数近似系统内的功率分布与反应性反馈权重，其优势在于计算效率高，然而由于大量乏燃料后处理的需求，后续后处理工业示范厂规模大幅增加，势必加大溶液处理设备容量。而大容量溶液处理设备为保证系统安全大量布置控制毒物，并控制容器形状保证几何安全，上述假设不再适用，势必引入误差，对乏燃料存储策略、临界事故预防措施的制定带来不确定性；此外，随着计算机科学的发展，三维计算模型的计算时间逐渐下降至可接受范围内，基于时空动力学的溶液系统瞬态分析程序变得可能，从而进一步提高安全分析结论的可靠性。因此，开发一套几何适用范围广、模型精度高的临界安全分析程序，对于提高核临界事故的分析精度、保障乏燃料后处理厂的安全具有重要的意义。

为此，本研究开发了一套用于溶液系统的确定论并行三维瞬态分析程序Hydra-TD，主要包含截面制作模块，中子动力学模块及热工-辐解气体模块，其中中子动力学模块采用离散纵标并行输运程序计算，在保持三维中子学计算几何适用性广、分辨率高的优点的同时保证中子学计算效率；热工水力-辐解气体模块在AGNES2模型基础上进行简化以加速计算；热工水力-辐解气体-中子学耦合中考虑了溶液的膨胀、辐解气体对热散射的影响等因素，增加耦合精度。将程序与SILENE实验数据进行对比验证的结果误差较小，显示程序具备临界安全分析的能力。

本文对研究工作分为三个部分介绍，首先为程序的理论模型，介绍中子学截面制作、时空动力学模型和热工-辐解气体模型，以及中子学-热工水力-辐解气体模型的耦合方法；随后利用法国SILENE实验装置对本研究开发程序的模型进行了检验，其裂变次数、裂变率等结果精度较高；最后所得结论为本研究开发程序能够正确模拟瞬发临界第一裂变峰的功率变化，但无法还原功率下降阶段的振荡，辐解气体模型可能需要进一步改进。

1 理论模型

1.1 截面制作与产生

本文采用多群近似进行中子学计算，需要获取溶液系统各材料各个状态的有效多群截面。多群截面制作模块流程如下：

首先，使用核数据处理程序NECP-Atlas［8］，从评价核数据库出发，进行共振重构及线性化处理、多普勒展宽、不可分辨区有效自屏截面处理，热中子散射计算和多群常数计算，制作问题无关的多群截面库；然后，在非共振段使用温度和背景截面插值获得多群参数，在共振能量段使用超细群方法求解均匀问题的慢化方程，使用通量归并获得问题相关的多群参数。

式中：σi，k，g是第i区、k核素、多群g的有效自屏截面；σi，k，h是第i区、核素k、超细群h的截面；φi，h是第i区、超细群h的中子通量密度。

共振计算中，针对超细群方法的特点，模块采用了散射源加速及碰撞概率表插值以提高计算效率。散射源加速方法下，为避免计算源项时进行大量超细群的累加，每次计算时从上次的散射源中减去高能群超过最高对数能降的散射源，加上低能群的散射源：

式中：Sjg为j区域内其他能群散射至g能群的散射源；Δuf为对数能群宽度；Σsjk，g为j区域内k核素在g能群宏观散射截面；PNk，k为中子穿越Nk个能群后与核素k发生核反应的概率；φjg是j区域g群通量。

为提高计算效率，多群截面制作模块产生WIMS-69群各材料微观截面后，使用NECP-Hydra程序［9］进行一次溶液系统的稳态计算，使用全问题能谱将各材料截面归并为适用于热谱堆芯的4群2阶散射截面。在STACY燃料溶液罐问题中的稳态计算结果表明这一归并能够保持对系统有效增殖因子keff及溶液内通量分布的计算精度，keff偏差为十万分之十四，如表1所示；将溶液区轴向分为三层，径向分为4层，分群统计其通量分布（按蒙特卡罗程序输出的一个源中子的统计结果进行归一），所得分群分布偏差最大为5%，如表2所示。

表1 STACY问题keff计算结果对比Table 1 Comparison on calculated keff of STACY problem

表2 STACY问题溶液区分群通量计算结果对比Table 2 Comparison on calculated neutron flux distribution of STACY problem

针对溶液系统的辐射分解现象，为考虑辐解气体对热散射的影响，上述计算时同时制作了游离氢的截面数据，在瞬态计算时将对应的辐解气体由水中氢转化为游离氢。

1.2 中子时空动力学计算

本文在中子时空动力学计算中采用预估-校正准静态方法［10］，将中子通量密度与缓发中子先驱核密度分解为幅值与形状函数分别求解，在大时间步上求解空间形状分布，从而减少输运求解次数节省计算时间。

在多群、预估校正准静态的离散下，幅函数满足精确点堆方程，其具体形式见文献［10］，求解时，中子密度方程采用全隐式向后差分离散，缓发中子先驱核方程对其右端产生项βˉd(t)n(t)/Λ(t)采用亚当斯方法离散，即作二阶插值多项式在［tn，tn+1］内积分，其余项采用解析积分；中子通量密度形状函数的求解采用全隐式向后差分，需要求解的方程形式与有源次临界问题类似。

在输运方程的求解上，为保证计算效率，本文采用圆柱几何r-θ-z六面体网格下的KBA（Koch-Baker-Alcouffe）并 行离 散纵 标（SN）方法，基 于NECP-Hydra程序［9］完成开发。在KBA并行算法下，输运问题被空间区域分解，根据各个SN角度上的空间依赖关系，依次在各个并行单元上逐层逐角度启动扫描，如图1所示。KBA并行算法能避免迭代格式退化，是结构几何下并行效率最高的算法之一。在500×500×500网格的输运问题中，NECP-Hydra在百核下能保持90%的并行效率。

图1 KBA并行算法示意图Fig.1 Diagram of KBA parallel sweep algorithm

除NECP-Hydra中已经采用的切比雪夫外推加速、扩散综合加速方法［11］外，针对上述时间相关固定源问题，本文还采用了如下加速方法［12］：

1）对固定源项归一：

式中：S0为外中子源的归一化形状分布；·为在问题相空间区域内积分值。

2）进行类似于幂法的迭代过程直至迭代收敛。

式中：M为包含泄漏、消失与散射贡献的输运算子；F为裂变源算子。在瞬态计算过程中，每次输运计算的通量矩分布、k特征值都被保留，作为下次计算的迭代初值，在通量形状变化不大时可大幅减少迭代次数。

对于精确点堆方程系数及形状函数方程时间相关项中的角度相关项，本文采用与离散纵标方法内部储存一致的PN通量矩展开表示，即：

为使形式简洁，方程（6）省去了空间变量，方程（7）、（8）省去了时间和空间变量。其中：φk，lg为g能群k，l阶通量矩；Ylk(Ω)为对应的球谐函数；ψ*g为求解系统稳态共轭输运方程得到的共轭角通量；ψ͂g为归一化后的中子角通量形状函数。

由于瞬态过程中系统各处的溶液组分随时变化，为便于考虑溶液的各种物理变化，本研究使用微观截面进行中子动力学计算，在每个时间步按当地温度插值获得各核素的微观截面，随后结合当前区域的核子密度计算得到宏观截面，代入中子输运模块求解。

1.3 热工-辐解气体计算

在溶液系统临界过程中，由于功率升高，溶液温度上升，体积、浓度也发生变化。同时，裂变反应产生的高能碎片与燃料溶液中的分子发生碰撞，产生辐照裂解气体（H2、O2、N2等），并形成微泡；此外，中子和伽马射线一并与分子的相互作用对辐解气体的产生亦有贡献［1］。上述因素在瞬态过程中影响热散射、多普勒效应及核子密度分布，其计算精度关系到中子输运计算精度。

对于溶液中辐解气体计算，经调研国内外相关研究，本文参考AGNES2程序［3］的模型，建立如下气泡输运方程：

式中：Fi，j为（i，j）控制体内的空泡体积份额；νt为能量转换系数，m6·J-1·mol-1，Pi，j为控制体内功率密度，W·m-3；Ci，j为控制体内辐解气体浓度，mol·L-1；C0为辐解气体临界释放浓度；θ为单位阶跃函数；vi，j为气泡轴向迁移速度，m·s-1；G为辐解气体产生率，mol·J-1；τ为气体溶解时间。

上述方程中的未知量，除空泡体积份额、辐解气体量外，还包括气泡迁移速度。然而，溶液中辐解气体气泡的迁移过程受溶液对流、相间曳力、浮升力、气泡的生长、合并、溶解、破裂等多种因素影响，精确模拟其运动需进行精细建模的两相流数值计算，将极大影响计算效率。为此，本文中的辐解气体迁移速度采用近似关系式，通过一系列预先选取的两相流工况数值计算结果对当地功率密度、当地溶液动力学黏度与当地压力等因变量拟合得到。

对于系统热工计算，本文在溶液区及容器区分别计算导热，并在区域交界面上将对流、辐射换热等效为导热系数进行计算：

式中：ρ为材料密度；Cp为定压比热；T为温度；λeff为导热系数或等效导热系数；Φ̇为体积热源即功率密度。

1.4 耦合方法

在耦合计算中，中子动力学模块求解获得各个时刻的功率分布，热工-辐解气体模块求解获得各个时刻的温度、质量密度及辐解气体浓度分布，两者通过耦合接口按体积或质量权重转化为所需网格上的分布；质量密度与辐解气体浓度分布通过如下计算转化为中子学计算所需的核子密度。

式中：Ninuc为核素密度，l表示溶液初始包含核素的集合；εl为网格内溶液体积份额；G为辐解气体产生率，mol·J-1；Φ̇为功率密度；NA为阿伏伽德罗常数；δ为Kronecker符号，仅当下标相同时为1，其余时刻为0；l to gas表示溶液中辐解为气体的核素集合；nigas∈gas为辐解气体平均每分子内各核素的原子数；cgas为辐解气体的摩尔浓度分布。上式的物理含义为溶液在升温膨胀的同时被辐解使一部分核素析出并转移，其核素密度按物质的量混合至当地网格内，如图2所示。

此外，随溶液膨胀与辐解气体析出，溶液液位将发生变化。耦合计算中，热工-辐解气体模块根据溶液受热膨胀后密度与辐解气体所占体积得出总体积，从而得出溶液液位，耦合接口中结合液位在上溢网格内获得当地溶液与空气的体积份额，随后利用当地溶液质量密度获得当地核素密度，以考虑该效应：

式中：Voverflow为发生溶液上溢的网格；air表示气体包含核素的集合；εair为网格内空气体积份额；ρair为当地空气密度。式（14）和（15）的物理含义为：在溶液由于膨胀上溢处将溶液与空气按各自的体积份额混合至当地网格，如图2所示。

图2 中子学与热工水力-辐解气体耦合中的体积混合示意图Fig.2 Volume mixing operation when coupling neutronics with thermal-hydraulics / radiolysis gas

本文在耦合中采用算子分裂方法，即将温度、辐解气体场与中子场分离解耦，先后求解后更新至下一时间步。算子分裂方法是对多物理非线性耦合问题的线性化，其误差随时间步减小而减小。

本研究采用的计算流程如下：首先进行一次稳态前向计算，以获取初始边界条件；随后进行共轭计算，以共轭通量作为准静态下分离形状函数与幅函数的权重函数；进入瞬态计算后：

1）先更新截面进行预估求解，获取中子通量及先驱核浓度的形状函数；

2）随后在大时间步内进行形状函数插值更新点堆参数；

3）进行点堆求解获取校正通量及功率分布；

4）从而进行热工-辐解气体计算；

5）根据温度分布进行截面插值，根据质量密度及辐解气体含量更新核素密度分布从而归并得到宏观截面。

瞬态计算的时间步进流程如图3所示。

图3 瞬态时间步进流程Fig.3 Transient time stepping process of Hydra-TD

2 数值验证

按照§1的理论模型，本文编制了圆柱几何乏燃料溶液系统临界安全时空瞬态分析程序。为了验证该程序的计算精度，对SILENE实验装置的部分实验进行了建模计算。SILENE实验装置是一个硝酸铀酰溶液的均匀实验堆（235U富集度达92.7%），位于法国Burgundy的Valduc中心，主要用于临界事故研究，但还可以用于核加热、剂量测量、临界训练。堆芯在一个很大的混凝土房间中央，是一个小的环形水槽，如图4所示。不锈钢容器高度为1 m，外径360/368 mm，外壁厚度4 mm，内径70/76 mm，内壁厚度3 mm，底部厚度36 mm，顶部厚度30 mm。反应性引入棒为环形镉棒，厚度为1 mm，内外均有1 mm厚的不锈钢包壳，直径60/66 mm，长度为1 130 mm［1］。

图4 SILENE实验装置纵截面示意图Fig.4 Diagram of axial section of SILENE facility

计算设备采用Linux 18.04 LTS 64位系统，处理器是Intel Xeon Gold 6230 2.10GHz，内存空间为8 G，采用7个CPU核心进行计算。中子学计算中，划分的网格为r-θ-z方向24×1×210，角度采用层对称S4求积组。热工辐解气体计算中，对容器壁的网格进行细化，对容器上部空气部分网格进行粗化，其网格为r-z方向54×134。设置临界浓度为10 mol·m-3，辐解气体产生率为1.5×10-7mol·J-1，能量转换系数为1×10-7mol·J-1［3］。

本文所计算实验为S2-300实验，其初始液位高度为38.91 cm，通过以2 cm·s-1的速度提升反应性引入棒，共引入约0.97 $的反应性。对该问题进行时间步长敏感性分析，分别在脉冲阶段取大时间步为62.5 ms、31.25 ms、20 ms下计算，其裂变率曲线对比如图5所示，时间步为31.25 ms与20 ms的裂变率曲线已重合，可见在脉冲阶段取大时间步为31.25 ms能够使结果收敛。计算结果与实验测量结果的对比如图6所示，溶液高度、温度及空泡份额的变化和反应性变化绘于图7。随着反应性引入，系统功率开始指数上升，使溶液温度上升，密度下降，溶液高度升高；同时辐解气体开始积累，到达临界浓度后析出，使反应性快速下降。

图5 S2-300实验各时间步长模拟计算裂变率结果对比Fig.5 Comparison of calculated results of fission rate for S2-300 experiment with different time step length

图6 S2-300实验裂变率与温度变化计算结果（彩图见网络版）(a) 脉冲期间裂变率与温度变化，(b) 全过程裂变率与温度变化Fig.6 Calculated results of fission rate and temperature for S2-300 experiment (color online)(a) The change of fission rate and temperature in the pulse duration, (b) The change of fission rate and temperature in the whole process

计算获得的裂变率曲线趋势在第一功率峰处与实验值曲线一致，平均温度曲线位于实验值波动区间内。计算值与实验值的详细对比列于表2，峰值裂变率误差为22%，总裂变次数、峰值前裂变次数、功率倍增时间与峰值时间的误差分别为10%、2%、4%和3%，均在10%以内，以上结果表明，程序正确模拟了瞬态过程中溶液系统的复杂的多物理过程。

然而，由图6（b）中与实验值的比较结果可知，现阶段程序还无法精确模拟第一裂变峰后的功率振荡，因此未能模拟第二裂变峰，计算所得溶液温度在第二裂变峰后偏低，功率下降偏慢。由图7中空泡反应性反馈与空泡份额的变化可知，可能的原因是程序高估了空泡的反馈，或对气泡迁移速度存在低估，使负反馈偏强，这将在后续重点研究和改进。

图7 S2-300实验反应性与溶液变化计算结果Fig.7 Calculated results of reactivity and variation of solution system for S2-300 experiment

表3 S2-300实验功率变化特性计算结果Table 3 Calculated power feature of S2-300 experiment

3 结语

本文基于超细群共振方法、预估-校正准静态方法及导热微分方程与辐解气体输运方程，开发了一套用于圆柱几何溶液系统的临界安全分析的并行三维瞬态分析程序，并对SILENE实验装置的S2-300实验进行了建模计算，结果表明，程序能够正确模拟瞬发临界第一裂变峰的功率变化；程序现阶段无法还原S2-300的在第一裂变峰后的功率振荡现象，由溶液空泡份额下降较慢得知可能是辐解气体相关模型的影响，这部分内容将在后续工作中重点研究和改进。

致谢本文研究内容受国防科工局乏燃料后处理科研专项中“后处理厂核临界安全事故研究”项目资助，特此致谢。

作者贡献声明银华北：程序开发、验证计算、结果分析、起草文章；王永平：程序开发及理论指导，文章整体设计和修改；苟军利：热工水力-辐解气体模块开发；刘国明：提供实验数据，结果分析；祖铁军：截面模块程序开发及理论指导；尹文：截面制作程序模块开发；郑友琦：三维时空动力学程序开发理论指导；杜夏楠：程序截面模块验证指导。