多物理耦合计算中动态输运问题高效蒙特卡罗模拟方法

在定量讨论各种核系统(如反应堆等)的动态行为时, 中子的时-空及能量分布是很关键的. 中子与原子核反应导致核能释放, 从而导致介质温度、密度和运动速度等物理量发生变化; 反过来, 因为中子与核的反应依赖于二者之间的相对运动能量, 从而介质中原子核热运动和流体运动将影响中子的输运过程. 基于物理上的考虑, 可以用Maxwell分布来描述核的热运动. 在对核的密度分布函数作为权重求平均以后, 中子输运方程中的系数将只依赖于中子与流体运动之间的相对速度. 基于相空间中子数守恒原理, 可得考虑流体宏观运动的中子输运方程[1]:

$ \begin{split} &\frac{{\partial N}}{{\partial t}} + {\boldsymbol{v}} \cdot \nabla N + \frac{{{\rm{d}}{\boldsymbol{V}}}}{{{\rm{d}}t}} \cdot {\nabla _V}N + V\varSigma ({\boldsymbol{r}},V)N({\boldsymbol{r}},{\boldsymbol{V}},t) \\ =\;& \int {{\varSigma '}{f'}{V'}{N'}{\rm{d}}{{\boldsymbol{V}}'}} + q({\boldsymbol{r}},{\boldsymbol{V}},t) {, } \\[-15pt]\end{split}$

其中$ {\boldsymbol{V}} = {\boldsymbol{v}} - {\boldsymbol{u}} $是相对速度, v是中子速度, u是流体速度, $ N({\boldsymbol{r}}, {\boldsymbol{V}}, t){\rm{d}}{\boldsymbol{r}}{\rm{d}}{\boldsymbol{V}} $是t时刻r处体积元${\rm{ d}}{\boldsymbol{r}} $内相对速度V处$ {\rm{d}}{\boldsymbol{V}} $范围内的中子可几数, ${\varSigma '}{f'}{V'}{N'} = $$ \varSigma ({\boldsymbol{r}}, {V'})f({\boldsymbol{r}}, {{\boldsymbol{V}}'} \to {\boldsymbol{V}}){V'}N({\boldsymbol{r}}, {{{V}}'}, t)$, $\varSigma({\boldsymbol{r}}, V)$是宏观总截面, $\varSigma({\boldsymbol{r}}, {V'})f({\boldsymbol{r}}, {{\boldsymbol{V}}'} \to {\boldsymbol{V}}) = \displaystyle\sum\nolimits_i \displaystyle\sum\nolimits_x \varSigma _x^i({\boldsymbol{r}}, {V'})\times $$ f_x^i({{\boldsymbol{V}}'} \to {\boldsymbol{V}})$, $\varSigma _x^i({\boldsymbol{r}}, {V'})$是一个中子在单位距离内与i类核经受某一由x标志的特定反应的概率, $ f_x^i({\boldsymbol{r}}, {{\boldsymbol{V}}'} \to {\boldsymbol{V}}) $表示r处速度为$ {{\boldsymbol{V}}'} $的中子与i类核发生x类反应产生速度为V的中子的概率, $ q({\boldsymbol{r}}, {\boldsymbol{V}}, t){\rm{d}}{\boldsymbol{r}}{\rm{d}}{\boldsymbol{V}} $表示r处$ {\rm{d}}{\boldsymbol{r}} $体积元内产生的速度为V处$ {\rm{d}}{\boldsymbol{V}} $范围内的中子数. 严格来说, 由于宏观截面依赖于中子通量, 所以上述方程中宏观截面是随时间变化的. 所以, 输运方程其实是一个非线性的微分-积分方程. 由于流体运动的时间尺度较中子运动的时间尺度长, 所以可选一个时间间隔, 其相对于中子输运相当长, 而对于流体运动来说相当短, 从而在该间隔内可以忽略核数密度改变的影响. 在时间间隔末通过解线性输运方程得到中子角通量以后, 就可以利用所得中子角通量、流体力学方程组和燃耗方程组求解新的物质密度和核数密度. 这样就把非线性输运方程转化为多时间步的线性输运问题, 形成了多物理耦合计算(本文仅考虑流体力学、输运和燃耗的多物理耦合计算).

随着科学技术的发展, 特别是并行计算能力的飞跃, 粒子输运问题的建模(包括几何建模、物理建模等)越来越逼近真实情况, 以前由于计算能力不足而不得不采取的各种近似逐渐被替换, 这就给蒙特卡罗粒子输运模拟带来了巨大的挑战. 以反应堆pin-by-pin模型蒙特卡罗粒子输运模拟领域著名的Kord-Smith挑战[2]为例, 由于其高达百万以上栅元, 几十亿数量级的计数规模的设置, 使得蒙特卡罗模拟面临存不下、算不快和算不准的难题. 对于多物理耦合计算中的动态蒙特卡罗粒子输运模拟而言, 其面临挑战的难度相较Kord-Smith挑战而言只高不低. 这主要是由动态输运问题高达千万以上网格的几何建模、由于极端反应条件造成的数十万以上的计算步、由于输运与燃耗耦合而必须求解的大量网格计数所造成的, 而且动态输运问题由于时间累积效应, 定态输运蒙特卡罗模拟所常用的许多技巧并不能直接应用, 需要根据动态问题的特点设计新的算法.

在公开发表的文献中, 反应堆大型pin-by-pin模型的高效蒙特卡罗模拟方法呈现出蓬勃发展的趋势. 区域分解[3,4]和数据分解[5,6]方法的提出主要是为了解决存不下的问题; 为了解决算不快的问题, 提出了若干高效样本数并行算法; 为了解决全局计数算不准的问题, 除了加大样本并高效并行以外, 还提出了UFS (uniform fission site method)[7, 8]和UTD (uniform tally density method)[9,10]等算法(这些算法都属于偏倚算法). 表1列出了常见的大型反应堆pin-by-pin模型的蒙特卡罗模拟方法对于存不下、算不快和算不准这3个难点各自的优缺点分析. 其中“√”表示左边的算法对上面难点的解决有促进作用, “ × ”表示表示左边的算法对上面难点的解决有阻碍作用, “—”表示没有影响. 框内的语句表示的是促进或阻碍作用最本质的原因. 这些研究对于动态输运问题的高效蒙特卡罗模拟也具有重要的借鉴作用.

表 1  常见算法优缺点比较

Table 1.  Relative merits of common algorithms.

从表1可以看出, 各种算法在解决一个问题的同时常引入新的问题. 对于本文关注的动态输运问题, 由于其不同的特点, 上述算法可以借鉴但不能照搬. 鉴于动态输运蒙特卡罗模拟极大的计算量, 设计高效算法的一个基本思路就是在基本不改变大样本直接模拟结果的基础上尽可能减少计算时间. 为此, 本文开展了针对输运燃耗耦合计算的计数规约优化算法和针对样本并行的样本数自适应算法研究. 这些研究是在北京应用物理与计算数学研究所自主开发的极端条件下三维多物理多介质问题数值模拟软件JUPITER (jointed numerical simulation software for multi-physics and multi-materials problems under extreme conditions)上进行的, 此平台对模型采取三维非结构网格建模, 在此基础上进行多物理耦合计算. 其中的蒙特卡罗粒子输运模拟模块(以下简称为JUPITER-MC程序)以区域分解加样本并行为基本架构, 可支持数千万三维非结构多面体网格、数百亿样本的中子输运及燃耗耦合计算, 具备数万处理器核的高效并行计算能力.

本文第2节介绍了一种基于输运燃耗耦合计算模式的计数规约新算法; 第3节基于蒙特卡罗模拟分批计算提出了一种样本数自适应算法; 第4节的数值结果表明两种算法都可以在基本不改变计算结果的前提下大幅减少计算时间, 从而提高了效率; 最后是总结.