2006核データ研究会

核データチュートリアル（２）MVPの利用

日本原子力研究開発機構
森 貴正
1

モンテカルロ法の基礎- モンテカルロ法による中性子輸送シミュレーション- 乱数とサンプリング- 物理量の評価法- MC法による固有値計算の困難

MVPの概要- 機能- 幾何形状表現法- 物理モデル・核データ表現法

MVPの利用例
2

モンテカルロ法とは？
z 乱数(random number)を用いて数学的問題を解く
数値解法
z 統計的手法(確率論的手法) ⇔ 決定論的手法
例：数値積分

ξ1 , ξ 2 は[0,1]の乱数

y
1

001x

ξ12 + ξ2 2 > 1

カウントせずξ12 + ξ2 2 ≤ 1

カウントする

カウント数
青色の π= ≈
面積4
総試行回数

3

モンテカルロ法による中性子輸送シミュレーション
飛行する個々の中性子を物理法則と核データに従う確率を用いて追跡する（ランダムウォーク）
散乱
核分裂
捕獲

捕獲

手 順
核分裂
散乱

捕獲

1. 体系を正確にモデル化(計算機の中に炉心を構築)
2. 大量の中性子の飛行追跡
3. 統計量から集団的振舞いを評価
本質的な近似が少ない

4

ランダムウォーク：自由飛行過程
次の衝突位置までの距離 s の確率密度関数（確率分布）s⎡
P( s )ds = Σ t ( s, E )exp - ∫ ds ′Σ t (r ′ + s ′Ω, E )⎤ ds⎢⎣ 0⎥⎦

手 順(1) 無限体系を仮定
P ( s )ds = Σ t exp[- Σ t s ]ds1
s = − ln(ξ )Σt(2) 領域境界までの距離 d の計算

r′

例：右図（中心r0、半径Rの球の内部）

l = −Ω ⋅ (r ′ − r0 ) ± R − (| r ′ − r0 | − | Ω ⋅ (r ′ − r0 ) | )2

簡潔なアルゴリズム

22Ωr0d(3) sとdの比較Æ今の領域内で衝突するかどうか決定(4) 今の領域から出た粒子が入射する領域の探索

5

ランダムウォーク：衝突過程
衝突後のエネルギーEと飛行方向 Ωの分布（連続エネルギー法）
C ( E ← E ′, Ω ← Ω′) = ∑m

ρ mσ tm ( E ′)σ rm ( E ′) m m
ν r f r ( E ← E ′, Ω ← Ω′)∑mΣ t ( E ′) r σ t ( E ′)

A

手 順

BC

D

(1) 衝突核種の決定

A：核種 m と衝突する確率(2) 反応タイプの決定

B： 核種 m との衝突で、反応タイプが r である確率(3) 衝突から現れる粒子数の決定
C:核種 m 、反応タイプ r の衝突から出てくる粒子数の期待値(4) 衝突後のエネルギーと飛行方向の決定
D：核種、反応毎の確率分布
弾性散乱等 ： 散乱角分布＋エネルギー・運動量保存則6
その他 ： 散乱角・エネルギー分布

モンテカルロ法の特徴
z 確率的な事象を忠実に模擬しているので、近似が少ない。
利
点 z 位相空間を離散化する必要がない。
¾ 体系を正確に(as bulit)に模擬できる。
¾ 連続エネルギー的取り扱いが可能。

ほぼ核データの精度で計算できる。
z 統計的手法であるので解にばらつきがある。
欠 z 解の精度を上げるために試行回数を増やさなけ
点ればならない。
z 計算時間がかかる。
7

乱数とサンプリング積算確率を用いる方法+∞∫ dx f ( x) = 1,確率密度関数 f (x )：f ( x) ≥ 0,積算確率関数 F (x )：dF ( x)= f ( x).
F ( x) = ∫ dx f ( x),
−∞
dx−∞x

確率変数 ρ = F (x ), その確率密度関数を p ( ρ )dF
dx = f ( x ) dx
p ( ρ ) dρ = f ( x ) dx, dρ =dx
p ( ρ ) = 1 = 一定 (0 ≤ ρ ≤ 1).
一様乱数 ξ (0 ≤ ξ ≤ 1).
ξ = F (x ) を満たす xk −1

離散的事象

f(x)からのサンプリングk∑ P < ξ ≤ ∑ P , ( k ≤ n)
i =1

i

i =1i

事象kを選択8その他の種々のサンプリング法
参考文献３）：”A Monte Carlo Sampler”,LA-5061-MS(1972)

例１：Evaporation spectrumからのサンプリングf ( E ) ∝ E exp(− E / T ),E = −T ln(ξ1 ⋅ ξ 2 ),T : 核温度

ξ1 , ξ 2 : 一様乱数

例２：Simple fission spectrumからのサンプリングf ( E ) ∝ E exp(− E / T ), T : 核温度⎧⎫2 π
E = −T ⎨ln ξ 3 ) + ln ξ1 ⋅ cos ( ξ 2 )⎬, ξ1 , ξ 2 , ξ 3 :一様乱数29⎩⎭

粒子輸送モンテカルロ法で用いられる乱数：
合同乗算法(Linear Congruential Generator)• Lehmerによって提案（1951）• 乱数ξは次の漸化式により生成

Si = (aSi −1 + c) mod mSi
ξi =m• c = 0のとき乗算合同法（multiplicative CG）• c ≠ 0のとき混合合同法（mixed CG）
10

粒子輸送MCコードの線形合同法
コード名
MCNP5 (default)(optional)

ma519
248
263 9,219,741,426,499,971,445
2,806,196,910,506,780,709
3,249,286,849,523,012,805

MCNP (Version 4) 248519c

周期0111246263
263
263246263
231
245
246
2450

MVP (Version 2)

263 9,219,741,426,499,971,445

MVP (Version 1)

231

32771

1,234,567,891

RACER

247

84,000,335,758,957

0

VIM

248

5190MORSE

247

515

0

RCP

248

29+1

59,482,192,516,946

148211物理量の評価法（エスティメータ）中性子束の評価によく用いられる方法
飛程長エスティメータ
小領域 V 内の各トラックからの寄与 Wi ⋅ l i / V の総和

計算体系

領域 A の平均中性子束

衝突エスティメータ
小領域 V 内の各衝突からの寄与 Wi / Σ t ,i / V の総和
領域 A の平均中性子束面交差エスティメータ

面Al1l2

小領域 Vl4l5l3

Wi
面 A を横切る毎の寄与
A n ⋅ Ωiの総和
面Aでの平均中性子束（ カレント評価にも利用可）

12

固有値（keff）計算に用いられる方法
核分裂源反復世代毎の追跡

核分裂で中性子の一生は終了
発生する核分裂中性子を保存
核分裂源分布の収束後

次世代の線源

実効増倍率等の評価

実効増倍率keffの評価法
中性子増倍による評価
νΣ f φdEd r
核分裂により生成され る中性子の重み∫=keff =
source neutrons
源の中性子の総重み

中性子バランスによる評価
νΣ f φdEd r
核分裂により生成され る中性子の重み∫keff ==
漏れ＋正味の吸収
leakage + ∫ Σ loss φdEd r

断面積を含む項は飛程長エスティメータ等で評価
漏れは面交差エスティメータで評価

13

誤差評価（期待値と分散の推定）• n回の試行 Æ n個のサンプル Xi• サンプル間の独立性を仮定

1 n
期待値の推定： X = ∑ Xi
n i=1()

n21
分散の推定： σ ( X ) =Xi − X∑n(n −1) i=12n()2
12X −X ,=n −12X∑ iX 2 = i=1n14

MC法による固有値問題の困難（１）z 固有値の分散のバイアス¾ 世代（バッチ）間の相関があるために固有値の分散は過小評価される。
¾ スーパーヒストリー法¾ Ueki等の方法（MVPで採用している方法）¾ 時系列解析法
z 収束判定¾ 線源分布（核分裂源分布）が何世代（バッチ）で収束するのか分からない。
¾ 収束の速さはドミナンス比に依存する。
¾ 核分裂源分布の収束は固有値の収束より遅い。
15数値計算例：分散のバイアス

Case#1 一次元平板60cm

σ /σ2R

Case#1

1.43

Case#2

2.13

Case#3

3.07

Case#42A

5.62

k ∞ = 1.0
Case#22.4cm

k ∞ = 0.8
Case#31.2cm

k ∞ = 1.3

k ∞ = 0.8

σ : 500ラン毎の σ S2,i2Aの単純平均

σ R2 : 500ラン毎の k i のCase#4 3次元円筒z35cm

k ∞ = 0.65

ばらつきからk ∞ = 1.01

1639.5cmr

数値計算例：固有値と核分裂源分布の収束性
z Keff is an integral quantity - converges faster than source shape
Keff calculation for 2 nearly symmetric slabs,
with Dominance Ratio = .9925
value divided by true mean

1.11keff

0.90.80.7

source at right fuel0.60.50.40100200300

4005006007008009001000cycle
17
Cited from F. Brown’s viewgraph at Monte Carlo 2005. Calculated by T. Ueki, UNM

MC法による固有値問題の困難（２）z アンダーサンプリング¾ Whitesidesによって提案された”k-eff of the world”の問題¾ 結合の弱い体系¾ 中性子が入射しないユニットがあることにより、固有値が正しく評価されない。
¾ バッチサイズを大きくする。
¾ Stratified Sampling

未臨界臨界超過18

¾ MVP : (Monte Carlo code for Vector Processors)• 連続エネルギー法に基づく汎用モンテカルロコード¾ GMVP : (Group-wise Monte Carlo code for Vector Processors)• 多群法の基づく汎用モンテカルロコード
z 1980年代後半より日本原子力研究機構で開発されているモンテカ
ルロコード。
z 1994年に第1版を国内公開。2005年に第2版を世界に向けて公開。
z MVP/GMVPコードではベクトル計算に適したアルゴリズムを採用。
（事象駆動型アルゴリズム ⇔ ヒストリー駆動型アルゴリズム）z ベクトル計算機上ではスカラーモンテカルロコードより10倍以上の高速化を達成。（スカラー計算機では計算速度はほぼ同じ。）
z 汎用並列化ライブラリーMPI, PVMを用いた並列化にも対応。
19

MVPの機能z 解く問題： 中性子・光子輸送に対する固有値・固定源問題。
時間依存問題も可能。
z エネルギー範囲： 中性子10-5eV – 20MeV, 光子1keV – 100MeV
z 幾何形状表現： 組み合わせ幾何形状表現（CG）、多重格子形状。
z 粒子源表現： 多数のサンプリング関数により柔軟な指定が可能。
z 断面積：
9 LICEMコードによって評価済み核データから処理された専用のライブラリを使用。（JENDL-3.3, END/B-VII, JEF-3.1, etc.）
9 任意温度における計算が可能。
z 分散低減法： ロシアン・ルーレット、スプリッティリング、インポータ
ンス、ウエイト・ウインドウ、パス・ストレッチング、線源バイアス
z エスティメータ： 飛程長、衝突、面検出器、点検出器が利用可能。
固有値はいくつかの評価法より最尤法により評価。Uekiの方法によ20
る真の分散の推定。

MVPの幾何形状表現法
組合わせ形状表現（CG）
Bodyと呼ぶ単位形状の組合せ（集合演算）でCG幾何形状を表現利用できるBody+

直方体、円柱、球、円錐台、三角柱、平行六面体、正六角柱、トーラス、一般二次曲面など
Pre-defined bodies

多重格子表現Multiple Lattice

複雑な繰り返し形状の表現
多重の四方及び六方格子21

確率論的幾何形状モデル（１）• 高温ガス炉（HTR）の開発• 膨大な数のランダムに分布する被覆燃料粒子（CFP）、燃料ぺブルの取り扱い
Coated fuelparticle(CFP)
Fuelcompact

Prismaticblock

Fuel pebble

個々の粒子位置を知る事は困難
入力データ作成が困難
非均質性は炉心特性に重要

Moderator
pebblePebble bed
幾何形状の確率論的取り扱い22確率論的幾何形状モデル（２）個々の中性子に対して燃料粒子を確率分布（最近接セル分布、NND）を用いて順次配置しながら追跡する
燃料粒子の存在する領域

2６４

注目する中性子の飛行に関係する燃料粒子のみを配置する
使用するNND(1) モンテカルロ最密充填法（MCRDFコード）５
1

3
原子核との衝突(2) 統計的一様性を仮定した理論的近似式⎛ 3 fp⎞dNND( r ) 3 f p= ⋅⋅ exp⎜ − ⋅⋅ r⎟⎜⎟dr
2 1− fp⎝ 2 1− fp ⎠f p : 充填率23計算例：HTTR全炉心解析

0.98975±0.09%

確立論的幾何形状ﾓﾃﾞﾙ
1.00092±0.14%―――

1.13±0.1493 万47 万10 時間

12 時間

均質ﾓﾃﾞﾙ
実効増倍率(k)非均質効果(%Δk/kk’)
追跡粒子数計算時間(VPP-500)

３０燃料カラム装荷炉心*JENDL-3.2 を使用
J

Cf. ２１カラム装荷炉心における非均質効果
0.52±0.04 %Δk/kk’
24（HTTR／安藤氏等)

MVPの物理モデル・核データ表現法
ENDF-6形式で表現されたデータのほぼ全てを利用

zS.C.TzFree gas model
25

MVPライブラリー作成システム：LICEM
ENDF-6selected
tape

LINEAR

ENDF-6
thermal
tape

U3R-J

THERMJ-B6

U3RPDS01~ 11THERMALMVPLIBMKWork file

RECENT
SIGMA1NJOY
HEATR

MVP libraryMaterial file26

中性子反応断面積¾ 分離共鳴領域＋スムーズ断面積
線形内装可能なpoint-wise断面積¾非分離共鳴領域

断面積確率テーブル27

断面積確率テーブル法の効果
NEA/NSC MOX燃料PWRベンチマーク

U-238捕獲断面積
By MVPMOX

PWRセル――――――――――――――――――――――――――
k-infinity
非ボイドセルボイドセル――――――――――――――――――――――――――
確率テーブル1.2160±0.03% 1.2694±0.02%
無限希釈断面積1.2155±0.03% 1.2612±0.02%
Δk0.0004
0.0082±4.4%――――――――――――――――――――――――――28

ν値（MF=1）
評価済み核データに従う¾多項式表示¾テーブル表示Æ線形内挿可能なエネルギー点

遅発中性子のデータ（MF=1、MF=5）
評価済み核データに従い、 νd値、λi, βi, χd,i を格納¾νd値： ν値と同じ表現¾βi, χd,i：エネルギー分布（サブセクション）29

角度分布（MF=4）¾テーブル表示Æ32個の等確率ビンで表現（熱中性子散乱を除く、核分裂は等方を仮定）

30

エネルギー分布（MF=5)
評価済み核データに従う¾テーブル形式（LF=1、5）
入射エネルギー毎に、線形内挿可能2次エネルギービンに対するテーブル（内挿法に従い使用するテーブルを選択）¾関数形式（LF=7(Maxwellian)、9(Evaporation)、11（WATT））¾Madland and Nix fission spectrum（LF=12）
テーブル形式に変換（LF=1と同じ）

31

角度エネルギー分布（MF=6）
評価済み核データに従い、次の形で取り扱う¾Kalbach-87分布¾多体フェーズスペース分布¾実験室系角度エネルギー分布¾角度分布のみ（2体反応角度分布）32熱中性子散乱データ¾干渉性弾性散乱（MF=7、MT=2、LTHR=1）
評価済み核データ（ブラッグエッジエネルギー等）をそのまま利用¾非干渉性弾性散乱（MF=7、MT=2、LTHR=2）
評価済み核データ（Debye-Waller integral）Æ32の等確率ビン（ μ1 ≤ μLab ≤ μN+1）、エネルギーは変化せず。
¾非干渉性非弾性散乱（MF=7、MT=4）
162のエネルギー点（10-5～12eV）に対する散乱マトリックス
各エネルギートランスファーに対する角度分布：Double P1近似：
5つの事象（前方(μ >0)等方、後方(μ <0)等方、μ=1、0、-1）に対する確率を与える¾Free gas model33
マックスウェル分布の標的原子核と中性子が重心系等方散乱任意温度における計算の準備
New Data Processing Scheme

JENDL-3.3

New LICEMNew MVP library

Interpolation parameters• Thermal scattering data
• Data in the unresolved resonance region
MVP code
ART code$CROSS SECTION& IDMAT( 1 )• Doppler broaden data in theTEMPMT( 600.0 )resolved resonance region
U02350J33( 5.070E-5 )• Interpolate data for thermal
U02380J33( 1.685E-2 )scattering/unresolved resonance$END CROSS SECTION
34

MVP-BURNにおける燃焼計算手法
ART
Doppler Broadening, etc.

MVP Library(RoomTemperature)

ENDF/B-VI
JEF-2.2

LICEM

JENDL-3.2

MVP-BURN
User’s LibraryR(Arbitrary (Predictor-Corrector Method)
Temperature)
Microscopic
Depletion Calculation
Reaction Rates
by Bateman’s method
newMonte CarloCalculation

R( i ) + R old ( i + 1)
( i + 1) =2

MVP
Input datafor MVP

BURN

sub-step

MVP step(i)NuclideDensities

MVP step(i+1)

35核データの積分テスト
詳細な実験解析
決定論的手法の検証と補完

36

核データの積分テスト
by K. Okumura (2006)37

遅発中性子スペクトルの効果 (by Y. Nagaya)
ケース

燃料漏れ

Δk（DNあり-なし）(1σ)

Godiva

高濃縮U

57%

ORNL-10高濃縮U

7%0.00012 (0.00003)
0.00039 (0.00003)

BIGTEN

中濃縮(10% 235U)

11%-0.00183 (0.00005)

TCA(1.83U)

低濃縮U
19×19

4%

0.00079 (0.00008)

TCA(3.00U)

低濃縮U
19×19

6%

0.00069 (0.00008)

JezebelPu

67%Pu(NO3)4Pu

7%-0.00023 (0.00003)
0.00008 (0.00004)

38
ICSBEPベンチマークから選択。DNあり、なしでそれぞれ10回ずつ実行。

遅発中性子(by Y. Nagaya)

U-235 spheren (t ) = l∑j =1

Point source

attenuationof promptneutrons

ωtm +1e jβ i λi
2i (ω j + λi )l+∑l : Generation time

βi : Delayed neutron fraction

λi : Decay constantω j : Root of inhour equation
Delayed neutrons

39

FCAボイド反応度価値測定実験の詳細解析
FCA II-1炉心Blanket(DUB)

SS

ボイド反応度価値測定で置換された領域

Half height (cm)66.04

Nat-Ublock(NUB)
45.72

0.00

9Z
8Z
7Z
6Z
5Z
4Z3Z2Z
1Z0.001-9Z (65%→95%)
4-9Z (65%→95%)
7-9Z (65%→95%)
1Z (65%→95%)
1Z (65%→80%)

Blanket(SB)

ポリスチレン板

CoreDriver

Al2O3 Depleted-UO
2Pu

24.84

35.88

63.48

Distance from core center (cm)80.04

Depleted-UO2

Nat-U
40

ボイド反応度価値（微分演算子サンプリング法、JENDL-3.3）

C/E=1.034±2.3%
10億h, (2007年）

41空間依存共鳴吸収計算法の計算精度の検証
By T. Kugo燃料ペレットを径方向分割した燃料棒

U-238吸収断面積（barn)（詳細解：MVP）
詳細解80

従来法

新手法

60
40
20000.1

0.20.3半径(cm)

0.40.5

中性子増倍率詳細解 1.3719(MVP)
従来法 1.3612(－0.8%)
新手法 1.3706(－0.1%)42

散乱マトリックスタリー
17-group macroscopic scattering matrices
Al clad
2.6%
UO2

MVP

H 2OSRAC(collision probability)43最後に、MVPの最近の話題
z 摂動計算機能¾ 実効増倍率への密度摂動、核種数密度摂動¾ 核分裂源分布の変化に対する摂動効果¾ 微分演算子サンプリング法に対する高次の摂動効果¾ 断面積摂動、温度摂動
z 多群定数計算機能¾ 散乱行列タリー¾ 拡散係数計算機能（等方、非等方拡散係数）¾ 輸送断面積用定数（高次のルジャンドル展開散乱係数）44

45