«上一篇
文章快速检索     高级检索
下一篇»
  哈尔滨工程大学学报  2021, Vol. 42 Issue (12): 1762-1768  DOI: 10.11990/jheu.202010033
0

引用本文  

康乐, 郝琛. 精细化中子输运计算程序HNET动力学计算初步验证[J]. 哈尔滨工程大学学报, 2021, 42(12): 1762-1768. DOI: 10.11990/jheu.202010033.
KANG Le, HAO Chen. Preliminary verification of the kinetics calculation for the high-fidelity neutron transport program HNET[J]. Journal of Harbin Engineering University, 2021, 42(12): 1762-1768. DOI: 10.11990/jheu.202010033.

基金项目

国家重点研发计划项目(2018YFE0180900)

通信作者

郝琛,E-mail: haochen.heu@163.com

作者简介

康乐,男,博士研究生;
郝琛,男,副教授,博士生导师

文章历史

收稿日期:2020-10-21
网络出版日期:2021-09-23
精细化中子输运计算程序HNET动力学计算初步验证
康乐 , 郝琛     
哈尔滨工程大学 核安全与仿真技术国防重点学科实验室,黑龙江 哈尔滨 150001
摘要:针对精细化瞬态中子输运计算的高耗时问题,本文开展了瞬态中子输运计算方法研究及程序开发。基于二维/一维方法开发了三维全堆芯一步法精细化中子输运计算程序,着重开展了时空中子动力学的算法研究及程序模块开发。针对三维瞬态输运问题,将其转化为二维与一维的瞬态固定源问题,分别采用特征线法与节块展开法进行求解,同时,采用两级基于广义等价理论的粗网有限差分方法作为三维计算框架并加速全局的迭代收敛。最后,采用C5G7-TD非均匀瞬态基准题对该程序的中子动力学模块开展计算精度方面的验证,并与国际知名高保真中子学程序进行对比分析。结果表明: 自主开发的时空动力学模块计算精度高,能够满足三维高保真瞬态中子输运计算的需求。
关键词时空动力学    中子输运方程    二维/一维方法    特征线法    广义等价理论    粗网有限差分    两级加速    C5G7-TD    高保真    
Preliminary verification of the kinetics calculation for the high-fidelity neutron transport program HNET
KANG Le , HAO Chen     
Fundamental Science on Nuclear Safety and Simulation Technology Laboratory, Harbin Engineering University, Harbin 150001, China
Abstract: To solve the time consumption problem in high-fidelity transient neutron transport calculations, the development of a transient transport calculation method and program is presented in this work. Based on the 2D/1D method, a 3D whole-core one-step high-fidelity neutron transport calculation program HNET (high-fidelity neutron transport program) is developed. This study focuses on the calculation method and the program development of neutron kinetics. The 3D transient transport problem is translated into 2D and 1D transient fixed source problems, which are solved by the method of characteristics and nodal expansion method, respectively. Then, a two-level generalized equivalence theory-based coarse mesh finite difference method is adopted at the same time as a 3D calculation framework and accelerator for the convergence of global iteration. To verify the calculation accuracy of the transient calculation module of HNET, C5G7-TD heterogeneous transient benchmark problems are used. The obtained results are analyzed and compared with the well-known high-fidelity neutronics program MPACT. The numerical results indicate that the space-time kinetics module of HNET is accurate and can satisfy the requirements of 3D high-fidelity transient transport calculations.
Keywords: space-time kinetics    neutron transport equation    2D/1D method    method of characteristics    generalized equivalence theory    coarse mesh finite difference    two-level acceleration    C5G7-TD    high-fidelity    

随着高性能计算能力的快速提升,反应堆数值模拟技术在精细化方面的需求逐渐提高。然而,直接采用标准三维中子输运方法开展全堆芯计算仍不现实,主要是因为典型反应堆堆芯的实际未知数总数过大,对于稳态计算接近1015数量级,而对于瞬态计算来说则要更大。目前,一种实用且切实可行的方案是采用二维/一维方法[1-4],即通过定义轴向泄漏项和径向泄漏项,将全三维输运计算转为径向二维输运计算耦合轴向一维低阶输运计算。同时,精细化全三维瞬态输运计算依赖于有效的全局加速方法的应用。其中,粗网有限差分方法(coarse mesh finite difference, CMFD)是一种流行的全局加速方法,特别针对二维/一维聚合算法,CMFD还为径向二维输运计算和轴向一维低阶输运计算提供了一种天然聚合计算框架。但是,传统CMFD有收敛慢的缺陷,同时还存在收敛发散及耦合系数为负的潜在缺陷[5]。为了有效降低多群CMFD带来的计算负担,在研发HNET程序(high-fidelity neutron trasport program for 3D nuclear reactor core,HNET)过程中,提出了基于广义等价理论的两级瞬态CMFD加速方案,通过建立等价的单群CMFD来加速多群CMFD瞬态固定源问题的收敛,并进一步加速精细化全三维瞬态输运计算。基于CMFD并行计算框架,径向二维输运计算采用流行的特征线方法(method of characteristics, MOC),轴向一维低阶输运计算采用节块展开法(nodal expansion method, NEM)。

为了验证HNET的时空中子动力学计算功能,本文采用C5G7-TD非均匀瞬态基准题对HNET的计算精度进行验证,计算结果与高保真中子学程序MPACT进行对比。C5G7-TD系列基准题[6-7]包括一系列非均匀几何时空动力学基准题,用于校验动力学程序对非均匀几何问题的计算能力和计算精度。目前该基准问题已经发布第一阶段的系列问题,用于验证时空中子学输运算法与程序。

1 HNET瞬态计算方法

HNET采用二维MOC/一维NEM的聚合算法,并通过基于广义等价理论的两级瞬态CMFD加速全局三维瞬态输运计算,同时还采用基于MPI的区域分解并行技术,进一步缩减运行时间,更多细节见文献[5, 8-11]。

1.1 三维瞬态输运固定源方程

为了简要说明三维瞬态固定源方程的推导过程,本文从瞬态多群中子输运方程出发,在各向同性散射近似的条件下,输运方程为:

$ \begin{array}{c} \frac{1}{{{\nu _g}(\mathit{\boldsymbol{r}})}} \cdot \frac{{\partial {\varphi _g}(\mathit{\boldsymbol{r}}, \mathit{\boldsymbol{ \boldsymbol{\varOmega} }}, t)}}{{\partial t}} = - \mathit{\boldsymbol{ \boldsymbol{\varOmega} }} \cdot \nabla {\varphi _g}(\mathit{\boldsymbol{r}}, \mathit{\boldsymbol{ \boldsymbol{\varOmega} }}, t) - \\ {\mathit{\mathit{\varSigma} }_{t, g}}(\mathit{\boldsymbol{r}}, t){\varphi _g}(\mathit{\boldsymbol{r}}, \mathit{\boldsymbol{ \boldsymbol{\varOmega} }}, t) + {S_{s, g}}(\mathit{\boldsymbol{r}}, \mathit{\boldsymbol{ \boldsymbol{\varOmega} }}, t) + \\ \frac{1}{{4{\rm{ \mathit{ π} }}}}{S_F}(\mathit{\boldsymbol{r}}, t)\left[ {{\chi _{p, g}} + \sum\limits_k {{\beta _k}} \left( {{\chi _{d, k, g}} - {\chi _{p, g}}} \right)} \right] + \\ \frac{1}{{4{\rm{ \mathit{ π} }}}}\sum\limits_k {{\chi _{d, k, g}}} \left[ {{\lambda _k}{C_k}(\mathit{\boldsymbol{r}}, t) - {\beta _k}{S_F}(\mathit{\boldsymbol{r}}, t)} \right] \end{array} $ (1)

缓发中子先驱核方程为:

$ \frac{\partial C_{k}(\boldsymbol{r}, t)}{\partial t}=\beta_{k} S_{F}(\boldsymbol{r}, t)-\lambda_{k} C_{k}(\boldsymbol{r}, t) $ (2)

式中:v为中子飞行速度;φ为角通量密度;Ck为第k缓发群的先驱核浓度;SF为经过稳态修正后的裂变源项;Ss, g为散射源项;χpχd分别为瞬发和缓发中子能谱;βk为缓发中子份额;λk为缓发中子先驱核衰变常数;r为空间向量;Ω为角度向量;g为能群;Σt, g为总截面。

对瞬态角通量做指数变换:

$ {{\mathit{ φ }}_g}(\mathit{\boldsymbol{r}}, \mathit{\boldsymbol{ \boldsymbol{\varOmega} }}, t) = {\tilde \varphi _g}(\mathit{\boldsymbol{r}}, \mathit{\boldsymbol{ \boldsymbol{\varOmega} }}, t){{\rm{e}}^{{\alpha ^n}(\mathit{\boldsymbol{r}}, \mathit{\boldsymbol{ \boldsymbol{\varOmega} }})\left( {t - {t_{n - 1}}} \right)}}, t \in \left[ {{t_{n - 1}}, {t_n}} \right] $ (3)

其中,系数αn为:

$ {\alpha ^n}(\mathit{\boldsymbol{r}}, \mathit{\boldsymbol{ \boldsymbol{\varOmega} }}) = \frac{1}{{\Delta {t_{n - 1}}}}{\rm{ln}}\frac{{p\left( {{t_{n - 1}}} \right)}}{{p\left( {{t_{n - 2}}} \right)}} $ (4)

将式(4)代入式(1)~(3)中,并且对时间导数项采用全隐式方法,则瞬态输运方程可转化为瞬态固定源方程为:

$ \begin{array}{l} \mathit{\boldsymbol{ \boldsymbol{\varOmega} }} \cdot \nabla \varphi _g^n(\mathit{\boldsymbol{r}},\mathit{\boldsymbol{ \boldsymbol{\varOmega} }}) + \left( {\mathit{\Sigma }_{t,g}^n(\mathit{\boldsymbol{r}}) + \frac{{{\alpha ^n}}}{{{\nu _g}(\mathit{\boldsymbol{r}})}}} \right)\varphi _g^n(\mathit{\boldsymbol{r}},\mathit{\boldsymbol{ \boldsymbol{\varOmega} }}) = \\ S_{s,g}^n(\mathit{\boldsymbol{r}},\mathit{\boldsymbol{ \boldsymbol{\varOmega} }}) + \frac{{{{\tilde \chi }_g}(\mathit{\boldsymbol{r}})}}{{4{\rm{ \mathsf{ π} }}}}S_F^n(\mathit{\boldsymbol{r}}) + \frac{1}{{4{\rm{ \mathsf{ π} }}}}\sum\limits_k {{\chi _{d,k,g}}} {\lambda _k}\hat C_k^{n - 1}(\mathit{\boldsymbol{r}}) - \\ \frac{{\varphi _g^n(\mathit{\boldsymbol{r}},\mathit{\boldsymbol{ \boldsymbol{\varOmega} }})}}{{\Delta {t_n}{\nu _g}(\mathit{\boldsymbol{r}})}} + \frac{{\varphi _g^{n - 1}(\mathit{\boldsymbol{r}},\mathit{\boldsymbol{ \boldsymbol{\varOmega} }})}}{{E_g^n(\mathit{\boldsymbol{r}},\mathit{\boldsymbol{ \boldsymbol{\varOmega} }})\Delta {t_n}{\nu _g}(\mathit{\boldsymbol{r}})}} \end{array} $ (5)

其中:

$ E_g^n(\mathit{\boldsymbol{r}}, \mathit{\boldsymbol{ \boldsymbol{\varOmega} }}) = {{\rm{e}}^{ - \alpha _g^n\left( {{t_n} - {t_{n - 1}}} \right)}} $ (6)
$ \tilde{\chi}_{g}^{n}(\boldsymbol{r})=\chi_{p, g}(1-\beta)+\sum\limits_{k} \chi_{d, k, g} \lambda_{k} \beta_{k} \Delta t_{n} \kappa_{1}\left[\tilde{\lambda}_{k}^{n}(\boldsymbol{r})\right]_{k} $ (7)
$ \begin{aligned} &\hat{C}_{k}^{n-1}(\boldsymbol{r})=\mathrm{e}^{-\hat{\lambda}_{k}^{n}} C_{k}^{n-1}(\boldsymbol{r})+ \\ &\mathrm{e}^{\alpha^{n} \Delta t_{n}} \beta_{k} \Delta t_{n} S_{F}^{n-1}(\boldsymbol{r})\left[\kappa_{0}\left(\tilde{\lambda}_{k}^{n}\right)-\kappa_{1}\left(\tilde{\lambda}_{k}^{n}\right)\right] \end{aligned} $ (8)
$ \boldsymbol{\kappa}_{m}(x)=\mathrm{e}^{x}-\frac{m \boldsymbol{\kappa}_{m-1}(x)}{x} $ (9)

至此,式(6)可以使用任何标准的中子固定源问题求解器来解决。在HNET程序中,径向二维输运计算采用流行的MOC方法,轴向一维低阶输运计算采用节块展开NEM方法,并通过基于广义等价理论的两级瞬态CMFD加速全局三维瞬态输运计算。

1.2 基于广义等价理论的两级瞬态粗网有限差分加速方法

CMFD不仅为二维径向MOC输运计算和一维轴向NEM计算提供了天然的计算框架,同时还对全堆芯输运计算起到加速作用,三维多群瞬态扩散方程为:

$ \begin{array}{l} \sum\limits_q {J_{g, c, q}^{n, l}} {A_{c, q}} + {V_c}\left[ {\mathit{\boldsymbol{ \boldsymbol{\bar \varSigma} }}_{t, g, c}^n\bar \phi _{g, c}^{n, l} - \sum\limits_{{g^\prime }} {\mathit{\boldsymbol{ \boldsymbol{\bar \varSigma} }}_{{g^\prime }, g, c}^n} \bar \phi _{{g^\prime }, c}^{n, l}} \right] = \\ {V_c}\left( {\frac{{\bar \chi _g^n}}{{k_{{\rm{eff }}}^{ss}}}\sum\limits_{{g^\prime }} \nu \mathit{\boldsymbol{ \boldsymbol{\bar \varSigma} }}_{f, {g^\prime }, c}^n\bar \phi _{{g^\prime }, c}^{n, l - 1} + \bar S_{tr, g, c}^n} \right) \end{array} $ (10)

式中:c为中心节块;q为相邻节块;n为时间步;l为迭代次数;Sdcn-1Sdtn-1为上一时间步提供的缓发中子源项;νΣfnΣg, gnχnΣt, gn分别为节块内的平均裂变截面、散射截面、裂变谱和总截面;keffss为稳态输运计算提供的特征值。

根据基于广义等价理论的粗网有限差分方法,中心节块c到相邻节块q的界面流为:

$ J_{g, c, q} \cdot A_{c, q}=\tilde{D}_{g, c, q} \cdot \bar{\phi}_{g, c}-\tilde{D}_{g, q, c} \cdot \bar{\phi}_{g, q} $ (11)

其中$\tilde D$的定义为:

$ \tilde{D}_{g, c, q}=A_{c, q} \frac{2 f_{g, c, q}^{\mathrm{dis}}}{\frac{f_{g, c, q}^{\mathrm{dis}} h_{c, q}}{f_{g, c, q}^{\mathrm{dif}} D_{g, c}}+\frac{f_{g, q, c}^{\mathrm{dis}} h_{q, c}}{f_{g, q, c}^{\mathrm{dif}} D_{g, q}}} $ (12)
$ \tilde{D}_{g, q, c}=A_{c, q} \frac{2 f_{g, q, \mathrm{c}}^{\mathrm{dis}}}{\frac{f_{g, c, q}^{\mathrm{dis}} h_{c, q}}{f_{g, c, q}^{\mathrm{dif}} D_{g, c}}+\frac{f_{g, q, c}^{\mathrm{dis}} h_{q, c}}{f_{g, q, c}^{\mathrm{dif}} D_{g, q}}} $ (13)

对于边界节块,其对应边界界面的$\tilde D$为:

$ \tilde{D}_{g, c, q}=\frac{2 f_{g, c, q}^{\mathrm{dis}}}{\frac{f_{g, c, q}^{\mathrm{dis}} h_{c, q}}{f_{g, c, q}^{\mathrm{dif}} D_{g, c}}+\frac{2}{\tilde{\alpha}}} $ (14)

式中:hc, q为中心节块在对应方向上的厚度;Dg为扩散系数;$\tilde \alpha $为边界上的反射率。同时,fc, qdisfq, cdif分别为节块不连续因子和扩散系数修正因子,其来自于径向MOC与轴向NEM。

在HNET的二维/一维算法中,使用$\tilde D$强制低阶算法的界面流与高阶算法相等,因此在全局收敛之后CMFD计算得到的节块通量和界面流能够与输运计算的结果保持一致。同时,基于二维MOC计算获得的中子通量信息,采用通量-体积加权平均方法可计算得到不同节块内的均匀化截面。这样,便可以建立基于广义平衡理论的全局多群CMFD线性系统为:

$ \begin{aligned} &V_{c}\left(-\sum\limits_{g^{\prime}} \mathit{\bar{\Sigma}}_{g^{\prime}, g, c}^{n} \bar{\phi}_{g^{\prime}, c}^{n, l}+\mathit{\bar{\Sigma}}_{t, g, c}^{n} \bar{\phi}_{g, c}^{n, l}\right)+ \\ &\left(\tilde{D}_{g, c, \mathrm{~W}}^{n}+\tilde{D}_{g, c, \mathrm{~N}}^{n}+\tilde{D}_{g, c, \mathrm{E}}^{n}+\tilde{D}_{g, c, \mathrm{~S}}^{n}+\tilde{D}_{g, c, \mathrm{~T}}^{n}+\right. \\ &\left.\tilde{D}_{g, c, \mathrm{~B}}^{n}\right) \bar{\phi}_{g, c}^{n, l}-\left(\tilde{D}_{g, \mathrm{~W}, c}^{n} \bar{\phi}_{g, \mathrm{~W}}^{n, l}+\tilde{D}_{g, \mathrm{~N}, c}^{n} \bar{\phi}_{g, \mathrm{~N}}^{n, l}+\right. \\ &\left.\tilde{D}_{g, \mathrm{E}, c}^{n} \bar{\phi}_{g, \mathrm{E}}^{n, l}+\tilde{D}_{g, \mathrm{~S}, c}^{n} \bar{\phi}_{g, \mathrm{~S}}^{n, l}+\tilde{D}_{g, \mathrm{~T}, c}^{n} \bar{\phi}_{g, \mathrm{~T}}^{n, l}+\tilde{D}_{g, \mathrm{~B}, c}^{n} \bar{\phi}_{g, \mathrm{~B}}^{n, l}\right)= \\ &V_{c}\left(\frac{\bar{\chi}_{g}^{n}}{k_{\mathrm{eff}}^{s s}} \sum\limits_{g^{\prime}} \nu \mathit{\bar{\Sigma}}_{f, g^{\prime}, c}^{n} \bar{\phi}_{g^{\prime}, c}^{n, l-1}+\bar{S}_{\mathrm{tr}, g, c}^{n}\right) \end{aligned} $ (15)

式中:W、E、N、S、T、B代表节块在西、东、北、南、顶、底方向上的6个表面。

式(15)可以由任意固定源问题求解器来求解,且所用的特征值来自于稳态输运计算,其在整个瞬态过程中应保持不变。

由于多群CMFD线性系统的条件数很大,如若采用Wielandt shift方法加速逆幂迭代,其条件数会变得更大,因此,一级多群CMFD线性系统存在收敛慢的缺陷。一群CMFD线性系统条件数较低使得其容易求解,可利用一群CMFD获得中子通量及裂变源加速多群CMFD。基于多群CMFD获得的中子通量、中子流及多群截面信息,可获得一群CMFD线性系统为:

$ \begin{gathered} \left(\tilde{D}_{c, \mathrm{~W}}^{n}+\tilde{D}_{c, \mathrm{~N}}^{n}+\tilde{D}_{c, \mathrm{E}}^{n}+\tilde{D}_{c, \mathrm{~S}}^{n}+\tilde{D}_{c, \mathrm{~T}}^{n}+\tilde{D}_{c, \mathrm{~B}}^{n}\right) \bar{\phi}_{c}^{n, l}- \\ \left(\tilde{D}_{\mathrm{W}, c}^{n} \bar{\phi}_{\mathrm{W}}^{n, l}+\tilde{D}_{\mathrm{N}, c}^{n} \bar{\phi}_{\mathrm{N}}^{n, l}+\tilde{D}_{\mathrm{E}, c}^{n} \bar{\phi}_{\mathrm{E}}^{n, l}+\right. \\ \left.\tilde{D}_{\mathrm{S}, c}^{n} \bar{\phi}_{\mathrm{S}}^{n, l}+\tilde{D}_{\mathrm{T}, c}^{n} \bar{\phi}_{\mathrm{T}}^{n, l}+\tilde{D}_{\mathrm{B}, c}^{n} \bar{\phi}_{\mathrm{B}}^{n, l}\right)+ \\ V_{c}\left(\mathit{\bar{\Sigma}}_{a, c}^{n}+\frac{\alpha^{n}}{v_{g}}\right) \bar{\phi}_{c}^{n, l}=V_{c}\left(\frac{1}{k_{\mathrm{eff}}^{s s}} \nu \mathit{\bar{\Sigma}}_{f f}^{n} \bar{\phi}_{c}^{n, l-1}+\bar{S}_{\mathrm{tr}, c}^{n}\right) \end{gathered} $ (16)

一群CMFD用以加速多群CMFD的特征值和裂变源项,求解多群CMFD线性系统,以计算得到新的多群中子通量分布和特征值,进而用于更新一群宏观截面、中子通量和中子流。同时,计算得到的多群中子通量将进一步更新每一个节块的界面流,而界面流则会被用于更新一群CMFD的节块不连续因子和扩散系数修正因子,直到整个线性系统收敛。其中,当残差范数下降到指定标准或最大迭代次数达到时,一群CMFD迭代或多群CMFD迭代相应退出,整体收敛通过检查不连续因子的相对误差是否达到收敛标准,具体细节可参考文献[7]。

2 C5G7-TD基准题模型

C5G7-TD系列基准题二维算例包含TD0、TD1、TD2和TD3系列;三维算例包含TD4和TD5系列。每个系列都包含数个算例,这些算例包含控制棒移动和慢化剂密度变化等瞬态事件。控制棒移动是以组件为单位同时移动,即在同一组件内所有控制棒均保持相同的高度。二维与三维的瞬态问题,更多细节见文献[5]。

2.1 二维算例模型

1) TD0系列。

TD0系列算例描述了二维情况下控制棒瞬间插入与提出反应堆的过程。在初始时刻,控制棒突然插入到堆芯10%的位置,并停留1 s。随后,控制棒提出到堆芯5%的位置,并再次停留1 s。最后,在2 s结束时,控制棒从堆芯中完全撤出。在二维模型中,这一过程被模拟为控制棒区1区截面的变化。

2) TD1与TD2系列。

TD1与TD2系列描述了二维情况下控制棒线性插入与提出反应堆的过程,通过线性地将导向管截面变化成控制棒截面来模拟控制棒移动过程。在初始时刻,所有的控制棒均在堆芯外部,在0~1 s过程中,控制棒从堆芯顶部匀速插入到堆芯的特定位置。其中,TD1为堆芯5%的位置,TD2为堆芯10%的位置。然后以相同的速度提出控制棒,并在2 s结束时完全提出堆芯区域。TD1和TD2基准题各有5道算例。

3) TD3系列。

TD3系列描述了二维情况下反应堆内慢化剂密度变化过程。在起始时刻,燃料组件内的慢化剂密度(不包括轴向反射层)为正常密度,在0~1 s过程中,慢化剂密度以一定速度降低。在1 s结束时,密度达到最小值,即初始密度的w值。然后慢化剂密度以同样的速度增加,并在2 s结束时回到初始值。对于TD3中的不同算例,w的值有所区别,其中“TD3-1”的w为0.95;“TD3-2”的w为0.90;“TD3-3”的w为0.85;“TD3-4”的w为0.80。

2.2 三维算例模型

1) TD4系列。

TD4系列描述了三维情况下控制棒线性插入与提出反应堆的过程。在初始时刻,所有控制棒均处于顶部反射层内。随后在0~8 s过程中,不同组件内的控制棒以不同的速度匀速地插入与提出反应堆。TD4系列共包含5个算例,分别模拟不同组件内控制棒不同的移动情况。

2) TD5系列。

TD5系列描述了三维情况下反应堆内慢化剂密度变化过程。在整个瞬态过程中,所有控制棒均保持在燃料组件顶部的反射层中。随后在0~8 s过程中,不同组件内的慢化剂密度以组件为单位同步地进行变化。TD5系列共包含4个算例,分别模拟不同组件内慢化剂变化情况。

3 HNET数值结果与分析

对于二维算例,本文展示TD1、TD2与TD3的计算结果;对于三维算例,本文展示TD4与TD5的计算结果。对于二维算例,瞬态输运计算的步长为10 ms,对于三维算例,时间步长为25 ms。对于所有算例,径向网格的划分方式是一致的。对于燃料组件内的非均匀栅元,在燃料区域径向划分3圈且周向划分8区,在慢化剂区域径向划分1圈且周向划分8区;对于反射层内的慢化剂栅元,靠近燃料组件的部分栅元划分为6×6的平源区网格,远离燃料组件的部分栅元划分为1×1的平源区网格。对于所有三维算例,轴向网格尺寸为5.355 cm,特征线间距为0.03 cm,轴向求积组采用Yamamoto求积组[12],半空间内采用64个方位角与3个极角的组合。特征值的收敛准则为1.0×10-6,中子标通量采用无穷范数误差,其收敛准则为1.0×10-4。算例所使用的CPU型号均为2.60 GHz Intel Xeon E5-2690 v4,二维算例均采用径向9核的区域分解并行计算方案,三维维算例均采用轴向32核的区域分解并行计算方案。为了保证数值结果对比的有效性,上述参数除求积组和CPU型号外均与MPACT的参考解保持一致。

3.1 二维算例计算结果

对于C5G7-TD基准题的稳态特征,MPACT程序的参考解为1.186 673,HNET程序的计算结果为1.186 867,HNET结果与参考解的相对误差约为0.016%,表明HNET计算得到的初始时刻的堆芯状态是正确的。二维算例的堆芯相对功率历史、与MPACT参考解的堆芯功率误差和反应性历史如图 1~9所示。

Download:
图 1 TD1系列堆芯相对功率历史 Fig. 1 TD1 fractional core power history
Download:
图 2 对比MPACT程序的TD1系列堆芯功率相对误差 Fig. 2 Relative errors of TD1 fractional core power compared to MPACT results
Download:
图 3 TD1系列反应性历史 Fig. 3 TD1 reactivity history
Download:
图 4 TD2系列堆芯相对功率历史 Fig. 4 TD2 fractional core power history
Download:
图 5 对比MPACT程序的TD2系列堆芯功率相对误差 Fig. 5 Relative error of TD2 fractional core power compared to MPACT results
Download:
图 6 TD2系列反应性历史 Fig. 6 TD2 reactivity history
Download:
图 7 TD3系列堆芯相对功率历史 Fig. 7 TD3 fractional core power history
Download:
图 8 对比MPACT程序的TD3系列堆芯功率相对误差 Fig. 8 Relative error of TD3 fractional core power compared to MPACT results
Download:
图 9 TD3系列反应性历史 Fig. 9 TD3 reactivity history

对于TD1、TD2和TD3系列,在0~1 s和1~2 s的过程中,宏观截面以相同的速度增减变化,使得裂变率和反应性的增减也出现了对称的趋势。同时,TD2系列的控制棒插入深度大于TD1系列,导致TD2系列引入了数值更大的负反应性。图 1~9表明相对功率和反应性的变化规律与MPACT的结果相符。同时,表 1列出了每个二维算例的10 s瞬态过程的运行时间并与MPACT参考解对比。

表 1 TD1~TD3系列算例运行时间对比 Table 1 Comparison of total run times for TD1~TD3 

对于TD1系列,功率峰值的最大偏差为0.48%;功率的最大偏差为0.48%,发生在“TD1-5”的1 s时刻;反应性峰值最大偏差为0.26%,发生在“TD1-1”。对于TD2系列,功率峰值的最大偏差为0.21%;功率的最大偏差为0.14%,发生在“TD2-5”的10 s时刻;反应性峰值的最大偏差为1.26%,发生在“TD2-3”。对于TD3系列,功率峰值的最大偏差为-0.72%;功率的最大偏差为-0.74%,发生在“TD3-1”的1.05 s时刻;反应性峰值的最大偏差为1.14%,发生在“TD3-1”。

3.2 三维算例计算结果

对于三维问题TD4和TD5系列,堆芯相对功率变化和反应性变化见图 10~15。其中,“TD4-4”和“TD4-5”描述了不同控制棒同时插入和提出的过程,这导致了相对复杂的裂变率变化和反应性变化。表 2给出了所有三维算例的运行时间并与MPACT参考解对比。

Download:
图 10 TD4系列堆芯相对功率历史 Fig. 10 TD4 fractional core power history
Download:
图 11 对比MPACT程序的TD4系列堆芯功率相对误差 Fig. 11 Relative error of TD4 fractional core power compared to MPACT results
Download:
图 12 TD4系列反应性历史 Fig. 12 TD4 reactivity history
Download:
图 13 TD5系列堆芯相对功率历史 Fig. 13 TD5 fractional core power history
Download:
图 14 对比MPACT程序的TD5系列堆芯功率相对误差 Fig. 14 Relative error of TD5 fractional core power compared to MPACT results
Download:
图 15 TD5系列反应性历史 Fig. 15 TD5 reactivity history
表 2 TD4、TD5算例的运行时间 Table 2 Run time summaries for TD4, TD5 

对于TD4系列,功率峰值的最大偏差为-0.21%,出现在“TD4-1”;功率的最大偏差为-0.55%,发生在“TD4-5”的3.375 s时刻;稳定之后相对偏差为0.10%。反应性峰值最大相对偏差为2.21%,发生在“TD4-4”。对于TD5系列,功率峰值的最大偏差为-0.21%,出现在“TD5-1”;功率最大偏差为-0.23%,发生在“TD5-1”的2.3 s时刻;稳定之后相对偏差为-0.05%。反应性峰值最大相对偏差为0.66%,发生在“TD5-4”。

4 结论

1) 对于二维问题,在瞬态事件期间堆芯功率的相对误差有所波动,瞬态事件结束后相对误差逐渐趋于平稳。虽然计算模型的细微差别可能引入了一些偏差,如空间几何模型、离散策略、动力学参数的处理和加速方法等,但是总体而言,HNET的数值结果与MPACT基准解吻合得很好,堆芯功率的相对误差最大不超过0.8%,在保持足够精度的同时HNET程序算例拥有明显的效率优势。

2) 对于三维问题,在TD4问题上HNET的相对误差不超过0.55%,在TD5问题上相对误差小于0.23%,在输运时间步长相同的情况下,HNET程序的计算总耗时更少,与MPACT程序相比HNET程序只需要50%~85%的计算时间便可实现相同精度的瞬态输运计算。

3) 总体而言,HNET瞬态过程中堆芯总功率变化与反应性变化均与MPACT符合良好,已完全具备三维瞬态精细化计算的能力。

同时,用于瞬态计算的热工反馈功能正在研发中,未来将进一步测试HNET的时空动力学计算能力。

参考文献
[1]
JOO H G, CHO J Y, KIM K S, et al. Methods and performance of a three-dimensional whole-core transport code DeCART[C]. PHYSOR 2004-The Physics of Fuel Cycles and Advanced Nuclear Systems, Chicago: Global Developments, 2004. (0)
[2]
CHO N Z, LEE G S, PARK C J. Fusion of method of characteristics and nodal method for 3-D whole-core transport calculation[J]. Transactions of the American nuclear society, 2002, 86: 322-324. (0)
[3]
COLLINS B, STIMPSON S, KELLEY B W, et al. Stability and accuracy of 3D neutron transport simulations using the 2D/1D method in MPACT[J]. Journal of computational physics, 2016, 326: 612-628. DOI:10.1016/j.jcp.2016.08.022 (0)
[4]
吴文斌. 基于并行技术的2D/1D耦合三维全堆输运方法研究[D]. 北京: 清华大学, 2014.
WU Wenbin. Research on 2D/1D Coupled 3D whole-core transport calculations based on parallel computation[D]. Beijing: Tsinghua University, 2014. (0)
[5]
HAO Chen, KANG Le, XU Yunlin, et al. 3D whole-core neutron transport simulation using 2D/1D method via multi-level generalized equivalence theory based CMFD acceleration[J]. Annals of nuclear energy, 2018, 122: 79-90. DOI:10.1016/j.anucene.2018.08.014 (0)
[6]
HOU J, IVANOV K, BOYARINOV V, et al, 2017. OECD/NEA benchmark for time-dependent neutron transport calculations without spatial homogenization. Nuclear Engineering and Design 317, 177-189. (0)
[7]
HOU J, IVANOV K N, BOYARINOV V F, et al. OECD/NEA benchmark for time-dependent neutron transport calculations without spatial homogenization[J]. Nuclear engineering and design, 2017, 317: 177-189. DOI:10.1016/j.nucengdes.2017.02.008 (0)
[8]
SHEN Qicang, WANG Yiran, JABAAY D, et al. Transient analysis of C5G7-TD benchmark with MPACT[J]. Annals of nuclear energy, 2019, 125: 107-120. DOI:10.1016/j.anucene.2018.10.049 (0)
[9]
HAO Chen, XU Yunlin, DOWNAR T J. Multi-level coarse mesh finite difference acceleration with local two-node nodal expansion method[J]. Annals of nuclear energy, 2018, 116: 105-113. DOI:10.1016/j.anucene.2018.02.002 (0)
[10]
KANG Le, HAO Chen, ZHAO Qiang, et al. Two-level time-dependent GET based CMFD acceleration for 3D whole core transient transport simulation using 2D/1D method[J]. Annals of nuclear energy, 2020, 142: 107405. DOI:10.1016/j.anucene.2020.107405 (0)
[11]
KANG Le, HAO Chen, ZHAO Qiang, et al. The advanced multilevel predictor-corrector quasi-static method for pin-resolved neutron kinetics simulation[J]. Frontiers in energy research, 2021, 9: 747148. DOI:10.3389/fenrg.2021.747148 (0)
[12]
YAMAMOTO A, TABUCHI M, SUGIMURA N, et al. Derivation of optimum polar angle quadrature set for the method of characteristics based on approximation error for the Bickley function[J]. Journal of nuclear science and technology, 2007, 44(2): 129-136. DOI:10.1080/18811248.2007.9711266 (0)