HTML
-
FLUKA是一款广泛应用于粒子加速器辐射源项和屏蔽计算的蒙特卡罗程序包,它包含了详细的强子和原子核相互作用模型、强子和电磁过程之间的耦合等物理作用原理,可以用来计算带电粒子输运和带电粒子与物质的相互作用[6]。
屏蔽参数计算采用的模型如图1所示,400 MeV/u碳离子束流打圆柱形靶体,以打靶位置为球心在半径30至36 m处设置球壳状屏蔽体。靶体与屏蔽体内表面之间及屏蔽体外的空间均设置为真空。整个屏蔽体被划分为多层来减小模拟计算的统计方差[10-11],同时也便于采用USRBDX选项卡统计不同屏蔽层的粒子信息。如图1中红色部分所示,取自球心向外延伸的圆锥与屏蔽体球壳相交的范围作为统计区域来统计发射角0°~10°、20°~30°、40°~50°、60°~70°、80°~90°、120°~130°、170°~180°范围内的粒子信息。
在使用FLUKA程序计算时,调用DPMJET、IONTRANS和LOW-NEUT选项卡模拟重离子和低能中子的输运,调用PHYSICS选项卡对重离子核反应中的蒸发碎片、电磁离解等过程进行设置,中子及其他粒子的输运阈值设置为10–3 eV。设置了AMB74选项卡来计算周围剂量当量,它采用的是ICRP74号出版物和Pelliccioni等[12-13]的报告中关于高能中子部分的通量-剂量转换系数。
靶的材料选择铁和水,靶厚度略大于400 MeV/u碳离子在该靶材中的射程,以确保打靶可以最大限度地产生次级辐射。靶体的材料、尺寸以及利用SRIM软件[14]算得的碳离子射程如表1所列。
靶体材料 半径/cm 厚度/cm 射程/cm 铁 3 6 5.093 水 15 30 27.867 混凝土造价便宜、屏蔽适应性广、与土建契合性高,是加速器屏蔽设计中最常用的材料。铁和铅等重材料也经常被用来与混凝土结合使用,以实现局部屏蔽的功能。本次计算选用了密度2.3 g/cm3的混凝土、7.874 g/cm3的铁及11.35 g/cm3的铅作为屏蔽材料。混凝土材料成分如表2所列。
元素 质量分数/% 元素 质量分数/% H 1.00 Al 3.39 C 0.10 Si 33.70 O 52.91 K 1.30 Na 1.60 Ca 4.40 Mg 0.20 Fe 1.40
-
加速器打靶产生次级辐射场具有以下两个特征:是一个类点源场,剂量水平与关注点距打靶位置距离r的平方成反比;剂量随屏蔽厚度的增加呈指数衰减。在束流集中损失的屏蔽计算中,通常使用双参数的指数式(1)来描述周围剂量当量在屏蔽体中的衰减[15]。
式中:H是每个粒子打靶产生的周围剂量当量,单位是Sv·ion–1;Ep是粒子能量;r是粒子打靶位置与统计位置之间的距离,单位是m;θ是束流方向与统计位置的夹角;H0是辐射源项,单位是Sv·m2·ion–1;d是屏蔽层的厚度,单位是cm;ρ是屏蔽体材料密度,单位是g/cm3;λ(θ)是周围剂量当量在屏蔽材料中的衰减长度,单位是g/cm2;α是r方向与屏蔽面法线之间的夹角,g(α)取值为cosα,本次计算使用的模型屏蔽体为球状,因此cosα恒为1。
下文的剂量衰减曲线中,将以H·r2作为纵轴,其单位是Sv·m2·ion–1,用来表征单一粒子带来的辐射剂量大小。以H·r2作为纵轴一方面便于直接观察一定角度范围内周围剂量当量随屏蔽厚度增加的衰减趋势,而不体现与距离r的平方成反比的规律;另一方面也便于对数据进行拟合处理从而得到H0与λ的值。
图4是400 MeV/u碳离子打靶产生的周围剂量当量在屏蔽体中的衰减曲线。从图中可以看出,前向的衰减曲线在刚到达屏蔽体表面时会出现爬高或衰减缓慢的现象。这是由于前向产生的次级粒子中,中高能成分占绝对优势,刚到达屏蔽体表面时会迅速损耗能量产生大量低能中子,从而出现明显的二次中子增殖效应。而侧向及后向的衰减曲线刚开始衰减较快,随着屏蔽厚度增加衰减速度逐渐减慢。这是由于侧向和后向产生的次级粒子中,低能中子占绝对优势,随着屏蔽厚度的增加,较低能量中子所占的比例越来越少,能谱逐渐硬化,最终达到能谱平衡状态。
从图4可以看出,三种材料的屏蔽效果具有很大差异,这是主要是由其中子反应截面的差异导致的。中高能中子在进入屏蔽材料后,主要通过非弹性散射反应迅速降低能量;当能量低于非弹性散射阈能时,主要通过弹性散射来降低能量;当能量足够低时,主要通过俘获反应被吸收[16]。图5给出了混凝土、铁和铅的中子反应总截面、非弹反应截面、弹性散射截面的对比图。图中纵轴为σ·ρ·NA·M –1,单位是cm–1,用来表征单位厚度的屏蔽材料与中子发生反应的概率大小。其中σ是单一原子核的中子反应截面,单位cm2,作图所用的俘获截面数据均取自ENDF数据库[17];ρ是屏蔽材料的密度,单位g/cm3;NA是阿伏伽德罗常数,单位mol–1;M是屏蔽材料的摩尔质量,单位是g·mol–1。混凝土的σ·ρ·NA·M–1计算方法如式(2)所示,
式中:ni是第i种元素归一化的原子数目占比;σi是第i种元素的中子反应截面;Mi是第i种元素的摩尔质量。
从图5可以看出,对于非弹性散射反应(主要针对几MeV以上能区的中子),铁和铅的反应截面差距不大,混凝土的反应截面要比铁和铅的小一个量级左右。对于弹性散射反应,铁的反应截面整体上明显大于铅的,除过中子能量约10–3 eV以下时铁的反应截面略小于混凝土的,其余能量段均大于混凝土的;在10–1 eV以下能区铅的反应截面明显小于混凝土的,在10–1 eV至105 eV能区铅和混凝土的反应截面基本一致,105 eV以上时,铅的反应截面比混凝土的大。对于俘获反应(主要针对低能中子),铁的反应截面比铅的高一个量级左右,铅的反应截面比混凝土的高一倍左右。整体来看,铁在各个能量段的中子反应截面都要大于铅和混凝土的;铅的非弹性散射反应截面明显大于混凝土的,其低能中子俘获反应截面也稍大于混凝土的,但其弹性散射截面比混凝土的小;另外铅、铁等重材料对俘获反应伴生的
$\gamma $ 射线的屏蔽能力也要比混凝土的更强[13]。 -
在混凝土、铁屏蔽体中,前向的周围剂量当量在屏蔽厚度0到100 cm时衰减较为缓慢,屏蔽厚度达到100 cm以后衰减速度会增快。结合Agosteo等[18-19]之前的研究中使用的拟合方法,选用屏蔽厚度100 cm以上的范围对混凝土和铁屏蔽体的前项数据进行拟合,将式(1)简化为拟合式(3):
在混凝土、铁屏蔽体中,侧向和后向周围剂量当量在屏蔽厚度0到100 cm时衰减较快,100 cm以上时衰减变慢。如果仅利用式(3)进行拟合,会造成辐射源项的严重低估[10]。因此采用双指数式(4)对混凝土、铁屏蔽体的侧、后向数据进行分段拟合,拟合区间分别选取0到100 cm以及100到600 cm。
铅屏蔽体厚度达到400 cm以上时,周围剂量当量的衰减趋势加快,但考虑到厚度400 cm以上时周围剂量当量非常低,实际屏蔽设计工作中几乎不可能采用厚度400 cm以上的铅屏蔽体,故选取100 到400 cm的屏蔽厚度,采用拟合式(5)对铅屏蔽体中的衰减数据进行拟合。
400 MeV/u碳离子打靶产生的次级辐射场具有很强的前冲性,在束流正前向(0°方向)屏蔽厚度估算时,如果采用0°~10°数据的拟合参数,会对辐射源项值H0产生约2到5倍的低估,因此补充了发射角0°~1°的相关数据。根据前文提到的拟合方法,使用Origin软件对衰减曲线进行了拟合,拟合结果如表3至表8所列。表中,式(4)的拟合结果列入H1、λ1、H0、λ0四栏;式(3)和式(5)的拟合结果列入H0、λ0两栏;
$R_0^2 $ 、$R_1^2 $ 分别表示两次拟合数据的相关系数。发射角 H1/(Sv·m2·ion–1) λ1/(g·cm–2) R12 H0/(Sv·m2·ion–1) λ0/(g·cm–2) R02 拟合公式 0°~1° — — — (2.48±0.08)×10–12 121.05±0.74 0.997 1 (3) 0°~10° — — — (9.10±0.27)×10–13 125.41±0.77 0.997 1 (3) 10°~20° — — — (3.04±0.09)×10–13 124.73±0.75 0.996 8 (3) 20°~30° — — — (1.34±0.03)×10–13 120.42±0.53 0.998 4 (3) 40°~50° (5.28±0.16)×10–14 95.24±2.09 0.993 8 (4.02±0.02)×10–14 110.00±0.09 0.999 9 (4) 60°~70° (3.98±0.20)×10–14 67.47±2.08 0.986 7 (1.45±0.04)×10–14 99.01±0.55 0.998 3 (4) 80°~90° (3.40±0.26)×10–14 49.85±1.91 0.977 2 (5.60±0.28)×10–15 87.99±0.89 0.994 0 (4) 120°~130° (3.08±0.17)×10–14 34.72±0.91 0.991 0 (1.09±0.11)×10–15 76.44±1.47 0.981 2 (4) 170°~180° (2.97±0.16)×10–14 31.61±0.77 0.9917 (4.44±0.61)×10–16 78.77±1.86 0.968 4 (4) 发射角 H1/(Sv·m2·ion–1) λ1/(g·cm–2) R12 H0/(Sv·m2·ion–1) λ0/(g·cm–2) R02 拟合公式 0°~1° — — — (3.25±0.23)×10–11 162.99±1.44 0.990 4 (3) 0°~10° — — — (9.55±0.22)×10–12 172.15±0.54 0.999 3 (3) 10°~20° — — — (3.22±0.06)×10–12 172.26±0.44 0.999 6 (3) 20°~30° — — — (1.37±0.02)×10–12 169.44±0.36 0.999 7 (3) 40°~50° — — — (3.58±0.15)×10–13 164.59±0.85 0.997 5 (3) 60°~70° (1.83±0.05)×10–13 154.91±1.66 0.997 5 (1.17±0.10)×10–13 162.75±1.85 0.992 5 (4) 80°~90° (1.57±0.02)×10–13 137.25±0.62 0.999 5 (4.89±0.57)×10–14 164.42±2.58 0.987 9 (4) 120°~130° (1.39±0.01)×10–13 128.39±0.46 0.999 7 (2.23±0.24)×10–14 173.86±2.61 0.990 8 (4) 170°~180° (1.35±0.02)×10–13 127.89±0.60 0.999 5 (2.55±0.23)×10–14 169.01±2.00 0.987 1 (4) 发射角 H1/(Sv·m2·ion–1) λ1/(g·cm–2) R12 H0/(Sv·m2·ion–1) λ0/(g·cm–2) R02 拟合公式 0°~1° — — — (1.01±0.04)×10–10 340.53±2.09 0.995 9 (5) 0°~10° — — — (3.53±0.09)×10–11 356.13±1.53 0.999 0 (5) 10°~20° — — — (1.28±0.02)×10–11 355.47±1.09 0.999 5 (5) 20°~30° — — — (5.79±0.05)×10–12 354.02±0.49 0.999 9 (5) 40°~50° — — — (1.97±0.05)×10–12 351.39±1.59 0.999 3 (5) 60°~70° — — — (1.12±0.04)×10–12 347.09±1.95 0.998 5 (5) 80°~90° — — — (8.43±0.39)×10–13 345.09±2.47 0.996 0 (5) 120°~130° — — — (7.50±0.34)×10–13 343.11±2.45 0.995 9 (5) 170°~180° — — — (7.08±0.34)×10–13 345.72±2.73 0.997 6 (5) 发射角 H1/(Sv·m2·ion–1) λ1/(g·cm–2) R12 H0/(Sv·m2·ion–1) λ0/(g·cm–2) R02 拟合公式 0°~1° — — — (7.28±0.23)×10–12 115.58±0.70 0.995 9 (3) 0°~10° — — — (1.67±0.05)×10–12 124.19±0.78 0.996 7 (3) 10°~20° — — — (3.58±0.13)×10–13 125.55±0.85 0.995 9 (3) 20°~30° — — — (1.49±0.04)×10–13 122.02±0.63 0.998 0 (3) 40°~50° — — — (3.72±0.04)×10–14 111.16±0.26 0.999 5 (3) 60°~70° (1.94±0.04)×10–14 80.73±0.99 0.998 1 (1.14±0.03)×10–14 98.33±0.45 0.998 6 (4) 80°~90° (1.22±0.05)×10–14 59.85±1.37 0.993 5 (3.52±0.17)×10–15 89.11±0.85 0.995 2 (4) 120°~130° (9.40±0.57)×10–15 40.72±1.22 0.988 0 (9.46±0.68)×10–16 74.82±0.96 0.988 9 (4) 170°~180° (8.00±0.69)×10–15 37.47±1.39 0.979 0 (4.17±0.43)×10–16 84.50±1.70 0.978 6 (4) 发射角 H1/(Sv·m2·ion–1) λ1/(g·cm–2) R12 H0/(Sv·m2·ion–1) λ0/(g·cm–2) R02 拟合公式 0°~1° — — — (7.30±0.27)×10–11 164.94±0.78 0.997 6 (3) 0°~10° — — — (1.65±0.03)×10–11 172.00±0.43 0.999 5 (3) 10°~20° — — — (4.05±0.06)×10–12 172.83±0.39 0.999 7 (3) 20°~30° — — — (1.54±0.03)×10–12 170.14±0.40 0.999 6 (3) 40°~50° — — — (3.67±0.12)×10–13 162.72±0.74 0.998 7 (3) 60°~70° (1.24±0.03)×10–13 157.17±1.75 0.997 8 (8.85±0.78)×10–14 161.92±1.72 0.987 6 (4) 80°~90° (1.02±0.00)×10–13 137.15±0.30 0.999 9 (2.76±0.24)×10–14 169.48±2.04 0.991 2 (4) 120°~130° (9.53±0.08)×10–14 129.15±0.5 0.999 7 (1.96±0.21)×10–14 169.30±2.62 0.991 4 (4) 170°~180° (9.59±0.08)×10–14 128.87±0.39 0.999 8 (1.85±0.25)×10–14 170.65±2.98 0.985 3 (4) 发射角 H1/(Sv·m2·ion–1) λ1/(g·cm–2) R12 H0/(Sv·m2·ion–1) λ0/(g·cm–2) R02 拟合公式 0°~1° — — — (2.07±0.07)×10–10 351.5±2.32 0.998 9 (5) 0°~10° — — — (6.32±0.14)×10–11 353.8±1.28 0.999 4 (5) 10°~20° — — — (1.56±0.04)×10–11 356.25±1.51 0.998 7 (5) 20°~30° — — — (6.14±0.12)×10–12 356.47±1.18 0.999 5 (5) 40°~50° — — — (1.96±0.05)×10–12 348.91±1.42 0.998 6 (5) 60°~70° — — — (1.10±0.06)×10–12 341.76±2.88 0.994 1 (5) 80°~90° — — — (7.41±0.13)×10–13 349.12±1.03 0.999 7 (5) 120°~130° — — — (8.47±0.40)×10–13 339.41±2.30 0.996 1 (5) 170°~180° — — — (5.95±0.40)×10–13 361.46±2.97 0.997 8 (5)
4.1. 周围剂量当量的衰减曲线
4.2. 屏蔽参数的获取
-
福建莆田市即将建造一台医用重离子加速器,下面以该装置水平治疗室的屏蔽设计为例,介绍本文给出的屏蔽参数的使用方法。该装置加速粒子为12C,最大能量为400 MeV/u,治疗室内最大束流损失参数为1×108 pps。治疗室的FLUKA计算模型如图6所示,屏蔽体材料选用密度2.3 g/cm3的混凝土,靶体为厚度30 cm的水靶。
FLUKA模拟计算得到的治疗室周围剂量当量率分布如图7所示。蒙特卡罗计算结果显示,当屏蔽体外周围剂量当量率限值设定为2.5 μSv/h时,图6中A位置所在方向需屏蔽体厚度438.21 cm;B位置所在方向(约45°发射角处)是侧向场强最大的位置,所需屏蔽体厚度为168.31 cm;C位置所在方向需屏蔽体厚度128.91 cm;D位置所在方向所需屏蔽体厚度70.6 cm。
下面采用本文得到的屏蔽参数计算所需的厚度。周围剂量当量率与单个粒子带来的周围剂量当量的换算关系如式(6)所示。式中,H*是周围剂量当量率,单位是μSv/h;H是单个粒子产生的周围剂量当量,单位是Sv·ion–1;I是平均束流损失,设计值是1×108 pps。当H*=2.5 μSv/h时,计算得到H=6.944 4×10–18 Sv·ion–1。
选取图5中A、B、C、D四个位置进行计算。A位置选用0°~1°的屏蔽参数,将下列数据代入式(3)中:
计算得到A位置所需屏蔽体厚度d为447.13 cm。
与A位置同样的方法,选用40°~50°的屏蔽参数进行计算,计算得到B位置所需的屏蔽体厚度d为170.29 cm。
C、D位置的计算方法与A位置稍有不同,选用80°~90°、170°~180°的数据进行计算。计算时先通过H1、λ1(厚度100 cm以下的拟合数据)计算得到厚度d1,若d1<100 cm,则d1即为所求厚度;若d1>100 cm,则使用H0、λ0(厚度100 cm以上的拟合数据)进行计算,得到的厚度d0即为所求厚度。经计算,C位置所需的屏蔽厚度为126.25 cm,D位置所需屏蔽厚度为70.36 cm。
使用FLUKA直接模拟计算和通过本文给出的屏蔽参数得到的结果如表9所列。可以看出两种方法的计算结果符合很好,相对误差小于3%,说明本文给出的屏蔽参数是可靠的,能够直接应用于医用重离子加速器的屏蔽设计工作。
位置 厚度/cm(模拟结果) 厚度/cm(屏蔽参数估算) 相对误差/% A 438.21 447.13 2.04 B 168.31 170.29 1.18 C 128.91 126.25 2.06 D 70.60 70.36 0.34