厚靶轫致辐射在粒子输运模拟中的实现与验证

李华锋，盛 伟，韩 斌，卓 俊，陶应龙

（西北核技术研究所，西安710024）

厚靶轫致辐射在粒子输运模拟中的实现与验证

李华锋，盛 伟，韩 斌，卓 俊，陶应龙

（西北核技术研究所，西安710024）

为提升自行开发的粒子输运蒙特卡罗模拟软件PHEN对电子输运模拟的处理能力，根据厚靶轫致辐射（thick－target bremsstrahlung，TTB）模型，在PHEN中加入了自编的TTB模块，分析了加入TTB模块后PHEN的模拟计算速度、TTB模块生成的轫致辐射光子能量及TTB模块模拟电子输运的适用范围。通过与MCNP中轫致辐射截面及电子射程的对比，验证了TTB模块模拟电子输运的正确性。计算了加入TTB模块前后正电子湮没峰计数的变化，用于判断使用TTB模块模拟的合理性。计算结果表明：当材料尺寸小于电子在材料中的射程时，TTB模型不考虑电子输运的计算结果是不合理的，应采用压缩历史算法处理电子输运；当材料尺寸大于电子在材料中的射程时，TTB模型计算速度快，计算结果与压缩历史算法的计算结果吻合较好，可用厚靶轫致辐射处理电子输运。

厚靶轫致辐射；蒙特卡罗；粒子输运；正电子湮没峰

随着现代计算机技术的发展，蒙特卡罗方法已发展成为粒子辐射输运研究领域中的一种重要方法。在粒子输运的蒙特卡罗模拟中，光子只需经历很少的碰撞就会被介质完全吸收。但对电子而言，由于其在每次碰撞过程中一般只损失很少的能量（约几十电子伏），电子在被介质吸收前要经历多次碰撞。与光子模拟相比，电子模拟的时间相当长，受计算机计算能力的限制，往往无法实现对电子输运过程的直接蒙特卡罗模拟。为了克服对电子直接模拟的困难，Berger提出了电子的压缩历史模拟算法，MCNP（Monte Carlo N－particle transport code）基于该算法实现了对电子的输运模拟。

在不考虑电子／正电子输运的光子模拟中，MCNP不采用压缩历史算法模拟电子输运，而是利用厚靶轫致辐射（thick－target bremsstrahlung，TTB）模型对次级电子近似模拟，计算速度得到了明显提升。西北核技术研究所在射线与物质相互作用数值模拟仿真平台研制中，开发了粒子输运蒙特卡罗模拟软件PHEN（the transportation program of photon，hadron，electron，and neutron）。该软件能够模拟光子、质子、电子、中子在常见介质中的输运过程。现有的PHEN在模拟电子时只可采用压缩历史算法，计算速度较慢。为了进一步完善电子输运模拟算法，提升PHEN对电子的处理能力，本文基于简化TTB模型，在PHEN中加入了TTB模块，并进行了验证计算和相关讨论。

1 TTB模型及算法

在粒子输运蒙特卡罗模拟中，电子输运模拟计算缓慢，模拟时间相当长。在不考虑电子／正电子输运的光子模拟中，可认为电子／正电子在其生成位置被局部吸收，生成轫致辐射光子，这种处理模型即为TTB模型。

在TTB模型中，生成的轫致辐射光子数、轫致辐射光子的能量及角度和空间分布由预先计算好的截面确定。该过程仅需用光子输运方程描述为

其中，E＝Eω；v＝vω；E为光子能量；ω为光子运动方向矢量；v为光子速率；v为光子速度矢量；f（r，E，t）为光子的分布函数；σ（r，E）为总截面；K（r，E′→E）表示点r处，能量为E′、运动方向为ω′的光子碰撞后的能量、角度分布函数，其物理意义为光子在r处经碰撞由能量E′变为能量E的概率；g（r，E，t）为光子源的分布函数；Nbrem（E′）表示能量为E′的光子生成的轫致辐射光子数r，E）表示点r′处，能量为E′、运动方向为ω′的光子单次相互作用事件生成的轫致辐射光子的能量、角度和空间分布函数［1］。

MCNP在光子输运模拟mode p模式中，当phys：p卡中ides＝0时，使用TTB模型模拟次级电子／正电子［2］。MCNP通过子程序brem对厚靶轫致辐射截面初始化，通过子程序ttbr（d）完成轫致辐射光子的生成，其中，d为电子／正电子的径迹长度。ttbr（d）在MCNP中的调用情况如图1所示，其中，ttbr（0）为厚靶轫致辐射。从图1可以看出，在不考虑电子输运的模拟中，MCNP在所有次级电子生成后均通过TTB模型将次级电子转化为光子。

值得注意的是，MCNP采用的是简化TTB模型，仅考虑了厚靶轫致辐射光子的能量分布，而没有考虑厚靶轫致辐射光子的角度和空间分布。轫致辐射光子的位置和运动方向用电子的位置和运动方向近似模拟。

下面结合PHEN使用的轫致辐射截面数据结构，介绍简化TTB模型的计算方法。

轫致辐射生成光子的能量与初始电子的能量比为k，k＝W／E。能量Ee的电子被完全吸收，生成的轫致辐射光子数为

其中，Ea为材料的吸收能；μ（E，k）是能量为E的电子在单位路程上发生轫致辐射生成能量kE光子的概率；σ（E，k）是能量为E的电子通过轫致辐射生成能量kE光子的微观截面；n为材料的原子数密度；R（E）为电子的CSDA（continuous slowing down approximation）射程，单位为b－1。电子射程在电子截面初始化时已计算，可作为已知函数。

PHEN的电子截面数据库与PENELOPE的电子截面数据库相似，其存储的电子折合轫致辐射微分截面（“scaled”bremsstrahlung differential cross sections）为［3-4］

对于由多种元素组成的材料，其电子轫致辐射微分截面为

其中，Zeq为材料的等效原子序数；χeq（Zeq，E，k）为等效的电子折合轫致辐射微分截面。

令

则

轫致辐射光子电子能量比k分布的累积概率密度为

轫致辐射光子能量Ep分布的累积概率密度为

式（8）右端分子中第一项表示电子能量大于Ep时，轫致辐射生成的能量小于等于Ep光子的数目；第二项表示电子能量小于Ep时生成的光子数目；分母表示电子轫致辐射生成的光子总数。

简而言之，轫致辐射光子数通过函数Np（Ee）计算，轫致辐射光子能量Ep由累积概率密度FE（Ep；Ee）抽样得到，轫致辐射光子的位置和运动方向与电子的位置和方向相同。

2 TTB算法在PHEN中的实现

由式（2）至式（8），可实现在PHEN中采用TTB模型模拟次级电子／正电子。程序实现时需要注意的是，与电子射程有关的积分中，能量为降序。MCNP中可计算的电子能量为1keV～100MeV，电子能量网格为降序；PHEN中可计算的电子能量为50eV～1GeV，电子能量网格为升序。

通过对比MCNP和PHEN中采用TTB模型涉及的物理量，可分析采用TTB模型的可行性。在电子截面初始化的基础上，采用TTB模型使用的物理量为轫致辐射微分截面和电子射程。

2.1 轫致辐射微分截面

MCNP中涉及电子的数据库有el和el03。钨在el和el03中的电子数据表格分别为74000.01e和74000.03e。现有MCNP程序在计算FE（Ep；Ee）时使用的是ITS3.0数据格式表格74000.03e中的电子轫致辐射微分截面［5-6］。PHEN中电子折合轫致辐射微分截面存放在pdebrZZ.p06文件中［3］。截面数据为Bethe－Heitler公式的计算结果，该公式通过波恩近似推导得出［7］。MCNP的截面数据库中的电子能量网格和PHEN的截面数据库中的电子能量网格一致，均为57个节点。以钨为例，两个程序的数据库中，轫致辐射截面一致，如图2所示。

2.2 电子射程

PHEN中使用的电子射程单位为cm，MCNP中使用的电子射程单位为b－1。若材料的原子数密度为n，单位cm－3，则PHEN中使用的电子射程和MCNP中使用的CSDA射程间的关系式为

图3给出了电子在钨中的射程，包括PHEN、MCNP和经验公式［8］的计算结果。

从图3可以看出，PHEN和MCNP中的电子射程在100keV以上完全符合。在100keV以下，PHEN中的电子射程较MCNP中的电子射程偏小。100keV以下电子射程的差异对TTB计算结果的影响分析如下：

从式（6）和式（8）可知，Np（Ee）和FE（Ep；Ee）均与电子射程有关。程序输出的电子在钨中的Np（Ee）如图4所示。从图4可以看出，Np（Ee）在高能区完全一致，在低能区略有差异。不同能量的电子的轫致辐射截面有明显的差异，低能电子生成的光子数远少于高能电子生成的光子数。从图4（b）可以看出，PHEN和MCNP每模拟1 000个5keV的电子，轫致辐射生成的光子数分别是7和8，即两程序中电子射程在低能区的差异对轫致辐射光子数的影响可忽略不计。

程序输出的电子在钨中的轫致辐射光子能量分布的累积概率密度FE（Ep；Ee）如图5所示。从图5（a）可以看出，PHEN和MCNP在100MeV时的轫致辐射光子能量分布完全相同。从图5（b）可以看出，虽然5keV时的电子射程有细微差异，但其轫致辐射光子能量分布无明显差异。从图5（c）可以看出，两程序在1keV处的轫致辐射光子能量分布有明显差异，PHEN抽样得到的光子能量偏低。即便如此，这种差异也仅对能量为1keV以下的光子模拟有影响，在计算时可以忽略。

综上所述，低能电子射程的细微差异对TTB模型计算结果的影响可忽略不计，这也为在PHEN中实现TTB模型奠定了基础。

3 验证与讨论

3.1 模型验证

为了确保PHEN中加入的TTB模块的正确性，需进行验证计算。验证计算模型为半径5cm的钨球，密度19.3g·cm－3，光子源能量5MeV，位于钨球球心，记录穿出钨球球面的光子能谱。MCNP在mode p模式下计算，phys：p卡取默认值（ides＝0）。PHEN用加入TTB模块的程序计算。计算结果如图6所示。

图6中，从低能到高能的3个峰依次为钨的特征X射线峰（57.98，59.32，66.95，67.25keV 4个能峰重叠）、电子对湮没峰（511keV）和光子源峰（5 MeV）。两程序计算得到的光子能谱完全一致，验证了PHEN中加入TTB模块的正确性。

3.2 TTB对PHEN计算速度的影响

在PHEN中加入TTB后，原程序包中的压缩历史算法模拟电子输运可由厚靶轫致辐射模型代替，由于忽略了电子输运过程，可大幅提高模拟计算速度。以验证计算模型为例，使用压缩历史算法，1s内可模拟计算138个初级光子；使用TTB模型，1s内可模拟计算9 150个初级光子，速度提升了约65倍。对更高能量的光子源和更大尺寸的介质材料，使用TTB模型的计算速度可提升更多。

3.3 轫致辐射光子能量

能量Ee的电子在材料中生成的轫致辐射光子的总能量为

经计算得，50MeV的电子在钨中经厚靶轫致辐射产生的光子数为41.38，每个光子的平均能量为0.801MeV，轫致辐射光子总能量占入射电子能量的66.3%。3.125MeV的电子在钨中经厚靶轫致辐射产生的光子数为7.88，每个光子的平均能量为0.056MeV，厚靶轫致辐射光子总能量占入射电子能量的14.1%。

从上述计算可以看出，虽然电子通过TTB模型直接转化为光子，但并非将电子的全部能量转化成了轫致辐射光子能量。这是由于TTB模型并非只考虑了电子的轫致辐射过程，还考虑了电子与物质的弹性与非弹性碰撞。TTB模型中使用的电子射程为CSDA射程，计算公式为

其中，Sin（E′）为碰撞阻止本领；Sbr（E′）为辐射阻止本领［3］。可认为电子能量未全部转化为轫致辐射光子能量，是由于电子受弹性或非弹性碰撞后，能量在局部沉积造成的。若计算电子射程时不考虑碰撞阻止本领，则经TTB模型后电子的能量会全部转化为轫致辐射光子能量。

3.4 TTB的适用条件

以验证计算模型为例，分别计算了不同钨球半径时的光子能谱。图7给出了使用PHEN＋TTB和MCNP两种计算模式下得到的钨球表面的出射光子能谱。由图7可见，钨球半径为0.01cm时，PHEN＋TTB模式和MCNP模式给出的能谱仅在3MeV以上符合较好。钨球半径为0.2cm或0.6 cm时，两种模式的结果在1MeV以上符合较好。

从图6可知，电子经轫致辐射生成的光子中，80%的光子能量小于电子能量的1／10。由图8给出的光子平均自由程可知，光子在500keV以下主要发生光电效应，在1～4MeV主要发生康普顿散射，在4～5MeV主要发生康普顿散射和电子对效应。因此，图7中低能段主要反映了两种模式对电子轫致辐射的计算差异，高能段主要反映了两种模式对康普顿散射和电子对效应的计算差异。

除湮没外，正电子的输运过程与电子在材料中的输运过程相似，因此，可通过正电子的输运情况大致判断电子在材料中的输运情况。表1给出了不同钨球半径时的正电子湮没峰计数。正电子湮没峰计数与正电子出射钨球表面的数量及湮没光子在钨球中的衰减有关。5MeV的电子／正电子在钨中的射程约为0.20cm，5MeV的光子在钨中的平均自由程为1.26cm。当钨球半径小于正电子在钨中的射程（0.2cm）时，湮没光子在钨球中的衰减可忽略，正电子湮没峰计数主要受正电子出射钨球表面的数量影响，两种模式计算得到的正电子湮没峰计数差别较大。当钨球半径大于正电子在钨中的射程（0.2 cm）时，两种模式计算得到的正电子湮没峰计数差别小于10%。钨球半径越大，正电子湮没峰计数差异越小。因此，可通过对计算模型预先试算，比较两种计算模式下的正电子湮没峰计数，若计数偏差不超过10%，可认为用TTB代替压缩历史算法的计算结果是合理的。

当材料尺寸小于电子在材料中的射程时，TTB模型与压缩历史算法的计算结果差异很大，这是因为此时不考虑电子在材料中的输运是不合理的。由于此时压缩历史算法的计算速度并不慢，不采用TTB模型也能得到计算结果。当材料尺寸大于电子在材料中的射程时，TTB模型与压缩历史算法的计算结果基本一致，这种情况下，用压缩历史算法计算特别耗时，采用TTB模型计算效果更好。

4 下一步工作

可考虑在简化TTB模型中加入轫致辐射光子的空间分布和角度分布。模拟流程如图9所示。

可通过电子轫致辐射后的角度分布抽样确定光子方向，根据光子方向和轫致辐射平均自由程计算得到光子位置，认为电子方向在轫致辐射过程中保持不变。

5 结论

本文根据TTB模型，在PHEN中加入了TTB模块，通过与MCNP中轫致辐射截面及电子射程的对比，验证了TTB模块模拟电子输运的正确性；分析了PHEN加入TTB模块后的模拟计算速度、TTB模块生成的轫致辐射光子能量及TTB模块模拟电子输运的适用范围。正电子湮没峰计数可用于判断TTB模型计算结果的合理性。计算结果表明：当材料尺寸小于电子在材料中的射程时，TTB模型不考虑电子输运的计算结果是不合理的；当材料尺寸大于电子在材料中的射程时，TTB模型计算速度快，计算结果与压缩历史算法的计算结果吻合较好。下一步可根据轫致辐射的空间及角度分布对简化TTB模型进行进一步完善，同时还可加入轫致辐射能量偏倚抽样等相关功能。

［1］DONSKOY E N，ZALYALOV A N.Bremsstrahlung account in photon transport［C］／／Book of Abstracts of International Conference“Parallel Computational Fluid Dynamics”.Moscow，2003：249 252.

［2］X－5MONTE CARLO TEAM.MCNP－A general Monte Carlo N－particle transport code：Version 5［R］.LA－CP－03－0284.Los Alamos National Laboratory，2003.

［3］SALVAT F，FERNANDEZ－VAREA J M，SEMPAU J.PENELOPE－2008：A Code System for Monte Carlo Simulation of Electron and Photon Transport［M］.Barcelona：NZS.OZCD，2009.

［4］SALVAT F，FERNANDEZ－VAREA J M，SEMPAU J，et al.Monte Carlo simulation of bremsstrahlung emission by electrons［J］.Radiation Physics and Chemistry，2006，75（10）：1 201 1 219.

［5］ADAMS K J.Electron upgrade for MCNP 4B［R］.LA－UR－00－3581.Los Alamos National Laboratory，2000.

［6］HALBLEIB J A，KENSEK R P.ITS：The integrated TIGER series of electron／photon transport codes：Version 3.0［J］.IEEE Trans Nucl Sci，1992，39（4）：1 025 1 030.

［7］HAUG E，NAKEL W.The Elementary Process of Bremsstrahlung［M］.New Jersey：World Scientific，2004.

［8］卢希庭.原子核物理［M］.修订版.北京：原子能出版社，2000.（LU Xi－ting.Nuclear Physics［M］.Rev.ed.Beijing：Atomic Energy Press，2000.）

The Implementation and Verification of the Thick－Target Bremsstrahlung in the Particle Transportation Simulation

LI Hua－feng，SHENG Wei，HAN Bin，ZHUO Jun，TAO Ying－long
（Northwest Institute of Nuclear Technology Xi'an，710024，China）

To enhance the processing capability of electron transportation simulation with PHEN（the transportation program of photon，hadron，electron，and neutron），the TTB module was compiled for PHEN based on the TTB（thick－target bremsstrahlung）model.The simulating speed of TTB module，the energy of the bremsstrahlung photon generated by TTB module，and the scope of application of TTB module were analyzed.The TTB module was verified by comparing the bremsstrahlung cross－section and electron range with those given by MCNP.The count of positron annihilation peak was calculated to estimate the simulation results.The calculation result shows that TTB model is irrational as it does not include the electron transport when the medium size is less than the electron range in the medium，and the simulation result of TTB module fits well with the simulation result of the compression history method，and the simulating speed of TTB module is faster than that ofthe compression history method when the medium size is greater than the electron range in the medium.

thick－target bremsstrahlung；Monte Carlo；particle transport；positron annihilation peak

O571.1

A

2095 6223（2015）03 155 08

2015 03 23；

2015 05 19

李华锋（1989－），男，湖北松滋人，助理工程师，硕士，主要从事电子、光子与物质相互作用的蒙特卡罗模拟方法研究。

E－mail：lihuafeng＠nint.ac.cn