4.12. 例題: 包括的なDCHAIN-PHITS計算¶
本例題では、DCHAINの主要な機能を一連のDCHAIN-PHITS計算を通して説明する。使用するPHITS入力ファイルは<PHITS-install>/dchain-sp/sample/simple/phits_simple.inpにあり、付録6.1 にも掲載している。この入力からDCHAIN入力ファイルを自動生成する。
直径1 cmのペンシル状250 MeV陽子ビームを、タングステン、水、鉄の3区間で構成される円柱状ロッドへ入射する。ロッドの側面と後端はベリリウムで覆われている。ビーム電流100 nAで6分間照射した後、出力を50 %（50 nA）へ下げて4分間照射し、その後50分間冷却する。
PHITS入力パラメータicntlを一時的にicntl= 8とすると、図 4.12.1 の計算体系を描画できる。
[ T i t l e ]
sample input for dchain: 250MeV protons incident on W/H2O/Fe rod surrounded in Be
[ P a r a m e t e r s ]
icntl = 0 # (D=0) 3:ECH 5:NOR 6:SRC 7,8:GSH 11:DSH 12:DUMP
maxcas = 1000 # (D=10) number of particles per one batch
maxbch = 10 # (D=10) number of batches
igamma = 3 # (D=0) 0:No, 1:Old, 2:EBITEM, 3:EBITEM+Isomer
file(1) = c:/phits # (D=c:/phits) PHITS install folder name
$ required options for DCHAIN
jmout = 1 # (D=0) Density echo, 0:input, 1:number density
e-mode = 0 # (D=0) Event generator mode is not recommended for DCHAIN calculation
[ S o u r c e ]
s-type = 1 # mono-energetic axial source
proj = proton # kind of incident particle
e0 = 250.00 # energy of beam [MeV/u]
r0 = 0.5000 # radius [cm]
x0 = 0.0000 # (D=0.0) center position of x-axis [cm]
y0 = 0.0000 # (D=0.0) center position of y-axis [cm]
z0 = -10.000 # minimum position of z-axis [cm]
z1 = -10.000 # maximum position of z-axis [cm]
dir = 1.0000 # z-direction of beam [cosine]
[ M a t e r i a l ]
$ Tungsten, Density (g/cm3) = 19.300000
M1 74000 1.000000 $ W
$ Water, Density (g/cm3) = 0.998207 ~= 1.0
M2 1000 0.666657 $ H
8000 0.333343 $ O
$ Iron, Density (g/cm3) = 7.874000
M3 26000 1.000000 $ Fe
$ Beryllium, Density (g/cm3) = 1.848000
M4 4000 1.000000 $ Be
[ S u r f a c e ]
set:c1[3] $ thickness of W portion
set:c2[10] $ thickness of H2O portion
set:c3[5] $ thickness of Fe portion
set:c4[5] $ thickness of Be back end
set:c5[4] $ radius of inner cylinder
set:c6[8] $ radius of outer cylinder
10 pz 0
11 pz c1
12 pz c1+c2
13 pz c1+c2+c3
14 pz c1+c2+c3+c4
20 cz c5
21 cz c6
99 so 80.0
[ C e l l ]
1 1 -19.3 10 -11 -20 $ W
2 2 -1.0 11 -12 -20 $ H20
3 3 -7.874 12 -13 -20 $ Fe
4 4 -1.848 (10 -13 20 -21):(13 -14 -21) $ Be
99 0 (-99 -10):(-99 14):(-99 21) $ vacuum
999 -1 99 $ particle graveyard
[ V o l u m e ]
$ required section for DCHAIN
reg vol
1 pi*c1*c5**2
2 pi*c2*c5**2
3 pi*c3*c5**2
4 pi*(((c1+c2+c3+c4)*c6**2)-((c1+c2+c3)*c5**2))
[ M a t N a m e C o l o r ]
mat name size color
0 void 1.00 lightgray
1 Tungsten 1.00 darkgreen
2 Water 1.00 pastelcyan
3 Iron 1.00 orange
4 Beryllium 1.00 pastelviolet
[ T - D C H A I N ]
$ required section for DCHAIN
title = 250MeV protons incident on Be-encased W/H2O/Fe column
mesh = reg # mesh type is region-wise
reg = 1 2 3 4
file = spallation.out # file name of dchain-phits input file
timeevo = 3 # time evolution / irradiation schedule
6.0 m 1.0
4.0 m 0.5
50.0 m 0.0
outtime = 8 # output times
1.0 m
3.0 m
6.0 m
10.0 m
15.0 m
20.0 m
30.0 m
40.0 m
$ beam current (nA)
set:c11[100.0]
$ beam power (source/sec)
set:c12[ c11 * 1.0e-9 / (1.602177e-19) ]
amp = c12 # (D=1.0) Source Intensity(source/sec)
[ T - T r a c k ]
title = proton and neutron distributions
mesh = xyz
x-type = 2
xmin = -10
xmax = 10
nx = 1
z-type = 2
zmin = -10
zmax = 30
nz = 200
y-type = 2
ymin = -20
ymax = 20
ny = 200
e-type = 2
emin = 0
emax = 2000
ne = 1
axis = yz
file = yz-track.out
part = proton neutron
unit = 1
2D-type = 7 # 1:Cont, 2:Clust, 3:Color, 4:xyz, 5:mat, 6:Clust+Cont, 7:Col+Cont
gshow = 3 # 0: no 1:bnd, 2:bnd+mat, 3:bnd+reg 4:bnd+lat
epsout=1
[ E n d ]
図 4.12.1 例題のPHITS計算体系¶
本例題では、タングステン部の核破砕反応による放射化と、鉄部の中性子反応による放射化の両方を扱う。
icntl= 0へ戻して通常モードでPHITSを実行すると、図 4.12.2 から 図 4.12.3 に示す陽子・中性子フラックス分布が得られる。一次陽子ビームは水領域内で飛程に達し、鉄領域までは到達しない。
このPHITS計算により、[T-Dchain]のfileで指定したDCHAIN入力ファイルspallation.outが生成される。この内容は 付録6.2 にも掲載している。併せて、中性子フラックスを記録したspallation.dtrk、20 MeV未満の中性子反応を除く核種生成収率を記録したspallation.dyld、DCHAINデータライブラリフォルダーのパスを渡すdch_link.datが生成される。PHITS実行時に、既定の簡略形式iredufmt= 1ではなく従来形式iredufmt= 0を使用した場合は、核種生成収率の統計誤差を記録するspallation_err.dyldも生成される。
htitle = 250MeV protons incident on Be-encased W/H2O/Fe column
! --- control parameters ---
imode = 2
jmode = 2
! --- calculation parameters ---
idivs = 50
ichain = 100
itdecs = 1
itdecn = 1
isomtr = 2
ifisyd = 0
ifisye = 0
ixsrall = 1
! --- data library parameters ---
inxslib = 100
idcylib = 5
! hnxslib = h-JP-US-EU
! hdcylib = J15-EB8-ES
! --- output parameters ---
iyild = 2
iggrp = 3
ibetap = 1
acmin = 1.0000E-20
istabl = 0
igsdef = 1
igsorg = 1
iwrtchn = 1
chrlvth = -1.0000E+00
iwrchdt = 0
iwrchss = 0
iwrchnuc = 0
idosecf = 1
idosunit = 333
imtcard = 1
imtcnum = 0
imtcmeta = 2
thmatnd = 1.0000E-06
! --- Proton beam current and neutron flux ---
amp = 1.0000E-04
ebeam = 3.0000E+00
prodnp = 1.0000E+00
! --- irradiation time ---
itstep = 3
6.0000E+00 m 1.0000E+00
4.0000E+00 m 5.0000E-01
5.0000E+01 m 0.0000E+00
! --- output time ---
itout = 8
1.0000E+00 m
3.0000E+00 m
6.0000E+00 m
1.0000E+01 m
1.5000E+01 m
2.0000E+01 m
3.0000E+01 m
4.0000E+01 m
! --- irradiation condition ---
inmtcf = 1
! ffnmtcf = 1.0
itdchout = 1
hhnmtcf = spallation.dyld
iertdcho = 1
henmtcf = spallation_err.dyld
iredufmt = 1
irdonce = 1
! --- calculation region details ---
iregon = 4
!1)HRGCMM 2)IREGS 3)ITGNCLS 4)FLUXS 5)HNFLUXS 6)VOLUMES
DUMMY001 1 5 2.3797E+10 spallation.dtrk 1.5080E+02
W-180 7.5882E-05
W-182 1.6756E-02
W-183 9.0514E-03
W-184 1.9376E-02
W-186 1.7975E-02
!1)HRGCMM 2)IREGS 3)ITGNCLS 4)FLUXS 5)HNFLUXS 6)VOLUMES
DUMMY002 2 3 7.3851E+09 spallation.dtrk 5.0265E+02
H-1 6.6880E-02
H-2 1.0033E-05
O-16 3.3446E-02
!1)HRGCMM 2)IREGS 3)ITGNCLS 4)FLUXS 5)HNFLUXS 6)VOLUMES
DUMMY003 3 4 9.4207E+08 spallation.dtrk 2.5133E+02
Fe-54 4.9642E-03
Fe-56 7.7928E-02
Fe-57 1.7997E-03
Fe-58 2.3951E-04
!1)HRGCMM 2)IREGS 3)ITGNCLS 4)FLUXS 5)HNFLUXS 6)VOLUMES
DUMMY004 4 1 1.6216E+09 spallation.dtrk 3.7196E+03
Be-9 1.2352E-01
! --- output option for ANGEL ---
phitsout = 1
angelout = 1
iangbpwr = 1
! angelout_act = 0
! angelout_heat = 0
! angelout_heat_alpha = 0
! angelout_heat_beta = 0
! angelout_heat_gamma = 0
! angelout_dose = 0
! angelout_region = all, sum
! angelout_nuclides = 0
end
生成されたspallation.outは修正せずにDCHAINへ入力できる。DCHAINの出力ファイルは、いずれもspallation.*という基底名に用途別の拡張子を付けて生成される。以下に各ファイルを説明する。
推奨されるバッチ／シェルスクリプトでDCHAINを実行すると、DCHAIN終了後にANGELが自動実行され、放射能や崩壊熱などの時間変化が描画される。図 4.12.4 に放射能のプロットを示す。この図はANGEL入力ファイルspallation.angから生成されたspallation.epsに基づく。図を変更する場合は、spallation.angを編集してANGELを再実行する。照射中はベリリウム領域（セル4）が全放射能へ大きく寄与するが、その寄与は照射終了後すぐに相対的に小さくなる。
図 4.12.4 各領域の放射能の時間変化（ANGELで作成）¶
DCHAIN入力ファイルの末尾には、*.angへ書き出す内容を制御するパラメータがある。IANGBPWRは、ビーム出力履歴を図へ重ねるかどうかを指定する。図 4.12.4 の赤い網掛け部分がこれに対応する。ANGELOUT_REGIONは出力する領域を指定し、ANGELOUT_NUCLIDESは領域別データとともに描画する個別核種を指定する。
たとえば、鉄領域（セル3）の\(^{56}\)Mn、\(^{53}\)Fe、\(^{53m}\)Feだけに注目する場合は、末尾のパラメータを次のように変更する。結果を 図 4.12.5 に示す。
angelout_region = 3
angelout_nuclides = 3
Mn-56
Fe-53
Fe-53m
ANGELOUT_NUCLIDESで指定した核種の値は、[T-Dchain]のregで指定した全領域の合計であり、ANGELOUT_REGIONで指定した領域だけの値ではないことに注意すること。
図 4.12.5 領域3の放射能と指定核種の放射能の時間変化¶
最も詳細な主要出力は*.actファイルであり、本例ではspallation.actである。各領域・時間ステップの核種インベントリ、放射能、崩壊熱、光子スペクトルなどを確認できる。一方、全放射能への寄与が無視できる核種も多数含まれるため、ファイルが長くなることがある。出力核種のしきい値はACMINで制御する。PHITSの既定値は\(10^{-20}\)であり、放射能が\(10^{-20}\)Bq以上の核種を含める。これは極めて小さいしきい値であるため、通常は生成した核種のほぼすべてが出力される。
*.actファイルの「Top 10」には、放射能、崩壊熱、光子線量への寄与が大きい10核種がまとめられる。ファイル全体から寄与の小さい核種を除外したい場合は、たとえばACMIN= -1.0E+4とし、原子数密度が\(10^4\)atoms/cm\(^3\)以上の核種だけを出力する。本例では、同程度の収率を持つ多数の核破砕生成物が存在するため、タングステン領域の記載量はあまり減らない。一方、鉄領域（領域3）では、まれな中性子反応で生成される核種が除外され、1分時点の出力核種数が45から17へ減少する。
後述するDCHAIN Toolsで*.actファイルを解析する場合、通常は人間が直接読む際の長さを気にする必要はない。ただし、領域数や出力時刻数が多く、ファイル容量自体を抑える必要がある計算では、ACMINによる制限が有効である。
*.actファイルの「Top 10」と同様に、spallation.alrには、各出力時刻outtime/ITOUTにおいて、全領域で積算した全放射能と全崩壊熱への寄与が大きい上位30核種が記載される。spallation.yldには、各領域における各核種の最終インベントリが記載される。ただし、単位体積当たりの全原子数に対する割合が\(10^{-12}\)以上の核種だけを対象とする。IYILDは、この全原子数に初期標的核種のインベントリを含めるかどうかを指定し、既定では含める。
DCS（崩壊系列スキーム）ファイルspallation.dcsには、各時間ステップにおける各核種のインベントリ変化と、その変化を生じさせた崩壊系列が記載される。既定では変化量と系列だけを示す。IWRCHDT= 1（既定値0）とすると、系列内の各過程がインベントリへ与える寄与と、その起源が高エネルギー反応か、崩壊・中性子反応かを詳細に出力する。
照射期間では、既定でIDIVSサブステップの最後に生じた変化だけを出力する。IWRCHSS= 1（既定値0）とすると、すべての照射サブステップの変化を出力する。両機能、特にIWRCHSS= 1はDCSファイルを大きくするが、併用すればDCHAINが計算した各インベントリ変化の原因を詳しく追跡できる。
CHRLVTHは、単一の崩壊・反応系列によるインベントリ変化がどの程度以上であればDCSファイルへ出力するかを制御する。詳細な追跡が必要な核種が限られる場合は、IWRCHNUCで対象核種だけを指定できる。
本例では、DCS出力を制御する入力を次のように変更し、鉄領域（領域3）で生成される8核種に関係する系列だけを表示できる。さらに、各系列のどの過程がインベントリへ寄与したかを詳細に出力する。ここではIWRCHSS= 0のままとする。通常は各照射期間の最後のサブステップだけでも、系列内の各過程の相対寄与を十分に把握できる。IWRCHSS= 1が必要となるのは、照射期間内に標的核種のインベントリが大きく変化すると予想される場合や、すべてのインベントリ変化を完全に追跡する場合である。
iwrtchn = 1
chrlvth = -1.0000E+00
iwrchdt = 1 ! write more detailed information on each link in decay chain file
iwrchss = 0
! add special "iwrchnuc" parameter for tracking of specific nuclides in this DCS file
iwrchnuc = 8
Cr-49
Mn-56
Fe-53
Cr-55
Cr-51
V-48
Mn-54
Fe-55
この入力で再計算すると、指定核種を含まない他の3領域には核種情報が出力されず、各崩壊系列には生成機構別の寄与を表す3行が追加される。DCHAINは構築した系列全体を表示するが、多くの場合、対象核種のインベントリへ大きく寄与するのは系列末端付近の数過程だけである。
spallation.phtには、各時間ステップ・領域の光子スペクトルを表すPHITS形式の[Source]セクションが記載される。既定では全時間ステップを1ファイルへ出力するが、PHITSのinfl:機能で挿入ファイルとして使用する場合などは、PHITSOUTを1から2へ変更し、出力時刻ごとにファイルを分割できる。
この出力は、冷却時間や遮蔽条件を変えた二次線量評価に利用できる。ただし、各領域内の放射能分布は一様と仮定されるため、大きな領域では実際の分布を適切に表さない場合がある。二次評価を予定している場合は、タリー領域を小さく分割するかxyzメッシュを使用し、必要な空間分解能を確保すること。
従来機能として、spallation.gsdには各時間ステップ・領域の光子スペクトルを表すMCNP形式のSDEF（線源定義）カードが記載される。spallation.gsoには、各光子エネルギー群に含まれる核種と、その群への相対寄与が領域・出力時刻ごとに記載される。
spallation.matには、各領域・出力時刻の組成を表すPHITS形式の[Material]セクションが記載される。この材料カードはMCNPとも互換性がある。[Source]セクションと同様に、既定では全時刻を1ファイルへ出力するが、IMTCARDで時間ステップごとに分割できる。
DCHAINは、標的核種の減少によって核破砕反応による生成率が低下する標的燃焼の効果を近似的に扱う。材料組成がほとんど変化しない一般的な計算では、この近似の影響は小さい。一方、組成が大きく変化する計算で高精度な結果を得るには、PHITSとDCHAINを交互に繰り返す反復計算が必要である。特に、中性子フラックススペクトルへ影響する核種が蓄積する中性子場では、DCHAINが一定と仮定するフラックス自体が時間とともに変化し得る。その場合、spallation.matはPHITS-DCHAIN反復燃焼計算の入力更新を簡略化する。設定の詳細は カード4 を参照すること。
spallation.lstには、DCHAIN実行時の端末出力の大部分が転送される。入力エコー、入力パラメータの説明、警告、診断情報などが含まれる。
DCHAIN Toolsモジュールを使用すると、DCHAIN出力を短時間で解析できる。現在はPHITS Toolsのサブモジュールとして提供されており、pip install PHITS-Toolsによる導入を推奨する。詳細は DCHAIN Tools documentation を参照すること。旧リポジトリには PDFマニュアルと出力辞書構造の概要 も用意されている。以下のPythonスクリプトで本例題の出力を読み込み、結果の一部を確認する。併せて実行時の端末出力も示す。
from dchain_tools import *
import numpy as np
import matplotlib.pyplot as plt
# path to the folder containing output files
simulation_folder_path = r'C:\phits\lecture\advanced\DCHAIN3\\'
# the filename part before the extension of all of the output files
simulation_base_name = 'spallation'
# Provide this path information to the main DCHAIN parsing function
# It parses *.act and, if found, *.dtrk and *.dyld too (and *.dcs if enabled).
dchain_output = process_dchain_simulation_output(simulation_folder_path,simulation_base_name,process_DCS_file=False)
# Output can be accessed in 'dictionary' or 'class/attribute' styles.
# print list of all nuclides (as text strings) found in the second region over all times
print(dchain_output['nuclides']['names'][1]) # dictionary-style access
# print total activity and its absolute error in the second region and first time step
print (dchain_output.nuclides.total.activity.value[1][0], dchain_output.nuclides.total.activity.error[1][0]) # attribute-style access
# find the activity of Fe53 in the 3rd region at the 5th time step
ri = 2 # region index
ti = 4 # time index
# get desired nuclide name string formatted in DCHAIN's specific 6-character syntax
Fe53_Dname = nuclide_plain_str_to_Dname('Fe-53')
# determine index of desired nuclide among all nuclides in this region
Fe53_index = dchain_output['nuclides']['names'][ri].index(Fe53_Dname)
# now extract the activity value
A_Fe53 = dchain_output.nuclides.activity.value[ri][ti,Fe53_index]
print('A(Fe-53) in [Bq/cc] in region 3 at 5th output time: ',A_Fe53)
Parsing DCHAIN activation file... (0.01 seconds elapsed)
Restructuring nuclide data table array... (0.05 seconds elapsed)
['H 3 ', 'He 6 ', 'Li 5 ', 'Li 8 ', 'Li 9 ', 'Be 7 ', 'Be 8 ', 'Be 10 ', 'Be 11 ', 'B 8 ', 'B 9 ', 'B 12 ', 'B 13 ', 'C 10 ', 'C 11 ', 'C 14 ', 'C 15 ', 'N 13 ', 'N 16 ', 'N 17 ', 'N 18 ', 'O 14 ', 'O 15 ', 'O 19 ', 'Hf178n', 'Re188 ', 'Re188m']
15921500.0 563170.0
A(Fe-53) in [Bq/cc] in region 3 at 5th output time: 76742.0
水で構成される領域2の核種一覧に、\(^{178}\)Hfの第2準安定状態と\(^{188}\)Reの2つの核異性体が含まれている。これは不具合ではなく、DCHAIN出力とPHITSが生成した*.dyldの両方に実際に存在する。[T-Dchain]が自動生成する[T-Yield]タリーは、既定のoutput= productではなくoutput= cutoffを使用する。この設定では、粒子・原子核が生成されたセルではなく、そのヒストリーが終了したセルで計数する。放射化計算には後者が必要なためである。通常、この違いはトリトンのような軽い粒子で重要となるが、本例のように、セル境界付近で生成された重い反跳核が隣接セルへ移動する場合にも影響する。
Pythonへ読み込んだ結果は容易に可視化できる。図 4.12.6 から 図 4.12.7 に、先にANGELで作成した放射能プロットの再現と、照射終了時における領域1の光子スペクトルを示す。これらの図を作成するため、先ほどのスクリプトへ追加するコードも以下に示す。
# Recreate ANGEL activity plot
plt.figure(0,(4.5,3.5)) # figure index and dimensions
t_minutes = np.array(dchain_output.time.from_start_sec)/60
total_activity , total_activity_error = 0 , 0
for ri in range(len(dchain_output.region.numbers)): # for each region
region_volume = dchain_output.region.volume[ri]
plt.errorbar(t_minutes,
dchain_output.nuclides.total.activity.value[ri][:]*region_volume,
dchain_output.nuclides.total.activity.error[ri][:]*region_volume,
label='Cell {}'.format(ri+1))
total_activity += dchain_output.nuclides.total.activity.value[ri][:]*region_volume
total_activity_error += (dchain_output.nuclides.total.activity.error[ri][:]*region_volume)**2
total_activity_error = total_activity_error**0.5
# Total activity line
plt.errorbar(t_minutes,total_activity,total_activity_error,
label='All cells',linestyle='--',color='k')
plt.yscale('log')
plt.legend(loc='best')
plt.xlabel('time [min]')
plt.ylabel('activity [Bq]')
plt.title('activity')
plt.grid(b=True, which='major', linestyle='-', alpha=0.25)
plt.grid(b=True, which='minor', linestyle='-', alpha=0.10)
plt.tight_layout()
# Plot photon spectrum in region 1 at 3rd time step
# Note: as in the *.act file, photon spectra are binned by energy in descending order
plt.figure(1,(4,3)) # figure index and dimensions
bin_widths = dchain_output.gamma.spectra.E_upper[0][2,:]-dchain_output.gamma.spectra.E_lower[0][2,:]
plt.bar(x=dchain_output.gamma.spectra.E_lower[0][2,:],align='edge',
width=bin_widths,
height=dchain_output.gamma.spectra.flux.value[0][2,:]/bin_widths,
edgecolor='k',linewidth=0.5)
plt.xscale('log')
plt.yscale('log')
plt.xlabel('photon energy [MeV]')
plt.ylabel(r'$\gamma$ flux [#/(sec$\cdot$cm$^3$$\cdot$MeV)]')
plt.title('Photon flux in cell 1 at t = {}'.format(seconds_to_dhms(dchain_output.time.from_start_sec[2])))
plt.grid(b=True, which='major', linestyle='-', alpha=0.25)
plt.grid(b=True, which='minor', linestyle='-', alpha=0.10)
plt.tight_layout()
plt.show()
DCHAIN Toolsには、時刻や領域ごとに重要な核種を短時間で把握できるよう、核種を自動的に順位付けして描画する機能もある。領域3についての例を 図 4.12.8 に示す。
plot_top10_nuclides(dchain_output,region_indices=2)
plt.show()
図 4.12.8 各時間インデックスにおける領域3の核種の放射能順位¶
DCHAINでは、複数の崩壊ライブラリと中性子反応断面積ライブラリから使用するデータを選択できる。20 MeV未満の中性子反応に用いる断面積ライブラリは、INXSLIBだけを変更して切り替えられる。ここでは、既定のハイブリッドライブラリとFENDL/A-3.0を使用した場合の、最終時刻における\(^{55}\)Fe放射能を比較する。既定のハイブリッドライブラリは、\(^{54}\)Fe(n,\(\gamma\))\(^{55}\)Fe反応と\(^{56}\)Fe(n,2n)\(^{55}\)Fe反応にJENDL/AD-2017のデータを使用する。
*.dcsファイルで\(^{55}\)Fe生成に寄与する系列を調べると、この2反応が主要な中性子反応であることを確認できる。本例では、両反応が\(^{55}\)Fe生成量の約3分の2を占める。残りの約3分の1は高エネルギー反応による生成であり、その大部分は20 MeVを超える\(^{56}\)Fe(n,2n)\(^{55}\)Fe反応と考えられる。このエネルギー領域はDCHAINの断面積ライブラリではなく、[T-Yield]タリーによって扱われる。以下のDCS出力例では、対応する値が青色で強調表示される。
... Ga 56 --( p)-> Zn 55 --(B+)-> Cu 55 --(B+)-> Ni 55 --(B+)-> Co 55 --(B+)-> Fe 55
... 0.00000E+00 1.49003E+06
... 1.35589E-14 -2.17501E+00
... 1.35589E-14 1.49003E+06
... Co 53m --(B+)-> Fe 53m --(IT)-> Fe 53 --(B+)-> Mn 53 --(nx)-> Mn 54 --(B-)-> Fe 54 --(nx)-> Fe 55
... -1.15353E-23 -1.04382E-23 -2.15030E-23 -1.80614E-16 0.00000E+00
... 8.14926E-24 1.03221E-23 2.20627E-23 5.29997E-22 1.95226E+06
... -3.38604E-24 -1.16058E-25 5.59759E-25 -1.80614E-16 1.95226E+06
... Ca 56 --(B-)-> Sc 56 --(B-)-> Ti 56 --(B-)-> V 56 --(B-)-> Cr 56 --(B-)-> Mn 56 --(B-)-> Fe 56 --(nx)-> Fe 55
... 0.00000E+00 0.00000E+00 -2.24854E-11
... 2.61382E-24 1.31490E-14 1.17913E+06
... 2.61382E-24 1.31490E-14 1.17913E+06
比較用に、PHITSが生成したspallation.outをspallation_FENDL.outという名前でコピーする。コピーしたファイル内のinxslib= 100をinxslib= 50へ変更し、FENDL/A-3.0を選択してDCHAINを実行する。生成された2つの*.actファイルはテキストエディターでも比較できるが、ここではDCHAIN Toolsを使用する。
既存のPythonスクリプトへ次のコードを追加し、既定ライブラリとFENDL/A-3.0の出力から\(^{55}\)Fe放射能を抽出する。生成される核種の種類はライブラリによって異なるため、領域3の核種一覧から\(^{55}\)Feのインデックスをそれぞれ取得する必要がある。
# Obtain activity of Fe-55 in final time step
nuclide_name = 'Fe-55'
ri, ti = 2, -1 # region 3 (region index 2), final time index
reg_num = dchain_output.region.number[ri]
time_str = seconds_to_dhms(dchain_output.time.from_start_sec[ti])
ni = dchain_output['nuclides']['names'][ri].index(nuclide_plain_str_to_Dname(nuclide_name)) # nuclide index (in reg 3 of spallation.act)
A_default_lib = dchain_output.nuclides.activity.value[ri][ti,ni]
simulation_base_name_2 = 'spallation_FENDL'
dchain_output_2 = process_dchain_simulation_output(simulation_folder_path,simulation_base_name_2,process_DCS_file=False)
ni_2 = dchain_output_2['nuclides']['names'][ri].index(nuclide_plain_str_to_Dname(nuclide_name)) # nuclide index (in reg 3 of spallation_FENDL.act)
A_FENDL = dchain_output_2.nuclides.activity.value[ri][ti,ni_2]
print("Activity of {} in region {} at {} was {} Bq/cc using default libraries and {} Bq/cc using FENDL/A-3.0".format(nuclide_name,reg_num,time_str,A_default_lib,A_FENDL))
Could not find default .dtrk file spallation_FENDL.dtrk, using spallation.dtrk in same directory instead.
Could not find default .dyld file spallation_FENDL.dyld, using spallation.dyld in same directory instead.
Parsing DCHAIN activation file... (0.95 seconds elapsed)
Restructuring nuclide data table array... (0.98 seconds elapsed)
Activity of Fe-55 in region 3 at 40m 0.00s was 14.796 Bq/cc using default libraries and 15.032 Bq/cc using FENDL/A-3.0
鉄の放射化データは十分に評価されているため、両ライブラリによる予測放射能は互いに近く、他の多くの核種ではさらによく一致する。
端末出力には、spallation_FENDLと同じ基底名を持つPHITSの[T-Track]・[T-Yield]出力が見つからなかったことが示されている。この場合、DCHAIN Toolsは同じフォルダー内にある、対応する拡張子の別ファイルをフラックス・収率データとして自動選択する。多数のファイルが同じフォルダーにあり、自動選択されたファイルが正しくない場合は、関数へ正しいパスを明示的に渡すこと。この点が結果へ影響するのは、読み込んだフラックスや収率を後続の計算に使用する場合である。
DCHAIN Toolsには、以下のように、各ライブラリについてDCHAINを個別に実行せず、特定反応の1群断面積を計算して相対的な反応確率を比較する機能もある。詳細は DCHAIN Tools documentation を参照すること。
# Obtain single-group cross section for various libraries
lib_folder_path = r'C:\phits\dchain-sp\data\\'
lib_dchain_names = ['ENDF-B-8-0','JENDL-AD17','JENDL-4-0-','FENDL-A-30','JEFF-3-3--','TENDL-2017','BROND-3-1-','CENDL-3-1-']
targets = ['Fe54','Fe56'] # targets
n_flux = dchain_output.neutron.spectra.flux.value[2]
n_flux_abs_err = dchain_output.neutron.spectra.flux.error[2]
for rxni in range(len(targets)):
print('\nCross sections (in millibarns) for {}(n,x)Fe55'.format(targets[rxni]))
for libi in range(len(lib_dchain_names)):
libfile = lib_folder_path+lib_dchain_names[libi]+'_n_act_xs_lib'
xs, xs_err = calc_one_group_nrxn_xs_dchain(n_flux,n_flux_abs_err,libfile,targets[rxni],product='Fe-55')
print('{}\t{:g} +/- {:g}'.format(lib_dchain_names[libi],1000*xs,1000*xs_err))
Cross sections (in millibarns) for Fe54(n,x)Fe55
ENDF-B-8-0 348.229 +/- 16.6916
JENDL-AD17 347.876 +/- 16.6426
JENDL-4-0- 347.876 +/- 16.6426
FENDL-A-30 347.308 +/- 16.6399
JEFF-3-3-- 348.611 +/- 16.6863
TENDL-2017 348.297 +/- 16.681
BROND-3-1- 439.381 +/- 21.0636
CENDL-3-1- 337.008 +/- 16.0373
Cross sections (in millibarns) for Fe56(n,x)Fe55
ENDF-B-8-0 14.1537 +/- 3.03331
JENDL-AD17 13.3846 +/- 2.86822
JENDL-4-0- 13.3846 +/- 2.86822
FENDL-A-30 14.2596 +/- 3.07603
JEFF-3-3-- 14.2596 +/- 3.07603
TENDL-2017 14.567 +/- 3.10122
BROND-3-1- 13.3852 +/- 2.86822
CENDL-3-1- 14.4761 +/- 3.09249