电子输运过程的确定论计算方法研究

李晓英,李云召,邵睿智,徐 宁

(西安交通大学 能源与动力工程学院,陕西 西安 710049)

电子输运过程在电子元件对辐射环境的响应、医疗辐射治疗、食品杀菌等领域应用广泛。目前,功能强大的蒙特卡罗程序,如MCNP[1]、Geant4[2]等,已经能准确模拟许多场景中的电子输运过程。但是,由于电子输运过程中的短平均自由程和高度前向峰等特点,蒙特卡罗方法计算效率比较低、计算成本非常大。相比之下,在解决分布式计算、深穿透过程模拟和低概率事件概率估计上,确定论计算方法的优势更加明显[3]。

已有研究表明,电子输运的确定论方法可有效代替蒙特卡罗方法[4]。由桑迪亚国家实验室Lorence等[5]开发的截面生成程序CEPXS可生成用于电子-光子耦合输运计算所需的多群勒让德截面,基于这些截面,可直接利用一维确定论粒子输运计算程序ONEDANT[6]进行电子-光子的耦合输运计算。CEPXS/ONELD[7]程序包已达到了很高的复杂精度,并已广泛用于一维几何的问题。研究人员在CEPXS的基础上后续开发了截面生成程序CEPXS-BFP[8],适用于使用离散纵标法的确定论电子输运程序,用于求解Boltzmann-Fokker-Planck(BFP)[9-10]形式的电子输运方程。确定论栅格计算程序DRAGON-5.0.8也增加了电子输运计算功能,基于CEPXS-BFP生成的FMAC-M格式电子库,利用高阶菱形差分法求解一维、二维[11]、三维[12]笛卡尔几何下的电子输运问题,Hébert等还在NJOY12[13]中开发了与DRAGON-5.0.8中电子输运计算相匹配的ELECTR[14-15]模块。

考虑到电子输运过程的蒙特卡罗模拟方法计算效率低,确定论方法的研发尚未成熟,本文在已知电子发射源分布的条件下,基于Boltzmann方程,设计一套针对EPICS2017中的电子截面数据(EEDL)[16]的多群电子截面库处理流程,以NECP-Atlas[17]程序的多群光子截面数据库处理模块为平台基础,对电子反应多群总截面和散射矩阵进行加工,直接利用压水堆堆芯分析软件NECP-Bamboo的栅格计算程序Bamboo-Lattice[18]中基于特征线方法的中子输运计算核心,构建确定论电子输运计算程序。

1 理论模型

1.1 电子输运方程

稳态线性Boltzmann输运方程可以描述中子、光子、电子等粒子在介质中的输运过程。由于电子输运过程和中子输运过程在输运方程形式上的一致性,Bamboo-Lattice中的中子输运计算核心可以在不改变计算模型的基础上直接用于电子输运过程的计算,只是需要为其提供多群电子截面数据作为输入数据。对于不同类型的粒子,Boltzmann输运方程碰撞项的具体形式有所不同,它敏感地依赖于粒子的类型以及粒子与介质原子相互作用的微观机制与特点。针对电子输运过程,散射产生项主要考虑电离、激发、轫致辐射、弹性散射等的影响,相应的多群形式Boltzmann方程如下:

(1)

(2)

1.2 电子截面数据库

电子-光子相互作用截面(EPICS)是ENDF/B[19]系统的一部分,用于补充ENDF/B的中子数据,开展工程应用中的耦合计算。EEDL是其中以ENDF-6格式存储的电子截面数据,具体内容列于表1,可以作为计算所需的核数据。表1中MF为数据类型标号,MT为反应道类型标号。

表1 EEDL中的数据列表

1.3 多群电子截面数据处理

NECP-Atlas是由西安交通大学核工程计算物理实验室(NECP)自主研发的多功能核数据处理程序,可将评价核数据库转化为多种形式的核数据库,为中子-光子耦合输运、燃耗计算、辐照损伤计算、屏蔽计算、释热计算以及不确定度量化等分析提供高精度的基础数据[20]。本文基于NECP-Atlas程序,处理EEDL[21],生成电子多群总截面和散射矩阵。

多群电子反应总截面的表达式为:

(3)

式中:σ501,g为MT501的第g群微观截面;σ501(E)为MT501中能量E下的微观点截面;φ0(E)为0阶典型电子权重谱。

电子的弹性散射不向外辐射能量,只改变其运动方向,可由EEDL中的525反应道获得馈送函数如下:

f525,l(E→g′)=

(4)

式中:f525,l(E→g′)为大角度弹性散射的馈送函数;Pn(μLAB)为勒让德多项式;f525(E,μLAB)为电子与原子发生大角度弹性散射后的实验室坐标系下的角度分布,在26/525中给出;μLAB为实验室坐标系下的散射余弦。

激发和轫致辐射的处理类似,假定自由电子与介质原子发生激发或轫致辐射后仍沿着原来的运动方向前进,则馈送函数形式如下:

(5)

式中:fx(E→E′)为x反应道中入射能量为E的自由电子与介质原子反应后能量变为E′的概率分布;ΔE为反应后的电子能损,其分布在MF26中给出;反应道x的取值为527或528。

入射能量为E的快电子与原子核外电子发生非弹性碰撞,使核外电子摆脱束缚能Ei成为自由电子的过程被称为电离。假定自由电子和反冲电子都沿入射电子原来的运动方向继续前进,且入射电子的能量不会转移给剩余原子。一般认为电离后能量较大的电子((E-Ei)/2

(6)

式中:i为电离壳层的索引;Er为反冲电子的能量;Eg′与Eg′+1为g′能群的上下限;fi(E→Er)为电子与介质原子的第i层亚壳层发生反应后的反冲电子能量分布,在26/534～26/572中依次给出。

由馈送函数fx,l(E→g′),可进一步求得多群电子散射矩阵:

(7)

式中:σx,l,g→g′为x反应道第g群到第g′群的l阶多群转移矩阵;φl(E)为入射能量为E的l阶注量率;σx(E)为x反应道下入射能量为E时的微观截面。

2 数值结果

本文通过对单核素均匀问题、化合物均匀问题、化合物非均匀问题等算例的计算,使用蒙特卡罗程序对该确定论程序进行验证,代表性算例问题的几何均如图1所示。所采用的源均匀分布在半径0.235 6 cm的圆域内,能量分布在0.001～10 MeV,共分为60个能群,能群的数据结构如图2所示,在低能区域时分布较多,在高能区域时较少。本文的电子源项是通过蒙特卡罗程序计算AP1000堆芯中D类组件[22]内SPND发射体处产生的电子源项。在确定论程序Bamboo-Lattice计算过程中,特征线条数设置为50,收敛准则设置为10-4。

图1 代表性算例问题的几何

图2 源项能群的数据结构

2.1 单核素均匀问题

所有单核素均匀问题中,核素的核数密度均为4.293 54×10-2barn-1·cm-1,只统计圆域内的通量分布,计算的归一化通量分布曲线及其与蒙特卡罗程序的对比如图3所示。计算各类单核素均匀问题时,Bamboo-Lattice程序和蒙特卡罗程序所需时间列于表2。由图3和表2结果可知,H、B、O、Al、S、Cl的Bamboo-Lattice结果与蒙特卡罗程序结果符合得较好,但Bamboo-Lattice的计算耗时明显比蒙特卡罗程序计算耗时少,计算效率可提高83～571倍。

a——H;b——B;c——O;d——Al;e——S;f——Cl

表2 单核素均匀问题计算时长

2.2 化合物均匀问题

化合物均匀问题计算了Al2O3与含硼水两类,其中Al2O3中O、Al元素的核数密度分别为3.396 898×10-2和2.264 599×10-2barn-1·cm-1,含硼水中H、O、B 3种元素的核数密度分别为4.624 176×10-2、2.315 550×10-2、2.307 683×10-5barn-1·cm-1。该类问题只统计了圆域内的通量分布,归一化通量分布如图4所示。计算化合物均匀问题的时长比较列于表3。可以看出,化合物均匀问题中Al2O3和含硼水的Bamboo-Lattice结果与蒙特卡罗程序结果符合得较好,但计算效率提高了82～109倍。

a——Al2O3;b——含硼水

表3 化合物均匀问题计算时长

2.3 化合物非均匀问题

化合物非均匀问题为1个两层结构栅元问题,内层为Al2O3,外层为含硼水,材料的元素构成列于表4。该问题的归一化通量谱对比如图5所示,计算时长列于表5。可以看出,仅由较小核素构成的Al2O3层和含硼水层的Bamboo-Lattice结果与蒙特卡罗程序结果的趋势相一致,且在计算时长上效率提高约169倍。

a——Al2O3层;b——含硼水层

表4 非均匀问题构成元素及密度

表5 化合物非均匀问题计算时长

3 结论与展望

本文基于核数据处理程序NECP-Atlas将EEDL中连续能量的电子截面数据加工为多群形式,利用栅格计算程序Bamboo-Lattice中的中子输运计算核心进行电子输运计算。对单核素均匀问题、化合物均匀问题和化合物非均匀问题进行计算,并与蒙特卡罗程序计算结果进行对比验证。结果表明,该确定论程序可正确模拟电子输运过程,且计算效率比蒙特卡罗程序快80倍以上。后续研究将利用该确定论电子输运计算软件计算实际压水堆中自给能中子探测器(SPND)的响应电流,用以替换现有的蒙特卡罗程序[23]。