-
加速器驱动次临界系统(ADS)是一种由加速器、散裂靶、次临界堆三个模块组成的装置,通过质子束轰击重金属散裂靶产生的外源中子驱动次临界反应堆运行。ADS专门用于嬗变乏燃料中的长寿命次锕核素和裂变产物,增殖核燃料和产能,提高核燃料的利用率和降低乏燃料对环境的影响[1]。
有效增殖因数
$ k_{\mathrm{eff}} $ 是关系到反应堆安全运行的重要监测参数。虽然ADS具有次临界下的固有安全性,但仍需确保意外反应性引入下不会诱发超临界事故。次临界系统反应性的控制和调节必须确保反应堆始终处于一个相对稳定的次临界度。随着ADS系统的运行,核燃料的燃耗会使堆芯次临界度出现涨落,系统需要通过改变加速器束流流强(即外中子源强度)来调节核燃料的裂变率以补偿反应性变化带来的功率漂移[2]。因此,反应性的精确测量对于ADS安全运行非常重要。以ADS在次临界状态
$ k_{\mathrm{eff}} $ =0.98情况下为例,为维持运行功率保持不变,当$ k_{\mathrm{eff}} $ 降低1%,加速器束流强度就要增加50%[3]。核反应堆实验领域存在多种不同的测量方法探索和使用。针对次临界堆的反应性测量技术,目前普遍使用的是脉冲中子源法(PNS)[4],即向反应堆连续注入一定频率脉冲中子,测量堆内中子衰减规律从而计算出次临界度。其中,基于脉冲中子源的肖斯特兰德积分面积比方法可以给出比较精确的次临界度测量值。该方法已广泛用于KUCA[5],VENUS-F[6]和YALINA[7]等零功率次临界反应堆的在线反应性测量实验,测量值和蒙卡程序[8]的计算值显示出了较好的一致性。先前针对PNS方法的实验测量和模拟计算大多是建立在热中子反应堆的基础上,使用D-T、D-D反应的中子管来模拟外中子源。为了研究较为真实ADS装置下PNS方法的适用性,本文参考CiADS拟定建设方案,使用蒙卡程序建立次临界铅铋快堆模型。利用蒙卡程序模拟500 MeV了脉冲质子束轰击铅铋靶与次临界堆耦合下的粒子输运过程并得到了中子信息。随后,我们利用PNS面积比方法实验原理编写了处理程序,展开了次临界堆的反应性测量的模拟研究。
蒙卡程序无法直接模拟脉冲质子重复注入堆内达到准稳态平衡时中子通量密度的时间响应曲线,以及脉冲中子法所需要的稳定缓发中子本底。但可以给出单个脉冲质子束轰击堆芯后瞬发及缓发中子时间响应数据,使用程序对数据在时间上进行累计叠加,就能得到重复脉冲束流注入下堆芯内中子通量随时间的实际变化以及缓发中子的本底。本文使用 Python[9]编程语言编写程序处理脉冲数据积分与实验测量过程,实现了脉冲中子源法对 CiADS 次临界度的模拟测量。本文的研究目标:第一,给出CiADS不同次临界度下脉冲中子源面积比法的测量精度和实际适用性;第二,由于次临界堆内空间效应的影响,比较堆内探测器位置摆放选择不同而导致的测量差异,寻找测量结果较为准确的位置用于摆放探测器。
-
通过研究次临界堆芯在脉冲中子作用下的响应用来分析反应堆动态变化规律,以此来计算出反应性。向堆内注入一束脉冲中子,将引起堆芯燃料中的易裂变物质发生链式裂变反应,产生瞬发中子和缓发中子,在次临界中,裂变链式反应无法稳定维持,中子会在堆内持续衰减,其规律与反应堆次临界度有关。在周期性的脉冲中子作用下,当重复频率
$ R $ 远远大于缓发中子的衰减常数$ \lambda_{i} $ ,远远小于瞬发中子的衰减常数$ \alpha_{\rm{np}} $ ,即$ \alpha_{\rm{np}} \gg R $ $ \gg \lambda_{i} $ ,在准稳态平衡条件下,将反应堆中子分为瞬发中子和缓发中子两部分,次临界中缓发中子几乎不随时间改变形成一个常数本底,分别求出瞬发中子数对时间的面积积分和缓发中子数对时间的面积积分,并由它们的比值求出反应性$ \rho $ ,如图1所示。PNS面积比方法基于点堆动力学的模型推导,将中子通量密度
$ N $ 分为瞬发中子通量和缓发中子通量两部分[10-11]:式中:
$ N_{\rm{p}} $ 为瞬发中子通量密度;$ N_{\rm{d}} $ 为缓发中子通量密度;$ C_{i}(t) $ 为第$ i $ 组缓发中子先驱核浓度;$ \lambda_{i} $ 为第$ i $ 组缓发中子先驱核衰减常数;$ \beta_{{i}} $ 为第$ i $ 组缓发中子总有效份额;$ \Lambda $ 为中子代时间;$ \rho $ 为反应性;$ S(t) $ 为脉冲中子源。脉冲中子源
$ S(t) $ 的时空分布函数公式如下:其中
$ S_{0} $ 为总的脉冲中子数。当时间
$ t $ 处于两个脉冲之间时,$ S(t) = 0 $ ,将式(1)简化并从$ t = 0^{-} $ 到$ T $ 进行积分,$ T $ 为脉冲周期,并假设变量的整个周期的积分为零,即并令:
$ A_{\rm{p}} $ 和$ A_{\rm{d}} $ 分别为瞬发中子和缓发中子在$ T $ 时间内的积分面积,最终简化导出瞬发中子通量密度对时间的积分面积$ A_{\rm{p}} $ 以及缓发中子通量密度对时间的积分面积$ A_{\rm{p}} $ ,得到反应性$ \rho $ 的表达式: -
在蒙卡程序上构建的CiADS反应堆堆芯布局如图2所示,堆芯中心圆柱结构为散裂靶区,此处加速器质子束流通过靶窗处轰击铅铋散裂靶产生外源中子。堆芯组件类型为六角形,燃料组件包含61根燃料棒,每根燃料棒由95根长10 mm、半径5.4 mm的圆柱体二氧化铀燃料芯块组成,235U的富集度为19.75%。反射组件包含7根不锈钢棒,半径17.5 mm,屏蔽组件包含7根碳化硼棒,半径17.5 mm。冷却剂和散裂靶组件材料为铅铋合金(Pb-46%, Bi-54%)。燃料组件主要设计参数如表1所列。通过调整堆芯中燃料组件的数目配置了8种不同次临界深度的堆芯布局,如图3所示。
参数名称 数值 参数名称 数值 组件类型 有盒六角形 包壳内径 10.8 mm 组件长度 2 700 mm 芯块外径 10.4 mm 燃料棒数量 61根 活性区燃料 UO2芯块 六角管材料 高硅铁马钢 富集度 19.75% 包壳外径 12.8 mm 活性区高度 950 mm 燃料棒中心距 13.4 mm 气腔长度 600 mm 蒙卡程序中使用F4栅元计数卡记录中子通量,设定燃料棒中心处芯块和反射组件中心出处芯块的圆柱体栅元,记录该栅元处的中子通量随时间的响应。模拟设定单次
$ 10^6 $ 个500 MeV的质子注入到散裂靶中,记录各探测处体栅元从0到500 s内对数插值200道的中子通量计数。为了计算不同位置探测器摆放对计算结果的影响及差异,选择10个组件位置放置探测器,标号1到10,分为两组,标号1–5为第一组,标号6–10为第二组。根据离堆中心位置, 由内圈到外圈到依次排布5个探测器,1号和6号为第一圈组件的两个探测器,2号和7号为第二圈组件的两个探测器,同理依次递推,位置标号见图2。 -
假设堆芯处于稳态,材料温度和冷却剂流量分布不变,中子核截面不变,因此, 反应性不随着脉冲周期变化。MCNPX模拟粒子输运过程中,每个质子注入的反应过程中是相互独立的,在时间上可以线性叠加;重复注入质子束流引发的中子通量衰变曲线可以用单质子注入的中子通量变化曲线通过叠加来得到[12]。
在MCNPX的模拟过程中,截面、材料密度和几何模型没有变化。
$ N $ 个脉冲的中子通量与时间$ \varPhi(t) $ 的函数可以表示为$ N $ 个第一个脉冲后的中子通量衰减曲线函数$ \varPhi(t) $ (随时间移位$ T $ ,即脉冲周期)之和,如下所示:MCNPX给出了从0到500 s,按照对数插值取200道的时间谱数据,编写Python程序时,将数据通过样条插值拟合为函数,获得时间间隔1 μs,长度
$ 5 \times 10^{8} $ 个数据点线性列表$ X $ 。设定脉冲周期$ T = n $ ms,中子通量衰减曲线为500 s,单位时间间隔1 μs,达到准稳态平衡需叠加$ 5 \times 10^{5}/n $ 次,核心代码过程如下。#所需的循环叠加次数5E5/n for i in range{[int(5E5)/n]}: #n*1000为脉冲周期长度 j = (i+1)* n* 1000 for k in range[j, int(5E8)]: #每单位脉冲周期叠加 X_pns[k] = X_pns[k]+X[k-j] 以
$ n = 5 $ 为例,前18次脉冲叠加过程如图4所示。当完整叠加中子通量衰减曲线后,根据图中得到准稳态平衡曲线和面积比方法计算公式,分别计算瞬发和缓发区域的面积,得到反应性$ \rho $ 。 -
蒙卡程序计算给出了8种不同次临界深度下单脉冲质子束注入下的中子通量随时间的衰减曲线,对数据进行归一化处理,选取了次临界深度差别最大的两组数据比较,如图5所示。
由图5可知,在不同次临界度下,瞬发中子和缓发中子通量变化呈现两个峰,二者变化趋势基本一致。次临界度较深时[图5(a)],探测器位置分别由近到远排布,由于中子的传播速度有限,所以在不同位置到达峰值的时刻存在细微差别;在次临界度较浅时[图5(b)],未表现出峰值的时间随位置的变化,且相对与次临界度较深时峰值时间延后。
另外,两种反应性下中子通量曲线在
$ 10^{-9} $ 到$ 10^{-6} $ s内的瞬发中子通量前端变化有明显的凸起,越靠近散裂中心区域的位置中子通量凸起越高。考虑到散裂外源中子通量达到峰值的时间在$ 10^{-7} $ s左右,并且由于探测器无法区分外源中子和裂变中子,我们认为该部分的凸起是受到外源中子的影响所造成的。外源中子首先穿过内圈组件,所以内圈相对外圈的中子计数要高,由于外源中子的穿透距离有限,随着距离的增加, 外源中子数量就越来越少。因此表现出曲线在$ 10^{-9} $ 到$ 10^{-6} $ s内不同探测位置处的通量变化的分散。由于散裂反应产生的外源中子通量时间谱相对中子管产生的中子通量时间谱更加展宽,外源中子的计数也更多包含在探测器计数中,会对模拟结果产生一定的影响。图5(c)中,a1为
$ k_{\mathrm{eff}} $ =0.846 74时1号位置处探测器计数,b1为$ k_{\mathrm{eff}} $ =0.986 10时1号位置处探测器计数,由此可知,在有效增殖因数$ k_{\mathrm{eff}} $ 较低时,散裂中子所占的瞬发时刻总中子的份额较高,而随着有效增殖因数$ k_{\mathrm{eff}} $ 升高,由裂变产生的中子增多,散裂中子的所占中子份额有所下降。 -
将蒙卡程序给出的数据通过Python叠加计算模拟了8种不同深度次临界下周期脉冲质子束(
$ T=5 $ ms)重复注入下的中子通量随时间的衰减曲线,选取了次临界深度差别最大的两组数据比较,如图6所示。由图6(a)和6(b)可知,由于探测器所处的位置不同,叠加后形成的缓发中子常数本底在燃料组件中比反射组件中要高。由于反射区位置不发生裂变反应,只能探测到活性区向外泄露的中子。另外,在外源强度不变的情况下,随着增殖因数的提高[图6(b)],堆内裂变产生的中子增多,中子通量更高,缓发中子本底也相应提高。
分别从数据中提取出缓发区域中子面积
$ A_{\mathrm{d}} $ 和瞬发区域中子面积$ A_{\mathrm{p}} $ ,如图1所示,根据面积比式(5)和反应性公式:计算出反应性
$ \rho $ 和有效增殖因数$ k_{\mathrm{eff}} $ 。缓发中子份额$ \beta_{\mathrm{eff}} $ 是通过反应堆动态参数的蒙特卡罗计算研究[13]方法给出的计算结果,不同反应性下$ \beta_{\mathrm{eff}} $ 的计算结果如表2所列。$ k_{\mathrm{eff}}$ $ \beta_{\mathrm{eff}} $ $ k_{\mathrm{eff}}$ $ \beta_{\mathrm{eff}} $ 0.846 74 0.007 421 0.959 58 0.007 429 0.862 58 0.007 428 0.974 50 0.007 432 0.900 13 0.007 416 0.981 93 0.007 435 0.939 38 0.007 421 0.986 10 0.007 432 表3中给出不同次临界深度下,PNS面积比法的计算结果。堆内各标号位置处由面积比方法计算出的增殖因数
$ k_{\mathrm{eff}} $ 与基准值$ k_{\mathrm{eff}} $ 的偏差如表4所示。$ k_{\mathrm{eff}} $ 基准值是在蒙卡程序中kcode临界源卡给出的计算结果,计算设定100 000个源中子,在两处位置,跳过前50次迭代,累计迭代250次。$ k_{\mathrm{eff}} $基准值 第一组探测器$ k_{\mathrm{eff}} $结果 第二组探测器$k_{\mathrm{eff}} $结果 1 2 3 4 5 6 7 8 9 10 0.846 74 0.846 52 0.864 32 0.864 08 0.863 32 0.861 24 0.837 01 0.859 12 0.860 27 0.861 44 0.859 32 0.862 58 0.862 26 0.877 29 0.877 42 0.875 81 0.875 26 0.854 55 0.872 54 0.875 42 0.875 88 0.875 70 0.900 13 0.897 72 0.907 41 0.908 76 0.908 17 0.907 03 0.893 07 0.905 50 0.909 60 0.908 86 0.908 34 0.939 38 0.936 85 0.941 83 0.942 80 0.942 69 0.942 29 0.935 00 0.940 51 0.942 02 0.942 47 0.942 29 0.959 58 0.958 39 0.960 97 0.961 58 0.961 42 0.961 53 0.957 27 0.960 36 0.961 62 0.961 62 0.961 54 0.974 50 0.973 28 0.974 59 0.975 08 0.975 03 0.974 90 0.972 67 0.974 21 0.974 77 0.974 87 0.974 68 0.981 93 0.980 41 0.981 19 0.981 53 0.981 46 0.981 44 0.980 05 0.980 98 0.981 37 0.981 41 0.981 37 0.986 10 0.985 26 0.985 78 0.986 00 0.985 99 0.985 94 0.985 03 0.985 65 0.985 89 0.985 94 0.985 93 $ k_{\mathrm{eff}} $基准值 第一组探测器反应性偏差 第二组探测器反应性偏差 1 2 3 4 5 6 7 8 9 10 0.846 74 $-30$ 2 402 2 369 2 268 1 988 −1 372 1 701 1 857 2 015 1 728 0.862 58 $-43$ 1 943 1 960 1 751 1 679 −1 089 1 323 1 700 1 760 1 736 0.900 13 $-298$ 891 1 055 983 845 $-878$ 658 1 156 1 067 1 004 0.939 38 $-287$ 276 386 373 328 $-498$ 127 298 349 328 0.959 58 $-129$ 150 216 199 211 $-251$ 84 221 221 212 0.974 50 $ -128$ 9 61 55 42 $-193$ −30 28 38 18 0.981 93 $-157$ $ -76 $ $-41$ $-48 $ $-50$ $-195$ −98 $-58 $ $-53$ $-58$ 0.986 10 $-86$ $-32$ $-10$ $-11$ $-16$ $-110$ $-46$ $-21$ $-16$ $-17$ 表4中反应性偏差的计算方式为
单位pcm。
从表4中可以看出,随着有效增殖因数
$ k_{\mathrm{eff}} $ 提高,考虑两组探测器位置之间的差别时,第一组探测器与第二组探测器给出的有效增殖因数,与基准值的偏差越来越小。在次临界度较深下,由于空间效应的和点堆模型的不适用,两组探测器之间的偏差明显,随着增殖因数提高,随位置差别的影响也越来越小。不同次临界度下,PNS面积比法计算结果与基准值的偏差的大小决定了该方法是否可用,以及适用的测量增殖因数$ k_{\mathrm{eff}} $ 范围,根据EFIT反应堆设计给定的反应性测量偏差[14]范围约为600 pcm,以此确定PNS面积比法可达到用于测量次临界深度的范围。图7和图8给出了第一组和第二组探测器位置的
$ k_{\mathrm{eff}} $ 计算结果与基准值的对比,有效增殖因数$ k_{\mathrm{eff}} $ 的基准值为0.846 74、0.900 13、0.939 38、0.986 10,给出基于基准值的600 pcm上下限。图7和图8结果显示,最靠近堆中心的第一圈组件位置探测器(标号1)和(标号6)获取数据的计算结果相对外圈组件的$ k_{\mathrm{eff}} $ 结果有明显的偏低。由于探测器不能区分外源中子和裂变中子,次临界度越深,第一圈组件受到外源中子的影响越显著,瞬发中子通量计数中外源中子的贡献越高。这实际上变相提高了瞬发中子面积$ A_{\rm{p}} $ , 使得模拟的增殖因数$ k_{\mathrm{eff}} $ 偏低。从图2看出,探测器6位置比探测器1距离散裂靶更近,外源中子份额贡献更高,因此探测器6的模拟数值比探测器1更加偏低。在低次临界度下($ k_{\mathrm{eff}} <0.9 $ )距离散裂靶较远的探测器位置PNS法的模拟的数值都偏高;靠近散裂靶的探测器由于受外源中子影响,PNS模拟结果会偏低,两者趋势相互抵消使得实际测量反应性偏低。通过分析表4中的反应性偏差,结合图2中探测位置分布可知,1号和6号探测器都处于中心第一圈,当
$ k_{\mathrm{eff}} $ =0.946 74时,1号测量结果偏差(−30 pcm)相较于2、3、4号的测量偏差(2 042、2 369、2 268 pcm)较低,但6号的测量偏差较大(−1 372 pcm)。6号偏差较大是因为它比1号的位置更靠近散裂靶中心。图9给出相同有效增殖因数下1号和6号探测位置的中子通量随时间演化,可以看出,6号位置处外源中子的占比更大。我们认为,1号位置的模拟结果相比与其它位置更接近基准值的原因,是受到外源中子的影响和PNS系统方法误差共同导致。随着次临界堆反应性升高,外源中子占比逐渐降低,探测位置1、6的偏低程度也在变小。另外,当$ k_{\mathrm{eff}} $ 升到0.97时,从表4中可知,临近散裂靶内圈组件放置探测器所模拟的PNS偏差反而比外圈要高。随着次临界堆在有效增殖因数
$ k_{\mathrm{eff}} $ 的提高,面积比法测量结果越接近基准值,结合表4可以看出,当$ k_{\mathrm{eff}} $ = 0.939 38时,偏差小于600 pcm。据此认为在有效增殖因数$ k_{\mathrm{eff}} $ >0.94时,PNS面积比方法的模拟计算结果可以满足精度要求,可以作为次临界堆有效的反应性测量手段。 -
本文利用蒙卡软件模拟脉冲单质子脉冲注入次临界反应堆的粒子输运过程,采用 Python 编写程序处理中子数据实现脉冲束流下中子的叠加过程,实现了重复脉冲质子束周期注入下的次临界堆中子学模拟。最终实现了在CiADS次临界反应堆上使用PNS面积比方法进行反应性测量的模拟,并与
$ k_{\mathrm{eff}} $ 基准值进行对比,基于8种的不同次临界深度与测量位置处的PNS面积比法模拟,结果显示:(1) 随着堆芯有效增殖因数
$ k_{\mathrm{eff}} $ 增大,PNS方法的计算结果越相对准确,有效增殖因数的偏差越小,理论上最低可达到10 pcm左右。证明PNS方法在高有效增殖因数下($ k_{\mathrm{eff}} >0.94$ ),可以较为准确地测量反应性。(2) 由于受到外源中子的影响,探测位置对反应性测量结果有较大影响。紧邻散裂靶的燃料组件处的模拟结果显示出较强的空间效应,外层燃料组件和反射组件中的模拟值基本一致。考虑燃料组件内的中子通量数量比反射组件高出几个数量级,对探测器寿命和设计的影响较大。探测器放置在燃料外区附近的反射组件位置处较为合适,后续可以考虑在CiADS装置开展实验测量进行验证。
Monte Carlo Simulation Analysis of CiADS Reactivity Measurement Based on Pulsed Neutron Source Method
doi: 10.11804/NuclPhysRev.39.2021077
- Received Date: 2021-10-15
- Rev Recd Date: 2021-12-23
- Publish Date: 2022-03-01
-
Key words:
- neutron /
- reactivity /
- area method /
- CiADS /
- pulsed neutron source
Abstract: Pulsed neutron source method(PNS) is an important technique for the reactivity measurement of accelerator driven sub-critical systems. The CiADS sub-critical reactor model is established using Monte Carlo code to simulate the neutron transport process bombarded by pulsed proton beam, and the time spectrum of neutron flux is obtained. A code is developed by Python programming language to finish the pulsed neutron superposition process. It can obtain a stable delayed neutron background and the variation of neutron flux at different positions of the core under the continuously pulsed period proton injection. The reactivity measurement simulation of the sub-critical reactor by PNS method come true. Subsequently, the PNS method is used to simulate the reactivity of the CiADS core under different sub-criticalities and detection positions, and the results is compared with the reference values. It shows that this method can accurately predict the reactivity of sub-critical reactor at high effective multiplication factor (
Citation: | Xiang LI, Xunchao ZHANG, Huan JIA, Yuan HE. Monte Carlo Simulation Analysis of CiADS Reactivity Measurement Based on Pulsed Neutron Source Method[J]. Nuclear Physics Review, 2022, 39(1): 127-133. doi: 10.11804/NuclPhysRev.39.2021077 |