点探测器敏感性关联抽样计算方法.pdf

第58 卷第3期2024年3月原子能科学技术Atomic Energy Science and TechnologyVol.58,No.3Mar.2024点探测器敏感性关联抽样计算方法李瑞1-2，付元光1-2，邓力，许海波2（1.中物院高性能数值模拟软件中心，北京10 0 0 8 8;2.北京应用物理与计算数学研究所，北京10 0 0 94)摘要：为避免点探测器敏感性计算中随机统计涨落的影响，研究了在蒙特卡罗关联抽样算法框架下，点探测器敏感性的计算方法。在算法中考虑了输运过程偏倚以及存在次级粒子的情况，给出了输运扰动权乘子与计数扰动权乘子的更新方案。针对点探测器敏感性中次级粒子的模拟特点，引入了单因素扰动系统，并基于粒子扩展属性实现了对单因素扰动系统微扰的模拟，从而避免了对次级粒子历史进行回溯。利用脉冲球形实验装置，测试了散射截面扰动下点探测器的敏感性计算功能，并与微分算符方法的计算结果进行了比较，以验证其计算精度。利用NUREG/CR-6115屏蔽基准题测试了权重技巧下算法的适用性。测试表明，算法与微分算符方法计算精度相当，并适用于屏蔽问题点探测器敏感性模拟。关键词：点探测器敏感性；关联抽样；蒙特卡罗中图分类号：TL99;0571.53；0 2 11.9doi:10.7538/yzk.2023.youxian.0428Correlated Sampling Method for Point Detector Sensitivity Calculation(1.CAEP Software Center for High Performance Numerical Simulation,Beijing 100088,China;2.Institute of Applied Physics and Computational Mathematics,Beijing 100094,China)Abstract:In reactor radiation shielding design,the shielding scheme often needs to beoptimized by changing shielding body size,material density and isotope atomic frac-tions.The knowledge of radiation dose rates difference can help a lot.Besides,theuncertainty of the nuclear data itself affects the accuracy of the radiation shielding simu-lating results,which also needs to be considered in the evaluation of the shieldingscheme.When the difference of radiation dose rate is small,perturbation method can beused.Perturbation is especially useful in Monte Carlo calculation in which statistic errormasks the inherent difference between two schemes.In Monte Carlo radiation shieldingcalculation,point detector tally has great efficiency by using next event estimator.Inorder to calculate the sensitivity of point detector tally,correlated sampling algorithmwas introduced in JMCT software.Random walk was assumed the same for both pertur-bated and original systems.Instead,dose difference was expressed only in particleweight.Perturbation caused by material and nuclear data propagated in the entirerandom walk.General particle weight correction was calculated in colliding and flying文献标志码：ALI Ruil-2,FU Yuanguangl-2,DENG Li?,XU Haibo?文章编号：10 0 0-6 9 31(2 0 2 4)0 3-0 6 38-0 7收稿日期：2 0 2 3-0 6-15；修回日期：2 0 2 3-0 9-13基金项目：国家自然科学基金（12 2 7 50 30）第3期李瑞等：点探测器敏感性关联抽样计算方法events,and virtual particle weight correction was calculated in source sampling andcolliding events.Weight correction was deduced for every process considering transpor-tation biasing,often used in deep penetrating problem simulation.System with singleperturbation was introduced to encapsulate weight correction calculation.Within oneperturbated system,only perturbation in one energy region,isotope and reaction typewas considered.Meanwhile,the perturbated system had the same number of tallies withoriginal system,and the sensitivity tallies could be built by cloning from original sys-tems point detector tallies.This feature avoids a lot of underlying mistakes in memoryallocation and tally initialization.Multiply perturbation can be considered in one MonteCarlo simulation without interference with each other.Weight correction of descendantis managed in form of particle extended property so that no history retrieving is needed.This makes the algorithm compatible with variance reduction techniques such as weightwindow or cell importance.Validation was performed on Livermore pulsed sphereexperiment and NURGE/CR-6115 benchmark.Good agreement is found in time-of-flight leakage spectra in 2.9 mean free paths carbon sphere with results published by R.L.Perel which uses differential operator method.For the perturbation of inelastic scat-tering cross section in specified energy range,difference below 5%is achieved in bothsource direct leakage peak and single scattering leakage peak compared with referenceresults.In NUREG/CR-6115 benchmark,consistent adjoint driven importance samplingwas used to reduce the statistical variance.Comparing with the result from analogyMonte Carlo,qualitative agreement is found.However,time consuming in correlatedsampling method is 4.1 CPU hours compared with 4 286.2 CPU hours in analogy MonteCarlo method,which makes the correlated sampling method is applicable in shieldingdesign optimization.Key words:point detector sensitivity;correlated sampling;Monte Carlo在反应堆辐射屏蔽设计中，常需要通过改变屏蔽体尺寸、材料密度与组分实现屏蔽方案的优化。在此过程中，掌握这些因素导致的辐射剂量率变化规律可有效提高屏蔽方案优化效率。另外，核数据本身的不确定性会对辐射屏蔽模拟结果的准确性造成影响，在屏蔽方案的评估中也需要加以考虑。当不同屏蔽方案的差异较小时，上述计算辐射剂量率变化的问题可视为微扰问题。反应堆辐射屏蔽模拟一般采用离散纵标(SN）方法与蒙特卡罗（MC)方法。其中SN方法属于确定论方法，在计算微扰问题时可直接采用两次独立计算并做差的方式进行。但是，由于SN方法对于复杂几何处理能力有限，并存在射线效应等问题，其计算精度一般低于MC方法。与S方法相比,MC方法可处理比较复杂的几何体，采用的连续能量核反应截面639也更适用于核数据扰动分析。但MC方法的计算量巨大，特别是对于辐射屏蔽这一类深穿透问题，需要使用降方差技巧才能得到统计收敛结果。然而,在微扰问题计算中，由于辐射剂量率本身的改变就很小，如果采用直接MC模拟方法计算，则对于其统计收敛性要求更加苛刻,导致时间开销进一步增加。因此，为排除随机涨落的影响,MC方法一般在与参考系统相同的随机序列下完成微扰系统的模拟，常见的方法包括关联抽样1与微分算符2。近年来，采用伴随通量作为重要性的降方差技巧在MC辐射屏蔽模拟中被广泛研究，如针对单一目标的CADIS算法 3 以及多目标的FW-CADIS算法 4。这些算法的实质是基于粒子权重的赌/分裂以及偏倚抽样的组合。另一方面，在MC辐射屏蔽模拟中，常采用点探测器计数方法提高计数率。640MCNP软件L51基于微分算符法，提供了体平均通量及其衍生计数的敏感性计算功能。Perel等 6 在此基础上研究了点探测器计数的敏感性计算方法。与体平均通量这一类径迹长度估计方法不同，点探测器计数采用了下次事件估计方法7，需要额外考虑虚粒子的模拟过程。在Perel等 6 1展示的算法中，需多次遍历粒子的径迹历史贡献。在采用权重赌/分裂技巧的屏蔽问题模拟中，次级粒子的数目很多且平均粒子随机游动链较长，这种算法的性能开销会很大。本工作基于关联抽样算法架构，利用单因素扰动系统将系统扰动因素局限在可通过一个扰动权重描述的简单情形，在此基础上以随输运过程实时更新的方式代替频繁的粒子历史遍历。另外，将扰动权作为粒子扩展属性也可避免对次级粒子径迹的回溯需求，从而降低算法的实现难度。1算法介绍关联抽样算法的重点在于采用与原系统相同的随机序列模拟扰动系统的输运与计数统计过程。在形式上，扰动系统粒子随机游动过程与原系统是完全相同的，统计差异仅体现在粒子权重上。本工作基于三维粒子输运模拟软件JMCT开展算法研究 1。对于辐射屏蔽问题，点探测器计数估计量 6 1可表示为：1j=1式中：EstTa为点探测器计数估计量；Ta为点探测器计数；N为模拟的粒子历史总数；i为粒子历史索引；J；为第i个粒子历史的径迹段总数；i为粒子历史中累计模拟的径迹段数目；k为每个径迹段的索引;lin(j)为从源抽样到第j径迹段随机序列中所有径迹段索引的集合；;为粒子经过个累计径迹段之后触发的点探测器计数响应，对应的能量取决于飞人探测器的虚粒子能量Ea；Pk（E）为每个真实粒子随机游动径迹段产生的概率，相应径迹段内粒子的能量表示为E。式（1)中仅在粒子随机游动链中的碰撞事件或源抽样事件中触发。具体来说，是在确原子能科学技术第58 卷定粒子的反应类型之后，且在反应模拟之前，即进行下次事件估计7 时。；可表达为式（2)所示形式。由于点探测器计数的触发特点，在径迹段的定义中，不单独统计粒子飞行到几何边界的事件，而是将其作为某个径迹段的一部分。8(2dR&-R7;(Ea)=Rd-R2Ra-RP,(2a,Ea)e/o式中：为克罗内克函数，表示仅当虚粒子飞行方向朝向点探测器时触发计数；2 a为虚粒子飞行方向矢量；Ra为点探测器的位置矢量；R为随机游动链中粒子的位置矢量；P。为粒子发生反应时在点探测器方向发射粒子的概率，即指向概率；Z.为总反应截面；s为虚粒子到达探测器过程中经历的径迹段长度。在触发计数时从碰撞点创建虚粒子，使之飞人点探测器，因此虚粒子的能量及飞行方向与随机游动链中碰撞后的粒子是不同的。为显示这一区别，式（2)中以E。表示虚粒子的能量，2 a表示虚粒子的飞行方向矢量，（2 aRd-R一1)描述了虚粒子飞行方向与点探测TR:-RT器位置Ra和碰撞点R之间满足的约束条件。P（2 d,Ea)为碰撞点处产生的飞入点探测器的虚粒子概率，需要根据粒子碰撞的次级角分布计算，其中次级角分布的扰动不在本工作研究kelin(i)范围之内。指数项表示虚粒子从碰撞点飞入点(1)探测器过程中的衰减因子，其数值取决于Ea能量虚粒子在飞行径迹上的，与飞行距离。式(1)中的P（E)为粒子随机游动链中径迹段的产生概率，按照Perel 等 6 1 对于径迹段指标的约定，每个径迹段以碰撞或源抽样事件开始，因此点探测器计数是在径迹段开始处触发的。由于虚粒子与实际碰撞或源抽样过程是独立的，因此影响第k个径迹段点探测器计数的所有径迹段kElin(i)并不包括k径迹段本身 6。因此，对于首个点探测器计数，仅依赖径迹段0，为保证形式上的一致性，将径迹段指标定义进行简单扩展，扩展之后的径迹段概率为：Rd-R1)(2)第3期李瑞等：点探测器敏感性关联抽样计算方法1k=0T(R,Ei,2i;R)Q(R,Ei,2)k=1Pk=T(R,Ek,2k;R-1):C(Rk-1,Ek,2k;Ek-1,2k-1)Q为辐射屏蔽计算中的外源分布，源粒子的位置属性指标为0，能量与方向指标为1。由于粒子从R。飞行到R1过程中,其能量与方向是不变的,因此E1与2 1分别对应了第2 个径迹段开始之前粒子的能量与飞行方向。对于k1的情况，位置属性指标与径迹段的关系以此类推。T为输运核，表示粒子飞行至下一次碰撞的概率，其形式为：T(R,Ek,2;R-1)=C，为各类反应的碰撞核，如式(5)所示。C,(Rk-1,Ek,2k;E-1,2k-1)=V(E-1)Z,(E,2;Ek-1,2k-1)Z.(E-1,52k-1)式中：Z为宏观总截面；Z，为宏观反应截面；V为反应发生后产生的次级粒子数目，主要用于中子诱发裂变、（n，n)以及中子产光等反应过程。在MC辐射屏蔽模拟计算中,使用的权重技巧包括源抽样偏倚、权窗、几何重要性等，这些技巧都可以描述为对粒子随机游动过程中输运核与碰撞核的偏倚。以P录表示偏倚后粒子径迹段的概率，可将式(1)改写为：NExLET.1-22(E.).1j-1I(P(E)e(E2)-Pik(E)kelinGi)=1j-1式中，w为第k径迹段对应的粒子权重，根据统计无偏性的要求，其定义如式（7)所示。由于赌/分裂等技巧的存在，粒子径迹段的数目也会发生改变，采用lin(j)表示。W=kelinG)Pi(Ek)641T(R,Ek,2;Rk-1)C,(Re-1,E,2k;Ek-1,2k-1)T(Re,Eb,2;R-1)C(R-1,Ek,2;E-1,2-1)kelin(j)T(R1,E1,21;R。)Q(R,E,2)k1T(R1,Ei,21;Ro)Q(Ro,Ei,2)(3)式（6)所示为权重技巧下的点探测器计数形式。当系统存在扰动时，基于关联抽样计算点探测器计数对应变化量（T.）)的过程为：(T a)=Es t*T a-Es t T a =NJ7;(Ea)II P(Ek)kelin(i)7;(Ea)IIPi(E)=1 j=1KElin(i)W,(Ea)II Pu(Er)=kelin(i)NJ=1j=1式中，上标*表示扰动系统相关项。(4)式（8)对应的随机游动过程与式（6）类似，唯一的区别就是引入了扰动权乘子W：7(Ea)II Pi(E.)(5)kelin(j)W=7(Ea)II Pa(E)kelin(i)wdwi-1式中：w为仅由于虚粒子计数部分导致的计数扰动权乘子,其形式如式（10)所示；w为由于随机游动过程导致的输运扰动权乘子，其形式如式(11)所示。(Ea)7;(Ea)PA(EL)W=Pi(Ek)kelin(i)根据式（8）（11)可构建点探测器扰动量的关联抽样算法。仅需要在原系统随机游动模拟过程中更新并统计扰动权对应的计数即可，算法流程如图1所示。(6)其中，W仅对应虚粒子的扰动贡献，在触kelin(j)发点探测器计数的同时更新即可。与之不同，W包含了从源粒子到虚粒子产生之前所有径迹段的扰动贡献，需要在输运过程中同步更新。简洁起见，图1中未画出次级粒子产生情况下wi的更新流程。根据JMCT软件次级粒子的Pi(E,)处理流程，本工作将W作为粒子扩展属性进行管理。当产生1个以上次级粒子时，将当前(7)kelin(j)一1=(9)(10)(11)642粒子对应的与次级粒子一同存入粒子库。相应地，当对存库粒子进行模拟时，恢复数值即可。值得注意的是，尽管裂变截面、(n,n)截面或裂变平均释放中子数等数据的扰动会影响次级粒子数目，但根据式（5)碰撞核的定义，这一扰动也会在输运扰动权乘子中体现，不需要进行额外处理。更新源抽样输运扰动权乘子W*Q(Ro,E1,2,)0(Ro,E,2.)图1点探测器扰动关联抽样计算流程Fig.1Flowchart of point detector perturbationbased on correlated sampling另一方面，为完成点探测器扰动的模拟，对于每一个扰动因素，都需要建立与原系统结构相同的统计计数。结构相同是指点探测器的数目、能量区间划分、时间区间划分都相同。具有相同结构的扰动系统计数可通过对原系统计数进行克隆得到。所有扰动系统克隆计数的总内存使用量M如式（12）所示。原子能科学技术NdM=N,ZMi一1式中：N为扰动系统总数；Na为点探测器计数数目；M.为每个点探测器计数的内存使用量。JMCT中允许在一次模拟中同时分析多个扰动因素，这些扰动因素包括核素截面与材料密度的变化，以及针对扰动因素的能量区间划分。在本工作中,将某一能量区间内具有单一开始核素与反应道扰动或材料密度扰动的系统定义递增源粒子索引i-+1为单因素扰动系统。因此，一次MC模拟过程中,可存在多个单因素扰动系统。扰动系统对初始化径迹段索引k=1初始化径迹段概率Po-1初始化输运扰动权乘子！触发点探测器计数碰撞事件？是更新碰撞输运扰动权乘子C(Rk-1,Ek,2k;E-1,2k-1)W*C(R-1,Ex,2k;Ex-1,2k-1)计算点探测器虚粒子扰动权乘子元(Ea)累计点探测器敏感性计数(wawi-1)wki(Ea)递增径迹段索引k-k+1粒子飞行，更新输运扰动权乘子T(R,Ex,2k;RR-1)W*T(Rk,Ex,2k;R-1)径迹结束？是否完成所有粒子模拟？工是点探测器敏感性计数归一化与输出结束第58 卷(12)应的计数内存与内部状态是独占的，有效地避免了不同扰动因素在计算流程上相互干扰。对于每个单因素扰动系统，可通过扩展固定数目为1的粒子属性即可完成输运扰动权乘子的管理。此外,JMCT允许定义扰动所在几何体与材料，从而限制扰动的空间范围。2计算结果利用劳伦斯利弗莫尔国家实验室（LaW-rence Livermore National Laboratory,LLNL)的脉冲球形实验装置对算法进行正确性检验。该实验探测D-T脉冲中子在半径为2 0.96 cm的石墨球中输运后的时间泄漏谱，装置如图2所示。其中D-T中子源发生器通过孔道经一方向放置于=一0.131cm处。D-T中子能否量介于13.2 0 15.11MeV区间内，具体分布依赖于出射角度。实验中，NE-213探测器与一轴呈30 角，距离石墨球心7 7 6 cm。文献 9给出了该模型相应的MCNP程序 5输入文件。X图2 脉冲球形实验装置模型Fig.2Model of pulsed sphere experiment device第3期李瑞等：点探测器敏感性关联抽样计算方法JMCT首先基于该模型进行了泄漏谱的计算，作为对比的实验值来自文献 10。JMCT整体计算结果与Perel等 6 工作符合良好，如图3所示。-1.0-1.5F(_Su/率源+1)8 I-2.0-2.5F-3.0F-3.5-4.051o0150200250300350400时间/ns图3脉冲球形实验装置时间泄漏谱计算结果Fig.3Calculation result of time leakage spectrumfor pulsed sphere experiment device在此基础上，对14.9 18 15.2 9 6 MeV区间内的（n，n）反应截面引人十1%的扰动，计算时间谱的变化。由于该实验装置中D-T中子在石墨球中平均运行2.9 个自由程，因此扰动引起的时间泄漏谱变化主要集中在0 2 次散射对应的时间区间上。JMCT计算结果如图4所示。其中由于(n,ni)反应截面的增加导致扰动系统的总截面变大，进而造成点探测器直穿项计数降低。另外，散射截面的贡献会由于散射截面的增加而变大,其中一次散射的贡献最为显著。将JMCT点探测器直穿项和一次散射峰的变化与Perel等 6 工作中微分算符计211-Su,s-01/率:F骤10-1-2-3卜-4-5-61001502000250300350400时间/ns图4脉冲球形实验装置时间泄漏谱扰动计算结果Fig.4Calculation result of time leakage spectrumperturbation for pulsed sphere experiment device643算结果进行定量对比，相对偏差均在5%以内，如表1所列。表1时间泄漏谱扰动结果对比Table 1Comparison of time leakage71.0spectrum perturbation一JMCT计算结果x实验值10.8统计标准差10.60.40.20.071.0一JMCT计算结果.统计标准差0.80.60.40.20.0泄漏率/10-5nm-1扰动计数JMCT计算结果参考结果 6 偏差/%直穿6.675一次散射1.020针对反应堆辐射屏蔽问题，采用NUREG/CR-6115屏蔽基准题 11I验证使用降方差技巧的点探测器敏感性计算结果。通过在压力容器材料中引入十1%的核子密度扰动，计算压力容器外中子能谱的扰动。为保证计算收敛，该模拟采用了CADIS降方差方法 31,其中伴随通量来自JSNT12计算。在4X108粒子历史规模下,JMCT压力容器外表面中子能谱与文献11中的DORT结果对比如图5所示。1011一DORT参考结果1010+JMCT计算结果10%10810710610-7图5压力容器外表面中子能谱计算结果Fig.5Calculation result of neutron energy spectrumon outer surface of pressure vessel这里采用直接MC模拟作为参考对比，关联抽样的计算结果如图6 所示。结果表明，两者在趋势上是符合的，即12 MeV能量区间的中子通量下降，而1MeV能量以下的中子通量增加。导致这一现象的原因主要有2 个：一方面是由于压力容器材料对快中子的慢化作用导致高能中子慢化为低能中子；另一方面,压力容器内的Fe等核素对快中子的辐射俘获较低能中子更加显著。在该算例中，为避免两次直相对7.0004.641.0002.0010-510-3能量/MeV10-1101644接MC模拟的统计误差掩盖真实扰动，要求各能群统计标准差均低于0.0 2,两次MC计算共耗时42 8 6.2 CPU小时，关联抽样计算耗时4.1CPU小时。201510100246810 12 14 16能量/MeV图6 压力容器外表面中子能谱扰动计算结果Fig.6 Calculation result of neutron energy spectrumperturbation on outer surface of pressure vessel3结论脉冲球形实验装置计算结果表明，基于关联抽样的点探测器敏感性计算方法可达到与Perel等 6 工作中微分算符方法相当的计算精度，说明了算法对输运与计数过程的微扰描述是正确的。在NUREG/CR-6115屏蔽基准题模拟中,CADIS降方差技巧的使用一方面对源抽样与输运过程进行了偏倚，另一方面也会在权重轮盘赌过程中产生大量的次级粒子。本工作利用单因素扰动系统的概念将扰动因素进行了规范化，在每个单因素扰动系统内，以粒子扩展属性的方式实现了扰动权乘子的更新。这一方法解决了权重技巧下输运与计数过程的微扰描述问题，并且利用次级粒子的存库与获取顺序避免了敏感性统计中对粒子历史进行频繁的回溯，从而降低了算法的运行时间开销。参考文献：1 KITADA T，Y A MA NE A,T A K ED A T.Improvement of the correlated sampling methodin Monte Carlo perturbation theoryC/Proc.Int.Conf.on the Physics and Reactors(PHY-SOR96).Japan:s.n.J,1996.2RIEF H.Generalized Monte Carlo perturbationalgorithms forcorrelated sampling and a second-原子能科学技术第58 卷order taylor series approachJ.Annals of Nu-clear Energy,1984,11(9):455-476.3WAGNER J C,PEPLOW D E,EVANS T M.Automated variance reduction applied to nuclearwell-logging problemsJ.Nuclear Technology,2009,168(3):799-809.1.0+JMCT直接模拟结果一JMCT关联抽样结果.关联抽样统计标准差+0.20.0184WAGNER J C,PEPLOW D E,MOSHER S W.0.8FW-CADIS method for global and regional vari-ance reduction of Monte Carlo radiation transport0.6calculationsJ.Nuclear Science and Engineer-0.4ing,2014,176(1):37-57.5MCKINNEY G W,IVERSON J L.Verificationof the Monte Carlo differential operator techniquefor MCNP,LA13098R.USA:Los AlamosNational Laboratory,1996.6PEREL R L,WAGSCHAL J J,YEIVIN Y.Monte Carlo calculation of point-detector sensi-tivities to material parametersLJJ.Nuclear Sci-ence and Engineering,1996,124(1):197-209.7HENDRICKS J S,PRAEL R E.Monte Carlonext-event estimates from thermal collisionsJ.Nuclear Science and Engineering，1991,10 9:150-157.8DENG Li,LI Gang,ZHANG Baoyin,et al.Ahigh fidelity general purpose 3-D Monte Carloparticle transport program JMCT3.oJI.Nucle-ar Science and Techniques,2022.doi:10.1007/$41365-022-01092-0.9WHALEN DJ,CARDON D A,UHLEJ L,etal.MCNP:Neutron benchmark problemsJ.Nuclear Physics&.Radiation Physics,1991.doi:10.2172/10103487.1o WONG C,ANDERSON J D,BROWN P,et al.Livermore pulsed sphere program:Programsummary through July 1971,UCRL-51144R.S.1.:s.n.,1992.11 CAREW J F,HU K,ARONSON A,et al.PWR and BWR pressure vessel fluence calculationbenchmark problems and solutions,NUREG/CR-6115,BNL-NUREG-52395R.S.1.:s.n.,2001.12 YANG Chao,CHENG Tangpei,DENG Li,etal.Development of 3-D parallel first-collisionsource method for discrete ordinate code JSNT-SJJ.Annals of Nuclear Energy，2 0 2 0,135:106942.