多物理一体化临界安全程序INSL-UniFoam的研究

叶思远; 李泓龙; 李玥航; 柴翔; 刘晓晶; 何东豪

doi:10.11884/HPLPB202537.240369

摘要: 快中子脉冲堆（FBR）是临界安全分析研究的重要对象。它们具有几何形状不规则、强瞬态过程、多物理紧密耦合以及复杂的反馈特性等特点。为了精确模拟并分析快中子脉冲堆在瞬发超临界过程中各物理场的变化情况，开发了一门基于OpenFOAM的多物理核临界安全分析程序，名为INSL-UniFoam。该程序集成了离散纵坐标中子输运求解器、传热和应力-应变求解器，能够模拟快中子脉冲堆的瞬态超临界脉冲过程。程序在Godiva-I基准题中进行了验证，对稳态条件下的多个物理参数进行了敏感性分析。同时，程序还对Godiva-I的瞬态脉冲场景进行了计算并与实验结果进行了比对。结果表明，程序在中子学计算方面具有较高的精度，能精确反映脉冲堆的功率、中子通量的分布情况。同时瞬态耦合计算所得的脉冲功率曲线、峰值功率、裂变产额等方面与参考解符合良好，能够较好地反映脉冲过程并且能够完整地输出脉冲过程的功率、温度、应力应变在内多个物理场的分布情况并与实验结果较好地匹配。

关键词:

Abstract: Fast Burst Reactors (FBRs) are important subjects for criticality safety analysis research. They are characterized by irregular geometry, strong transient processes, tight multi-physics coupling, and complex feedback characteristics. This paper introduces an OpenFOAM based multi-physics nuclear criticality safety analysis code named INSL-UniFoam. It integrates discrete ordinate neutron transport solver, heat transfer and stress-strain solvers to detailly model the prompt super-critical burst pulse of FBRs. The UniFoam is first verified in the Godiva-I benchmark under both the steady-state condition and several transient scenarios. The results demonstrate that the program aligns well with the reference solution in terms of K_eff calculation, peak power, and fission yield. Furthermore, it is capable of comprehensively outputting the distributions of power, temperature, and stress-strain throughout the pulse process.

Key words:

快中子脉冲堆（Fast Burst Reactors, FBRs）是一类具有热膨胀自熄灭机制的反应堆^[1]，在功率脉冲过程中通过燃料热膨胀带来的负反馈效应使功率回落进而完成脉冲过程。自20纪50代以来，以Godiva-I为代表的快中子脉冲堆为对象开展了大量精细临界安全实验^[2]。快中子脉冲堆具有不规则几何、强瞬态过程、多物理耦合紧密、反馈特性复杂等特征^[3]，是核临界安全分析实验研究的重要对象。

针对Godiva-I一类的金属固体材料反应堆的瞬发超临界过程的分析涉及到中子学模拟、传热分析、应力应变模拟等多学科交叉集成，对程序的复杂几何处理能力以及耦合计算能力具有较高要求。Aufiero^[4]等人将蒙特卡罗学程序Serpent与CFD平台OpenFOAM相耦合，对Godiva-I的瞬发超临界过程进行了模拟，结果与实验数据符合良好。Wilson^[5]等人利用弹性力学方程以及中子学点堆方程，对脉冲堆进行了耦合计算。高辉^[6]等人利用蒙特卡罗程序与有限元程序ANSYS程序获得了Godiva-I的多物理耦合关系，并将耦合关系带入点堆方程从而获得了脉冲过程的波形。Fiorina^[7]等人基于有限体积法开发了多物理耦合程序，并对Godiva-I的脉冲过程进行了计算。郭树伟^[8]等人基于点堆方程、蒙卡程序和有限元程序，对Godiva-I脉冲堆开展了核热力耦合计算分析研究，并获得了表面位移等变化情况。Wang^[9]等人在MOOSE平台上基于蒙特卡罗输运计算以及点堆方程对Godiva-I脉冲堆的温度反馈系数进了计算，并在此基础上进行了耦合计算。

INSL-UniFoam是由上海交通大学先进核能系统实验室（INSL）开发的多物理临界安全分析程序。相较于传统耦合程序，INSL-UniFoam能够进行一体化的耦合计算，无需针对不同的温度、膨胀量进行多次蒙特卡罗计算，以获得反应性反馈情况，或在模拟过程中采用蒙特卡罗准静态方法并进行大量的蒙卡计算^[10]。同时，由于程序基于OpenFOAM平台，具有强大的几何处理能力，能够适应于不同的几何结构。在先前的研究中，INSL-UniFoam已经基于中子学扩散方程对Godiva-I基准题的稳态以及瞬态过程进行了模拟计算，其瞬态结果能够模拟出Godiva-I的脉冲波形。但由于扩散方程本身存在精度限制，在具有强各向异性散射特征的快中子脉冲堆中误差较大。本文建立了以离散纵标法中子输运方程为基础的多物理求解器，对Godiva-I基准题的稳态以及瞬态工况进行了精确的模拟。

1. 多物理方程以及求解器结构

1.1 中子输运方程

根据中子守恒关系，连续能量的中子输运方程可以写为

$\dfrac{1}{{v(E)}}\dfrac{{\partial \psi (E,{{\boldsymbol{\varOmega}} },{\boldsymbol{r}},t)}}{{\partial t}} + {{\boldsymbol{\varOmega}} } \cdot \nabla \psi (E,{{\boldsymbol{\varOmega}} },{\boldsymbol{r}},t) + {{{\mathit{\Sigma}} }_{\text{t}}}(E,{\boldsymbol{r}},t)\psi (E,{{\boldsymbol{\varOmega}} },{\boldsymbol{r}},t) = S(E,{{\boldsymbol{\varOmega}} },{\boldsymbol{r}},t)$

(1)

式中：$ v(E) $表示能量E的中子速度；$ \psi (E,{{\boldsymbol{\varOmega}} },{\boldsymbol{r}},t) $表示能量E、方向角$ {\boldsymbol{\varOmega }} $、位置r以及时间t下所对应的角通量大小；$ {{{\mathit{\Sigma }}}_{{\mathrm{t}}}} $代表当前位置材料的总截面；S代表总的中子源项，包括有散射、裂变以及外部的源项。

在实际应用当中需要对中子输运方程进行近似处理。对于能量，采用分群近似的方法将能量划分为多个能群从而进行离散化处理。对于空间坐标，本文程序采用有限体积法进行离散处理。方向角则采用离散纵标法将方向角划分为多组离散的方向，其中每个方向对应有不同的权重。

采用离散纵标法后，中子输运方程可以被表示为

$\dfrac{1}{{{v^g}}}\dfrac{{\partial \psi _\alpha ^g({\boldsymbol{r}},t)}}{{\partial t}} + {{{\boldsymbol{\varOmega}} }_\alpha } \cdot \nabla \psi _\alpha ^g({\boldsymbol{r}},t) + {{\mathit{\Sigma}} }_{\text{t}}^g({\boldsymbol{r}},t)\psi _\alpha ^g({\boldsymbol{r}},t) = Q_{\mathrm{s}}^g({\boldsymbol{r}},t,{{{\boldsymbol{\varOmega}} }_\alpha }) + \dfrac{{Q_d^g({\boldsymbol{r}},t)}}{{4\pi }} + \dfrac{{(1 - \beta ){\chi ^g}}}{{4\pi }}Q_{\text{f}}^g({\boldsymbol{r}},t) + \dfrac{{Q_{{\text{ex}}}^g({\boldsymbol{r}},t)}}{{4\pi }}$

(2)

式中：$ {v^g} $表示第g群的中子速度；$ \psi _\alpha ^g({\boldsymbol{r}},t) $表示第g群、在位置r处、在角度$\alpha $处、时间t时的中子通量密度；上标g代表对应的中子能群的编号；$ {\chi ^g} $为第g能群的中子裂变能谱；$ \; \beta $为总的缓发中子先驱核份额。缓发中子源项$ Q_d^g({\boldsymbol{r}},t) $以及裂变中子源项$ Q_{\text{f}}^g({\boldsymbol{r}},t) $可以表示为

$Q_d^g({\boldsymbol{r}},t) = \displaystyle\sum\limits_d^{{G_d}} {\chi _d^g} {\lambda _d}{C_d}({\boldsymbol{r}},t)$

(3)

$Q_{\text{f}}^g({{\boldsymbol{r}}},t) = \displaystyle\sum\limits_{g' = 1}^G v \Sigma _f^{g'}\psi _\alpha ^{g'}({{\boldsymbol{r}}},t)$

(4)

式中：$ {G_d} $表示缓发中子群的总数；G表示能群的总数；v表示每次裂变释放的平均中子数；$ {C_d} $表示第d组缓发中子先驱核的密度。

对应的缓发中子先驱核守恒方程则如下

$\dfrac{{\partial {C_d}({\boldsymbol{r}},t)}}{{\partial t}} = {\beta _d}Q_{\text{f}}^{\text{g}}({\boldsymbol{r}},t) - {\lambda _d}{C_d}({\boldsymbol{r}},t)$

(5)

式中：$ {C_d} $为第d组的缓发中子先驱核浓度；$ \;{ \beta _d} $和$ {\lambda _d} $分别为第d组的缓发中子先驱核份额和第d组的缓发中子先驱核衰变常数。

1.2 各向异性散射源的计算

中子输运方程的源项包括有散射源项、裂变源项以及外部的中子源项，其中散射源项在现实的反应堆中通常是各向异性的^[11]。在以快中子脉冲堆为代表的具有强泄露的堆型中，各向异性散射源项带来的影响往往不可忽视。因此为了使离散纵标法在快中子脉冲堆问题中得到更好的应用，需要引入各向异性散射源以提高计算精度。

在中子输运方程中，散射源项可以被写为

${Q_{\mathrm{s}}}({\boldsymbol{r}},E,{{\boldsymbol{\varOmega}} }) = \displaystyle\int_0^\infty {\displaystyle\int_{4\pi } {{{{\mathit{\Sigma}} }_{\mathrm{s}}}} \left( {{\boldsymbol{r}},E'} \right)f\left( {{\boldsymbol{r}},E' \to E,{{\boldsymbol{\varOmega}} '} \to {{\boldsymbol{\varOmega}} }} \right)\psi \left( {{\boldsymbol{r}},E',{{\boldsymbol{\varOmega}} '}} \right){\rm{d}} {{\boldsymbol{\varOmega}} '}} {\rm{d}} E'$

(6)

式中：$ f\left( {{\boldsymbol{r}},E' \to E,{{\boldsymbol{\varOmega}} '} \to {{\boldsymbol{\varOmega}} }} \right) $为散射函数，表示能量为$ E' $、方向为$ {{\boldsymbol{\varOmega}} '} $的中子经过散射后能量变为E、方向变为$ {{\boldsymbol{\varOmega}} } $的概率。

若令$ \;{ \mu _0} = \cos \left( {{{\boldsymbol{\varOmega}} '},{{\boldsymbol{\varOmega}} }} \right) = \cos {\theta _0} $，同时将能量进行离散化，使方程中的能量按能群进行划分，则有

$Q_{\mathrm{s}}^g({\boldsymbol{r}},{{\boldsymbol{\varOmega}} }) = \dfrac{1}{{2\pi }}\displaystyle\sum\limits_{g' = 1}^G {\displaystyle\int_{4\pi } {{{\mathit{\Sigma}} }_{\mathrm{s}}^{g' \to g}} \left( {{\boldsymbol{r}},{\mu _0}} \right){\psi _{g'}}\left( {{\boldsymbol{r}},{{\boldsymbol{\varOmega}} '}} \right){\rm{d}} {{\boldsymbol{\varOmega}} '}}$

(7)

在有限体积法的一个体积元中，可以认为截面等信息是均匀分布的，即$ {{\mathit{\Sigma}} }_{\mathrm{s}}^{g' \to g}\left( {{\boldsymbol{r}},{\mu _0}} \right) = {{\mathit{\Sigma}} }_{\mathrm{s}}^{g' \to g}\left( {{\mu _0}} \right) $。进而对于方向角$\varphi $，散射源项$Q_{\mathrm{s}}^g$可以被表示为

$Q_{\mathrm{s}}^g\left( {{\boldsymbol{r}},{\mu ^m},{\eta ^m},{\xi ^m}} \right) = \dfrac{1}{{2\pi }}\displaystyle\sum\limits_{g' = 1}^G {\displaystyle\int_{ - 1}^1 {\displaystyle\int_0^{2\pi } {{{\mathit{\Sigma}} }_{\mathrm{s}}^{g' \to g}\left( {{\mu _0}} \right){\psi _{g'}}\left( {{\boldsymbol{r}},\mu ',\eta ',\xi '} \right)} } {\rm{d}} \varphi '{\rm{d}} \mu '}$

(8)

式中：$\left( {{\mu ^m},{\eta ^m},{\xi ^m}} \right)$代表中子运动方向角与坐标轴的余弦。因此可以利用勒让德多项式将散射截面${{\mathit{\Sigma}} }_{\text{s}}^{g' \to g}\left( {{\mu _0}} \right)$展开为

$\Sigma _{\text{s}}^{g' \to g}\left( {{\mu _0}} \right) = \displaystyle\sum\limits_{l = 0}^L {\dfrac{{2l + 1}}{2}\Sigma _{{\mathrm{s}},l}^{g' \to g}{{\mathrm{P}}_l}\left( {{\mu _0}} \right)}$

(9)

式中：l表示各向异性散射的阶数；$ {{\rm{P}}_l}\left( {{\mu _0}} \right) $为勒让德多项式，可以被写为

${{\mathrm{P}}_l}\left( {{\mu _0}} \right) = {{\mathrm{P}}_l}(\mu ){{\mathrm{P}}_l}\left( {\mu '} \right) + 2\displaystyle\sum\limits_{k = 1}^l {\dfrac{{(l - k)!}}{{(l + k)!}}{\mathrm{P}}_l^k(\mu ){\mathrm{P}}_l^k\left( {\mu '} \right)\cos \left[ {k\left( {\varphi - \varphi '} \right)} \right]}$

(10)

通过结合以上方程，最终中子散射源项可以被写为

$Q_{{\text{s}},{\text{g}}}^m = \displaystyle\sum\limits_{l = 0}^L {(2l + 1)} \displaystyle\sum\limits_{g' = 1}^G {\Sigma _{s,l}^{g' \to g}{{\mathrm{P}}_l}\left( {{\mu ^m}} \right)\psi _{g'}^l} + 2\displaystyle\sum\limits_{l = 0}^L {(2l + 1)} \displaystyle\sum\limits_{g' = 1}^G {\Sigma _{s,l}^{g' \to g}} \displaystyle\sum\limits_{k = 1}^l {\dfrac{{(l - k)!}}{{(l + k)!}}{\mathrm{P}}_l^k\left( {{\mu ^m}} \right)\left( {\cos \left( {k{\varphi ^m}} \right)\psi _{g'}^{l,k,1} + \sin \left( {k{\varphi ^m}} \right)\psi _{g'}^{l,k,2}} \right)}$

(11)

$\psi _{g'}^l = \displaystyle\sum\limits_{m = 1}^M {{\omega ^m}} {{\mathrm{P}}_l}\left( {{\mu ^m}} \right)\bar \psi _g^m$

(12)

$\psi _{g'}^{l,k,1} = \displaystyle\sum\limits_{m = 1}^M {{\omega ^m}} {\mathrm{P}}_l^k\left( {{\mu ^m}} \right)\cos \left( {k{\varphi ^m}} \right)\bar \psi _g^m$

(13)

$\psi _{g'}^{l,k,2} = \displaystyle\sum\limits_{m = 1}^M {{\omega ^m}} {\mathrm{P}}_l^k\left( {{\mu ^m}} \right)\sin \left( {k{\varphi ^m}} \right)\bar \psi _g^m$

(14)

式中：${\omega ^m}$表示求积组中方向m的权重系数；${{\mathrm{P}}_l}\left( {{\mu ^m}} \right)$表示$ \left( {n,m} \right) $阶的勒让德多项式。

1.3 传热以及应力计算

在模拟快中子脉冲堆的瞬态脉冲过程时，需要计算短时间内由堆芯功率快速变化带来的热应力冲击。带有内热源的传热微分方程为

$\rho c\dfrac{{\partial T}}{{\partial t}} = \nabla \cdot (k\nabla T) + Q$

(15)

式中：$ \rho $是物质的密度；c是物质的比热容；T是温度；t是时间；k是物质的热导率；$ \nabla $是梯度算子；Q是热源。通过将中子学计算求得的热功率分布代入方程中，即可求得小型快中子脉冲堆在某一时刻的温度分布情况。

在求得温度分布情况之后，使用基于位移的瞬态分离有限体积方程进行热应力计算，即

$\dfrac{{{\partial ^2}{\boldsymbol{D}}}}{{\partial {t^2}}} = \nabla \cdot \left( {\mu \left( {\nabla {\boldsymbol{D}} + \nabla {{\boldsymbol{D}}^{\mathrm{T}}}} \right) + \lambda (\nabla \cdot {\boldsymbol{D}})I} \right) + \nabla \dfrac{E}{{1 - 2v}}T$

(16)

式中：$ {{\boldsymbol{D}}} $代表位移量；E代表杨氏模量；$ \; \mu $与$ \lambda $分别为拉梅系数的第一和第二个参数，这两个参数的计算公式分别为

$\; \mu = \dfrac{E}{{2(1 + v)}}$

(17)

$\lambda = \dfrac{{vE}}{{(1 + v)(1 - 2v)}}$

(18)

进而可以求得应变张量$ {{\boldsymbol{U}}} $为

${{\boldsymbol{U}}} = \dfrac{1}{2}\left( {\nabla {{\boldsymbol{D}}} + \nabla {{{\boldsymbol{D}}}^{\mathrm{T}}}} \right)$

(19)

求解以上方程，能够获得几何网格每个体积元所对应的应变张量，进而能够获得每个体积元的体积变化量为

$\varepsilon = |{{\mathrm{grad}}} {{\boldsymbol{D}}} - {\boldsymbol{E}}| - 1$

(20)

在计算由于堆芯温度上升以及膨胀引起的中子学反应性反馈时，仅考虑由热膨胀效应引起的反应性变化。因而材料对中子的宏观截面随着体积膨胀发生改变可表示为

${{{\mathit{\Sigma}} }^{t + 1}} = \dfrac{1}{{1 + \varepsilon }}{{{\mathit{\Sigma}} }^t}$

(21)

式中：$ {\Sigma ^t} $表示时间步t对应使用的中子各项截面。

1.4 边界条件

程序对脉冲堆的各个物理场建立方程并进行计算前，需要指定物理场的初始边界条件。本求解器的边界条件设定遵循OpenFOAM的格式，通过设定算例在开始时间点下的各个物理场的初始数值以及边界信息来确立模拟过程中使用的边界条件。

对于所有的中子通量场，由于模拟过程中不考虑反应堆外的环境对于中子的反射，程序在模拟中将中子的外边界设定为泄露边界条件，即所有的中子角通量能够自由穿出外边界，但边界外部没有中子进入反应堆内。

对于固体材料外部自由膨胀的情况，描述材料位移量的位移场$ {{\boldsymbol{D}}} $的边界条件设置为零牵引力边界条件，即初始条件下边界处无位移值，同时边界上的牵引力以及所受的压力均为0。每个体积元的体积变化量对应的物理场$ \varepsilon $设定为初始值为0的零体积变化条件，即边界处的体积元对应的体积变化量为0。考虑到所模拟的脉冲堆的脉冲过程时间短，同时外部空气导热能力差，在模拟中将温度场设定为零梯度（绝热）边界条件。

1.5 程序结构

本文的程序INSL-UniFoam基于OpenFOAM （Open Field Operation And Manipulation）^[12]计算平台建立，基于有限体积法进行求解。基于C++的OpenFOAM平台是一款完全开源并已被广泛应用于传热、应力以及流体力学等领域求解的平台。通过采用有限体积法（Finite Volume Method，FVM），将要计算的包括中子角通量、温度场、应力场等物理实体切分为网格单元来描述，并使用发散定理，将所有包含发散项的偏微分方程中的体积积分转换为表面积分，基于物理量的守恒性，将每个网格的项加总成每个体积表面的通量从而进行计算^[13]。

INSL-UniFoam程序采用内耦合的方式，利用公共变量共享数据，在各个物理场对应的偏微分方程之间建立联系。程序在求解器内部的每个时间步内分别求解中子输运方程、热传导方程以及热应变方程，从而得出每个体积元的体积变化量，并根据整个反应堆的体积变化量对应的标量场生成截面变化的系数分布标量场，进而在后续的中子输运方程的计算中引入截面变化系数，实现核-热-力的耦合计算。

相较于传统使用点堆方程的耦合程序，INSL-UniFoam无需在进行耦合计算前针对不同温度以及膨胀情况进行多次蒙卡计算以获得膨胀产生的反应性反馈。这样极大提高了针对不同物性参数以及反应堆构型进行临界安全分析验证的效率。

整个程序的计算流程如图1所示。

图 1 INSL-UniFoam程序计算流程图

Figure 1. Calculation flowchart of INSL-UniFoam

下载: 全尺寸图片幻灯片

2. 数值验证

2.1 敏感性分析

INSL-UniFoam的计算验证针对Godiva-I基准题进行。Godiva-I脉冲堆为洛斯阿拉莫斯实验室（LANL）于1953年设计建造的块中子脉冲堆^[14]。Godiva-I的堆芯是一个由三块金属拼接而成的半径8.740 7 cm，密度18.739 8 g/cm³，富集度93.8%的铀金属球。为了验证计算精度并确定瞬态计算时的参数，分别对瞬态过程中使用的包括散射矩阵阶数、网格模型、能群结构以及方向角数量等参数进行了敏感性分析。

首先对使用散射矩阵的阶数对计算精度的影响进行了分析。在通过blockMesh软件对Godiva-I进行建模以及网格划分，网格数量27 001、能群数量33群、采用S8全对称求积组（每个卦限10个离散角）的角度离散情况下，分别对散射矩阵不同阶数进行了测试，分别对比了其与OpenMC参考解（有效增殖因子k_eff=0.994 56）的误差。测试结果如表1所示。

表 1 不同散射阶数下k_eff计算结果

Table 1. Calculated k_eff of Godiva-I in different orders of scattering matrices

scattering matrix	k_eff	error/10⁻⁵
OpenMC	0.99456	−
P0	1.09455	10003
P1	0.97971	1487
P3	0.98583	869

下载: 导出CSV

| 显示表格

从表1中可见，在对Godiva-I这种具有强泄露的模型的模拟中，采用高阶散射能够大大降低误差。在后续计算验证中，均采用P3散射矩阵进行各向异性散射源的计算。

由于OpenFOAM平台具有支持多种几何模型文件的能力，能够测试不同的网格生成方式以及网格数量对计算结果的影响。分别利用OpenFOAM自带的blockMesh程序以及开源网格划分软件Gmsh^[15]构建六面体网格以及四面体网格进行测试计算，两个软件生成的网格结构如图2所示。测试在能群数量8群、采用S4全对称求积组、使用P3阶散射矩阵的情况下进行。计算结果如表2所示。可以看到在相似数量的网格下，六面体网格在Godiva-I中的表现优于四面体网格。这是由于六面体网格总体网格质量相对较高，而四面体网格相对较为杂乱，使得四面体网格需要更多的网格数量才能达到同样的精度。因此，本文后续工作将采用六面体网格（blockMesh）进行Godiva-I建模，模型总共56 000个网格。

图 2 不同软件生成的Godiva-I网格模型

Figure 2. Meshing model of Godiva-I by different tools

下载: 全尺寸图片幻灯片

表 2 不同网格下的k_eff计算结果

Table 2. Calculated k_eff in different meshes

meshing tool	cell number	k_eff	error/10⁻⁵
OpenMC	−	0.99456	−
blockMesh	56000	0.99272	184
Gmsh	60187	0.99971	515
Gmsh	384803	0.99529	73

下载: 导出CSV

| 显示表格

由于Godiva-I 是一种无慢化剂的快谱装置，需要较为精细的能群划分。在均采用全队成求积组S8、P3阶散射矩阵的情况下，对不同的能群结构分别进行了计算。所使用的能群数分别为1群、8群以及33群。能群结构如表3所示：

表 3 不同能群数量对应的能量网格划分

Table 3. Energy grids of different energy group structures

33 group energy grid/MeV				8 group energy grid/MeV
1.964033E+01	1.831564E-01	2.034684E-03	2.260329E-05	1.00E+01
1.000000E+01	1.110900E-01	1.234098E-03	1.370959E-05	8.21E-01
6.065307E+00	6.737947E-02	7.485183E-04	8.315287E-06	5.53E-03
3.678794E+00	4.086771E-02	4.539993E-04	4.000000E-06	4.00E-06
2.231302E+00	2.478752E-02	3.043248E-04	5.400000E-07	6.25E-07
1.353353E+00	1.503439E-02	1.486254E-04	1.000000E-07	2.80E-07
8.208500E-01	9.118820E-03	9.166088E-05	1.000010E-11	1.40E-07
4.978707E-01	5.530844E-03	6.790405E-05		5.80E-08
3.019738E-01	3.354626E-03	4.016900E-05		1.00E-11

下载: 导出CSV

| 显示表格

从表4中可以看到，使用8个能群时k_eff计算结果的误差为174×10⁻⁵。将能群数量增加到33个后，与OpenMC的k_eff误差减少到79×10⁻⁵。然而，在求解瞬态中子输运方程时，33个能群显著增加了计算成本，因此我们选择了能群结构如表3所示的8群划分的截面用于后续计算。

表 4 不同能群数量下k_eff计算结果

Table 4. Calculated k_eff in different energy groups

energy group number	k_eff	Error/10⁻⁵
OpenMC	0.99456	−
1	0.95974	3482
8	0.99282	174
33	0.99377	79

下载: 导出CSV

| 显示表格

为了正确模拟快中子脉冲堆中的强各向异性散射，INSL-UniFoam允许用户自定义求积组。对于Godiva-I，由于其球形几何结构，选择使用全对称求积组。在Godiva-I采用56 000个几何网格单元、P3阶散射矩阵和总共8个能群的条件下，使用全对称S2、S4和S8求积组进行了角度离散化计算。图3显示了计算中使用的S4和S8全对称求积组在其中一个卦限上的示意图。计算结果如表5所示。

图 3 S4和S8全对称求积组示意图

Figure 3. S4 and S8 fully symmetric quadrature sets

下载: 全尺寸图片幻灯片

表 5 选取不同求积组的k_eff计算结果

Table 5. Calculated k_eff in different quadrature sets

quadrature set	k_eff	error/10⁻⁵
OpenMC	0.99456	−
S2	0.94014	5442
S4	0.99272	184
S8	0.99282	174

下载: 导出CSV

| 显示表格

从表中可以看出，当使用S2全对称求积组时，此时每个卦限只有一个方向角，k_eff计算结果与参考解相差5442×10⁻⁵。当使用S4全对称求积组后，每个卦限对应三个方向角，k_eff误差减少到184×10⁻⁵。此后继续提升方向角数量对于计算精度的改进相对有限。

最终，瞬态计算使用的网格由blockMesh生成，网格数量为56000个。能群采用8群，结构如表3所示。各向异性散射计算采用P3阶散射矩阵。在这种情况下，Godiva-I的稳态计算结果k_eff误差为184×10⁻⁵。

2.2 瞬态验证

瞬态过程中热-应力子求解器所使用到的部分物性参数如表6所示。

表 6 热应力求解使用的物性参数

Table 6. Thermo-mechanical parameters of Godiva-I

density/ (kg·m⁻³)	Poisson’s ratio/GPa	Young’s modulus	specific heat capacity/(J·kg⁻¹·K⁻¹)	thermal conductivity/ (W·m⁻¹·K⁻¹)	coefficient of linear expansion/K⁻¹
18740	0.23	208	117.7	27.5	1.39$ \times $10⁻⁵

下载: 导出CSV

| 显示表格

通过采用上文的各项参数，对Godiva-I进行了瞬发超临界工况的模拟。瞬发超临界起始条件下，反应堆周期为29.5 µs，计算过程中时间步步长为10⁻⁶ s。将计算获得的功率变化曲线与实验数据^[16]进行对比，结果如图4所示。可以看到，功率曲线在上升阶段能够较好地吻合实验数据，但在功率脉冲到达峰值后的下降阶段，UniFoam计算所得的功率数值低于实验数据。进一步对比UniFoam与实验数据的峰值功率以及半峰宽结果，对比情况如表7所示。程序计算得到的反应堆的峰值功率为7.397 3×10⁸ W，对应的峰值裂变率为2.36×10⁸ s⁻¹，而实验数据的峰值裂变率为2.62×10¹⁹ s⁻¹。UniFoam计算得到的半峰宽大小为8.4×10⁻⁵ s，而实验测得数据对应的半峰宽约为11.1×10⁻⁵ s。因此，UniFoam计算所得功率峰值为实验数据的90.1%，而半峰宽则为实验数据的75.7%。

图 4 29.5 μs初始反应堆周期下反应堆功率随时间变化情况

Figure 4. Total reactor power over time with initial reactor periods of 29.5 μs

下载: 全尺寸图片幻灯片

表 7 进一步对比时热应力求解使用的物性参数

Table 7. Additional thermo-mechanical parameters of Godiva-I

	peak power/W	full width at half maximum/s
INSL-UniFoam	7.397E+8	8.4E-5
experiment result	8.223E+8	11.1E-5

下载: 导出CSV

| 显示表格

造成此误差的主要原因和实验装置的外部环境有关。部分学者^[17]分析认为，在临界实验期间，Godiva-I被放置在建筑物内，这导致一些中子从墙壁反射回来并继续参与核反应。这种部分反射的边界条件显然与本文耦合程序采用的泄漏边界条件不一致，进而使得计算结果的峰值以及半峰宽数值偏低。同时，反应堆膨胀以及温度上升对材料中子学截面产生反馈的计算方式也会对计算结果造成影响。此外，本研究中Godiva-I模型未指定内部的力学边界也是模拟结果与实验值产生偏差的另一个原因。Godiva-I的实验装置的球形结构由三个金属块拼接而成，在脉冲过程中每个金属块接触面的力学性质不能等同于本研究中作为一个整体的固体金属内部的力学性质。

除了总功率曲线之外，还对平均温升以及平均体积膨胀率曲线进行了绘制。图5（a）为平均温度曲线，图5（b）为平均体积膨胀率曲线。可以看到，在脉冲过程的前期，由于中子通量较低导致功率较低，反应堆温度上升速度较慢。随后随着功率的迅速上升，且越靠近反应堆中心功率数值越大，大量热能累积在金属球内部，使得金属发生膨胀。体积的膨胀导致材料密度减小，进而对反应性产生了负反馈，使得脉冲功率下降并达到一个相对稳定的数值。随着功率的下降，温度以及体积膨胀率曲线也趋于平稳。

图 5 反应堆平均温度以及平均体积膨胀率随时间变化情况

Figure 5. Changes in average temperature and average volume expansion rate of the reactor over time

下载: 全尺寸图片幻灯片

图6（a）为功率分布情况，图6（b）为应力场的分布情况。可以看到经历了脉冲后，虽然功率集中在球体中心位置，但是中心热量的堆积引起的体积的膨胀在球体外部更为显著，这导致了金属球的外部产生了更大的位移，使大量应力集中在球体的外层区域。

图 6 功率最大点的功率场分布情况以及脉冲峰结束后的应力分布情况

Figure 6. Distribution of power field at the point of maximum power and the stress distribution after the pulse peak

下载: 全尺寸图片幻灯片

3. 结　论

本文基于OpenFOAM建立了临界安全分析计算程序INSL-UniFoam，对Godiva-I基准题进行了稳态以及瞬态求解，分别对比了不同散射阶数以及网格条件下的计算结果并将其与OpenMC参考解比较。稳态计算结果表明程序能够较好地对k_eff值进行计算，与参考解误差小于200×10⁻⁵。计算的精度相较于传统使用点堆方程以及扩散方程来模拟瞬发超临界过程的耦合程序有较大提升，能更真实地反应并输出脉冲堆内的功率、中子通量等物理场的分布情况。

同时程序能较好地对快Godiva-I基准题的瞬发超临界工况进行模拟，在峰值功率、裂变产额等方面与实验参考解符合较好，能够完整地输出整个脉冲过程中功率场、应力场等物理场的变化过程。同时对于温度场、膨胀量等安全分析中较为关注的信息，能够绘制出脉冲过程中相关物理场在整个反应堆中的分布情况。

程序计算结果表明，程序的精度较高，能够适应于具有强泄露的快中子脉冲堆的多群中子瞬态输运问题的求解以及多物理场耦合计算。同时，由于OpenFOAM平台具有强大几何处理能力，INSL-UniFoam能够对非结构化网格的复杂几何模型进行处理计算，为后续开展针对具有更复杂几何结构的快中子脉冲堆的计算建立了基础。

图 1 INSL-UniFoam程序计算流程图

Figure 1. Calculation flowchart of INSL-UniFoam

下载: 全尺寸图片幻灯片

图 2 不同软件生成的Godiva-I网格模型

Figure 2. Meshing model of Godiva-I by different tools

下载: 全尺寸图片幻灯片

图 3 S4和S8全对称求积组示意图

Figure 3. S4 and S8 fully symmetric quadrature sets

下载: 全尺寸图片幻灯片

图 4 29.5 μs初始反应堆周期下反应堆功率随时间变化情况

Figure 4. Total reactor power over time with initial reactor periods of 29.5 μs

下载: 全尺寸图片幻灯片

图 5 反应堆平均温度以及平均体积膨胀率随时间变化情况

Figure 5. Changes in average temperature and average volume expansion rate of the reactor over time

下载: 全尺寸图片幻灯片

图 6 功率最大点的功率场分布情况以及脉冲峰结束后的应力分布情况

Figure 6. Distribution of power field at the point of maximum power and the stress distribution after the pulse peak

下载: 全尺寸图片幻灯片

表 1 不同散射阶数下k_eff计算结果

Table 1. Calculated k_eff of Godiva-I in different orders of scattering matrices

scattering matrix	k_eff	error/10⁻⁵
OpenMC	0.99456	−
P0	1.09455	10003
P1	0.97971	1487
P3	0.98583	869

下载: 导出CSV

表 2 不同网格下的k_eff计算结果

Table 2. Calculated k_eff in different meshes

meshing tool	cell number	k_eff	error/10⁻⁵
OpenMC	−	0.99456	−
blockMesh	56000	0.99272	184
Gmsh	60187	0.99971	515
Gmsh	384803	0.99529	73

下载: 导出CSV

表 3 不同能群数量对应的能量网格划分

Table 3. Energy grids of different energy group structures

33 group energy grid/MeV				8 group energy grid/MeV
1.964033E+01	1.831564E-01	2.034684E-03	2.260329E-05	1.00E+01
1.000000E+01	1.110900E-01	1.234098E-03	1.370959E-05	8.21E-01
6.065307E+00	6.737947E-02	7.485183E-04	8.315287E-06	5.53E-03
3.678794E+00	4.086771E-02	4.539993E-04	4.000000E-06	4.00E-06
2.231302E+00	2.478752E-02	3.043248E-04	5.400000E-07	6.25E-07
1.353353E+00	1.503439E-02	1.486254E-04	1.000000E-07	2.80E-07
8.208500E-01	9.118820E-03	9.166088E-05	1.000010E-11	1.40E-07
4.978707E-01	5.530844E-03	6.790405E-05		5.80E-08
3.019738E-01	3.354626E-03	4.016900E-05		1.00E-11

下载: 导出CSV

表 4 不同能群数量下k_eff计算结果

Table 4. Calculated k_eff in different energy groups

energy group number	k_eff	Error/10⁻⁵
OpenMC	0.99456	−
1	0.95974	3482
8	0.99282	174
33	0.99377	79

下载: 导出CSV

表 5 选取不同求积组的k_eff计算结果

Table 5. Calculated k_eff in different quadrature sets

quadrature set	k_eff	error/10⁻⁵
OpenMC	0.99456	−
S2	0.94014	5442
S4	0.99272	184
S8	0.99282	174

下载: 导出CSV

表 6 热应力求解使用的物性参数

Table 6. Thermo-mechanical parameters of Godiva-I

density/ (kg·m⁻³)	Poisson’s ratio/GPa	Young’s modulus	specific heat capacity/(J·kg⁻¹·K⁻¹)	thermal conductivity/ (W·m⁻¹·K⁻¹)	coefficient of linear expansion/K⁻¹
18740	0.23	208	117.7	27.5	1.39$ \times $10⁻⁵

下载: 导出CSV

表 7 进一步对比时热应力求解使用的物性参数

Table 7. Additional thermo-mechanical parameters of Godiva-I

	peak power/W	full width at half maximum/s
INSL-UniFoam	7.397E+8	8.4E-5
experiment result	8.223E+8	11.1E-5

下载: 导出CSV

[1]	Kadioglu S Y, Knoll D A, De Oliveira C. Multiphysics analysis of spherical fast burst reactors[J]. Nuclear Science and Engineering, 2009, 163(2): 132-143. doi: 10.13182/NSE09-07
[2]	Briggs J B, Scott L, Nouri A. The international criticality safety benchmark evaluation project[J]. Nuclear Science and Engineering, 2003, 145(1): 1-10. doi: 10.13182/NSE03-14
[3]	Wimett T F. Los Alamos: Los Alamos Scientific Lab, 1965.
[4]	Aufiero M, Fiorina C, Laureau A, et al. Serpent–OpenFOAM coupling in transient mode: simulation of a Godiva prompt critical burst[C]//Proceedings of Joint International Conference on Mathematics and Computation (M&C), Supercomputing in Nuclear Applications (SNA) and the Monte Carlo (MC) Method. 2015.
[5]	Wilson S C, Biegalski S R, Coats R L. Computational modeling of coupled thermomechanical and neutron transport behavior in a Godiva-like nuclear assembly[J]. Nuclear science and engineering, 2007, 157(3): 344-353. doi: 10.13182/NSE06-28
[6]	高辉, 钟力晗, 梁文峰, 等. 基于反应性温度系数的金属型脉冲堆波形计算[J]. 原子能科学技术, 2017, 51(5):798-802 doi: 10.7538/yzk.2017.51.05.0798 Gao Hui, Zhong Lihan, Liang Wenfeng, et al. Waveform calculation of metal burst reactors based on reactivity temperature coefficient[J]. Atomic Energy Science and Technology, 2017, 51(5): 798-802 doi: 10.7538/yzk.2017.51.05.0798
[7]	Fiorina C, Aufiero M, Pelloni S, et al. A time-dependent solver for coupled neutron-transport thermal-mechanics calculations and simulation of a Godiva prompt-critical burst[C]//Proceedings of the 2014 22nd International Conference on Nuclear Engineering. 2014.
[8]	郭树伟, 陈珍平, 江新标, 等. 金属核燃料快中子脉冲堆核-热-力耦合计算方法研究[J]. 核动力工程, 2022, 43(4):31-37 Guo Shuwei, Chen Zhenping, Jiang Xinbiao, et al. Study on neutronic/thermal-mechanical coupling calculation method for fast-neutron pulse reactor with metallic nuclear fuel[J]. Nuclear Power Engineering, 2022, 43(4): 31-37
[9]	Wang Lipeng, Guo Shuwei, Hu Tianliang, et al. Transient simulation and parameter sensitivity analysis of Godiva experiment based on MOOSE platform[J]. Energies, 2023, 16(18): 6575. doi: 10.3390/en16186575
[10]	Blanco J A. Neutronic, thermohydraulic and thermomechanical coupling for the modeling of criticality accidents in nuclear systems[D]. Grenoble: Université Grenoble Alpes, 2020.
[11]	杨波. 离散纵标法求解含有各向异性散射的输运方程[D]. 绵阳: 中国工程物理研究院, 2005 Yang Bo. Solution of transport equations with anisotropic scattering using the discrete ordinates method[D]. Mianyang: China Academy of Engineering Physics, 2005
[12]	Jasak H, Jemcov A, Tuković Z. OpenFOAM: a C++ library for complex physics simulations[C]//Proceedings of the International Workshop on Coupled Methods in Numerical Dynamics. 2007.
[13]	Eymard R, Gallouët T, Herbin R. Finite volume methods[J]. Handbook of Numerical Analysis, 2000, 7: 713-1018.
[14]	Peterson R, Newby G. Lady Godiva: an unreflected uranium-235 critical assembly[R]. Los Alamos: Los Alamos National Lab, 1953.
[15]	Geuzaine C, Remacle J F. Gmsh: a 3-D finite element mesh generator with built-in pre-and post-processing facilities[J]. International Journal for Numerical Methods in Engineering, 2009, 79(11): 1309-1331. doi: 10.1002/nme.2579
[16]	Wimett T, Engle L, Graves G, et al. Time behavior of Godiva through prompt critical[R]. Los Alamos: Los Alamos Scientific Lab, 1956.
[17]	张驰, 周琦, 朱庆福, 等. 金属核燃料系统瞬态特性分析研究[J]. 原子能科学技术, 2016, 50(12):2170-2174 doi: 10.7538/yzk.2016.50.12.2170 Zhang Chi, Zhou Qi, Zhu Qingfu, et al. Transient characteristic analysis of nuclear metallic fuel system[J]. Atomic Energy Science and Technology, 2016, 50(12): 2170-2174 doi: 10.7538/yzk.2016.50.12.2170

姓名
邮箱
手机号码
标题
留言内容
验证码

留言板