粒子輸送モンテカルロ計算の基礎原理とそこに潜む不確かさ
日本原子力研究開発機構

1

長家

康展

はじめに
モンテカルロ法は乱数 (random number) を用いて数学的問題を解く数値解析手法である。方程式を離散化

して数値的に解く決定論的手法と対比して、モンテカルロ法は確率論的手法 (統計的手法) とも呼ばれる。
モンテカルロ法に関連した統計手法の起源は古く、18 世紀にまでさかのぼる。18 世紀後半、Buﬀon は、平行線に何回も針を投げ、針が平行線と交わる確率から実験的に円周率 π を求める方法を編み出した [1] 。この実験
は”Buﬀon の針実験” として知られている。19 世紀後半になると、Lord Rayleigh は、吸収妨害物 (absorbing

barriers) のない 1 次元のランダムウォークが 2 次微分方程式の近似解を与えることを示した [2]。20 世紀に入ると統計的手法に関する研究は盛んになり、Lord Kelvin は、気体の運動理論に現れる時間積分を実行するためにランダムサンプリングを用いた [3]。Courant ら [4]、Kolmogorov[5]、Petrowsky[6] はランダムウォーク過程やマルコフ連鎖過程と微分方程式の関係を明らかにしていった。また、Fermi は、1930 年代前半にローマにおいて中性子拡散に関する研究で統計的手法を用いていた。
しかし、まだこの当時、電子計算機はなく、モンテカルロ法の基礎も確立していなかった。それが確立したのは 1940 年代に入ってからである。1945 年に世界初の電子計算機 ENIAC が完成し、von Neumann と Ulamは、その電子計算機上でランダムウォーク過程に基づくサンプリング実験を行うことを提案した。また、von

Neumann、Ulam、Metropolis はモンテカルロ法の数学的な基礎を確立した [7, 8]。モンテカルロ法という名前
はこの頃 Metropolis によって提案され、この統計的手法を表わす言葉として用いられるようになった [9]。
モンテカルロ法は、乱数を用いた統計的な手法であるので、確率的な事象に関する広範な問題を解くことができる。現在では、様々な分野で用いられ、原子力の分野 (炉心設計、臨界安全、遮蔽解析) はもちろんのこと、統計物理、素粒子物理、プラズマ物理、放射線防護、宇宙物理、更には石油探査、金融工学の分野でも用いられるようになっている。このようにモンテカルロ法の利用範囲が広まった一因は計算機技術の進歩にある。
モンテカルロ法は統計的手法であるので、その結果にはばらつきがあり、解の精度を上げようとすると必然的に試行回数を多くしていかなければならない。これはモンテカルロ法の欠点であるが、ベクトル化や並列化といった計算機技術の進歩と最近の低価格で高性能のパーソナルコンピュータの出現により、現実的な時間で十分な精度の結果が得られるようになってきている。
モンテカルロ法が広範な分野で用いられるようになったのは、計算機技術の進歩の他にも注目すべき利点がいくつかあるからである。モンテカルロ法は、確率的な事象を忠実に模擬することにより、興味のある物理量を評価する手法であるので、決定論的手法のように物理現象を記述する方程式を差分化して解く必要はない。
それゆえ、決定論的手法では複雑すぎて解けないような問題もモンテカルロ法では解くことが可能となる。また、位相空間 (空間、エネルギー、時間等) の離散化をする必要がないので、計算体系を正確に記述できるという利点がある。従って、物理事象を正確にモデル化しておけば、近似を伴わない正確な物理量を求めることができる。
様々な分野で用いられているモンテカルロ法であるが、ここでは原子炉分野で用いられる粒子輸送モンテカ
ルロ法について述べる。即ち、取り扱う粒子は中性子で、エネルギー範囲は 10−5 eV から 20 MeV であるものを対象とする。また、中性子と原子核との反応は評価済み核データに記述されているものとする。このような粒子輸送モンテカルロ法に基づく汎用のモンテカルロコードは既に数多く開発されており、最も有名なものではロスアラモス国立研究所 (LANL) で開発された MCNP[10]、日本原子力研究開発機構 (JAEA) で独自に開発された MVP[11] などがある。これらコードをブラックボックスとして使用しても、求める解は得られるかもしれないが、その解が意味のあるものになっているのかどうか判断するのはユーザ次第である。本稿では、モンテカルロコードの仕組みを理解できるように粒子輸送モンテカルロ計算の基礎原理について述べ、モンテカ
ルロコードで用いられている計算手法について解説する。また、モンテカルロ計算に潜む不確かさについても言及し、モンテカルロコードのユーザが入力データを作成したり、計算結果を解釈する際に注意しなければな
らない点について述べる。

モンテカルロ法の原理

2
2.1

乱数

乱数はモンテカルロ法の基礎であり、モンテカルロ計算において非常に重要な役割を果たしている。モンテ
カルロ計算では、どの事象が起こるのかを確率論的に決定し、忠実に対象となる現象をシミュレーションしているが、乱数を用いてある事象に対する確率分布からのサンプリングを行い、事象を決定している。
乱数を発生する仕組みは乱数ジェネレータと呼ばれ、通常、乱数はソフトウェア的に (計算機上のアルゴリ
ズムにより) 生成される。しかし、計算機上の変数が有限ビットで表わされるため、ソフトウェア的に真の乱数を生成することは不可能であり、乱数列は周期を持つことになる。即ち、同じ乱数が現れると、それ以降、前の乱数列と全く同じ乱数列になってしまうことになる。それゆえ、このようなソフトウェア的な乱数ジェネ
レータで生成した乱数は擬似乱数 (pseudorandom number) と呼ばれる。
擬似乱数を生成するアルゴリズムはこれまで数多くのものが提案されてきているが、粒子輸送モンテカルロ
計算では Lehmer によって提案された線形合同法 (linear congruential generator)[12] が最も多く用いられている。この方法では次の漸化式により、0 から 1 の間の一様乱数を生成する。

Si
ξi

=

(aSi−1 + c) mod m
Si
=
m

(1)
(2)

ここで、a, c, m, Si は正の整数である。また、m を法 (modulus) とする剰余演算から分かるように Si のとりう
る範囲は 0 ≤ Si < m であり、ξi のとりうる範囲は 0 ≤ ξi < 1 である。しばしば、c は 0 とされることが多いが、このときのアルゴリズムは乗算合同法 (multiplicative congruential generator) と呼ばれる。これに対し、

c 6= 0 のときは混合合同法 (mixed congruential generator) と呼ばれる。
線形合同法は原理が単純で、高速に乱数を生成できるという利点がある。一方、他のアルゴリズムに比べて周期が短く、連続した乱数をいくつの組にして座標にプロットすると規則的に並ぶ [13] といった欠点もある。
しかしながら、通常の粒子輸送モンテカルロコードで用いられる線形合同法は周期が 245 ≈ 3.5 × 1013 以上のものが多く、通常の計算において周期は十分に長いので、ランダム性 (randomness) は保証されていると考えても良い。
最近では、計算機の進歩により、モンテカルロコードによる大規模計算がよく行われるようになってきている。乱数の周期以上の個数の乱数を使ったからといって、すぐに計算結果がおかしくなるということはないが、このような大規模計算のために MCNP5 のオプショナル乱数ジェネレータや MVP Version 2 の乱数ジェネレー
タのように周期が 263 ≈ 9.2 × 1018 である線形合同法も用いられるようになりつつある。

2.2

サンプリング法

サンプリング法はモンテカルロ法の基礎となる手法であり、乱数を用いてある確率分布からどのような事象が起こるかを決定する。今、事象 E1 , · · · , En は独立で、それぞれ確率 p1 , · · · , pn で起こるものとする。このとき、p1 + · · · + pn = 1(3)

p1 + · · · + pi−1 ≤ ξ < p1 + · · · + pi

(4)であり、乱数 ξ がの範囲にあるとすると、事象 Ei が起こったと決定する。
一例として、中性子が原子核と衝突し、どのような事象が起こるかを決定することを考えてみる。簡単のため、吸収、弾性散乱、非弾性散乱しか起こらないとするし、それらの断面積をそれぞれ Σ a , Σ el , Σ ie とする。
全断面積は

Σ t = Σ a + Σ el + Σ ie

(5)であり、各反応が起こる確率は p1 = Σ a /Σ t , p2 = Σ el /Σ t , p3 = Σ ie /Σ t である。このとき、乱数 ξ の値により次のように反応を決定することができる。

0≤ξ<

Σa
Σt

→

吸収

(6)

Σa
Σa
Σ el
≤ξ<
+
Σt
Σt
Σt
Σa
Σ el
+
≤ξ<1
Σt
Σt

→ 弾性散乱(7)→ 非弾性散乱(8)

図 1はこのサンプリングに対する乱数と各反応の起こる確率を示している。

Σa
Σt

Σ el
Σt

Σ ie
Σt
1.0

ξ

0.0

for various types of neutron
Fig. 2 Probabilities
図 1: 中性子反応に対する確率
reactions
この例では、確率密度関数 (probability density function, PDF) がステップ関数で与えられる場合のサンプ
リングを示している。このときの確率密度関数は一般的に次のように定義することができる。

p(x) = pi , (i − 1 ≤ x ≤ i) for i = 1, 2, · · · , n(9)

この関数を図で表わすと図 2(a) のようになる。そこで、積算密度関数 (cumulative probability function, CDF)
を次のように定義する。

∫ xp(t)dt, (0 ≤ x ≤ n)

P (x) =(10)0すると、図 2(b) に示すように P (0) = 0, P (n) = 1 の単調増加関数となる。
1
p

n-1

p

P(x)

1

p(x)

p

p

2

n

p +p
1

p

p
0
0

1

2

2

ξ

3

n-1

3

1

0

n

0

x

(a) 確率密度関数

1

2

3

n-1

n

x

(b) 積算密度関数

図 2: 離散確率密度関数と積算密度関数

ここで、乱数 ξ と積算密度関数を次のように結びつけると、x は ξ の関数として一意的に決定することができる。

ξ = P (x)(11)

即ち、x が i − 1 ≤ x < i に入る確率は Pi であり、それに対応した事象を決定することができる。従って、乱数を使って事象を決める際に必ず (11) 式の逆関数

x = P (ξ)−1(12)

を求める必要がある。
これまで確率密度関数が離散的な (ステップ関数で表わされる) 場合について述べてきたが、一般的には連続的な確率密度関数が与えられる場合がある。この場合についても同様にして、乱数を用いてある事象をサンプリングすることができる。確率密度関数が区間 a ≤ x < b で定義され、図 3(a) に示されるのもとすると、積算
密度関数は次式で定義される。

∫ x
P (x) =

p(t)dt
a

(13)

1

P(x)

p(x)

ξ

0

b

a

a

x

(a) 確率密度関数

x(ξ)

x

b

(b) 積算密度関数

図 3: 連続確率密度関数と積算密度関数

積算密度関数は図 3(b) に示すように単調増加関数になるので、乱数が 0 から 1 の間で一様分布していれば、乱数 ξ に対応した x(ξ) が一意的に決定できる。即ち、x(ξ) = P (ξ)−1

(14)

である。
連続確率密度関数の例として、乱数を用いて飛行距離を決定することを考える。今、中性子が原点 (x = 0)にあり、最初の衝突を距離 x と x + dx の間で起こす確率 p(x)dx はp(x)dx = Σ t exp(−Σ t x)dx
で与えられる。ここで、0 ≤ x ≤ ∞ である。(15) 式を x のとりうる範囲で積分すると∫ ∞∫ ∞p(x)dx =Σ t exp(−Σ t x)dx = 1
0

(15)

(16)

0

であり、p(x) は確率密度関数とみなすことができる。乱数 ξ から飛行距離 x を決定するために積算密度関数を
計算すると

∫ x
Σ t exp(−Σ t x)dx = 1 − exp(−Σ t x)

P (x) =

(17)0であり、次式より飛行距離 x を決定することができる。

ξx1 − exp(−Σ t x)1
= −ln(1 − ξ)
Σt
=

(18)
(19)

また、(1 − ξ) の分布と ξ の分布は同じであるから、次式で飛行距離 x を決定してもよい。

x=−

1
ln ξ
Σt

(20)

次に、2 変数の連続確率密度関数の例として、等方的に散乱された粒子の飛行方向を決定することを考える。
微小立体角 dΩ は極角 (polar angle)θ と方位角 (azimuthal angle)ω を用いて次式で表わされる。

dΩ = dθ sin θdω(21) 式を全立体角で積分すると

∫

∫ π
dΩ =

(21)

∫ 2π
dθ sin θ0dω = 4π(22)0であるから、単位立体角に粒子が散乱する確率 p(Ω)dΩ はp(Ω)dΩ =dθ sin θ dω
dΩ
=
4π
2
2π

(23)

で与えられる。極角と方位角はそれぞれ独立であるから、

p(µ)dµ =p(θ)dθ

=

1
dµ
2
1
dθ
2π

と分離することができる。よって、乱数 ξ1 , ξ2 を用いて
∫ µ
dµ
=
−1 2
∫ ω
dω
=
0 2π

(24)
(25)

ξ1

(26)

ξ2 ,

(27)

µ

= 2ξ1 − 1(28)ω

= 2πξ2(29)から、極角と方位角を次のように決定することができる。

モンテカルロ計算ではこれまでに述べてきたような離散または連続確率密度関数からのサンプリングが随所に用いられる。これら確率密度関数は評価済み核データに基づいて決められるが、様々な関数からのサンプリングが必要となり、既に上で述べたように積算密度関数の逆関数を求めなければならない。積算密度関数の逆関数は上で述べたような簡単に求められるものから、非常に高度な数学的手法を必要とするものまで様々である。これまで開発されてきた積算密度関数の逆関数を求める方法は参考文献 [14] にまとめられており、実用上必要となるほぼすべての確率密度関数からのサンプリングが可能である。積算密度関数の逆関数が求めることができない確率密度関数は、棄却法 (rejection method)[15] によりサンプリングが可能であり、逆関数からの直接サンプリングよりも速い場合もある。
サンプリング法をコードへ実装する際に様々なサンプリング手法が工夫されている。例えば、ベクトル化な
どの高速化を計るために、discrete conditional sampling method[16] が提案されており、この手法では離散確率密度関数について、積算確率密度を用いずに確率密度関数から直接サンプリングすることができる。

2.3

誤差評価

モンテカルロ法では確率密度関数から事象や値をサンプリングし、物理現象を模擬してゆく。即ち、確率過程により物理現象を模擬するので、求めたいパラメータは平均値として評価される。それとともにパラメータの値は平均値の回りにばらつき、そのばらつきは統計誤差として評価される。モンテカルロ法では誤差評価は重要であり、この節では簡単に統計理論の用いられる定義と結果を引用し、モンテカルロ法の統計誤差につい
て考える。

2.3.1

期待値

X を確率変数とするとその期待値は次式で定義される。
∫ ∞E[X] =
xf (x)dx ≡ x̄,

(30)

−∞

ここで、f (x) は確率密度関数で、x̄ は確率変数 X の真の平均値とも呼ばれる。もし、f (x) から独立な n 個のサンプルをとり、それらを確率変数 X1 , X2 , · · · , Xn で表わすと、それぞれの確率変数の期待値はE[Xi ] = E[X] = x̄で、(30) 式の真の平均値と等しい。また、これらの確率変数の平均値 X =とができ、その期待値は

(31)
∑n

i=1 Xi /n も確率変数とみなすこ

[

]
n
n
1∑
1∑
E[X] = EXi =E [Xi ] = E[X] = x̄
n i=1
n i=1

(32)

で、X は真の平均値 x̄ の不偏推定量 (unbiased estimator) になっている。以上の式が示していることは、Xiや X の期待値は真の平均値 x̄ に等しいが、サンプリングした値は x̄ とはならないことを示している。つまり、サンプリングした n 個の平均値は x̄ の周りに拡がりを持って分布する。
そこで、必要となるのは拡がりを表わす尺度である。この尺度を導入するためにここでいくつかの式を定義しておく。確率変数 X の実関数 g(X) に対する期待値を次のように定義する。
∫ ∞E[g(X)] =g(x)f (x)dx ≡ ḡ

(33)

−∞

そのとき、(31) 式に対応して

E[g(Xi )] = E[g(X)] = ḡ

(34)

である。

2.4

分散

x̄ の周りの X の拡がりを評価するために (33) 式の g(X) に x̄ の周りの 2 次モーメントg(X) = (X − x̄)2を代入する。このとき、分散 (variance) は次のように定義される。
∫ ∞
[
]
2
2
σ (X) ≡ E (X − x̄) =(x − x̄)2 f (x)dx

(35)

(36)

−∞

(36) 式は良く知られた次式に書き直すことができる。
σ 2 (X) = x2 − x̄2(37)

この分散もしくは標準偏差 (standard deviation) はしばしば x̂ の周りの拡がりを表わす尺度として用いられる。
標準偏差は分散の平方根であり、次式で定義される。

√
σ(X) = x2 − x̄2(38)

さて、ここで確率変数 X の分散を確率変数 X の分散で表わすことを考える。

[]
σ 2 (X) = E (X − x̄)2 .

(39)

を定義すると、

σ 2 (X) =

=
=

(
)2 
n
∑
1
E
Xi − x̄ 
n i=1
(
)2 
n
∑
1
E
(Xi − x̄) 
n i=1
1 2
σ (X)
n

(40)

(41)
(42)

となり、従って、標準偏差は次式で与えられる。

σ(X)
σ(X) = √n(43)

この式はモンテカルロ計算における統計誤差の重要な性質を表わしている。n 個のサンプルを確率密度関数√
f (x) から取ってきて求めた平均の標準偏差は 1/ n の割合で減少していくことを示す。
期待値の場合と同様にして、真の分散 σ 2 (X) は知ることができない。そこで、不偏推定量から真の分散を推定することになる。この不偏推定量は次式
n
)
1 ∑(
U =
Xi − Xn − 1 i=12(44)で定義され、その期待値は[ ]E U 2 = σ 2 (X)(45)

である。実際には分散の不偏推定量は (44) 式からではなく次の式によって計算される。
)
n ( 2
2
U2 =
X −X
n−1
ここで、

(46)

1∑ 2
X
n i=1 i
n

X2 =

(47)

である。このとき標本標準偏差 (sample standard deviation) は次のようになる。
v
( n
)
u
∑
u n
1
2
U =t
Xi2 − X
n−1 n

(48)

i=1

2.5

中心極限定理

前節において分散の評価の仕方について述べた。しかし、分散はただサンプル平均値 x̂ =∑ni=1 xi やサンプ

ル xi の拡がりを示しているに過ぎない。ここでは、分散と推定された平均値の信頼度を結びつける関係式を導く。
定量的な信頼性を評価するために、よく知られている中心極限定理 (central limit theorem) を用いる。この定理によると、どのような組の有限の n 個のサンプルに対しても、サンプル平均値 x̂ の分布について記述する
確率分布関数 fn (x) が存在する。また、n が無限大に近づくにつれて x̂ には特別な極限分布が存在し、それは次のような正規分布となる。

[]
1−(x̂ − x̄)2
fn (x̂) ≈ √exp,n → ∞
2σ 2 (X)
2πσ(X)(49) 式は (43) 式を用いて、次のように書き直すこともできる。
√
[
]
n 1
−n(x̂ − x̄)2
fn (x̂) ≈exp,n → ∞
2π σ(X)2σ 2 (X)

(49)

(50)

この式を用いて、x̂ が x̄ −  と x̄ +  の間に入る確率を知ることができる。その確率を次式で定義する。
∫ x̄+
P {x̄ −  < x̂ < x̄ + } =fn (x̂)dx̂
(51)
x̄−

√

2/nt = (x̂ − x̄)/σ(X) と変換すると][√n 
P {x̄ −  < x̂ < x̄ + } = erf
2 σ(X)(51) 式に (50) 式を代入し、変数をとなる。ここで、erf(x) は誤差関数である。
さて、この  に適当な値を入れることにより、x̂ の拡がりに対応した信頼性を評価することができる。

σ(X) = σ(X) = √nσ(X) = 2σ(X) = 2 √nσ(X) = 3σ(X) = 3 √n→P = 0.683→

P = 0.954→

P = 0.997(52)

1st step

determine a neutron source position

2nd step

determine the flight direction

3rd step

determine the next collision point

escape fromthe reactor?

yes

no

4th step

determine the reaction

scatteredabsorbed or
scattered?

absorbed

図 4: 簡単な中性子ランダムウォークのシミュレーション

Fig. 4 Algorithm for the neutron random walk

中性子輸送のモンテカルロシミュレーション

3
3.1

はじめに

中性子の振る舞いはそれ自身ランダムで統計的なものであり、個々の中性子と物質の相互作用を確率と物理法則に従って忠実にシミュレーションすることにより、物質中での中性子の集合としての物理量を求めることができる。中性子輸送のシミュレーション、即ち個々の中性子のランダムウォークのシミュレーションは、前節で述べたサンプリング法を用いて、図 4に示すようなアルゴリズムに従って行うことができる。
最初のステップとして、中性子の発生点を決める。点源の場合はそこから発生させればよいし、空間的な拡がりをもつ中性子源の場合は乱数を用いて決定する。2 つ目のステップでは、飛行方向を (28) 式と (29) 式から決定する。3 つ目のステップでは、(19) 式もしくは (20) 式から次の衝突点を決定する。もし、飛行距離が体系の境界までの距離よりも長い場合は、中性子が体系から漏れたとみなし、追跡 (トラッキング) をやめ、最初のステップに戻り、次の中性子を追跡していく。もし、飛行距離が体系の境界までの距離よりも短い場合は、中性子が原子核と衝突したとして、(6) 式から (8) 式を用いてどの反応が起きたかを決定する。吸収反応が起きた
ときは中性子の追跡をやめ、最初のステップに戻る。散乱反応が起きたときは 2 つ目のステップに戻り、散乱の方向を決定する。
基本的には中性子輸送モンテカルロシミュレーションは上で述べたようなアルゴリズムに基づいているが、実際のモンテカルロコードでは幾何形状の取り扱い、エネルギーについての連続的な取り扱い、分散低減のための非アナログ的取り扱いがコード内でなされている。以下の節ではこれらについて簡単に説明していくこと
にする。

3.2

幾何形状

モンテカルロ計算は、基本的にある有限の空間領域内のある体系a について行われ、それぞれ異なる情報を必要とするいくつかの領域に分割される。これらの基本領域は、物質組成 (核種組成) が同じであるとか、その領域である物理量の評価量を求めたいとかの理由で分割される。その領域は直交座標系で表わされる任意の数の関数 f (x, y, z) = 0 で囲まれており、その各面は f (x, y, z) > 0 側と f (x, y, z) < 0 側とで 3 次元空間を二分割する。通常、この関数が表わす面の多くは 2 次関数

ax2 + by 2 + cz 2 + dxy + eyz + f zx + gx + hy + iz + j = 0(53)で表現されることがほとんどであるが、トーラスのように４次式になることもある。また、実用上、我々が対象とする体系を構成する領域は直方体、球、円柱等のように比較的単純な立体であることが多く、その場合、(53) 式の係数のいくつかは 0 となる。計算コードでは、簡略化のため、そのような単純な面は一般的な式から
ではなく、個別にとり使われるようになっていることが多い。
実際のモンテカルロコードでは領域を定義する入力方法は大まかに 2 通りに分けられる。一つは上に示したよ
うに面の方程式を直接定義する方法で、もう一つは組み合わせ幾何形状 (combinatorial geometry, CG) である。
組み合わせ幾何形状入力ではあらかじめコード内で定義されたボディと呼ばれる形状と交わり (intersection) と和 (union) のような集合演算を用いて領域を定義する。ボディは通常、球、直方体、円柱等、単純な形状があらかじめコード内で用意されている。前者の方法 (直接法) ではフレキシブルに領域を定義することができるが、3 次元空間内において関数 f (x, y, z) = 0 のどちらが正でどちらが負になるのか判定することが難しいこともある。後者の方法 (CG) は直感的な幾何形状の指定方法であるので、ユーザは簡単に計算体系を構成することができる。よく用いられているモンテカルロコードで分類すると MCNP は直接法を用いており、MVP、MORSE、
KENO は組み合わせ幾何形状を用いている。MCNP においても最近のバージョンでは、MACROBODY と呼ばれる入力で組み合わせ幾何形状のような入力ができるようになっている。逆に、MVP などにおいても、平面、無限円柱、一般 2 次曲面などの面を定義することもできる。
原子炉の炉心計算や臨界集合体の解析、また核物質貯蔵容器の配列などをモデル化する際には非常に多くの同じ領域を定義する必要が出てくる。もちろん、すべての領域を個別に定義しても良いが、入力データは複雑で膨大になる。また、計算速度の観点からそのような入力に対するモンテカルロ計算は効率的ではない。このような場合、1 つの単位領域b を定義し、繰り返し格子状に配置することにより全体の領域を定義することができる。このような入力法は繰り返し形状や格子形状と呼ばれ、ほとんどすべてのモンテカルロコードで利用可能である。
最後に、モンテカルロ計算体系の一例として、MVP コードにおける高温工学試験研究炉 (HTTR)[17] の全炉心を図 5(a) と (b) に示す。これらの図から分かるように、燃料ブロックの燃料棒から、減速材ブロックの制御棒挿入孔まで as built に模擬されており、このようなモンテカルロ計算により、空間の離散化による近似の入らない非常に精度の高い解を得ることができる。

3.3

トラッキング

一般に、モンテカルロ計算は、複雑な形状の体系中で輸送方程式を解く場合によく用いられる。この場合、領域や物質の境界を横切ってが自由飛行した後に衝突を起こす中性子も当然多い。モンテカルロ計算では、このような中性子を追跡 (トラッキング) し、どの面のどの点を横切り、どの領域に入ったかを逐次求めていかな
ければならない。そのアルゴリズムは単純であれ、複雑であれほとんど変わることはない。実際には、この中性子の追跡がモンテカルロ計算における計算時間を消費し、かつこの部分のプログラムは複雑になることが多い。ここでは中性子の追跡の手法について簡単に説明する。
図 6に示すような体系における中性子の自由飛行を考える。

1. まず、中性子は領域 R1 中の点 (x0 , y0 , z0 ) にあり、飛行方向 Ω であるとする。これは中性子源で発生したか、衝突がこの点で起こり、Ω 方向へ散乱されたと考えてよい。この中性子が次の衝突を起こすまで
a 無限については境界条件で取り扱う。
b セルと呼ばれることがある。

(a) 全炉心モデル(b) 拡大図
図 5: HTTR 全炉心計算モデル

R3

R2

R1

Ω

next
crossing
pointcrossing
point

u3

L
u2
(x0,y0,z0)

collision
point

u1

図 6: Tracking in a three-dimensional geometry.

の飛行距離 L を次式によって決定する。

L=−

1
ln ξ.
Σt

(54)

もし、中性子が物質の境界を越えて自由飛行するときは全断面積 Σ t が変化するので、一般に平均自由行
程 (mean free path) を単位として、飛行距離は次式で表わされる。

s = Σ t L = − ln ξ.

(55)

2. 飛行距離 L が決まると、領域 R1 において点 (x0 , y0 , z0 ) から Ω 方向へ直線を伸ばし、領域境界と交わる点を求める。一般には、飛行方向を表わす直線の方程式と領域 R1 を構成する面の方程式群 fi (x, y, z) =0, i = 1, 2, · · · , N との交点を求め、どの交点と点 (x0 , y0 , z0 ) の距離が最短でかつ +Ω にある点が境界を横切る点となる。

3. 飛行距離 L と境界までの距離 u1 を比べ、もし、L が u1 よりも小さければ中性子は同一領域内 R1 で衝突を起こす。もし、L が u1 よりも大きければ中性子は領域 R1 を出て隣の領域に入ったとして、その領域を同定する。領域 R2 が同定されると横切った点から Ω 方向の次の境界を横切る点までの距離 u2 を求める。再び、飛行距離 L と u1 + u2 を比較し、衝突が領域 R2 で起こるかどうかを判別する。このようにして中性子が衝突を起こすか、体系から漏れるまでこの過程を繰り返し、衝突点を決める。

3.4

衝突過程

中性子輸送のシミュレーションにおいては、中性子が自由飛行した後、媒質中のある原子核と衝突したとき、その入射中性子のエネルギーと飛行方向から散乱後のエネルギーと飛行方向を決める必要がある。ここでは、多群法の場合と連続エネルギー法の場合に分けてその決定方法について説明する。
多群法の場合、散乱後のエネルギーと飛行方向を決める方法は単純である。中性子の散乱前のエネルギー群
を g とすると、2.2節で述べた手法を用いて

Σ s,g→1 + · · · + Σ s,g→g0 −1
Σ s,g→1 + · · · + Σ s,g→g0
(56)
≤ξ<
Σ s,g
Σ s,g
∑
で決定される。ここで、Σ s は巨視的散乱断面積を表わしており、Σ s,g = g00 Σ s,g→g00 である。また、散乱角については、等方散乱の場合、(28) 式と (29) 式からµ= 2ξ1 − 1(57)ω

= 2πξ2(58)のように決定される。非等方散乱の場合は、通常群定数として低次のルジャンドル展開係数で実験室系での µの非等方性が与えられているので、これを用いて µ をサンプリングする。しかし、低次のルジャンドル展開形

式で与えられる散乱角分布は負の値をとる可能性があり、このような確率密度分布からサンプリングするには多少の工夫が必要である。古くからよく用いられるのは、離散角法 (discrete angle method)[18] で、µ についてのモーメントを保存するようにして、離散的な方向から µ をサンプリングする。その他には、散乱角分布を
等確率ステップ関数で表現する方法 [19] や、比較的最近では、非等確率ステップ関数で表現する方法 [20] が提案されている。
連続エネルギー法の場合、多群法では巨視的取り扱いしかできなかったのに対し、微視的な取り扱いが可能となる。散乱前のエネルギー E 0 、飛行方向 Ω0 の中性子がある点で衝突し、エネルギー E 、飛行方向 Ω へ散乱される確率は次式で表わされる。

p(E ← E 0 , Ω ← Ω0 ) =Σ s (E ← E 0 , Ω ← Ω0 )
Σ t (E 0 , Ω0 )(59) 式を構成核種と散乱のタイプで分けて表現すると次のようになる。
∑ Nk σt,k (E 0 ) ∑ σj,k (E 0 )f (E ← E 0 , Ω ← Ω0 )p(E ← E 0 , Ω ← Ω0 ) =0 ) j,kΣ t (E 0 )
σ
(E
t,k
j

(59)

(60)

k

ここで、k は核種を表わすインデックスで、j は弾性散乱、非弾性散乱等の散乱のタイプを表すインデックスである。また、Nk は核種 k の数密度、fj,k (E ← E 0 , Ω ← Ω0 ) は核種 k 、散乱タイプ j のエネルギー角度分布を表わしている。右辺最初の因子 Nk σt,k (E 0 )/Σ t (E 0 ) は k 番目の核種と衝突する確率を表わしており、この確率に従って衝突核種が決定される。右辺の第 2 因子 σj,k (E 0 )/σt,k (E 0 ) は散乱反応 j を起こす確率であり、この確率に従って散乱反応のタイプを決定する。エネルギー角度分布 fj,k は中性子と原子核の散乱の運動学(kinematics) から決まる場合もあるし、一般的には散乱タイプ毎に評価済み核データにテーブル形式や関数形式で与えられている。エネルギー角度分布から散乱後のエネルギーと角度を決定する方法は、次節にて中性子の反応毎に述べる。

3.5

連続エネルギー法における中性子の反応

3.5.1

弾性散乱

弾性散乱の場合、散乱後のエネルギーと角度は散乱の運動学から決定することができる。まず、重心系での散乱角余弦 µcm を決める必要があるが、重心系等方の場合、µcm は −1 ≤ µcm ≤ 1 で一様分布するので等確率密度分布から決定できる。一方、非等方散乱の場合には、解析的に求めるか、前もってテーブル化された分布から決定する。実際には、これらの元のデータは評価済み核データファイルに与えられているが、限られたエ
ネルギー点に対してのみ角度分布が与えられているので、それぞれのエネルギー点で角度分布の確率をテーブ
ルにしてモンテカルロコード用のライブラリに持たせておくのが一般的であり、かつ計算も速い。この場合、入射エネルギー E 0 とテーブルを与えているエネルギー点は通常異なるので、E 0 を挟む 2 点のエネルギー点に対し、用意されているテーブルを内挿することにより E 0 に対する値を求める。例えば、その 2 点のエネルギー点を E1 と E2 (E1 ≤ E 0 ≤ E2 ) とすればf (µcm , E 0 ) =E2 − E 0
E 0 − E1f (µcm , E1 ) +f (µcm , E2 )E2 − E1
E2 − E1(61)のような直線内挿式を用いる。

µcm が求められると、エネルギーと運動量保存則から散乱後のエネルギー E を次のように決定することができる。

E 0 (A2 + 2Aµcm + 1)
(62)
(A + 1)2
ここで、A は中性子に対する原子質量比である。また、µcm を実験室系に変換するためには次式を用いればE=よい。

1 + Aµcmµlab = √1 + 2Aµcm + A2
方位角については、散乱則に依存しないので
ω = 2πξ

(63)

(64)

により決定することができる。このようにして求めた µlab 、ω は、粒子の飛行方向に対する相対角度であるから、更に (x, y, z) 座標に対する散乱角余弦に変換する。

3.5.2非弾性散乱

多くの場合、非弾性散乱は重心系等方であり、この場合は µcm 、ω は (28) 式と (29) 式とからサンプリング
する。また、非等方の散乱角度分布が核データに与えられている場合は、弾性散乱と同様の方法でサンプリン
グする。このとき、散乱後のエネルギー E は次式で決まる。

E=

[
]
E0
1 + 2γ(E 0 )µcm + γ 2 (E 0 )2(A + 1)

ここで、√0

γ(E ) = A

1+

A+1 Q
A E0

(65)

(66)

で、Q は Q 値を表わし、負の値を取る。
また、実験室系での散乱角 Elab は次式で求められる。

1 + γ(E 0 )µcm
µlab = √1 + 2γ(E 0 )µcm + γ 2 (E 0 )(67)

これを用いて実験室系の散乱後のエネルギー E を表わすと次のようになる。

[
]2
√
E0
2
2
0
E=
µlab ± γ (E ) − 1 + µlab
(A + 1)2(68)

以上のように非弾性散乱の場合においても、ある入射エネルギーと入射角を持って非弾性散乱を起こすとある確率分布で散乱角が決まり、それに従って散乱後の中性子のエネルギーも一意的に決定できることが分かる。
しかし、以上の議論は離散レベル非弾性散乱に対し成立するものであり、核データの表現においては、すべての非弾性散乱が離散レベルとして表現されてはいない。分離できない離散レベルは連続レベル (continuum level)として表現されており、この場合の散乱エネルギー分布の与え方は、それぞれの核データによって異なるが、テーブル形式か蒸発モデルで与えることが多い。前者の場合は、与えられた確率テーブルから E を乱数で決定する。後者の場合は、次のような関数

f (E ← E 0 ) = CEe

E− T (E
0)(69)

で与えられ、

E = −T (E 0 ) ln(ξ1 ξ2 ),0 < E < E0 − U(70)

から E を決定する。ここで、T (E 0 ) は核温度であり、U は評価済み核データに与えられている定数である。

3.5.3(n,2n) 反応(n,2n) 反応のエネルギー分布は連続レベル非弾性散乱と同じ蒸発モデルで与えられていることが多く、この場合は前節で述べた方法で散乱後のエネルギーをサンプリングすることができる。角度については等方角度分
布であることがほとんどであるが、例外的に 2 つの中性子のエネルギー分布関数が与えられていることもある
ので、その場合にはそれを使ってサンプリングする。更に、(n,2n) 反応の場合には、放出された 2 つの中性子に相関を持たせずに衝突前と同じ重みを持たせてランダムに放出する場合と、放出する粒子を 1 つにして、重みを 2 倍にする方法がある。

3.5.4

核分裂反応

核分裂により放出される平均の中性子数 ν̄ はテーブル形式または解析式により与えられており、その分布か
ら乱数を使って発生中性子数を決定する。発生する中性子数は常に整数個でなければならないので、発生する
核分裂中性子数の期待値が ν̄ となるように確率的に発生数は決められる。
核分裂中性子の角度分布はほとんどの場合等方であり、この場合は (28) 式と (29) 式から発生した中性子の飛行方向を決定することができる。最近の核データでは角度分布がテーブル形式で与えられていることもあるが、この場合でもその確率テーブルから乱数を用いて中性子の飛行方向を決定することができる。

核分裂中性子のエネルギー分布は、核データにテーブル形式やマクスウェル分布で与えられている。マクス
ウェル分布の場合は 3 つの乱数 ξ1 、ξ2 、ξ3 を用いて次式から直接エネルギーを決定することができる。
[( π )]
E = −R(E 0 ) ln ξ1 + ln ξ2 cos2
ξ3
(71)
2
ここで、R(E 0 ) は入射エネルギーに依存する定数である。

3.6

熱中性子散乱熱中性子散乱を扱うモデルもいくつかあるが、化学結合の効果や原子核の熱運動の影響も取り入れた散乱則

S(α, β) に基づくのが最も厳密である。通常、モンテカルロコードでは散乱則データから直接サンプリングす
るのではなく、それを基にエネルギー遷移確率と対応する確率分布テーブルを作成し、そこから散乱後のエネ
ルギーと角度を決定する。
散乱則データが与えられていない場合は、自由ガスモデル (free gas model) と呼ばれる方法で熱中性子散乱を取り扱う。このモデルでは物質を構成する原子核は単一原子のガスとして存在し、運動は等方的で、速度についてマクスウェル分布していると仮定する。このとき、標的原子核の速度 (速さと飛行方向) を確率分布から決定する。その後、二体弾性衝突の運動学 (kinematics) を用いて、中性子の散乱後のエネルギーと方向を決定
する。

3.7

アナログモンテカルロ法と非アナログモンテカルロ法

モンテカルロ法で中性子を追跡するとき、2 つの手法がある。１つは、物質中での中性子の振る舞いを現実と全く同じ方法で追跡する方法である。この方法では、中性子が原子核と衝突したときに起こった反応を乱数で決定し、それが吸収であれば、そのヒストリーの追跡は終了し、中性子源から新たに中性子を発生させる。
反応が吸収でなければ追跡し続け、吸収されるか体系から漏れるまでその中性子を追跡する。この方法をアナログモンテカルロ法 (analog Monte Carlo) と呼ぶ。この方法の欠点は、問題によっては分散が極めて大きくな
り、信頼できる解が得られないことである。例えば、吸収確率の大きな体系での透過確率を求める問題や、共鳴領域の吸収が大きい媒質に対し、それ以下のエネルギー領域での情報を得たい場合などである。
このような欠点を克服するために、多くのモンテカルロコードでは、非アナログモンテカルロ法 (non-analog
Monte Carlo) を採用しているc 。非アナログモンテカルロ法では、重み (weight) という概念を導入し、中性子の存在割合をこれで表わす。アナログモンテカルロ法では粒子は存在するかしないか、即ち、重みは 1 か 0 かしかないが、非アナログモンテカルロ法ではその中間の状態を重みによって記述するのである。従って、この方法においては、衝突において吸収という事象は起こさず、その代わりに吸収確率だけ中性子の重みを小さくする。即ち、W の重みを持った中性子が衝突を起こすと、衝突後の重みは W (1 − Σ a /Σ t ) に変えられる。係数1 − Σ a /Σ t は生き残り確率 (survival probability) に対応している。このようにすると、中性子の数は減ることがなく、後に述べるエスティメータ (estimator) に寄与する機会も多くなり、分散を小さくすることができる。
非アナログモンテカルロ法で中性子を追跡すると、体系から漏れる以外は追跡が終了しないので、計算時間は増加することになる。特に大きな体系を解いている場合は、中性子が漏れる確率も小さくなるので、極めて小さな重みしか持たない粒子をいつまでも追跡していることになり、計算時間の観点から効率の悪い計算になってしまう。このように分散を小さくすることと、計算時間を短くすることは相反することになるので、ある種の尺度を用いて判断することが多い。モンテカルロ計算においては次式で定義される FOM が用いられることが多く、これを最大にする方法が最適と判断できる。

FOM =

1
σ2 t

(72)

ここで、σ 2 は分散、t は計算時間を表わす。そこで、非アナログモンテカルロ法では、重みの小さくなった中性子の追跡を終了させるいくつかの技法が導入されている。それについては以下の節で述べる。

c コードによっては、オプションによりアナログと非アナログモンテカルロ法を使い分けることができる。

エスティメータ

4

これまで、中性子のランダムウォークを模擬する手法を説明してきたが、ただ模擬するだけでは中性子束や反応率といった物理量を求めることはできない。モンテカルロ法でこれらの物理量を計算する手段のことをエスティメータ (estimator)d と呼ぶ。具体的には、ランダムウォークの過程で、衝突、自由飛行、漏れなどの事象が起こるごとにエスティメータへの寄与をスコアリングし、それらを蓄積していき、統計処理をすることにより物理量を評価する。また、エスティメータへの寄与をスコアリングすることを一般にタリー (tally) と呼ぶ。
ここでは、中性子に対するエスティメータをいくつか説明していくことにする。反応率は中性子束に反応断面積を係数としてかけてやることにより同様に評価することができる。

4.1

衝突エスティメータ

衝突エスティメータ (collision estimator) では平均の衝突回数から中性子束を評価する。ある体積 V における全断面積が Σ t で、その体積の平均中性子束が φ̄、平均の衝突回数が c̄ であったとすると、次式が成り立つ。

c̄ = Σ t φ̄V

(73)

1
c̄
ΣtV

(74)

この式は次のように書きなおすことができ、φ̄ =衝突エスティメータはこの式を根拠に中性子束を評価する。非アナログモンテカルロ法では、中性子束に対する衝突当りの寄与は

W
ΣtV

(75)

となる。ここで、W は衝突前の中性子の重みである。
衝突エスティメータは衝突が起こらないと中性子束を評価することができないという欠点がある。それゆえ、ボイド領域は媒質が希薄な体系では有効ではない。

4.2

飛程長エスティメータ

中性子束はすべての中性子が横切る単位体積当り全飛程長 (track length) としても定義される。即ち、

φ̄ =

1¯
`
V

(76)

のように定義される。ここで、`¯ は平均の飛程長である。これに基づき、飛程長エスティメータ (track length

estimator) ではある領域を横切る中性子の飛程長をスコアリングすることにより中性子を評価する。非アナロ
グモンテカルロ法では、中性子束に対する各飛程からの寄与は
`
W
V

(77)

飛程長エスティメータではボイド領域の中性子束も評価でき、領域がある程度以上の大きさをもてば、信頼度の高い結果が得られる。

4.3

面交差エスティメータ

面交差エスティメータ (surface crossing estimator) は面平均の中性子束を評価するもので、ある面の領域を中性子が通過するたびに次の値をスコアリングする。

W
A |n · Ω|(78)

この評価法では、中性子の飛行方向が面にほぼ平行に入射したとき (n · Ω → 0) には、寄与が発散する問題が生じる。このとき、各コードにおいては、|n · Ω| に対してある下限値を設定し、それ以下になると下限値に設定するようにしていることが多い。
d 日本語では評価法とも呼ばれるようであるが、あまりなじみがないのでここではエスティメータと記述することにする。

このエスティメータは体系からの漏れなど中性子流れ (current) の評価にも用いられる。この場合には分母の n · Ω がなくなるので、信頼性の高い結果が得られる。

4.4

点検出器エスティメータ

炉心計算や遮蔽計算の解析ではある点における中性子束を求めたいことがある。しかし、ある点を中性子が横切る可能性はほとんど 0 であるので、衝突エスティメータや飛程長エスティメータのように中性子が飛んでくるのを待つような方法では、点における中性子束を評価することができない。そこで、点検出器エスティメー
タ (point detector estimator) では、図 7に示すように中性子源と散乱点から仮想的な粒子を直線的に飛ばし、評価したい点への寄与を評価することにより点における中性子束を評価する。具体的には次のように中性子源
からの直接寄与と散乱源からの寄与とに分けて評価する。

Φ(r d ) =∑i+

Wi Pi (θ, φ)∑1
exp [−β(r i , r d )]|r i − r d |2

Wj Pj (θ, φ)jΣ s (r j )1
exp [−β(r j , r d )] ,Σ t (r j ) |r j − r d |2(79)

ここで、Wi は中性子源での重み、Wj は j 回目の衝突前の重み、P は線源または衝突した中性子が方向 (θ, φ)の周りの立体角 dΩ に入る確率、β は検出器点 r d と線源/衝突点 r i /r j の光学距離である。

Real particle path

Pseudo-particle path

Collision point
Collision point

SourceContribution
from source
W 0 P0 ( θ , φ )r1r0e − β0
r02

r2

Contribution fromcollision
W 2 P 2 (θ , φ )qDetector pointΣ s e − β2
Σ t r 22

図 7: 点検出器の概念図(79) 式から分かるように、点検出器エスティメータは 1/r2 の特異点を持ち、r j = r d のときは発散してしまい、無限大の分散を持つことになる。これを避けるために、各コードでは検出器からある距離の範囲では中性子束を一定と仮定して、検出器近傍における散乱からの寄与を計算するなどの近似的手法が用いられる。従って、散乱媒質中で点検出器エスティメータを用いるときは注意が必要である。

5

実効増倍率の評価
炉心解析において実効増倍率 keﬀ は最も重要なパラメータの一つである。実効増倍率は固有値問題の固有値に相当するものであり、固有関数に対応する核分裂源分布もあらかじめ分かっていない。このような固有値問
題を決定論的手法で解く場合は源反復を用いるが、モンテカルロ法においても世代 (generation) の概念を導入
し、べき乗法 (power iteration method) を用いる。モンテカルロ法におけるべき乗法では世代をバッチまたはサイクルと呼び、各バッチ (サイクル) において同じ数の中性子を発生させるか、または中性子の総重みが同じになるように発生させる。最初のバッチでは核分裂源分布が分かっていないので、初期分布 (initial guess) として点源か平坦分布を仮定し、中性子を追跡する。追跡過程において中性子が衝突した点が次の世代の中性子の出発点となり、そこから次の世代の中性子を追跡している。これを何バッチも繰り返し、核分裂源分布が基本モードに収束したと思われるバッチ以降の情報を処理することにより、実効増倍率をはじめとする中性子束、反応率などを評価する。
実効増倍率の評価法について述べる。実効増倍率はある世代に生まれる核分裂中性子数と次の世代に生まれ
る核分裂中性子数の比と定義される。即ち、

kj =

Fj
Fj−1

(80)

である。ここで、kj は第 j 世代の実効増倍率、Fj は第 j 世代に核分裂で生まれる中性子数 (総重み) を表わす。

(80) に従って実効増倍率を評価するには、基本的に以下の 3 つの方法がある。
1. 衝突エスティメータ
kj =

[
]
1 ∑ ∑ νΣ f
Wn,i
N n i
Σ t n,i(81)

ここで、N は 1 バッチ当りの源中性子の総重み、n はヒストリーを表わすインデックス、i は衝突を表わすインデックス、Wn,i はヒストリー n の i 回目の衝突前の重みを表わす。

2. アナログまたは吸収エスティメータ
kj =

[
]
1 ∑ ∑ νσf,kW0
N n i
σa,k n,i n,i(82)

Σc + Σf
Wn,i
Σt

(83)

0
Wn,i
=

ここで、k は衝突核種、Σ c 、Σ f はそれぞれ巨視的捕獲、核分裂断面積である。アナログエスティメータと衝突エスティメータの違いは、衝突核種だけの寄与をカウントするか構成核種すべての寄与を平均的にカウントするかである。

3. 飛程長エスティメータ
kj =

1 ∑∑
[νΣ f `]n,i Wn,i
N n i

(84)

ここで、` は飛程長、i は飛行を表わすインデックスである。
上では実効増倍率を (80) 式で定義したが、実効増倍率を

kj =

Fj
Aj + Lj(85)と定義し、これに基づいて、上と同様の衝突、アナログ、飛程長エスティメータにより実効増倍率を評価することもできる。ここで、Aj は第 j 世代に吸収でなくなる中性子数 (総重み)、Lj は第 j 世代に漏れてなくなる
中性子数 (総重み) を表わす。
体系の実効増倍率は、上記のようにして得られた各バッチの実効増倍率 kj を用いて次のように評価することができる。

k̄ =1
M − Mskip

M
∑

kj

(86)

j=Mskip +1

ここで、Mskip は統計処理から除くバッチ数、M は総バッチである。その標準偏差は

√
σk̄

=

k2

=

(
)
1
k 2 − k̄ 2
M − Mskip − 1

1
M − Mskip

M
∑

kj2

(87)

(88)

j=Mskip +1

以上のように複数個の方法で実効増倍率を評価するのは、体系によって有利な評価法が異なるためである。
更に、一般には 2 つ以上の評価法によって得た値の平均をとると、その分散は各々の分散より小さい。平均のとり方は、単純な算術平均や相関度を考慮して共分散を重みとして行う場合がある。

6

分散低減法

モンテカルロ計算の精度を上げるには、ヒストリー数 N を増加させなければならないが、統計誤差 (標準偏√差) は 1/ N でしか減少しない。同じ計算時間でより精度を上げる (つまり、FOM を大きくする) 方法が、いろ
いろ工夫されている。ここでは、一般的なコードでよく用いられる分散低減法 (variance reduction technique)についていくつか説明する。一般には、以下に述べる分散低減法を効率よく使うにはかなりの経験を必要とす
るが、実際には何回か試行を行うことに分散低減法のパラメータを決めることがよく行われる。

1. ロシアンルーレット (Russian roulette)
ロシアンルーレットは最も一般的で有効な方法である。ある中性子が位相空間のある領域境界を横切ったり、原子核と衝突したりしたとき、指定された重みに従った確率で粒子を消滅 (kill) させる。そのとき、粒子の重みはバイアスがかからないように調整される。ロシアンルーレットは非アナログモンテカ
ルロ法で粒子のヒストリーを終了させるのにしばしば用いられ、最も単純な使用はウェイトカットオフ(weight cutoﬀ) で用いられる場合である。ウェイトカットオフでは重みを W の中性子がある指定した重み以下になったとき、その中性子にロシアンルーレットを適用する。このとき、生き残り確率

p=

W
WSRV

(89)

で中性子を生き残らせ、確率 1 − p で中性子を消滅させる。生き残った中性子には重み WSRV を与える。
このようにすると、

pWSRV + (1 − p) × 0 = W(90)

で、元の中性子の重みは保存されていることになる。

2. スプリッティング (splitting)
スプリッティングは中性子がある位相空間を横切ったとき、同じ重みを持つ ν 個の中性子に分割する方法である。スプリッティングは重要な (ユーザが情報を得たい) 領域の粒子数を増やす方法で、通常、ロシアンルーレットと組み合わせて用いられる。

3. インポータンス (importance)
インポータンスは形状スプリッティング/ロシアンルーレット (geometrical splitting/ Russian roulette)とも呼ばれ、各領域に対してインポータンスを設定し、インポータンスに従ってロシアンルーレットとスプリッティングを行う方法である。中性子がインポータンス In の領域 n からインポータンス In+1 の領域 n + 1 に入ったとする。もし、In+1 /In ≥ 1 なら、スプリッティングを行い、ν = In+1 /In 個の中性子に分割する。もし、In+1 /In < 1 なら、ロシアンルーレットを行い、確率 ν で中性子を生き残らせる。
生き残った中性子には重み W In /In+1 を与える。インポータンスは一般的で分散を低減させる有力な方法であるが、最適なインポータンスを見つけることが容易でないことが多い。

4. ウェイトウインドウ (weight window)
ウェイトウインドウは、エネルギー、空間依存のスプリッティング/ロシアンルーレットと呼べるものである。この手法では、上限の重み WU と下限の重み WL を決めておき、衝突後の粒子の重み W が W > WUの時はスプリッティングを行い、W > WL のときはロシアンルーレットを行う。これにより中性子の重みは WU と WL の間に揃うことになる。分散は各中性子の重みが揃っているほど小さくなることから、分散の低減に有効な手法である。

7

モンテカルロ計算に潜む不確かさ

この節では、モンテカルロ計算に潜む不確かさ (uncertainty) を考える。ここで言う不確かさとは、推定誤差(estimated error) のことであり、特に誤解のないときには、単に誤差 (error) と呼ばれることが多い。
一般論として、不確かさには大きく分けて、偶然誤差 (random error) と系統誤差 (systematic error) の 2 種類がある [21]。(正確には、偶然不確かさ (random uncertainty) と系統不確かさ (systematic uncertainty) と呼ぶべきであるが、そのような日本語は聞いたことがないので、ここでは誤差と呼ぶことにする。) 偶然誤差は、
読んで字のごとく、偶然に起因したばらつきの誤差である。現象が本質的にランダムなものである場合、統計的手法により推定することができるので、統計誤差 (statistical error) と呼ばれる。系統誤差は、偶然によらない、一定の傾向を持った誤差で、バイアス (bias) とも呼ばれる。系統誤差は、統計的に扱えるものではないが、多くの場合、真値の確率分布を与えることによって考慮される。以下では、モンテカルロ計算に潜む統計誤差とバイアスについて述べる。

統計誤差

7.1

モンテカルロ計算でまず第一に考量しなければならないのが統計誤差である。これについては、既に 2.3 節で述べた。統計誤差を評価した例として、小林によって提案された 3 次元輸送計算ベンチマーク問題 [22] をモン
テカルロコードで計算した結果を示す。図 8は、小林ベンチマーク問題 1 の計算体系であり、3 つの立方体領域で構成される。xy 平面、yz 平面、zx 平面は反射境界条件で取り扱われ、それ以外の外側境界面は真空境界条
件として取り扱われる。左下の 10 cm×10 cm×10 cm の領域は中性子源領域で、それを 50 cm×50 cm×50 cmのボイド領域が取り囲んでおり、更に 100 cm×100 cm×100 cm の遮蔽体領域が取り囲んでいる。

Flux (1/sec/cm^2)z (cm)
100y (cm)
100
50

Void

0
0

Source

100

200
300
Number of bacthes

Standard deviation

Shield

50

0.080
0.075
0.070
0.065
0.060
0.055
0.050
0.045
0.0400
0.014
0.012
0.010
0.008
0.006
0.004
0.002
0.000050

100

x (cm)

図 8: 小林ベンチマーク問題 1 の体系

100

200
300
Number of bacthes

400

500

Standard deviation
Fit with a/sqrt(n)

400

500

図 9: 飛程長エスティメータで評価した (x, y, z)=(45 cm,

45 cm, 45 cm) における中性子束

図 9は、GMVP コード [11] で (x, y, z)=(45 cm, 45 cm, 45 cm) における中性子束を飛程長エスティメータで計算した結果である。1 バッチ当り 50,000 ヒストリーで計算した。解析的に求めた参照解は、0.0525132 cm−2 sec−1
であり、バッチ数 (ヒストリー数) が増加するにつれて GMVP で計算した結果は参照解に近づいていくことが分かる。図 9の下の図は、計算値の統計誤差 (1 標準偏差) のバッチ数に対する依存性を示している。点線は計√算で評価された統計誤差であり、実線はそれを a/ n(n:バッチ数) でフィッティングしたものである。2.3 節で述べたように、統計誤差は、ヒストリ数の平方根の逆数に比例して減少していくことが分かる。
図 10は、小林ベンチマーク問題 1 に対し、点検出器エスティメータで中性子束を評価した結果である。(a)は、中性子源領域 ((x, y, z)=(5 cm, 5 cm, 5 cm)) における中性子束を評価したものであり、80 バッチ目と 180
バッチ目において、中性子束と標準偏差が突然大きくなる現象が見られる。これは点検出器エスティメータに特有の現象であり、4.4 節で述べたように、点検出器エスティメータが 1/r2 の特異点を持つために生じている。
即ち、検出器位置の近傍で散乱 (もしくは中性子の発生) が起きたために、非常に大きなエスティメータへの寄与があったということである。このような場合、統計誤差は、ヒストリ数の平方根の逆数に比例して減少していくことはなく、統計誤差の大きさ自体も信用できないかもしれない。検出器近傍からの寄与が十分に考慮されていない場合には、中性子束とその統計誤差にバイアスがかかるので、点検出器エスティメータを使い際には注意が必要である。MCNP マニュアル [10] では、点検出器エスティメータに対するガイドラインとして、相対誤差が 5%以下であれば、一般的に評価が信用できるとしている。
図 10(b) は、ボイド領域 ((x, y, z)=(45 cm, 45 cm, 45 cm)) における中性子束を点検出器エスティメータで

0.25

100

200
300
Number of bacthes

400

500

100

200
300
Number of bacthes

400

500

0.000120.000100.00008
0.00006
0.00004
0.000020.000000

100

200
300
Number of bacthes

400

500

100

200
300
Number of bacthes

400

500

Standard deviation

Standard deviation

0.20
0.15

0.10
0.05

0.000

0.05260
0.05255
0.05250
0.052450.052400.05235
0.05230
0.052250

Flux (1/cm^2/sec)

Flux (1/cm^2/sec)

6.20
6.10
6.00
5.90
5.80
5.70
5.600

(a) (x, y, z)=(5 cm, 5 cm, 5 cm) における中性子束(b) (x, y, z)=(45 cm, 45 cm, 45 cm) における中性子束

図 10: 小林ベンチマーク問題 1 に対して点検出器エスティメータで中性子束を評価した結果評価した結果である。(a) の結果と異なり、検出器近傍からの寄与はないため、大きな寄与が発生して、中性子束とその統計誤差が突然大きくなることはない。このように検出器近傍からの寄与がないことがあらかじめ分かっている場合には、点検出器エスティメータを用いて、容易に信頼できる中性子束や反応率を計算することが可能である。

7.2

固有値問題における核分裂源の規格化によるバイアス

5 節で述べたように、粒子輸送モンテカルロ法では、実効増倍率 (固有値) をべき乗法により計算し、ある世代に発生した核分裂中性子の情報 (発生位置など) をすべて保存しておき、それを次の世代の中性子源 (スターター) として中性子を発生させる。ほとんどすべての場合において、実効増倍率がちょうど 1 となることはな
いので、すべての核分裂中性子を次の世代のスターターとすると、中性子の数は発散してしまう (keﬀ > 1) か、

0 になってしまう (keﬀ < 1) ことになる。そのため、1 世代当りのヒストリー数 (バッチサイズ) もしくは 1 世
代当りの総重みを固定し、べき乗反復を行っている。即ち、核分裂中性子数 (総重み) は、世代毎に規格化されることになり、これによりバイアスが生じる。このバイアスについての理論的研究は既に過去行われており[23, 24, 29]、Brissenden と Garlick によると、keﬀ のバイアス ∆k は、バッチサイズの逆数 1/N に比例することが示されている。
1
∆k ∝
(91)
N
(91) 式を数値的に検証するために、バッチサイズの値を変えて、Godiva 炉心に対する実効増倍率を計算した結果が、図 11である。MVP コードで JENDL-4.0 を用い、バッチサイズを 10,000, 1,000, 100, 50, 25, 20 と変化させて計算した。すべての計算において、スキップバッチ数は 100 とし、統計処理に用いる総ヒストリー
数は一定 (100 万ヒストリー) となるようにした。即ち、バッチサイズ 10,000 の場合は、有効バッチ数 100、ス
キップバッチ数 100 の計算を行い、バッチサイズ 1,000 の場合は、有効バッチ数 1,000、スキップバッチ数 100の計算を行った。この図で示されるように、バッチサイズが小さくなるほど、実効増倍率の評価値に対するバ
イアスは直線的に大きくなることが分かる。
図 11から、バッチサイズを 1,000 以上にしておけば、世代間の規格化によって生じるバイアスは非常に小さいことが分かる。最近の計算機環境では、通常、バッチサイズは 10,000 以上にすることが多く、ベクトル計算
や並列計算では、計算効率の点からもっと大きいバッチサイズを用いるため、このバイアスはほとんど無視することができる。

7.3

固有値問題における統計誤差のバイアス

5 節と前節で述べたように、粒子輸送モンテカルロ法で固有値問題を解く場合 (実効増倍率を計算する場合)には、べき乗法を用いる。即ち、ある世代に発生した核分裂中性子の情報 (発生位置など) をすべて保存してお0.9990.998 N=10000
N=1000
0.997

N=100

0.996
keff

N=50

0.995
0.994

N=25

0.993

N=20

0.992
0.00

0.01

0.02

0.03
1/N

0.04

0.05

0.06

図 11: 実効増倍率のバッチサイズ N に対する依存性

き、それを次の世代の中性子源 (スターター) として中性子を発生させる。これにより、世代間に相関が生じ、実効増倍率などの評価する量の統計誤差を正確に求めるにはこの相関を考慮する必要がある。例えば、実効増倍率に対する正確な分散 (真の分散) は次式 [25] で定義され、この式の右辺第 2 項が世代間の相関から生じる分散を表している。
2=
σR

M −1 M12 ∑ ∑cov [ki , ki ] + 2
cov [ki , kj ]M
M i=1 j=i+1(92)2
ここで、σRは実効増倍率に対する真の分散、M は統計処理に用いた (アクティブ) バッチ数、ki は第 i バッ

チの増倍率で、cov [ki , kj ] は ki と kj の共分散である。右辺第 1 項は、第 i 世代の増倍率に対する分散を表しており、平均値と同様に二乗平均値を計算すれば容易に評価できるので、すべてのモンテカルロコードにおいて世代間の相関を無視した第一次近似として評価されているものである。この近似は、以前より統計誤差に対してバイアスを生じさせることが知られており、正確な統計誤差を推定する手法もいくつか提案されている[26, 27, 28, 29]。しかしながら、世代間の相関から生じる分散を評価するのは難しく、現在でも様々な手法が提案されている。最近では、時系列解析を用い、より洗練された真の分散推定法が提案されている [30, 31, 32, 33]。
これらの手法の詳細については、論文を参照していただくことにして、ここでは、実際に真の統計誤差と世代間の相関を無視して評価した統計誤差の間にどれくらいの差を生じるか文献調査の結果を示すことにする。真の分散推定の研究のため、世代間の相関が大きくなるような極端な体系の結果は除き、なるべく現実的な体系について評価しているものを選択した。
植木らの実効増倍率に対する計算例

植木らは、17×17 燃料集合体の PWR を簡略化した 3 領域の円柱体系

2(内側炉心、外側炉心、水反射体) について、多群モンテカルロ計算で実効増倍率に対する真の分散 σRと世代2間の相関を無視した見かけの分散 σAを評価している [25]。真の分散は、乱数列を変えて実行した独立なモン

テカルロ計算を 500 回実行することによって評価している。表 1は、彼らが評価した統計誤差のバイアスをまとめたものである。実効増倍率に対する統計誤差は、世代間の相関を無視することにより 15%程度過小評価されていることが分かる。

Shim らの計算例

Shim らは、小型の PWR 炉心と中型の PWR 炉心に対する実効増倍率と集合体出力の統計

誤差のバイアスを連続エネルギーモンテカルロコード McCARD によって評価している [34]。真の統計誤差 σR
は、植木らの例と同様に乱数列を変えて実行した 500 回の独立したモンテカルロ計算から求めている。表 2は、彼らが評価した実効増倍率に対する統計誤差のバイアスをまとめたものである。Shim らの計算例においても、実効増倍率に対する統計誤差は、世代間の相関を無視することにより 15%程度過小評価されることが分かる。

表 1: 簡略化 PWR 炉心 3 次元円柱体系に対する統計誤差のバイアス
見かけの分散と真の分散の比

標準偏差の比

2
2
(σA
/σR
)

(σA /σR )4 群飛程長エスティメータ
4 群衝突エスティメータ

0.690
0.676

0.830
0.822

1 群飛程長エスティメータ
1 群衝突エスティメータ

0.714
0.690

0.845
0.830

表 2: 小型と中型 PWR 炉心の実効増倍率に対する統 表 3: 小型と中型 PWR 炉心の集合体出力に対する統計誤差のバイアス

計誤差のバイアス
標準偏差の比標準偏差の比(σR /σA )(σR /σA )小型 PWR

0.833

小型 PWR(集合体 11 体の平均)

0.685

中型 PWR

0.858

中型 PWR(集合体 26 体の平均)

0.802

彼らはまた、集合体出力に対する統計誤差のバイアスも評価しており、集合体出力に対する統計誤差は、世代間の相関を無視することにより 20%-30% 過小評価されることが分かる (表 3)。
植木の出力分布に対する計算例

植木は、1,100 MWe 級 PWR 全炉心体系について、出力分布に対する標準偏差のバイアスを連続エネルギーモンテカルロコード MCNP を用いて評価している [33]。真の分散は、乱数列を変えて実行した独立なモンテカルロ計算を 25 回実行することによって評価している。表 4は、彼が評価した統計誤差のバイアスをまとめたものである。出力分布の代表値として、集合体の中央ノード部の出力を計算し、炉心中心集合体 (タリーセル 0)、x 軸上集合体 (タリーセル 2, 5)、非対称位置集合体 (タリーセル 1, 3, 4)
における出力をモンテカルロ計算により計算している。この表より分かるように、世代間の相関を無視することにより、60%-80%も統計誤差は過小評価され、対称な位置のセルより、非対称な位置のセルの方が大きく過小評価されていることが分かる。この原因については、核分裂源分布の固有関数の非対称性から理論的に説明できることが指摘されている [35]。
表 4: 1,100 MWe 級 PWR 全炉心体系における出力分布に対する統計誤差のバイアス
標準偏差の比(σA /σR )炉心中心集合体 (タリーセル 0)0.356

x 軸上集合体 (タリーセル 2)
x 軸上集合体 (タリーセル 5)
非対称位置集合体 (タリーセル 1)

0.377
0.318
0.213

非対称位置集合体 (タリーセル 3)

0.225
0.250

非対称位置集合体 (タリーセル 4)7.4

固有値問題における収束判定

固有値問題では、実効増倍率 (固有値) とともに核分裂源分布 (固有関数) も未知であり、これらを同時に求める必要がある。これまで既に、粒子輸送モンテカルロ法で固有値問題を解く場合、実効増倍率を求めるためにべき乗法という反復法が用いられることを述べたが、核分裂源分布については最初から分かっているわけではないので、最初はある分布を仮定し、反復毎に求めるべき定常状態の核分裂源分布に収束させていく。このとき、収束していない核分裂源分布から評価された値をサンプルとして用いると平均値の評価にバイアスがかかる可能性があるので、定常状態に収束していない核分裂源分布から評価された平均値のサンプル値は統計処
理から除外する必要がある。ほとんどのモンテカルロコードでは、核分裂源分布が収束するまでの世代数を統計処理から除外する世代数e として入力値として与えてやる必要があり、ユーザが何世代で核分裂源分布が収束するか前もって判断しておかなければならない。
統計処理から除外する世代数は、安全のため大きいほうが無難であるが、大きい値に設定すると計算時間もかかり、適切な値を設定することが求められる。最近のコンピュータは処理も早いので、この世代数は 100 程度とされることが多いように思われる。この値は、過去の実効増倍率に対する計算結果からの工学的判断で、核設計における一般的な問題に対して十分であると考えられる値である。しかしながら、臨界安全分野では、結合の弱い体系などを扱うことが多く、核分裂源分布の収束が非常に遅くなることがあることが指摘されている [36]。最近では、OECD/NEA の臨界安全に関するワーキングパーティ(WPNCS) で、モンテカルロ計算の核分裂源収束に関する専門家グループ (expert group) が立ち上げられ、収束に関する研究が精力的に行われている [37]。
核分裂源分布の収束が遅い例として、植木らによって提案された結合の弱い体系に対するベンチマーク問題[38] を取り上げる。体系は、図 12に示されるような 1 次元平板状体系で、2 つの燃料領域があり、その間に吸収体があるような体系である。燃料領域の幅は、左が 1 cm で右が 1.01 cm で微妙に異なっており、核分裂源
分布は非対象になることが予想される。この体系について、GMVP コードを用いて、実効増倍率と左右の燃料領域における核分裂率を計算した。各世代のヒストリー数 (バッチサイズ) は 5,000 とし、初期中性子源分布は全領域で一様とした。図 2 領域の燃料領域で構成される結合の弱い平板上体系における収束は、横軸にバッ
チ数をとり、実効増倍率と核分裂源分布の収束の様子を示したものである。縦軸は、参考文献 [38] にある参照
解 (実効増倍率 0.427425, 右側核分裂率割合 0.9744) で規格化した値を示している。この図から分かるように、実効増倍率は早く収束するが、核分裂源分布は 400 世代以降にならないと収束しない。従って、実効増倍率の傾向だけで核分裂源分布の収束を判断することは難しく、大きな体系で出力分布などを計算するときは、注意が必要である。

1 cm

1 cm

1.1
Normalize values divided by true mean

VacuumV

1.01 cm

Scatterer
Fuel
F

5 cm

Ab
bsorberF
FuelScatterer

Vacuum
V

1 cm

1.0
0.9
0.8
0.7
0.6
0.5

k-eff
source at right fuel

0.4
0.30

200

400
600
Number of bacthes

800

1000

図 12: 2 領域の燃料領域で構成される結合の弱い平板

図 13: 2 領域の燃料領域で構成される結合の弱い平板

上体系

上体系における収束

これまでの研究から、核分裂源分布の収束の早さは、ドミナンス比 (dominance ratio, 基本モード固有値と 1
次モード固有値の比) と密接に関連しており、ドミナンス比が 1 に近いほど収束が遅くなることが知られている。上の例では、ドミナンス比は、0.9925 と計算されている [38]。ドミナンス比を前もって知ることは難しいが、それが分かれば収束の速さを予測することができ、非常に役に立つ。最近では、モンテカルロ法でドミナンス比を計算する手法も開発されている [39]。
核分裂源分布の収束判定についても、近年研究が進んでおり、情報理論に基づく相対エントロピー (relative

entropy) やシャノン・エントロピー (Shannon entropy) を用いた収束判定法が提案されている [38]。また、2e MCNP コードなどでは、inactive cycle とも呼ばれる。

回の独立なモンテカルロ計算から実効増倍率の収束値を挟み撃ちにするサンドウィッチ法なども提案されている [40]。
核分裂源分布の収束を加速する方法もこれまでにいくつか提案されている。核分裂行列 (fission matrix) を用いる加速法は、古くからよく知られており [41]、最近でも適用した例 [42] がある。また、最近では、Wielandt法をモンテカルロ計算に適用し、核分裂源分布の収束を加速する方法も提案されている [43]。

7.5

アンダーサンプリング

粒子輸送モンテカルロ法で解くのが解くのが困難な問題としてよく知られているのが、Whitesides によって提案された “k-eﬀective of the world” 問題 [44] である。この問題は、臨界安全の分野において、複数の核燃料物質が規則的に離れて配置されている体系の未臨界度を計算するような場合に起こりうる。半径約 4 cm の金属プルトニウム球を 60 cm ピッチで 9 × 9 × 9 で規則上に並べ、その周りを厚い水反射体で囲った体系があり、この体系の実効増倍率は約 0.93 である。この体系の中心にあるプルトニウム球を、裸の単体でちょうど臨界であるプルトニウム球で置き換えると、全体系の実効増倍率は 1 以上になるはずであるが、通常のモンテカルロ
計算では 0.93 の値となるというが、“k-eﬀective of the world” 問題である。
この問題は、中心の臨界であるプルトニウム球に中性子が集まらず、周りの未臨界のプルトニウム球ばかり
で中性子が核分裂反応を起こしている (衝突している) のが原因である。このような現象は、アンダーサンプリ
ング (undersampling) もしくはロスト・ユニット (lost unit) と呼ばれる。プルトニウム球の結合が弱い場合、一度次の世代のスターターがなくなってしまった球は、スターターを得るのが難しくなってしまい、その間、球からの寄与は統計処理にカウントされないことになる。
実際に、“k-eﬀective of the world” 問題を再現するために、MVP コードで図 14に示すような体系をモデル化した。金属プルトニウム球の密度は 19.84 g/cm3 とし、組成は重量割合で Pu-239 が 95.5%、Pu-240 が 4.5%とした。中心以外のプルトニウム球の半径は 4 cm で、中心のプルトニウム球のみ半径 5 cm とした。(組成はどちらの球も同じである。) プルトニウム・アレイの外側は、厚さ 15 cm の水で囲った。すべてのプルトニウム
球を半径 4 cm とし、JENDL-4.0 で計算すると実効増倍率は、0.933781(1σ= 0.0261%) となる。(バッチサイ
ズ 50,000 ヒストリー、スキップバッチ数 200、有効バッチ数 200 で計算した。)(a) 第 5 層 (中央層) の断面図 (中心の燃料球のみ半径 5 cm(b) 第 5 層以外の断面図 (燃料球の半径はすべて 4 cm で、で、超臨界)

未臨界)

図 14: Whitesides の k-eﬀective of the world の問題

中心のプルトニウム球が半径 5 cm の場合の体系に対して、バッチサイズを 100 から 50,000 まで徐々に増やした結果を、図 15 から図 21 に示す。初期の中性子源分布は、水反射体の内側で一様として与えた。各図とも、(a) は、飛程長エスティメータで評価した各バッチにおける増倍率のバッチ数に対する推移を表わしており、(b) は、各バッチ以降の積算平均を表わしている。(b) については、バッチ数が大きくなるほど、平均を評

1.3

1.01

1.2

1.00
k-effective (cumulative)

1.1
k-effective

1.0

0.9

0.8

0.99
0.98
0.97

0.7
0.96

0.6
0.50

200

400
600
Number of bacthes

800

1000

(a) 各バッチにおいて飛程長エスティメータで評価した増

0.950

200

400
600
Number of bacthes

800

1000

(b) 各バッチ以降の増倍率の平均値 (飛程長エスティメータ)倍率

図 15: 飛程長エスティメータで評価した増倍率とその平均値のバッチ数に対する依存性 (バッチサイズ 100 の場合)価するサンプル数が少なくなるので統計誤差が大きくなる。
図 15は、バッチサイズ 100 に対する結果であるが、(a) の図を見ると変動幅が非常に大きく、収束しているかどうかを判定するのは困難である。積算平均値の (b) の図からも、1000 バッチまで計算しても実効増倍率は1 以上となることはなく、バッチサイズ 100 では正しい計算はできないことが分かる。
図 16は、バッチサイズを 500 に増やした結果であるが、依然として各バッチ毎に評価した増倍率の変動幅は大きく、収束を判定することは困難である。(b) の図を見ると 800 バッチ付近でようやく 1 になっているが、その後は大きく 1 を過小評価しており、バッチサイズ 500 でも正確に実効増倍率を評価することは難しい。
図 17は、バッチサイズを 1,000 に増やした結果であり、(a) の図からは大まかに 150 バッチぐらいで一定値に落ち着く傾向が見られるが、まだ変動幅が大きく、収束判定をするのは困難である。(b) の図からは、小さなバッチ数の段階から積算平均値は 1 を超えており、合理的な値が得られている。しかし、150 バッチ目における積算平均値の統計誤差 (1 標準偏差) は 0.1%程度であり、現在の計算精度から考えるともう少し精度のよい結果が望まれる。
図 18は、バッチサイズを 5,000 に増やした結果であり、(a) の図の変動幅もだいぶ小さく抑えられ、100 バッ
チぐらいで一定値に落ち着く傾向がある。(b) の図からも 100 バッチ目ぐらいで一定値となっており、統計処
理から除外する世代数は 100 バッチとすれば十分であることが分かる。
図 19 から図 21 の結果も、バッチサイズ 5,000 の結果の傾向と同様であり、更に変動幅が抑えられ、増倍率の収束の様子が明確に分かる。以上の結果より、この体系については、バッチサイズは最低 5,000 以上 (できれば 10,000 以上) にし、統計処理から除外する世代数は 100 バッチとすることが推奨される。

Whitesides の “k-eﬀective of the world” 問題は、7.4節で述べた収束判定の問題と混同されることがあるが、Blomquist によって指摘されたようにこれはアンダーサンプリングの問題である [45]。また、植木により、この体系のドミナンス比が約 0.7 になることから、この問題はアンダーサンプリングの問題であることが定量的に示されており、情報理論を用いてアンダーサンプリングを検知する手法が提案されている [46]。

7.6核データ表現における不確かさ

原子炉の炉心計算で必要となる断面積は、主に 20 MeV 以下の断面積であり、このエネルギー領域の断面積については評価済み核データを利用する。一般によく用いられている評価済み核データとしては、日本で整備されている JENDL、米国で整備されている ENDF/B、ヨーロッパで整備されている JEFF がある。これらは、それぞれ独立した機関で評価されたものであるので、評価も異なる。20 MeV 以下の中性子輸送計算は評価
済み核データが基礎となっているので、どんなに幾何形状を正確にモデル化し、莫大なヒストリー数で統計精度のよい計算を行っても用いる評価済み核データが異なると計算結果も異なり、評価済み核データ自体に起因する不確かさが生じる。逆にいえば、連続エネルギーモンテカルロ法で体系を正確にモデル化し、精度のよい

1.3

1.01

1.2

1.00
k-effective (cumulative)

1.1
k-effective

1.0

0.9

0.8

0.99
0.98
0.97

0.7
0.96

0.6
0.50

200

400
600
Number of bacthes

800

0.950

1000

(a) 各バッチにおいて飛程長エスティメータで評価した増

200

400
600
Number of bacthes

800

1000

(b) 各バッチ以降の増倍率の平均値 (飛程長エスティメータ)倍率

図 16: 飛程長エスティメータで評価した増倍率とその平均値のバッチ数に対する依存性 (バッチサイズ 500 の
場合)
1.3

1.01

1.2

1.00
k-effective (cumulative)

1.1
k-effective

1.0

0.9

0.8

0.99
0.98
0.97

0.7
0.96

0.6
0.50

200

400
600
Number of bacthes

800

0.950

1000

(a) 各バッチにおいて飛程長エスティメータで評価した増

200

400
600
Number of bacthes

800

1000

(b) 各バッチ以降の増倍率の平均値 (飛程長エスティメータ)倍率

図 17: 飛程長エスティメータで評価した増倍率とその平均値のバッチ数に対する依存性 (バッチサイズ 1,000
の場合)
1.3

1.01

1.2

1.00
k-effective (cumulative)

1.1
k-effective

1.0

0.9

0.8

0.99
0.98
0.97

0.7
0.96

0.6
0.50

200

400
600
Number of bacthes

800

1000

(a) 各バッチにおいて飛程長エスティメータで評価した増

0.950

200

400
600
Number of bacthes

800

1000

(b) 各バッチ以降の増倍率の平均値 (飛程長エスティメータ)倍率

図 18: 飛程長エスティメータで評価した増倍率とその平均値のバッチ数に対する依存性 (バッチサイズ 5,000
の場合)

1.3

1.01

1.2

1.00
k-effective (cumulative)

1.1
k-effective

1.0

0.9

0.8

0.99
0.98
0.97

0.7
0.96

0.6
0.50

200

400
600
Number of bacthes

800

0.950

1000

(a) 各バッチにおいて飛程長エスティメータで評価した増

200

400
600
Number of bacthes

800

1000

(b) 各バッチ以降の増倍率の平均値 (飛程長エスティメータ)倍率

図 19: 飛程長エスティメータで評価した増倍率とその平均値のバッチ数に対する依存性 (バッチサイズ 10,000
の場合)
1.3

1.01

1.2

1.00
k-effective (cumulative)

1.1
k-effective

1.0

0.9

0.8

0.99
0.98
0.97

0.7
0.96

0.6
0.50

200

400
600
Number of bacthes

800

0.950

1000

(a) 各バッチにおいて飛程長エスティメータで評価した増

200

400
600
Number of bacthes

800

1000

(b) 各バッチ以降の増倍率の平均値 (飛程長エスティメータ)倍率

図 20: 飛程長エスティメータで評価した増倍率とその平均値のバッチ数に対する依存性 (バッチサイズ 20,000
の場合)
1.3

1.01

1.2

1.00
k-effective (cumulative)

1.1
k-effective

1.0

0.9

0.8

0.99
0.98
0.97

0.7
0.96

0.6
0.50

200

400
600
Number of bacthes

800

1000

(a) 各バッチにおいて飛程長エスティメータで評価した増

0.950

200

400
600
Number of bacthes

800

1000

(b) 各バッチ以降の増倍率の平均値 (飛程長エスティメータ)倍率

図 21: 飛程長エスティメータで評価した増倍率とその平均値のバッチ数に対する依存性 (バッチサイズ 50,000の場合)

計算をすれば、評価済み核データの不確かさのみを抽出することができ、炉物理実験の実験値と比較することにより評価済み核データの精度向上を図ることができる。このような妥当性評価 (validation) は、積分ベンチ
マークと呼ばれており、これまで数多くの研究がなされ、論文が発表されているので、ここでは取り上げない。
評価済み核データ自体の不確かさとは別に、評価済み核データをモンテカルロ計算用のデータに変換したときに不確かさが生じる。多群モンテカルロ計算では、決定論的手法の輸送計算と同じ多群断面積データを読み込み、その断面積データを元に粒子輸送モンテカルロ計算を行うので、多群断面積に起因する不確かさ (群数の効果、共鳴の取り扱い方法など) は、そのままモンテカルロ計算にも当てはまる。一方、連続エネルギーモ
ンテカルロ計算では、モンテカルロコード毎に中で用いられている物理モデルやサンプリング方法が異なるため、それぞれのモンテカルロコードに特有のフォーマットの断面積データが必要となる。例えば、MCNP コー
ドは、核データ処理システム NJOY[47] を用いて作成される ACE(A Compact ENDF) と呼ばれる形式の断面積データを利用し、MVP コードは、NJOY と同様の核データ処理システム LICEM[48] を用いて作成される

MVP 専用の断面積データを利用する。これらのコードで用いられる断面積データは、ポイントワイズ断面積(pointwise ENDF) が基本となっており、これに熱中性子散乱データ、非分離共鳴データ、その他粒子輸送計算
で必要なデータを加えることによって構成されている。評価済み核データで定義されている様々な反応断面積
は、一般的に反応毎に評価されているエネルギー点やその間の内挿方法が異なっているが、ポイントワイズ断面積を作成する際、線形化 (linearization) することにより、共通のエネルギー点を持つように変換される。この処理では、通常、0.1%の精度で元の断面積と一致するように線形化される。また、モンテカルロコードに依存するが、熱中性子散乱データや非分離共鳴データは、テーブル化されることが多く、有限個数の離散化されたビンで断面積や散乱方向などが表現される。
以上で述べたように、評価済み核データに基づく粒子輸送モンテカルロ計算では、評価済み核データをそのまま用いるのではなく、処理された断面積データを用いている。その処理の際、線形化やテーブル化による不確かさが生じているが、そのような断面積データ処理による不確かさが計算結果に与える影響を調べた例はな
い (公表されていない) と思われる。また、検証のため、処理された断面積データと元の評価済み核データを可視化して比較することがよく行われるが、見た目には非常によく一致しており、断面積データ処理による不確かさが計算結果に与える影響はほとんどないものと考えられる。

7.7

物理モデルにおける不確かさ

粒子輸送モンテカルロコードでは、様々な物理モデルを用いて粒子の輸送をシミュレーションしているが、実際の物理モデルを忠実にモデル化しているわけではなく、影響は小さいと仮定し、近似を導入している場合がある。このような近似によるバイアスは、以前は統計誤差と比べて小さかったが、計算機の発達とともに計算精度が向上し、徐々に見えてくることがある。以下では、比較的最近になって改良された物理モデルの例を取り上げる。

7.7.1

遅発中性子スペクトル

中性子輸送モンテカルロ計算で、遅発中性子を陽にシミュレーションする必要があるのは、主に原子炉動特性計算をするときで、そのようなシミュレーションはまだ研究段階にある。中性子輸送モンテカルロ計算が用いられるのは、主に定常状態の計算においてであり、通常、遅発中性子は即発中性子として取り扱われてきた。
即ち、遅発中性子スペクトルが即発中性子スペクトルとして取り扱われてきた。この近似の影響を MVP コー
ドで調べたのが、表 5である [49]。評価済み核データは JENDL-3.3 を用い、遅発中性子スペクトルを考慮した計算としない計算ともそれぞれ 10 回ずつ乱数列を変えて独立な計算を実行し、その平均から差を算出した。
ほとんどのケースにおいて、遅発中性子スペクトルが実効増倍率に与える効果は、0.1%∆k 以下である。ただし、BIGTEN のような U-238 の核分裂の寄与の大きな体系では 0.2%∆k ほど固有値は小さくなる。
最近のモンテカルロコードでは、ほとんど遅発中性子スペクトルは考慮されているが、オプションを立てないと計算されない場合があるので、注意する必要がある。(正式公開版の MVP 2.0 の場合、DELAYED-NEUTRON
オプションを指定する必要がある。)

表 5: 遅発中性子スペクトルの効果

7.7.2

ケース

燃料

漏れ

∆k （DNS 考慮-なし）(1σ)

GodivaORNL-10BIGTEN

高濃縮 U

57%
7%
11%

0.00012 (0.00003)
0.00039 (0.00003)-0.00183 (0.00005)
0.00079 (0.00008)
0.00069 (0.00008)-0.00023 (0.00003)
0.00008 (0.00004)

高濃縮 U
中濃縮 (10% U-235)

TCA(1.83U)
TCA(3.00U)

低濃縮 U(19 × 19 pins)
低濃縮 U(19 × 19 pins)

4%
6%

JezebelPu(NO3 )4

Pu
Pu

67%
7%

共鳴弾性散乱モデル

ほとんどの中性子輸送モンテカルロコードにおいて、共鳴領域における弾性散乱は、漸近減速散乱モデル(asymptotic slowing-down scattering model, ここでは簡単に漸近モデルと呼ぶことにする) によって取り扱われる。このモデルは、3.5.1節で既に述べたモデルであり、標的核は静止しているものと仮定し、標的核と入射
中性子の運動学 (kinematics) で、弾性散乱を記述するモデルである。熱外領域 (epithermal region) と熱領域においては、3.6節で述べた自由ガスモデルが用いられており、標的核の運動を考慮しているものの、0 K における散乱断面積はエネルギーに依存しないという仮定を用いている。通常、自由ガスモデルと呼ばれるモデル
は、断面積一定を仮定した自由単一原子モデル (free-monatomic-gas model with constant cross section, ここでは簡単に一定断面積モデルと呼ぶことにする) である。最近の研究では、共鳴領域における重い核種に対する弾性散乱モデルとして、漸近モデルや一定断面積モデルを用いると、計算結果に影響を与えることが指摘されており [50]、特に、軽水炉におけるドップラーフィードバック反応度が大きな過小評価となることが指摘されている [51, 52, 53, 54]。
ここでは、MVP コードに実装した共鳴断面積を用いた自由単一原子モデル (free-monatomic-gas model withresonance cross section, ここでは簡単に厳密モデルと呼ぶことにする) を用いて、Mosteller によって提案されたピンセル体系ドップラー反応度ベンチマーク [55] を計算した結果 [56] を図 22に示す。漸近モデルの結果は、厳密モデルの結果と比較すると、絶対値でドップラー反応度を 7%-11%過小評価することが分かる。また、一定断面積モデルも漸近モデルとほぼ同じ結果を与えており、単に標的核の熱運動を考慮するだけではだめなことが分かる。

-1.5

Dopper coefficient (pcm/K)

-2.0
-2.5
-3.0
-3.5
-4.0
-4.5
-5.0
-5.50

1

2

4
3
Enrichment (wt %)Asymptotic
Constant XS
Exact
5

6

図 22: UO2 セルに対するドップラー反応度の計算結果

ここでは示していないが、実効増倍率についても差異があり、漸近モデルの結果は、厳密モデルの結果と比較すると、0.1%-0.2%程度過大評価する。

7.8

幾何形状モデリングから生じる不確かさ

ここでは、モンテカルロ計算において幾何形状モデリングから生じる不確かさについていくつか考えてみる。
製作公差

3.2節で述べたように、粒子輸送モンテカルロ計算における幾何形状のモデリングは、2 次曲面で表現できる体系であれば、正確にモデル化することができ、体系のモデリングによる不確かさは生じない。しかしながら、実験値と計算値を比較する際には、体系の製作公差 (manufacturing tolerance) を考慮しなければならない。国際臨界安全ベンチマークハンドブック (ICSBEP)[57] のように、最近の臨界実験ベンチマーク問題
では、製作公差を正確に評価し、問題に与えることが求められており、製作公差の不確かさが実験値に与える影響が評価されるようになりつつある。このような影響を定量的に評価する手法として、サンプリングによる感度解析 (sampling-based sensitivity and uncertainty analysis)[58] をモンテカルロ計算と組み合わせて行われることがある。
被覆粒子燃料

高温ガス炉では、燃料として被覆粒子燃料が用いられ、JAEA にある高温工学試験研究炉 (HTTR)のようなプリズム型ブロック炉心の場合、大量の被覆粒子燃料がランダムに燃料コンパクトに閉じ込められ、燃料コンパクトはピン状に並べられて、六角形の燃料ブロックに配置される。ぺブル・ベッド炉の場合、大量の被覆粒子燃料はランダムに燃料球に閉じ込められ、燃料球は炉心内にランダムに装荷される。炉心解析においては、このような二重非均質性を考慮することは重要であり、連続エネルギーモンテカルロ計算でこの効果を計算する場合、規則的な格子形状表現を用いて近似をしたり、ランダム性を考慮できる確率論的幾何形状モデル[59] を用いたりして、この効果を評価する。規則的な格子形状表現を用いた近似でも、非常によく二重非均質性を考慮することはできるが、実際の体系とはやはり異なるので、不確かさを生じている可能性がある。また、
確率論的幾何形状モデルも非常によく二重非均質性を考慮することはできるが、このモデルにおいても、最近
接粒子分布にどのような分布を与えるかという問題があり、これにより不確かさを生じている可能性がある。
ボイド分布

BWR 体系の炉心計算や集合体計算においては、軸方向のボイド率分布を考慮する必要がある。

実際には、ボイド率分布は連続的に変化しているが、モンテカルロコードでこのような体系を計算する場合には、軸方向をいくつかのノードに分割し、各ノード内でボイド率が一定であるとして計算する。このとき、実際のボイド率分布と計算で用いるボイド率分布は異なり、結果に差異を生じる可能性がある。このような差異
は、モデリングによる不確かさとなる。
ボイド率分布に適用した例ではないが、物質の断面積が連続的に変化するような体系についてモンテカルロ
法で粒子の飛行パスを正確にサンプリングする手法が提案されており [60]、このような手法を用いれば、物質の密度が連続的に変化するも同様に取り扱える可能性がある。
燃焼計算

最近では、モンテカルロ法を用いた燃焼計算もよく行われるようになっているが、この場合もボイド率分布と同様に、ノード分割による不確かさを生じる。実際の核種数密度分布は連続的に変化するが、モン
テカルロ燃焼計算では、燃焼領域を設定し、その領域では数密度は一定であると仮定して計算を行う。
モンテカルロ燃焼計算における不確かさに関するもう一つの問題は、燃焼に伴う統計誤差の伝搬である。モンテカルロ燃焼計算の統計誤差としては、燃焼方程式を解く時の反応率及び原子数密度の誤差の伝播と輸送方
程式を解く時の原子数密度誤差の伝播を評価する必要がある。これまで既に、この問題についての論文がいく
つか発表されている [61, 62, 63, 64]。

8

おわりに

2節から 6節において粒子輸送モンテカルロ計算の基礎原理と中性子輸送モンテカルロコードを使用するにあたり理解しておいたほうがよいと思われる事項についてまとめた。この部分については、本稿は次の参考文献

中川正幸, “モンテカルロ法による輸送方程式の解法,” 第 22 回炉物理夏期セミナーテキスト 日本
原子力学会 (1990)
長家康展, “モンテカルロ計算の基礎理論及び実験解析への適用,” 第 38 回炉物理夏期セミナーテキ
スト

日本原子力学会 (2006)

をベースに加筆修正したものである。

7節では、モンテカルロ計算で生じうる不確かさ (偶然誤差、バイアス) について述べた。これまでモンテカル
ロ法の固有値計算で課題となっていた、統計誤差のバイアス、収束判定、アンダーサンプリングの問題については、近年、精力的に研究が行われ、有望な評価方法や診断方法が提案されてきている。統計誤差のバイアスについては、よく実際の計算値のばらつきより統計誤差が小さいと言われることが多かったが、比較的正確な
統計誤差が評価できるようになってきている。収束判定については、提案されている診断方法を用いれば、モンテカルロコードのユーザが入力データの設定ミスによるバイアスをなくすことができるようになりつつある。

7節では、また、計算機の発達とともにモンテカルロ計算の精度が向上してきており、これまで統計誤差に隠れて問題にならなかった物理モデルの近似が問題になることがある例を示したが、モンテカルロコードで用いられる物理モデルも改良されてきており、改良された物理モデルを用いれば、物理モデルから生じるバイアスはほぼないものと考えてよい。即ち、適切な入力値で、詳細な物理モデルを用いれば、評価済み核データに基づく粒子輸送計算では、評価済み核データに起因する誤差、統計誤差、幾何形状モデリングから生じる誤差のみであると考えてよいだろう。
最後に、もっと進んで粒子輸送モンテカルロ法を学習するときに役立つであろうと思われる文献をいくつか挙げて本稿を終わりにする。

関連図書
入門書

• E. E. Lewis and W. F. Miller, Jr. : “Computational Method of Neutron Transport,” Chapter 7, John
Wiley & Sons (1984).
• L. L. Carter and E. D. Cashwell : ”Particle-Transport Simulation with the Monte Carlo Method”,
TID-26607 (1975).
実用書• E. D. Cashwell and C. J. Everett : ”A Practical Manual on the Monte Carlo Method for Random Walk
Problems”, Pergamon Press, Inc., New York (1959).
• C. J. Everett and E. D. Cashwell : ”A Third Monte Carlo Sampler (A Revision and Extension of
Samplers I and II)”, LA-9721-MS (1983).
• S. A. Dupree and S. K. Fraley : “A Monte Carlo Primer, A Practical Approach to Radiation Transport,”
Kluwer Academic/Plenum Publishers (2002).
• 「モンテカルロ計算法高度化」研究専門委員会編, “モンテカルロ計算ハンドブック” 日本原子力学会(2006).
理論書• Herman Kahn :”Application of Monte Carlo,” RAND Corporation (1954).
• James Wood : ”Computational Methods in Reactor Shielding”, Chapter 7.1 The Monte Carlo Method,
Pergamon Press (1982).
• J. M. Hammersley and D. C. Handscomb : ”Monte Carlo Methods”, John Wiley & Sons, Inc., NewYork (1964).

• J. Spanier and E. M. Gelbard : ”Monte Carlo Principles and Neutron Transport Problems”, AddisonWesley Publishing Company, Inc., Reading, Mass. (1969).
• I. Lux and L. Koblinger : ”Monte Carlo Particle Transport Methods: Neutron and Photon Calculations”, CRC press (1991).
• Fishman : “Monte Carlo, Concepts, Algorithms, and Applications,” Springer (1996).
乱数についての本• D. E. Knuth : “ The Art of Computer Programming,” Volume2, Third edition, Addison-Wesley (1998).
• J. E. Gentle : “Random Number Generation and Monte Carlo Methods,” Springer (1998).

参考文献[1] G. Buﬀon, “Essai d’arithmétique morale.” Histoire naturelle, générale er particulière, Supplèment 4,
pp.100-104 (1777).
[2] Lord J. W. S. Rayleigh, “ On James Bernoulli’s theorem in probabilities,” Philosophical Magazine, 47,
pp.246-251 (1899).
[3] Lord Kelvin, “Nineteenth Century Clouds over the Dynamical Theory of Heat and Light,” Philosophical
Magazine, Vol. 2, No. 6, pp.1-40 (1901).
[4] R. Courant, K. Friedrichs and H. Lewy, “Über die Partiellen Diﬀerenzengleichungen der Matematischen
Physik, Math. Annalen, 100, pp.32-74 (1928); translated into English by P. Fox, NYO-7689, ABC
Computing Facility, Institute of Mathematical Sciences, New York University (1956).
[5] A. Kolmogorov, “Über die Analytischen Methoden in der Wahrscheinlichkeitsrechnung, Math. Annalen,
104, pp.415-458 (1931).
[6] I. Petrowsky, “Über das Irrfahrtproblem,” Math. Annalen, 109, pp.425-444 (1933).
[7] S. Ulam, R. D. Richtmyer and J. von Neumann, “Statistical methods in neutron diﬀusion,” LAMS-551(1947).
[8] N. Metropolis and S. Ulam, “The Monte Carlo Method,” Journal of the American Statistical Association,
44, pp.335-341 (1949).
[9] N. Metropolis, “The beginning of the Monte Carlo Method,” Los Alamos Science Special Issue (1987).
[10] X-5 Monte Carlo Team, “MCNP – A General Monte Carlo N-Particle Transport Code, Version 5,”
LA-UR-03-1987 (2003).
[11] Y. Nagaya, et al., “MVP/GMVP II : General Purpose Monte Carlo Codes for Neutron and Photon
Transport Calculations based on Continuous Energy and Multigroup Methods,” JAERI 1348 (2005).
[12] D. H. Lehmer, “Mathematical methods in large-scale computing units,” Proceedings of the Second
Symposium on Large Scale Digital Computing Machinery, Harvard University Press, Cambridge, Massachusetts, pp.141-146 (1951).
[13] G. Marsaglia, ”Random Numbers Fall Mainly in the Planes,” Proc. National Academy of Sciences,
U.S.A., 61, pp.25-28 (1968).
[14] C. J. Everett and E. D. Cashwell, “A Third Monte Carlo Sampler,” LA-9721-MS (1983).
[15] L. L. Carter and E. D. Cashwell, “Particle Transport Simulation with the Monte Carlo Method,” TID26607 (1975).

[16] F. B. Brown, W. R. Martin and D. A. Calahan, Trans. Am. Nucl. Soc., 39, 755 (1981).
[17] K. Yamashita, et al., “Nuclear design of the High-Temperature Engineering Test Reactor (HTTR),”
Nucl. Sci. Eng., 122, pp.212-228 (1996).
[18] M. B. Emmett, “The MORSE Monte Carlo Radiation Transport Code System,” ORNL-4972 (1975).
[19] L. L. Carter and C. A. Forest, “Transfer Matrix Treatments for Multigroup Monte Carlo Calculations— The Elimination of Ray Eﬀects,” Nucl. Sci. Eng., 59, pp.27-45 (1976).
[20] L. Mao and J. C. Nimal, “Semianalytical Procedures for Obtaining Low-Order Non-Equally-Probable
Step-Function Representations in Multigroup Monte Carlo,” Nucl. Sci. Eng., 136, pp.409-414 (2000).
[21] 吉沢 康和, “新しい誤差論 実験データ解析法,” 共立出版 (1989).
[22] K. Kobayashi, N. Sugimura, “3D radiation transport benchmark problems and results for simple geometries with void region,” Prog. Nucl. Energy, 39, pp.119-144 (2001).
[23] E. M. Gelbard, R. E Preal, “Monte Carlo Work at Argonne National Laboratory,” in Proc. NEACRP
Meeting on a Monte Carlo Study Group, ANL-75-2 (1974).
[24] V. G. Zolotukhin, L. V. Maiorov, “An estimate of the systematic errors in the calculation of criticality
by the Monte Carlo method,” Atom. Energy, 55, 173-175 (1983).
[25] T. Ueki, T. Mori, M. Nakagawa, “Error Estimations and Their Biases in Monte Carlo Eigenvalue
Calculations,” Nucl. Sci. Eng., 125, 1-11 (1997).
[26] D. B. MacMillan, “Monte Carlo Confidence Limits for Iterated Source Calculations,” Nucl. Sci. Eng.,50, 73 (1973).
[27] E. M. Gelbard, R. E Preal, “Computation of Standard Deviations in Eigenvalue Calculations,” Prog.
Nucl. Energy, 24, 237-241 (1990).
[28] T. J. Urbatsch, R. A. Forster, R. E. Prael, R. J. Beckman, “Estimation and Interpretation of keﬀ
Confidence Intervals in MCNP,” Nucl. Technol., 111, 169-182 (1995).
[29] R. J. Brissenden, A. R. Garlick, “Biases in the estimation of Keﬀ and its error by Monte Carlo methods,”
Ann. Nucl. Energy, 13, 63-83 (1986).
[30] L. Demaret, A. Nouli, L. Carraro, O. Jaquet, “Accurate determination of confidence intervals in Monte
Carlo eigenvalue calculations,” ICNC’99, Sixth Int. Conf. on Nuclear Criticality Safety, Versailles,
France, Sep. 20-24, 1999, Vol. 1, 66-80 (1999).
[31] O. Jaquet, R. Chajari, X. Bay, A. Nouri, L. Carraro, “Eigenvalue uncertainty evaluation in MC calculations, using time series methodologies,” Proc. Monte Carlo 2000 Conf., Lisbon, Portugal, Oct. 23-26,
2000 (2000).
[32] T. Ueki, B. R. Nease, “Time series analysis of Monte Carlo fission sources: II. Confidence interval
estimation,” Nucl. Sci. Eng., 153, 184-191 (2006).
[33] T. Ueki, “Standard Deviation of Local Tallies in Global Monte Carlo Calculation of Nuclear Reactor
Core,” J. Nucl. Sci. Technol., 47, 739-753 (2010).
[34] H. J. Shim, C. H. Kim, “Real Variance Estimation Using an Intercycle Fission Source Correlation forMonte Carlo Eigenvalue Calculations,” Nucl. Sci. Eng., 162, 98-108 (2009).
[35] T. Ueki, F. B. Brown, D. K. Parsons, D. E. Kornreich, “Autocorrelation and dominance ratio in Monte
Carlo criticality calculations,” Nucl. Sci. Eng., 145, 279-290 (2003).

[36] T. Yamamoto, T. Nakamura, Y. Miyoshi, “Fission Source Convergence of Monte Carlo Criticality
Calculations in Weakly Coupled Fissile Arrays,” J. Nucl. Sci. Technol., 37, 41-52 (2000).
[37] http://www.oecd-nea.org/science/wpncs/convergence/[38] T. Ueki, “Stationary Modeling and Informatics-Based Diagnostics in Monte Carlo Criticality Calculations,” Nucl. Sci. Eng., 149, 38-50 (2005).
[39] T. Ueki, F. Brown, D. K. Parsons, “Time Series Analysis of Monte Carlo Fission Sources – I: DominanceRatio Computation,” Nucl. Sci. Eng., 148, 374-390 (2004).
[40] Y. Naito, J. Yang, “The Sandwich Method for Determining Source Convergence in Monte Carlo Calculations,” J. Nucl. Sci. Technol., 41, 559-568 (2004).
[41] K. W. Morton, “Criticality Calculations by Monte Carlo Methods,” AERE T/R 1903, Harwell (1956).
[42] T. Kitada, T. Takeda, “Eﬀective Convergence of Fission Source Distribution in Monte Carlo Simulation,” J. Nucl. Sci. Technol., 38, 324-329 (2001).
[43] T. Yamamoto, Y. Miyoshi, “Reliable Method for Fission Source Convergence of Monte Carlo Criticality
Calculation with Wielandt’s Method,” J. Nucl. Sci. Technol., 41, 99-107 (2004).
[44] G. E. Whitesides, “A Diﬃculty in Computing the k-eﬀective of the World,” Trans. Am. Nucl. Soc., 14,
680 (1971).
[45] R. N. Blomquist, “Monte Carlo Source Convergence and the Whitesides Problem,” Proc. Int. Topl. Mtg.
Advances in Reactor Physics and Mathematics and Computation into the Next Millennium (PHYSOR
2000), Pittsburgh, Pennsylvania, May 7-12, 2000, American Nuclear Society (2000) [DVD-ROM].
[46] T. Ueki, “Information Theory and Undersampling Diagnostics for Monte Carlo Simulation of Nuclear
Criticality,” Nucl. Sci. Eng., 151, 283-292 (2005).
[47] R. E. MacFarlane, A. C. Kahler, “Methods for Processing ENDF/B-VII with NJOY,” Nucl. Data Sheets,
111, 2739-2890 (2010).
[48] 森 貴正、長家 康展、奥村 啓介、金子 邦男, “最新の評価済み核データに基づく MVP 中性子断面積
ライブラリーの整備,” JAERI-Data/Code 2004-011 (2004).

[49] 長家 康展、奥村 啓介、森 貴正, “MVP を用いた遅発中性子のシミュレーション,” 日本原子力学会
2004 年秋の年会, B33 (2004).
[50] M. Ouisloumen, R. Sanchez, “A Model for Neutron Scattering oﬀ Heavy Isotopes that Accounts for
Thermal Agitation Eﬀects,” Nucl. Sci. Eng., 107, 189-200 (1991).
[51] O. Bouland, V. Kolesov, J. L. Rowlands, “The Eﬀect of Approximations in the Energy Distributions of
Scattered Neutrons on Thermal Reactor Doppler Eﬀects,” JEF/DOC-486 (1994).
[52] V. V. Kolesov, V. F. Ukraintsev, “Temperature Eﬀects and Resonance Elastic Cross Section Influence on
Secondary Energy Distributions of Scattered Neutrons in the Resolved Resonance Region,” PHYSOR2006, Vancouver, Canada (2006).
[53] R. Dagan, C. H. M. Broeders, “On the Eﬀect of Resonance Dependent Scattering-kernel on Fuel Cycle
and Inventory,” PHYSOR-2006, Vancouver, Canada (2006).
[54] D. Lee, K. Smith, J. Rhodes, “The impact of 238 U resonance elastic scattering approximations on
thermal reactor Doppler reactivity,” Int. Conf. on Reactor Physics, Nuclear Power: A Sustainable
Resource, Casino-Kursaal Conference Center, Interlaken, Switzerland, Sep. 14-19, 2008 (2008).

[55] R. D. Mosteller, “Computational Benchmarks for the Doppler Reactivity Defect,” LA-UR-06-2968(2006).
[56] T. Mori, Y. Nagaya, “Comparison of Resonance Elastic Scattering Models Newly Implemented in MVP
Continuous-Energy Monte Carlo Code,” J. Nucl. Sci. Technol., 46, 793-798 (2009).
[57] Nuclear Energy Agency, “International handbook of evaluated criticality safety benchmark experiments,” NEA/NSC/DOC(95)03, September 2008 Edition (2008), [CD-ROM].
[58] D. G. Cacuci, “Sensitivity and Uncertainty Analysis of Models and Data,” Section 6.3.1, in (Eds.) Y.
Azmy, E. Sartori, Nuclear Computational Science, A Century in Review, Springer (2010).
[59] I. Murata, T. Mori, M. Nakagawa, “Continuous Energy Monte Carlo Calculations of Randomly Distributed Spherical Fuels in High-Temperature Gas-Cooled Reactors Based on a Statistical Geometry
Model,” Nucl. Sci. Eng., 123, 96 (1996).
[60] F. B. Brown, W. R. Martin, “Direct sampling of Monte Carlo flight paths in media with continuously
varying cross-sections,” Nuclear Mathematical and Computational Sciences: A Century in Review, A
Century Anew, Gatlinburg, Tennessee, April 6-11, 2003 (2003), [CD-ROM].
[61] T. Takeda, N. Hirokawa and T. Noda, “Estimation of Error Propagation in Monte-Carlo Burn-up
Calculations,” J. Nucl. Sci. Technol. 36, 738 (1999).
[62] M. Tohjoh, T. Endo, M. Watanabe, A. Yamamoto, “Eﬀect of error propagation of nuclide number
densities on Monte Carlo burn-up calculations,” Ann. Nucl. Energy, 33 1424-1436 (2006).
[63] H. J. Shim, H. J. Park, C. H. Kim, “Error propagation in Monte Carlo depletion analysis,” Trans. Am.
Nucl. Soc., 96, 582-584 (2007).
[64] E. Dumonteil, C. M. Diop, “Biases and Statistical Errors in Monte Carlo Burnup Calculations: An
Unbiased Stochastic Scheme to Solve Boltzmann/Bateman Coupled Equations,” Nucl. Sci. Eng., 167,
165-170 (2011).