基于蒙特卡罗均匀化理论与有限体积方法的溶液系统临界事故分析方法

孙 旭，周 琦,*，余慧莺，朱庆福，夏兆东，宁 通，马骁笛

(1.中国原子能科学研究院 反应堆工程技术研究所，北京 102413；2.国家电投集团科学技术研究院有限公司，北京 102209)

在核燃料循环领域，易裂变核材料在某些环节上是以溶液的形式存在，溶液中的水是中子的优良慢化材料，使得溶液中具有较少的易裂变核材料便可达到临界。同时溶液具有易流动、易变形的特点，以及多相铀钚混合体系、多体相互作用等使得核材料的临界安全控制更加复杂[1]。据文献报道，各国在核燃料循环领域共发生过22起临界事故，其中有21起发生在溶液或浆液中[2]。临界事故不仅会引起国际社会的广泛关注，严重影响核能发展，还会造成大量放射性物质的释放、人员伤亡和环境污染。因此，有必要对溶液系统临界事故进程进行研究，从而对临界事故的预防、屏蔽设计及后果评价提供指导。

自2011年以来，中国原子能科学研究院一直承担核燃料系统瞬发临界事故分析过程与方法研究的任务，研制了包括GETAC-M(固体)、GETAC-S(溶液)、GETAC-WP(湿粉末)及GETAC-DP(干粉末)等一系列核燃料系统瞬态分析程序[3-4]。包括GETAC-S在内的大多溶液系统瞬态分析程序使用点堆模型对事故进程进行估计，然而核燃料循环体系中溶液系统存在多种几何形态，且存在多相溶液的不均匀混合等特征，点堆模型难以适应复杂的计算条件。本文基于蒙特卡罗均匀化理论与有限体积方法，提出适用于瞬态分析问题的三维扩散时空动力学模型，将该模型与非稳态传热模型、辐照裂解气泡模型相耦合，对GETAC-S进行升级。

1 基础理论

核溶液系统发生临界事故后，系统功率将在短时间内迅速上升得到一个功率峰值，而后温度反馈效应和辐照裂解气泡反馈效应会使系统功率回到较低水平。随温度下降和气泡消失，这些反馈效应的减小和消失，系统功率会重新达到一个比第一个功率峰值低的值，如此反复形成功率震荡，直至部分水被蒸发或易裂变材料溅出而使系统变为次临界状态停止功率震荡[5]。为了模拟临界事故的发展过程，需要从中子物理和热工水力反馈两方面进行分析。物理上采用三维扩散时空动力学模型可为热工计算提供更加准确的功率分布，从而对事故进程的发展给出更为精确的模拟结果。

三维扩散时空动力学的求解通常需要将系统均匀化近似后得到多群常数作为输入条件，本文采用蒙特卡罗方法产生均匀化多群常数，并采用有限体积方法进行三维扩散时空动力学的求解。

1.1 蒙特卡罗均匀化理论

蒙特卡罗方法一般使用连续能量截面，对三维实际问题建立精确的几何模型，并考虑各向异性散射，相较于确定论输运计算方法，在能量、空间和角度等方面的近似更少，能够准确描述非均匀性几何、边界条件、共振等物理特性，是较为精确的输运方法。

蒙特卡罗均匀化的核心是将求解输运问题得到的各种事件统计结果转化为均匀化常数，这些事件主要包括散射、裂变、俘获吸收、(n，xn)等反应，均匀化常数包括各种反应的群截面、输运修正、散射矩阵、各阶勒让德项，次级中子产生矩阵、裂变产额、裂变能谱，以及动力学计算需要的中子速度倒数、缓发中子先驱核产额与衰变常量等。

为简化表达式，引入内积算符〈·，·〉来代表能量、空间或角度的积分形式，这样第x类核反应的反应率计数估计可表示为：

(1)

式中：Σx为第x类核反应的截面；φ为中子通量密度；V为体积；S为角度；E为中子能量；r为空间位置；Ω为中子运动的方向。

内积函数的统一性，如空间均匀化与能量积分通量可表示为：

〈φ〉k,g≡〈φ,1〉k,g=

(2)

式中：k为空间网格编号；g为能群编号；Vk为第k网格体积；Eg为第g能群的能量上限。

(3)

1.2 多群时空中子动力学扩散方法

有限体积方法是一种处理偏微分方程的数值方法，在热工、流体、力学等领域有着广泛应用，将有限体积方法用于中子扩散方程的求解，能够更加有效地完成中子物理参数与其他参数的交互，实现耦合与反馈计算，提高反应堆物理分析的准确性。

多群时空中子动力学扩散方程如式(4)所示，本文使用有限体积方法对其进行离散求解[6]。

g=1,…,G

(4)

i=1,…,ND

(5)

式中：vg为第g群中子运动速度；φg为第g群中子通量密度；t为时间；Dg为第g群中子的扩散系数；Σr，g为第g群中子的移出截面；Σs，g′→g为第g′群中子到第g群的散射截面；χg为第g群中子的裂变能谱；β为缓发中子份额；(νΣf)g′为第g′群中子的平均裂变中子数与裂变截面的乘积；i为缓发中子先驱核组编号，组数为ND；λi为第i组先驱核衰变常量；Ci为第i组先驱核浓度；G为能群数。

2 软件实现

使用蒙特卡罗均匀化方法进行多群常数的生成，然后采用有限体积方法进行三维扩散中子时空动力学行为研究，采用非稳态传热模型与辐照裂解气泡模型进行热工反馈分析。

首先需要进行稳态中子扩散方程的求解(临界源问题)，得到初始状态下求解区域内的中子通量密度分布；然后进行非稳态扩散方程的求解(固定源问题)，以得到随时间变化的中子通量密度与缓发中子先驱核浓度，再通过将热工反馈计算得到的温度、空泡系数反馈到非稳态中子扩散方程的求解中，反复迭代，直到模拟时间结束。

2.1 基于蒙特卡罗方法的多群常数生成

蒙特卡罗程序可通过计数功能进行多群常数的生成。这些多群截面可提供给各种能群结构或空间网格划分的确定论程序。从能量上看，蒙特卡罗程序可计算单独的核素或元素的宏观或微观核反应率，可得到一个或多个任意能群结构的多群截面，便于研究能群结构对多群截面精度的影响；从空间上看，蒙特卡罗程序还能对各种空间划分方式进行多群截面的计算，包括材料、栅元以及区域相关的笛卡尔坐标系下的空间划分，还支持重复结构系统，如组件或堆芯的计算。利用蒙特卡罗方法生成多群常数，满足软件适用于复杂几何、材料以及能谱结构的目标要求。采用OpenMC程序计算了多群截面、瞬发中子常数以及缓发中子常数等参数[7]。

2.2 基于OpenFOAM的三维扩散时空动力学程序

三维扩散时空动力学程序是在使用有限体积方法求解偏微分方程的OpenFOAM的基础上进行开发[8]。OpenFOAM是一个完全面向对象的开源CFD类库，这些库全部是由C++语言编写，OpenFOAM针对偏微分方程具有多种离散求解方式，这些求解方法同样适用于三维扩散时空动力学方程的求解。

基于OpenFOAM求解三维扩散时空动力学方程的基本流程如图1所示。主要分为3部分，分别是软件输入、网格生成以及扩散方程的求解。其中输入主要分为6部分：neutronProperties、topoSetDict、controlDict、decomposeParDict、blockMeshDict及phi.group，用于初始条件及控制条件等参数的输入。网格生成用于将整个计算空间离散为网格，每个网格内选择1个点作为此网格的代表。由于蒙特卡罗均匀化理论已将复杂的几何模型了进行简化，这里使用正方形网格划分即可满足要求。离散后的扩散方程求解使用OpenFOAM提供的函数库实现。

图1 基于OpenFOAM求解时空动力学方程的基本流程Fig.1 Basic process of solving spatiotemporal dynamics equation based on OpenFOAM

三维扩散时空动力学程序的输出主要分为3部分，分别为keff、通量分布以及功率分布。瞬态过程的计算首先需要求解稳态扩散方程，将稳态计算通量分布用于瞬态的计算条件。因此，内核方程的求解分为两种情形：不含时间项的稳态中子扩散方程的求解与含时间项的中子扩散时空动力学方程(式(4)、(5))的求解。

2.3 热工反馈与耦合的实现

式(6)、(7)分别为热工反馈采用非稳态传热模型和辐照裂解气泡模型。

(6)

式中：ρ为介质密度，kg/m3；c为比热容，J·kg-1·K-1；T为温度，K；λ为介质导热系数，W·m-1·K-1；S为内热源，对核燃料溶液系统即为功率密度，W/m3；h为轴向上的坐标；r为径向上的坐标。

(7)

式中：V(t，z)为t时刻、轴向位置z处的空泡体积份额；ζ为氢气产量与其他气体产量之比；G(H2)为单位能量的氢气产额，mol/J；P(t，z)为t时刻、z位置处的功率密度，W/m3；R为摩尔气体常数；Tg为气泡内部温度，K；p0为环境大气压强，Pa；g为重力加速度，m2/s；H(z)为z位置上方溶液的高度，m；σ为溶液的表面张力系数，N/m；r为气泡半径，m；v(z)为z位置处的气泡向上运动速度，m/s，参考CRITEX的两速度模型，设定功率上升时，气泡的上升速度为4.0×10-2m/s，功率下降时，气泡的上升速度为1.5×10-2m/s。

基于OpenFOAM平台，对上述两项热工反馈模型进行离散求解[9]，并实现物理计算与热工反馈计算的网格匹配。由于热工反馈计算得到的参数为温度与空泡份额，GETAC-S通过设置温度与空泡反馈系数将其传递到反应性的变化中。在升级的程序中，由于使用三维扩散时空动力学模型代替点堆模型，而多群截面、动力学参数等作为扩散软件的输入代替作为点堆模型的反应性输入。因此，使用原来的反馈系数传递反应性的方法不能充分体现出三维扩散时空动力学的优势。

采用以温度与空泡份额为自变量，将截面预先制成相关函数，形成函数库(或拟合成曲线)，在每一时间步长内热工反馈及计算结束后，得到相应的温度与空泡份额，然后调用相关函数库(或拟合曲线)得到对应温度与空泡份额下的少群截面与动力学参数，将其代入到时空动力学方程中，实现每一时间步长内的迭代求解。

3 TRACY实验验证

使用日本TRACY瞬态实验装置的模拟临界事故实验数据对升级的GETAC-S程序进行可靠性验证[10]。TRACY装置的主体是一个带有中间孔道的环形不锈钢容器，如图2所示。实验时可通过拔出中间孔道内的碳化硼控制棒从而模拟不同反应性引入所引起的瞬时功率上升事故。TRACY装置可记录瞬态事故发展过程中功率、温度等参数的变化，本文对其中的RUN72实验[11-12]进行了验证，表1列出RUN72实验模拟参数，表2列出RUN72实验计算结果及其相对偏差。由表2可看到，程序计算的功率峰值、总能量和最终温度与实验结果的相对偏差均在10%以内，结果均符合良好，验证了升级后的GETAC-S在计算结果上有较大改进(改进前GETAC-S计算结果与实验的相对偏差均在10%以上)。图3示出RUN72实验功率与温度模拟结果，二者的程序计算值与RUN72实验值均符合较好。

图2 TRACY装置结构[12]Fig.2 TRACY device structure[12]

表1 RUN72实验模拟参数Table 1 RUN72 experiment simulation parameter

表2 RUN72实验计算结果及相对偏差Table 2 RUN72 experiment calculation result and relative deviation

图3 RUN72实验功率与温度模拟结果Fig.3 Power and temperature simulation results of RUN72 experiment

4 瞬态分析程序在事故反演上的应用

核燃料瞬态分析程序的一个重要作用在于对已发生的事故进行反演，了解事故发生的过程，从而对事故的预防、屏蔽设计等提供指导。本文利用改进后的GETAC-S对日本JCO临界事故进行了反演[13]。JCO事故模拟参数列于表3[14-15]。

表3 JCO事故模拟参数Table 3 JCO accident simulation parameter

JCO事故的裂变功率震荡持续时间很长，由于事故发生后工厂采取了多种干预措施，对事故全过程的模拟非常困难。临界事故发生后人员受到的辐照程度主要由临界事故的第一裂变功率峰产额决定，因此本文主要对第一裂变功率峰裂变次数进行了模拟验证，其计算值和估计值分别为4.54×1016和5×1016，二者符合较好。

图4示出程序计算得到的功率和温度。由图4a可见，程序模拟结果能较好地体现临界事故的发展过程，即功率在短时间内上升到一较高峰值后迅速下降，最终趋向于一稳定的较低值。在这一过程中，由于辐照裂解气体逸出水面等因素的影响，在第一个功率峰值之后还可能会出现几个较低的峰值。

图4 JCO事故功率与温度的变化Fig.4 Power and temperature variations of JCO accident

5 结论

本文针对传统的基于点堆模型的溶液瞬态分析程序难以适应复杂几何及复杂物料情况的事实，基于蒙特卡罗均匀化理论与有限体积方法建立了适用于瞬态分析问题的三维扩散时空动力学模型，然后将该模型与非稳态传热模型、辐照裂解气泡模型相耦合，对GETAC-S程序进行升级。改进后的GETAC-S从截面生成端到时空动力学的求解，均实现了对复杂几何、材料及能谱条件下模型的适配，从而使程序可针对任意条件下的溶液系统进行瞬态特性分析。采用日本TRACY瞬态实验装置的临界事故模拟实验数据对改进后的GETAC-S的可靠性进行了验证，结果在10%内符合。针对日本JCO事故进行了事故进程反演，取得了较好结果。验证结果表明改进后的GETAC-S具备对复杂条件下溶液系统的临界事故后果进行评价与反演的能力，可为今后溶液核系统的设计起到一定的指导作用。