ISSN 2434-1088

KURNS-EKR- 2

第7回「炉物理専門研究会」
Proceedings of 7th Reactor Physics Workshop(RPW2018)

編集： 卞 哲浩
Edited by : Cheol Ho Pyeon

京都大学複合原子力科学研究所
Institute for Integrated Radiation and Nuclear Science, Kyoto University要旨

本研究会は、京都大学臨界集合体実験装置（KUCA）で行われた共同利用研究者による実験および解析結果を内外に広く公表し、その成果を多くの研究者たちと議論することによって、KUCA の共同利用実験の発展に資することを目的としている。京都大学複合原子力科学研究所では、KUCA（A 架台）と FFAG 加速器を組み合わせて加速器駆動システム（以下 ADS: Accelerator-Driven System）
を構成し、核変換技術への適用性に関する基礎研究を行っている。特に、KUCA で行われている
Pb-Bi に関連する ADS 実験は、日本原子力研究開発機構（JAEA）の J-PARC 施設のひとつとして建設が検討されている核変換実験施設 TEF（Transmutation Experimental Facility）を用いた ADS 研
究に対して、炉物理研究および核データ研究の基盤基礎強化に大きく貢献することが期待されている。これらの実験結果が外部の研究者たちによって客観的に評価され、意見交換を積極的に行うことによって、ADS 研究のさらなる発展が研究会を通して行われている。
原子炉物理実験の解析を精度良く行うためには、計算科学および核データ分野との連携は極めて重要である。核計算および核データ関連の研究者たちによる広範な視点から、これまで得られた研究成果を活発に議論し、ADS 研究における計算科学および核データ分野の研究課題を互いに共有することが本研究会において可能になっている。原子炉物理の研究成果を国内で議論する機会が原子力学会および炉物理夏期セミナーなどに限られていることから、参加者の研究成果が第三者により評価される機会として、また、原子炉物理研究をさらに発展させる場としてこの研究会が大いに活用され、さらに、産学官の研究機関の若手研究者および学生たちのスキルアップの機会となれば幸いである。
最後に、本研究会の開催に向けてご尽力いただいた名古屋大学・山本章夫教授、北海道大学・
千葉豪准教授、名古屋大学・遠藤知弘助教、東北大学・相澤直人助教、大阪大学・竹田敏助教、JAEA・
山根祐一氏および京都大学複合原子力科学研究所・山中正朗氏に心より感謝申し上げます。
卞 哲浩
2018 年 12 月iPreface

The objective of this workshop is to open all the results of experiments carried out at the
Kyoto University Critical Assembly (KUCA) and develop all future activities of joint use at
KUCA through the discussion about the experimental topics together with all researchers andengineers. In the Institute for Integrated Radiation and Nuclear Science, Kyoto University(KURNS), the accelerator-driven system (ADS) is composed of the KUCA core and thefixed-field alternating gradient (FFAG) accelerator, and the research and development of
ADS are being conducted to examine the feasibility of the application of ADS to the nuclear
transmutation techniques.
It is very important to share the experimental field with the mathematical andcomputational (M&C), and nuclear data (ND) fields in terms of the analyses of reactor
physics experiments. From this context, another purpose of this workshop is to share theresults of experimental data with the researchers in the M&C and nuclear data fields throughthe discussion with them.
Further, it is expected that this workshop could be contributed to the human resource
training for young researchers and students in domestic, through their researchpresentations.
Finally, we would like to give special thanks for their support and patience, by Prof. Akio
Yamamoto of Nagoya University, Prof. Go Chiba of Hokkaido University, Prof. Tomohiro
Endo of Nagoya University, Prof. Naoto Aizawa of Tohoku University, Prof. Satoshi Takeda,
Dr. Yuichi Yamane of JAEA, and Dr. Masao Yamanaka of KURNS, to hold this workshop.

Cheol Ho Pyeon
December 2018

Keywords:
Reactor physics, KUCA, M&C, Nuclear data, ADS

ii

目
1.

次
人生 100 年時代の炉物理キャリアパス
山崎正俊（Studsvik Japan）

2.

Adaptive smooth-lasso を用いた感度係数評価に関する検討
方野量太（JAEA）････････････････････････････････････････････････････････････････

3.

1

水素の熱中性子散乱則の違いによる軽水炉体系での影響
竹田 敏（大阪大学）････････････････････････････････････････････････････････････ 13

4.

KUCA 固体減速架台の数値計算ベンチマーク問題
千葉 豪（北海道大学）･･････････････････････････････････････････････････････････

5.

22

Effect of Differences Fuel Materials on Neutronic Parameters in Kartini Research Reactor
Argo Satrio Wicaksono（大阪大学大学院）････････････････････････････････････････････

41

未臨界施設での炉物理実験の検討
中込宇宙（東芝 ESS）

6.

KUCA における外部中性子源を用いた逆動特性解析の検討
山中正朗（京都大学複合研）

7.

様々な核種の崩壊や生成などを考慮した福島第一原子力発電所内部の放能インベントリ解析
石井亮憲（北海道大学大学院）･･････････････････････････････････････････････････

8.

49

臨界シナリオに基づくデブリ取り出し作業を想定した臨界挙動解析
木村 礼（東芝 ESS）

9.

Supercritical transient analysis using Multi-region Integral Kinetic code: Basics and applications
Delgersaikhan Tuya（東京工業大学）･･････････････････････････････････････････････

62

10. ブートストラップ法を活用した Feynman-α法の実験共分散評価
遠藤知弘（名古屋大学）････････････････････････････････････････････････････････

73

11. 臨界近接監視システムの開発(1) KUCA を用いた臨界近接監視システムの検証
加納慎也（東芝 ESS）

12. 臨界近接監視システムの開発(2) MVP 2.0.33 - Time-list モードの検証
和田怜志（東芝 ESS）

13. KUCA における鉛とビスマスの Al 置換反応度価値の測定
大泉昭人（JAEA）

14. KUCA 実験データを利用した粒子フィルタ法による未臨界度推定池田卓弥（名古屋大学大学院）････････････････････････････････････････････････････ 84

15. 核破砕中性子源によるウラン-鉛ゾーン炉心 ADS 体系における放射化反応率への中性子スペ
クトルの影響相澤直人（東北大学）････････････････････････････････････････････････････････････ 97iii

Contents
1.

Nuclear Physics Career Path: Working in the Age of 100-year Life
(Studsvik, Japan) Masatoshi Yamasaki

2.

A Study for the Estimation of the Sensitivity Coefficients using Adaptive Smooth-Lasso(JAEA) Ryota Katano ･･････････････････････････････ 1

3.

Study on Impact of Thermal Scattering Law Data Improvements on Nuclear Calculation forLight Water Reactor(Osaka Univ.) Satoshi Takeda ･･････････････････････ 13

4.

Benchmark Problems for Neutron Transport Calculations about KUCA Solid-Moderated
Cores(Hokkaido Univ.) Go Chiba ････････････････････････ 22

5.

Effect of Differences Fuel Materials on Neutronic Parameters in Kartini Research Reactor(Osaka Univ.) Argo Satrio Wicaksono ･･･････････････ 41

6.

Nuclear Facility for Subcritical Reactor Experiment(Toshiba ESS) Hiroshi Nakagomi

7.

Reactivity Measurement by Inverse Point Kinetics Method with External Neutron Source at
KUCA(Kyoto Univ.) Masao Yamanaka

8.

Radioactivity Inventory Analysis in the Fukushima Daiichi NPP Considering such as Various
Nuclides Generation and Decay(Hokkaido Univ.) Akinori Ishii ･････････････････････ 49

9.

Criticality Behavior Analyses based on Criticality Scenario during Fuel Debris Retrieval
Process(Toshiba ESS) Rei Kimura

10. Supercritical transient analysis using Multi-region Integral Kinetic Code: Basics andapplications(Tokyo Inst. Technol.) Delgersaikhan Tuya ･･･････････ 62
11. Estimation of Experimental Covariance in the Feynman-α Method using the Bootstrap
Method(Nagoya Univ.) Tomohiro Endo ････････････････････ 73
12. Development of Sub-Criticality Monitoring System for Fukushima Daiichi Fuel Debris(1) Verification of Sub-Criticality Monitoring System using KUCA
(Toshiba ESS) Shinya Kanoiv

13. Development of Sub-Criticality Monitoring System for Fukushima Daiichi Fuel Debris(2) Verification of Time-List Mode of MVP 2.0.33(Toshiba ESS) Satoshi Wada
14. Sample Worth Measurement with Systematically Changed Mixing Ratio of Lead and Bismuth
in KUCA(JAEA) Akito Oizumi
15. Estimation of Subcriticality and Neutron Generation Time using the Particle Filter Method(Nagoya Univ.) Takuya Ikeda ･･････････････････････ 84
16. Effect of Neutron Spectrum on Reaction Rate in ADS Experiment with U-Pb Zoned Core andSpallation Source(Tohoku Univ.) Naoto Aizawa ･････････････････････ 97v

第 7 回「炉物理専門研究会」のご案内

日時：2018 年 11 月 26 日（月）および 11 月 27 日（火）場所：京都大学 複合原子力科学研究所 事務棟大会議室
参加費：11/26 の懇親会費 1,500 円、11/27 の昼食代 500 円
2018 年 11 月 26 日（月）
12:30

受付

13:00 - 13:05

開会（敬称略、名大・山本章夫）

13:05 - 14:05

Session I：Special session（阪大・竹田 敏）

山崎正俊（Studsvik Japan）「人生 100 年時代の炉物理キャリアパス」
14:15 - 15:45

Session II：核データおよび不確かさ解析（北大・千葉 豪）

方野量太（JAEA）「Adaptive smooth-lasso を用いた感度係数評価に関する検討」
竹田 敏（阪大）「水素の熱中性子散乱則の違いによる軽水炉体系での影響」
千葉 豪（北大）「KUCA 固体減速架台の数値計算ベンチマーク問題」
16:00 - 17:30

Session III：炉物理一般（名大・遠藤知弘）

Argo Satrio Wicaksono（阪大）「Effect of Differences Fuel Materials on Neutronic Parameters in Kartini Research Reactor」
中込宇宙（東芝 ESS）
「未臨界施設での炉物理実験の検討」
山中正朗（京大複合研）「KUCA における外部中性子源を用いた逆動特性解析の検討」
17:45 - 20:00

懇親会

vi

2018 年 11 月 27 日（火）9:00 - 10:30

Session IV：臨界安全・燃料デブリ（JAEA・山根祐一）

石井亮憲（北大）「様々な核種の崩壊や生成などを考慮した福島第一原子力発電所内部の放能インベント
リ解析」
木村 礼（東芝 ESS）
「臨界シナリオに基づくデブリ取り出し作業を想定した臨界挙動解析」
Delgersaikhan Tuya（東工大）「 Supercritical transient analysis using Multi-region Integral Kinetic code: Basics andapplications」
10:45 - 12:15

Session V：KUCA 実験 I（東北大・相澤直人）

遠藤知弘（名大）「ブートストラップ法を活用した Feynman-α法の実験共分散評価」
加納慎也（IRID／東芝 ESS）
「臨界近接監視システムの開発(1) KUCA を用いた臨界近接監視システムの検証」
和田怜志（IRID／東芝 ESS）
「臨界近接監視システムの開発(2) MVP 2.0.33 - Time-list モードの検証」
13:15 - 14:45

Session VI：KUCA 実験 II（京大複合研・山中正朗）

大泉昭人（JAEA）「KUCA における鉛とビスマスの Al 置換反応度価値の測定」
池田卓弥（名大）「KUCA 実験データを利用した粒子フィルタ法による未臨界度推定」
相澤直人（東北大）「核破砕中性子源によるウラン-鉛ゾーン炉心 ADS 体系における放射化反応率への中性
子スペクトルの影響」
14:45

閉会（京大複合研・卞 哲浩）

vii

Adaptive smooth-lasso を用いた感度係数評価に関する検討
A study for the estimation of the sensitivity coefficients using adaptive smooth-lasso＊

方野 量太 1
1

JAEA

Forward 計算のみで核特性の感度係数を評価する手法として、筆者は罰則化線形回帰手法である adaptive
smooth-lasso を提案した[1]。本検討は、本手法の更なる適用性評価として、加速器駆動システム(ADS)のビーム電流値の感度係数を評価したものである。

キーワード：感度係数、ランダムサンプリング、罰則化線形回帰、L1 ノルム

1. 緒言

通常、核反応数に比べ微視的多群断面積数が大きい場合、感度係数は(一般化)摂動論を用いて評価される。しかし、複雑な解析システムに対する一般化随伴方程式ソルバーの実装は多大な労力を要し、また核特性数の数が増えるほど摂動論による計算コストは増大する。そこで筆者は、Forward 計算のみを用いる感度係数評価手法として、adaptive smooth-lasso を提案した[1]。本手法はランダムサンプリングで得られる断面積摂動と核特性との罰則化回帰分析を通じて感度係数を評価するもので、感度係数の特徴である
1)多くは非常に小さい(スパースである)、2) 炉型に応じてあるエネルギー領域にまとまるが共鳴や閾値反応によってそのエネルギー依存性は複雑である、ということを考慮して罰則項を設定するものであり、従来の罰則化回帰分析法である lasso や smooth-lasso を改良したものである。本検討は、出力 800MW の ADS
初期炉心のビーム電流値の感度係数評価を通じ、本手法の更なる適用性評価を行ったものである。
2. 計算条件

ADS 未燃焼炉心(keff = 0.97)について、固定源計算を実施し、炉心熱出力が 800MW となる陽子ビーム加速器の電流値の感度係数を評価した。考慮した微視的多群断面積数は 13286 個（73 群、26 核種、
7 反応）であり、ランダムサンプリングは断面積間の相関を考慮せず±5%の範囲で一様にサンプリングした。参照解は直接法で算出し、本手法、lasso、smooth-lasso を用いて感度係数を推定した。
3. 計算結果

図 1 は 750 サンプルによる Pu-239 の核分裂断面積に対する感度係数の推定結果である。参照

解と比較して、lasso や smooth-lasso はやや絶対値を過小評価する傾向にあるが、提案手法の結果は同じサンプル数で絶対値と形状をよく再現しており、ビーム
電流値の感度係数に対しても従来法と比べてより良い性能を持つことを示した。但し、提案手法は 4 つのユーザ定義の調整パラメータを有しており、パラメ
ータの自動最適化が今後の課題である。
参考文献[1] R. Katano, et al., J. Nucl. Sci. Technol., 55,
1099-1109, (2018).
*Ryota Katano1
1JAEA

図 1 Pu-239 核分裂断面積に対する感度係数

Adaptive smooth-lassoを用いた
感度係数評価に関する検討(JAEA) 方野 量太
第7回炉物理専門研究会
2018年11月26日

1

背景 感度係数 不確かさ評価、断面積調整法に利用 (一般化)摂動論を用いて評価

 核特性毎に Forward/Adjoint 方程式の2つを解けばよい
 実装は煩雑...

 燃焼感度→複雑 非線形フィードバック→???
 複数コードによる炉心解析への適用は簡単ではない ランダムサンプリング (RS) Forward計算のみで不確かさ評価

 計算コスト < 直接法 断面積調整法も提案 [1] 複雑な処理の不要、実行が簡単[1] T. Watanabe, et al. J. Nucl. Sci. Technol., 51, (2014).
2

背景 RSによる感度係数評価
 いくつか提案されている

 一般化逆行列法[2], L1ノルム最小化法[3] 新たに線形回帰分析法を提案 [4] Adaptive smooth-lasso と命名、後程説明 シンプルながら従来法より高精度
 未燃焼keffは検証済み、更なる検討の必要性

 本検討
 加速器ビーム電流値の感度係数を評価

 固定源＋拘束条件（定格出力）→ 提案手法の更なる検証[2] G. Chiba, et al. Ann. Nucl. Energy., 75, (2015).
[3] T. Watanabe, et al. Trans. Am. Nucl. Soc., 111, (2014).
[4] R. Katano, et al. J. Nucl. Sci. Technol., 55, (2018).

3

理論
 RSによる感度係数評価
 断面積摂動と核特性摂動の連立方程式

 核特性摂動 = 断面積摂動 ・ 感度係数

𝛥𝑅
𝛥𝑅
⋮
𝛥𝑅

𝛥𝑇
𝛥𝑇
⋮
𝛥𝑇

𝛥𝑇
𝛥𝑇
⋮
⋯

⋯ 𝛥𝑇
⋯
⋮
⋱
⋮
⋯ 𝛥𝑇

𝑔
𝑔
⋮
𝑔

4

理論
 RSによる感度係数評価
 断面積摂動と核特性摂動の連立方程式

 核特性摂動 = 断面積摂動 ・ 感度係数

𝛥𝑅
𝛥𝑅
⋮
𝛥𝑅

𝛥𝑇
𝛥𝑇
⋮
𝛥𝑇

𝛥𝑇
𝛥𝑇
⋮
⋯

⋯ 𝛥𝑇
⋯
⋮
⋱
⋮
⋯ 𝛥𝑇

𝑔
𝑔
⋮
𝑔

断面積数

5

理論
 RSによる感度係数評価
 断面積摂動と核特性摂動の連立方程式

サンプル数

 核特性摂動 = 断面積摂動 ・ 感度係数

𝛥𝑅
𝛥𝑅
⋮
𝛥𝑅

𝛥𝑇
𝛥𝑇
⋮
𝛥𝑇

𝛥𝑇
𝛥𝑇
⋮
⋯

⋯ 𝛥𝑇
⋯
⋮
⋱
⋮
⋯ 𝛥𝑇

𝑔
𝑔
⋮
𝑔

6

理論
 RSによる感度係数評価
 断面積摂動と核特性摂動の連立方程式

 核特性摂動 = 断面積摂動 ・ 感度係数

𝛥𝑅
𝛥𝑅
⋮
𝛥𝑅

𝛥𝑇
𝛥𝑇
⋮
𝛥𝑇

𝛥𝑇
𝛥𝑇
⋮
⋯

⋯ 𝛥𝑇
⋯
⋮
⋱
⋮
⋯ 𝛥𝑇

𝑔
𝑔
⋮
𝑔

 サンプル数 < 断面積数なら計算コスト削減
 劣決定系となり一意に決定できない
 拘束条件 or 罰則化回帰分析が必要

 L2 ノルム最小化 [2] L1 ノルム最小化 [3] ridge, lasso, smooth-lasso
7

理論
 解をどうやって決める？
 例）断面積が２つ、計算回数１回の場合
 解は直線上の点のどこか

解空間

8

理論
 解をどうやって決める？
 例）断面積が２つ、計算回数１回の場合
 解は直線上の点のどこか

 例：”一番小さい”ものを選ぶ

円との交点！
解空間

感度係数の特徴をとらえているか？
9理論
 感度係数の特徴
Relative sensitivity coefficient [‐] 例) 未燃焼ADS keffの感度 (下図：Np-237 cap.)
0.04
0.03
0.02

ほとんどゼロ、スパース

0.01
0.00
‐0.01
‐0.02
0

1022

2044

3066

4088

5110 6132 7154 8176
Cross section index

9198 10220 11242 12264 13286

Relative sensitivity coefficient[/lethergy]

0.0E+0‐5.0E‐3

高エネルギーにまとまる
滑らかでなはい、急峻な変化‐1.0E‐2‐1.5E‐2‐2.0E‐2‐2.5E‐2
1.E‐5

1.E‐2

1.E+1
Energy [eV]

1.E+4

1.E+7

10

理論
 罰則化線形回帰
 残差 + 罰則項を最小化

 罰則項：不自然に大きな解は選択しない

1
2

𝑟 𝑔⃗

𝛥𝑅

Δ𝑇 𝑔

疎な解
or
𝑔
𝑔

疎な解

密な解

ridge

lasso

𝜆
2

𝑝 𝑔⃗

𝑝 𝑔⃗

𝑝 𝑔⃗

𝑔

smooth‐lasso

𝜆

𝑔

𝑝 𝑔⃗

𝜆

𝑔

𝜆
2

𝑔

𝑔
11

理論
 Adaptive smooth-lasso
 smooth-lasso の 𝑝 𝑔
⃗

𝑝 𝑔⃗

𝜆

に重みづけ
𝜆
2

𝑢 𝑔

𝑣 𝑔

𝑔

 罰則を弱める、大きい値、急峻な変化を許す

 大きい感度 or 変化が大きい感度で小さい値

𝑢

1

1
𝑐𝑔

,𝑣

1

1
𝑐 0.5 ⋅ 𝑔𝑔

 smooth-lasso で粗く推定→重み決定

 重み: (0, 1] の範囲
12

検証計算
 ADS加速器ビーム電流値
 keff = 0.97, 未燃焼 熱出力: 800MW
 断面積数: 13286

 26核種、7反応、73群 サンプリング：±0.5% で一様分布
 参照値: 中央差分近似

 13286 x 2 回 forward 計算
 重みパラメータ: 参照値と合うように設定

𝜆0.01, 𝜆0.4, 𝑐 50.0, 𝛾
 自動最適化は今後の課題

0.8

 lasso, smooth-lasso も比較 lasso: 𝜆0.0, 𝑐 smooth-lasso: , 𝑐

0.0
0.0
13

検証計算
 相対感度係数 (750サンプル)

Relative sensitivity coefficient [‐]Relative sensitivity coefficient [‐] 上：lasso, 下：参照値
6.0E‐1
4.0E‐1
2.0E‐1
0.0E+0‐2.0E‐1‐4.0E‐1‐6.0E‐1‐8.0E‐1‐1.0E+0‐1.2E+0
6.0E‐1 0

1022

2044

3066

4088

5110

6132
7154
Cross section index

8176

9198

10220

11242

12264

13286

1022

2044

3066

4088

5110

6132
7154
Cross section index

8176

9198

10220

11242

12264

13286

4.0E‐1
2.0E‐1
0.0E+0‐2.0E‐1‐4.0E‐1‐6.0E‐1‐8.0E‐1‐1.0E+0‐1.2E+0
0

14

検証計算
 相対感度係数 (750サンプル)

Relative sensitivity coefficient [‐]Relative sensitivity coefficient [‐] 上：smooth-lasso, 下：参照値
6.0E‐1
4.0E‐1
2.0E‐1
0.0E+0‐2.0E‐1‐4.0E‐1‐6.0E‐1‐8.0E‐1‐1.0E+0‐1.2E+0
6.0E‐1 0

1022

2044

3066

4088

5110

6132
7154
Cross section index

8176

9198

10220

11242

12264

13286

1022

2044

3066

4088

5110

6132
7154
Cross section index

8176

9198

10220

11242

12264

13286

4.0E‐1
2.0E‐1
0.0E+0‐2.0E‐1‐4.0E‐1‐6.0E‐1‐8.0E‐1‐1.0E+0‐1.2E+0015

検証計算
 相対感度係数 (750サンプル)

Relative sensitivity coefficient [‐]Relative sensitivity coefficient [‐] 上：提案法, 下：参照値
6.0E‐1
4.0E‐1
2.0E‐1
0.0E+0‐2.0E‐1‐4.0E‐1‐6.0E‐1‐8.0E‐1‐1.0E+0‐1.2E+0
6.0E‐1 0

1022

2044

3066

4088

5110

6132
7154
Cross section index

8176

9198

10220

11242

12264

13286

1022

2044

3066

4088

5110

6132
7154
Cross section index

8176

9198

10220

11242

12264

13286

4.0E‐1
2.0E‐1
0.0E+0‐2.0E‐1‐4.0E‐1‐6.0E‐1‐8.0E‐1‐1.0E+0‐1.2E+0
0

16

検証計算
 Pu-239 Fis. (750サンプル)

17

検証計算
 Pu-239 Fis. (750サンプル) lasso: 特に裾で過小評価

18

検証計算
 Pu-239 Fis. (750サンプル) lasso: 特に裾で過小評価
 s. lasso : 不自然に滑らか、過小評価

19

検証計算
 Pu-239 Fis. (750サンプル) lasso: 特に裾で過小評価
 s. lasso : 不自然に滑らか、過小評価
 提案法 : 形状、大きさともによく再現

20

まとめと今後の課題
 RSによる感度係数評価の改良手法を提案
 核特性の感度係数の特徴を考慮
 論文[1]では keff のみ検討

 ADSビーム電流値の感度係数評価
 固定源+拘束条件(定格出力)
 提案手法の性能を示した。

 今後の課題 パラメータ自動最適化
 複雑な炉心解析、様々な核特性、核特性以外の検討

21

軽水炉における水の熱中性子散乱則データの違いによる核計算結果への影響評価
Study on the Impact of Thermal Scattering Law Data Improvements on Nuclear Calculation for Light
Water Reactor
1

敏 ，牟田 幹彦 ，山口 光 ，北田 孝典 1，大岡 靖典 2，松岡 正悟 2，長野 浩明 2

＊

竹田

1

1

1

大阪大学，2 原子燃料工業株式会社

JENDL-4.0、ENDF/B-VIII、JEFF-3.2、阿部氏の作成された水の熱中性子散乱則データ[1]を用いて、水の熱中
性子散乱則データを変更したことによる軽水炉での核計算結果への影響を評価した。
キーワード：軽水炉，熱中性子散乱則，S(α, β)，JENDL-4.0，JEFF-3.2，ENDF/B-VIII
1. 緒言
水の熱中性子散乱則データ S(α, β)は軽水炉において重要な核データである。JENDL-4.0 に含まれる水のS(α, β)は 1960 年代から検討されてきた GA モデルに基づいて作成され、JEFF-3.2 に含まれる水の S(α, β)は 1980 年代に発案された IKE モデルを改良して作成されている。ENDF/B-VIII では、近年検討されている分子動力学に基づいて水の S(α, β)が整備されている。日本においても、阿部氏により分子動力学に基づいた水の S(α, β)(以下、Abe)が整備されている[1]。上記に示したように、様々なモデルに基づいて水の S(α,β)が整備されていることをふまえ、本研究では水が多く含まれる軽水炉体系において、水の S(α, β)を変更したことによる核計算結果への影響を評価した。
2. 計算条件
ベースとなる評価済み核データライブラリは JENDL-4.0 とし、水の S(α, β)のみを変更して SRAC 詳細群
構造の 107 群断面積を作成した。評価対象は UO2 燃料ピンセル体系と MOX 燃料ピンセル体系とし、GMVPにより計算を実施した。
表 1 ピンセル体系における無限増倍率の比較
3. 計算結果と考察
UO2 燃料MOX 燃料Abe を用いた場合の無限増倍率と、他のAbe との
Abe との
k-infinity*
k-infinity*
相対差[%]
相対差[%]
S(α, β)を用いた場合の無限増倍率を比較した
1.39093―
1.17039―(表 1)。UO2 燃料ピンセル体系では無限増倍率 Abe
JENDL-4.0
1.39050-0.0311.16796-0.208の相対差は 0.05%以下となった。一方、MOX
JEFF-3.2
1.390760.019
1.17015-0.021
燃料ピンセル体系では、JENDL-4.0 との相対 ENDF/B-VIII
1.39097
0.015
1.170850.039
差が約 0.2%となることが確認された。
*統計誤差(1σ)は 0.00002 以下
MOX 燃料ピンセル体系での無限増倍率の差異をふまえ、MOX 燃料ピンセル体系で無限増倍率の差異を特に生じさせる入射エネルギー1eV の散乱断面積について、Abe から JENDL-4.0 の S(α, β)に変更した場
合の反応度を厳密摂動により評価した。散乱断面積の差を図 1 に、投入される反応度を図 2 に示す。図 1
及び図 2 より、1eV から 0.1eV 付近に下方散乱される断面積の差異は比較的小さいものの、反応度としては有意な影響を及ぼすことが明らかになった。
Reactivity per lethargy(⊿k/kk'lethargy)

Difference in crosssection (1/cm)

6.0E-03
1.2E-02 (1/cm)

4.0E-032.0E-03
0.0E+00-2.0E-03

1.0E-04
0.0E+00-1.0E-04-4.0E-03-6.0E-03
1.0E-03

1.0E-02

1.0E-01-2.0E-04
1.0E-03

1.0E+00

入射エネルギー1eV の散乱断面積差(JENDL-4.0 ― Abe)

1.0E-02

1.0E-01

1.0E+00

Energy of scattered neutron (eV)Energy of scattered neutron (eV)

図 1

2.0E-04

図 2

入射エネルギー1eV の散乱断面積変更により投入される反応度(Abe → JENDL-4.0)

4. まとめ
本検討では UO2 燃料及び MOX 燃料ピンセル体系において S(α, β)の違いによる影響を評価した。
他の S(α,β)と比べ、JENDL-4.0 では MOX 燃料体系での無限増倍率が約 0.2%低くなることが確認された。
参考文献[1] Y. Abe, T. Tsuboi, S. Tasaki, “Evaluation of the neutron scattering cross-section for light water by molecular dynamics,”
Nuclear Instruments and Methods in Physics Research Section A, Vol. 735, Pages 568-573, (2014).
*Satoshi TAKEDA1, Mikihiko MUTA1, Hikaru YAMAGUCHI1, Takanori KITADA1, Yasunori OHOKA2, Shogo MATSUOKA2,

Hiroaki NAGANO2. 1Osaka University, 2Nuclear Fuel Industries, Ltd.

第7回「炉物理専門研究会」
2018年11月26日, 京都大学 複合原子力科学研究所

水素の熱中性子散乱則の違いによる軽水炉体
系での影響*竹田 敏1、牟田 幹彦1、山口 光1、北田 孝典1、
大岡 靖典2、松岡 正悟2、長野 浩明2（1. 大阪大学、2. 原子燃料工業）
1

背景～水の熱中性子散乱則データの開発～


熱中性子散乱則データは評価済み核データライブラリで整備されてきている。
 ENDF/B： ENDF/B-IIIから熱中性子散乱則のデータが整備されており、
ENDF/B-VIIIでは分子動力学のシミュレーション結果に基づき整備(CABモデル)
 JENDL： JENDL-4.0ではENDF/B-VI.8から引用(GAモデル)
 JEFF： JEFF-3.1.1にて熱中性子散乱則データを整備(IKEモデル)
 近年では京大安部先生が分子動力学のシミュレーション結果に基づき整備*
ENDF/B

JENDL

熱中性子散乱則データ
ENDF/B-III
水の分子振動を考慮するためのモデルを整備(GAモデル)
ENDF/B-VI.2 実験値をふまえてモデルに用いるパラメータを整備(GAモデル)
ENDF/B-VI.8 ↑ENDF/B-VII.0 水分子のクラスタを考慮するなど、モデルを大幅に変更(IKEモデル)
ENDF/B-VII.1 ↑ENDF/B-VIII
分子動力学に基づきパラメータを整備するモデルを採用(CABモデル)

熱中性子散乱則データ
JENDL-3.3 評価済み核データライブラ
リの中には存在しない。
JENDL-4.0 ENDF/B-VI.8から引用(元はENDF/B-VI.2)

JEFF
熱中性子散乱則データ
JEFF-3.1.1 IKEモデルに基づき整備
JEFF-3.2↑*Y.Abe, T.Tsuboi,S.Tasaki, “Evaluation of the neutron scattering cross-section for light water,” 735, 568–573, Nuclear
Instruments and Methods in Physics Research A, 2014.

2

背景～分子動力学の応用の違い～
ENDF/B-VIIIと安部先生の熱中性子散乱則データはどちらも分子動力学に基づくシ
ミュレーション結果を応用している。



ENDF/B-VIII

散乱則を非干渉性散乱近似で表し、ガウス近似により表された自己中間散乱
関数を並進運動，分子回転，分子内振動に分離してモデル化している。

分子回転

並進運動

分子内振動

各モデルの入力パラメータとして、分子動力学のシミュレーション結果を活用
安部先生


散乱則を非干渉性散乱近似で表し、ガウス近似により表された自己中間散乱関数で用いられる速度の自己相関関数を分子動力学のシミュレーションにより直接求める。
S 𝜅, 𝜔

1
2𝜋

𝑒

exp

𝜅 ℏ
2𝑚

𝑔 𝜔
𝜔

𝑐𝑜𝑡ℎ

ℏ𝜔
2𝑘𝑇

1

cos 𝜔𝑡𝑖𝑠𝑖𝑛 𝜔𝑡𝑑𝜔𝑔 𝜔 :速度の自己相関関数を余弦変換させた関数
𝜅：散乱前後の運動量変化, 𝜔:散乱前後のエネルギー変化

3

本検討の目的
軽水炉において核計算結果は水の熱中性子散乱則に強く依存
安部先生、ENDF/B‐VIII、JEFF‐3.2、JENDL‐4.0の熱中性子散乱則データを用いて、水の熱中性子散乱則データを変更したことによる核計算結果への影響を評価

ベースとなる評価済み核データライブラリはJENDL‐4.0として、水の熱中性子散乱則
データのみ変更させて定数を用意し、結果を比較する。
安部先生S(α, β)

ENDF/B‐VIIIJENDL‐4.0

JEFF‐3.2

JENDL‐4.0
H‐1, O‐16断面積

JENDL‐4.0
H‐1, O‐16断面積

JENDL‐4.0
H‐1, O‐16断面積

JENDL‐4.0
H‐1, O‐16断面積

安部先生
H‐1(in H2O)
S(α, β)ENDF/B‐VIII
H‐1(in H2O)
S(α, β)JENDL‐4.0
H‐1(in H2O)
S(α, β)JEFF‐3.2
H‐1(in H2O)
S(α, β)

4

計算方法


GMVPにより107群のBWR燃料セル体系を計算する。



燃料・被覆管の断面積についてはJENDL-4.0に基づきSRAC2006を用いて作成する。



水の断面積については、NJOY(ver. 99.396)で多群断面積を作成し、断面積変換プロ
グラムTRANSX2.15を用いてGMVP用の断面積形式に変換させる。
燃料・被覆管等の断面積
SRAC2006
 PijモジュールでGMVPと同様の体系で計算
 PEACO
水の断面積
NJOY + TRANSX2.15

燃料セル体系

減速材温度 [K]燃料温度 [K]エネルギー群数
有効バッチ数スキップバッチ数
ヒストリー数/バッチ

GMVP

0.480
：Fuel
：Cladding 0.550：Coolant 0.725

600
900
107
1万
0.1万
4万

評価する燃料 4.45wt%UO2燃料
 MOX燃料(富化度10wt%)(cm)

Pu‐239Pu‐240Pu‐241
Pu‐242

45%
30%
15%
10%

5

結果～ウラン燃料体系での無限増倍率～


ウラン燃料体系におけるボイド率0, 40, 70%での無限増倍率を求めた。各ボイド率について、安部先生の熱中性子散乱則データを用いる場合（Abe）を基準として、無限増倍率の比較を行った。



ENDF/B-VIIIとJEFF-3.2についてはAbeとの相対差は0.010%程度以下であった。



JENDL‐4.0との相対差は比較的大きく0.030%程度以下となった。

① Abe② ENDF/B‐VIII③ JENDL‐4.0④ JEFF‐3.2(② ‐ ①) / ① (%)(③ ‐ ①) / ① (%)(④ ‐ ①) / ① (%)

Void = 0%1.39093
1.39097
1.39050
1.39076
0.003‐0.031‐0.012

Void = 40%
1.25200
1.25205
1.251631.251880.005‐0.030‐0.010

Void = 70%
1.04144
1.04153
1.04123
1.041450.009‐0.020+0.001※無限増倍率の標準偏差σは0.002%以下

6結果～MOX燃料体系での無限増倍率～


MOX燃料体系についても同様に、Abeを基準として、無限増倍率の比較を行った。



ENDF/B-VIIIとJEFF-3.2についてはAbeとの相対差の大きさは0.040%程度以下となり、ウラン燃料よりも相対差は大きくなる傾向が確認された。



JENDL-4.0については特に相対差が大きくなり、ボイド率0%の場合、Abeとの相対差の大きさは0.200%を超えた。

① Abe② ENDF/B‐VIII③ JENDL‐4.0④ JEFF‐3.2(② ‐ ①) / ① (%)(③ ‐ ①) / ① (%)(④ ‐ ①) / ① (%)

Void = 0%1.17039
1.17085
1.16796
1.17015
0.039‐0.208‐0.021

Void = 40%
1.09226
1.092631.09040
1.09204
0.034‐0.170‐0.020

Void = 70%
1.03691
1.03715
1.036111.03689
0.023‐0.077‐0.002※無限増倍率の標準偏差σは0.002%以下
AbeとJENDL-4.0の結果の差異について考察

7

考察～摂動論～

AbeとJENDL-4.0の無限増倍率を比較した結果、ウラン燃料体系では相対差が小さく、MOX燃料体系では相対差が大きくなる傾向が確認された。



散乱断面積の違いが反応度に与える影響を考察するため、摂動論により反応度を評価した。
AbeからJENDL-4.0に散乱断面積を変更することを考えると、反応度の散乱成分は下記のように示すことができる。

∑

𝑑𝑟 𝜙 𝑟 ∑
∑

Σ
𝑑𝑟 𝜒

∗ →

𝑟

Σ

→

𝑟 𝜑∗ 𝑟 ∑ υΣ

𝜑∗

𝑟

𝑟 𝜙

𝑟

𝑟

𝜑∗ 𝑟

𝜙：中性子束𝜑：随伴中性子束*：摂動後(=JENDL‐4.0の場合)のパラメータを示す。


散乱成分を分析することで、MOX体系で実効増倍率の相対差が大きくなった
原因を調査し、反応度に与える影響が大きい散乱断面積差異を明らかにする。

8

考察～ウラン燃料/MOX燃料の反応度差異～摂動論により散乱断面積の変更(Abe→JENDL-4.0)により投入される反応度を求め、投入される反応度を中性子の入射エネルギーで整理した。



ウラン燃料と比較して、MOX燃料では1eV付近以上の散乱断面積を変更すること
により大きな負の反応度が投入されることが確認された。
Reactivity per lethargy(⊿k/kk'lethargy)2.0E-03
1.5E-03
1.0E-035.0E-04
0.0E+00-5.0E-04-1.0E-03-1.5E-03-2.0E-03
1.0E-03

UOX

MOX

1.0E-02 1.0E-01
Energy (eV)

1.0E+00

散乱断面積の変更(Abe→JENDL-4.0)により投入される反応度9

考察～反応度の算出式～


MOX燃料では、1eV付近以上の散乱断面積が変更されることにより負の反応度が投入されることが分かった。この要因を考察するため、反応度の算出式に立ち戻る。



規格化された条件下では、反応度は中性子束(𝜙 𝑟 )、散乱断面積の差
∗ →

(Σ
𝑟
求まる。

Σ

→

𝑟 )、散乱前後の随伴中性子束の差(𝜑 ∗中性子束∑

𝑑𝑟 𝜙 𝑟 ∑∑

散乱断面積の差
Σ
𝑑𝑟 𝜒

∗ →

→

𝜑∗

𝑟

𝑟 𝜑 ∗ 𝑟 ∑ υΣ

𝑟 𝜙

𝑟

Σ

𝜑∗ 𝑟 )の積で随伴中性子束の差
𝑟

𝑟

𝑟

𝜑∗ 𝑟

上記の三つのうち、散乱断面積の差(AbeとJENDL-4.0の散乱断面積差)について
は、ウラン燃料とMOX燃料で同等となる。このため、ウラン燃料とMOX燃料で反応度が異なる要因としては、中性子束と随伴中性子束が考えられる。

10

考察～中性子束の確認～

ウラン燃料体系とMOX燃料体系の水領域の中性子束を比較した。


比較のため、中性子束は体系における生成率が1になるよう規格化している。
υΣ

𝑟 𝑑𝑟

1

ウラン燃料と比較して、MOX燃料では1eV付近以上の中性子束が同程度以下となることが確認された。このため、中性子束の大きさの違いは、反応度の違いを引き起こす主要因でないことがわかった。
Normalized forward flux(/lethargy)



𝑟 𝜙

1
UOX
MOX

0.8
0.6
0.4
0.2
0
1.0E-03

1.0E-02

1.0E-01
Energy (eV)

1.0E+00

水領域における中性子束の比較

11

考察～随伴中性子束の確認～

ウラン燃料体系とMOX燃料体系の水領域の随伴中性子束を比較した。


比較のため、随伴中性子束は以下の式により規格化している。
𝜒

𝑟 𝜑

𝑟 𝑑𝑟

1



ウラン燃料と比較して、MOX燃料では1eV付近の前後で随伴中性子束が変化する
ことが確認できる。

𝜑∗ 𝑟 )により計算できること
反応度は散乱前後の随伴中性子束の差(𝜑 ∗ 𝑟から、散乱前後での随伴中性子束の差により、MOX燃料ではウラン燃料と比較して大きな反応度が投入されることがわかった。
Normalized adjoint flux

2.5
2
1.5
1
0.5
0
1.0E-03

UOX
MOX
1.0E-02

1.0E-01
Energy (eV)

1.0E+00

水領域における随伴中性子束の比較12

考察～散乱後のエネルギーでの分解①～AbeからJENDL-4.0へ散乱断面積を変更することにより投入される反応度を分析した結果、MOX燃料では1eV付近の随伴中性子束の変化傾向により反応度が投入されることが確認できた。



散乱前後の随伴中性子束の差が反応度に寄与することをふまえ、入射エネルギー
1eV付近(0.99 - 1.13 eV)の散乱断面積を変更した場合に投入される反応度について、散乱後のエネルギーで整理して考察を深める。
Reactivity per lethargy(⊿k/kk'lethargy)2.0E-03
1.5E-03
1.0E-035.0E-04
0.0E+00-5.0E-04-1.0E-03-1.5E-03-2.0E-03
1.0E-03

反応度を散乱後のエネルギーで整理

1.0E-02 1.0E-01
Energy (eV)

1.0E+00

散乱断面積の変更(Abe→JENDL-4.0)により投入される反応度13

考察～散乱後のエネルギーでの分解②～

散乱断面積の差(JENDL-4.0 - Abe)としては、1eV付近で大きく、1→0.1eVの散乱断面積の差は相対的に小さいことが確認できる。
1.8E-01
1.6E-01
1.4E-01
1.2E-01
1.0E-01
8.0E-02
6.0E-02
4.0E-022.0E-02
0.0E+00
1.0E-03

Difference in cross section(1/cm)

入射エネルギーが1eVの散乱断面積を確認した。断面積の大きさとしては1eV付近
でピークをもつことが確認できる。

Cross section (1/cm)

Abe
JENDL-4.0

1.0E-02

1.0E-01

1.0E+00

Energy of scattered neutron (eV)

入射エネルギー1eV付近の散乱断面積

4.0E-033.0E-03

1.2E-02 (1/cm)

2.0E-03
1.0E-03
0.0E+00-1.0E-03-2.0E-03

JENDL-4.0 - Abe-3.0E-03-4.0E-03
1.0E-03

1.0E-02

1.0E-01

1.0E+00

Energy of scattered neutron (eV)

入射エネルギー1eV付近の散乱断面積の差

14

考察～散乱後のエネルギーでの分解③～
入射エネルギー1eV付近の散乱断面積を変更した場合に投入される反応度を評価すると、 1→0.1eV 付近の散乱断面積の変更による反応度が負側に大きいことが確認できる。



散乱断面積の差が小さい場合でも、散乱前後で随伴中性子束の差が大きくなる場
合は、投入される反応度が大きくなることが確認された。
1.5E-04

Difference in cross section(1/cm)

Reactivity per lethargy(⊿k/kk'lethargy)

1.0E-04
5.0E-05
0.0E+00-5.0E-05-1.0E-04-1.5E-04
1.0E-03

1.0E-02

1.0E-01

1.0E+00

Energy of scattered neutron (eV)

入射エネルギー1eV付近の散乱断面積を変更(Abe→JENDL-4.0)した場合に投入される反応度

4.0E-033.0E-03

1.2E-02 (1/cm)

2.0E-03
1.0E-03
0.0E+00-1.0E-03-2.0E-03

JENDL-4.0 - Abe-3.0E-03-4.0E-03
1.0E-03

1.0E-02

1.0E-01

1.0E+00

Energy of scattered neutron (eV)

入射エネルギー1eV付近の散乱断面積の差

15

まとめ

安部先生、ENDF/B‐VIII、JEFF‐3.2、JENDL‐4.0の熱中性子散乱則を用いて、水の熱中性子散乱則データを変更したことによる核計算結果への影響を評価した。



MOX燃料体系において、安部先生の熱中性子散乱則データを用いた場合とJENDL-4.0のデータを用いた場合の結果を比較すると、ボイド率0%の場合では無限増倍率の相対差の大きさは0.200%を超えた。



摂動論により、用いる熱中性子散乱則データを安部先生のデータからJENDL-4.0に変更した場合の反応度を分析した。分析の結果、MOX燃料では1eV付近の随伴中性子束の変化傾向により比較的大きな反応度が投入されることが確認できた。



さらに、散乱先のエネルギーで整理して分析した結果、1eVから0.1eV付近への散乱断面積の差異について、差異の大きさは比較的小さいものの、反応度へは有意な影響を及ぼすことが明らかになった。

16

KUCA 固体減速架台の数値計算ベンチマーク問題
Benchmark problems for neutron transport calculations about KUCA solid-moderated cores＊1

千葉 豪 1 , 遠藤 知弘 2 ,北海道大学，2 名古屋大学

KUCA 固体減速架台における濃縮ウラン燃料を用いた臨界実験体系に関する数値計算ベンチマーク問題の整備状況について報告する。
キーワード：KUCA、固体減速架台、濃縮ウラン、ベンチマーク問題
1.

背景：京都大学臨界実験装置（KUCA）の固体減速架台では近年、加速器駆動システム（ADS）の核特性に関わる基礎データが精力的に取得されており、そのいくつかは公開され ADS 開発に関わる研究者に広く利用されている（文献[1]など）。公開データには、余剰反応度や反応率分布に加えて、ゾーン物質置換反応度や動特性に関わるデータも含まれており、その実験解析においては統計論的手法と同様に決定論的手法も有用となる。決定論的手法による KUCA 固体減速架台の実験解析においては、実験体系のモデル化・簡略化が必要となるため、その影響を定量的に把握しておくことが重要となる。そのような観点から、現在、KUCA 固体減速架台の数値計算のためのベンチマーク問題の整備を進めている。
2.

ベンチマーク問題：1/16 インチの濃縮ウラン（EU）板 2 枚と 1/8 イ

ンチのポリエチレン（PE）板 1 枚からなる単位燃料セル、及びポリエチ
レン板で構成される単位反射体セルの 2 種類を考え、それらがアルミニ
ウム製の鞘管に格納されて炉心に配置されるものとした。なお、固体減速架台特有の鞘管間の空隙も考慮した。XY 平面上での炉心配置図を図1 に示す。この体系を Phase III とした。また、これに加えて、鞘管を EU
板、PE 板に希釈させて単位燃料セルを 1 次元でモデル化可能としたPhase II、Phase II において鞘管間の空隙を無視した Phase I、Phase I において単位燃料セルを体積荷重で平均化した Phase 0 の体系を設定した。
各 Phase について連続エネルギーモンテカルロコードにより参照解を与えることで、ベンチマーク問題の利用者が自らの決定論的手法の妥当性

図1

Phase III の炉心配置図

をより詳細に吟味することが可能となる。
3.

解析例：連続エネルギーモンテカルロコード MVP による解析結果を表 1 に示す。使用した核データは

JENDL-4.0 である。単位燃料セルを均質モデルから 1 次元モデルにすることにより実効増倍率は 0.6%程度低下した。これは断面積の空間的な自己遮蔽効果によりセル平均の核分裂断面積が低下することに由来する。
また、鞘間の空隙を考慮することにより 0.8%程度の実効増倍率の減少が見られた。2mm 程度の空隙であるが、中性子のストリーミングパスとなる
ため、中性子の漏洩が大きくなり実効増
倍率が減少したものと考えられる。単位
燃料セルを 1 次元から 2 次元モデルにすることによる実効増倍率への影響は 0.1%

表1

Phase
0
I
II
III

モンテカルロコードによる解析結果

MVP1.03457
1.02859
1.020701.02028未満であった。当日は決定論コードによる解析例も併せて示す。
参考文献：[1] C.H.Pyeon, KURRI-TR-444 (2012).
* Go Chiba1, Tomohiro Endo2
1 Hokkaido Univ., 2 Nagoya Univ.,

Homogeneous cell model, no gap
1D cell model, no gap
1D cell model(Reference)

KUCA固体減速架台の数値計算ベンチマーク問題
北大 千葉 豪名大 遠藤 知弘

1

KUCA固体減速架台（A・B架台）

今回の発表はA架台に関するものです。
2

KUCA軽水減速架台（C架台）

3

A架台における種々のADS核特性に関する実験データ

100MeV陽子入射の未臨界実験により以下の核特性パラメータが取得：
・絶対反応率・Cd比・反応率空間分布・即発減衰定数・未臨界度
4A架台における種々のADS核特性に関する実験データ

鉛・ビスマスの置換反
応度実験

*C.H.Pyeon, et al., JNST, 53, p.602‐612 (2016)5

実験解析を行うにあたって

遠藤知弘、2018年炉物理夏期セミナー講演資料

6

実験解析を行うにあたって

連続エネルギーのモンテカルロ計算でこれらは（統計誤差付きで）精度良く考慮できるので、それで良い？

遠藤知弘、2018年炉物理夏期セミナー講演資料

7

実験解析を行うにあたって

遠藤知弘、2018年炉物理夏期セミナー講演資料

8実験解析を行うにあたって

小回りの効く決定論的手法の役割はまだまだ残っており、決定論的手法でKUCAの実験解析を行う人のためのベンチマーク問題があると便利と考えた。

遠藤知弘、2018年炉物理夏期セミナー講演資料

9

今回設定したベンチマーク問題

燃料と反射体からなる体系を設定。
10

今回設定したベンチマーク問題

単位燃料格子は、アルミニウムの鞘管に93%濃縮ウラン・Al合金板とポリエチレン板が規則的に配列しているものとする。
11

U‐235の反応断面積

単位燃料格子の均質モデルで評価すると、核分裂反応率は捕獲反応率の4倍程度。
12

HEU燃料を用いた実験での随伴中性子束エネルギースペクトル

濃縮ウラン燃料使用、かつ高漏洩体系であることにより、エネル
ギーに対して単調に減少する形状となっている。

13

今回設定したベンチマーク問題
Phase 0Phase I

Phase II

Phase III

0D cell

1D unit cell

1D unit cell

2D unit cell

Phase IIIに対して、モデルを簡略化したPhase 0～IIを追加で設定（Phase IVでは制御棒を導入する予定）。
14

連続エネルギーモンテカルロコードによる参照解Phase 0Phase I

Phase II

Phase III

0D cell

1D unit cell

1D unit cell

2D unit cell1.03486
1.03384

1.028881.02794

1.021111.02004

1.01951
1.01863

上段はMVP‐3、下段はMCNP‐6.2計算値。JENDL‐4.0使用。
統計誤差は0.00004程度。予稿の数値と異なることに注意。
15

連続エネルギーモンテカルロコードによる参照解Phase 0Phase I

Phase II

Phase III

0D cell

1D unit cell

1D unit cell

2D unit cell1.03486
1.03384

1.028881.02794

1.021111.02004

1.01951
1.01863

MVP‐3とMCNP‐6.2とで統計誤差（0.00004）を超える有意な差が見られた → 現在確認中。

16

連続エネルギーモンテカルロコードによる参照解Phase 0Phase I

Phase II

Phase III

0D cell

1D unit cell

1D unit cell

2D unit cell1.03486
1.03384

1.028881.02794

1.021111.02004

1.01951
1.01863

均質格子モデルから非均質格子モデルに変更することで実効増倍率が低下。
17

連続エネルギーモンテカルロコードによる参照解Phase 0Phase I

Phase II

Phase III

0D cell

1D unit cell

1D unit cell

2D unit cell1.03486
1.03384

1.028881.02794

1.021111.02004

1.01951
1.018632mmの空隙を考慮することで0.8%弱もの実効増倍率の低下が見られた。
18

連続エネルギーモンテカルロコードによる参照解Phase 0Phase I

Phase II

Phase III

0D cell

1D unit cell

1D unit cell

2D unit cell1.03486
1.03384

1.028881.02794

1.021111.02004

1.01951
1.01863

格子モデルで鞘管を陽に考慮することで実効増倍率が低下。主に反射体領域のモデル化の影響の模様。
19

CBZによる計算

・JENDL‐4.0に基づく107群のライブラリを使用。
・単位燃料セルは１次元平板でモデル化し、実効断面積は東捻の方法で計算。衝突確率法で中性子束分布を計算し均質化を行う。
・KUCAのA架台特有のギャップは炉心計算で直接考慮。
・炉心計算は三次元Sn輸送計算で実施（CBZ/SNT）。

20

計算結果：Phase 0

非等方散乱についてはP2程度まで考慮する必要がある。
21

計算結果：Phase 0

角度求積背セットについてはLS12程度まで考慮すればよい。P2S12で0.4%程度の過少評価が見られる（後で議論）。
22

計算結果：Phase I

P2計算では参照値との差異が0.6%程度まで拡大→均質化の影響？

23

計算結果：Phase I

均質化を行わない場合は、参照値との差異は0.4%程度となった。ただし、かなり詳細な角度分点が必要。
24

角度求積セットLevel symmetric（黒点がEven moment condition、白点
はLevel moment condition）

Legendre‐Chebyshev productquadrature

・Product quadratureについては、Double Legendre‐Chebyshev
quadratureを使用。「DPnTn」という呼称とし、この場合は1/8象限に対してn×nの離散点が配置されることとする。
*C.Rabiti, at al., ANS transaction, 96, 2007. 25

計算結果：Phase I

空間に関しても、詳細なメッシュ分割が必要となる。
26

計算結果：Phase II

CBZのバイアスはPhase Iと同様に0.6%程度のため、0.2mmの間隙による中性子ストリーミング効果は決定論で考慮できていると考
えられる。
27

計算結果：Phase 0

このシンプルな問題での0.4%の過小評価の原因は何だろうか？
28

計算結果：Phase 0

実効断面積計算（共鳴処理）が第一に考えられるので、まずは共鳴処理の影響がどの程度あるか、
無限希釈断面積で計算を行ったところ、‐0.06%程度の影響しかなかった。
→ HEU燃料でU‐235しか含まれていないので、共鳴処理はあまり影響しない？

このシンプルな問題での0.4%の過小評価の原因は何だろうか？
29

U‐235の反応断面積

単位燃料格子の均質モデルで評価すると、核分裂反応率は捕獲反応率の4倍程度。
30

計算結果：Phase 0
実効断面積計算（共鳴処理）が第一に考えられるので、まずは共鳴処理の影響がどの程度あるか、
無限希釈断面積で計算を行ったところ、‐0.06%程度の影響しかなかった。
念のため、実効断面積を無限希釈断面積に置換したときの反応度を摂動計算により要因別に求めた
ところ、以下の結果を得た。
・生成項：+1.09%
・吸収項：‐1.60%
・散乱項＋漏洩項：+0.45%31

摂動計算結果：生成項と吸収項

U‐235の共鳴計算が重要である。

・生成項：+1.09%
・吸収項：‐1.60%
・散乱項＋漏洩項：+0.45%32

摂動計算結果：散乱項と漏洩項・生成項：+1.09%
・吸収項：‐1.60%
・散乱項＋漏洩項：+0.45%
33

Al‐27の全反応断面積
Al‐27の共鳴計算も重要である。

34

まとめ・KUCAの固体減速架台（A架台）における高濃縮
ウランを用いた実験の数値解析のためのベンチ
マーク問題を作成している。
・決定論コードを用いた実験解析にとって有用となるものと考える。
・高濃縮ウラン燃料を使用した高漏洩体系であり、他の原子炉とは異なる特性が観察された。
・高エネルギー領域の共鳴計算方法が重要となることが示唆された。
35

Effect of Differences Fuel Material on Neutronic Parameters in Kartini Research Reactor*Argo Satrio Wicaksono1, Takanori Kitada1, Satoshi Takeda1
1

Osaka University.

This research offers a possibility fuel material conversion from TRIGA standard fuel assembly, U-ZrH1.65 to others
fuel materials for solving the operation problem in many TRIGA reactors. In this work, neutronic parametersassessment of TRIGA Kartini reactor’s core is loaded by U3Si2-Al, U3O8-Al, and UO2 then simulated using
SRAC2006. The neutronic parameters of each fuel materials were compared based on the calculation of parameters
for utilization of reactor and safety operation.
Keywords: TRIGA reactor, SRAC2006, neutronic parameters
1. Introduction: Fuel fabrication of TRIGA standard fuel assembly has been suspended, it is important to get asolution for TRIGA reactor with alternative fuel assembly. A reliable estimation of neutronic parameters, for each offuel materials will play important role in decision making of selecting suitable alternative fuel materials. In order to
get exact information about core characteristic, the standard reactor analysis code SRAC[1] was chosen because of its
general geometry modeling capability, and provide multi-dimensional diffusion code. The objective in this research is
to study the possibility of TRIGA Kartini reactor’s core that loaded by others fuel materials and to select the optimum
based on the calculation of neutronic parameters for utilization of reactor and safety operation.
2. Methodology: The collision probability method lattice transport
code SRAC-PIJ and the SRAC-CITATION code are used to perform
global core calculations of Kartini reactor. The TRIGA standard fuelgeometry and original core configuration are used for U-ZrH1.65,
U3Si2-Al, U3O8-Al, and UO2 respectively. Since the alternative fuelmaterials are in the under-moderated situation, fuel geometry changed
by increasing the Moderator To Fuel Ratio (MTFR) from 0.226-4.291(Fig.1). Neutronic analysis of the U3Si2-Al and U3O8-Al are carriedout when the reactor core is in a critical condition, MTFR = 0.330.

Fig.1 Relation between MTFR and K-eff

3. Result and Discussion: Neutronic parameters of U3Si2-Al andU3O8-Al, as alternative fuel materials have been compared with
U-ZrH1.65 as the TRIGA standard fuel material. Safety parameters suchas power peaking factor, core excess, and shutdown margin have been
compared with the limit value in the safety analysis report. Total flux
distribution at irradiation facilities (Fig.2) shows the higher result in
U3Si2-Al and U3O8-Al than U-ZrH1.65 materials. However, at the centercore region and edge core position still needs moderation to change
fast to be thermal neutron by increasing the MTFR at fuel regions.

Fig.2 Flux distribution at irradiation facilities

4. Conclusion: The parameters for utilization of reactor and safety
operation for each fuel materials are evaluated using a comprehensive neutronic code, SRAC2006. U3Si2-Al andU3O8-Al, can be utilized as TRIGA Kartini reactor alternative fuel materials by considering MTFR.
References [1] K. Okumura, T. Kugo, K. Kaneko, and K. Tsuchihashi, SRAC2006: A Comprehensive Neutronics Calculation
Code System, JAEA-Data/Code 2007-004, 2007.

Reactor Physics Workshop
November 26-27, 2018
Kyoto University Critical Assembly

Effect of Differences Fuel Material on Neutronic
Parameters in Kartini Research Reactor〇 Argo S. Wicaksono, Takanori Kitada, Satoshi Takeda
Division of Suistanable Energy and Environmental Engineering,
Graduate School of Engineering, Osaka University, Japan

1

Contents
I. Introduction
II. Calculation techniques and tools
III. Result and Discussion
IV. Summary
V. Future Work

2

I.Introduction

Kartini Research Reactor Profile

1 of 70 TRIGA reactors in the world (35 in US, 3 in JPN).

Commissioned at the National Nuclear Energy Agency
of Republic Indonesia (BATAN), Yogyakarta on March
01, 1979.

Operates steady state at 100 kwth and since 1994
using TRIGA stadandard fuel elements U‐ZrH (type 104
and 204) manufactured by General Atomics (US).



Used for neutron activation analysis, sample
irradiation, and nuclear physics experimentutilizations.
TRIGA Kartini Reactor Core

3

I.Introduction

Reactor Fuel Materials and TRIGA Kartini Problems~ Fuel Shortage ~Fuel fabrication of TRIGA standard fuel assembly has been suspended.
~ Requirements for alternative fuel assembly ~Fabricated in domestically or many are produced globally.
Material Compound

U‐ZrH1.65

U3Si2‐AlU3O8‐Al

UO2

Utilizing

TRIGA

MTR

MTR

LWR

FabricatedUS (suspended) INA, etc. INA, etc.

US, JPN, etc.

~ Desirable alternative fuel assembly ~I. Sufficient excess reactivity.
II. Produces the high neutron fluxes.
III. Power peaking factor for thermal-hydraulicIV. Nuclear fuel life time for the long operation period.

4

I.Introduction

Kartini Reactor Core Conf. and Fuel Assembly Dim.
Fuel Assembly Cross Section
Zr‐Rod [Outer radius=0.288 Cm]
Fuel Meat [Outer radius= 1.780 Cm]
Helium Gap [Outer radius=1.790 Cm]
SS‐305 Cladding [Outer radius= 1.850 Cm]
Water coolant [Outer radius = 2.029 Cm]
TRIGA Std. Fuel Assembly

Kartini Reactor Core Cofiguration (2017)

Reactor Core Configuration
Fuel Assemblies [71 pcs]
Control Reactivity Rods [3 pcs]
Graphite Dummy Element [14 pcs]
Neutron Source [1 pcs]
Irradiation Facilities [2 pcs]
Water Coolant and Moderator
Graphite Annular Reflector [Outer radius= 30Cm]

5

II. Calculation techniques and Tools

Fuel Materials and Core parameters
Material Specifications

U‐ZrH1.65 U3Si2‐Al U3O8‐Al

UO2

U235 enrich. (%)

19.75%

19.75%

19.75%

3.0%

Material density (gr/cc)

5.96

12.2

8.3

10.96

Core Parameters
Thermal Power

100 kW

Active core dia.

43.81 Cm

Active core height

38 Cm

Annular reflector

30 Cm

Fuel Temperature

403 K

Water Temperature

313 K

6

II. Calculation techniques and Tools

Simulation Methodology
SRAC-PIJSRACCITATIONJENDL-4.0

E‐ring model
1 control rod
23 fuel rods

F‐ring model
14 dummy rods
13 fuel rods
IGT=3 Cylindrical cell

Central thimble, B and C‐ring model
IGT=11 Annular assembly with asymmetric pin rods

Reactor Core ½ Cylindrical R‐Z Geometry

Fuel assembly cell models

D‐ring model
all fuel rods

7

III. Result and Discussion

Multiplication factor (Reactor Criticality)

0.226

0.330

U‐ZrH1.65

1.0346‐

U3Si2‐Al

0.94691.0438

U3O8‐Al

0.95871.0540UO2

0.81000.8864

Reactor core K‐eff

Moderator to
Fuel Ratio

Under moderated

1.2

Fuel Material Multiplication Factor

Over moderated

1.1
Doppler (pcm/K) = ‐1.91 ; ‐2.04
MTC (pcm/K)= ‐11.46 ; ‐10.88

1

0.9

0.8
0

0.330
0.226

1

U3Si2‐Al

2

U3O8‐Al3

4

UO2

Moderator to fuel ratio

8III. Result and Discussion

Neutron Flux at Irradiation Facilities
U3Si2‐Al

2E+12

Fast

0
UZrH

Edge U3Si2Al
core position
U3O8AlEdge core
position
U3Si2AlU3O8AlUZrH

U3Si2‐AlU3O8‐Al

5E+11

U‐ZrH1.65

U‐ZrH1.65

1E+12

Kartini Reactor Core Whole Model

Epithermal

U3O8‐AlU3Si2‐Al

1.5E+12

ThermalU‐ZrH1.65

Edge core position

Rotary rack region

U3O8‐Al

Center core region

RotaryU3Si2Al
rack regionU3O8Al

UZrH

9

III. Result and Discussion

Power peaking analysis[b] Axial power peaking factor
fZ = Ppeak/ Pave= axial peak to average power ratio.
The total peaking factor fT defined as a product of hot
rod factor and axial peaking factor are found to be:
fT = fHR x fZType ofpeaking

Limit
value

U‐ZrH1.65

U3Si2‐AlU3O8‐AlfHR 1.7908

1.536

1.510

1.529

fZ 1.4050

1.281

1.155

1.166

fT 2.5161

1.968

1.744

1.783

2.05
Power Distribution (watts/cc)

The power peaking factors have been defined andanalyzed in terms of the following parameters:[a] Hot rod power peaking factor
fHR = (Prod.ave) max / Pcore.ave= hot rod factor determines the power
generation in the hottest fuel element

UZrH

U3O8‐AlU3Si2‐Al

16

22

1.85
1.65
1.45
1.25
1.05
0.85
0.65
0.45
0.25

4

10

Distance in Cm (radius)10III. Result and Discussion

Nuclear fuel life time and burn up
Fresh fuel1.06

Irradiated
fuel

1.05

UZrH (0.5 gU/cc)

SRAC-PIJ

Maximum Burn Up(K‐eff)Irradiation
simulation

K-eff
<1

U3Si2‐Al (3.26 gU/cc)

1.04

SRACCITATION

SWAT

U3O8‐Al (3.26 gU/cc)

1.03
1.02
1.01

Maximumburn up

1.00

Fuel Material

U‐ZrH1.65

U3Si2‐AlU3O8‐Al

Maximum burn up

18

54

68

0.99
0

15

30
45
Burn Up (GWD / ton of U)

60

75

11

IV. Summary

Neutronic parameter calculations, such as muliplication factor, power peaking factor, neutron flux,
nuclear life time and burn up have been performed for alternative fuel materials on TRIGA Kartini
reactor core:Especially for UO2 the reactor core is still in a subcritical situation on the all geometry changes.
For the excess reactivity U3O8‐Al > U3Si2‐Al > U‐ZrH1.65.
For total neutron flux U3Si2‐Al > U3O8‐Al > U‐ZrH1.65 . However, U‐ZrH1.65 can produce the highest
thermal neutron flux at all irradiation facilities.
All fuel materials produce power densities that are still in accordance with the peaking factor safety
standards.
The high level of uranium content makes the U3O8‐Al and U3Si2‐Al fuel life longer than U‐ZrH1.65.
Based on utilization and safety criteria, U3O8‐Al and U3Si2‐Al can be utilized as fuel materials
alternative but with different fuel density and different MTFR.

12

V. Future WorkNew core design of kartini reactor using alternative material fuel element plate type.
‐ Study the fuel type conversion, from TRIGA fuel rod type to the fuel plate type with the alternative fuel materials, U3O8‐
Al and U3Si2‐Al from the neutronic point of view, as a solution to the fuel shortage.

13

End of this presentation14

様々な核種の崩壊や生成などを考慮した
福島第一原子力発電所内部の放射能インベントリ解析
Radioactivity inventory analysis in the Fukushima Daiichi NPP
considering such as various nuclides generation and decay＊

亮憲 1，千葉

石井
1

豪1

北海道大学

汎用炉物理解析コード CBZ を用いて、福島第一原子力発電所内部に存在する放射能インベントリ解析を行った。この解析では様々な領域へ放射性物質が移行することを簡易的なモデルによって扱いつつ、自発核分裂による短半減期 FP の生成の影響なども詳細に想定し、定量的な評価を行った。
キーワード：廃炉、福島第一原子力発電所、自発核分裂、未臨界中性子増倍
1.

緒言：原子力発電所の廃炉作業では、作業員の被曝防止が重要である。本研究で取り扱う福島第一原子力発電所(以下 1F)は過酷事故を起こした発電所であり、健全な原子炉には存在しない箇所にも放射性物質が漏洩している。そこで現在では、SAMPSON や MAAP といった SA 解析コードによって、1F 内部の放射能イン
ベントリ解析が進められている。本研究の主目的は従来炉物理解析コードとして使用されていた CBZ を活用
し、発電所内部をいくつかの領域に簡易的に分割し、なおかつ実際の事故状況を踏まえ放射性物質の移行条件などを詳細に設定しインベントリ解析を行うことである。また、インベントリ解析を行う中で、自発核分
裂由来の中性子による未臨界核分裂連鎖反応によって生成された短半減期 FP の影響評価などにも取り組む。
2.

研究の内容：本解析では、事故以降に各領域中の放射能インベントリがそれぞれの領域間を移行することを想定している。 発電所内部を𝐽個の空間領域に分割した場合、ある時刻𝑡で各領域𝑗における核種𝑖の単位
体積当たりの数密度𝑁 , t が従う微分方程式を以下のように定義する。ここで𝑉は体積、𝜆は崩壊定数、𝐶は各領域間の移行係数、P は崩壊による核種の遷移確率を表す。本研究では、移行係数𝐶を元素ごとに与えて評価を行うが、この移行係数については、事故後報告されている様々な測定結果や研究結果を踏まえて設定する。

3.

計算結果：本解析では原子炉内部を燃料、ペデスタル、

ドライウェル、ウェットウェル、建屋、オペレーションフ
ロア、外部環境の 7 領域に分割して評価した。本解析はセシウムの移行挙動をエネ総研の解析結果や、大気中への放出の推定結果などに基づいて再現[1]し、その後に他の元素については、フランスで行われた過酷事故実験結果を引用した[2]。図 1 に各領域でのセシウムの存在割合の解析結果を示す。本結果から外部環境への放出と、またオペレーシ
ョンフロア内への漏洩も少ないことが確認できた。
4.

図 1 セシウムの初期インベントリ
に対する各領域での存在割合

結言：今後は、更に事故進展状況や、事故後の放射性物質の大気中への拡散のデータなどを統合して解析を行い、より詳細なインベントリ解析を行う。
参考文献：[1] 白井博嗣(エネルギー総合工学研究所)、日本原子力学会 2015 年春の大会予稿集 H26[2] 日本原子力学会「水化学」部会、Phébus FP プロジェクトにおける核分裂生成物挙動のまとめ* Akinori Ishii1, Go Chiba2
1Hokkaido Univ.

様々な核種の崩壊や⽣成などを考慮した
福島第⼀原⼦⼒発電所内部の放射能インベントリ解析
Radioactivity inventory analysis in the Fukushima Daiichi NPP
considering such as nuclides generation and decay

平成30年11⽉27⽇
Session Ⅳ:臨界安全・燃料デブリ
北海道⼤学⼤学院⼯学院
エネルギー環境システム専攻
原⼦炉⼯学研究室○⽯井亮憲 千葉豪⽯井 亮憲1⽯井5

研究背景• 福島第⼀原⼦⼒発電所の廃炉作業において、
放射能インベントリの解析が必要である。
• そこで、過酷事故解析コードを⽤いて事故後から、
原⼦炉内部の放射能インベントリの解析が各所で⾏われている。
• 事故直後から現在までに、放射性物質の漏洩のしかたなどについては、多くの調査結果が存在しており、これらも活⽤することによって、より信頼できるインベントリ解析が⾏えると考えられる。
• 本研究室で使⽤している、1400核種程度の解析が可能な
炉物理解析コードCBZを⽤いて、発電所内部の領域分けなどを設定しインベントリ評価が⾏えるようになった。

1

スライド 2
⽯井 亮憲1 ⼀番最後について、前々から1400核種程度のインベントリ解析はできているので、それをノーディングすることによって、発電所内部の解析ができるようになったというイメージ⽯井 亮憲, 2018/09/01

⽯井5

編集済み

⽯井 亮憲, 2018/09/02

2
CBZコードシステムで⽤いられている核種の崩壊と移⾏の計算モデル1
ある解析結果における放射性物質の移動の流れ
ウェットウェル

燃料

オペフロ
ドライウェル

ペデスタル

環境※⽮印は放射性物質の流れ

放射性物質がどのぐらいの時間にどの程度移動するかの設定が重要

3
CBZコードシステムで⽤いられている核種の崩壊と移⾏の計算モデル2
ある時刻t、あるノード𝑗における、原⼦核𝑖のインベントリの従う微分⽅程式

他の核種𝑖′が崩壊したことによる⽣成

半減期に従う崩壊

他領域𝑗′へ核種𝑖が移⾏することによる減少分

他領域𝑗′からの移⾏による核種𝑖の増加分

4
CBZコードシステムで⽤いられている核種の崩壊と移⾏の計算モデル3
ある時刻t、あるノード𝑗における、原⼦核𝑖のインベントリの従う微分⽅程式

この移⾏係数Cを、元素ごとに設定する。
このCはさまざまな実験結果や事故報告を踏まえて設定する

5
CBZ内で⾏われる⾃発核分裂ならびに、中性⼦による未臨界核分裂連鎖反応の計算1
本計算に使う⽅程式(⾃発核分裂を考えない場合)

この⽅程式が核種𝑖ごとにあるので、それをまとめたものをA⾏列とすると
↓↓

nは数密度

この𝒏を解いていくことによって、数密度(インベントリ)がわかる↓続いて、この計算に⾃発核分裂の影響を考慮する↓

6
CBZ内で⾏われる⾃発核分裂ならびに、中性⼦による未臨界核分裂連鎖反応の計算2
1回の核分裂によって⽣成される平均の中性⼦数

⾃発核分裂による⽣成
未臨界中性⼦増倍による⽣成
⾃発核分裂性核種が、⾃発核分裂を起こすことによって⽣じる平均の中性⼦数※今回は⾃発核分裂性核種としてCm-242, Cm-244を採⽤

上記2つの影響を𝑨⾏列に組み込み、計算を⾏う

研究⽬的

7

本研究のメインテーマ
福島第⼀原⼦⼒発電所内部の放射能インベントリの定量的な評価
また、この研究を進めるにあたり、以下のようなデータを活⽤する
1. 東電や学会などが提⽰している1Fの事故報告書
2. 各地で⾏われている過酷事故実験などの調査結果
3. 他の過酷事故解析コードで出⼒された結果
4. そのほか、事故後に1Fに関して報告されたデータ(例えば事故後のCs-137の⼤気中への拡散状況など)

さらに、このようにインベントリ解析を進めながら⾃発核分裂の影響などについても評価する。

⻑期間解析(3領域解析)-1
1. 燃料の損傷度合いについては、95%が損傷していると仮定する。
2. まず、エネ総研での解析結果に基づいて汚染⽔中のCs-137の移⾏を再現する。
3. Cs-137の挙動を再現するために設定したセシウムの移⾏係数を踏まえて、他の元素の移⾏係数を設定する。他の元素については、
Phébus原⼦炉のFPT3の実験結果をもとにして設定する。
4. FPT3の結果を踏まえて設定するが、FPT3で移⾏が確認された元素の中で⽔溶性のある物質のみが移⾏すると仮定した場合と、⾮⽔溶性のものもすべて移⾏すると仮定した場合の2種類で解析

89⻑期間解析(3領域解析)-2セシウム(Cs-137)の移⾏係数(移⾏半減期)設定の仕⽅• 事故発⽣~10⽇経過 ⇒ 5時間• 事故後のデブリ⼤量放出期間• 11⽇~30⽇⇒ 12時間• 3/23より冷却開始を受けて、それによる漏えい量の低下を考慮• 31⽇~400⽇⇒ 10年• 注⽔等が完了し、また燃料が⼗分冷却されたとみなす• 400⽇~1000⽇⇒ 1.5年• 燃料の劣化などによって新たに放射性物質が漏洩• 1001⽇以降

⾼
↑
移⾏
↓

⇒ 1.35年......

低• 劣化の進⾏

移⾏のしやすさのイメージ図

セシウム以外の元素については、過酷事故実験結果の放出量を踏まえて、セシウムの移⾏係数に⽐例する形で設定。

⻑期間解析(3領域解析)-3

時刻 t⻑期間解析のイメージ図

燃料

浄化系

Fuel

Purify system

汚染⽔
Water

解析箇所
移⾏の仕⽅にPhébusの実験結果を適⽤

10⻑期間解析(3領域解析)-411

汚染⽔中のFPT3(含⾮⽔溶性)の解析結果
事故後(⻑期)

汚染⽔中ではストロンチウムやセシウムが⽀配的な物質となる

⻑期間解析(3領域解析)-5

⾃発核分裂の影響評価1

12

セシウムの放射能Cs-134やCs-137の影響が強い

短半減期同位体
Cs-136の影響は弱い

事故から⻑期間が経過すると汚染⽔中では、⾃発核分裂による短半減期FPや、⾃発核分裂によって発⽣した中性⼦による核分裂反応の影響は少ない⻑期間解析(3領域解析)-6

⾃発核分裂の影響評価213

⾃発核分裂無し

短半減期核種の多いキセノン等は事故後⾃発核分裂等が発⽣すると放射能濃度が増加する。
↓↓

しかしながら、汚染⽔中の放射能濃度に与える影響は極めて⼩さい⻑期間解析(3領域解析)-7

3領域解析での結果

14• これまでの、過去の事故調査結果や放射性物質の移⾏の仕⽅を取り⼊れ、特に汚染⽔について解析をした。
• 今回は特に、エネ総研Cs-137の汚染⽔中の放射能濃度や、Phébus原⼦炉での実験結果を参考にして、解析した。
→今回の解析では、⾃発核分裂による短半減期FPが汚染⽔の放射能
インベントリに与える影響は、現時点では少ないと考えられる。

短期間解析(7領域解析)-1

解析体系の設定

15• エネ総研の発表内容を踏まえて7領域で解析
燃料、ペデスタル、ドライウェル、ウェットウェル、原⼦炉建屋、オペフロ、外部環境
CBZで新たに付け加えた箇所• 燃料からペデスタルとウェットウェルに漏れ出ると仮定して解析• 移⾏係数なども事故後想定されることなどを考えて
設定した。

燃
料

引⽤:原⼦炉過酷事故における放射性核分裂⽣成物挙動の評価
6.過酷事故解析コードによる福島事故時のFP放出・移⾏の検討(2015年春の年会 H26)

短期間解析(7領域解析)-2

解析シナリオ・モデルの設定

16• 事故後6時間程度

1. 溶融燃料落下開始
2. SR弁からウェットウェルに蒸気移動• 事故後11時間程度

格納容器から漏洩が発⽣(東電報告)• 事故後18時間程度

格納容器ベント開始(東電報告)• 事故後24時間程度

⽔素爆発(1号機)(※↑エネ総研のデータでは考慮されていないと推測)

引⽤:原⼦炉過酷事故における放射性核分裂⽣成物挙動の評価
6.過酷事故解析コードによる福島事故時のFP放出・移⾏の検討(2015年春の年会 H26)

短期間解析(7領域解析)-3

解析条件

17

• 今回は全燃料から移⾏が発⽣すると仮定• 左図のエネ総研の解析結果を基にして、セシウムの移⾏係数を設定。
• セシウムの移⾏係数に基づいて、他の元素についても移⾏係数を設定。これは⻑期間での解析と同様にPhébus
FPT3の試験結果を引⽤• エネ総研の解析結果では考慮されていない、⽔素爆発による外部環境放出の影響なども併せて考慮する。

引⽤:原⼦炉過酷事故における放射性核分裂⽣成物挙動の評価
6.過酷事故解析コードによる福島事故時のFP放出・移⾏の検討(2015年春の年会 H26)

短期間解析(7領域解析)-4

解析モデルの再現結果

18

CBZでの再現↓続いて、⽔素爆発の影響も考慮する↓

引⽤:原⼦炉過酷事故における放射性核分裂⽣成物挙動の評価
6.過酷事故解析コードによる福島事故時のFP放出・移⾏の検討(2015年春の年会 H26)

短期間解析(7領域解析)-519

⽔素爆発の影響を考慮(Cs-137)

初期量の2.2%放出

Cs-137における初期インベントリに対する各領域での存在割合

⽔素爆発の影響により、環境中の放射能の結果が可視化→環境中に数PBq程度放出

Cs-137の各領域における放射能

短期間解析(7領域解析)-6

⽔素爆発の影響を考慮(I-131)

初期量の1.8%放出

I-131における初期インベントリに対する各領域での存在割合

20

セシウムと同様に
⽔素爆発による放出量が可視化→環境中に数⼗PBq程度放出

I-131の各領域における放射能短期間解析(7領域解析)-7

他の核種(元素)の解析結果

キセノンは移⾏しないと仮定

キセノンにおける
初期インベントリに対する各領域での存在割合

21

ヨウ素の崩壊などによって環境中のキセノンの放射能も新たに解析→環境中に数PBq程度放出

キセノンの各領域における放射能

結論・今後の研究⽅針

22• ⻑期間解析の結果(再度解析を⾏った結果)→今回の解析では、⾃発核分裂による短半減期FPが汚染⽔の放射能
インベントリに与える影響は、現時点では少ないと考えられる。
• 短期間解析の結果→事故後⽔素爆発の影響によって放出される
放射性物質とエネ総研の解析結果と組み合わせることにより他元素の外部環境中の放射能の変化などを確認することができた。
今後は、事故直後から現時点までの⼀連の移⾏モデルを作成してインベントリ評価を⾏っていく。

Supercritical transient analysis using Multi-region Integral Kinetic code: Basics andapplications
Delgersaikhan Tuya and Toru Obara
Laboratory for Advanced Nuclear Energy, Institute of Innovative Research, Tokyo Institute of Technology
Supercritical transient analysis in arbitrary nuclear fissile system generally requires a model that describes a time- andspace-dependent power while taking related feedback phenomena in the considered system into account. Multi-region
Integral Kinetic (MIK) code is a Monte Carlo neutron transport method- and integral kinetic model-based numerical
tool for a supercritical transient analysis in an arbitrary fissile system. MIK code allows complicated and complex
feedback mechanisms in the system to be treated through its time-dependent feedback modeling capability provided
that the governing equations for feedback mechanisms are given.
Keywords：supercritical transient; criticality accident; Godiva; MIK code; integral kinetic model;
1. Introduction： MIK code [1] essentially consists of three parts: 1) an integral kinetic model (IKM), 2) Monte
Carlo neutron transport code MVP2.0/3.0, and 3) a time-dependent feedback model. The IKM [2] describes a timeand region-dependent fission rate in arbitrary fissile system using secondary fission probability density functions. The
Monte Carlo neutron transport code MVP2.0/3.0 [3][4] is utilized to
calculate the secondary fission probability density functions at an
instantaneous state of the system. The time-dependent feedback
model calculates time-dependent feedback variables (e.g.,
region-wise temperature and density) via user-provided governing
equations.
2. Methodology：MIK code has been applied to supercriticaltransient in Godiva core with step-wise reactivity insertion of about
8 cents above prompt-criticality as an example of application of
MIK code. Furthermore, the obtained simulation results have been
compared with experimental results. In the simulation, Godiva coreis divided into two regions of inner sphere and outer concentric
sphere. As for feedback, adiabatic temperature increase due toFig. 1. Total fission ratefission energy deposit and approximate thermal expansion are taken
into account.
3. Results： The obtained results in terms of total fission rate,
region-wise temperature, and region-wise density profiles areshown in Fig. 1, Fig. 2, and Fig. 3, respectively. Also shown in Fig.
1 is the total fission rate obtained from the experiment [5]. From
Fig. 1, it is seen that the simulation reproduced the total fission ratequite well near the peak, but diverged after the peak. The reason isthat the tail region of the experimental results was due to ashort-lived delayed neutron precursors, however, this effect was nottaken into account in the simulation. On the other hand, Fig. 2 andFig. 3 showcase the time-dependent feedback capability of MIK
code in case of time- and region-dependent temperature and density.
Fig. 2. Region-wise temperature
As heat loss from the system was neglected, the region-wisetemperature kept increasing leading to increasing thermal expansion(see Fig. 3), which in turn led to increasing negative feedback effectthat eventually stopped the supercritical transient in the system.
References: [1] D. Tuya and T. Obara. “Development of Monte Carloneutron transport method-based supercritical transient code with
time-dependent feedback capability”. Nucl. Sci. Eng. (submitted) (2018)[2] H. Takezawa and T. Obara. “New approach to space-dependent kinetic
analysis by the integral kinetic model”. Nucl. Sci. Eng. 171. p.1-12. (2012)[3] Y. Nagaya et al,. “MVP/GMVP II: General purpose Monte Carlo codes
for neutron and photon transport calculations based on continuous energy
and multigroup method”, JAERI-1348. (2005)[4] Y. Nagaya et al,. “MVP/GMVP Version 3: General purpose Monte
Carlo codes for neutron and photon transport calculations based oncontinuous energy and multigroup method”, JAEA-Data/Code 2016-018. (2016)[5] TF. Wimett et al,. “Time behavior of Godiva through prompt critical”, LA-2029. (1956)

Fig. 3. Region-wise density

7th Reactor Physics Special Meeting
Kyoto University
November 26 - 27, 2018

Supercritical transient analysis using
Multi-region Integral Kinetic code:
Basics and applications
Delgersaikhan TUYA and Toru OBARA

Contents
I.

Introduction

II.

Multi-region Integral Kinetic (MIK) code

III. Example application
IV. Conclusions
V.

Future study2

Contents
I.

Introduction

II.

Multi-region Integral Kinetic (MIK) code

III. Some applications
IV. Conclusions
V.

Future study

3

Introduction• Supercritical transient analysis aims to quantify a power profile, energy release, and dose
release etc in a supercritical excursion.
• Supercritical transient analysis in a given system is often complicated due to a presenceof complex and nonlinear feedback mechanisms, nonuniform distribution of material andother variables/parameters, and complicated geometry.
• Generally, conventional deterministic methods such as point-reactor kinetics model andmulti-point kinetics model have been used.
• Recently, Monte Carlo neutron transport based methods1 for transient/dynamic analysis
are getting more practical.
• Multi-region Integral Kinetic (MIK) code is a Monte Carlo neutron transport based toolthat can be used to calculate a power profile and energy release in a supercriticaltransient excursion, taking the relevant feedback mechanisms into account via userdefined governing equations.
1 Sjenitzer, Bart L., and J. Eduard Hoogenboom. "Dynamic Monte Carlo method for nuclear reactor kinetics calculations." Nuclear Scienceand Engineering 175.1 (2013): 94‐107.
Mylonakis, A. G., Varvayanni, M., Grigoriadis, D. G. E., & Catsaros, N. (2017). “Developing and investigating a pure Monte‐Carlo module fortransient neutron transport analysis”. Annals of Nuclear Energy, 104, 103‐112.
J. Leppanen, “Development of a dynamic simulation mode in SERPENT 2 Monte Carlo code,” in M&C 2013, pp. 117–127, Sun Valley, Idaho(2013).

4

Contents
I.

Introduction

II.

Multi-region Integral Kinetic (MIK) code

III. Example application
IV. Conclusions
V.

Future study

5

Multi-region Integral Kinetic code - Introduction• MIK code has been developed in Obara lab., Tokyo Institute of Technology.
• MIK code essentially calculates a time- and region-dependent fission rate(i.e., power) in arbitrary fissile system.
• MIK code allows a treatment of time- and region-dependent feedback via
Ruser-defined governing equations for update and non-update variables.
• MIK code essentially consists of 3 parts:• Integral Kinetic Model (IKM)• Monte Carlo neutron transport based method (based on MVP2.0/MVP3.0 codes)• Time-dependent feedback model• MIK code can run in parallel mode using Open MPI 3.0 or higher.
• MIK code is an open source code (as of now, MIK code, Version 1.0).
• Request to use MIK code can be made to Prof. Toru OBARA of Tokyo Tech.

6

Multi-region Integral Kinetic code - Integral kinetic model (IKM)• The Integral Kinetic Model (IKM) for coupled system of R fissile region is:
Here, 𝑵𝒊 𝒕 is total fission rate in region i at time t (fissions/ s); 𝜶𝒊𝒋 𝒕 𝒕′ is secondary prompt‐fission probability
density function in region i provided by the first fission in
the source region j with time difference t‐t’ (fissions in i∙s‐1/ source fissions in j )• Discretized form (with initial critical condition):𝜏

where, 𝐶𝑖𝑗 𝜏

0𝛼𝑖𝑗 𝜏′ 𝑑𝜏′• The key kinetic functions 𝛼

 Cumulative distribution function

𝜏′ or 𝐶𝜏 :• How to obtain?

7

Multi-region Integral Kinetic code – Monte Carlo neutron transport based
methodi
S2
F

Arbitrary material containing fissile
material(s)

S1
B

j

B – point in region j, at which a neutron wasborn at τ=0
S1 and S2 – points at which neutron was
scatteredF – point in region i, at which neutron induced
a fission τ= τ’

MVP2.0/3.0 code can be used to calculate a neutron transport time‐ and region‐
dependent fission reaction per time interval/bin as:

8

Multi-region Integral Kinetic code – Monte Carlo neutron transport based
methodCumulative fission by time k∆t

Cumulative fission by time entire transport time

Total fission in source region at end of

Total fission in source region at

generation

beginning of generation

Cumulative secondary fission distribution function by time k∆t

9Multi-region Integral Kinetic code – Time-dependent feedback model• In MIK code, a state of a system is partially characterized by set of variables related
to feedback - they are called update variables in the code
• Example of update variables: temperature, density etc• These update variables change during a supercritical transient• Accordingly, kinetic function 𝐶
Fission energy
deposition in thesystem𝜏 changes too!

Change in update
variables

Change in kinetic
functions

Change indynamics of thetransient• Need to recalculate/update the kinetic functions at the instantaneous state of thesystem• How often/when do we need to update the kinetic functions?

10Multi-region Integral Kinetic code – Time-dependent feedback model• The criteria to update the kinetic function is:U

𝑣‐th update variable in region 𝑟𝑘current time step𝑘time step at which the last update wasdoneδupdate criterion for 𝑣‐th update variable
in region 𝑟• Update when the change in one of update
variables in any region exceeds the userdefined criteria!
• User can control the trade-off between
computation cost and calculation accuracy
via update criteria

11

Contents
I.

Introduction

II.

Multi-region Integral Kinetic (MIK) code

III. Example application
IV. Conclusions
V.

Future study

12

Example application – Godiva supercritical transient• A series of supercritical critical experiments1 had been performed at Godiva facility at Los
Alamos Scientific Laboratory in 1950s• MIK code was used to simulate supercritical transient with 8.12 ¢ step-wise reactivity
insertion above prompt criticality (delayed neutron effect was ignored in the simulation).
• In Godiva core, the dominant feedback is thermal expansion (although in high reactivity
insertion case, there is a delay of thermal expansion)• The thermal expansion feedback without the delay was considered in the MIK simulation

Godiva core in its
disassembled state

Godiva core as a two‐region system in
MIK simulation (XY cross‐sectional view)

1 T. F. WIMETT et al., “Time behavior of Godiva through prompt critical, Report

No. LA‐2029” (1956).

13

Example application – Godiva supercritical transient• Some basic parameters:

Enrichment of 235U, wt%

93.7

Outer radius, cm

8.85

𝜀, MeV/f

178• Time-dependent feedback:• Update variables: temperature (K), density (g/cm3), radius (cm)• Non-update variable: energy release (J)• Update criteria matrix:
δ

2
2

0.05
0.05

0.005
0.005

• User-defined governing equations:Temperature
J/K/mol
Volume expansion1

1 D. H. GURINSKY and G. J. DIENES, Nuclear Fuels. The Geneva Series on the Peaceful Uses ofAtomic Energy,D. Van Nostrand Inc. (1956).
14

Example application – Godiva supercritical transient• Obtained results:

As delayed neutrons wereignored, the tail region of
overall fission rate profile
was not simulated!

Experiment

Simulation

Peak fission rate,f/s

2.9×1020

3.2×1020 f/s(1 %)Fission release, f

1.3×1016

1.2×1016(1×10‐3 %)15

Example application – Some other applications
For some other applications of MIK code including TRACY and SILENE, please take look at:
Delgersaikhan Tuya and Toru Obara, “Development of Monte Carlo neutron transport method‐based
supercritical transient code with time‐dependent feedback capability”, Nuclear Science and Engineering(2018)
DOI: https://doi.org/10.1080/00295639.2018.1540209

16

Contents
I.

Introduction

II.

Multi-region Integral Kinetic (MIK) code

III. Example application
IV. Conclusions
V.

Future study

17

Conclusions• Multi-region Integral Kinetic (MIK) code has been developed• MIK code can calculate time- and region-dependent power in arbitrary fissile systemfollowing a reactivity insertion above prompt criticality• MIK code allows a treatment of time-dependent feedback mechanisms via user-defined
governing equations• MIK code can run in single or parallel mode on shared-memory system using Open MPI
3.0 or higher• MIK code is an open source code and can be requested from Prof. Toru Obara of Tokyo
Institute of Technology• MIK code was applied to Godiva supercritical transient and the obtained results were in
reasonably good agreement with the experimental results• MIK code was also applied to supercritical transient in TRACY and SILENE (not reported
here)
18

Contents
I.

Introduction

II.

Multi-region Integral Kinetic (MIK) code

III. Example application
IV. Conclusions
V.

Future study

19

Future study• Currently a work is underway to include the effect of delayed neutron in MIK code
• A work is also underway to investigate the possibility of coupling MIK code with an open
source thermal hydraulics (CFD / thermo-mechanics) code to improve feedback modelingcapability• As MIK code is an open source code, if you have interest and need, please request anduse it; modify/improve it; or find bugs and criticize us (not in a mean way!)20

ブートストラップ法を活用した Feynman-𝛂法の測定共分散評価
Estimation of experimental covariance in the Feynman-α method using the bootstrap method＊1

遠藤 知弘 1

名古屋大学

バンチング法による𝑌 𝑇 の共分散を推定できるようブートストラップ法を利用した Feynman-α 法を改良した。
キーワード：Feynman-α 法、バンチング法、共分散、相関、ブートストラップ法
1.

緒言：実験の測定時間は限られるため、有限の測定データに基づいて測定結果の不確かさ評価を行うことが重要となる。炉雑音解析手法の一種である Feynman-α 法の場合、中性子検出時刻情報が記録された一連の時系列データに対してバンチング法を利用することで、任意の検出時間幅𝑇に対する二次相関量𝑌 𝑇 を求め、解析式のフィッティングにより即発中性子減衰定数𝛼の推定を行う。バンチング法では同じデータを再利用するため、異なる𝑇間の𝑌 𝑇 に相関(共分散𝚺𝒀 )が生じる。そこで本研究では𝚺𝒀 の推定法を新たに考案した。
2.

改良ブートストラップ法：図 1 で示したようにリサンプリング法を改良した[1]。元の測定結果から再標本化する際に、𝑁

個の連続した時系列データを並べることでバンチング法による共分散𝚺𝒀 を推定する。

(a) 元の中性子検出数時系列の測定データ(b) 改良ブートストラップ法
共分散行列の推定

図 1 𝑌 𝑇 共分散推定のためのブートストラップ法の適用例(𝑁
3.

5の場合)

実験結果：停止中 KUCA 未臨界炉心において 10 分間炉雑音測定を 93 回実施し、𝑌 𝑇 の共分散の参照値を求めた。また、93 回の測定結果のうち 1 回の 10 分間炉雑音測定を選び、提案手法により共分散𝚺𝒀 を推定した。図 2 より、提案手法により 1 回の炉雑音測定から共分散𝚺𝒀 の概算が可能であることを確認できた。

図 2 バンチング法に起因した𝑌 𝑇 の測定共分散行列𝚺𝒀
参考文献：[1] T. Endo, A. Yamamoto, Ann. Nucl. Energy. 124, pp.606–615, DOI: 10.1016/j.anucene.2018.10.032 (2019).
謝辞：本研究は JSPS 科研費(15K18317, 17K14909)の助成を受けたものです。
* Tomohiro Endo1
1 Nagoya Univ.

ブートストラップ法を活用した
Feynman-α法の測定共分散評価
名古屋大学
 遠藤 知弘
要約：バンチング法による2次中性子相関量𝑌 𝑇 の共分散を推定できるようブートストラップ法を利用したFeynman-α法を改良した
謝辞：本研究はJSPS科研費(15K18317,17K14909)の助成による第7回炉物理専門研究会, 京都大学 複合原子力科学研究所, 2018/11/26–27

1

本日の発表内容
 炉雑音解析手法
 Feynman-α法

 バンチング法に起因した測定共分散 ブートストラップ法

 KUCA実験による検証
 参照値(複数回測定の共分散)との比較
 共分散が𝛼フィッティング結果に与える影響
第7回炉物理専門研究会, 京都大学 複合原子力科学研究所, 2018/11/26–27

2

背景：未臨界実験の重要性
1. 未臨界監視技術の開発
exp

 加速器駆動システム
 1F燃料デブリ取り出し時

𝛼𝑡

2. 未臨界実験結果を活用したデータ同化
 停止状態の実験施設の有効活用→核データ更新
 測定しやすい量として即発中性子減衰定数𝛼に注目
 臨界近傍： 𝛼
 パルス中性子源を用いたパルス中性子法
 だけでなく炉雑音測定でも測定可能
3

第7回炉物理専門研究会, 京都大学 複合原子力科学研究所, 2018/11/26–27

Feynman-α法(1/2)
 検出時間幅𝑇の中性子検出数𝐶 𝑇 の時系列データ
𝐶

𝐶

𝐶

𝐶

𝐶

中性子
検出
𝐶

𝐶

𝐶

𝐶

𝐶

𝐶

t

𝐶

𝑇
 𝐶 𝑇 の平均𝜇 と分散𝜎 の推定

𝜇 𝑇 ≡ 𝐶 𝑇

𝜎 𝑇 ≡

𝐶 𝑇

𝜇 𝑇

第7回炉物理専門研究会, 京都大学 複合原子力科学研究所, 2018/11/26–27

4Feynman-α法 (2/2)
𝜎 𝑇
𝑌 𝑇 ≡
𝜇 𝑇

1
Y(T) [-]

 𝑌 𝑇 測定結果に対して解析式をフィッティング→即発中性子減衰定数
𝛼を推定

𝛼

𝜌$

1

𝛼

T [sec]𝛼crit𝛽 /Λ: 臨界時の𝛼5

第7回炉物理専門研究会, 京都大学 複合原子力科学研究所, 2018/11/26–27

Feynman-α法における測定共分散 検出時間幅𝑇を変えた測定のやり直し
 𝑇を変えるたびに長い測定時間を要する

𝑇
𝐶

𝐶

𝐶

𝐶

𝐶

𝐶

𝐶

𝐶

𝐶

𝐶

𝐶

𝐶

𝐶

𝐶

𝐶

𝐶

𝐶

𝐶

𝐶

𝐶

𝐶

𝐶

𝐶

𝐶

2𝑇
𝐶

𝐶

𝐶

𝐶

𝐶

𝐶

𝐶

𝐶

𝐶

𝐶

𝐶

𝐶

3𝑇
𝐶

𝐶

𝐶

𝐶

𝐶

𝐶

𝐶

第7回炉物理専門研究会, 京都大学 複合原子力科学研究所, 2018/11/26–27

𝐶
6

Feynman-α法における測定共分散 検出時間幅𝑇を変えた測定のやり直し
 𝑇を変えるたびに長い測定時間を要する

 バンチング法
 1つの時系列データを再活用して𝑌 𝑇 の変化を推定

𝑇
𝐶

𝐶

𝐶

𝐶

𝐶 𝐶

𝐶

𝐶

𝐶 𝐶

𝐶

𝐶

𝐶

𝐶 𝐶

𝐶

𝐶

𝐶 𝐶

𝐶

𝐶

𝐶
𝐶
𝐶

𝐶

𝐶
𝐶
𝐶

𝐶

𝐶
𝐶
𝐶

𝐶

𝐶

𝐶

𝐶
𝐶

𝐶

𝐶

𝐶
𝐶

𝐶

𝐶

𝐶
𝐶

𝐶

𝐶
𝐶

同一の時系列データに基づいて𝑌 𝑘𝑇 , 𝑌 𝑘 𝑇 を推定するため異なる時間幅𝑇間の𝑌値に相関(共分散)が生じざるを得ない
7

第7回炉物理専門研究会, 京都大学 複合原子力科学研究所, 2018/11/26–27

課題：Feynman-α法の測定共分散推定
時系列データ1

𝐶,

𝐶,

𝐶,

𝐶,

𝐶,

𝐶,

𝑌 𝑇

𝐶,

𝐶,

・・
・

𝐶,

∑

𝑌

𝑌
𝑀

時系列データ2

𝐶,

𝑌 𝑇 共分散行列𝚺𝒀

𝐶,

𝐶,

𝑌

𝑌

1

𝑌 𝑇

 多数回測定によるバンチング法起因の共分散推定
 長い測定時間
 各時系列データ長さ&回数𝑀の設定が悩ましい

目的：1回の時系列データから共分散
𝚺𝒀 推定手法の開発
第7回炉物理専門研究会, 京都大学 複合原子力科学研究所,
2018/11/26–27

8

ブロック・ブートストラップ法(a) 元の中性子計数・時系列データ𝐶 𝑇
𝐶

𝐶

𝐶

𝐶

𝐶

𝐶

𝐶

𝐶

𝐶

𝐶

𝐶

𝐶

𝐶

𝐶

𝐶

𝐶

𝐶

𝐶

𝐶

𝐶

𝐶

𝐶

𝐶

𝐶

(b) ブートストラップ法による再標本化
𝐶

𝐶

𝐶

𝐶

𝐶

𝐶

𝐶

𝐶

𝐶

𝐶

𝐶

𝐶

𝐶

𝐶

𝐶

𝐶

𝐶

𝐶

𝐶

𝐶

𝐶

𝐶

𝐶

𝐶

𝐶

𝐶

𝐶

𝐶

𝐶

𝐶

𝐶

𝐶

𝐶

𝐶

𝐶

𝐶

𝐶

𝐶

𝐶

𝐶

𝐶

𝐶

𝐶

𝐶

𝐶

𝐶

𝐶

𝐶

𝐶

𝐶

𝐶

𝐶

𝐶

𝐶

𝐶

𝐶

𝐶

𝐶

𝐶

𝐶

𝐶

𝐶

𝐶

𝐶

𝐶

𝐶

𝐶

𝐶

𝐶

𝐶

𝐶

𝐶

𝐶

𝐶

𝐶

𝐶

𝐶

𝐶

𝐶

𝐶

𝐶

𝐶

𝐶

𝐶

𝐶

𝐶

𝐶

𝐶

𝐶

𝐶

𝐶

𝐶

𝐶

𝐶

𝐶

𝐶

𝐶

𝐶

𝐶

𝐶

𝐶

𝐶

𝐶

𝐶

𝐶

𝐶

𝐶

𝐶

𝐶

𝐶

𝐶

𝐶

𝐶

𝐶

𝐶

𝐶

𝐶

𝐶

𝐶

𝐶

9

第7回炉物理専門研究会, 京都大学 複合原子力科学研究所, 2018/11/26–27

再帰バンチング法
𝐶

𝐶

𝐶

𝐶

𝐶

𝐶

𝐶

𝐶

𝐶

𝐶

𝐶

𝐶

𝐶

𝐶

𝐶

𝐶

𝐶

𝐶

𝐶

𝐶

𝐶

𝐶

𝐶

𝐶

𝐶

𝐶

𝐶

𝐶

𝐶

𝐶

𝐶

𝐶

𝐶

𝐶

𝐶

𝐶

𝐶

𝐶

𝐶

𝐶

𝐶

𝐶

𝐶

𝐶 𝐶

𝐶

𝐶

𝐶

𝐶
𝐶

𝐶
𝐶

𝐶
𝐶

𝐶
𝐶

𝐶
𝐶

𝐶
𝐶

𝐶
𝐶

𝐶
𝐶

𝐶
𝐶

𝐶
𝐶

𝐶
𝐶

𝐶
𝐶

𝐶

𝐶
𝐶

𝐶

𝐶
𝐶

𝐶

𝐶
𝐶

𝐶
𝐶

𝐶
𝐶

𝐶
𝐶

𝐶

𝐶
𝐶

𝐶
𝐶

𝐶
𝐶

𝐶
𝐶

𝐶
𝐶

𝐶
𝐶

𝐶
𝐶

𝐶
𝐶

𝐶
𝐶

𝐶
𝐶

𝐶
𝐶

𝐶
𝐶

ある時系列データを起点として再帰的に隣接する2個の和をとることで加算回数を削減→処理時間の短縮
10
第7回炉物理専門研究会, 京都大学 複合原子力科学研究所, 2018/11/26–27

KUCAでの炉雑音測定実験
 運転停止状態

全制御棒 全挿入 中心架台落下 外部中性子源無し


p

p

p

p

p

p

p

p

p

p

p

p

4

p

p

p

燃料内の固有中性子源のみ、 p
主に238U 自発核分裂p

 ゆえに低計数率

p

p

p

p

p

p

1

p

p

p

p

p

p

p

p

p

p

p

p

p

p

p

p

p

p

p

p

p

p

2

p

p

p

p

p

p

p

p

F

F

F

F

F

F

F

F

F

F

F

F

F

F

F

F

F

F

F

F

F

F

F

F

F

p

p

p

p

p

p

F

F

F

p

p

p

p

p

p

3

p

p

p

p

p

p

p

p

enriched U plateA3/8” p36EU-NU core

#2位置のHe3検出器で4.444 0.011 [count/sec] Ave.235U wt% ： 5.4wt%
H/235U ratio

unit cell



p

： 270

natural U plate11

第7回炉物理専門研究会, 京都大学 複合原子力科学研究所, 2018/11/26–27

測定共分散の参照値
 (10分間の炉雑音測定)

×93回

𝚺𝒀

∑

𝑌

𝑌
93

𝑌

𝑌

1

相関係数
12

参照解

ブートストラップ法による測定共分散の推定
 10分間の炉雑音測定1回のみ


リサンプリング数𝐵1000

the 47th trial

13の共分散が与える影響
 最小二乗法によるフィッテング
 「各点が独立な正規分布に従う」と仮定する場合が多い カイ二乗値𝜒

∑を最小化

 誤差𝜎 の相対値は推定できるが絶対値が不明な場合
 「𝜒

自由度」となるように誤差𝜎 の絶対値を規格化
してフィッティング誤差を推定する場合あり 例) scipy.optimize.curve_fitの初期設定：
absolute_sigma=False

第7回炉物理専門研究会, 京都大学 複合原子力科学研究所, 2018/11/26–27

14

𝑌 𝑇 の共分散を無視、相対誤差を重みとした𝛼フィッティング結果
3500±1.96σ内の被覆確率39%

[1/s]

3000
2500
2000
1500
1000
500
0

20

40
60
trial number

80

100

𝛼のフィッティング誤差を過少評価
第7回炉物理専門研究会, 京都大学 複合原子力科学研究所, 2018/11/26–27

15

𝑌 𝑇 の共分散を無視、絶対誤差を与えた
𝛼フィッティング結果
3500±1.96σ内の被覆確率80%

[1/s]

3000
2500
2000
1500
1000
500
0

20

40
60
trial number

80

100

被覆確率が改善→共分散の影響がある場合
𝜒 に基づき規格化すべきでない
第7回炉物理専門研究会, 京都大学 複合原子力科学研究所, 2018/11/26–27

16

𝑌 𝑇 の共分散を考慮、絶対誤差を与えた
𝛼フィッティング結果
3500±1.96σ内の被覆確率89%

[1/s]

3000
2500
2000
1500
1000
500
0

20

40
60
trial number

80

100

共分散𝚺𝒀∗ を考慮することで被覆確率は95%により近くなる
第7回炉物理専門研究会, 京都大学 複合原子力科学研究所, 2018/11/26–27

17

ブートストラップ法による 頻度分布
 リサンプリングした各時系列データで𝑌 ∗𝑇 だけでなく即発中性子減衰定数𝛼 ∗ も併せて推定(共分散考慮)

第7回炉物理専門研究会, 京都大学 複合原子力科学研究所, 2018/11/26–27

18

ブートストラップ法による𝛼信頼区間
推定結果(共分散を考慮)
3500

[1/s]

3000
2500
2000
1500
1000
95%bootstrap信頼区間内の被覆確率94%

500
0

20

40
60
trial number

80

100

正規分布を仮定せずに統計誤差推定が可能
第7回炉物理専門研究会, 京都大学 複合原子力科学研究所, 2018/11/26–27

19まとめ バンチング法に起因した𝑌 𝑇 測定共分散推定方法
 ブートストラップ法により概算可能
 参照解(多数回測定の結果)とも比較的良く一致

 共分散𝚺𝒀 は即発中性子減衰定数𝛼の誤差推定に影響 共分散𝚺𝒀 を考慮することで妥当な誤差推定が可能

 今後の課題
 バンチング法に起因した共分散𝚺𝒀 の理論式導出∵ブートストラップ法には多くの処理時間が必要→実時間での推定には不向き
第7回炉物理専門研究会, 京都大学 複合原子力科学研究所, 2018/11/26–27

20粒子フィルタ法を用いた未臨界度と中性子生成時間の推定
Estimation of subcriticality and neutron generation time using the particle filter method＊

池田 卓弥 1, 遠藤 知弘 1，山本 章夫 1
1

名古屋大学

未臨界度および一点炉動特性パラメータの同時推定が可能な手法として、粒子フィルタ法に注目している。
本手法の適用可能性検討の第一段階として、未臨界度と中性子生成時間が同時変化すると仮定した場合について、京都大学臨界集合体(KUCA)の実験で得られた中性子計数の時系列データから未臨界度を推定した。
キーワード：未臨界度、一点炉動特性方程式、粒子フィルタ法、京都大学臨界集合体(KUCA)

1. 緒言

福島第一原子力発電所のような過酷事故後の炉心を対象に、従来の逆動特性法による未臨界度推定を行う場合、一点炉動特性パラメータ(Λ, 𝛽 )の不確かさが大きいため、その不確かさが伝播することにより、未臨界度𝜌 の推定結果にも大きな不確かさが生じるという問題がある。その問題を解決する手法

として本研究では粒子フィルタ法[1]に着目した。粒子フィルタ法では、パラメータ 例: 𝜌, Λ, 𝛽が従う確率分布を多数の”粒子”で近似する。個々の粒子における中性子数の時間変化を支配方程式(例:一点炉動特性方程式)に基づいて予測し、観測値(例:中性子計数)との差異が小さい(尤度の値が大きい)粒子をベイズ推定
により優先的に残すことで、パラメータの予測結果を動的に改善する。粒子フィルタ法に着目した理由としては、行列演算を必要とせず実装が容易であり、数値モデルの非線形性が強い場合でも適用可能なことが挙げられる。

𝜌, Λ, 𝛽の同時推定を行おうとする場合、数値モデルの非線形性が強くなることが考えられるが、粒子フィルタ法によって推定が行える可能性がある。本発表では、未臨界度と中性子生成時間が同時に変化すると仮定した場合について、粒子フィルタ法により未臨界度と中性子生成時間の逆推定を実施した結果を報告する。

0.00

10/30-11/2 に KUCA-A 架台において DT 中性子源を用いた未臨界度
過渡変化実験を実施した。制御棒(C1,C2,C3)操作および中心架台落下による中性子計数の時間変化を測定し、粒子フィルタ法を適用した。粒子フィルタ法
では遅発中性子 6 群一点炉動特性方程式に基づき中性子計数を予測し、未臨界度

Reactivity[dk/k]

2. 実験解析

C1挿⼊

-0.05
PF＋2σ

-0.10

PF-2σ
-0.15

Reference-0.20
0と中性子生成時間の逆推定を行った。
MCNP の計算結果を参照値として、粒子C3落下 C2落下中⼼架台落下

500

1000

1500
2000
Time[s]

2500

3000

図 1 粒子フィルタ法による未臨界度の推定結果

フィルタ法による推定結果と比較する
ことで、粒子フィルタ法による推定の妥当性を確認する。
3. 推定結果

未臨界度推定を行った結果を図 1 に示す。ここで、𝜎は 10000 個の粒子から得られた推定結果の標準偏差である。図 1 から、粒子フィルタ法による推定結果 2𝜎の範囲内に参照値が含まれていることが確認できる。中性子生成時間の推定結果については当日の発表において報告する。
参考文献
謝辞[1] M.S. Arulampalam, et al., IEEE Trans. Signal Processing, 50[2], pp.174-188, (2002).

本研究は KUCA 共同利用研究(CA3004)の成果であり、JSPS 科研費(17K14909)の助成による。

* Takuya Ikeda1, Tomohiro Endo1, Akio Yamamoto1
1Nagoya Univ.

KUCA実験データを利⽤した粒⼦フィルタ法による未臨界度と推定○池⽥卓弥、遠藤知弘、⼭本章夫
名古屋⼤学
第7回炉物理専⾨研究会, 京都⼤学複合原⼦⼒研究施設, 2018年11⽉27⽇

謝辞
本研究はKUCA共同利⽤研究(CA3004)の成果である
本研究はJSPS科研費(17K14909)の助成による内容
背景と⽬的粒⼦フィルタ法
未臨界度と中性⼦⽣成時間の相関
KUCA実験解析条件
粒⼦フィルタ法による推定結果
まとめ＆今後の課題

2

背景 燃料デブリ取出し3

福島第⼀原⼦⼒発電所(1FNPS)の燃料デブリ取出し
現在の推定: 炉内は未臨界状態
未臨界状態の燃料デブリ取出しにより、H/U⽐の変化が⽣じる可能性あり

H/U⽐の変化が予期せぬ正の反応度添加を引き起こす可能性あり予期せぬ放射線被ばくを防ぐために、臨界近接の
リアルタイム監視が必要[1]IRID, 2016, http://irid.or.jp/_pdf/201509to10_10.pdf?v=3

背景 未臨界度推定における問題
臨界近接検知⼿法の⼀つ
従来の逆動特性法
𝜌 𝑡

𝛽

Λ
𝑛 𝑡

𝑑𝑛 𝑡
𝑑𝑡

𝜆 𝐶 𝑡

𝑆

リアルタイムでの未臨界度監視が可能

問題点
適切な⼀点炉動特性パラメータが必要
1FNPS燃料デブリ取出しに逆動特性法の利⽤を考える
 H/U⽐の変化により未臨界度と⼀点炉動特性パラメータの同時変化が考えられる

適切な⼀点炉動特性パラメータを⽤いなければ適切な未臨界度が推定できない

4

⽬的

5

本研究の⽬的 適切な⼀点炉動特性パラメータが分からない状態で未臨界度が推定できる⼿法を検討すること“粒⼦フィルタ法”に注⽬

発表内容

t6
 ロボット⼯学、経済学等で利⽤される⼿法

 粒⼦フィルタ法による未臨界度推定の妥当性検討
 未臨界度と中性⼦⽣成時間が同時変化すると仮定してKUCA実験データを利⽤し実験解析を実施
 推定結果と参照値とを⽐較

スライド 5t6

具体的な分野追加しました
t-ikeda, 2018/11/23

粒⼦フィルタ法

6利⽤するモデル:⼀点炉動特性⽅程式
𝑑𝑛 𝑡
𝑑𝑡

𝛽
𝜌 𝑡
Λ 𝑡

𝑛 𝑡

𝜆 𝐶 𝑡

𝑆

𝛽
𝑑𝐶 𝑡
𝑛 𝑡
𝜆 𝐶 𝑡 𝑘 1,2, ⋯ , 6
Λ 𝑡𝑑𝑡
可観測な中性⼦数𝑛に対して 𝜌, Λ, 𝛽 が⾮線形な関係
粒⼦フィルタ法[2]⾮線形性の強いモデルでも利⽤可能パラメータの確率分布を複数の“粒⼦”で近似
粒⼦︓中性⼦数𝑛や反応度𝜌など推定したいパラメータを成分としてもつベクトル𝑛 𝐶 ⋯ 𝐶 𝜌 Λ𝑥[2]M.S. Arulampalam, et al., IEEE Trans. Signal Processing, 50[2], pp.174-188, (2002).

粒⼦フィルタ法
2.7E+03

2.7E+03

中性⼦計数率[count/0.5 s]

中性⼦計数率[count/0.5 s]

中性⼦数

2.8E+03予測
𝑡 2.8E+03
1

2.9E+03

2.6E+03
2.5E+03
2.4E+03
0

50

頻度

0

予測

25
頻度

𝑡

2.6E+03
2.5E+03

2.4E+03100150 0

反応度[Δk/k]

反応度[Δk/k]

反応度1.6E-03
1.5E-03
1.4E-03
1.3E-03
1.2E-03
1.1E-03
9.7E-04
8.7E-04
7.7E-04
6.7E-04

2.9E+03

予測

1.6E-03
1.5E-03
1.4E-03
1.3E-03
1.2E-03
1.1E-03
9.7E-04
8.7E-04
7.7E-04
6.7E-04
50 0

50

頻度

25
頻度

100

150

50

 予測したいパラメータ(反応度や中性⼦⽣成時間)に”システムノイズ”を与える
 システムノイズ︓予測の不確かさ
 ⼀点炉動特性⽅程式に基づき中性⼦数を予測

7

粒⼦フィルタ法
2.7E+03

2.7E+03

2.6E+03
2.5E+03
2.4E+03
0

50

頻度

0

2.6E+03
2.5E+03

2.4E+03100150 0

反応度[Δk/k]

反応度[Δk/k]

反応度1.6E-03
1.5E-03
1.4E-03
1.3E-03
1.2E-03
1.1E-03
9.7E-04
8.7E-04
7.7E-04
6.7E-04

2.9E+0325
頻度

中性⼦計数率[count/0.5 s]

中性⼦計数率[count/0.5 s]

中性⼦計数率[count/0.5 s]

中性⼦数

2.8E+03予測
𝑡 2.8E+03
1

2.9E+0350

1.6E-03
1.5E-03
1.4E-03
1.3E-03
1.2E-03
1.1E-03
9.7E-04
8.7E-04
7.7E-04
6.7E-04
50 0

頻度

尤度計算
2.9E+03

𝑡

8

尤度計算

2.8E+03
2.7E+03
2.6E+03
2.5E+03
2.4E+03
0.E+00
150

100

25
頻度

5.E-03
重み1.E-0250

尤度の計算 尤度関数には正規分布を仮定
 観測値に近い粒⼦ほど尤度が⼤きくなる

2.7E+03

2.7E+03

2.5E+03
2.4E+03
0

50

頻度

0

25
頻度

2.5E+03

2.4E+03100150 0

反応度[Δk/k]

反応度[Δk/k]

反応度1.6E-03
1.5E-03
1.4E-03
1.3E-03
1.2E-03
1.1E-03
9.7E-04
8.7E-04
7.7E-04
6.7E-04

2.6E+03

1.6E-03
1.5E-03
1.4E-03
1.3E-03
1.2E-03
1.1E-03
9.7E-04
8.7E-04
7.7E-04
6.7E-04
50 0

リサンプリング50

頻度

2.9E+03𝑡尤度計算

2.8E+03
2.7E+03
2.6E+03
2.5E+03
2.4E+03
0.E+00150100

5.E-03重み反応度[Δk/k]

2.6E+03

2.9E+03

中性⼦計数率[count/0.5 s]

中性⼦計数率[count/0.5 s]

中性⼦計数率[count/0.5 s]

中性⼦数

2.8E+03予測
𝑡 2.8E+03
1

2.9E+03

リサンプリング 9
中性⼦計数率[count/0.5 s]

粒⼦フィルタ法25
頻度50

 尤度の⼤きさを確率とみなして復元抽出 尤度の⼤きいものほど多数復元
 粒⼦が0に近いものは消滅

2.9E+03

リサンプリング

2.8E+03
2.7E+03
2.6E+03
2.5E+03
2.4E+03
1.E-02 0

50

頻度

100

150

1.6E-03
1.5E-03
1.4E-03
1.3E-03
1.2E-03
1.1E-03
9.7E-04
8.7E-04
7.7E-04
6.7E-04
0

25
頻度

50

未臨界度と中性⼦⽣成時間の相関(1/2) 10
未臨界度と中性⼦⽣成時間の相関関係のモデル化
連続エネルギーモンテカルロ計算コードにより
KUCA-A架台EE1炉⼼体系を計算
計算条件
計算コード
MCNP6.2[3]核データ

JENDL-4.0[4]

ヒストリー数

100000

バッチ数

1300

捨てバッチ数300

燃料枚数(3008/3000/2520/2280/2240枚)と状態を変えて計算
p

p

p

p

p

p

p

p

p

p

p

p

p

p

p

p

p

p

p

p

p

p

p

p

p

p

p

p

p

p

p

p

p

p

p

p

p

p

p

p

p

p

p

p

p

C3

F

F

F

F

F

S5

p

p

C3

F

F

F

F

F

S5

p

p

C3

p

F

F

F

p

S5

p

p

C3

p

F

F

F

p

S5

p

p

C3

p

F

F

F

p

S5

p
p

p

p

F

F

F

F

F

p

p

p

p

F

F

F

F

F

p

p

p

p

F

F

F

F

F

p

p

p

p

p

F

F

F

p

p

p

p

p

p

F

F

F

p

p

p

p

p

F

F

F

F

F

p

p

p

p

F

F

F

F

F

p

p

p

p

F

F

F

F

F

p

p

p

p

F

F

F

F

F

p

p

p

p

p

F

F

F

p

p

p

p

S4

F

F

F

F

F

C1

p

p

S4

F

F

F

F

F

C1

p

p

S4

F

F

F

F

F

C1

p

p

S4

F

F

F

F

F

C1

p

p

S4

F

F

F

F

F

C1

p

p

F

F

F

F

F

p

p

p

p

F

F

F

F

F

p

p

p

p

p

F

F

F

p

p

p

p

p

p

F

F

F

p

p

p

p

p

p

F

F

F

p

p

p

p

C2

p

p

4

p

p

S6

p

p

C2

p

p

p

p

p

S6

p

p

C2

p

p

p

p

p

S6

p

p

C2

p

p

p

p

p

S6

p

p

C2

p

p

p

p

p

S6

p

p

p

p

p

p

p

p

p

p

p

p

p

p

p

p

p

p

p

p

p

p

p

p

p

p

p

p

p

p

p

p

p

p

p

p

p

p

p

p

p

p

p

p

p

p

p

[2]Christopher. J Werner, LA-UR-17-29981, (2017).
[3]K. Shibata, O. Iwamoto, et al., J. Nucl. Sci. Tecnol., 48[1], pp1-30, (2011).

未臨界度と中性⼦⽣成時間の相関(2/2) 11
未臨界度と中性⼦⽣成時間の相関関係のモデル化
主成分分析により相関関係をモデル化
2.0
1.5
Normalized Λ

t4

1.0

MCNP result

0.5

PC1

0.0
-0.5
-1.0
-1.5
-2.0

-1.5

-1.0

-0.5
0.0
Normalized ρ

0.5

1.0

1.5

中性⼦⽣成時間は反応度に対して負の相関スライド 11t4

載せる図間違えました。
今回の結果を計算した際の使ったデータはこちらです。
t-ikeda, 2018/11/23

実験概要

12

実験内容
KUCA-A架台
未臨界度過渡変化実験
ADS未臨界炉⼼体系
中性⼦源︓100 HzのDT外部中性⼦源
中性⼦計数の時系列データを測定→未臨界度を推定
p

炉⼼体系(2018年10⽉30⽇撮影)

p

p

p

p

p

p

p

p

p

p

p

p

p

p

p

p

p

p

p

C3

p

p

p

p

p

p

p

F

F

F

p

S5

p

p

p
p

F

F

F

F

F

p

p

F

F

F

F

F

p

p

p

S4

F

F

p

p

p

F

F

F

F

C1

p

F

F

p

p

p

p

C2

p

p

p

p

p

p

p

p

S6

p

p

p

p

p

p

炉⼼体系の模式図

実験体系
p

p

p

p

p

p

p

p

p

p

p

C3

p

13
p

p

p

p

p

p

p

p

p

p

p

p

p

p

p

p

F

F

F

p

S5

p

p

F

F

F

F

F

p

p

p

F

F

F

F

F

p

p

p

S4

F

F

F

F

F

C1

p

p

p

p

F

F

F

p

p

p

p

C2

p

p

p

p

p

S6

p

p

p

p

p

p

p

p

p

F

燃料体

p

C

制御棒

S

BF3検出器未臨界度過渡変化
⼩
𝜌
⼤

①定常状態
②C1制御棒挿⼊③C3制御棒落下④C2制御棒落下⑤中⼼架台落下(燃料・反射体引抜)

ポリエチレン減速材
安全棒
ファイバー検出器

実験解析条件(1/2)
観測値
 BF3検出器から得られた
中性⼦計数率p

C1挿⼊

C3落下 C2落下中⼼架台落下

600

p

p

p

p

p

p

p

p

p

p

p

p

p

p

p

p

p

p

p

p

p

p

C3

p

F

F

F

p

S5

p

p

p

F

F

F

F

F

p

p

p

F

F

F

F

F

p

p

S4

F

F

F

F

F

C1

p

p

p

p

F

F

F

p

p

p

p

C2

p

p

p

p

p

S6

p

p

p

p

p

p

p

p

p

Count Rate[s]

p

700

p

p

p

14

500
400
300
200
100
0
0

500

1000

1500
Time[s]

2000

2500

中性⼦計数率の時間変化

未臨界度と中性⼦⽣成時間の参照値
 MCNPによる計算結果を利⽤

3000

実験解析条件(2/2)

15

粒⼦フィルタ法の条件
粒⼦数︓10000 サンプル
時間ステップ幅︓1.0 s
観測ノイズ分散︓中性⼦計数に等しいと仮定
外部中性⼦源強度︓⼀定と仮定
システムノイズ︓𝜌は正規分布を仮定
Λは相関関係のモデルから算出後、モデルとのずれに基づき正規分布で摂動
2.0

Normalized Λ

1.5
1.0

MCNP result

0.5

PC1

0.0
-0.5
-1.0
-1.5
-2.0

-1.5

-1.0

-0.5
0.0
Normalized ρ

0.5

1.0

1.5

未臨界度の推定結果
Reactivity[dk/k]

0.00

C1挿⼊

16

C3落下 C2落下中⼼架台落下

-0.05
PF＋2σ

-0.10

PF-2σ
-0.15

Reference

-0.20
0

500

1000

1500
2000
Time[s]粒⼦10000サンプルの標準偏差𝜎
 2𝜎以内に参照値が含まれる

2500

3000

t8

中性⼦⽣成時間の推定結果
Generation time[s]5.0E-05

C1挿⼊

4.5E-05

17

C3落下 C2落下中⼼架台落下

PF＋2σ
PF-2σ
Reference

4.0E-05
3.5E-05
3.0E-05
0

500

1000

1500
Time[s]

2000

2500

3000

現在は熱中性⼦の計数率のみを計測しデータ同化している
 C1挿⼊前、中⼼架台落下後のように𝜌の変化によりΛが変化する場合は対応できる 制御棒挿⼊のように中性⼦スペクトルが変化するような場合には対応できないスライド 17t8

下の⽂章を書き換えました。
t-ikeda, 2018/11/23

18

まとめ & 今後の課題まとめ 粒⼦フィルタ法を⽤いた未臨界度推定の妥当性を検討
 未臨界度と中性⼦⽣成時間の同時変化を仮定して適⽤性検討
 KUCA未臨界過渡変化実験で得られた計測データを利⽤し未臨界度逆推定を実施
 未臨界度は不確かさ 2𝜎以内に参照値を含む
 未臨界度変化に伴う中性⼦⽣成時間変化には対応できるが制御棒挿⼊のような中性⼦スペクトル変化の場合には対応できない

今後の課題
反応度と中性⼦⽣成時間のモデルの改良
実効遅発中性⼦割合𝛽 や外部中性⼦源𝑆の時間変化推定

19

尤度関数

20

尤度(正規分布を仮定)
y
Hx
1
exp
𝑙
2R
2𝜋R
R :観測不確かさ分散、y ︓観測値
x ︓粒⼦、H ︓x をy と同じ物理量とするための⾏列重み(規格化された尤度)
𝑙
∝ exp
𝑤
∑ 𝑙

y

Hx
2R

核破砕中性子源によるウラン-鉛ゾーン炉心 ADS 体系における放射化反応率への中性子スペクトルの影響
Effect of neutron spectrum on reaction rate in ADS experimentwith U-Pb zoned core and spallation neutron source＊1

相澤直人 1, 山中正朗 2, 卞

哲浩 2

東北大学，2 京都大学複合原子力科学研究所

京都大学臨界集合体実験装置（KUCA）において、加速器中性子源によるウラン―鉛(U-Pb)ゾーン炉心 ADS（加速器駆動システム）体系を用いて、放射化反応率測定を実施した。その結果、通常燃料ゾーンと U-Pb
ゾーン間にて放射化反応率に関して中性子スペクトルの違いによる差異が測定された。
キーワード：加速器駆動システム（ADS）、中性子スペクトル、放射化反応率1.

緒言：加速器駆動システム(ADS)において、放射化反応率の測定は、炉内のスペクトル情報の把握ならびに ADS の中性子増倍の指標を表す中性子未臨界増倍率の推定を行うために重要である。本研究では、KUCA-A 架台において、ウラン―鉛(U-Pb)ゾーンを有する炉心と核破砕中性子源を組み合わせた ADS 体系(Fig. 1) [1]において In 線を炉内に設置し、熱中性子に感度を持つ
115

In(n,γ)116mIn および高速中性子に感度を持つ 115In(n,n’)115mIn 反応による放射化反応率を測定し、未臨界増倍率および熱／高速中性子比のスペ
クトル指標である In ratio[2]の測定を行った。
2.

実験結果と数値解析との比較：実験では、ポリエチレン減速通常燃料の本数を変更することにより、実効増倍率 0.98~0.95 の範囲内で 4 つの未臨界度を設定し、In 放射化反応率の測定を行った。115In(n,γ)116mIn

Fig. 1: U-Pb ゾーン炉心
赤線部は In 線設置位置である。

反応率分布に関しては、U-Pb ゾーンにおいてスペクトル硬化の影響に
より反応率が低くなる傾向が全測定結果にて得られ、数値解析結果とも概ねの一致を示した。この結果を用いて従来法[2]により未臨界増倍率を求めたところ、Fig. 2 に示すように実験値と計算値での不一致が見られ、原因を調査したところ、従来法では未臨界増倍率の導出に必要な換算係数の計算において、U-Pb ゾ
ーンのスペクトル硬化の影響を考慮できていないことが判明した。なお、スペクトル硬化の影響を考慮した換算係数を使用したところ、未臨界増倍率は実験値と計算値にて良い一致を示した。一方、115In(n,n’)115mIn 反応率については、計算値が実験値を大幅に下回る測定結果となった。これらの結果よりFig. 2: 未臨界増倍率の比較

In ratio を求めたところ、115In(n,n’)115mIn 反応率の差異が影響し、計算値と実験値の間に大幅な乖離が見られた。これらの要因について分析した内容については、発表にて報告の予定である。
参考文献[1] C. H. Pyeon and M. Yamanaka, KURNS-EKR-001, Inst. for Integrated Radiat. and Nucl. Sci., Kyoto Univ. (2018).
[2] C. H. Pyeon, et al., J. Nucl. Sci. Technol, 55(2), 190-198, 2018
* Naoto Aizawa1, Masao Yamanaka2 and Cheol Ho Pyeon2
1Tohoku Univ., 2Institute for Integrated Radiation and Nuclear Science, Kyoto University

核破砕中性子源によるウラン‐
鉛ゾーン炉心ADS体系における放射化反応率への中性子スペクトルの影響
東北大学 ○相澤直人，
京都大学複合研 山中正朗、卞哲浩
本研究は、平成29年度京都大学原子炉実験所（現 京都大学
複合原子力科学研究所）臨界集合体実験装置共同利用研究による成果である。

研究の背景・目的• 加速器駆動システム(ADS)における炉心特性の測定• 静特性実験による反応率測定• 炉内の中性子スペクトル情報の把握• 中性子束分布（反応率分布）の測定• 中性子未臨界増倍率𝑘 の推定• 中性子スペクトルが局所的に異なるウラン‐鉛(U‐Pb)ゾーン炉心と核破砕中性子源を組合せたADS実験体系において、放射化反応率測定を実施し、中性子スペクトルが放射化反応率及び反応率から求まる炉心特性に与える影響を明らかにする• 反応率分布• 未臨界増倍率 𝑘• In ratio（スペクトルインデックス）
第7回炉物理専門研究会，2018年11月27日

2

実験の概要・
炉心特性の導出手法

3
第7回炉物理専門研究会，2018年11月27日

実験の概要
Fig. 燃料集合体 ”F”• 実験体系• 炉心・燃料：• “F”: 1/8”p60EUEU• (1/8”PE+1/16”EU×2)を60セル重ねた燃料• “f”: 1/8”p10EUEU+ 1/8”Pb40EUEU+ 1/8”p10EUEU

Fig. 燃料集合体 ”f”• ”F”の中央40セルのPEをPbに置換した燃料• 制御棒・安全棒は全引き抜き、“F”の本数により、keff=0.98～0.95に設定• 中性子源：100MeV陽子ビーム+Pb‐Biターゲットによる核破砕中性子源• 放射化反応率測定• 炉内にIn線(φ1.0mm)、ターゲットにIn箔
(1×1×0.1mm)を設置、
In線は小片にカットし、放射化反応率を測定• 115In(n,γ)116mIn : 熱中性子の感度大• 115In(n,n’)115mIn: 高速中性子の感度あり（閾エネルギー0.3MeV）
第7回炉物理専門研究会，2018年11月27日

Fig. ADS実験体系図4

実験の概要• 実験体系• 炉心・燃料：• “F”: 1/8”p60EUEU• (1/8”PE+1/16”EU×2)を60セル重ねた燃料• “f”: 1/8”p10EUEU+ 1/8”Pb40EUEU+ 1/8”p10EUEU• ”F”の中央40セルのPEをPbに置換した燃料(a) Case A (≈0.98) (b) Case B (≈0.965)• 制御棒・安全棒は全引き抜き、“F”の本数により、keff=0.98～0.95に設定• 中性子源：100MeV陽子ビーム+Pb‐Biターゲットによる核破砕中性子源• 放射化反応率測定• 炉内にIn線(φ1.0mm)、ターゲットにIn箔
(1×1×0.1mm)を設置、
In線は小片にカットし、放射化反応率を測定• 115In(n,γ)116mIn : 熱中性子の感度大(c) Case C (≈0.955) (d) Case D (≈0.95)• 115In(n,n’)115mIn: 高速中性子の感度あり

5（閾エネルギー0.3MeV）
第7回炉物理専門研究会，2018年11月27日

未臨界増倍率の導出• 未臨界増倍率𝑘 の導出[1‐3]𝑘
𝐹: 核分裂中性子の総数, 𝑆: 源中性子の総数• 核分裂中性子•𝐹

𝑃𝜙 𝑟, 𝐸 𝑑𝑟 𝑑𝐸
𝐴·𝐶

𝑅𝑅

,

𝑥 , 𝑦, 𝑧 𝑑𝑦

• 源中性子•S

𝑠 𝑟, 𝐸 𝑑𝑟 𝑑𝐸

𝐶

𝑅𝑅

,

𝑟 𝑑𝑟

• 𝐶 ： 235U(n,f)と115In(n,γ)の反応率の換算係数• 従来手法：熱中性子領域(<0.1 eV)での235U(n,f)と115In(n,γ)の断面積の比例関係から導出• 比較手法：解析による235U(n,f)と115In(n,γ)の反応率比から導出[1] H. Shahbunder, et al., Ann. Nucl. Ene., 37, (2010), 592‐597, [2] H. Shahbunder, et al., Ann. Nucl. Ene., 37, (2010), 1214‐1222[3] C. H. Pyeon, et al., Nucl. Technol., 192, (2015), 181‐190
第7回炉物理専門研究会，2018年11月27日

6

In ratioの導出• KUCAにおけるADS炉心体系は235U高濃縮燃料により構成• 235Uにおける高速中性子(>0.1 MeV)および熱中性子(<0.1 eV)による核分裂比により炉心スペクトル指標を定義[4], ,235U fission ratio, ,• 235Uの高速/熱核分裂率の直接測定は不可能なことから、115Inの2つの反応から間接的に導出• 115In(n,Υ)116mIn: 熱中性子の感度大• 115In(n,n’)115mIn: 高速中性子の感度あり
,

, ,

,

, ,

235U fission ratio

7[4] C. H. Pyeon, et al., J. Nucl. Sci. Technol., 55(2), (2018), 190‐198
第7回炉物理専門研究会，2018年11月27日

数値解析による換算係数の導出・評価• PHITSとMVPにより実験値を評価、換算係数の導出を行う• PHITS: 核破砕反応と20MeV以上の中性子輸送計算• MVP: 20MeV以下の中性子輸送計算および反応率計算• ライブラリ:ENDF/B‐VII.1（輸送計算）、

JENDL/D‐99（中性子線量測定用特殊ライブラリ、反応率計算）• PHITS+MVPの解析の妥当性の評価（2018年秋の大会で発表）• 固有値計算（MVP単独）：

実験ベンチマークと概ね一致• 固定源計算（PHITS+MVP、実験の模擬）：実験値を概ね再現（後述）• 未臨界増倍率の導出のための換算係数の計算等に使用

8
第7回炉物理専門研究会，2018年11月27日

放射化反応率・
炉心特性の結果

9
第7回炉物理専門研究会，2018年11月27日

115In(n,γ)反応率分布の比較（全体）

Fig. Case Aにおける反応率分布

Fig. Case Bにおける反応率分布

Fig. Case Cにおける反応率分布

Fig. Case Dにおける反応率分布

第7回炉物理専門研究会，2018年11月27日

10

115In(n,γ)反応率分布の比較（炉心）U‐Pbゾーンの放射化反応率が比較的小さい⇒ U‐Pbゾーンの中性子スペクトルの硬化
Fig. Case Aにおける反応率分布
Fig. Case Bにおける反応率分布
炉心領域の局所的な中性子スペクトルの違いの影響が放射化反応率分布に反映されている

Fig. Case Cにおける反応率分布

Fig. Case Dにおける反応率分布11

第7回炉物理専門研究会，2018年11月27日

未臨界増倍率への影響• 領域毎に未臨界増倍率を導出• 𝐶 の導出手法の影響• 従来手法：熱中性子領域での235U(n,f)と115In(n,γ)の断面積の比例関係（≈0.309）

Fig. 従来手法による𝑘 の比較• 実験値と計算値の不一致の程度が大きい• U‐PEゾーンが比較的計算値に近い• 比較手法：235U(n,f)と115In(n,γ)の反応率比• 実験値と計算値が概ね一致• 燃料領域の違いに関係なく、いずれも同様の値

Fig. 比較手法による𝑘 の比較

12
第7回炉物理専門研究会，2018年11月27日

未臨界増倍率に対する中性子スペクトルの影響• In線における中性子束・反応率スペクトル• U‐Pbゾーンにおける熱中性子束が低い• 反応率について、2eV付近の共鳴による反応が顕著• 断面積の比例関係が成り立つ0.1eV以下の熱中性子による反応が支配的とはいえない• U‐PE領域にて30~40%@In / 60~70%@U• U‐Pb領域にて20~30%@In / 40~50%@U

Fig. CaseAにおける中性子スペクトル

Fig. CaseAにおける 115In(n,γ)の反応率スペクトル
13

第7回炉物理専門研究会，2018年11月27日

In ratioの導出結果①• In線の小片では115In(n,n’)によるγ線②スペクトルの検出が困難だったため、広範囲でのIn ratioを導出①• Case B: 20サンプル分(50cm分)• Case C: U‐Pb部 6サンプル（黄）,Case B

Case CU‐PE部 7サンプル（黄緑）• 高速領域： E > 0.111 MeV , 熱領域： E < 5.04 eV (文献[4]ではE < 0.1 eV)
ENDF/B‐VII.1

Case Bcase C (U‐Pb)case C (U‐PE)

case C (All)Calculation(2.54 ± 0.03)E+2 (1.56 ± 0.07)E+1 (3.80 ± 0.10)E+1 (2.91 ± 0.07)E+1

Experiment(2.00 ± 0.10)E+2 (7.91 ± 0.66)E+0 (1.29 ± 0.14)E+1 (1.13 ± 0.09)E+1

C/E

1.27 ± 0.06

1.97 ± 0.19

2.94 ± 0.34

2.57 ± 0.22

• 定性的には、熱／高速中性子比としては妥当な傾向• PE反射体を含むCase BでIn ratioが大きく、Case C(U‐Pb)で小さな値となる• 実験値と計算値に大きな差が存在• Case Cでは特に大きな隔たりが見られる
第7回炉物理専門研究会，2018年11月27日

14

In ratioに関する考察

235U fission ratio
,

, ,

,

, ,

• 計算値と実験値の不一致について、熱中性子/高速中性子の反応率を反応毎に分けて検討• 115In(n,γ) （熱中性子が主）は、反応率分布が概ね一致のため影響小• 115In(n,n’)：

高速中性子が支配的• この反応に限定してC/Eを比較したところ、1を下回る結果• Case B (炉心領域~PE領域の50cm分): 0.74~0.79• Case C, U‐Pb region: 0.53~0.56, Case C, U‐PE region: 0.36~0.38• 過去のPyeonの研究[4]において、放射化箔実験において、未臨界度が深くなるにつれてC/Eが減少する結果が得られており、測定数は少ないが同様の傾向が本研究でも得られた• 計算値はあくまでも20MeV以下の中性子による反応率であり、20MeV以上の中性子による反応が寄与している可能性がある[4] C. H. Pyeon, et al., J. Nucl. Sci. Technol., 55(2), (2018), 190‐19815

第7回炉物理専門研究会，2018年11月27日

結論

16
第7回炉物理専門研究会，2018年11月27日

結論• U‐Pbゾーン炉心＋核破砕中性子源によるADS静特性実験にて反応率分布測定を実施し、測定結果を元に未臨界増倍率・In
ratioの導出を行い、中性子スペクトルの影響を調査した• 115In(n,γ)放射化反応率分布• 炉心の中性子スペクトルの違いを定性的に示す分布が得られた• 未臨界増倍率• 𝐶 の導出に従来手法(熱領域の断面積比)と比較手法（反応率比）を適用• 従来手法では中性子スペクトル硬化の影響により適切な未臨界増倍率の導出不可、中性子スペクトルの違いを考慮した導出が必要• In ratio• 導出した値は、定性的には適切な熱/高速中性子の比を表す数値• 全体的なC/Eの不一致が見られた⇒20MeV以上の中性子による反応の寄与を明らかにする必要がある

17
第7回炉物理専門研究会，2018年11月27日

KURNS REPORT OF
INSTITUTE FOR INTEGRATEDRADIATION AND NUCLEAR SCIENCE,
KYOTO UNIVERSITY
発行所
発行日
住所

京都大学複合原子力科学研究所
平成 30 年12月
大阪府泉南郡熊取町朝代西二丁目
TEL（072）451- 2300

掲載された論文等の出版権、複製権および公衆送信権は原則として京都大学複合原子力科学研究所に帰属する。
本誌は京都大学学術情報リポジトリに登録・公開するものとする。 http://repository.kulib.kyoto-u.ac.jp/dspace/