Leakage-corrected fast reactor assembly calculation with Monte-Carlo code and its validation methodology
-
摘要: 一种基于B1均匀化方程的泄漏修正模型在连续能量蒙特卡罗程序TRIPOLI4中得以实现并且用于制作少群截面参数。此蒙卡泄漏修正模型通过在连续能量的蒙卡模拟以及求解B1均匀化方程之间迭代,最终得到蒙卡模拟下的临界状态。通过此方法得到的少群截面参数较其他蒙卡以及确定论方法有两点显著优势:用于求解B1均匀化方程的少群常数是用通过临界状态的通量谱得到的;考虑了泄漏效应的蒙卡模拟可以更真实地反映组件计算时的能谱状态。为验证此泄漏修正模型,一个由连续能量的TRIPOLI4模拟而得到的数值临界实验被用于分析与比较。通过与其他蒙卡程序SERPENT以及确定论程序ECCO进行结果对比,可证明此B1泄漏修正方法能够给出更精确的用于堆芯计算的少群截面参数。Abstract: A leakage model based on B1 homogeneous equations has been implemented in continuous-energy Monte Carlo code TRIPOLI4. This leakage model algorithm iterates between the point-wise Monte Carlo simulation and a B1 homogeneous equation solver till reaching a final critical state in Monte Carlo simulation. The two advantages of our leakage model compared with the others are: we use critical flux spectrum to generate the multi-group constants for solving the B1 homogeneous equation; the leakage coefficients calculated are considered in point-wise Monte Carlo simulation. This leakage model is validated by a pre-designed numerical experiment simulated with continuous-energy TRIPOLI4 and the results obtained by this leakage model are proved to be more accurate by comparison with those from SERPENT leakage model and deterministic leakage model in ECCO code.
-
Key words:
- leakage model /
- critical buckling /
- cross section
-
蒙特卡罗(简称“蒙卡”)方法在核反应堆堆芯设计与分析中起着越来越重要的作用。鉴于此方法能够对堆芯几何结构进行精确描述并且能够使用连续能量截面的优势,基于该方法的程序经常被用于三维全堆输运计算并给出参考结果。蒙特卡罗算法的优点亦可以用到现有的两步法堆芯计算中(燃料组件计算+堆芯计算),而少群参数的制作是两步法计算中的关键。目前,针对蒙卡方法计算截面参数的研究仍然比较受限,大部分研究还是集中在适用于扩散程序的少群均匀化参数上。
在组件计算时,边界条件经常设为全反射或者周期性边界条件。热中子的平均自由程在1~2 cm之间,基本可满足全反射边界条件的假设。相反,快中子的平均自由程远超过单个燃料组件的尺寸。相邻组件的联系在快中子之间更为紧密,故而全反射边界条件对其并不适用。因此,在考虑快中子堆组件计算时需要考虑泄漏效应。考虑泄漏效应的少群常数制作可以通过在计算中嵌入泄漏修正模型来实现,此方法期望通过理论模拟在即使不清楚反应堆堆芯真实运行条件下也可以构造出组件临界状态。基于蒙卡方法的组件泄漏模型大致分为三种。方法一来自于确定论中常用的基模修正法,该方法首次用于韩国KAERI的蒙卡程序McCARD中[1],芬兰VTT实验室开发的蒙卡燃耗程序SERPENT[2]也使用此方法做泄漏修正。第二种泄漏修正方法来自于韩国KAIST实验室,这种方法的主要思想是通过改变边界反照率来搜索组件临界状态[3],此方法已经在MCNP5中实现并且与MONTEBURNS程序耦合做燃耗计算。第三种泄漏模型由日本京都大学的Yamamoto教授提出,通过给中子权重赋予复数形式,新颖地将中子泄漏项与其传播方向结合起来[4]。与前两种方法相比,第三种方法能够给出每个中子的泄漏概率,更合理地应用于蒙卡程序的模拟过程中。尽管此方法占据了理论上的优势,但它的几何适用范围受限,为了消除裂变源项中的虚数部分,该方法只适用于对称的几何结构上。
本文介绍了法国蒙卡程序TRIPOLI4中嵌入的泄漏修正模型。此泄漏模型的开发主要用于少群截面参数的制作[5]。目前,TRIPOLI4程序中选择的泄漏修正模型为基模修正方法。另外,基于蒙卡方法制得的少群截面参数的验证也是本文的一个重点内容。本文设计了一套自然泄漏的数值计算模型用以验证TRIPOLI4程序中少群截面参数的制作功能。此外,该程序产生的截面参数也将与已经相对成熟的快堆组件程序ECCO[6-7]和另一套蒙卡程序SERPENT进行对比。最后,这些少群截面参数还将被用于堆芯计算中并分析比较其计算结果。
1. 模型与算法
主要叙述基模修正模型及其在程序中使用的具体算法,包括详细的针对基模方程的物理近似,以及对此模型在程序中的流程嵌入和代码实现。
1.1 B1均匀化方程
输运方程在均匀无穷介质结构下的解即通常所称的基态模式。组件计算一般在全反射或者周期性边界条件下进行。因而,泄漏模型必须要加进燃料组件计算中才能强制有效增殖系数Keff=1。基模的解通常由两部分组成,第一部分是基态通解φ(r, E, Ω)对于任意位置的组件通解保持不变;第二部分f(r)被用于堆芯中的通量幅值形状,并且,此形状函数f(r)是拉普拉斯方程(1)的解。
∇2f(r)+B2f(r)=0 (1) 式中:B2为堆芯的几何曲率。该方程对应的通解为f(r)=αeiB·r。对于每个独立组件,基态通解已经不再随空间变化。因此,组件内的角通量密度可写为
φ(r,E,Ω)=φ(E,Ω)eiB⋅r (2) 式中:B代表曲率向量,其值满足B·B=B2。
把式(2)变量分离后的通量密度表达式代入连续能量的Boltzmann方程中,并采用线性各向异性近似假设,即可得到均匀化的B1基模方程
[iΩ⋅B+Σt(E)]φ(r,E)=∫∞014π{Σs0(E′→E)φ(E′)+3Σs1(E′→E)Ω⋅J(E)}dE′+14πKeff∫∞0χ(E′→E)νΣf(E′)φ(E)dE′ (3) 式中:Σt,Σs0,Σs1分别代表总截面、零阶散射截面、一阶散射截面;ν为中子裂变增殖系数;χ为中子裂变谱;J (E)代表基态中子流,表示为
J(E)=∫4πΩφ(E,Ω)d2Ω (4) 为了消除曲率和中子流的方向性,假设J和B具有相同的方向,故而有J(E)·B =J(E)B。
通过对方程(3)进行角度积分以及乘以1iΩ⋅B+∑t(E)后再进行角度积分,可分别得到下面两组与标通量密度相关的方程
iBJ(E)+Σt(E)φ(E)=∫∞0Σs0(E′→E)φ(E′)dE′+1Keff∫∞0χ(E′→E)νΣf(E′)φ(E′)dE′ (5) iJ(E)B=1γ[B,Σt(E)]Σt(E){13φ(E)+∫∞0Σs1(E′→E)iJ(E′)B dE′} (6) 式中:γ[B, Σt(E)]为新引入的变量,其定义如下
γ[B,Σt(E)]={13B∑t(E)arctanB∑t(E)1−B∑t(E)arctanB∑t(E),B2<0B23∑t(E)ln∑t(E)+Im(B)∑t(E)−Im(B)2Im(B)−∑t(E)ln∑t(E)+Im(B)∑t(E)−Im(B),B2>01+415(B∑t(E))2−12175(B∑t(B))4+⋯,B2>0 (7) 需要特别指出的是,当曲率接近为零时,iJ(E)/B应该一直保持为实数并且有限。实际上,根据所在介质的不同,J(E)可以为虚数或为实数。当堆芯处于超临界状态时,J(E)为虚数并且B为正数;反之,当堆芯处于次临界状态时,J(E)为实数而且B为虚数。
于是,通过菲克定律,可以得到扩散系数为
d(B,E)=1BiJ(E)φ(E) (8) 将式(8)分别代入方程(5)与方程(6)中,可得到如下变形后的B1方程组
[Σt(E)+d(B,E)B2]φ(E)=∫∞0Σs0φ(E′)dE′+1Keff∫∞0χ(E′→E)νΣf(E′)φ(E′)dE′ (9) d(B,E)=1γ[B,Σt(E)]Σt(E){13+∫∞0Σs1(E′→E)d(B,E′)φ(E′)φ(E)dE′} (10) 有了上述连续能量的B1方程组之后,可以很容易通过并群得到其少群形式下的基态方程组
Σt,gφg+dgB2φg=∑Gg′=1∑g′→gs0φ′g+1Keff ∑Gg′=1χ(g′→g)νΣf,g′φ′g (11) dg=1∑t,gγg{13+∑Gg′=1∑g′→gs1φ′gφg d′g} (12) 式中:φg为g群内通量密度,φg=∫E∈gφ(E)dE;dg为g群内扩散系数,dg=1φg∫E∈gd(B,E)φ(E)dE;γg为系数,可简化为近似表达形式γg≈γ(B, Σt, g);Σi, g为一维少群截面,例如总截面、吸收截面、裂变截面等,∑i,g=∫E∈g∑i(E)φ(E)dE∫E∈gφ(E)dE;∑g′→gs0为多群散射截面的零阶项,∑g′→gs0=(Σt,g′−Σa,g′)∗P(g′→g),其中,P(g′→g)为群间散射概率;∑g′→gs1为多群散射截面的一阶项,∑g′→gs1=∑g′→gs0ˉμg′→g,其中,ˉμg′→g代表从g′到g群的平均散射角余弦。
值得指出的是,少群截面参数是由连续能量的TRIPOLI4模拟计算提供,而方程组(11)~(12)中未知量为dg,φg和γg,其具体求解方法将在下一节详细介绍。
1.2 B1算法流程
本节给出TRIPOLI4程序中求解B1均匀化方程的计算流程,其求解过程由内迭代和外迭代两部分交替完成。内迭代过程即在固定B2值下求出收敛后的通量密度以及Keff值;外迭代过程即通过方程(12)变化B2值来使得Keff=1。
B2n+1=B2n+(1Kexp eff −1Kneff )¯ν∑nfˉdn (13) 式中:n为内迭代小标数;Keffexp为用户设定的Keff=1值;νΣfn为第n次内迭代后,由通量密度权重得到的平均裂变产生截面;dn为第n次内迭代后,由通量密度权重得到的扩散系数。
图 1给出了完整的迭代流程。图中左边部分代表连续能量截面的蒙卡模拟过程,右边部分代表求解B1均匀化方程的过程。左右两部分通过两条线连接。其一为收敛后的少群截面参数,将作为已知量被用于求解B1均匀化方程中;其二为由求解B1方程后得到的泄漏截面dgB2,此泄漏系数将被代入蒙卡的连续截面模拟过程中。由于该泄漏修正模型涉及到基理完全不同的蒙卡程序和确定论程序的耦合,因此在收敛准则上需要注意匹配。此处主要涉及到三种收敛准则。第一,用户根据需求可以给出蒙卡模拟过程中的总粒子数(批数量×每批粒子数)以确保每个能群下计算得到的截面参数都是收敛值;第二,求解B1均匀化方程时对Keff的收敛标准也控制得比较严格,一般在|Keff-1| < 1×10-5;第三,最终的蒙卡模拟得到的Keff值与用户设定的期望值Keffexp可保持在如下收敛范围内:|Keff-Keffexp| < 100×10-5+3σ。公式中σ为蒙卡模拟时给出的自身不确定性。
2. 验证
介绍了TRIPOLI4程序中泄漏修正模块的验证工作。首先,通过蒙卡连续能量模拟得到的一个一维临界几何结构被设计出来用于给出各项参考量,例如临界通量密度谱、反应率以及少群截面参数等。与TRIPOLI4程序中泄漏模块进行对比的有蒙卡程序SERPENT以及已经较成熟的用于快堆截面制作的确定论程序ECCO。本验证工作选择了快堆组件Zero Power Plutonium Reactor Double Column Fuel(ZPPRDCF)的材料成分作为一维临界几何模型来展开计算以及分析对比工作。此组件模型的Kinf高达1.67,也就意味着此临界模型的中子泄漏率达40%。这种高泄漏结构也就对泄漏模型的处理要求更加严格。另外,为了使各种程序的计算结果有可比性,各软件统一使用JEFF3.1.1作为输入截面库。下面将分别给出计算结果参数的对比分析。
2.1 临界几何曲率B2
现今包含泄漏模型的组件计算程序给出几何曲率已经不再困难,而如何找到一套可靠的参考值却是值得推敲的问题。其实,该问题可以转移到如何确定方程(1)中空间分布函数f(r)。多数情况下,f(r)被认为是矢量通量密度分布φ(r)。分布函数f(r)实际上是两个不同通量密度的内积函数
f(r)=∫∞0φ+(E)φ(r,E)dE (14) 式中:φ+ (E)为基态模式下的共轭通量密度;φ(r, E)为堆芯计算下的直接通量密度。此两者通量密度皆由TRIPOLI4连续蒙卡模拟计算给出。因而,在已知f(r)轴向分布的情况下,可以拟合出几何曲率B2。
表 1给出了各种不同程序计算得到的几何曲率值以及他们的比较,其中,由通量密度内积拟合出的曲率值作为参考结果。可看出,TRIPOLI4的计算结果与参考值最为接近,相差1.09%,而其他两套程序的计算结果与参考值相比偏差稍大。此结果可归功于TRIPOLI4的泄漏修正模块中通过连续截面的蒙卡模拟过程与B1方程求解的交叉迭代直至达到真正临界状态,而这样的计算模型在其他两套程序中是不具备的。
表 1 临界曲率比较结果Table 1. Comparison of critical buckling valuesB2/(10-3cm-2) relative diff/% numerical curve 2.854 - TRIPOLI4 2.823 -1.09 SERPENT 2.807 -1.65 ECCO 2.747 -3.75 2.2 归一化临界通量密度谱
图 2给出了多种包含泄漏修正模型程序的归一化临界通量密度计算结果。其中的实线(Numerical Exp)代表由连续能量截面的TRIPOLI4模拟给出的参考结果。图 2(a)给出了各程序计算得到的归一化通量密度值,其中归一化原则为∑Gg=1φg=1,能群总数G可取33。可以看出,这三套程序中的泄漏修正模型都能够给出与参考值趋势一致的结论。图中灰色区域为重要能群区,故而将此区域数据单独描在图 2(b)中。此外,三套程序与参考值的相对误差也在图 2(c)中给出。可以看出,在高能区部分,TRIPOLI4,SERPENT以及ECCO程序中的泄漏修正模型均低估了临界通量密度水平。相反,在较低能量区域,三套程序中的泄漏修正模型却高估了通量密度水平。
2.3 泄漏率
本节给出了各程序中泄漏修正模型模拟得到的泄漏率(dgB2φg)及其与参考值的对比。参考值依然来自于TRIPOLI4连续能量截面模拟得到的中子泄漏率。图 3中描述了各软件计算得到的泄漏率与参考值的偏差。可以看到,这三套程序模拟的泄漏率与参考值相比在每个能群内偏差都控制在500×10-5以内。此外,这三套程序中都存在高能群区与低能群区的互补现象,这一现象主要由基态模型的近似造成。在所使用的基态模型中,认为中子谱不随组件在堆芯中的位置变化。实际上,在一个真实堆芯中,堆芯外缘的谱要比堆芯中部的谱稍硬,因而在堆芯外缘处,快中子泄漏的机率偏小,这也就解释了图 1中在2 MeV处中子通量密度估算值偏高的现象。
3. 结论
本文介绍了法国蒙卡程序TRIPOLI4中嵌入的泄漏修正模型以及利用此模型制作适用于快堆的少群截面参数的工作。此泄漏模型利用通过蒙卡产生的少群参数作为B1均匀化方程的已知量从而求得相应的少群泄漏参数。该泄漏模型的优点在于TRIPOLI4程序能够利用以上解得的泄漏系数在连续蒙卡模拟过程中直至迭代成功到真正临界状态。本文也给出了关于此泄漏模型的基础验证工作。通过对比几何曲率、临界中子谱、泄漏率等参数,可以看到TRIPOLI4程序中的泄漏修正模型能够给出参考值相符合的结果,甚至比同类程序(SERPENT,ECCO)的计算精度更高。
-
表 1 临界曲率比较结果
Table 1. Comparison of critical buckling values
B2/(10-3cm-2) relative diff/% numerical curve 2.854 - TRIPOLI4 2.823 -1.09 SERPENT 2.807 -1.65 ECCO 2.747 -3.75 -
[1] Park H J, Shim H J, Joo H G. Generation of few group diffusion theory constants by Monte Carlo code McCARD[J]. Nuclear Science and Engineering, 2012, 172(1): 66-77. [2] Leppänen J. SERPENT: a continuous-energy Monte Carlo reactor physics burnup calculation code[R]. Finland: VTT Technical Research Centre of Finland. 2013. [3] Cho N Z, Yun S H, Lee J. Generation of homogenized nodal parameters by Monte Carlo method with non-zero leakage spectra in global-local iteration framework[J]. Transactions of the American Nuclear Society, 2009, 101: 707-710. [4] Yamamoto T. Monte Carlo algorithm for buckling search and neutron leakage-corrected calculations[J]. Annals of Nuclear Energy, 2012, 47: 14-20. doi: 10.1016/j.anucene.2012.04.017 [5] Cai Li, Pénéliau Y, Diop C M, et al. P1 adaptation of TRIPOLI-4 code for the use of 3D realistic core multigroup cross section generation[C]//Joint International Conference on Supercomputing in Nuclear Application and Monte Carlo. 2013. [6] Rimpault G. Algorithmic features of the ECCO cell code for treating heterogeneous fast reactor assemblies[C]//International Conference on Mathematics and Computations, Reactor Physics, and Environmental Analyses. 1995. [7] Grimstone M J, Tullett J D, Rimpault G. Accurate treatments of fast reactor fuel assembly heterogeneity with the ECCO cell code[C]//International Conference on the Physics of Reactors: Operation, Design and Computation-PHYSOR. 1990. 期刊类型引用(2)
1. 吴宏春,杨红义,郑友琦,曹良志,杜夏楠,杨勇,刘一哲,胡赟. 快中子反应堆堆芯物理分析方法的研究现状与发展建议. 原子能科学技术. 2024(03): 513-527+507 . 百度学术
2. 李金洲,张滕飞,何东豪,潘清泉,刘晓晶. MORPHY程序在铅冷快堆中的应用. 强激光与粒子束. 2024(06): 142-153 . 本站查看
其他类型引用(0)
-