1 / 48

P HI T S

P HI T S. Multi-Purpose P article and H eavy I on T ransport code S ystem. PHITS 講習会　基礎実習（ III ） : 計算条件の設定. 2014 年 4 月改訂. title. 1. はじめに. PHITS 計算は， [Parameters] セクションで定義される様々なパラメータ（物理モデルの選択など）によりコントロールされています。 全てのパラメータには初期設定値があり，基本的には変更する必要はありません。

tricia
Download Presentation

P HI T S

An Image/Link below is provided (as is) to download presentation Download Policy: Content on the Website is provided to you AS IS for your information and personal use and may not be sold / licensed / shared on other websites without getting consent from its author. Content is provided to you AS IS for your information and personal use only. Download presentation by click this link. While downloading, if for some reason you are not able to download a presentation, the publisher may have deleted the file from their server. During download, if you can't get a presentation, the file might be deleted by the publisher.

E N D

Presentation Transcript

  1. PHITS Multi-Purpose Particle and Heavy Ion Transport code System PHITS 講習会　基礎実習（III）: 計算条件の設定 2014年4月改訂 title 1

  2. はじめに • PHITS計算は，[Parameters]セクションで定義される様々なパラメータ（物理モデルの選択など）によりコントロールされています。 • 全てのパラメータには初期設定値があり，基本的には変更する必要はありません。 • ただし，最適な計算結果を得るためには，いくつかのパラメータを変更する必要があります。 本実習では，パラメータの設定方法について学習します 2

  3. 本実習の目標 PHITSの便利な機能を活用し，統計数や核反応モデルを変化させ，最適な計算結果が得られるようになる 宿題体系における初期設定での 陽子（上）・中性子（下）フルエンス ヒストリー数，切断エネルギー，核データを適切 に設定した陽子（上）・中性子（下）フルエンス 3

  4. 実習内容 • 計算モードの選択 • 便利な入力補助機能 • 統計的に信頼できる結果を得るための設定 • モンテカルロ積分を使った体積・面積計算 • 試行回数と統計誤差 • 物理的に信頼できる結果を得るための設定 • 切断エネルギーの設定 • 断面積データライブラリの利用 • 物理モデルの選択 • まとめ 4

  5. 計算モードの選択 通常輸送計算 反応なしモード 体系確認モード 5

  6. 計算モードの選択 体系確認モード まずは、基礎実習(II)で行ったジオメトリを確認するタリーを用いて、サンプルインプットの体系を確認してみましょう。 6

  7. 体系の確認([t-3dshow]) [ T - 3 D s h o w ] output = 3 material = -1 6 x0 = 0 y0 = 0 z0 = 0 e-the = 70 $ eye e-phi = 20 e-dst = 80 l-the = 20 $ light l-phi = 0 l-dst = 100 w-wdt = 50 $ window w-hgt = 50 w-dst = 25 heaven = z line = 1 shadow = 2 file = 3dshow.out title = Check onion structure using [T-3dshow] tally epsout = 1 ・ ・ ・ ・ ・ ・ lec03.inp file = lec03.inp [ T i t l e ] ・ ・ ・ ・ ・ ・ [ P a r a m e t e r s ] icntl = 11 maxcas = 100 maxbch = 10 file(6) = phits.out [t-3dshow]を実行 1層5cm幅の玉ねぎ構造 3dshow.eps 7

  8. 体系の確認([t-gshow]) lec03.inp file = lec03.inp [ T i t l e ] ・ ・ ・ ・ ・ ・ [ P a r a m e t e r s ] icntl = 11 maxcas = 100 maxbch = 10 file(6) = phits.out file = lec03.inp [ T i t l e ] ・ ・ ・ ・ ・ ・ [ P a r a m e t e r s ] icntl = 7 maxcas = 100 maxbch = 10 file(6) = phits.out [ T - G s h o w ] mesh = xyz x-type = 2 nx = 100 xmin = -50. xmax = 50. y-type = 1 ny = 10 -25. -20. -15. -10. -5. 0. 5. 10. 15. 20. 25. z-type = 2 nz = 100 zmin = -50. zmax = 50. axis = xz output = 6 file = gshow.out title = Check onion structure using [T-gshow] tally epsout = 1 各領域のセル番号と満たされている物質名をチェック！ 修正して実行 gshow.eps 8

  9. 線源の確認([t-track]) lec03.inp 修正して実行 [ T - T r a c k ] mesh = xyz x-type = 2 nx = 100 xmin = -50. xmax = 50. y-type = 1 ny = 1 -5.0 5.0 z-type = 2 nz = 100 zmin = -50. zmax = 50. e-type = 1 ne = 1 0. 200. unit = 1 axis = xz file = track_xz.out title = Check source direction using [T-track] tally epsout = 1 file = lec03.inp [ T i t l e ] ・ ・ ・ ・ ・ ・ [ P a r a m e t e r s ] icntl = 5 maxcas = 100 maxbch = 10 file(6) = phits.out 9

  10. 線源の飛跡確認 Track-xz.eps 全ての物質がvoidとなるため，まっすぐに飛んでいく （線源の発生位置・方向が確認できる） 10

  11. 実習内容 • 計算モードの選択 • 便利な入力補助機能 • 統計的に信頼できる結果を得るための設定 • モンテカルロ積分を使った体積・面積計算 • 試行回数と統計誤差 • 物理的に信頼できる結果を得るための設定 • 切断エネルギーの設定 • 断面積データライブラリの利用 • 物理モデルの選択 • まとめ 11

  12. 外部ファイルの挿入 lec03.inp onion.inp file = lec03.inp [ T i t l e ] ・ ・ ・ ・ ・ ・ [ P a r a m e t e r s ] icntl = 5 maxcas = 100 maxbch = 10 file(6) = phits.out set: c1[10] [ S o u r c e ] ・ ・ ・ ・ ・ ・ infl: {onion.inp}[1-33] [ M a t e r i a l ] ・ ・ ・ ・ ・ ・ [ M a t N a m e C o l o r ] ・ ・ ・ ・ ・ ・ [ S u r f a c e ] 10 so 500. 11 so 5. 12 so 10. 13 so 15. 14 so 20. 15 so 25. [ C e l l ] 100 -1 10 101 1 -19.32 -11 102 2 -1. 11 -12 103 3 -8.93 12 -13 104 4 -1. 13 -14 105 5 -0.9 14 -15 106 6 -1.20e-3 15 -10 • インプットファイルとは別のファイルを取り込む場合は”infl:”コマンドを使う {}でファイル名，[n1-n2]で取り込む行の範囲を指定する • インプットファイルの一行目に“file = インプットファイル名”を書く(重要) • “off”したセクションの下にある場合、無視されるので注意 12

  13. 変数の利用 同じパラメータは，変数にしておくと便利！ lec03.inp file = lec03.inp [ T i t l e ] ・ ・ ・ ・ ・ ・ [ P a r a m e t e r s ] icntl = 5 maxcas = 100 maxbch = 10 file(6) = phits.out set: c1[10] [ S o u r c e ] s-type = 9 proj = proton x0 = 0. y0 = 0. z0 = 0. r1 = c1 r2 = c1 dir = 1.0 e0 = 150 • s-type=9, 10:　球および球殻分布の線源 • x0, y0, z0: 球の中心の座標 • r1: 内半径 • r2: 外半径（球面線源の場合，常にr1=r2） Z PHITSでは変数を定義して、インプットファイル内でその変数を参照することが可能。 “set: ci[x]” ＊書いた箇所より下の部分で有効となる ＊“off”したセクションの下にある場合、無視されるので注意 r1 Y r2 （x0, y0, z0) X 13

  14. 数式の利用 PHITS入力ファイルでは，FORTRAN形式の数式が利用可能 規格化について lec03.inp • PHITSのタリー結果は，発生する線源数*（試行回数）当たりで出力される • 測定値などと比較するには，Bqや/cm2/sで表される線源の発生頻度で規格化する必要がある • タリー結果を定数倍（totfactで定義）することにより，直接， /cm2/sやmGy/hなどの単位で出力できるようになる • この例の場合，規格化しないと線源球面のフルエンスは1/πr2 (/cm2/source)となるため，タリー結果にπr2を乗じることにより，球面におけるフルエンスが1(/cm2)となるように規格化している set: c1[10] [ S o u r c e ] s-type = 9 proj = proton x0 = 0. y0 = 0. z0 = 0. r1 = c1 r2 = c1 dir = 1.0 e0 = 150 totfact = pi*c1**2 • 円周率（pi）は定義しなくても使える • べき乗は「**」なので注意 　（C言語では「^」） • expやcosなどの関数も使える *厳密には，線源粒子のウェイト当たり 14

  15. 実習内容 • 計算モードの選択 • 便利な入力補助機能 • 統計的に信頼できる結果を得るための設定 • モンテカルロ積分を使った体積・面積計算 • 試行回数と統計誤差 • 物理的に信頼できる結果を得るための設定 • 切断エネルギーの設定 • 断面積データライブラリの利用 • 物理モデルの選択 • まとめ 15

  16. 体積、面積計算の原理 モンテカルロ積分: 乱数を用いて定積分の近似解を求める方法。積分範囲の内か外かを判定することにより、積分値を確率的に評価する。計算精度は試行回数に依存する。 体積計算にこの考え方を適用すると… 線の長さの合計 [cm] ＝　線の数(試行回数) [source] ×　平均track length [cm/source] 求める体積 [cm3] ＝　線の長さの合計÷線の面密度 ＝　平均track length÷平均flux ＝　[t-track]の結果　×線源面積 線の面密度(間隔) [1/cm2] ＝　線の数(試行回数) [source] ×　平均flux [1/cm2/source] あらかじめ 予測可能 試行回数を増やすことで精度が向上する 16

  17. 体積、面積計算の原理 モンテカルロ積分: 乱数を用いて定積分の近似解を求める方法。積分範囲の内か外かを判定することにより、積分値を確率的に評価する。計算精度は試行回数に依存する。 体積計算にこの考え方を適用すると… ただし、効率的に統計計算を行うために、線源の形状を球殻に取る。 この場合、線源面積は球を垂直な平面に射影したものとなる，すなわち πr2 試行回数を増やすことで精度が向上する 17

  18. 課題1: 体積、面積計算の実行 lec03.inp [ P a r a m e t e r s ] icntl = 5 maxcas = 100 maxbch = 10 file(6) = phits.out set: c1[30] [ S o u r c e ] s-type = 9 proj = proton x0 = 0. y0 = 0. z0 = 0. r1 = c1 r2 = c1 dir = -all e0 = 150 totfact = pi*c1**2 [ T - T r a c k ] mesh = reg reg = 101 102 103 104 105 ・ ・ ・ ・ ・ ・ ・ ・ ・ ・ ・ ・ file = volume.out [ T - C r o s s ] mesh = reg reg = 5 r-in r-out area (101 102) (101 102) 1.0 (102 103) (102 103) 1.0 (103 104) (103 104) 1.0 (104 105) (104 105) 1.0 (105 106) (105 106) 1.0 ・ ・ ・ ・ ・ ・ ・ ・ ・ ・ ・ ・ file = area.out 修正して実行 体系の体積や面積を計算する場合: • icntl=5, • s-type=9, • r1=r2（体系が全部入る大きさ） • dir=-all 18

  19. 体積、面積計算の結果 モンテカルロ積分の要領で体系の体積や面積を求めることが可能 volume.out area.out ・ ・ ・ ・ ・ ・ #num reg volume flux r.err 1 101 1.0000E+00 4.1594E+02 0.1881 2 102 1.0000E+00 2.9858E+03 0.0823 3 103 1.0000E+00 1.0534E+04 0.0451 4 104 1.0000E+00 2.0323E+04 0.0295 5 105 1.0000E+00 3.2595E+04 0.0211 ・ ・ ・ ・ ・ ・ #num area flux r.err 1 1.0000E+00 1.5700E+02 0.1901 2 1.0000E+00 1.1874E+03 0.1019 3 1.0000E+00 3.1593E+03 0.0948 4 1.0000E+00 5.5458E+03 0.0723 5 1.0000E+00 7.7191E+03 0.0469 各層の体積は？ • 4π(5)3/3=524 cm3 • 4π(10)3/3 - 524=3665 cm3 • 4π(15)3/3 - 4189=9948 cm3 • 4π(20)3/3 - 14137=19373 cm3 • 4π(25)3/3 - 33510=31940 cm3 各層の表面積は？ • 4π(5)2=314 cm2 • 4π(10)2=1257 cm2 • 4π(15)2=2827 cm2 • 4π(20)2=5027 cm2 • 4π(25)2=7854 cm2 内側に行くほど一致していない →　統計が足りていない 19

  20. 実習内容 • 計算モードの選択 • 便利な入力補助機能 • 統計的に信頼できる結果を得るための設定 • モンテカルロ積分を使った体積・面積計算 • 試行回数と統計誤差 • 物理的に信頼できる結果を得るための設定 • 切断エネルギーの設定 • 断面積データライブラリの利用 • 物理モデルの選択 • まとめ 20

  21. 試行回数の設定 モンテカルロ計算の統計精度は、シミュレーションの試行回数に依存しています。したがって、信頼できる計算結果を得るためには、十分な数の試行回数を設定する必要があります。 • 初期設定では，常に同じ乱数を使う • 毎回違う結果を得たい場合は，rseed < 0とする maxcas ×maxbch ＝全試行回数 21

  22. 試行回数と統計誤差の関係 モンテカルロ積分の要領で体系の体積や面積を求めることが可能 lec03.inp volume.out file = lec03.inp [ T i t l e ] ・ ・ ・ ・ ・ ・ [ P a r a m e t e r s ] icntl = 5 maxcas = 10000 maxbch = 10 file(6) = phits.out set: c1[30] [ S o u r c e ] s-type = 9 proj = proton x0 = 0. y0 = 0. z0 = 0. r1 = c1 r2 = 30 dir = -all e0 = 150 file = lec03.inp [ T i t l e ] ・ ・ ・ ・ ・ ・ [ P a r a m e t e r s ] icntl = 5 maxcas = 100 maxbch = 10 file(6) = phits.out set: c1[30] [ S o u r c e ] s-type = 9 proj = proton x0 = 0. y0 = 0. z0 = 0. r1 = c1 r2 = 30 dir = -all e0 = 150 file = lec03.inp [ T i t l e ] ・ ・ ・ ・ ・ ・ [ P a r a m e t e r s ] icntl = 5 maxcas = 1000 maxbch = 10 file(6) = phits.out set: c1[30] [ S o u r c e ] s-type = 9 proj = proton x0 = 0. y0 = 0. z0 = 0. r1 = c1 r2 = 30 dir = -all e0 = 150 ・ ・ ・ ・ ・ ・ #num reg volume flux r.err 1 101 1.0000E+00 4.1594E+02 0.1881 2 102 1.0000E+00 2.9858E+03 0.0823 3 103 1.0000E+00 1.0534E+04 0.0451 4 104 1.0000E+00 2.0323E+04 0.0295 5 105 1.0000E+00 3.2595E+04 0.0211 ・ ・ ・ ・ ・ ・ #num reg volume flux r.err 1 101 1.0000E+00 5.3772E+02 0.0508 2 102 1.0000E+00 3.5871E+03 0.0242 3 103 1.0000E+00 1.0055E+04 0.0146 4 104 1.0000E+00 1.9882E+04 0.0095 5 105 1.0000E+00 3.2469E+04 0.0067 ・ ・ ・ ・ ・ ・ #num reg volume flux r.err 1 101 1.0000E+00 5.2132E+02 0.0162 2 102 1.0000E+00 3.6695E+03 0.0076 3 103 1.0000E+00 9.8648E+03 0.0046 4 104 1.0000E+00 1.9316E+04 0.0031 5 105 1.0000E+00 3.1970E+04 0.0021 各層の体積は？ • 4π(5)3/3=524 cm3 • 4π(10)3/3 - 524=3665 cm3 • 4π(15)3/3 - 4189=9948 cm3 • 4π(20)3/3 - 14137=19373 cm3 • 4π(25)3/3 - 33510=31940 cm3 area.out 各層の表面積は？ • 4π(5)2=314 cm2 • 4π(10)2=1257 cm2 • 4π(15)2=2827 cm2 • 4π(20)2=5027 cm2 • 4π(25)2=7854 cm2 ・ ・ ・ ・ ・ ・ #num area flux r.err 1 1.0000E+00 3.0015E+02 0.0685 2 1.0000E+00 1.1869E+03 0.0349 3 1.0000E+00 2.7765E+03 0.0237 4 1.0000E+00 5.3456E+03 0.0218 5 1.0000E+00 7.8880E+03 0.0179 ・ ・ ・ ・ ・ ・ #num area flux r.err 1 1.0000E+00 3.2130E+02 0.0255 2 1.0000E+00 1.2291E+03 0.0122 3 1.0000E+00 2.8454E+03 0.0086 4 1.0000E+00 5.0156E+03 0.0063 5 1.0000E+00 7.8085E+03 0.0056 ・ ・ ・ ・ ・ ・ #num area flux r.err 1 1.0000E+00 1.5700E+02 0.1901 2 1.0000E+00 1.1874E+03 0.1019 3 1.0000E+00 3.1593E+03 0.0948 4 1.0000E+00 5.5458E+03 0.0723 5 1.0000E+00 7.7191E+03 0.0469 22

  23. 再開始計算 統計が足りない過去のタリー結果を読み込んで計算の再開が可能 lec03.inp volume.out file = lec03.inp [ T i t l e ] ・ ・ ・ ・ ・ ・ [ P a r a m e t e r s ] icntl = 5 maxcas = 10000 maxbch = 10 file(6) = phits.out set: c1[30] [ S o u r c e ] s-type = 9 proj = proton x0 = 0. y0 = 0. z0 = 0. r1 = c1 r2 = 30 dir = -all e0 = 150 ・ ・ ・ ・ ・ ・ #num reg volume flux r.err 1 101 1.0000E+00 5.2132E+02 0.0162 2 102 1.0000E+00 3.6695E+03 0.0076 3 103 1.0000E+00 9.8648E+03 0.0046 4 104 1.0000E+00 1.9316E+04 0.0031 5 105 1.0000E+00 3.1970E+04 0.0021 ・ ・ ・ ・ ・ ・ #num reg volume flux r.err 1 101 1.0000E+00 5.1768E+02 0.0115 2 102 1.0000E+00 3.6609E+03 0.0054 3 103 1.0000E+00 9.8629E+03 0.0033 4 104 1.0000E+00 1.9284E+04 0.0022 5 105 1.0000E+00 3.1893E+04 0.0015 ・ ・ ・ ・ ・ ・ #num reg volume flux r.err 1 101 1.0000E+00 5.2004E+02 0.0094 2 102 1.0000E+00 3.6612E+03 0.0044 3 103 1.0000E+00 9.8995E+03 0.0027 4 104 1.0000E+00 1.9299E+04 0.0018 5 105 1.0000E+00 3.1800E+04 0.0012 istdev = -1 area.out ・ ・ ・ ・ ・ ・ #num area flux r.err 1 1.0000E+00 3.2130E+02 0.0255 2 1.0000E+00 1.2291E+03 0.0122 3 1.0000E+00 2.8454E+03 0.0086 4 1.0000E+00 5.0156E+03 0.0063 5 1.0000E+00 7.8085E+03 0.0056 ・ ・ ・ ・ ・ ・ #num area flux r.err 1 1.0000E+00 3.1065E+02 0.0180 2 1.0000E+00 1.2476E+03 0.0090 3 1.0000E+00 2.8271E+03 0.0060 4 1.0000E+00 5.0014E+03 0.0045 5 1.0000E+00 7.8189E+03 0.0040 ・ ・ ・ ・ ・ ・ #num area flux r.err 1 1.0000E+00 3.1013E+02 0.0148 2 1.0000E+00 1.2660E+03 0.0075 3 1.0000E+00 2.8413E+03 0.0050 4 1.0000E+00 5.0008E+03 0.0037 5 1.0000E+00 7.7990E+03 0.0033 23

  24. バッチ計算の活用 PHITSでは、タリーの結果を各バッチが終了する度に出力させることができ、それを見ながらモンテカルロ計算を中断することができます。 batch.now 計算の途中で1を0に書き換え保存すると その次のバッチで計算をやめる。 1 <--- 1:continue, 0:stop ------------------------------------------------------------------------------- bat[ 560] ncas = 560. : rijk = 151264979546685. low neutron = 0. : ncall/s = 4.000000000E+00 cpu time = 0.288 s. date = 2012-05-02 time = 15h 08m 25 24

  25. バッチ計算の活用 lec03.inp file = lec03.inp [ T i t l e ] ・ ・ ・ ・ ・ ・ [ P a r a m e t e r s ] icntl = 5 itall = 1 maxcas = 10000 maxbch = 10 file(6) = phits.out $ istdev = -1 • 計算が終わらない場合は、 batch.now source-direction.eps 1 <--- 1:continue, 0:stop ------------------------------------------------------------------------------- 0 <--- 1:continue, 0:stop ------------------------------------------------------------------------------- 途中のバッチで終わっているか phits.outで確認してみましょう。 25

  26. 実習内容 • 計算モードの選択 • 便利な入力補助機能 • 統計的に信頼できる結果を得るための設定 • モンテカルロ積分を使った体積・面積計算 • 試行回数と統計誤差 • 物理的に信頼できる結果を得るための設定 • 切断エネルギーの設定 • 断面積データライブラリの利用 • 物理モデルの選択 • まとめ 26

  27. 切断エネルギーとは？ • 注目していない、或いは重要でないと思われる粒子の輸送や核反応の発生をストップし、数値計算の量を減らすことが可能です。 • 初期設定では輸送しない粒子（低エネルギー中性子，光子，電子など）を輸送するように変更できます。 27

  28. 通常の輸送計算を実行 lec03.inp 輸送計算を実行する [ T - C r o s s ] mesh = reg reg = 5 r-in r-out area 102 101 1.0 103 102 1.0 104 103 1.0 105 104 1.0 106 105 1.0 e-type = 2 ne = 100 emin = 0. emax = 200. axis = eng unit = 1 output = flux file = cross.out epsout = 1 [ T - C r o s s ]off mesh = reg reg = 5 r-in r-out area 102 101 1.0 103 102 1.0 104 103 1.0 105 104 1.0 106 105 1.0 e-type = 2 ne = 100 emin = 0. emax = 200. axis = eng unit = 1 output = flux file = cross.out epsout = 1 修正して実行 [ P a r a m e t e r s ] icntl = 5 itall = 1 maxcas = 10000 maxbch = 10 file(6) = phits.out c istdev = -1 set: c1[30] [ S o u r c e ] s-type = 9 proj = proton ・ ・ ・ ・ ・ ・ ・ ・ ・ ・ ・ ・ dir = -all e0 = 150 [ P a r a m e t e r s ] icntl = 0 itall = 1 maxcas = 1000 maxbch = 10 file(6) = phits.out c istdev = -1 set: c1[30] [ S o u r c e ] s-type = 9 proj = proton ・ ・ ・ ・ ・ ・ ・ ・ ・ ・ ・ ・ dir = -all e0 = 150 玉ねぎ構造の各層を通過する粒子をタリーする。 28

  29. 切断エネルギーの設定 lec03.inp 50MeV以下の中性子を輸送させない。 [ P a r a m e t e r s ] icntl = 0 itall = 1 maxcas = 1000 maxbch = 10 emin(2) = 50 file(6) = phits.out # istdev = -1 set: c1[30] 修正して実行 cross.eps 50MeV未満の中性子はカット → 輸送されず各層を通過していない 29

  30. 実習内容 • 計算モードの選択 • 便利な入力補助機能 • 統計的に信頼できる結果を得るための設定 • モンテカルロ積分を使った体積・面積計算 • 試行回数と統計誤差 • 物理的に信頼できる結果を得るための設定 • 切断エネルギーの設定 • 断面積データライブラリの利用 • 物理モデルの選択 • まとめ 30

  31. 断面積データライブラリとは？ • 光子・電子・陽電子並びに低エネルギー中性子（20MeV以下）の断面積は，ターゲットやエネルギーに複雑に依存するため，画一的なモデルでは表現できない → 断面積データライブラリが必要 • 断面積ライブラリを使った計算は，計算時間が掛かる場合があるため，デフォルトでは，それらの粒子の切断エネルギーは高く設定されている 正しい輸送計算結果を得るためには，断面積データライブラリを使うように[parameters]セクションでeminなどを設定する必要がある 113Cdの中性子反応断面積（JENDL-4.0より） 31

  32. データライブラリの設定 • データライブラリファイルの確認 c:/phits/XS/neu • アドレス（xsdir）ファイルの確認 c:/phits/data/xsdir.jnd • アドレスファイルの１行目の確認 datapath=c:/phits/XS • [Parameters]セクションにおいて“emin(i)” “dmax(i)” ，“file(7)”パラメータを設定 中性子用核データライブラリ アドレスファイル データライブラリを含むフォルダ名 32

  33. 課題2: 中性子用核データライブラリの利用 lec03.inp phits.out [ P a r a m e t e r s ] icntl = 0 itall = 1 maxcas = 1000 maxbch = 10 emin(2)= 1.0e-10 dmax(2) = 20 file(6) = phits.out file(7) = c:/phits/data/xsdir.jnd # istdev = -1 --------------------------------------------------------------- CPU Summary --------------------------------------------------------------- ・ ・ ・ ・ ・ ・ ・ ・ ・ ・ ・ ・ ・ ・ ・ ・ ・ ・ other = 0.00 0.00 hydro = 0.00 0.00 310. dklos = 0.00 0.00 0. elast = 0.00 0.00 1348. ncasc = 0.00 0.00 4701. n-data = 0.00 0.00 38974. p-data = 0.00 0.00 0. e-data = 0.00 0.00 0. h-data = 0.00 0.00 0. berti = 0.00 0.00 0. ・ ・ ・ ・ ・ ・ ・ ・ ・ ・ ・ ・ 中性子のデータを使う (emin < dmaxを確認) アドレスファイル名 33

  34. 課題3: 光子・電子用データライブラリの利用 • 電子・陽電子・光子の切断エネルギーemin(12~14) を設定する 　　 　→計算時間が長くなるので，ここでは1MeVに設定 • 電子・陽電子・光子のライブラリ上限エネルギーdmax(12~14)を設定する 　→　全て1GeVに設定 PHITSを実行 • phits.outにあるe-data，p-dataを確認 演習 [t-track](track_xz.out)のpartを変えて電子・光子・陽電子の飛跡を確認してみよう　 34

  35. 実習内容 • 計算モードの選択 • 便利な入力補助機能 • 統計的に信頼できる結果を得るための設定 • モンテカルロ積分を使った体積・面積計算 • 試行回数と統計誤差 • 物理的に信頼できる結果を得るための設定 • 切断エネルギーの設定 • 断面積データライブラリの利用 • 物理モデルの選択 • まとめ 35

  36. γ脱励起オプション igamma: 核反応で生成される放射性核種からのγ線放出を考慮するためのオプション。デフォルトはγ線放出を考慮しないので注意が必要 Version 2.64以降の奨励値は igamma = 2 36

  37. 荷電粒子のビームライン設計オプション nspred & nedisp: 荷電粒子の物質中での角度・エネルギー分散を考慮するためのオプション。加速器ビームのビームライン設計の際，不可欠となる 奨励値は nspred = 1 & nedisp = 1 37

  38. 核反応モデルの変更 PHITSには、Bertini, JAM, JQMD, INCL, INC-ELFという核反応モデルがあり、状況に応じてユーザーが使い分けることができます。 38

  39. 核反応モデルの切替エネルギー (=emin) dmax(i) (3.0GeV) einclmax (1MeV) emin(i) Nucleon 核データ INCL (inclg=1) JAM (1MeV) emin(i) (3.0GeV) einclmax Pion INCL (inclg=1) JAM (3.5GeV/u) ejamqmd (10MeV/u) eqmdmin Nucleus JQMD JAMQMD (d, t, 3He, α) INCL (inclg=1) Kaon, Hyperon JAM 39

  40. イベントジェネレーターモード 核データ＋統計崩壊モデルの利用により、残留核からの放出粒子まで考慮し、かつエネルギーと運動量の保存則を満たした核反応イベントを生成するモード。 イベントジェネレータモードを使った方がよい場合 • 低エネルギー中性子が放出する陽子やα粒子スペクトルが見たい • 残留核の情報（反跳エネルギーなど）が知りたい • イベント毎の情報が知りたい（検出器応答関数の計算など） イベントジェネレータモードを使わない方がよい場合 • 中性子と光子の情報だけ知りたい（遮へい計算など） • 中性子の透過率を計算したい • 熱中性子の正確な挙動が見たい 40

  41. イベントジェネレーターモードの設定 • 使用する方法は、核データを使用する設定を施した上で、[parameters]セクションにおいて”e-mode=1”とするだけです。(指定しない場合は”e-mode=0”となっています。) • ガンマ崩壊データ（trxcrd.dat）の設定が必要です(“igamma”パラメータは自動的に2になります） 41

  42. 課題4: イベントジェネレーターモードの利用 lec03.inp 修正して実行 [ T - T r a c k ] mesh = reg reg = (101 102 103 104 105) e-type = 2 ne = 100 emin = 0 emax = 20. axis = eng unit = 1 part = proton neutron photon alpha file = track_eng.out epsout = 1 [ P a r a m e t e r s ] icntl = 0 itall = 1 maxcas = 1000 maxbch = 10 emin(2)= 1.0e-10 dmax(2) = 20 emin(12) = 1.0 emin(13) = 1.0 emin(14) = 1.0 dmax(12) = 1000.00000 dmax(13) = 1000.00000 dmax(14) = 1000.00000 file(6) = phits.out file(7) = c:/phits/data/xsdir.jnd # istdev = -1 e-mode = 0 [ S o u r c e ] s-type = 9 proj = neutron x0 = 0. y0 = 0. z0 = 0. r1 = c1 r2 = c1 dir = -all e0 = 20.0 totfact = pi*c1**2 • 各領域の合計値を出力したい場合は，合計したい領域を（）で囲む • まずはイベントジェネレーターモードを使わずに計算してみる。（デフォルト設定） • 中性子の影響を見るため，線源は20MeV中性子とする • 球全体の粒子フルエンスを出力して，陽子・α粒子が検出されたか確認する 42

  43. 課題4: イベントジェネレーターモードの利用 lec03.inp 修正して実行 [ P a r a m e t e r s ] icntl = 0 itall = 1 maxcas = 1000 maxbch = 10 emin(2)= 1.0e-10 dmax(2) = 20 emin(12) = 1.0 emin(13) = 1.0 emin(14) = 1.0 dmax(12) = 1000.00000 dmax(13) = 1000.00000 dmax(14) = 1000.00000 file(6) = phits.out file(7) = c:/phits/data/xsdir.jnd # istdev = -1 e-mode = 1 file(14) = c:/phits/data/trxcrd.dat Track_eng.eps 陽子やα粒子のスペクトルが見えた！ 43

  44. 入出力ファイルの設定 44

  45. 実習内容 • 計算モードの選択 • 便利な入力補助機能 • 統計的に信頼できる結果を得るための設定 • モンテカルロ積分を使った体積・面積計算 • 試行回数と統計誤差 • 物理的に信頼できる結果を得るための設定 • 切断エネルギーの設定 • 断面積データライブラリの利用 • 物理モデルの選択 • まとめ 45

  46. まとめ PHITSでは、[Parameters]セクションの設定を変更することで、様々な計算モードや物理モデルの選択を行うことができる。 各計算モード（icntl）を使い分けることにより、ジオメトリの確認，ソースのチェック、粒子の輸送計算を行う。 統計精度と試行回数(maxcas, maxbch)は密接に関連しており、適切な数を設定する必要がある。 低エネルギーの中性子や光子，電子，陽電子等は，切断エネルギー（emin）や核データの上限エネルギー（dmax）を変更することで精度の高い計算結果を得ることができる。その際，file(7)でアドレスファイルを指定する必要がある 状況に応じてイベントジェネレーターモード(e-mode)など物理モデルの設定を変更する必要がある。 最適な設定が分からない場合は，奨励設定ファイル(recommendation)をご参照ください 46

  47. 宿題 宿題ファイルで，熱中性子まで考慮し，20MeV以下の中性子に対して核データライブラリを使う イベントジェネレータモードを使ってみる 線量-深さ分布の相対標準偏差２％以内を目指す（maxcas, istdev, batch.nowなどを活用する） 47

  48. 宿題（解答例） 陽子（上）・中性子（下）フルエンス 円柱中心（上）と端側（下）における 付与エネルギーの深さ分布 48

More Related