2. 中国科学技术大学 合肥 230026
2. University of Science and Technology of China, Hefei 230026, China
在反应堆设计和物理分析中,确定临界状态下的硼浓度、控制棒位置等参数非常重要,该过程叫做临界搜索[1]。由于蒙特卡罗法具有适应性强、计算结果准确、易于并行计算等优势,故随着计算机技术的不断发展,蒙特卡罗法在临界搜索中的应用越来越广泛。
传统临界搜索主要采用Newton-Raphson法,它通过不断尝试,进行多次临界计算,直到获得符合临界条件的参数[1]。由于蒙特卡罗临界计算较耗时,而传统方法需要反复进行临界计算,所以不适合蒙特卡罗临界搜索。为了提高临界搜索的效率,国内外学者对蒙特卡罗临界搜索法进行了一些研究。例如,北京应用物理与计算数学研究所邓力等基于中子平衡理论实现临界硼浓度搜索[2];而更多的学者基于微扰理论展开了相关研究,密歇根大学Aaron等[3]使用相关抽样法实现了临界硼浓度搜索;清华大学李泽光等[4]使用微分算符法实现了临界硼浓度搜索。
蒙特卡罗微扰计算方法主要有三类[4]:共轭通量法、微分算符法以及相关抽样法。其中,微分算符法和相关抽样法是传统微扰计算方法,这两种方法在计算微扰对有效增值因子影响时必须要考虑裂变源扰动效应,而裂变源扰动的处理非常困难[5]。共轭通量法不需要考虑裂变源扰动效应,使用时更加方便,并且共轭通量计算结果与考虑裂变源扰动效应的一阶微分算符法结果等效[6]。因此,本文选取共轭通量法作为计算微扰问题的方法进行临界硼浓度搜索的研究。
本文基于FDS凤麟核能团队自主研发的中子输运设计与安全评价软件系统SuperMC[7-11] (Super Multi-functional Calculation Program for Nuclear Design and Safety Evaluation)开展。FDS凤麟核能团队在中子学分析与程序开发、堆芯设计[12-15]与材料研究[16-19]、核能安全研究、实验及核技术应用[20-21]等方面具有丰富的研究经验和突出的研究成果。本文研究在蒙特卡罗程序中使用反复裂变几率计算共轭通量的问题,发展了一种基于共轭通量的蒙特卡罗临界硼浓度搜索方法。该方法只需对系统进行一次临界计算,就可以得到系统初始的共轭通量以及有效增殖因子对硼浓度变化的响应系数,进而计算得到临界硼浓度。
1 共轭通量计算方法基于共轭通量进行临界硼浓度搜索,关键在于共轭通量的求解。共轭通量的物理意义是中子价值,本身不具有“通量”的含义。中子价值的定义[22]为:在临界反应堆中,一个运动方向为Ω、能量为E的中子投放到位置r处所引起的对该系统稳定功率的贡献。1965年,Hurwitz[23]指出:系统在临界状态下,由一个中子引起的每代裂变反应次数在中子代数趋于无穷时会趋近于一个极限值,该极限值即为此中子引发的反复裂变几率,等同于中子价值,可以使用反复裂变几率替代共轭通量。
根据反复裂变几率的物理含义,蒙特卡罗程序模拟粒子输运的过程可以分为三个阶段[24]:第一阶段,统计源中子的个数;第二阶段,跟踪源中子发生裂变反应产生的后代中子;第三阶段,统计属于相同的源中子的后代中子的个数,得到反复裂变几率的大小。反复裂变几率示意图如图 1所示。
微扰理论的基本方程[22]:
$ \left\langle {{\phi ^ * }, \mathit{\boldsymbol{P}}\phi } \right\rangle {\rm{ = }}0 $ | (1) |
式中:ϕ=ϕ(r, E, Ω)为中子通量密度;ϕ*=ϕ*(r, E, Ω)为共轭中子通量密度;P为扰动算符,包括裂变源项算子F、除裂变源项算子以外的其他输运算子L和keff等在内的系统内各种参数的微小扰动。
假设由于某种原因在系统内引入某种扰动,算子L及F获得相应扰动。为保持系统临界,其特征值λ(λ=1/keff)亦给予相应的扰动δλ,即L、F、λ扰动后分别为:
$ \mathit{\boldsymbol{L'}} \to \mathit{\boldsymbol{L}} + {\rm{ \mathsf{ δ} }}\mathit{\boldsymbol{L}}{\bf{、}}\mathit{\boldsymbol{F'}} \to \mathit{\boldsymbol{F}} + {\rm{ \mathsf{ δ} }}\mathit{\boldsymbol{F}}{\bf{、}}\lambda ' \to \lambda + {\rm{ \mathsf{ δ} }}\lambda $ |
因而扰动算符P=δL-λδF-δλF[25],根据扰动方程(1),得到:
$ \left\langle {{\phi ^*},{\rm{ \mathsf{ δ} }}\mathit{\boldsymbol{L}}\phi } \right\rangle - \left\langle {{\phi ^*},\lambda {\rm{ \mathsf{ δ} }}\mathit{\boldsymbol{F}}\phi } \right\rangle - \left\langle {{\phi ^*},{\rm{ \mathsf{ δ} }}\lambda \mathit{\boldsymbol{F}}\phi } \right\rangle = 0 $ | (2) |
根据反应性与有效增殖因子的关系
$ {\rm{ \mathsf{ δ} }}\rho = - \frac{{\left\langle {{\phi ^*},{\rm{ \mathsf{ δ} }}\mathit{\boldsymbol{L}}\phi } \right\rangle - \left\langle {{\phi ^*},\lambda {\rm{ \mathsf{ δ} }}\mathit{\boldsymbol{F}}\phi } \right\rangle }}{{\left\langle {{\phi ^*},\mathit{\boldsymbol{F}}\phi } \right\rangle }} $ | (3) |
这便是中子输运方程的扰动方程,δL及δF扰动包括各种截面的变化,如:总的反应截面的变化量δ∑t、裂变反应截面的变化量δ∑f、散射反应截面的变化量δ∑s等。把这些扰动带入式(3)便得到扰动方程:
$ \begin{array}{*{20}{l}} {{\rm{ \mathsf{ δ} }}\rho = - \left\{ {\int_0^\infty {{\rm{d}}E} } \right.\int_\mathit{\boldsymbol{ \boldsymbol{\varOmega} }} {{\rm{d}}\mathit{\boldsymbol{ \boldsymbol{\varOmega} }}\int_V {{\phi ^*}\left[ {\phi {\rm{ \mathsf{ δ} }}} \right.} } {\mathit{\Sigma }_{\rm{t}}} - }\\ {\;\;\;\;\;\;\;\int_{\mathit{\boldsymbol{ \boldsymbol{\varOmega} '}}} {{\rm{d}}\mathit{\boldsymbol{ \boldsymbol{\varOmega} '}}} \int_0^\infty {\phi \left( {\mathit{\boldsymbol{r}},E',\mathit{\boldsymbol{ \boldsymbol{\varOmega} '}}} \right)} }\\ {\;\;\;\;\;\;\;{\rm{ \mathsf{ δ} }}{\mathit{\Sigma }_{\rm{s}}}\left( {\mathit{\boldsymbol{r}};E',\mathit{\boldsymbol{ \boldsymbol{\varOmega} '}} \to E,\mathit{\boldsymbol{ \boldsymbol{\varOmega} }}} \right){\rm{d}}E' - \frac{\lambda }{{4{\rm{ \mathsf{ π} }}}}\int_{\mathit{\boldsymbol{ \boldsymbol{\varOmega} '}}} {{\rm{d}}\mathit{\boldsymbol{ \boldsymbol{\varOmega} '}}} }\\ {\;\;\;\;\;\;\;\left. {\left. {\int_0^\infty {{\rm{d}}E'{\rm{ \mathsf{ δ} }}\left( {\chi \left( E \right)\nu {\mathit{\Sigma }_{\rm{f}}}} \right)\phi \left( {\mathit{\boldsymbol{r}},E',\mathit{\boldsymbol{ \boldsymbol{\varOmega} '}}} \right)} } \right]{\rm{d}}\mathit{\boldsymbol{r}}} \right\}/Q} \end{array} $ | (4) |
$ \begin{array}{*{20}{l}} {Q = \left\langle {{\phi ^*},\mathit{\boldsymbol{F}}\phi } \right\rangle = \frac{1}{{4{\rm{ \mathsf{ π} }}}}\int_V {{\rm{d}}\mathit{\boldsymbol{r}}} \int_0^\infty {{\rm{d}}E\int_\mathit{\boldsymbol{ \boldsymbol{\varOmega} }} {{\phi ^*}\left( {\mathit{\boldsymbol{r}},E,\mathit{\boldsymbol{ \boldsymbol{\varOmega} }}} \right)} } {\rm{d}}\mathit{\boldsymbol{ \boldsymbol{\varOmega} }}}\\ {\;\;\;\;\;\;\;\int_0^\infty {{\rm{d}}E'} \int_{\mathit{\boldsymbol{ \boldsymbol{\varOmega} '}}} {\chi \left( E \right)\nu {\mathit{\Sigma }_{\rm{f}}}\left( {\mathit{\boldsymbol{r}},E'} \right)} \times \phi \left( {\mathit{\boldsymbol{r}},E',\mathit{\boldsymbol{ \boldsymbol{\varOmega} '}}} \right){\rm{d}}\mathit{\boldsymbol{ \boldsymbol{\varOmega} '}}} \end{array} $ | (5) |
式中:Q为系统内裂变中子的总价值,对于给定稳定系统,它等于常数。
基于上述微扰理论,本文使用反复裂变几率法计算共轭通量,发展了临界硼浓度搜索的方法。对于式(4)左侧δρ,可以由目标值keff=1和初始keff计算得到;式(4)右侧,由于硼浓度的变化对散射截面的变化δ∑s和裂变截面的变化δ∑f以及其他截面的变化影响较小,可以忽略不计,所以δ∑t=δ∑a+δ∑s≈δ∑a,式(4)可以简化为:
$ - \left( {1 - \frac{1}{{{k_{{\rm{eff}}}}}}} \right) = - \left\{ {\int_0^\infty {{\rm{d}}E\int_\mathit{\boldsymbol{ \boldsymbol{\varOmega} }} {{\rm{d}}\mathit{\boldsymbol{ \boldsymbol{\varOmega} }}\int_V {{\phi ^*}\phi {\rm{ \mathsf{ δ} }}{\mathit{\Sigma }_{\rm{a}}}{\rm{d}}\mathit{\boldsymbol{r}}} } } } \right\}/Q $ | (6) |
通过反复裂变几率方法获取每个栅元、能群的共轭中子通量密度ϕ*(r, E)。同时可以获得系统中慢化剂中硼对中子的吸收率Ra(r, E)和keff。对于式(6)中的Q值,归一化后可以取渐进代裂变中子的权重w。将计算出来的各个变量带入式(6)中,可以得到:
$ \begin{array}{*{20}{l}} { - \left( {1 - \frac{1}{{{k_{{\rm{eff}}}}}}} \right) = - \Sigma {\phi ^*}\left( {\mathit{\boldsymbol{r}},E} \right)\phi \left( {\mathit{\boldsymbol{r}},E} \right)\frac{{{R_{\rm{a}}}\left( {\mathit{\boldsymbol{r}},E} \right)}}{{\phi \left( {\mathit{\boldsymbol{r}},E} \right)}}{\rm{ \mathsf{ δ} }}{N_{\rm{b}}}/w}\\ {\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\; = - \Sigma {\phi ^*}\left( {\mathit{\boldsymbol{r}},E} \right){R_{\rm{a}}}\left( {\mathit{\boldsymbol{r}},E} \right){\rm{ \mathsf{ δ} }}{N_{\rm{b}}}/w} \end{array} $ | (7) |
其中:δNb表示由初始状态到临界状态(即keff=1)下的硼核子密度的变化量,可得到:
$ {\rm{ \mathsf{ δ} }}{N_{\rm{b}}} = w\left( {1 - \frac{1}{{{k_{{\rm{eff}}}}}}} \right)/\left[ {\Sigma {\phi ^*}\left( {\mathit{\boldsymbol{r}},E} \right){R_{\rm{a}}}\left( {\mathit{\boldsymbol{r}},E} \right)} \right] $ | (8) |
系统初始状态的硼核子密度为Nb,达到临界状态硼核子密度需要的变化量为δNb,所以临界硼核子密度Nb'为:
$ {N'_{\rm{b}}} = {N_{\rm{b}}} + {\rm{ \mathsf{ δ} }}{N_{\rm{b}}} $ | (9) |
基于上文对临界硼浓度搜索方法的描述,本文实现的临界硼浓度搜索流程图如图 2所示。
为了验证本文发展方法的正确性,选取了西屋公司标准17×17组件模型[26]进行数值验证。
3.1 例题描述与计算基于SuperMC的层次结构建模功能[27-29]对模型的几何、材料、源项、计数等进行建模,如图 3所示。
组件长和宽均为21.42cm,高为20cm,四周采用反射边界条件,燃料材料UO2的密度为30.196g·cm-3,硼水作慢化剂。在KCODE模式下,源中子数为50000,计算200个非活跃代,400个活跃代。本次测试中均使用ENDF/B-Ⅶ.0数据库,采用栅元计数模式统计各个栅元中热群、中间能群、快群的共轭通量及硼对中子的吸收率。
通过计算得出临界硼浓度以后,改变输入文件中的硼浓度,进行临界搜索结果验证。
3.2 计算结果及分析为了验证该方法的正确性,选取多个初始硼浓度,进行临界搜索。初始硼浓度、初始keff值以及通过上述方法搜索得到的临界硼浓度和验证计算的keff的值如表 1所示。
通过表 1可知,当初始keff与1相差0.02之内时,通过本文方法得到的硼浓度,可以使keff与1相差小于0.0005,结果比较理想;由于该方法只考虑了一阶扰动的影响,当初始keff与1相差大于0.02时,用搜索到的硼浓度计算得到的keff与1尚有一定差异,需要进一步考虑高阶扰动的影响。
4 结语为了克服传统临界搜索方法效率低的缺点,本文基于超级蒙特卡罗核模拟软件系统SuperMC发展了一种高效的蒙特卡罗临界硼浓度搜索法。该方法选取共轭通量法作为计算微扰问题的方法,并使用反复裂变几率的统计结果作为共轭通量的估计,仅需一次临界计算,可以获得系统的共轭通量和初始有效增殖因子对硼浓度的响应系数,进而求得临界硼浓度。经验证,该方法计算结果准确,表明了本文方法的正确性与可行性。
致谢 本文工作得到了FDS凤麟核能团队其他成员的帮助和大力支持,在此深表感谢![1] |
Dall'Osso A. A neutron balance approach in critical parameter determination[J]. Annals of Nuclear Energy, 2008, 35(9): 1686-1694. DOI:10.1016/j.anucene.2008.02.008 |
[2] |
Li R, Zhang L Y, Shi D F, et al. Criticality search of soluble boron iteration in MC code JMCT[C]. International Youth Nuclear Congress (IYNC 2016), Hangzhou, China, 2016. DOI: 10.1016/j.egypro.2017.08.118. https://www.researchgate.net/publication/319950518_Criticality_search_of_soluble_boron_iteration_in_MC_code_JMCT
|
[3] |
Aaron M B, William R M. Correlated sampling Monte Carlo for critical boron search[C]. MC2015, Nashville, Tennessee, 2015.
|
[4] |
李泽光, 王侃, 邓景康. 考虑源扰动效应的高阶蒙特卡罗微扰研究[J]. 核动力工程, 2014, 35(6): 6-10. LI Zeguang, WANG Kan, DENG Jingkang. Research on high-order perturbation calculation method with perturbed source effects[J]. Nuclear Power Engineering, 2014, 35(6): 6-10. DOI:10.13832/j.jnpe.2014.06.0006 |
[5] |
Yasunobu N, Takamasa M. Impact of perturbed fission source on the effective multiplication factor in Monte Carlo perturbation calculations[J]. Journal of Nuclear Science & Technology, 2005, 42(5): 428-441. DOI:10.3327/jnst.42.428 |
[6] |
Shim H J, Kim C H. Adjoint sensitivity and uncertainty analyses in Monte Carlo forward calculation[J]. Journal of Nuclear Science & Technology, 2011, 48(12): 1453-1461. DOI:10.1080/18811248.2011.9711838 |
[7] |
Wu Y. Multi-functional neutronics calculation methodology and program for nuclear design and radiation safety evaluation[J]. Fusion Science and Technology, 2018, 74: 1475162. DOI:10.1080/15361055.2018.1475162 |
[8] |
Wu Y, Chen Z, Hu L, et al. Identification of safety gaps for fusion demonstration reactors[J]. Nature Energy, 2016, 1: 16154. DOI:10.1038/nenergy.2016.154 |
[9] |
Wu Y, Song J, Zheng H, et al. CAD-based Monte Carlo program for integrated simulation of nuclear system SuperMC[J]. Annals of Nuclear Energy, 2015, 82: 161-168. DOI:10.1051/snamc/201406022 |
[10] |
Wu Y, FDS Team. CAD-based interface programs for fusion neutron transport simulation[J]. Fusion Engineering and Design, 2009, 84(7-11): 1987-1992. DOI:10.1016/j.fusengdes.2008.12.041 |
[11] |
吴宜灿, 宋婧, 胡丽琴, 等. 超级蒙特卡罗核计算仿真软件系统SuperMC[J]. 核科学与工程, 2016, 36(1): 62-71. WU Yican, SONG Jing, HU Liqin, et al. Super Monte Carlo simulation program for nuclear and radiation process:SuperMC[J]. Nuclear Science and Engineering, 2016, 36(1): 62-71. |
[12] |
Wu Y. Conceptual design of the China fusion power plant FDS-Ⅱ[J]. Fusion Engineering and Design, 2008, 83(10-12): 1. |
[13] |
Wu Y, Jiang J, Wang M, et al. A fusion-driven subcritical system concept based on viable technologies[J]. Nuclear Fusion, 2011, 51(10): 103036. DOI:10.1088/0029-5515/51/10/103036 |
[14] |
Wu Y, Bai Y, Song Y, et al. Development strategy and conceptual design of China Lead-based Research Reactor[J]. Annals of Nuclear Energy, 2016, 87: 511-516. DOI:10.1016/j.anucene.2015.08.015 |
[15] |
Wu Y. Design and R & D progress of China Lead-Based Reactor for ADS research facility[J]. Engineering, 2016, 2(1): 124-131. DOI:10.1016/J.ENG.2016.01.023 |
[16] |
Huang Q, Baluc N, Dai Y, et al. Recent progress of R & D activities on reduced activation ferritic/martensitic steels[J]. Journal of Nuclear Materials, 2013, 442(1-3): S2-S8. DOI:10.1016/j.jnucmat.2012.12.039 |
[17] |
Huang Q, Li C, Li Y, et al. Progress in development of China low activation martensitic steel for fusion application[J]. Journal of Nuclear Materials, 2007, s367-370(10): 142-146. DOI:10.1016/j.jnucmat.2007.03.153 |
[18] |
Huang Q. Status and improvement of CLAM for nuclear application[J]. Nuclear Fusion, 2017, 57: 086042. DOI:10.1088/1741-4326/aa763f |
[19] |
Wu Y. Design status and development strategy of China liquid lithium-lead blankets and related material technology[J]. Journal of Nuclear Materials, 2007, 367(4): 1410-1415. DOI:10.1016/j.jnucmat.2007.04.031 |
[20] |
Wu Y, Liu C, Song G, et al. Development of high intensity D-T fusion neutron generator (HINEG)[C]. European Physical Journal Web of Conferences, European Physical Journal Web of Conferences, 2017: 03006. DOI: 10.1051/epjconf/201715303006.
|
[21] |
吴宜灿, 刘超, 宋钢, 等. 强流氘氚聚变中子源HINEG设计研究[J]. 核科学与工程, 2016, 36(1): 77-83. WU Yican, LIU Chao, SONG Gang, et al. Design study of high intensity D-T fusion neutron generator HINGE[J]. Nuclear Science and Engineering, 2016, 36(1): 77-83. |
[22] |
谢仲生, 邓力. 中子输运理论数值计算方法[M]. 西安: 西北工业大学出版社, 2005, 303-345. XIE Zhongsheng, DENG Li. Neutron transport theory numerical method[M]. Xi'an: Northwestern Polytechnical University Press, 2005, 303-345. |
[23] |
Hurwitz H. Physical interpretation of the adjoint flux: iterated fission probability[M]. Naval Reactor Physics Handbook I, U. S. Atomic Energy Commission, 1964: 864-869.
|
[24] |
丘意书, 梁金刚, 王侃. 基于反复裂变几率法的敏感性分析初步研究[J]. 核动力工程, 2014, 35(S2): 83-86. QIU Yishu, LIANG Jingang, WANG Kan. Preliminary study of sensitivity analysis based on iterated fission probability method[J]. Nuclear Power Engineering, 2014, 35(S2): 83-86. DOI:10.13832/j.jnpe.2014.S2.0083 |
[25] |
刘勇, 曹良志, 吴宏春, 等. 基于经典微扰理论的特征值灵敏度和不确定性分析[J]. 原子能科学技术, 2015, 49(7): 1247-1253. LIU Yong, CAO Liangzhi, WU Hongchun, et al. Eigenvalue sensitivity and uncertainty analysis based on classical perturbation theory[J]. Atomic Energy Science and Technology, 2015, 49(7): 1247-1253. DOI:10.7538/yzk.2015.49.07.1247 |
[26] |
佘顶, 梁金刚, 吴高晨, 等. 堆用蒙卡程序RMC-Beta2. 0用户使用手册[CP/DK]. 北京: 清华 REAL, 2013: 18-19. SHE Ding, LIANG Jingang, WU Gaochen, et al. Monte Carlo code for reactor analysis RMC-Beta2. 0 user manual[CP/DK]. Beijing: Tsinghua University REAL, 2013: 18-19. |
[27] |
甘佺, 吴斌, 宋婧, 等. 基于参数的可视化裂变堆芯蒙特卡罗自动建模方法[J]. 核技术, 2016, 39(6): 060501. GAN Quan, WU Bin, SONG Jing, et al. Rapid parameter-based and visual Monte Carlo modeling method of fission reactor core[J]. Nuclear Techniques, 2016, 39(6): 060501. DOI:10.11889/j.0253-3219.2016.hjs.39.060501 |
[28] |
Gan Q, Wu B, Yu S, et al. CAD-based hierarchical geometry conversion method for modeling of fission reactor cores[J]. Annals of Nuclear Energy, 2016, 94: 369-375. DOI:10.1016/j.anucene.2016.03.013 |
[29] |
Pan X L, Wang J Q, Yuan R, et al. Biasing transition rate method based on direct MC simulation for probabilistic safety assessment[J]. Nuclear Science and Techniques, 2017, 28(7): 91. DOI:10.1007/s41365-017-0255-2 |