空間均質化誤差低減手法と高速化手法
日本原子力研究開発機構
多田 健一
1. はじめに本資料では核計算コードで利用されている高速化手法について、その概要と特徴を理解してもらうことを目的としている。タイトルにある通り、本資料では最初に空間均質化誤差低減手法について説明した上で、拡散計算コードや輸送計算コード 1で利用されている高速化手法について説明する。
空間均質化誤差低減手法と高速化手法には一見関連がないように見えるが、Generalized
Coarse Mesh Rebalance (GCMR)法 4)など、現在集合体計算コードで広く使われている粗メッシュ
加速法と呼ばれる拡散計算を使った高速化手法には、空間均質化誤差低減手法が用いられている。GCMR 法を理解することは核計算コードの開発や入力値の設定の際に大きく役立つことから、本資料ではまず第 2 章で空間均質化誤差低減手法について説明し、続いて第 3 章で内部反復、外部反復それぞれの高速化手法について説明する。
なお、各手法には出来る限り参考文献をつけたので、各手法について詳しく知りたい場合には、そちらを参考にして欲しい。また、配列アクセスの最適化(ストライドアクセスから連続アクセスへの変更)やプロファイラを利用したボトルネックの調査とその改善といった一般的な
プログラミングの高速化テクニックについては、様々な書籍が刊行されており、またインター
ネット上でも簡単に調べることができることから、ここでは割愛する。

2. 空間均質化と空間均質化誤差低減手法について
2.1 空間均質化とは
軽水炉の炉心設計では、図 1 のように燃料棒、被覆管、減速材領域から構成される単ピンセ
ル体系において、共鳴干渉効果を考慮した実効断面積を計算するセル計算、集合体平均断面積を計算する二次元集合体計算、そして原子炉の実効増倍率や燃料棒単位の出力分布などを計算する三次元全炉心計算の三つの計算を組み合わせている 5-7)。これらの計算の内、空間均質化の問題が顕在化するのは全炉心計算である。二次元集合体計算では図 1 のように燃料棒やウォーターロッドなどの集合体内の詳細な構造を考慮した計算を行っているが、三次元全炉心計算
で集合体内の詳細な構造を考慮することは、現在の計算機資源をもってしても現実的ではない
2。そのため、図 2 や図 3 のように集合体内を全て同じ均質な物質と見なすことで高速な計算

1炉物理分野では、Characteristics 法(MOC)1)や衝突確率法 2)、SN 法 3)など拡散計算以外の手法を輸送計算と呼んで区別している。
近年ではスーパーコンピューターを利用すれば集合体内の詳細な構造を考慮した三次元全炉心計算ができるようになってきているが、炉心設計業務で利用することはできない。

2

1

を実現している 3。図 2 や図 3 のように集合体内を均質な物質とみなすことを均質化という。
なお、単に均質化と言った場合、エネルギー群数を縮約することと混同することから、幾何形
状を均質化することを空間均質化と呼んでいる 4。

セル計算(実効断面積の作成、
数百群)

二次元集合体計算(集合体平均断面積の作成、10群)

三次元全炉心計算(2～3群)

図 1 炉心解析に用いる断面積作成の流れ

図 2 BWR 燃料集合体の均質化の例

図 3 PWR 燃料集合体の均質化の例
エネルギー群数の縮約も同様であるが、空間均質化では反応率を保存するように均質化を行う。例えば(1)式に示すような一次元拡散方程式を考えると、方程式を解く上で重要となるのは巨視的断面積 Σ と中性子束𝜙𝜙の積である反応率である。
3

炉心計算コードによっては集合体内を分割しているものもあるが、集合体内の詳細な構造を考慮しないという点では同じである。
4 本章で説明する空間均質化誤差低減手法はエネルギー群数の縮約に伴う誤差の低減にも利用できるが、本資料では説明を簡略化するため、空間均質化誤差に焦点を絞って説明する。
2

−𝐷𝐷𝑔𝑔,𝑖𝑖𝜒𝜒𝑔𝑔,𝑖𝑖d2
𝜙𝜙𝑔𝑔,𝑖𝑖 + Σ𝑡𝑡,𝑔𝑔,𝑖𝑖 𝜙𝜙𝑔𝑔,𝑖𝑖 =� 𝜈𝜈Σ ′ 𝜙𝜙 ′ + � Σ𝑠𝑠,𝑔𝑔′ →𝑔𝑔,𝑖𝑖 𝜙𝜙𝑔𝑔′ ,𝑖𝑖2d𝑥𝑥
𝑘𝑘 ′ 𝑓𝑓,𝑔𝑔 ,𝑖𝑖 𝑔𝑔 ,𝑖𝑖
′
𝑔𝑔

(1)

𝑔𝑔

ここで、g はエネルギー群、i はメッシュ番号、D は拡散係数、Σt は巨視的全断面積、k は実効増倍率、χ は核分裂スペクトル、ν は核分裂当たりに放出される平均中性子数、Σf は核分裂断面積、Σs,g’→g は g’群から g 群に散乱する散乱断面積である 5。
エネルギー群を G 群に縮約し、メッシュを粗メッシュ I に均質化することを考える。各メッ
シュでの反応率を保存させるため、中性子束𝜙𝜙とメッシュの体積 V の積を重みにして次式を用いてエネルギー・空間を均質化している 6。
Σ𝐺𝐺,𝐼𝐼 =∑𝑔𝑔∈𝐺𝐺 ∑𝑖𝑖∈𝐼𝐼 Σ𝑔𝑔,𝑖𝑖 𝜙𝜙𝑔𝑔,𝑖𝑖 𝑉𝑉𝑖𝑖
∑𝑔𝑔∈𝐺𝐺 ∑𝑖𝑖∈𝐼𝐼 𝜙𝜙𝑔𝑔,𝑖𝑖 𝑉𝑉𝑖𝑖(2)(2)式では、(3)式が成立すると仮定している。
Σ𝐺𝐺,𝐼𝐼 𝜙𝜙𝐺𝐺,𝐼𝐼 𝑉𝑉𝐼𝐼 = Σ𝐺𝐺,𝐼𝐼 � � 𝜙𝜙𝑔𝑔,𝑖𝑖 𝑉𝑉𝑖𝑖 = � � Σ𝑔𝑔,𝑖𝑖 𝜙𝜙𝑔𝑔,𝑖𝑖 𝑉𝑉𝑖𝑖
𝑔𝑔∈𝐺𝐺 𝑖𝑖∈𝐼𝐼

すなわち、

𝜙𝜙𝐺𝐺,𝐼𝐼 𝑉𝑉𝐼𝐼 = � � 𝜙𝜙𝑔𝑔,𝑖𝑖 𝑉𝑉𝑖𝑖

(3)

𝑔𝑔∈𝐺𝐺 𝑖𝑖∈𝐼𝐼(4)

𝑔𝑔∈𝐺𝐺 𝑖𝑖∈𝐼𝐼である。(3)、(4)式が成立するのであれば、非均質な体系の反応率と均質な体系の反応率は一致する。しかし、(2)式で求めた均質化断面積を用いた均質体系の計算で(3)、(4)式が成立する保証はない。例えば全てのメッシュが同じ断面積を持った体系であれば当然ながら(3)、(4)式は成立する。一般的な非均質体系でも(3)、(4)式は成立するように見える。しかし、一般的な非均質体系では(4)式は成立しない。そのため、均質化に伴う誤差である空間均質化誤差が発生してしまう。
なぜ(4)式が成立しないのかを考えるため、図 4 のような体系の中心に吸収断面積が非常に大きい中性子強吸収体を置いた体系を考える。図 4 のように中性子が左下から右上の方向にのみ飛行する場合、非均質体系では中性子強吸収体を通過する中性子は吸収されてしまうが、それ以外の領域を通る中性子はそのまま通過する。このような体系を均質化した場合、均質体系の吸収断面積は非常に大きくなる。そのため、メッシュ内の大半の領域で中性子が吸収されてしまい、非均質体系と均質体系で中性子の吸収量の総和が異なってしまう。このように反応率が保存するように均質化を行ったとしても、均質化によって断面積が平均化されることで体系内の吸収量に差異が生じてしまう。この体系内の吸収量の差異が空間均質化誤差の要因となる。
さて、読者の皆さんは上記の説明で空間均質化誤差について納得して頂けただろうか？学生時代の筆者であったら何の疑いもなく納得したことと思う 7。しかし、上記の説明では空間均(1)式はあくまでイメージであり、(1)式中の中性子束はメッシュ中央の中性子束とメッシュ平均の中性子束が混在しているので注意する必要がある。
6 ここでは断面積のみを示しているが、拡散係数 D を均質化する場合も断面積と同じように(2)式の断面積を拡散係数に置き換えればよい。
7 実際、筆者はこの資料を書き始めるまで、この説明で空間均質化誤差について全て説明できていると思っていた。
53

質化誤差を完全に説明したことになっていない。

吸収体を通過する中性子
のみ吸収

非均質体系

セル内の大半の領域で中性子が吸収

均質体系

図 4 中性子強吸収体を含む体系の均質化の例 (1)
先ほどの説明の問題点について解説する。反応率が保存するように均質化を行っている以上、中性子の吸収量の総量は変わらないはずである。仮に図 4 のような中性子強吸収体があったとしても、中性子強吸収体が置かれた領域では中性子束はほとんどゼロになることから、反応率
でみれば他の領域とほとんど変わらないはずである。その場合、中性子の吸収量は非均質体系と均質体系で変わらないことから、(3)、(4)式は成立するように見える。
(3)、(4)式が成立しない原因は正確には吸収量の変化ではなく、体系境界での中性子の入射に対する放出の関係が変化するためである 8。図 4 では説明が分かりにくいので、図 5 のように
中性子強吸収体が中心ではなく左上に置かれた体系を考える。この場合、非均質体系では図 4の体系と同じように中性子強吸収体を通過する中性子が吸収されてしまうことから、体系の上の辺に到達する中性子の数が少なくなる一方、体系の右の辺には全ての中性子が到達する。それに対し均質体系では体系全体で中性子が吸収されてしまうことから、非均質体系に比べて体系の上の辺に到達する中性子が増える一方、体系の右の辺に到達しない中性子が出てくることになる。図 5 のような体系では、仮に体系内での中性子の吸収量が同じであったとしても、非均質体系と均質体系で体系を出る中性子の数、すなわち体系境界での中性子流が異なるのは明白である。
このように、非均質体系と均質体系とで境界の中性子流が異なることが、空間均質化誤差の要因である 9。また、今回の説明で用いたように、空間均質化誤差は Gd 入り燃料棒や制御棒など、中性子強吸収体が存在する体系で特に大きくなる。そのため、中性子強吸収体が存在する体系を均質化する際には、空間均質化誤差の影響を考慮することが特に重要となる。
近年の燃料集合体は高燃焼度化に伴い、集合体内出力分布を平坦化させるため、濃縮度が異なる様々な燃料棒で構成されており、従来の燃料集合体に比べて非均質性が高くなっている。

8

中性子の入射と流出の関係が中性子流であり、中性子流で考えると、非均質体系と均質体系
で中性子流が一致していないことを意味している。
9 なお、非均質体系と均質体系の中性子流が等しくなるように断面積を均質化したとしても、体系内の反応率が保存されず、非均質体系と均質体系で反応率差異が生じてしまうので解決策にはならない。
4

また、Gd 濃度が高い Gd 入り燃料棒が入っていることも多いことから、高精度な炉心解析を行うためには空間均質化誤差を無視することは出来ない。

全ての中性子が通過

一部の中性子がセル内で吸収

非均質体系

均質体系

図 5 中性子強吸収体を含む体系の均質化の例 (2)2.2 空間均質化誤差低減手法の適用
空間均質化誤差を低減する手法として、Superhomogenization (SPH)法 8)や Discontinuity factor(不連続因子、もしくは DF)9-11)が考案され、炉心計算コードで広く利用されている。これらの空間均質化誤差低減手法は反応率だけでなく、中性子流も保存させることで空間均質化誤差の低減を図っている。以降ではこれらの手法について説明する。
2.2.1 SPH 法について
SPH 法では(3)式が成立するように、SPH 因子 μ を断面積に掛ける。
𝜇𝜇𝐺𝐺,𝐼𝐼 Σ𝐺𝐺,𝐼𝐼 𝜙𝜙𝐺𝐺,𝐼𝐼 V𝐼𝐼 = � � Σ𝑔𝑔,𝑖𝑖 𝜙𝜙𝑔𝑔,𝑖𝑖 𝑉𝑉𝑖𝑖
𝑔𝑔

(5)

𝑖𝑖

SPH 法を(1)式の一次元拡散方程式に適用する場合、以下のように書き表せられる。
−𝜇𝜇𝐺𝐺,𝐼𝐼 𝐷𝐷𝐺𝐺,𝐼𝐼d2
𝜙𝜙 + 𝜇𝜇𝐺𝐺,𝐼𝐼 Σ𝑡𝑡,𝐺𝐺,𝐼𝐼 ϕ𝐺𝐺,𝐼𝐼d𝑥𝑥 2 𝐺𝐺,𝐼𝐼𝜒𝜒𝐺𝐺,𝐼𝐼=� 𝜇𝜇 ′ 𝜈𝜈Σ ′ 𝜙𝜙 ′ + � 𝜇𝜇𝐺𝐺 ′ ,𝐼𝐼 Σ𝑠𝑠,𝐺𝐺 ′ →𝐺𝐺,𝐼𝐼 𝜙𝜙𝐺𝐺 ′ ,𝐼𝐼𝑘𝑘 ′ 𝐺𝐺 ,𝐼𝐼 𝑓𝑓,𝐺𝐺 ,𝐼𝐼 𝐺𝐺 ,𝐼𝐼
′
𝐺𝐺

(6)

𝐺𝐺

(6)式では分かりやすくするために𝜇𝜇𝐺𝐺,𝐼𝐼 を陽に書いているため、複雑に見える。しかし、SPH 因

子 μ を掛けたものを断面積として与えれば、計算コードを書き換えることも余分な計算負荷をかけることもなく SPH 法を適用することが出来る。このように SPH 法は計算手法の制約が無く、適用範囲が広い手法である。
SPH 法で必要な SPH 因子の計算について説明する。SPH 因子は𝜇𝜇𝐺𝐺,𝐼𝐼 =∑𝑔𝑔 ∑𝑖𝑖 𝜙𝜙𝑔𝑔,𝑖𝑖 𝑉𝑉𝑖𝑖
𝜙𝜙𝐺𝐺,𝐼𝐼 𝑉𝑉𝐼𝐼(7)で計算される。(7)式のように SPH 因子の計算には均質化後の中性子束が必要となる。均質化
後の中性子束𝜙𝜙𝐺𝐺,𝐼𝐼 は(6)式から得られることから、SPH 因子を求めるためには(6)式を用いた均質化後の中性子束𝜙𝜙𝐺𝐺,𝐼𝐼 の計算と(7)式を用いた SPH 因子の計算を繰り返す反復計算が必要となる。
5

SPH 因子の計算手順は次の通りである。
1. (2)式を用いて非均質体系の断面積Σ𝑔𝑔,𝑖𝑖 、中性子束𝜙𝜙𝑔𝑔,𝑖𝑖 から均質体系の断面積Σ𝐺𝐺,𝐼𝐼 を計算する。

2. (7)式を用いて SPH 因子𝜇𝜇𝐺𝐺,𝐼𝐼 を計算する。なお、𝜇𝜇𝐺𝐺,𝐼𝐼 の初期値は𝜇𝜇𝐺𝐺,𝐼𝐼 = 1.0とする。

3. SPH 因子を乗じて補正した均質体系の断面積Σ′𝐺𝐺,𝐼𝐼 = 𝜇𝜇𝐺𝐺,𝐼𝐼 Σ𝐺𝐺,𝐼𝐼 を計算し、(6)式を用いて均質体系を計算する。

4. 2.と 3.を SPH 因子が(8)式のように収束条件を満たすまで繰り返す。なお、(8)式の n は反復回数で、𝜀𝜀は収束条件であり、一般的には𝜀𝜀 = 1.0 × 10−2程度とすることが多い。
𝑛𝑛
𝑛𝑛−1
𝜇𝜇𝐺𝐺,𝐼𝐼− 𝜇𝜇𝐺𝐺,𝐼𝐼< 𝜀𝜀𝑛𝑛𝜇𝜇𝐺𝐺,𝐼𝐼(8)

SPH 法を適用するには、非均質体系の中性子束が必要となる。図 1 に示すような一般的な炉心計算では、集合体計算の結果を再現するように同じ燃料集合体が無限に並んだ単一集合体体系で SPH 法を適用する。単一集合体体系では同じ燃料が隣接することから、集合体内の中性子束分布は(Gd 燃料棒や制御棒、水ロッドなどの変化を無視すれば)平坦な形となり、また集合体
境界の中性子流は零となる。しかし、実際の炉心計算では濃縮度分布や燃焼度の異なる燃料集合体が隣接しているなど、集合体計算時とは中性子束分布や中性子流分布が異なっている 10。
そのような場合、単一集合体体系において SPH 法を用いて空間均質化誤差を低減したとしても、異なる燃料と隣接する集合体境界で非均質体系との差異が見られる。しかし、その差異は十分に小さく、SPH 法を適用していない((2)式を使って均質化しただけの)場合に比べ、炉心計算において高い計算精度が期待できる。
ただし、SPH 法はあくまで SPH 因子を作成した際の非均質体系の結果を再現するように補正を行うものであり、UO2 燃料集合体と MOX 燃料集合体が隣接する境界など、SPH 因子を作成した際の非均質体系と中性子流が大きく異なる場合には、計算精度が悪化する場合がある 11。
SPH 因子の計算に必要な反復回数は十回～数十回程度である。一般的に均質体系の計算には拡散計算や Simplified P3 法 12)など、詳細輸送計算手法に比べて非常に高速な計算手法が用いられる。そのため、単一集合体体系で SPH 因子を求める計算時間は 1 分もかからない 12。また、集合体計算による核定数の準備には長い計算時間が許容されるため、SPH 法を適用するために
10

例として新しい燃料集合体と燃焼度の高い燃料集合体が隣接した場合を考える。燃焼度の高い燃料集合体は新しい燃料集合体に比べて核分裂性核種の量が少なく、出力(中性子束)が低くなる。そのため、新しい燃料集合体に比べて燃焼度の高い燃料集合体からの中性子の流れ込みが少なくなり、
集合体境界の中性子流は零ではない値となる(新しい燃料集合体から燃焼度の高い燃料集合体へと中性子が流出する)。このように、燃焼度の高い燃料集合体が隣接すると、集合体内の中性子束分布は傾いた分布となる。
11 濃縮度や燃焼度の異なる燃料集合体が隣接している場合に比べ、UO2 燃料集合体と MOX 燃
料集合体が隣接する場合は計算精度がやや悪化するものの、軽水炉の炉心計算に用いる範囲内であれば空間均質化誤差を低減する手法として非常に効果的な手法である。
12 燃料集合体内を燃料棒単位で均質化した pin-by-pin 体系でも SPH 因子の計算時間は 1 分程度である。一般的な炉心計算で取り扱う燃料集合体内を 1 メッシュとした計算では、数秒以内に終了するため、計算時間の増加はほとんどない。
6

は反復計算が必要というデメリットは、SPH 法を炉心解析コードに適用する場合には大きな問題にはならない。
2.2.2 不連続因子について(1)式の拡散方程式を解くことを考える。(1)式の右辺の中性子源の項を既知と見なし N 個の空間メッシュがあるとすると、(1)式は N 本の二階の微分方程式となる。二階の微分方程式を解くのに必要な条件は 2 つであり、一般的な拡散計算コードではメッシュ境界の中性子束(表面中性子束)とメッシュ境界の中性子流(表面中性子流)が連続であるという条件を用いて、この N 本の二階の微分方程式を解いている。
拡散方程式以外の計算手法の場合についても同様である。輸送計算で基本となる微積分型ボルツマン輸送方程式を用いた場合には(9)式の通り、一階の微分方程式となる。
��⃗ ∙ ∇𝜙𝜙𝑔𝑔 �𝑟𝑟⃗, Ω��⃗� + Σ𝑡𝑡,𝑔𝑔 (𝑟𝑟⃗)𝜙𝜙𝑔𝑔 �𝑟𝑟⃗, Ω��⃗� = 𝑄𝑄𝑔𝑔 �𝑟𝑟⃗, Ω
��⃗�
Ω

(9)

��⃗は中性子の角度方向、𝑟𝑟⃗は位置である。この場合、(9)式を解くのに必要な条件は 1
ここで、Ωつであり、メッシュ境界の片側の表面角度中性子束が連続であるという条件(これによって、中性子流も連続となる)を用いることが多い。
ここで、非均質体系の一部の空間領域をエネルギー1 群で 1 メッシュの均質体系に縮約・均質化する場合を考える。この時、非均質体系については、中性子束と中性子流は既知とする。
この縮約・均質化された領域（均質体系と呼ぶ）を対象として、(1)式を用いて計算するためには、2 つの境界条件が必要となる。この境界条件として、均質体系境界の表面中性子束が非均質体系のメッシュ境界の表面中性子束と一致するという条件を用いて計算を行うと、均質体系のメッシュ境界の中性子流は非均質体系のメッシュ境界の中性子流と一致しない。逆に均質体系のメッシュ境界の中性子流が非均質体系のメッシュ境界の中性子流と一致するという条件を用いて計算を行うと、今度は均質体系境界の表面中性子束が非均質体系のメッシュ境界の表
面中性子束と一致しなくなる。
このように境界条件が 2 つである以上、均質体系で非均質体系のメッシュ境界の表面中性子
束と中性子流を同時に一致させる計算を行うことは出来ない。これが空間均質化誤差の要因である。ではどうすれば空間均質化誤差を低減させることが出来るだろうか？SPH 法では、非均質体系と反応率が一致するように SPH 因子を用いて断面積を補正した。不連続因子では、メッシュ境界の均質体系の中性子流を非均質体系の中性子流と一致させる代わりに、均質体系と非均質体系のメッシュ境界の表面中性子束で差異が生じるのを許容している。そして、均質メッシュの境界において、不連続因子𝐷𝐷𝐷𝐷𝑠𝑠±,𝐺𝐺,𝐼𝐼 （今、一次元で考えているので、境界両側 𝑠𝑠 ±）を次式で定義する。

𝐷𝐷𝐷𝐷𝑠𝑠±,𝐺𝐺,𝐼𝐼 =∑𝑖𝑖∈𝐼𝐼 ∑𝑔𝑔∈𝐺𝐺 𝜙𝜙𝑠𝑠±,𝑔𝑔,𝑖𝑖 𝑆𝑆𝑠𝑠±,𝑖𝑖
𝜙𝜙𝑠𝑠±,𝐺𝐺,𝐼𝐼 𝑆𝑆𝑠𝑠±,𝐼𝐼(10)

ここで、𝑆𝑆𝑠𝑠±,𝐼𝐼 は均質メッシュ I の表面積である。(10)式では、非均質体系では均質メッシュ Iの両側の境界に対応する非均質メッシュ領域は異なるが、両者を i としている。すなわち、(10)式の𝜙𝜙𝑠𝑠±,𝑔𝑔,𝑖𝑖 は非均質体系における均質メッシュの両側境界に対応する表面中性子束である。な7

お、表面中性子流を一致させることによって、どのような縮約・均質化断面積 Σ𝑡𝑡,𝐺𝐺,𝐼𝐼 を使用しても、均質体系での反応率 Σ𝑡𝑡,𝐺𝐺,𝐼𝐼 𝜙𝜙𝐺𝐺,𝐼𝐼 𝑉𝑉𝐼𝐼 を非均質体系での反応率に一致させることができる。

次に、隣り合う 2 つの空間メッシュ(均質メッシュ𝐼𝐼 − 1と𝐼𝐼)で縮約・均質化する場合を考える。

非均質体系では 2 つのメッシュ境界で中性子束も連続なので、各メッシュの不連続因子を用い
ると、
𝐷𝐷𝐷𝐷𝑠𝑠+,𝐺𝐺,𝐼𝐼−1 𝜙𝜙𝑠𝑠+,𝐺𝐺,𝐼𝐼−1 =∑𝑖𝑖∈𝐼𝐼 ∑𝑔𝑔∈𝐺𝐺 𝜙𝜙𝑠𝑠+,𝑔𝑔,𝑖𝑖−1 𝑆𝑆𝑠𝑠+,𝑖𝑖−1 ∑𝑖𝑖∈𝐼𝐼 ∑𝑔𝑔∈𝐺𝐺 𝜙𝜙𝑠𝑠−,𝑔𝑔,𝑖𝑖 𝑆𝑆𝑠𝑠−,𝑖𝑖== 𝐷𝐷𝐷𝐷𝑠𝑠−,𝐺𝐺,𝐼𝐼 𝜙𝜙𝑠𝑠−,𝐺𝐺,𝐼𝐼
𝑆𝑆𝑠𝑠+,𝐼𝐼−1
𝑆𝑆𝑠𝑠−,𝐼𝐼(11)という関係が導き出される。メッシュ境界の表面積は同じであることから、
𝑆𝑆𝑠𝑠+,𝐼𝐼−1 = � 𝑆𝑆𝑠𝑠+,𝑖𝑖−1 = � 𝑆𝑆𝑠𝑠−,𝑖𝑖 = 𝑆𝑆𝑠𝑠−,𝐼𝐼
𝑖𝑖∈𝐼𝐼

(12)

𝑖𝑖∈𝐼𝐼

である。(11)、(12)式から、均質体系におけるメッシュ境界での表面中性子束の関係は
𝜙𝜙𝑠𝑠−,𝐺𝐺,𝐼𝐼 =𝐷𝐷𝐷𝐷𝑠𝑠+,𝐺𝐺,𝐼𝐼−1
𝜙𝜙
𝐷𝐷𝐷𝐷𝑠𝑠−,𝐺𝐺,𝐼𝐼 𝑠𝑠+,𝐺𝐺,𝐼𝐼−1(13)となる。ここでは、簡略化のため 𝜙𝜙𝑠𝑠±,𝐺𝐺,𝑖𝑖 = ∑𝑔𝑔∈𝐺𝐺 𝜙𝜙𝑠𝑠±,𝑔𝑔,𝑖𝑖 と表記した。(10)式は、均質体系においては、均質メッシュ境界において表面中性子束が境界両側の不連続因子で決まる量だけ不連続となることを意味している。従って、不連続因子を適用する場合、均質メッシュ境界における境界条件が、使用しない場合の①中性子流の連続＋②中性子束の連続に代えて、①中性子流の連続＋②中性子束の不連続(不連続量は(10)式で評価)となる。
ここで、不連続因子の計算方法について考える。不連続因子の計算では、均質体系のメッシ
ュ境界の中性子流が非均質体系でのメッシュ境界の中性子流と一致するという条件を用いて計算を行うと説明した。しかし、実際の拡散計算や輸送計算では、一般的に(1)式や(9)式のように中性子束を用いて計算を行っており、中性子流を陽に取り扱って計算を行うことは少ない 13。
もちろん、計算コードを書き換えて表面中性子流を境界条件として取り扱えるようにするという方法もあるが、炉心計算で不連続因子を適用する場合にはもっと簡単な方法で均質体系と非均質体系の中性子流を一致させた計算を行っている。
一般的な集合体計算では、同じ集合体が無限に並んでいる単一集合体体系で計算を行っている。単一集合体体系で計算を行うためには、体系を 1×1 集合体とし、その境界を完全反射境界
条件とすればよい。完全反射境界条件を用いた場合、メッシュ境界では燃料配置や水ギャップの幅、制御棒の有無などが全て同じ条件となっていることから、メッシュ境界での中性子の流入や流出はない。このように、完全反射境界条件を用いれば、メッシュ境界の中性子流は必ず零になる 14。そのため、不連続因子の計算に単一集合体体系のように完全反射境界条件を用いれば、メッシュ境界の中性子流を陽に取り扱わなくても非均質体系と均質体系のメッシュ境界
13

本セミナーの別の講義資料にある応答行列法のように、境界の部分中性子流を陽に取り扱う手法もある。
14 完全反射境界条件と同じように無限に並べる方法として、
周期境界条件を用いる場合がある。
完全反射境界条件は境界に鏡を置いたように反転した体系が隣接した体系であり、周期境界条
件では同じ体系が向きを変えずに並んだ体系である。周期境界条件では集合体境界での中性子流が零にならないなど、両者は明確に異なっているので自分が解析したい体系はどちらなのか
よく考えて境界条件を選択する必要がある。
8

の中性子流を一致させることが出来る。このように、炉心計算で使用する不連続因子を求める
場合に、単一集合体体系で計算を行うことで、計算手法や計算式を変えることなく、メッシュ
境界の中性子流を非均質体系と一致させた計算を行うことができる。
単一集合体体系における不連続因子の計算手順は次の通りである。
1. (2)式を用いて非均質体系の断面積Σ𝑔𝑔,𝑖𝑖 、中性子束𝜙𝜙𝑔𝑔,𝑖𝑖 から均質断面積Σ𝐺𝐺,𝐼𝐼 を計算する。

2. 1.で得られた均質体系の断面積Σ𝐺𝐺,𝐼𝐼 を用いて均質体系を計算する。

3.非均質体系のメッシュ境界の表面中性子束𝜙𝜙𝑠𝑠±,𝑔𝑔,𝑖𝑖 と 2.で得られた均質体系のメッシュ境界の表面中性子束𝜙𝜙𝑠𝑠±,𝐺𝐺,𝐼𝐼 から(10)式を用いて不連続因子𝐷𝐷𝐷𝐷𝑠𝑠±,𝐺𝐺,𝐼𝐼 を計算する。

このように(10)式中の均質体系のメッシュ境界の表面中性子束𝜙𝜙𝑠𝑠±,𝐺𝐺,𝐼𝐼 は、(2)式を用いて得られた均質断面積Σ𝐺𝐺,𝐼𝐼 を用いた均質体系の計算から得られた表面中性子束である。SPH 法では SPH

因子を収束させるために反復計算が必要であったが、不連続因子の計算では SPH 法と異なり、反復計算は必要ない。
ただし、SPH 法の場合と異なり、不連続因子を計算コードに適用する場合には、計算コードを書き換える必要がある。差分式を用いた拡散計算では一般的に、(1)式の左辺第一項の拡散項(漏洩項)を、中性子流を用いて 15−𝐷𝐷𝐺𝐺,𝐼𝐼d2
𝐽𝐽𝑠𝑠+,𝐺𝐺,𝐼𝐼 − 𝐽𝐽𝑠𝑠−,𝐺𝐺,𝐼𝐼
𝜙𝜙𝑔𝑔,𝑖𝑖 =
2
d𝑥𝑥
∆𝑥𝑥𝐼𝐼

(14)

で表す 13)。ここで、𝐽𝐽𝑠𝑠−,𝐺𝐺,𝐼𝐼 は負側のメッシュ境界での表面中性子流、𝐽𝐽𝑠𝑠+,𝐺𝐺,𝐼𝐼 は正側のメッシュ境

界での表面中性子流、∆𝑥𝑥𝐼𝐼 は x 軸方向のメッシュ幅である。表面中性子流を表面中性子束𝜙𝜙𝑠𝑠±,𝐺𝐺,𝐼𝐼とメッシュ中央の中性子束𝜙𝜙𝐺𝐺,𝐼𝐼 を用いて差分式で書くと、次式で書き表せる。
𝜙𝜙𝑠𝑠+,𝐺𝐺,𝐼𝐼 − 𝜙𝜙𝐺𝐺,𝐼𝐼
𝜙𝜙𝐺𝐺,𝐼𝐼+1 − 𝜙𝜙𝑠𝑠−,𝐺𝐺,𝐼𝐼+1
𝐽𝐽𝑠𝑠+,𝐺𝐺,𝐼𝐼 = −𝐷𝐷𝐺𝐺,𝐼𝐼= −𝐷𝐷𝐺𝐺,𝐼𝐼+1∆𝑥𝑥𝐼𝐼∆𝑥𝑥𝐼𝐼+1
2
2
𝜙𝜙𝐺𝐺,𝐼𝐼 − 𝜙𝜙𝑠𝑠−,𝐺𝐺,𝐼𝐼
𝜙𝜙𝑠𝑠+,𝐺𝐺,𝐼𝐼−1 − 𝜙𝜙𝐺𝐺,𝐼𝐼−1
𝐽𝐽𝑠𝑠−,𝐺𝐺,𝐼𝐼 = −𝐷𝐷𝐺𝐺,𝐼𝐼= −𝐷𝐷𝐺𝐺,𝐼𝐼−1∆𝑥𝑥𝐼𝐼∆𝑥𝑥𝐼𝐼−1
2
2

(15)

(16)

(15)、(16)式はメッシュ境界において両側のメッシュで決定される中性子流が一致する(連続である)ことを示しており、さらに不連続因子を適用しない場合には、両側のメッシュ境界の表面
中性子束は連続となる。すなわち、
𝜙𝜙𝑠𝑠−,𝐺𝐺,𝐼𝐼+1 = 𝜙𝜙𝑠𝑠+,𝐺𝐺,𝐼𝐼

𝜙𝜙𝑠𝑠−,𝐺𝐺,𝐼𝐼 = 𝜙𝜙𝑠𝑠+,𝐺𝐺,𝐼𝐼−1

(17)
(18)

が成立する。(17)、(18)式を(15)、(16)式に代入し、表面中性子束を消去すると、(15)、(16)式は𝐽𝐽𝑠𝑠+,𝐺𝐺,𝐼𝐼 = −2𝐷𝐷𝐺𝐺,𝐼𝐼 𝐷𝐷𝐺𝐺,𝐼𝐼+1
2�𝜙𝜙𝐺𝐺,𝐼𝐼+1 − 𝜙𝜙𝐺𝐺,𝐼𝐼 ��𝜙𝜙− 𝜙𝜙𝐺𝐺,𝐼𝐼 � = −∆𝑥𝑥𝐼𝐼+1 ∆𝑥𝑥𝐼𝐼
𝐷𝐷𝐺𝐺,𝐼𝐼 ∆𝑥𝑥𝐼𝐼+1 + 𝐷𝐷𝐺𝐺,𝐼𝐼+1 ∆𝑥𝑥𝐼𝐼 𝐺𝐺,𝐼𝐼+1
𝐷𝐷𝐺𝐺,𝐼𝐼+1 + 𝐷𝐷𝐺𝐺,𝐼𝐼(19)

3.4.1 節でも説明するが、(14)式は差分化しても導出できるが、ガウスの発散定理を用いても導出することが出来る。

15

9

𝐽𝐽𝑠𝑠−,𝐺𝐺,𝐼𝐼 = −となり、(14)式は−𝐷𝐷𝐺𝐺,𝐼𝐼

2𝐷𝐷𝐺𝐺,𝐼𝐼−1 𝐷𝐷𝐺𝐺,𝐼𝐼
2�𝜙𝜙𝐺𝐺,𝐼𝐼 − 𝜙𝜙𝐺𝐺,𝐼𝐼−1 ��𝜙𝜙𝐺𝐺,𝐼𝐼 − 𝜙𝜙𝐺𝐺,𝐼𝐼−1 � = −∆𝑥𝑥𝐼𝐼 ∆𝑥𝑥𝐼𝐼−1
𝐷𝐷𝐺𝐺,𝐼𝐼−1 ∆𝑥𝑥𝐼𝐼 + 𝐷𝐷𝐺𝐺,𝐼𝐼 ∆𝑥𝑥𝐼𝐼−1
𝐷𝐷𝐺𝐺,𝐼𝐼 + 𝐷𝐷𝐺𝐺,𝐼𝐼−1

d2
2𝐷𝐷𝐺𝐺,𝐼𝐼 𝐷𝐷𝐺𝐺,𝐼𝐼+1
𝜙𝜙 = −
�𝜙𝜙𝐺𝐺,𝐼𝐼+1 − 𝜙𝜙𝐺𝐺,𝐼𝐼 �
d𝑥𝑥 2 𝑔𝑔,𝑖𝑖
�𝐷𝐷𝐺𝐺,𝐼𝐼 ∆𝑥𝑥𝐼𝐼+1 + 𝐷𝐷𝐺𝐺,𝐼𝐼+1 ∆𝑥𝑥𝐼𝐼 �∆𝑥𝑥𝐼𝐼+

2𝐷𝐷𝐺𝐺,𝐼𝐼−1 𝐷𝐷𝐺𝐺,𝐼𝐼�𝐷𝐷𝐺𝐺,𝐼𝐼−1 ∆𝑥𝑥𝐼𝐼 + 𝐷𝐷𝐺𝐺,𝐼𝐼 ∆𝑥𝑥𝐼𝐼−1 �∆𝑥𝑥

(20)

(21)

�𝜙𝜙𝐺𝐺,𝐼𝐼 − 𝜙𝜙𝐺𝐺,𝐼𝐼−1 �とメッシュ中央の中性子束のみで書き表すことができる。差分式を用いた拡散計算では一般的に、(1)式に(21)式を代入して計算を行っている 13)。
一方、不連続因子を適用する場合には、(10)式に示したように表面中性子束に不連続性が現れ、その結果、(15)、(16)式を(11)式の関係を用いて次のように書き換える必要がある。
∑𝑖𝑖∈𝐼𝐼 ∑𝑔𝑔∈𝐺𝐺 𝜙𝜙𝑠𝑠+,𝑔𝑔,𝑖𝑖 𝑆𝑆𝑠𝑠+,𝑖𝑖∑𝑖𝑖∈𝐼𝐼 ∑𝑔𝑔∈𝐺𝐺 𝜙𝜙𝑠𝑠−,𝑔𝑔,𝑖𝑖+1 𝑆𝑆𝑠𝑠−,𝑖𝑖+1− ϕ𝐺𝐺,𝐼𝐼ϕ𝐺𝐺,𝐼𝐼+1 −𝐷𝐷𝐷𝐷𝑠𝑠+,𝐺𝐺,𝐼𝐼 𝑆𝑆𝑠𝑠+,𝐼𝐼
𝐷𝐷𝐷𝐷𝑠𝑠−,𝐺𝐺,𝐼𝐼+1 𝑆𝑆𝑠𝑠−,𝐼𝐼+1
𝐽𝐽𝑠𝑠+,𝐺𝐺,𝐼𝐼 = −𝐷𝐷𝐺𝐺,𝐼𝐼= −𝐷𝐷𝐺𝐺,𝐼𝐼+1∆𝑥𝑥𝐼𝐼∆𝑥𝑥𝐼𝐼+1
2
2∑𝑖𝑖∈𝐼𝐼 ∑𝑔𝑔∈𝐺𝐺 𝜙𝜙𝑠𝑠−,𝑔𝑔,𝑖𝑖 𝑆𝑆𝑠𝑠−,𝑖𝑖∑𝑖𝑖∈𝐼𝐼 ∑𝑔𝑔∈𝐺𝐺 𝜙𝜙𝑠𝑠+,𝑔𝑔,𝑖𝑖−1 𝑆𝑆𝑠𝑠+,𝑖𝑖−1
𝜙𝜙𝐺𝐺,𝐼𝐼 −− 𝜙𝜙𝐺𝐺,𝐼𝐼−1
𝐷𝐷𝐷𝐷𝑠𝑠−,𝐺𝐺,𝐼𝐼 𝑆𝑆𝑠𝑠−,𝐼𝐼
𝐷𝐷𝐷𝐷𝑠𝑠+,𝐺𝐺,𝐼𝐼−1 𝑆𝑆𝑠𝑠+,𝐼𝐼−1
𝐽𝐽𝑠𝑠−,𝐺𝐺,𝐼𝐼 = −𝐷𝐷𝐺𝐺,𝐼𝐼= −𝐷𝐷𝐺𝐺,𝐼𝐼−1∆𝑥𝑥𝐼𝐼∆𝑥𝑥𝐼𝐼−1
2
2

(22)

(23)

この場合、両側のメッシュ境界の非均質体系の表面中性子束は連続として計算する必要がある。
すなわち、

� � 𝜙𝜙𝑠𝑠−,𝑔𝑔,𝑖𝑖+1 𝑆𝑆𝑠𝑠−,𝑖𝑖+1 = � � 𝜙𝜙𝑠𝑠+,𝑔𝑔,𝑖𝑖 𝑆𝑆𝑠𝑠+,𝑖𝑖
𝑖𝑖∈𝐼𝐼 𝑔𝑔∈𝐺𝐺� � 𝜙𝜙𝑠𝑠−,𝑔𝑔,𝑖𝑖 𝑆𝑆𝑠𝑠−,𝑖𝑖 = � � 𝜙𝜙𝑠𝑠+,𝑔𝑔,𝑖𝑖−1 𝑆𝑆𝑠𝑠+,𝑖𝑖−1
𝑖𝑖∈𝐼𝐼 𝑔𝑔∈𝐺𝐺(24)𝑖𝑖∈𝐼𝐼 𝑔𝑔∈𝐺𝐺(25)

𝑖𝑖∈𝐼𝐼 𝑔𝑔∈𝐺𝐺が成立する。(24)、(25)式を(22)、(23)式に代入し、非均質体系の表面中性子束を消去すると、(22)、(23)式は𝐽𝐽𝑠𝑠+,𝐺𝐺,𝐼𝐼 = −2𝐷𝐷𝐺𝐺,𝐼𝐼 𝐷𝐷𝐺𝐺,𝐼𝐼+1

𝐷𝐷𝐷𝐷𝑠𝑠−,𝐺𝐺,𝐼𝐼+1 𝜙𝜙𝐺𝐺,𝐼𝐼+1 − 𝐷𝐷𝐷𝐷𝑠𝑠+,𝐺𝐺,𝐼𝐼 𝜙𝜙𝐺𝐺,𝐼𝐼
𝐷𝐷𝐺𝐺,𝐼𝐼 𝐷𝐷𝐷𝐷𝑠𝑠−,𝐺𝐺,𝐼𝐼+1 ∆𝑥𝑥𝐼𝐼+1 + 𝐷𝐷𝐺𝐺,𝐼𝐼+1 𝐷𝐷𝐷𝐷𝑠𝑠+,𝐺𝐺,𝐼𝐼 ∆𝑥𝑥𝐼𝐼

𝐷𝐷𝐷𝐷𝑠𝑠−,𝐺𝐺,𝐼𝐼+1 𝜙𝜙𝐺𝐺,𝐼𝐼+1 − 𝐷𝐷𝐷𝐷𝑠𝑠+,𝐺𝐺,𝐼𝐼 𝜙𝜙𝐺𝐺,𝐼𝐼= −2∆𝑥𝑥𝐼𝐼+1∆𝑥𝑥
𝐷𝐷𝐷𝐷+ 𝐼𝐼 𝐷𝐷𝐷𝐷
𝐷𝐷𝐺𝐺,𝐼𝐼+1 𝑠𝑠−,𝐺𝐺,𝐼𝐼+1 𝐷𝐷𝐺𝐺,𝐼𝐼 𝑠𝑠+,𝐺𝐺,𝐼𝐼(27)

𝐽𝐽𝑠𝑠−,𝐺𝐺,𝐼𝐼 = −2𝐷𝐷𝐺𝐺,𝐼𝐼−1 𝐷𝐷𝐺𝐺,𝐼𝐼

𝐷𝐷𝐷𝐷𝑠𝑠−,𝐺𝐺,𝐼𝐼 𝜙𝜙𝐺𝐺,𝐼𝐼 − 𝐷𝐷𝐷𝐷𝑠𝑠+,𝐺𝐺,𝐼𝐼−1 𝜙𝜙𝐺𝐺,𝐼𝐼−1
𝐷𝐷𝐺𝐺,𝐼𝐼−1 𝐷𝐷𝐷𝐷𝑠𝑠−,𝐺𝐺,𝐼𝐼 ∆𝑥𝑥𝐼𝐼 + 𝐷𝐷𝐺𝐺,𝐼𝐼 𝐷𝐷𝐷𝐷𝑠𝑠+,𝐺𝐺,𝐼𝐼−1 ∆𝑥𝑥𝐼𝐼−1
𝐷𝐷𝐷𝐷𝑠𝑠−,𝐺𝐺,𝐼𝐼 𝜙𝜙𝐺𝐺,𝐼𝐼 − 𝐷𝐷𝐷𝐷𝑠𝑠+,𝐺𝐺,𝐼𝐼−1 𝜙𝜙𝐺𝐺,𝐼𝐼−1∆𝑥𝑥𝐼𝐼∆𝑥𝑥𝐼𝐼−1
𝐷𝐷𝐺𝐺,𝐼𝐼 𝐷𝐷𝐷𝐷𝑠𝑠−,𝐺𝐺,𝐼𝐼 + 𝐷𝐷𝐺𝐺,𝐼𝐼−1 𝐷𝐷𝐷𝐷𝑠𝑠+,𝐺𝐺,𝐼𝐼−1

(28)

= −2
となり、(14)式は10

−𝐷𝐷𝐺𝐺,𝐼𝐼となる。

d2
𝐷𝐷𝐷𝐷𝑠𝑠−,𝐺𝐺,𝐼𝐼+1 𝜙𝜙𝐺𝐺,𝐼𝐼+1 − 𝐷𝐷𝐷𝐷𝑠𝑠+,𝐺𝐺,𝐼𝐼 𝜙𝜙𝐺𝐺,𝐼𝐼
𝜙𝜙 = −2𝐷𝐷𝐺𝐺,𝐼𝐼 𝐷𝐷𝐺𝐺,𝐼𝐼+1
d𝑥𝑥 2 𝑔𝑔,𝑖𝑖
𝐷𝐷𝐺𝐺,𝐼𝐼 𝐷𝐷𝐷𝐷𝑠𝑠−,𝐺𝐺,𝐼𝐼+1 ∆𝑥𝑥𝐼𝐼+1 + 𝐷𝐷𝐺𝐺,𝐼𝐼+1 𝐷𝐷𝐷𝐷𝑠𝑠+,𝐺𝐺,𝐼𝐼 ∆𝑥𝑥𝐼𝐼
𝐷𝐷𝐷𝐷𝑠𝑠−,𝐺𝐺,𝐼𝐼 𝜙𝜙𝐺𝐺,𝐼𝐼 − 𝐷𝐷𝐷𝐷𝑠𝑠+,𝐺𝐺,𝐼𝐼−1 𝜙𝜙𝐺𝐺,𝐼𝐼−1+ 2𝐷𝐷𝐺𝐺,𝐼𝐼−1 𝐷𝐷𝐺𝐺,𝐼𝐼
𝐷𝐷𝐺𝐺,𝐼𝐼−1 𝐷𝐷𝐷𝐷𝑠𝑠−,𝐺𝐺,𝐼𝐼 ∆𝑥𝑥𝐼𝐼 + 𝐷𝐷𝐺𝐺,𝐼𝐼 𝐷𝐷𝐷𝐷𝑠𝑠+,𝐺𝐺,𝐼𝐼−1 ∆𝑥𝑥𝐼𝐼−1(29)

このように不連続因子では、メッシュ境界の表面中性子束𝜙𝜙𝑠𝑠±,𝐺𝐺,𝐼𝐼 の不連続性を不連続因子

𝐷𝐷𝐷𝐷𝑠𝑠±,𝐺𝐺,𝐼𝐼 によって表すことで非均質体系のメッシュ境界の表面中性子束の連続性を再現し、空間均質化誤差を低減している。また、(21)式と(29)式を比べると、分母と分子に不連続因子が含まれているものの、式自体はほとんど同じであることが分かる。そのため、前述の通り不連続因子を計算コードに適用する場合には、計算コードを書き換える必要があるものの、不連続因子の適用に伴う計算コードの書き換えはそれほど大きな手間ではない。
2.2.3 SPH 因子と不連続因子の違い
SPH 法では断面積に SPH 因子を掛けるだけで計算コードを書き換える必要がなかったため、拡散計算だけでなく MOC や衝突確率法などの輸送計算にも適用することが可能であった。しかし、輸送計算に対する不連続因子は現時点では理論的に確立されておらず、不連続因子の適用範囲は拡散計算に限定されている 16。
また、(10)式のように不連続因子は各メッシュ境界で異なるため、三次元体系では一つの群、一つのメッシュ当たり x±、y±、z±と 6 個の不連続因子が必要となる。SPH 法では一つの群、一つのメッシュ当たりの SPH 因子は 1 個であり、不連続因子は SPH 法に比べて自由度が大きい。
しかし、不連続因子では一つのメッシュで 6 個×エネルギー群数のデータが追加で必要となる。
そのため、不連続因子の適用には、不連続因子を適用することで増加するメモリー量が許容範
囲かどうかを考える必要がある 17。
なお、SPH 法の場合と同様に不連続因子も不連続因子を作成するのに用いた単一集合体体系に近い中性子流の場合には高い計算精度を有するが、UO2 燃料集合体と MOX 燃料集合体が隣接する境界など、単一集合体体系と大きく異なる中性子流の場合には計算精度が悪化する場合があるので注意する必要がある。
2.2.4 SPH 因子と不連続因子の関係
SPH 因子と不連続因子の関係を考える。ここで、メッシュ表面の不連続因子が全て等しい場合を考える。この場合、メッシュ表面の中性子束が全て同じ割合で変化していることから、メッシュ内の中性子束分布は不連続因子の分だけ平行移動していることになる。これはメッシュ
ただし、近年は輸送計算や SP3 計算に対する不連続因子の適用に関する研究が報告されており 14,15)、この不連続因子の適用範囲の狭さの問題は克服されつつある。
17 例えば 764 体の BWR-5 の 1/4 炉心を pin-by-pin で解析する場合、11×11 メッシュ/集合体(9×9 燃料＋外側水領域一層)で、軸方向に 26 ノード、エネルギー群を 8 群とすると、
11×11 [メッシュ/集合体]×26 [ノード]×8 [群]×6 [表面/メッシュ]×764 [体/炉心]×4 [byte/データ]÷4= 115MB となる。一般的な炉心計算では 1 メッシュ/集合体であり、その場合に必要なメモリ
ー量は約 1MB であり、現在の計算機性能では大きな問題にはならない。
1611

平均の中性子束も同様であり、
� � 𝜙𝜙𝑔𝑔,𝑖𝑖 𝑉𝑉𝑖𝑖 = 𝐷𝐷𝐷𝐷𝐺𝐺,𝐼𝐼 𝜙𝜙𝐺𝐺,𝐼𝐼 𝑉𝑉𝐺𝐺,𝐼𝐼(30)

𝑔𝑔∈𝐺𝐺 𝑖𝑖∈𝐼𝐼となる。(30)式を整理すると𝐷𝐷𝐷𝐷𝐺𝐺,𝐼𝐼 =∑𝑔𝑔∈𝐺𝐺 ∑𝑖𝑖∈𝐼𝐼 𝜙𝜙𝑔𝑔,𝑖𝑖 𝑉𝑉𝑖𝑖
= 𝜇𝜇𝐺𝐺,𝐼𝐼
𝜙𝜙𝐺𝐺,𝐼𝐼 𝑉𝑉𝐺𝐺,𝐼𝐼(31)となる。このように、メッシュ表面の不連続因子が全て等しい場合、SPH 因子と不連続因子は等しくなる。不連続因子は均質体系と非均質体系の各メッシュ境界の中性子流を一致させていることを考えると、SPH 法は均質体系と非均質体系のメッシュ平均の中性子流を一致させていることになる。
前述したように、SPH 法は UO2 燃料集合体と MOX 燃料集合体が隣接する境界など、SPH 因子を作成した際の非均質体系と中性子流が大きく条件が異なる場合には計算精度が悪化する場合がある。これは不連続因子も同様であるが、メッシュ平均の中性子流のみを一致している
SPH 法の方がその影響が大きい。そのため、メッシュ境界の中性子流が大きく変わる場合には、メッシュ平均の中性子流ではなく、個々のメッシュ境界の中性子流を一致させる不連続因子の方が適している 18。
空間均質化誤差低減手法を適用する場合には、SPH 法や不連続因子の特徴を把握した上で、最適なものを選択することが重要である。

3. 高速化手法について
3.1 外部反復と内部反復(1)式の拡散方程式を解くことを考える。(1)式は固有値問題と呼ばれ 19、実効増倍率が固有値、中性子束が固有関数(固有ベクトル)となる。なお、固有値は重解、虚数解も含めると固有ベクトルの数だけ存在するが、一般的な核計算コードでは最大固有値である実効増倍率のみを計算している 20。
核計算コードでは、実効増倍率も含めた核分裂中性子源を固定して中性子束を計算する内部
反復と、実効増倍率を更新する外部反復の二つの反復計算で構成されている 21。例えば(1)式の拡散方程式を解く場合、(32)式のように右辺の核分裂源項と散乱源項の和を中性子源項𝑄𝑄𝑔𝑔,𝑖𝑖 として中性子束𝜙𝜙𝑔𝑔,𝑖𝑖 を計算する内部反復と、得られた中性子束を用いて(33)式を用いて実効増倍率 k
SPH 法と不連続因子の双方の特徴を生かし、pin-by-pin 炉心計算において、集合体境界については不連続因子を適用し、集合体内については SPH 法を適用する不連続因子と SPH 法を組み合わせた手法も提案されている 16)。
19 本資料では簡単のため、拡散方程式で説明するが、当然のことながら輸送計算も固有値問題である。
20 定常状態であれば最大固有値の評価のみで十分であるが、過渡時や未臨界の場合にはより高次の固有値・固有ベクトルが必要となる場合がある。
21 固定源計算では、実効増倍率が存在しないことから、外部反復を必要とせず、内部反復のみで計算することができる。
18

12

を更新する外部反復を実効増倍率が収束するまで交互に行う必要がある 22。
−𝐷𝐷𝑔𝑔,𝑖𝑖𝜒𝜒𝑔𝑔,𝑖𝑖d2
𝜙𝜙 + Σ𝑡𝑡,𝑔𝑔,𝑖𝑖 𝜙𝜙𝑔𝑔,𝑖𝑖 =� 𝜈𝜈Σ ′ 𝜙𝜙 ′ + � Σ𝑠𝑠,𝑔𝑔′ →𝑔𝑔,𝑖𝑖 𝜙𝜙𝑔𝑔′ ,𝑖𝑖 = 𝑄𝑄𝑔𝑔,𝑖𝑖
d𝑥𝑥 2 𝑔𝑔,𝑖𝑖
𝑘𝑘 ′ 𝑓𝑓,𝑔𝑔 ,𝑖𝑖 𝑔𝑔 ,𝑖𝑖′

𝑘𝑘 = � 𝑉𝑉𝑖𝑖 ′ � 𝜈𝜈Σ𝑓𝑓,𝑔𝑔′ ,𝑖𝑖 ′ 𝜙𝜙𝑔𝑔′ ,𝑖𝑖′
𝑖𝑖 ′

𝑔𝑔

𝑔𝑔

𝑔𝑔′

(32)

(33)

このように実効増倍率を計算するためには、内部反復と外部反復の二つの計算があり、それぞれに高速化手法が存在する。本資料では内部反復と外部反復のそれぞれの反復計算で用いられている高速化手法について説明する。
3.2 内部反復に用いられる高速化手法
前節で説明したように、内部反復では実効増倍率も含めた核分裂中性子源を固定して中性子束を計算する。中性子束の計算は計算手法によって大きく異なり、例えば MOC では(34)式に示す一次元微分方程式を解く 17)。
d𝜙𝜙𝑔𝑔,𝑖𝑖 (𝑠𝑠)+ Σ𝑡𝑡,𝑔𝑔,𝑖𝑖 𝜙𝜙𝑔𝑔,𝑖𝑖 (𝑠𝑠) = 𝑄𝑄𝑔𝑔,𝑖𝑖
d𝑠𝑠

(34)

ここで、s はメッシュ i の境界からの距離である。(34)式を解くにはメッシュ i の片側境界の中性子束を境界条件として用いて一階の微分方程式を解けばよく、(32)式を解く場合とは大きく異なっている。また、(32)式の二階の微分方程式も右辺の中性子源を定数、もしくは多項式関数と見なせばメッシュ i の境界両側の中性子束を境界条件として用いることで解析的に計算することが出来る 18)。
このように内部反復は中性子束の計算手法として様々な手法を採用することができるが、次の 3.2.1 節では高速化手法の説明として理解しやすいことから、(32)式を連立一次方程式として解く場合に有効な高速化手法について説明する。ここで紹介する手法は今では利用されていない手法も多いが、これらの手法に用いられている考え方は様々なところで応用することができる。
3.2.1 連立一次方程式で用いられる高速化手法(32)式を連立一次方程式として解くことを考える。三次元拡散方程式を差分化すると 23、次のように書き表せる 13)。

(33)式を用いることで、(32)式では核分裂中性子源の総和が 1.0 となるように規格化されている。
23 (35)式は、セル中央の中性子束がセル平均の中性子束と等しいという近似を用いて導出している。なお、(35)式の導出については本資料の 2.2.2 節や参考文献 13 を参考にして欲しい。また本資料では簡単のため、差分式で説明しているが、実際の炉心計算コードでは差分法ではな
く、セル内の中性子束分布を考慮した近代ノード法 10,11)を用いている場合が多い。なお、近代
ノード法と差分法は取り扱っている計算式が異なるだけであり、以降で説明する高速化手法は近代ノード法を用いた場合でも有効である。
22

13

𝑦𝑦−

𝑥𝑥−
𝑧𝑧−
0
𝐴𝐴𝑔𝑔,𝑖𝑖,𝑗𝑗,𝑘𝑘
𝜙𝜙𝑔𝑔,𝑖𝑖−1,𝑗𝑗,𝑘𝑘 + 𝐴𝐴𝑔𝑔,𝑖𝑖,𝑗𝑗,𝑘𝑘 𝜙𝜙𝑔𝑔,𝑖𝑖,𝑗𝑗−1,𝑘𝑘 + 𝐴𝐴𝑔𝑔,𝑖𝑖,𝑗𝑗,𝑘𝑘
𝜙𝜙𝑔𝑔,𝑖𝑖,𝑗𝑗,𝑘𝑘−1 + 𝐴𝐴𝑔𝑔,𝑖𝑖,𝑗𝑗,𝑘𝑘
𝜙𝜙𝑔𝑔,𝑖𝑖,𝑗𝑗,𝑘𝑘
𝑦𝑦+

𝑥𝑥+
𝑧𝑧+
+ 𝐴𝐴𝑔𝑔,𝑖𝑖,𝑗𝑗,𝑘𝑘
𝜙𝜙𝑔𝑔,𝑖𝑖+1,𝑗𝑗,𝑘𝑘 + 𝐴𝐴𝑔𝑔,𝑖𝑖,𝑗𝑗,𝑘𝑘 𝜙𝜙𝑔𝑔,𝑖𝑖,𝑗𝑗+1,𝑘𝑘 + 𝐴𝐴𝑔𝑔,𝑖𝑖,𝑗𝑗,𝑘𝑘
𝜙𝜙𝑔𝑔,𝑖𝑖,𝑗𝑗,𝑘𝑘+1
𝜒𝜒𝑔𝑔,𝑖𝑖,𝑗𝑗,𝑘𝑘=� 𝜈𝜈Σ𝑓𝑓,𝑔𝑔′ ,𝑖𝑖,𝑗𝑗,𝑘𝑘 𝜙𝜙𝑔𝑔′ ,𝑖𝑖,𝑗𝑗,𝑘𝑘 + � Σ𝑠𝑠,𝑔𝑔′ →𝑔𝑔,𝑖𝑖,𝑗𝑗,𝑘𝑘 𝜙𝜙𝑔𝑔′ ,𝑖𝑖,𝑗𝑗,𝑘𝑘
𝑘𝑘
′
′
𝑔𝑔

ここで、
𝑥𝑥±
𝐴𝐴𝑔𝑔,𝑖𝑖,𝑗𝑗,𝑘𝑘=−𝑦𝑦±

𝐴𝐴𝑔𝑔,𝑖𝑖,𝑗𝑗,𝑘𝑘 = −𝑧𝑧±
𝐴𝐴𝑔𝑔,𝑖𝑖,𝑗𝑗,𝑘𝑘
=−

(35)

𝑔𝑔

2𝐷𝐷𝑔𝑔,𝑖𝑖±1,𝑗𝑗,𝑘𝑘 𝐷𝐷𝑔𝑔,𝑖𝑖,𝑗𝑗,𝑘𝑘(36)

�𝐷𝐷𝑔𝑔,𝑖𝑖±1,𝑗𝑗,𝑘𝑘 ∆𝑥𝑥𝑖𝑖 + 𝐷𝐷𝑔𝑔,𝑖𝑖,𝑗𝑗,𝑘𝑘 ∆𝑥𝑥𝑖𝑖±1 �∆𝑥𝑥𝑖𝑖
2𝐷𝐷𝑔𝑔,𝑖𝑖,𝑗𝑗±1,𝑘𝑘 𝐷𝐷𝑔𝑔,𝑖𝑖,𝑗𝑗,𝑘𝑘(37)

�𝐷𝐷𝑔𝑔,𝑖𝑖,𝑗𝑗±1,𝑘𝑘 ∆𝑦𝑦𝑗𝑗 + 𝐷𝐷𝑔𝑔,𝑖𝑖,𝑗𝑗,𝑘𝑘 ∆𝑦𝑦𝑗𝑗±1 �∆𝑦𝑦𝑗𝑗
2𝐷𝐷𝑔𝑔,𝑖𝑖,𝑗𝑗,𝑘𝑘±1 𝐷𝐷𝑔𝑔,𝑖𝑖,𝑗𝑗,𝑘𝑘�𝐷𝐷𝑔𝑔,𝑖𝑖,𝑗𝑗,𝑘𝑘±1 ∆𝑧𝑧𝑘𝑘 + 𝐷𝐷𝑔𝑔,𝑖𝑖,𝑗𝑗,𝑘𝑘 ∆𝑧𝑧𝑘𝑘±1 �∆𝑧𝑧𝑘𝑘
𝑦𝑦−

(38)
𝑦𝑦+

0
𝑥𝑥−
𝑥𝑥+
𝑧𝑧−
𝑧𝑧+
𝐴𝐴𝑔𝑔,𝑖𝑖,𝑗𝑗,𝑘𝑘= Σ𝑡𝑡,𝑔𝑔,𝑖𝑖,𝑗𝑗,𝑘𝑘 − �𝐴𝐴𝑔𝑔,𝑖𝑖,𝑗𝑗,𝑘𝑘+ 𝐴𝐴𝑔𝑔,𝑖𝑖,𝑗𝑗,𝑘𝑘+ 𝐴𝐴𝑔𝑔,𝑖𝑖,𝑗𝑗,𝑘𝑘 + 𝐴𝐴𝑔𝑔,𝑖𝑖,𝑗𝑗,𝑘𝑘 + 𝐴𝐴𝑔𝑔,𝑖𝑖,𝑗𝑗,𝑘𝑘+ 𝐴𝐴𝑔𝑔,𝑖𝑖,𝑗𝑗,𝑘𝑘�(39)

であり、∆𝑥𝑥𝑖𝑖 、∆𝑦𝑦𝑗𝑗 、∆𝑧𝑧𝑘𝑘 はそれぞれ x 軸方向、y 軸方向、z 軸方向のメッシュ幅である。一般的

な拡散計算コードでは(35)式を用いて計算を行っている。しかし、今回は内部反復の高速化手法の例として用いることから、式を簡略化するため、一次元拡散方程式を取り扱う。一次元拡散方程式は(35)式中の z 軸と y 軸の項を削除すれば求められ、𝜒𝜒𝑔𝑔,𝑖𝑖𝑥𝑥−
0𝑥𝑥+
𝐴𝐴𝑔𝑔,𝑖𝑖
𝜙𝜙𝑔𝑔,𝑖𝑖−1 + 𝐴𝐴𝑔𝑔,𝑖𝑖
𝜙𝜙𝑔𝑔,𝑖𝑖 + 𝐴𝐴𝑔𝑔,𝑖𝑖
𝜙𝜙𝑔𝑔,𝑖𝑖+1 =� 𝜈𝜈Σ ′ 𝜙𝜙 ′ + � Σ𝑠𝑠,𝑔𝑔′ →𝑔𝑔,𝑖𝑖 𝜙𝜙𝑔𝑔′ ,𝑖𝑖
𝑘𝑘 ′ 𝑓𝑓,𝑔𝑔 ,𝑖𝑖 𝑔𝑔 ,𝑖𝑖′𝑔𝑔となる。ここで、

(40)

𝑔𝑔

0
𝑥𝑥−
𝑥𝑥+
𝐴𝐴𝑔𝑔,𝑖𝑖= Σ𝑡𝑡,𝑔𝑔,𝑖𝑖 − �𝐴𝐴𝑔𝑔,𝑖𝑖+ 𝐴𝐴𝑔𝑔,𝑖𝑖�(41)である。(40)式の変数は当然ながら中性子束𝜙𝜙𝑔𝑔,𝑖𝑖 であり、(40)式はエネルギー群数×メッシュ数の数だけ用意される。(40)式を(32)式のように右辺を中性子源項𝑄𝑄𝑔𝑔,𝑖𝑖 とおくと
𝑥𝑥−
0
𝑥𝑥+
𝐴𝐴𝑔𝑔,𝑖𝑖
𝜙𝜙𝑔𝑔,𝑖𝑖−1 + 𝐴𝐴𝑔𝑔,𝑖𝑖
𝜙𝜙𝑔𝑔,𝑖𝑖 + 𝐴𝐴𝑔𝑔,𝑖𝑖
𝜙𝜙𝑔𝑔,𝑖𝑖+1 = 𝑄𝑄𝑔𝑔,𝑖𝑖(42)となる。(42)式をエネルギー群数×メッシュ数分並べると、
0𝑥𝑥+⋯ + 𝐴𝐴𝑔𝑔,𝑖𝑖−2𝜙𝜙𝑔𝑔,𝑖𝑖−2 + 𝐴𝐴𝑔𝑔,𝑖𝑖−2𝜙𝜙𝑔𝑔,𝑖𝑖−1 + 0 ∙ 𝜙𝜙𝑔𝑔,𝑖𝑖+ 0 ∙ 𝜙𝜙𝑔𝑔,𝑖𝑖+1

𝑥𝑥−
0
𝑥𝑥+
⋯ + 𝐴𝐴𝑔𝑔,𝑖𝑖−1
𝜙𝜙𝑔𝑔,𝑖𝑖−2 + 𝐴𝐴𝑔𝑔,𝑖𝑖−1
𝜙𝜙𝑔𝑔,𝑖𝑖−1 + 𝐴𝐴𝑔𝑔,𝑖𝑖−1
𝜙𝜙𝑔𝑔,𝑖𝑖 + 0 ∙ 𝜙𝜙𝑔𝑔,𝑖𝑖+1⋯ + 0 ∙ 𝜙𝜙𝑔𝑔,𝑖𝑖−2⋯ + 0 ∙ 𝜙𝜙𝑔𝑔,𝑖𝑖−2⋯ + 0 ∙ 𝜙𝜙𝑔𝑔,𝑖𝑖−2⋯ + 0 ∙ 𝜙𝜙𝑔𝑔,𝑖𝑖−2𝑥𝑥−+ 𝐴𝐴𝑔𝑔,𝑖𝑖
𝜙𝜙𝑔𝑔,𝑖𝑖−1+ 0 ∙ 𝜙𝜙𝑔𝑔,𝑖𝑖−1+ 0 ∙ 𝜙𝜙𝑔𝑔,𝑖𝑖−1+ 0 ∙ 𝜙𝜙𝑔𝑔,𝑖𝑖−1+ 0 ∙ 𝜙𝜙𝑔𝑔,𝑖𝑖+2+ 0 ∙ 𝜙𝜙𝑔𝑔,𝑖𝑖+20𝑥𝑥++ 𝐴𝐴𝑔𝑔,𝑖𝑖
𝜙𝜙𝑔𝑔,𝑖𝑖 + 𝐴𝐴𝑔𝑔,𝑖𝑖
𝜙𝜙𝑔𝑔,𝑖𝑖+1 + 0 ∙ 𝜙𝜙𝑔𝑔,𝑖𝑖+2
𝑥𝑥−
0
𝑥𝑥+
+ 𝐴𝐴𝑔𝑔,𝑖𝑖
𝜙𝜙𝑔𝑔,𝑖𝑖+1 + 𝐴𝐴𝑔𝑔,𝑖𝑖+1
𝜙𝜙𝑔𝑔,𝑖𝑖+1 + 𝐴𝐴𝑔𝑔,𝑖𝑖+1
𝜙𝜙𝑔𝑔,𝑖𝑖+2+ 0 ∙ 𝜙𝜙𝑔𝑔,𝑖𝑖+ 0 ∙ 𝜙𝜙𝑔𝑔,𝑖𝑖+ ⋯ = 𝑄𝑄𝑔𝑔,𝑖𝑖−2+ ⋯ = 𝑄𝑄𝑔𝑔,𝑖𝑖−1+ ⋯ = 𝑄𝑄𝑔𝑔,𝑖𝑖+ ⋯ = 𝑄𝑄𝑔𝑔,𝑖𝑖+1

(43)

𝑥𝑥−
0
+ 𝐴𝐴𝑔𝑔,𝑖𝑖+2
𝜙𝜙𝑔𝑔,𝑖𝑖+1 + 𝐴𝐴𝑔𝑔,𝑖𝑖+2∙ 𝜙𝜙𝑔𝑔,𝑖𝑖+2 + ⋯ = 𝑄𝑄𝑔𝑔,𝑖𝑖+2+ 0 ∙ 𝜙𝜙𝑔𝑔,𝑖𝑖+1

𝑥𝑥−+ 𝐴𝐴𝑔𝑔,𝑖𝑖+3𝜙𝜙𝑔𝑔,𝑖𝑖+2 + ⋯ = 𝑄𝑄𝑔𝑔,𝑖𝑖+3といったようにメッシュ番号が一つ増える毎に係数𝐴𝐴の位置が右に一つずれていく。また、零でない係数は𝐴𝐴𝑥𝑥− 、𝐴𝐴0 、𝐴𝐴𝑥𝑥+ の三つのみであり、残りの係数は全てゼロとなる。このような連

立一次方程式を三点階差の連立一次方程式と呼ぶ 24。以降ではこの(43)式で記述される三点階差の連立一次方程式の解法について説明する 25。
二次元拡散方程式では、𝐴𝐴 𝑦𝑦− 、𝐴𝐴 𝑦𝑦+ の項が追加されて五点階差の連立一次方程式に、三次元拡散方程式ではさらに𝐴𝐴𝑧𝑧− 、𝐴𝐴𝑧𝑧+ の項が追加されて七点階差の連立一次方程式になる。
25 本資料で説明する連立一次方程式の解法は参考文献 19 に詳しい記述がある。また、本資料では簡単のため、三点階差の連立方程式を取り上げているが、五点階差、七点階差の連立方程
24

14

3.2.1.1 ガウスの消去法
連立一次方程式を解く方法として最も簡単なのは、ガウスの消去法である。例えば以下の三
元連立一次方程式を手計算で解くことを考える。
2𝑥𝑥 + 5𝑦𝑦 + 𝑧𝑧 = 12(44)

3𝑥𝑥 + 3𝑦𝑦 + 2𝑧𝑧 = 15(45)

4𝑥𝑥 + 𝑦𝑦 + 𝑧𝑧 = 12

(46)

まず、z を消去するため、(44)式-(46)式、(45)式－2×(46)式とすると、−2𝑥𝑥 + 4𝑦𝑦 = 0(47)−5𝑥𝑥 + 𝑦𝑦 = −9(48)となる。さらに y を消去するため(47)式－4×(48)式とすると18𝑥𝑥 = 36(49)となり、x が求められる。得られた x を(47)式もしくは(48)式に代入すれば y が得られ、得られた x と y を(44)式～(46)式に代入すれば z が得られる。このように手計算で連立一次方程式を解く手法をガウスの消去法と呼ぶ。
連立一次方程式を一般化し、行列形式で書き表すと次のようになる。
𝑨𝑨𝑥𝑥⃗ = 𝑏𝑏�⃗

(50)

ここで、

𝑎𝑎11 ⋯ 𝑎𝑎1𝑛𝑛⋱⋮ �
𝑨𝑨 = � ⋮𝑎𝑎𝑛𝑛1 ⋯ 𝑎𝑎𝑛𝑛𝑛𝑛
𝑥𝑥1𝑥𝑥⃗ = � ⋮ �
𝑥𝑥𝑛𝑛
𝑏𝑏1
𝑏𝑏�⃗ = � ⋮ �𝑏𝑏𝑛𝑛である。先ほどと同様に(50)式を解くため、𝑥𝑥1 を消去すると、
𝑎𝑎11 𝑎𝑎12⎡(2)
0 𝑎𝑎22⎢(2)⎢ 0 𝑎𝑎32⎢(2)⎢ 0 𝑎𝑎42⋮⎢ ⋮(2)⎣ 0 𝑎𝑎𝑛𝑛2となる。ここで、
(2)

𝑎𝑎13

(2)
𝑎𝑎23
(2)
𝑎𝑎33
(2)
𝑎𝑎43

⋮

(2)
𝑎𝑎𝑛𝑛3

𝑎𝑎𝑚𝑚𝑚𝑚 = 𝑎𝑎𝑚𝑚𝑚𝑚 −(2)

𝑏𝑏𝑛𝑛 = 𝑏𝑏𝑛𝑛 −

𝑎𝑎14

(2)
𝑎𝑎24
(2)
𝑎𝑎34
(2)
𝑎𝑎44

⋮

(2)
𝑎𝑎𝑛𝑛4

⋯

(51)

(52)

(53)

𝑎𝑎1𝑛𝑛

𝑏𝑏1
𝑥𝑥1
(2)
⎛𝑏𝑏2 ⎞⎥ 𝑥𝑥2(2)⋯ 𝑎𝑎3𝑛𝑛 ⎥ ⎛𝑥𝑥3 ⎞ ⎜𝑏𝑏3(2) ⎟⎟=⎜(2) ⎥ ⎜ 𝑥𝑥 ⎟⋯ 𝑎𝑎4𝑛𝑛 ⎥ ⎜ 4 ⎟ ⎜𝑏𝑏4(2) ⎟⋮⎜⎟⋱⋮ ⎥ ⎝𝑥𝑥 ⎠
⋮
𝑛𝑛
(2)
(2)
⋯ 𝑎𝑎𝑛𝑛𝑛𝑛 ⎦⎝𝑏𝑏𝑛𝑛 ⎠(2) ⎤⋯ 𝑎𝑎2𝑛𝑛𝑎𝑎𝑚𝑚1
𝑎𝑎𝑎𝑎11 1𝑛𝑛

𝑎𝑎𝑛𝑛1𝑏𝑏
𝑎𝑎11 1

(54)

(55)

(56)

である。同様の操作を繰り返し、𝑥𝑥2 、𝑥𝑥3 、𝑥𝑥4 …を消去していくと(54)式は式でも本資料で紹介する手法で解くことができる。
15

𝑎𝑎11
⎡
0
⎢
⎢ 0
⎢ 0
⎢
⎢ ⋮
⎣ 0

𝑎𝑎12

(2)
𝑎𝑎22

0
0
⋮
0

(2)
𝑎𝑎23
(3)
𝑎𝑎33

となる。ここで、
(𝑙𝑙)

𝑎𝑎13

(𝑙𝑙−1)

0
⋮
0

𝑎𝑎14

(2)
𝑎𝑎24
(3)
𝑎𝑎34
(4)
𝑎𝑎44

⋮
0

⋯

𝑎𝑎1𝑛𝑛

𝑏𝑏1
𝑥𝑥1
(2)
⎛ 𝑏𝑏2 ⎞
⎥ 𝑥𝑥2
(3)
⋯ 𝑎𝑎3𝑛𝑛 ⎥ ⎛𝑥𝑥3 ⎞ ⎜ 𝑏𝑏3(3) ⎟⎟=⎜(4) ⎥ ⎜ 𝑥𝑥 ⎟⋯ 𝑎𝑎4𝑛𝑛 ⎥ ⎜ 4 ⎟ ⎜ 𝑏𝑏4(4) ⎟⋮⎜⎟⋱⋮ ⎥ ⎝𝑥𝑥 ⎠
⋮
𝑛𝑛
(𝑛𝑛)
(𝑛𝑛)
⋯ 𝑎𝑎𝑛𝑛𝑛𝑛 ⎦⎝𝑏𝑏𝑛𝑛 ⎠(2) ⎤⋯ 𝑎𝑎2𝑛𝑛(𝑙𝑙−1)

𝑎𝑎𝑚𝑚,𝑙𝑙−1(𝑙𝑙−1)
𝑎𝑎𝑙𝑙−1,𝑛𝑛(𝑙𝑙−1)
𝑎𝑎𝑙𝑙−1,𝑙𝑙−1(𝑙𝑙−1)
𝑎𝑎𝑛𝑛,𝑙𝑙−1 (𝑙𝑙−1)(𝑙𝑙)(𝑙𝑙−1)
𝑏𝑏𝑛𝑛 = 𝑏𝑏𝑛𝑛
− (𝑙𝑙−1) 𝑏𝑏𝑙𝑙−1
𝑎𝑎𝑙𝑙−1,𝑙𝑙−1

𝑎𝑎𝑚𝑚𝑚𝑚 = 𝑎𝑎𝑚𝑚𝑚𝑚 −

(57)

(58)

(59)

である。ガウスの消去法は数式上では、(57)式のように行列 A を対角線から上の領域にのみ係数のある上三角行列 A’に変形することに等しい。(57)式が得られたら、後は手計算の場合と同様に𝑥𝑥𝑛𝑛 から順番に計算していけばよい。

以上が一般的な連立一次方程式の解法の説明である。それではこの連立一次方程式を早く解くにはどうすればよいかを考えよう。核計算で取り扱う連立一次方程式には次の特徴がある。
・ 左辺の係数行列 A は変化しない
・ 大規模疎行列である。
(40)式の左辺の係数行列 A は𝐴𝐴𝑥𝑥−、𝐴𝐴0、𝐴𝐴𝑥𝑥+ で構成されている。(36)～(39)式の通り、これら

はメッシュ幅と拡散係数、巨視的断面積からなる係数である。そのため、内部反復や外部反復の途中で変わることはない。
また、(40)式は三点階差の連立一次方程式であり、中性子束の係数はほとんどが零である。
このように係数のほとんどが零の行列を疎行列と呼ぶ。核計算では非常に多くの元数の連立一次方程式を取り扱うことから、
大規模疎行列と呼ばれている。
例えば、110 万 kWe 級 26の BWR-5の 1/4 炉心を集合体均質の三次元体系で計算する場合、行列 A の元数は 756 体×3 群 27×軸方向
24 ノード×1/4=13,608 と元数が非常に大きい。行列 A は 13,608×13,608 の二次元行列であるが、そのほとんどは零であり、零でない係数の数は 13,608×7 と非常に少ない 28。
このように核計算で取り扱う連立一次方程式には一般的な連立一次方程式には無い特徴がある。そのため、これらの特徴を利用することで、内部反復の高速化を実現している。以降でkWe の”e”は単位ではなく、発電量(electricity)であることを示す。発熱量(thermal)は kWthと”th”を付けて区別している。
27 エネルギー群数を増やすと連立一次方程式の元数は大幅に増加するため、計算時間に大きな影響がある。そのため、いかに少ないエネルギー群数で高い計算精度を得るかという研究が長年続けられている。例えば古い BWR の炉心解析コードでは高速化のため、エネルギー群数を減らした修正一群計算(参考文献 20 及び Appendix A 参照)を用いている場合もある。
28 厳密には境界条件を考慮する必要があるため、13,608×7 よりも小さくなる。例えば一番底𝑧𝑧−
部(1 ノード目、k=1)では、𝐴𝐴𝑔𝑔,𝑖𝑖,𝑗𝑗,𝑘𝑘は含まれないため、零ではない係数の数は 6 個となる。
26

16

はこれらの特徴を利用した連立一次方程式の計算の高速化について説明する。
3.2.1.2 LU 分解
まず、左辺の係数行列が変化しないという特徴を利用した内部反復の高速化について説明する。ガウスの消去法を用いた連立一次方程式の計算で長い計算時間が必要なのは、(57)式のように上三角行列の形に変形する操作である。そのため、この上三角行列の形に変形する操作を一度で済ますことが出来れば、計算時間の大幅な高速化が実現できる。
行列 A を(60)式のようにあらかじめ下三角行列 L と上三角行列 U に分解しておき、右辺のベ
クトル𝑏𝑏�⃗が変わってもあらかじめ作成した L と U で計算することで複数の連立一次方程式の計算時間を高速化する手法を LU 分解と呼ぶ。
𝑨𝑨 = 𝑳𝑳𝑳𝑳(60)

以降では下三角行列 L と上三角行列 U の計算について説明する。(55)式、(56)式の操作を詳細に記述すると、
1⎡− 𝑎𝑎 ⁄𝑎𝑎
⎢ 21 11⎢− 𝑎𝑎31 ⁄𝑎𝑎11⎢− 𝑎𝑎41 ⁄𝑎𝑎11⋮⎢⎣− 𝑎𝑎𝑛𝑛1⁄𝑎𝑎11

0
1
0
0
⋮
0

0
0
1
0
⋮
0

0
0
0
1
⋮
0

⋯
⋯
⋯
⋯
⋱
⋯

0 𝑎𝑎11⎤0 ⎡ 𝑎𝑎21⎥⎢0⎥ ⎢ 𝑎𝑎310⎥ ⎢ 𝑎𝑎41⋮⎥⎢ ⋮1⎦ ⎣𝑎𝑎𝑛𝑛1

𝑎𝑎11 𝑎𝑎12⎡(2)
0 𝑎𝑎22⎢(2)⎢ 0 𝑎𝑎32=⎢(2)⎢ 0 𝑎𝑎42⋮⎢ ⋮(2)⎣ 0 𝑎𝑎𝑛𝑛2

𝑎𝑎13

𝑎𝑎12
𝑎𝑎22
𝑎𝑎32
𝑎𝑎42
⋮
𝑎𝑎𝑛𝑛2

𝑎𝑎13
𝑎𝑎23
𝑎𝑎33
𝑎𝑎43
⋮
𝑎𝑎𝑛𝑛3

𝑎𝑎14
𝑎𝑎24
𝑎𝑎34
𝑎𝑎44
⋮
𝑎𝑎𝑛𝑛4

𝑎𝑎23

(2)

𝑎𝑎14

𝑎𝑎24

⋯

⋮

⋮

⋯
⋱
⋮ ⎥
(2)
⋯ 𝑎𝑎𝑛𝑛𝑛𝑛 ⎦

(2)
𝑎𝑎33
(2)
𝑎𝑎43
(2)

𝑎𝑎𝑛𝑛3

(2)

(2)
𝑎𝑎34
(2)
𝑎𝑎44
(2)

𝑎𝑎𝑛𝑛4

𝑎𝑎1𝑛𝑛

(2) ⎤

⋯ 𝑎𝑎2𝑛𝑛
⋯

⋯
⋯
⋯
⋯
⋱
⋯

𝑎𝑎1𝑛𝑛
𝑎𝑎2𝑛𝑛 ⎤
𝑎𝑎3𝑛𝑛 ⎥
𝑎𝑎4𝑛𝑛 ⎥⎥
⋮ ⎥𝑎𝑎𝑛𝑛𝑛𝑛 ⎦

(61)

⎥

(2)
𝑎𝑎3𝑛𝑛 ⎥
(2) ⎥𝑎𝑎4𝑛𝑛 ⎥となる。同様に、(58)式、(59)式の操作を行列形式で記述すると、
1
⎡⋮
⎢0
⎢
⎢0
⎢0
⎣0

⋯
⋱
⋯
⋯
⋯
⋯

0
⋮
1
0
0
0

0
0 0
⋮
⋮ ⋮
0
0 0
1
0 0
1
0
(𝑛𝑛−1)(𝑛𝑛−1)
0 − 𝑎𝑎𝑛𝑛,𝑛𝑛−1�𝑎𝑎𝑛𝑛−1,𝑛𝑛−1

⎤
⎥
⎥⋯
⎥
0⎥
1⎦

(62)

17

1
⎡− 𝑎𝑎 ⁄𝑎𝑎
⎢ 21 11⎢− 𝑎𝑎31 ⁄𝑎𝑎11⎢− 𝑎𝑎41 ⁄𝑎𝑎11⋮⎢⎣− 𝑎𝑎𝑛𝑛1⁄𝑎𝑎11

0
1
0
0
⋮
0

0
0
1
0
⋮
0

0
0
0
1
⋮
0

𝑎𝑎11
⎡
0
⎢
⎢
=⎢ 0
0
⎢
⎢ ⋮
⎣ 0

⋯
⋯
⋯
⋯
⋱
⋯

𝑎𝑎12

0 𝑎𝑎11⎤0 ⎡ 𝑎𝑎21⎥⎢0⎥ ⎢ 𝑎𝑎310⎥ ⎢ 𝑎𝑎41⋮⎥⎢ ⋮1⎦ ⎣𝑎𝑎𝑛𝑛1

(2)
𝑎𝑎22

0
0

⋮
0

𝑎𝑎13

(2)
𝑎𝑎23
(3)
𝑎𝑎33

0
⋮
0

𝑎𝑎12
𝑎𝑎22
𝑎𝑎32
𝑎𝑎42
⋮
𝑎𝑎𝑛𝑛2

𝑎𝑎14

となる。ここで、
𝑚𝑚𝑖𝑖𝑖𝑖 =

0
0
0
0
1
0
−𝑚𝑚43
1
⋮
⋮
−𝑚𝑚𝑛𝑛3 −𝑚𝑚𝑛𝑛4

𝑎𝑎11
⎡
0
⎢
⎢
=⎢ 0
0
⎢
⎢ ⋮
⎣ 0

𝑎𝑎12

(2)

𝑎𝑎22
0
0
⋮
0

⋯
⋯
⋯
⋯
⋱
⋯

𝑎𝑎13

(2)

𝑎𝑎23

(3)

𝑎𝑎33
0
⋮
0

𝑎𝑎14
𝑎𝑎24
𝑎𝑎34
𝑎𝑎44
⋮
𝑎𝑎𝑛𝑛4

⋯

𝑎𝑎1𝑛𝑛

0 𝑎𝑎11⎤0 ⎡ 𝑎𝑎21⎥⎢0⎥ ⎢ 𝑎𝑎310⎥ ⎢ 𝑎𝑎41⋮⎥⎢ ⋮1⎦ ⎣𝑎𝑎𝑛𝑛1

𝑎𝑎12
𝑎𝑎22
𝑎𝑎32
𝑎𝑎42
⋮
𝑎𝑎𝑛𝑛2

(2)
𝑎𝑎24
(3)
𝑎𝑎34
(4)
𝑎𝑎44

⋮
0

となる。(62)式の左辺を整理すると、(62)式は
1
0
⎡ −𝑚𝑚
1
⎢ −𝑚𝑚 21 −𝑚𝑚
31
32⎢⎢ −𝑚𝑚41 −𝑚𝑚42⋮⎢ ⋮⎣−𝑚𝑚𝑛𝑛1 −𝑚𝑚𝑛𝑛2

𝑎𝑎13
𝑎𝑎23
𝑎𝑎33
𝑎𝑎43
⋮
𝑎𝑎𝑛𝑛3

𝑎𝑎14

(2)

𝑎𝑎24

(3)

𝑎𝑎34

(4)

𝑎𝑎44
⋮
0

⋯
⋯
⋯
⋯
⋱
⋯

𝑎𝑎1𝑛𝑛
𝑎𝑎2𝑛𝑛 ⎤
𝑎𝑎3𝑛𝑛 ⎥
𝑎𝑎4𝑛𝑛 ⎥⎥
⋮ ⎥𝑎𝑎𝑛𝑛𝑛𝑛 ⎦

𝑎𝑎13
𝑎𝑎23
𝑎𝑎33
𝑎𝑎43
⋮
𝑎𝑎𝑛𝑛3

𝑎𝑎14
𝑎𝑎24
𝑎𝑎34
𝑎𝑎44
⋮
𝑎𝑎𝑛𝑛4

(2) ⎤

⋯ 𝑎𝑎2𝑛𝑛⎥(3)⋯ 𝑎𝑎3𝑛𝑛 ⎥(4) ⎥⋯ 𝑎𝑎4𝑛𝑛 ⎥
⋱
⋮ ⎥
(𝑛𝑛)
⋯ 𝑎𝑎𝑛𝑛𝑛𝑛 ⎦

⋯

𝑎𝑎1𝑛𝑛

(2) ⎤

⋯
⋯
⋯
⋯
⋱
⋯

𝑎𝑎1𝑛𝑛
𝑎𝑎2𝑛𝑛 ⎤
𝑎𝑎3𝑛𝑛 ⎥
𝑎𝑎4𝑛𝑛 ⎥⎥
⋮ ⎥𝑎𝑎𝑛𝑛𝑛𝑛 ⎦(63)⋯ 𝑎𝑎2𝑛𝑛⎥(3)⋯ 𝑎𝑎3𝑛𝑛 ⎥(4) ⎥⋯ 𝑎𝑎4𝑛𝑛 ⎥
⋱
⋮ ⎥
(𝑛𝑛)
⋯ 𝑎𝑎𝑛𝑛𝑛𝑛 ⎦(𝑗𝑗)𝑎𝑎𝑖𝑖,𝑗𝑗

(64)

(𝑗𝑗)

𝑎𝑎𝑗𝑗,𝑗𝑗である。(63)式の左辺の左の行列の逆行列は
1
0
⎡ −𝑚𝑚
1
⎢ −𝑚𝑚 21 −𝑚𝑚
31
32
⎢
−𝑚𝑚
−𝑚𝑚
41
42
⎢
⋮
⋮
⎢
⎣−𝑚𝑚𝑛𝑛1 −𝑚𝑚𝑛𝑛2

0
0
1
−𝑚𝑚43
⋮
−𝑚𝑚𝑛𝑛30001⋮−𝑚𝑚𝑛𝑛4

⋯
⋯
⋯
⋯
⋱
⋯

0
⎤
0
⎥
0⎥
0⎥
⋮⎥
1⎦

−𝟏𝟏

1
⎡ 𝑚𝑚
⎢ 𝑚𝑚2131= ⎢ 𝑚𝑚⎢ 41⎢ ⋮⎣𝑚𝑚𝑛𝑛1

0
1
𝑚𝑚32
𝑚𝑚42
⋮
𝑚𝑚𝑛𝑛2

となる。(65)式の行列を(63)式の両辺の左側から掛けると、

0
0
0
0
1
0
𝑚𝑚43 1⋮⋮
𝑚𝑚𝑛𝑛3 𝑚𝑚𝑛𝑛4

⋯
⋯
⋯
⋯
⋱
⋯

0
⎤
0
⎥
0⎥
0⎥
⋮⎥
1⎦

(65)

18

𝑎𝑎11
⎡ 𝑎𝑎21⎢ 𝑎𝑎31
⎢ 𝑎𝑎
⎢ 41
⎢ ⋮
⎣𝑎𝑎𝑛𝑛1

𝑎𝑎12
𝑎𝑎22
𝑎𝑎32
𝑎𝑎42
⋮
𝑎𝑎𝑛𝑛2

1

𝑎𝑎13
𝑎𝑎23
𝑎𝑎33
𝑎𝑎43
⋮
𝑎𝑎𝑛𝑛3

0
1
𝑚𝑚32
𝑚𝑚42
⋮
𝑚𝑚𝑛𝑛2

⎡𝑚𝑚
⎢ 𝑚𝑚2131= ⎢ 𝑚𝑚⎢ 41⎢ ⋮⎣𝑚𝑚𝑛𝑛1

𝑎𝑎14
𝑎𝑎24
𝑎𝑎34
𝑎𝑎44
⋮
𝑎𝑎𝑛𝑛4

⋯
⋯
⋯
⋯
⋱
⋯

0
0
0
0
1
0
𝑚𝑚43 1⋮⋮
𝑚𝑚𝑛𝑛3 𝑚𝑚𝑛𝑛4

𝑎𝑎1𝑛𝑛
𝑎𝑎2𝑛𝑛 ⎤
𝑎𝑎3𝑛𝑛 ⎥
𝑎𝑎4𝑛𝑛 ⎥⎥
⋮ ⎥𝑎𝑎𝑛𝑛𝑛𝑛 ⎦
⋯
⋯
⋯
⋯
⋱
⋯

𝑎𝑎11
0 ⎡
⎤ 0
0 ⎢
⎥
0⎥ ⎢ 00⎥ ⎢ 0
⋮⎥⎢
⋮
1⎦ ⎢
⎣ 0

𝑎𝑎12

(2)
𝑎𝑎22

0
0
⋮
0

𝑎𝑎13

(2)
𝑎𝑎23
(3)
𝑎𝑎33

0
⋮
0

𝑎𝑎14

(2)
𝑎𝑎24
(3)
𝑎𝑎34
(4)
𝑎𝑎44

⋮
0

⋯

𝑎𝑎1𝑛𝑛

(2) ⎤

⋯ 𝑎𝑎2𝑛𝑛

(66)

⎥
(3)
⋯ 𝑎𝑎3𝑛𝑛 ⎥(4) ⎥⋯ 𝑎𝑎4𝑛𝑛 ⎥
⋱
⋮ ⎥
(𝑛𝑛)
⋯ 𝑎𝑎𝑛𝑛𝑛𝑛 ⎦となる。(66)式は(60)式と同じであり、以上の操作により、行列 A を下三角行列 L と上三角行
列 U に分解することが出来る。あとはガウスの消去法と同様である。(50)式は LU 分解を行うと、
𝑨𝑨𝑥𝑥⃗ = 𝑳𝑳𝑳𝑳𝑥𝑥⃗ = 𝑏𝑏�⃗

(67)となる。ここで、���⃗
𝑼𝑼𝑥𝑥⃗ = 𝑥𝑥′(68)とおくと、(67)式は���⃗ = 𝑏𝑏�⃗
𝑳𝑳𝑥𝑥′

(69)

となる。(69)式を詳細に記述すると、
1
⎡𝑚𝑚
21
⎢ 𝑚𝑚
⎢ 31
⎢ 𝑚𝑚41⎢ ⋮⎣ 𝑚𝑚𝑛𝑛1

0
1
𝑚𝑚32
𝑚𝑚42
⋮
𝑚𝑚𝑛𝑛2

0
0
0
0
1
0
𝑚𝑚43 1⋮⋮
𝑚𝑚𝑛𝑛3 𝑚𝑚𝑛𝑛4

⋯
⋯
⋯
⋯
⋱
⋯

𝑏𝑏1
0 𝑥𝑥′1⎤ 𝑥𝑥′𝑏𝑏
0 ⎛ 2⎞⎛ 2⎞⎥0⎥ ⎜𝑥𝑥′3 ⎟ ⎜𝑏𝑏3 ⎟=0⎥ ⎜𝑥𝑥′4 ⎟ ⎜𝑏𝑏4 ⎟⋮⋮⎥ ⋮1⎦ ⎝𝑥𝑥′𝑛𝑛 ⎠ ⎝𝑏𝑏𝑛𝑛 ⎠(70)

���⃗を計算することなり、上三角行列と下三角行列の違いはあるものの、(57)式と同様に簡単に𝑥𝑥′���⃗を用いると、(68)式の関係からとが出来る。得られた𝑥𝑥′
𝑎𝑎11 𝑎𝑎12 𝑎𝑎13 𝑎𝑎14 ⋯ 𝑎𝑎1𝑛𝑛
𝑥𝑥′1
⎡
(2)
(2)
(2)
(2) ⎤ 𝑥𝑥1
0 𝑎𝑎22 𝑎𝑎23 𝑎𝑎24 ⋯ 𝑎𝑎2𝑛𝑛
⎢
⎥ 𝑥𝑥2
𝑥𝑥′
⎛ 2⎞
(3)
(3)
(3)
⎢ 0
𝑎𝑎33 𝑎𝑎34 ⋯ 𝑎𝑎3𝑛𝑛 ⎥ ⎛𝑥𝑥3 ⎞ ⎜ 𝑥𝑥′3 ⎟
0
(71)
⎢ 0
(4) ⎥ ⎜
(4)
⎜𝑥𝑥4 ⎟⎟ = ⎜ 𝑥𝑥′4 ⎟
0
⋯
𝑎𝑎
0
𝑎𝑎
4𝑛𝑛
44
⎢
⎥ ⋮
⋮
⋱
⋮ ⎥ ⎝𝑥𝑥 ⎠
⋮
⋮
⋮
⎢ ⋮
𝑛𝑛
(𝑛𝑛)
𝑥𝑥′
⎝
𝑛𝑛 ⎠
⋯ 𝑎𝑎𝑛𝑛𝑛𝑛 ⎦
0
0
0
⎣ 0
となり、𝑥𝑥⃗が計算できる。

行列 A を下三角行列 L と上三角行列 U にさえ分解してしまえば、(70)式、(71)式を用いた𝑥𝑥⃗の計算はほとんど計算時間がかからない。そのため LU 分解は、拡散方程式のような行列 A が変化せず、右辺のベクトル𝑏𝑏�⃗のみが変化する複数の連立一次方程式を解く場合には、非常に効果的な手法である。核計算では、(35)式の連立一次方程式を数十から数百回も解く必要がある 29。

そのため、LU 分解を用いることで、ガウスの消去法に比べて数十倍から数百倍の高速化が期
29

もちろん、後述する外部反復の加速法の利用によっても外部反復の数は変わってくるが、一
般的には数十～数百回である。
19

待できる。
3.2.1.3 ヤコビ法
次に、大規模疎行列であるという特徴を利用した内部反復の高速化について説明する。本手法と以降で説明するガウス・ザイデル法及び Successive Over Relaxation 法(SOR 法、逐次過緩和法)は反復法と呼ばれ、連立一次方程式を反復計算で解く手法である。なお、ガウスの消去法やLU 分解のように反復計算ではなく、四則演算を用いて直接計算する手法のことを直接法と呼ぶ。
反復法では、次式のように連立一次方程式の係数行列 A を対角部分 D、下三角部分 L、上三角部分 U 30に分ける。
𝑨𝑨 = 𝑫𝑫 + 𝑳𝑳 + 𝑼𝑼

(72)

ここで、

𝑎𝑎11
⎡ 0
⎢
0
𝑫𝑫 = ⎢
⎢ 0
⎢ ⋮
⎣ 0

0
𝑎𝑎22
0
0
⋮
0

0
0
⎡ 𝑎𝑎
0
⎢ 𝑎𝑎 21 𝑎𝑎
31
32
𝑳𝑳 = ⎢ 𝑎𝑎⎢ 41 𝑎𝑎42⋮⎢ ⋮⎣𝑎𝑎𝑛𝑛1 𝑎𝑎𝑛𝑛2

0
0
𝑎𝑎33
0
⋮
0

0
0
0
𝑎𝑎44
⋮
0

0 0
0 0
0
0
𝑎𝑎43 0⋮⋮𝑎𝑎𝑛𝑛3 𝑎𝑎𝑛𝑛4

⋯ 0
⎤
⋯ 0
⎥
⋯ 0 ⎥
⋯ 0 ⎥
⋱
⋮ ⎥
⋯ 𝑎𝑎𝑛𝑛𝑛𝑛 ⎦
⋯
⋯
⋯
⋯
⋱
⋯

(73)

0
⎤
0
⎥
0⎥
0⎥
⋮⎥
0⎦

(74)

0 𝑎𝑎12 𝑎𝑎13 𝑎𝑎14 ⋯ 𝑎𝑎1𝑛𝑛⎡0
0 𝑎𝑎23 𝑎𝑎24 ⋯ 𝑎𝑎2𝑛𝑛 ⎤⎢0
0 𝑎𝑎34 ⋯ 𝑎𝑎3𝑛𝑛 ⎥⎥0
𝑼𝑼 = ⎢⋯ 𝑎𝑎4𝑛𝑛 ⎥
0
0
0
⎢0
⎢ ⋮
⋮
⋮
⋮
⋱
⋮ ⎥
⎣0
0
0
0 ⋯ 0 ⎦である。(72)式を(50)式に代入すると、(75)(𝑫𝑫 + 𝑳𝑳 + 𝑼𝑼)𝑥𝑥⃗ = 𝑏𝑏�⃗

(76)となる。(76)式を次式のように下三角部分 L と上三角部分 U に分割したような形に書き換える。
𝑥𝑥⃗ = −𝑫𝑫−𝟏𝟏 (𝑳𝑳 + 𝑼𝑼)𝑥𝑥⃗ + 𝑫𝑫−𝟏𝟏 𝑏𝑏�⃗

(77)

ここで、

⎡
⎢
𝑫𝑫−𝟏𝟏 = ⎢⎢⎢⎣1⁄𝑎𝑎1100
0⋮001⁄𝑎𝑎2200⋮
0001⁄𝑎𝑎330⋮0

0
001⁄𝑎𝑎44
⋮
0

⋯
0
⎤
⋯
0
⎥
⋯
0
⎥
⋯
0
⎥
⋱
⋮
⎥
⋯ 1⁄𝑎𝑎𝑛𝑛𝑛𝑛 ⎦(78)

ここでの下三角部分 L と上三角部分 U は LU 分解の上三角行列、下三角行列と異なり、対角成分が 0 であることに注意すること。
30

20である。(77)式が反復法の基礎式となる。反復法は、𝑥𝑥⃗に適当な初期値を与え、𝑥𝑥⃗が設定した誤差未満になるまで(77)式の計算を繰り返す手法である。ヤコビ法と後述するガウス・ザイデル
法及び SOR 法では、(77)式の解き方が異なっているだけで、(77)式を用いるという点では同じである。
ヤコビ法は(77)式を解く上で最もシンプルな手法である。反復回数 l 回目の(77)式は次式となり、𝑥𝑥⃗ (𝑙𝑙) が(80)式や(81)式のような収束条件を満たすまで反復計算を繰り返す。
𝑥𝑥⃗ (𝑙𝑙) = −𝑫𝑫−𝟏𝟏 (𝑳𝑳 + 𝑼𝑼)𝑥𝑥⃗ (𝑙𝑙−1) + 𝑫𝑫−𝟏𝟏 𝑏𝑏�⃗(𝑙𝑙−1) 2(𝑙𝑙)� �𝑥𝑥𝑘𝑘 − 𝑥𝑥𝑘𝑘

� < 𝜀𝜀1

�

< 𝜀𝜀2

𝑘𝑘

(𝑙𝑙)

(𝑙𝑙−1)𝑥𝑥 − 𝑥𝑥� 𝑘𝑘 (𝑙𝑙)𝑘𝑘
𝑥𝑥⃗𝑘𝑘

𝑚𝑚𝑚𝑚𝑚𝑚

(79)
(80)

(81)

ヤコビ法の計算手順は簡単で、次の通りとなる。
1. (74)、(75)式および(78)式を用いて下三角部分 L と上三角部分 U 及び対角逆行列 D-1 を計算する。
2. 𝑥𝑥⃗に初期値𝑥𝑥⃗ (0) を与える。

3. (79)式を用いて𝑥𝑥⃗ (𝑙𝑙) を計算する。

4. (80)式や(81)式の条件を満たすまで 3.の計算を繰り返す。
このようにヤコビ法は(79)式の計算を繰り返すだけで𝑥𝑥⃗が計算できるという非常に単純な手法である。なお、ヤコビ法や後述するガウス・ザイデル法及び SOR 法では、(82)式に示すように
対角成分𝑎𝑎𝑖𝑖𝑖𝑖 が他の成分よりも十分に大きい対角優位な行列でなければ𝑥𝑥⃗ (𝑙𝑙) は発散し、正しい結果が得られない場合があることに注意する必要がある。
|𝑎𝑎𝑖𝑖𝑖𝑖 | > ��𝑎𝑎𝑗𝑗𝑗𝑗 �

(82)

𝑗𝑗≠𝑖𝑖

なお、拡散計算では(35)式のように𝑎𝑎𝑖𝑖𝑖𝑖 が巨視的全断面積を含む項となっており、多くの場合で(82)式の条件を満たす対角優位な行列であることから、計算が発散する恐れはほとんどない。

さて、ここで(79)式のヤコビ法を用いた連立一次方程式を早く収束させるためにはどうすればいいだろうか？ヤコビ法の反復回数は初期値𝑥𝑥⃗ (0) に依存する。例えば、初期値𝑥𝑥⃗ (0) が全て 1.0の場合と、その連立一次方程式の解を初期値として与えた場合とでは当然ながら反復回数に差が出るのは容易に想像が付くと思う。このように、ヤコビ法の計算時間を高速化させるためには、初期値𝑥𝑥⃗ (0) の設定が重要となる。

拡散計算の内部反復の解法としてヤコビ法を用いる場合には、m 回目の外部反復時の連立一次方程式を解く際の初期値に、m-1 回目の外部反復で得られた𝑥𝑥⃗を用いることが多い 31。外部反復の計算では中性子源項である右辺𝑏𝑏�⃗が変化するが、𝑏𝑏�⃗は外部反復によって大きく変化する訳で

外部反復 1 回目については全て 1.0 とする方法や、核分裂断面積で重みを付けた値を使う場合が多い。当然ながら外部反復 1 回目の初期値は全体の計算時間に大きな影響を与えることか
ら、外部反復 1 回目の初期値の設定は高速化の観点から重要である。

31

21

はない。そのため、𝑏𝑏�⃗がやや変化した時の連立一次方程式の解𝑥𝑥⃗も前回の外部反復の𝑥𝑥⃗に近い値となることが期待できる。このように初期値𝑥𝑥⃗ (0) を前回の外部反復の𝑥𝑥⃗とすることで、初期値を全て 1.0 などとする場合に比べて反復回数の低減が期待できる。
3.2.1.4 ガウス・ザイデル法
ガウス・ザイデル法も連立一次方程式の解に近い値を使うことで反復回数の低減を図る手法である。前述した反復回数の低減では、初期値𝑥𝑥⃗ (0) を連立一次方程式の解に近い値に設定することで反復回数の低減を図ったが、ガウス・ザイデル法では計算途中の𝑥𝑥⃗ (𝑙𝑙) をより連立一次方
程式の解に近づけることで反復回数の低減を図っている。

ヤコビ法を用いた連立一次方程式の計算が収束する場合、𝑥𝑥⃗ (𝑙𝑙) と𝑥𝑥⃗ (𝑙𝑙−1) を比べると𝑥𝑥⃗ (𝑙𝑙) の方が連立一次方程式の解に近くなる。そのため、ガウス・ザイデル法では、(79)式の右辺の𝑥𝑥⃗ (𝑙𝑙−1) のうち、すでに𝑥𝑥⃗ (𝑙𝑙) が得られているものは𝑥𝑥⃗ (𝑙𝑙−1) の代わりに𝑥𝑥⃗ (𝑙𝑙) を使うことで反復回数の低減を図っている。

ここで、反復回数 l 回目の𝑥𝑥⃗を計算することを考える。(76)式を用いて𝑥𝑥1 、𝑥𝑥2 、𝑥𝑥3 …と順番に計算していく場合、𝑥𝑥𝑛𝑛 を計算する際には𝑥𝑥1 ～𝑥𝑥𝑛𝑛−1 はすでに(76)式を用いて更新されている。(73)～(76)式を見ると、𝑥𝑥1 、𝑥𝑥2 、𝑥𝑥3 …の係数はそれぞれ𝑎𝑎𝑖𝑖1、𝑎𝑎𝑖𝑖2、𝑎𝑎𝑖𝑖3…となる。(76)式を用いて𝑥𝑥3 を計算する場合、下三角部分 L の要素である 𝑎𝑎31、𝑎𝑎32 の係数𝑥𝑥1 、𝑥𝑥2 にはすでに更新された値が存在する。このように、(76)式を用いて計算する場合、下三角部分 L に掛かる𝑥𝑥⃗は全て更新された
𝑥𝑥⃗ (𝑙𝑙) を用いることが出来る。

そこでガウス・ザイデル法では更新された𝑥𝑥⃗ (𝑙𝑙) を用いるように、(76)式を変形する。反復回数

l 回目のガウス・ザイデル法の計算式は(83)式の通りとなる。
𝑥𝑥⃗ (𝑙𝑙) = −𝑫𝑫−𝟏𝟏 𝑳𝑳𝑥𝑥⃗ (𝑙𝑙) − 𝑫𝑫−𝟏𝟏 𝑼𝑼𝑥𝑥⃗ (𝑙𝑙−1) + 𝑫𝑫−𝟏𝟏 𝑏𝑏�⃗(83)

なお、(83)式は上から(𝑎𝑎11 から)順番に計算していく場合の式である。このようにガウス・ザイ

デル法では新しい𝑥𝑥⃗ (𝑙𝑙) を用いることで、ヤコビ法に比べてより早く𝑥𝑥⃗が収束することが期待できる。

3.2.1.5 Successive Over Relaxation (SOR)法
ガウス・ザイデル法では、より連立一次方程式の解に近い𝑥𝑥⃗ (𝑙𝑙) を用いることで、反復回数の低減を図っていた。反復法を用いた連立一次方程式の計算では、𝑥𝑥⃗ (𝑙𝑙−1) から𝑥𝑥⃗ (𝑙𝑙) となって急に解に近づくのではなく、徐々に解に近づいていく。ガウス・ザイデル法では、𝑥𝑥⃗ (𝑙𝑙−1) の代わりにより解に近い𝑥𝑥⃗ (𝑙𝑙) を用いることで反復回数を低減させていたが、SOR 法ではさらに緩和係数 ωと呼ばれるパラメータを用いて𝑥𝑥⃗ (𝑙𝑙) を人為的に調整し、より解に近づけることで反復回数の低減を図っている。

𝑥𝑥⃗ (𝑙𝑙−1) から𝑥𝑥⃗ (𝑙𝑙) となる際に、どの程度値が変化するかを考える。(83)式のガウス・ザイデル法

を次式のように書き換える。

𝑥𝑥⃗ (𝑙𝑙) − 𝑥𝑥⃗ (𝑙𝑙−1) = −𝑫𝑫−𝟏𝟏 𝑳𝑳𝑥𝑥⃗ (𝑙𝑙) − (𝑫𝑫−𝟏𝟏 𝑼𝑼 + 𝑬𝑬)𝑥𝑥⃗ (𝑙𝑙−1) + 𝑫𝑫−𝟏𝟏 𝑏𝑏�⃗(84)

ここで、𝑬𝑬は単位行列である。(84)式は l 回目の反復計算において、(84)式の右辺の分だけ変化することを意味している。SOR 法では(85)式のように(84)式の右辺に緩和係数 ω を掛けること
22

で、𝑥𝑥⃗ (𝑙𝑙−1) から𝑥𝑥⃗ (𝑙𝑙) への変化量を人為的に調整し、反復回数の低減を図っている。
(𝑙𝑙)

���⃗′
𝑥𝑥

= 𝑥𝑥⃗ (𝑙𝑙−1) + 𝜔𝜔�−𝑫𝑫−𝟏𝟏 𝑳𝑳𝑥𝑥⃗ (𝑙𝑙) − (𝑫𝑫−𝟏𝟏 𝑼𝑼 + 𝑬𝑬)𝑥𝑥⃗ (𝑙𝑙−1) + 𝑫𝑫−𝟏𝟏 𝑏𝑏�⃗�(85)

ここで、ω が 1.0 の時はガウス・ザイデル法と同じである。
ω が 1.0 よりも大きい場合には、
𝑥𝑥⃗ (𝑙𝑙−1)から𝑥𝑥⃗ (𝑙𝑙) への変化量が大きくなり、より早く収束することが期待できる。なお、ω が 1.0 よりも大きい場合には前進加速と言い、ω が 1.0 よりも小さい場合には後退加速と呼ぶ。前進加速は解に対して一次関数状に近づいていく場合に、後退加速は解に対して振動している場合に有効な手法である。(85)式では𝑥𝑥⃗ (𝑙𝑙) が一次関数状に変化すると仮定しているが、𝑥𝑥⃗ (𝑙𝑙) が二次関数状に変化する場合には𝑥𝑥⃗ (𝑙𝑙−1) だけでなく𝑥𝑥⃗ (𝑙𝑙−2) を用いて補正することで、効果的に収束させることができる。

SOR 法では、緩和係数 ω を
1.0 < 𝜔𝜔 < 2.0(86)とすることが多い。これは ω の値が大きい場合、解に収束させるつもりが逆に解から離れてしまう場合があるためである。このように SOR 法では ω の値によっては、収束が遅くなる場合
や、計算が発散してしまう場合があるので、ω の値の設定には十分に注意する必要がある。一般的には、反復回数が少ない時には連立一次方程式の解から離れており、反復回数が増えるにつれて解に近くなる。(85)式では ω は反復回数によって一定であるが、(80)式や(81)式の左辺の値を指標として、解から遠い場合(すなわち(80)式や(81)式の左辺の値が大きい場合)には ω を大きな値とし、解に近づくにつれて ω も 1.0 に近い値にすることで効率的に収束させることが期待できる 32。
SOR 法は実装が簡単なこともあり、連立一次方程式の解法として広く利用されている。連立一次方程式の解法としては、他にも交互方向陰解法(Alternating-Direction Implicit Iterative method、ADI 法)21)、共役勾配法(Conjugate Gradient method、CG 法)22)や Krylov 部分空間法 23,24)を用いた手法などがあり、解析したい連立一次方程式の特徴によって手法を選ぶことが重要である。
3.2.1.6 反復法を内部反復に用いるときの収束条件
一般的に、外部反復も内部反復も収束条件が厳しくなればなるほど、反復計算回数は指数関数状に増加する。そのため、反復法を用いた連立一次方程式の解法では、厳しすぎる収束条件を適用すると非常に長い計算時間が必要となる。
ここで、内部反復の収束条件を厳しくした場合(例えば 1.0×10-3)と、緩くした場合(例えば1.0×10-2)での拡散計算の計算時間を比較する。外部反復の回数だけで考えると、当然収束条件を厳しくした方が、すなわちそれぞれの内部反復で中性子を完全に収束させた方が外部反復の回数は少なくなる。しかし、計算時間で考えると、内部反復の収束条件を緩くした方が早くなることがある。
但し、緩和係数 ω を細かく調整することで余計に計算時間が掛かってしまう場合もあり、解析する条件によって ω を一定とした方がいいのか、反復回数によって変えた方がいいのかは変化する。SOR 法を用いる場合には、最初は ω を固定としておき、収束の仕方などを見ながら ω の値の設定や ω を収束の途中で変えるなどの対策を考えることが望ましい。

32

23

外部反復の回数が少ない時点では、実際の中性子束分布とは大きく異なる分布になっている。
そのような状況では、中性子の収束条件が 1.0×10-3 でも 1.0×10-2 でも実際の中性子束からの差異という意味では大きな違いはない。しかし、内部反復の計算時間は大きく違うことから、外部反復回数が少ない時点では中性子束を完全に収束させることに時間をかけるよりも、中性子源項を更新して次の外部反復に移ってしまった方が、外部反復回数はやや増えるものの、全体の計算時間が短くなることが多い。
そのため、反復法を用いた拡散計算では、内部反復の上限回数を 10～20 回程度と少なくしたり 33、外部反復(実効増倍率)の収束条件に比べて内部反復の収束条件を緩くしたり、外部反復の回数によって内部反復の収束条件を変えたりすることで計算時間の高速化を図っている。
3.2.2 その他の高速化手法
3.2.2.1 除去断面積 Σr の導入(32)式に示したように、教科書などでは拡散方程式の左辺は拡散係数𝐷𝐷𝑔𝑔,𝑖𝑖 の二次微分項と

Σ𝑡𝑡,𝑔𝑔,𝑖𝑖 𝜙𝜙𝑔𝑔,𝑖𝑖 の二つの項で表記される。しかし、一般的な炉心解析コードでは全断面積Σ𝑡𝑡,𝑔𝑔,𝑖𝑖 を用いずに除去断面積(removal cross section)を用いることが多い。

(32)式の拡散方程式は右辺の中性子源の項にΣ𝑠𝑠,𝑔𝑔→𝑔𝑔,𝑖𝑖 𝜙𝜙𝑔𝑔,𝑖𝑖 が含まれていることから、両辺に g群の中性子束𝜙𝜙𝑔𝑔,𝑖𝑖 が含まれる形になっている。右辺は外部反復で得られることから、右辺の𝜙𝜙𝑔𝑔,𝑖𝑖

は一つ前の外部反復で得られた𝜙𝜙𝑔𝑔,𝑖𝑖 を用いることになり、収束性が悪化する原因となる。ガウス・ザイデル法の説明で述べたことと同様に、外部反復でも新しい中性子束𝜙𝜙𝑔𝑔,𝑖𝑖 を用いること

で収束性がよくなる。
そこで一般的な拡散計算コードでは(32)式の右辺のΣ𝑠𝑠,𝑔𝑔→𝑔𝑔,𝑖𝑖 𝜙𝜙𝑔𝑔,𝑖𝑖 を左辺に移し、(87)式のように変形することで収束性の向上を図っている。
−𝐷𝐷𝑔𝑔,𝑖𝑖𝜒𝜒𝑔𝑔,𝑖𝑖d2
𝜙𝜙 + Σ𝑟𝑟,𝑔𝑔,𝑖𝑖 𝜙𝜙𝑔𝑔,𝑖𝑖 =� 𝜈𝜈Σ ′ 𝜙𝜙 ′ + � Σ𝑠𝑠,𝑔𝑔′ →𝑔𝑔,𝑖𝑖 𝜙𝜙𝑔𝑔′ ,𝑖𝑖
d𝑥𝑥 2 𝑔𝑔,𝑖𝑖
𝑘𝑘 ′ 𝑓𝑓,𝑔𝑔 ,𝑖𝑖 𝑔𝑔 ,𝑖𝑖′

ここで、Σ𝑟𝑟,𝑔𝑔,𝑖𝑖 は除去断面積で、

𝑔𝑔

Σ𝑟𝑟,𝑔𝑔,𝑖𝑖 = Σ𝑡𝑡,𝑔𝑔,𝑖𝑖 − Σ𝑠𝑠,𝑔𝑔→𝑔𝑔,𝑖𝑖(87)

𝑔𝑔 ≠𝑔𝑔(88)である。

𝜒𝜒𝑔𝑔,𝑖𝑖

散乱源と同様に(87)式において核分裂源項の 𝑘𝑘 𝜈𝜈Σ𝑓𝑓,𝑔𝑔,𝑖𝑖 𝜙𝜙𝑔𝑔,𝑖𝑖 も左辺に移せば更なる高速化に繋がるのではないかと考える人がいるかもしれない。しかし、核分裂源項には実効増倍率 k が掛かっており、散乱源と同じようにはいかない。
実効増倍率 k は一つ前の外部反復の中性子束を用いて
𝑘𝑘 = � � 𝜈𝜈Σ𝑓𝑓,𝑔𝑔′ ,𝑖𝑖 𝜙𝜙𝑔𝑔′ ,𝑖𝑖 𝑉𝑉𝑖𝑖
𝑖𝑖

𝑔𝑔′

(89)

で計算される。すなわち(89)式において、実効増倍率 k は核分裂で発生する中性子の総和が体系全体で 1.0 となるように規格化している。そのため、核分裂中性子源項の中性子束を現在の外部反復の中性子束に変えると、実効増倍率を求めた際の中性子束と異なり、(89)式が成立し外部反復の回数が少ない時はこの程度の内部反復回数では中性子束が十分に収束することはない。しかし、十分に実効増倍率が収束してくると、外部反復による中性子束の変化が少な
くなり、この程度の内部反復回数でも中性子束が収束するようになる。

33

24

ない。このように核分裂中性子源項の中性子束を現在の外部反復の中性子束に変えると、中性子バランスが崩れてしまい、収束性が悪化する可能性があるため、注意が必要である。
3.2.2.2 散乱源の更新(35)式を見ると分かる通り、(35)式の左辺は 1 つのエネルギー群数のみの式となっている。そのため、(35)式の連立一次方程式を計算する場合、全エネルギー群をまとめて一つの連立一次方程式として解くのではなく、エネルギー群ごとに独立な連立一次方程式として計算している。
ここで中性子の散乱反応を考える。中性子の散乱は主に下方散乱と上方散乱に分けられる。
散乱前の中性子のエネルギーよりも散乱後の中性子のエネルギーの方が低くなる散乱を下方散乱と呼び、散乱後の中性子のエネルギーの方が高くなる散乱を上方散乱と呼ぶ 34。上方散乱
3は、原子核の熱運動によるエネルギー(最確値：k B 𝑇𝑇、平均値：2 k B 𝑇𝑇)をもらって中性子のエネ

ルギーが高くなる現象である。そのため、原子核の熱運動エネルギー(293.6K でk B 𝑇𝑇 = 0.0253eV、
3
2

k B 𝑇𝑇 = 0.0380eV)以下の低いエネルギー領域でのみ発生する 35。

3.2.2.1 節で述べたように、散乱源項の中性子束は新しい中性子束を用いることで収束性がよくなる。ここで、n 回目の外部反復の時の内部反復での中性子束𝜙𝜙𝑔𝑔,𝑖𝑖 を計算することを考える。

散乱源の更新を考慮すると、g 群の計算式は(90)式のように書き表せる。
−𝐷𝐷𝑔𝑔,𝑖𝑖𝜒𝜒𝑔𝑔,𝑖𝑖d2
𝜙𝜙𝑔𝑔,𝑖𝑖 + Σ𝑟𝑟,𝑔𝑔,𝑖𝑖 𝜙𝜙𝑔𝑔,𝑖𝑖 =� 𝜈𝜈Σ ′ 𝜙𝜙 𝑛𝑛−1
+ � Σ𝑠𝑠,𝑔𝑔′ →𝑔𝑔,𝑖𝑖 𝜙𝜙𝑔𝑔𝑛𝑛′ ,𝑖𝑖 + � Σ𝑠𝑠,𝑔𝑔′ →𝑔𝑔,𝑖𝑖 𝜙𝜙𝑔𝑔𝑛𝑛−1
′
′ ,𝑖𝑖
2
d𝑥𝑥
𝑘𝑘 ′ 𝑓𝑓,𝑔𝑔 ,𝑖𝑖 𝑔𝑔 ,𝑖𝑖′′
𝑔𝑔

𝑔𝑔 <𝑔𝑔(90)

𝑔𝑔 >𝑔𝑔

ここで、右辺第二項は下方散乱による散乱源項、右辺第三項は上方散乱による散乱源項である。
(90)式のように各エネルギー群の連立一次方程式を計算する際には、高エネルギーの群 36から中性子束の計算を行い、下方散乱による散乱源を n-1 回目の外部反復での中性子束𝜙𝜙𝑔𝑔𝑛𝑛−1′ ,𝑖𝑖 の代わ

𝑛𝑛
りに新しい中性子束𝜙𝜙𝑔𝑔,𝑖𝑖を用いる、すなわち自群よりもエネルギーの高い中性子については新

しい中性子束を用いることで、収束性が向上する。

下方散乱による散乱源だけでなく、上方散乱による散乱源も新しい中性子束を用いる方が更なる高速化に繋がるのではないかと考える人がいるかもしれない。実際、下方散乱による散乱源だけでなく、上方散乱による散乱源も新しい中性子束を用いる方が、外部反復回数は少なくなる。しかし、上方散乱による散乱源を計算するためには、内部反復の中で上方散乱による中性子源を更新するための反復計算が追加で必要となり、内部反復の計算時間が増加するため、あまり有効ではない場合が多い 37。下方散乱による散乱源のみを新しい中性子束とする場合には、高いエネルギー群から低いエネルギー群に向けて順番に各群の連立一次方程式を解くだけでよいため、連立一次方程式の計算回数はエネルギー群数と等しくなる。しかし、上方散乱に
34

散乱しても別の群に移らず、同じ群のままの場合の散乱を自群散乱と呼ぶ。
原子核の速度分布はマクスウェル・ボルツマン分布に従っており、その運動エネルギーの最3
3
確値はk B 𝑇𝑇、平均値は k B 𝑇𝑇である。そのため、厳密に言えば k B 𝑇𝑇以上のエネルギーでも上方散22
乱は発生するが、最大でもk B 𝑇𝑇の数倍程度である。
36 一般的には第 1 群が最もエネルギーが高く、群番号が大きくなるにつれてエネルギーは低くなっていく。
37 中性子スペクトルが非常に軟らかく、
上方散乱の割合が大きい場合には高速化に繋がる場合があるかもしれないが、多くの場合ではむしろ計算時間が増加してしまう。
35

25

よる散乱源も更新しようとすると、注目している群よりも低いエネルギー群の中性子束は決まっていないため、低いエネルギー群の連立一次方程式を計算した後に、再度高いエネルギー群の連立一次方程式を計算する必要がある。このように上方散乱による散乱源を更新するためには、反復計算が必要となる。そのため、連立一次方程式の計算回数はエネルギー群数よりも多くなる。この連立一次方程式の計算回数の増加によって、一回の外部反復の計算時間が長くな
り、結果的に全体の計算時間が増加してしまう。
このように外部反復回数の削減に繋がる手法でも、内部反復の計算時間次第では全体の計算時間が増加してしまう場合がある。そのため、高速化手法を適用する場合には、適用しようとしている手法が本当に全体の計算時間の高速化に繋がるのかどうかをよく考えることが重要である。
3.3 外部反復に用いられる高速化手法
3.3.1 SOR 法
3.2.1.5 章では SOR 法を連立一次方程式の解法として紹介したが、SOR 法は外部反復にも適𝑛𝑛
は n-1 回目の外部反復の中性子束用可能である。その場合、n 回目の外部反復の中性子束𝜙𝜙𝑔𝑔,𝑖𝑖
𝑛𝑛−1を用いて次のように補正される 13)。
𝜙𝜙𝑔𝑔,𝑖𝑖
𝑛𝑛

𝑛𝑛
𝑛𝑛
𝑛𝑛−1
𝜙𝜙 ′𝑔𝑔,𝑖𝑖 = 𝜙𝜙𝑔𝑔,𝑖𝑖+ 𝜔𝜔�𝜙𝜙𝑔𝑔,𝑖𝑖− 𝜙𝜙𝑔𝑔,𝑖𝑖
�

(91)

𝑛𝑛
外部反復の SOR 法では、(91)式で得られた中性子束𝜙𝜙 ′𝑔𝑔,𝑖𝑖 を n+1 回目の外部反復の実効増倍率

や中性子源項の計算に用いている。3.2.1.5 章で説明したのと同様に、外部反復での SOR 法も一つ前の中性子束と現在の中性子束を用いてより解に近い値に人為的に調整することで反復回数の低減を図っている。
3.3.2 Chebyshev 外挿法
𝑛𝑛𝑛𝑛−1
SOR 法では真の解が n 回目の外部反復の中性子束𝜙𝜙𝑔𝑔,𝑖𝑖と n-1 回目の外部反復の中性子束𝜙𝜙𝑔𝑔,𝑖𝑖を用いて得られる一次関数の先にあると仮定している。例えば真の解が二次関数の先にあるとすれば、中性子束を(91)式を用いて補正するよりも、二次関数を仮定した補正式を用いた方が効率的である。
Chebyshev 外挿法では、一次関数の代わりに、Chebyshev 多項式を用いて中性子束を補正する手法である 25)。Chebyshev 多項式は𝑇𝑇𝑛𝑛 (𝑥𝑥) = cos(𝑛𝑛𝑛𝑛) , 𝑥𝑥 = cos(𝜃𝜃)

(92)

もしくは

𝑇𝑇𝑛𝑛 (𝑥𝑥) =1(−1)𝑛𝑛𝑑𝑑 𝑛𝑛
�1 − 𝑥𝑥 2 𝑛𝑛 (1 − 𝑥𝑥 2 )𝑛𝑛−2(2𝑛𝑛 − 1)‼
𝑑𝑑𝑥𝑥

(93)

で定義され、その定義域は[-1, 1]である 19)。4 次までの Chebyshev 多項式は𝑇𝑇0 (𝑥𝑥) = 1
𝑇𝑇1 (𝑥𝑥) = 𝑥𝑥

𝑇𝑇2 (𝑥𝑥) = 2𝑥𝑥 2 − 1

𝑇𝑇3 (𝑥𝑥) = 4𝑥𝑥 3 − 3𝑥𝑥

(94)
(95)
(96)
(97)
26

𝑇𝑇4 (𝑥𝑥) = 8𝑥𝑥 4 − 8𝑥𝑥 2 + 1

(98)

となる。

Chebyshev 外挿法では n 回目、n-1 回目、n-2 回目の中性子束を用いて次式で表される。
𝑛𝑛

𝑛𝑛−2

𝑛𝑛−1

𝑛𝑛−2

𝑛𝑛−1

𝑛𝑛
𝜙𝜙 ′𝑔𝑔,𝑖𝑖 = 𝜙𝜙 ′𝑔𝑔,𝑖𝑖 + 𝛼𝛼𝑛𝑛 �𝜙𝜙 ′𝑔𝑔,𝑖𝑖 − 𝜙𝜙 ′𝑔𝑔,𝑖𝑖 + 𝛽𝛽 �𝜙𝜙𝑔𝑔,𝑖𝑖− 𝜙𝜙 ′𝑔𝑔,𝑖𝑖 ��(99)

𝜙𝜙は n 回目の外部反復後の中性子束で、
𝛼𝛼𝑛𝑛 、
ここで、𝜙𝜙 ′は Chebyshev 外挿法で補正した中性子束、

𝛽𝛽は

2

𝛽𝛽 =

2 − 𝑘𝑘𝑚𝑚𝑚𝑚𝑚𝑚 − 𝑘𝑘𝑚𝑚𝑚𝑚𝑚𝑚

𝛼𝛼𝑛𝑛 = 2𝛾𝛾𝛼𝛼0 = 1γ=

𝑇𝑇𝑛𝑛−1 (𝛾𝛾)1=𝛼𝛼𝑇𝑇𝑛𝑛 (𝛾𝛾)1 − 𝑛𝑛−14𝛾𝛾 2

2 − 𝑘𝑘𝑚𝑚𝑚𝑚𝑚𝑚 − 𝑘𝑘𝑚𝑚𝑚𝑚𝑚𝑚
𝑘𝑘𝑚𝑚𝑚𝑚𝑚𝑚 − 𝑘𝑘𝑚𝑚𝑚𝑚𝑚𝑚

(100)

(101)
(102)
(103)

であり、𝑘𝑘𝑚𝑚𝑚𝑚𝑚𝑚 、𝑘𝑘𝑚𝑚𝑚𝑚𝑚𝑚 は n 回目までの実効増倍率の内、最大と最小のものである 25)。Chebyshev

外挿法では、𝛼𝛼𝑛𝑛 、𝛽𝛽の設定が計算時間に大きな影響を与えることから、高速化のために 3.3.3

節で説明するドミナンス比を用いて𝛼𝛼𝑛𝑛 、𝛽𝛽を計算する手法 26)も提案されている。

Chebyshev 外挿法は SOR 法に比べて一般的に反復回数を低減することが出来るが、n-1 回目の中性子束だけでなく、n-2 回目の中性子束も保存しておかなければならないため、メモリー量の負担が大きくなる。近年の計算機性能の向上に伴い、メモリー容量の制約はあまり大きな
デメリットではないが、Chebyshev 外挿法を適用する場合にはメモリー量の増加に留意する必要がある。
3.3.3 Wielandt 法
実効増倍率の収束性は(104)式のように最大固有値𝜆𝜆0 と二番目の固有値𝜆𝜆1 の比(𝜆𝜆1⁄𝜆𝜆0 )である

ドミナンス比に依存する。ドミナンス比を用いると、予想される外部反復回数 n は次式のように底をドミナンス比とする収束条件𝜀𝜀の対数で書き表される。
𝑛𝑛 = log𝜆𝜆1⁄𝜆𝜆0 𝜀𝜀(104)

例えば、ドミナンス比が 0.1 のように極めて小さい場合、収束判定基準を 10-5 とすると、予想される外部反復回数は(104)式より 5 回となる。このようにドミナンス比が小さい条件では少ない外部反復回数で十分な計算精度を得ることが出来る。しかし、大型の軽水炉ではドミナンス
比は 0.97～0.99 程度になることがある 13)。この場合、収束判定基準を 10-5 とすると、予想される外部反復回数は 380～1150 回と非常に多くの外部反復が必要となる。
Wielandt 法は、(105)式のようにドミナンス比を小さくなるように拡散方程式の形を変えることにより、外部反復回数を低減する手法である 27)。

27

−𝐷𝐷𝑔𝑔,𝑖𝑖𝜒𝜒𝑔𝑔,𝑖𝑖𝜒𝜒𝑔𝑔,𝑖𝑖d2 𝑛𝑛
𝑛𝑛
𝜙𝜙𝑔𝑔,𝑖𝑖 + �Σ𝑟𝑟,𝑔𝑔,𝑖𝑖 −𝜈𝜈Σ𝑓𝑓,𝑔𝑔,𝑖𝑖 � 𝜙𝜙𝑔𝑔,𝑖𝑖−� 𝜈𝜈Σ𝑓𝑓,𝑔𝑔′ ,𝑖𝑖 𝜙𝜙𝑔𝑔𝑛𝑛′ ,𝑖𝑖
2
d𝑥𝑥
𝜆𝜆𝑒𝑒
𝜆𝜆𝑒𝑒 ′𝑔𝑔 ≠𝑔𝑔1 1𝑛𝑛= � − � 𝜒𝜒𝑔𝑔,𝑖𝑖 � 𝜈𝜈Σ𝑓𝑓,𝑔𝑔′ ,𝑖𝑖 𝜙𝜙𝑔𝑔𝑛𝑛−1′ ,𝑖𝑖 + � Σ𝑠𝑠,𝑔𝑔′ →𝑔𝑔,𝑖𝑖 𝜙𝜙𝑔𝑔′ ,𝑖𝑖
𝑘𝑘 𝜆𝜆𝑒𝑒′′(105)

𝑔𝑔 <𝑔𝑔𝑔𝑔+ � Σ𝑠𝑠,𝑔𝑔′ →𝑔𝑔,𝑖𝑖 𝜙𝜙𝑔𝑔𝑛𝑛−1′ ,𝑖𝑖
𝑔𝑔′ >𝑔𝑔

ここで、𝜆𝜆𝑒𝑒 は実効増倍率の推定値である。(105)式を用いると、ドミナンス比は11𝜆𝜆′1 1⁄𝜆𝜆′0 𝜆𝜆0 − 𝜆𝜆𝑒𝑒 𝜆𝜆1 𝜆𝜆𝑒𝑒 − 𝜆𝜆0
=
=
= ∙
1
1
𝜆𝜆′0 1⁄𝜆𝜆′1
𝜆𝜆0 𝜆𝜆𝑒𝑒 − 𝜆𝜆1−
𝜆𝜆1 𝜆𝜆𝑒𝑒
𝜆𝜆 −𝜆𝜆0となり、元のドミナンス比の 𝑒𝑒

𝜆𝜆𝑒𝑒 −𝜆𝜆1

𝜆𝜆𝑒𝑒 > 𝜆𝜆0(106)倍となる。なお、𝜆𝜆𝑒𝑒 は(107)でなければドミナンス比が負の値となり収束しない場合があるので、𝜆𝜆𝑒𝑒 の設定には注意が必要である。

ここで例として、𝜆𝜆0 = 1.0、𝜆𝜆1 = 0.97、𝜆𝜆𝑒𝑒 = 1.05とすると、Wielandt 法を適用した拡散方程式のドミナンス比は元のドミナンス比の 0.625 倍となり、外部反復回数は 380 回から 23 回へと大幅に低減される。このように Wielandt 法はドミナンス比を小さくしてやることにより、外部
反復回数を低減することが出来る。但し、(105)式のように除去断面積の項が小さくなり、行列の対角成分の項が𝜒𝜒𝑔𝑔,𝑖𝑖 𝜈𝜈Σ𝑓𝑓,𝑖𝑖 ⁄𝜆𝜆𝑒𝑒 の分だけ小さくなることに注意が必要である。内部反復の計算に

SOR 法などの反復法を用いる場合、内部反復の収束性は対角成分の大きさに依存する。Wielandt法では対角成分が小さくなることから、内部反復の反復回数が増加してしまう。また(105)式では一般的な多群の拡散方程式に Wielandt 法を適用したが、右辺に中性子源の項が入るため、元々の式よりも左辺に 0 ではない中性子束の係数が増え、内部反復の収束性が悪化することが懸念される 38。そのため、計算体系によっては全体の計算時間はあまり改善しない可能性がある。
また、この内部反復の収束性が悪化する問題を解決するため、次式のように(105)式の左辺第三項を右辺の散乱項に移項することを考える。
−𝐷𝐷𝑔𝑔,𝑖𝑖𝜒𝜒𝑔𝑔,𝑖𝑖d2 𝑛𝑛
𝑛𝑛
𝜙𝜙𝑔𝑔,𝑖𝑖 + �Σ𝑟𝑟,𝑔𝑔,𝑖𝑖 −𝜈𝜈Σ𝑓𝑓,𝑔𝑔,𝑖𝑖 � 𝜙𝜙𝑔𝑔,𝑖𝑖
2
d𝑥𝑥
𝜆𝜆𝑒𝑒

𝜒𝜒𝑔𝑔,𝑖𝑖
1 1= � − � 𝜒𝜒𝑔𝑔,𝑖𝑖 � 𝜈𝜈Σ𝑓𝑓,𝑔𝑔′ ,𝑖𝑖 𝜙𝜙𝑔𝑔𝑛𝑛−1𝜈𝜈Σ𝑓𝑓,𝑔𝑔′ ,𝑖𝑖 � 𝜙𝜙𝑔𝑔𝑛𝑛′ ,𝑖𝑖′ ,𝑖𝑖 + � �Σ𝑠𝑠,𝑔𝑔′ →𝑔𝑔,𝑖𝑖 +𝜆𝜆
𝑘𝑘 𝜆𝜆𝑒𝑒𝑒𝑒′′𝑔𝑔+ � �Σ𝑠𝑠,𝑔𝑔′ →𝑔𝑔,𝑖𝑖 +𝑔𝑔′ >𝑔𝑔(108)

𝑔𝑔 <𝑔𝑔𝜒𝜒𝑔𝑔,𝑖𝑖𝜈𝜈Σ𝑓𝑓,𝑔𝑔′ ,𝑖𝑖 � 𝜙𝜙𝑔𝑔𝑛𝑛−1
′ ,𝑖𝑖
𝜆𝜆𝑒𝑒

この場合、内部反復の収束性の悪化は抑えることが出来るものの、外部反復の収束性が悪化し

38

例えば(35)式では七点階差の連立一次方程式であったが、(105)式を用いると全てのエネル
ギー群の𝜙𝜙𝑔𝑔𝑛𝑛′ ,𝑖𝑖 が左辺に含まれることになる。

28

てしまう。3.2.2.2 節で述べたように、散乱反応では上方散乱はエネルギーの低い領域でしか起きないため、高いエネルギー群から順番に連立一次方程式を解くことで外部反復の収束性を向上させることが出来る。しかし、(108)式のように散乱源に核分裂の項が含まれると、自群よりも低いエネルギー群からの中性子の発生を全てのエネルギー群で考慮する必要が出てくるため、外部反復の収束性が悪化する。
Wielandt 法は後述する粗メッシュ加速法のように、一群の拡散方程式を解く場合に効果的である。(105)式より、一群の拡散方程式は−𝐷𝐷𝑖𝑖𝜈𝜈Σ𝑓𝑓,𝑖𝑖 𝑛𝑛
d2 𝑛𝑛
1 1
𝜙𝜙𝑖𝑖 + �Σ𝑎𝑎,𝑖𝑖 −� 𝜙𝜙𝑖𝑖 = � − � 𝜈𝜈Σ𝑓𝑓,𝑖𝑖 𝜙𝜙𝑖𝑖𝑛𝑛−1
2
d𝑥𝑥
𝜆𝜆𝑒𝑒
𝑘𝑘 𝜆𝜆𝑒𝑒(109)となる。ここで、Σ𝑎𝑎 は巨視的吸収断面積である。一群拡散方程式では(105)式のような中性子源の項が左辺に無いため、外部反復の収束性が悪化せずに Wiealdt 法によるドミナンス比の低減の効果だけを得ることができる。なお、この場合でも(109)式を用いると行列の対角成分が𝜈𝜈Σ𝑓𝑓,𝑖𝑖 ⁄𝜆𝜆𝑒𝑒 の分だけ小さくなるため、内部反復の収束性は悪化する。

3.4 粗メッシュ加速法
今まで説明してきた高速化手法は、Wielandt 法を除き、以前の反復計算での中性子束を用いて一回の反復計算での中性子束の更新量を大きくすることで、収束性を向上させていた。粗メッシュ加速法では、解析対象としている詳細なメッシュ数よりも少ない(粗い)メッシュ数の計
算結果を用いて中性子束を更新している。
拡散計算や輸送計算の計算時間はメッシュ数やエネルギー群数にほぼ比例する。例として、
BWR において、集合体単位で均質化した場合の炉心計算と、燃料棒単位で均質化した場合の炉心計算の計算時間を比較する。炉心内の燃料が全て 9×9 燃料集合体で構成されていると仮定する。チャンネルボックスとウォーターギャップ領域をまとめて 1 メッシュとすると、集合体内のメッシュは(1+9+1)×(1+9+1)=11×11=121 メッシュとなる。そのため、集合体単位で均質化した場合の炉心計算は燃料棒単位で均質化した場合の炉心計算に比べて 100 倍以上の高速化が期待できる。
粗メッシュ加速計算では、メッシュ数だけでなくエネルギー群数も縮約していることが多い。
そのため、粗メッシュ計算と元々の詳細メッシュ計算の計算時間の差異はメッシュ数の差以上に大きくなることが多く、非常に効果的な手法として多くの核計算コードで利用されている。
以降では粗メッシュ加速法として広く使われている、Coarse Mesh Rebalance(CMR)法 28-30)、
Coarse Mesh Finite Difference (CMFD)法 29-31)と Generalized Coarse Mesh Rebalance (GCMR)法 4)について説明する 39。なお、前述の通り、実際の粗メッシュ加速法はエネルギー群も縮約することが多いが、ここでは簡単のため、エネルギー群数については詳細メッシュと同じとする 40。

39

粗メッシュではないが、同様の手法として、SN 計算で拡散計算の結果を用いて加速する拡散合成加速(Diffusion Synthetic Acceleration)法がある 32)。
40 式の上は詳細群の g 群から縮約群の G 群に代わるだけであり、ほとんど変わらない。
29

3.4.1 CMR 法
𝑜𝑜𝑜𝑜𝑜𝑜
CMR 法では、エネルギーg 群、メッシュ i の表面流入中性子流及び表面流出中性子流𝐽𝐽𝑠𝑠,𝑔𝑔,𝑖𝑖、

𝑖𝑖𝑖𝑖を用いた以下のバランス方程式を利用する。
𝐽𝐽𝑠𝑠,𝑔𝑔,𝑖𝑖

𝑜𝑜𝑜𝑜𝑜𝑜𝑖𝑖𝑖𝑖
𝐽𝐽𝑠𝑠,𝑔𝑔,𝑖𝑖− 𝐽𝐽𝑠𝑠,𝑔𝑔,𝑖𝑖
𝜒𝜒𝑔𝑔,𝑖𝑖+ Σ𝑟𝑟,𝑔𝑔,𝑖𝑖 𝜙𝜙𝑔𝑔,𝑖𝑖 =� 𝜈𝜈Σ ′ 𝜙𝜙 ′ + � Σ𝑠𝑠,𝑔𝑔′ →𝑔𝑔,𝑖𝑖 𝜙𝜙𝑔𝑔′ ,𝑖𝑖∆𝑥𝑥𝑖𝑖
𝑘𝑘 ′ 𝑓𝑓,𝑔𝑔 ,𝑖𝑖 𝑔𝑔 ,𝑖𝑖
′
𝑔𝑔

(110)

𝑔𝑔 ≠𝑔𝑔𝑜𝑜𝑜𝑜𝑜𝑜𝑖𝑖𝑖𝑖、𝐽𝐽𝑠𝑠,𝑔𝑔,𝑖𝑖は部分表面中性子流を用いて次式で書き表される。
𝐽𝐽𝑠𝑠,𝑔𝑔,𝑖𝑖
𝑜𝑜𝑜𝑜𝑜𝑜+−+−𝐽𝐽𝑠𝑠,𝑔𝑔,𝑖𝑖= 𝐽𝐽𝑠𝑠+,𝑔𝑔,𝑖𝑖+ 𝐽𝐽𝑠𝑠−,𝑔𝑔,𝑖𝑖= 𝐽𝐽𝑠𝑠−,𝑔𝑔,𝑖𝑖+1+ 𝐽𝐽𝑠𝑠+,𝑔𝑔,𝑖𝑖−1
𝑖𝑖𝑖𝑖
+
−
+
−
𝐽𝐽𝑠𝑠,𝑔𝑔,𝑖𝑖= 𝐽𝐽𝑠𝑠−,𝑔𝑔,𝑖𝑖+ 𝐽𝐽𝑠𝑠+,𝑔𝑔,𝑖𝑖= 𝐽𝐽𝑠𝑠+,𝑔𝑔,𝑖𝑖−1+ 𝐽𝐽𝑠𝑠−,𝑔𝑔,𝑖𝑖+1

(111)
(112)

+
−
ここで、𝐽𝐽𝑠𝑠+,𝑔𝑔,𝑖𝑖
はエネルギーg 群、メッシュ i の正側の表面の表面放出中性子流であり、𝐽𝐽𝑠𝑠+,𝑔𝑔,𝑖𝑖
はエネルギーg 群、メッシュ i の正側の表面の表面入射中性子流である。また中性子流の連続性から、(111)、(112)式の部分表面中性子流には次の関係がある。
±±
𝐽𝐽𝑠𝑠+,𝑔𝑔,𝑖𝑖= 𝐽𝐽𝑠𝑠−,𝑔𝑔,𝑖𝑖+1(113)(110)～(113)式中の中性子束、部分表面中性子流の関係は言葉だけでは理解しにくいと思う。そこで(110)～(113)式中の中性子束、部分表面中性子流の関係をイメージしやすくするため、図 6に中性子束、表面中性子束、部分表面中性子流の関係を図示する。

図 6 中性子束、表面中性子束、部分表面中性子流の関係(110)式は(87)式をメッシュ i について体積積分し、拡散項をガウスの発散定理を用いて面積分に変えることによって得られ、次式のように(1)式の右辺第一項の拡散項(漏洩項)を表面流入中性子流と表面流出中性子流で置き換えただけである。
−𝐷𝐷𝑔𝑔,𝑖𝑖

𝑜𝑜𝑜𝑜𝑜𝑜𝑖𝑖𝑖𝑖
𝐽𝐽𝑠𝑠,𝑔𝑔,𝑖𝑖− 𝐽𝐽𝑠𝑠,𝑔𝑔,𝑖𝑖
d2
𝜙𝜙=
d𝑥𝑥 2 𝑔𝑔,𝑖𝑖
∆𝑥𝑥𝐼𝐼

(114)

n 回目の外部反復終了時を考えると、(110)式は𝑜𝑜𝑜𝑜𝑜𝑜,𝑛𝑛𝑖𝑖𝑖𝑖,𝑛𝑛
𝐽𝐽𝑠𝑠,𝑔𝑔,𝑖𝑖− 𝐽𝐽𝑠𝑠,𝑔𝑔,𝑖𝑖

∆𝑥𝑥𝐼𝐼

𝑛𝑛
+ Σ𝑟𝑟,𝑔𝑔,𝑖𝑖 𝜙𝜙𝑖𝑖,𝑔𝑔=

𝜒𝜒𝑔𝑔,𝑖𝑖� 𝜈𝜈Σ ′ 𝜙𝜙 𝑛𝑛−1
+ � Σ𝑠𝑠,𝑔𝑔′ →𝑔𝑔,𝑖𝑖 𝜙𝜙𝑔𝑔𝑛𝑛′ ,𝑖𝑖′
𝑘𝑘 𝑛𝑛−1 ′ 𝑓𝑓,𝑔𝑔 ,𝑖𝑖 𝑔𝑔 ,𝑖𝑖
′
𝑔𝑔

(115)

𝑔𝑔 ≠𝑔𝑔となる。3.2.1.6 節で述べたように、一般的に内部反復を完全に収束させると、外部反復の回数は少なくなるものの、内部反復の計算時間が大幅に増加するため、全体の計算時間は増加することが多い。このように外部反復の収束途中では内部反復は完全に収束していないため、(115)30

式は成立しない。
そこで CMR 法では、収束途中でも(115)式が成立するように、中性子源項以外にリバランス
因子𝑓𝑓𝑔𝑔,𝑖𝑖 を掛ける。

+,𝑛𝑛
−,𝑛𝑛
+,𝑛𝑛
−,𝑛𝑛
𝑛𝑛
�𝐽𝐽𝑠𝑠+,𝑔𝑔,𝑖𝑖+ 𝐽𝐽𝑠𝑠−,𝑔𝑔,𝑖𝑖+ Σ𝑟𝑟,𝑔𝑔,𝑖𝑖 𝜙𝜙𝑔𝑔,𝑖𝑖
�𝑓𝑓𝑔𝑔,𝑖𝑖 − 𝐽𝐽𝑠𝑠+,𝑔𝑔,𝑖𝑖−1
𝑓𝑓𝑔𝑔,𝑖𝑖−1 − 𝐽𝐽𝑠𝑠−,𝑔𝑔,𝑖𝑖+1
𝑓𝑓𝑔𝑔,𝑖𝑖+1𝜒𝜒𝑔𝑔,𝑖𝑖𝑛𝑛= 𝑛𝑛−1 � 𝜈𝜈Σ𝑓𝑓,𝑔𝑔′ ,𝑖𝑖 𝜙𝜙𝑔𝑔𝑛𝑛−1′ ,𝑖𝑖 + � Σ𝑠𝑠,𝑔𝑔′ →𝑔𝑔,𝑖𝑖 𝜙𝜙𝑔𝑔′ ,𝑖𝑖
𝑘𝑘
′
′
𝑔𝑔

(116)

𝑔𝑔 ≠𝑔𝑔𝑜𝑜𝑜𝑜𝑜𝑜,𝑛𝑛𝑖𝑖𝑖𝑖,𝑛𝑛
ここで、(115)式中の𝐽𝐽𝑠𝑠,𝑔𝑔,𝑖𝑖、𝐽𝐽𝑠𝑠,𝑔𝑔,𝑖𝑖
は(111)、(112)式を用いて部分表面中性子流に置き換えた。

(116)式で実効増倍率、中性子束、断面積等をリバランス因子の係数と見なし、リバランス因子を計算することが出来れば、次式を用いて中性子束を更新することで内部反復を完全に収束させた場合と同等の中性子束を得ることが出来る。
𝑛𝑛

𝑛𝑛
𝜙𝜙 ′𝑔𝑔,𝑖𝑖 = 𝜙𝜙𝑔𝑔,𝑖𝑖
𝑓𝑓𝑔𝑔,𝑖𝑖(117)

しかし、詳細メッシュ体系で(116)式を用いてリバランス因子を計算しても、通常の内部反復計
算を完全に収束させるのと計算時間はあまり変わらないため、CMR 法を用いるメリットは無い。CMR 法では、(116)式の計算を、(118)、(119)式のように粗メッシュ体系で実施する。
+
𝑛𝑛
+
−
−
�𝐽𝐽𝑠𝑠+,𝑔𝑔,𝐼𝐼+ 𝐽𝐽𝑠𝑠−,𝑔𝑔,𝐼𝐼+ Σ𝑟𝑟,𝑔𝑔,𝐼𝐼 𝜙𝜙𝑔𝑔,𝐼𝐼�𝑓𝑓𝑔𝑔,𝐼𝐼 − 𝐽𝐽𝑠𝑠+,𝑔𝑔,𝐼𝐼−1
𝑓𝑓𝑔𝑔,𝐼𝐼−1 − 𝐽𝐽𝑠𝑠−,𝑔𝑔,𝐼𝐼+1
𝑓𝑓𝑔𝑔,𝐼𝐼+1𝜒𝜒𝑔𝑔,𝐼𝐼𝑛𝑛= 𝑛𝑛−1 � 𝜈𝜈Σ𝑓𝑓,𝑔𝑔′ ,𝐼𝐼 𝜙𝜙𝑔𝑔𝑛𝑛−1′ ,𝐼𝐼 + � Σ𝑠𝑠,𝑔𝑔′ →𝑔𝑔,𝐼𝐼 𝜙𝜙𝑔𝑔′ ,𝐼𝐼
𝑘𝑘
′
′
𝑛𝑛
𝑛𝑛
𝜙𝜙 ′𝑔𝑔,𝑖𝑖 = 𝑓𝑓𝑔𝑔,𝐼𝐼 × 𝜙𝜙𝑔𝑔,𝑖𝑖, 𝑖𝑖 ∈ 𝐼𝐼𝑔𝑔

𝑔𝑔 ≠𝑔𝑔

(118)

(119)

ここで粗メッシュ体系の断面積や中性子束は詳細メッシュ体系の断面積、中性子束を用いて(2)、(4)式で計算される 41。粗メッシュ体系の計算は詳細メッシュ体系の計算に比べて圧倒的に高速なため、粗メッシュ体系で(118)式のバランス方程式を計算することで、ほとんど計算時間を増加させることなく、リバランス因子を計算することが出来る。このように CMR 法では粗メッシュ体系でバランス因子を計算することで、内部反復を完全に収束させた場合に近い中性子束を得ることができるため、外部反復回数の低減が可能になる。
拡散計算で CMR 法を用いる場合、(110)式ではなく、表面中性子流𝐽𝐽𝑠𝑠±,𝑔𝑔,𝐼𝐼 を用いる方が簡単である。(114)式と同様に中性子の漏洩項は表面中性子流𝐽𝐽𝑠𝑠±,𝑔𝑔,𝐼𝐼 を用いて−𝐷𝐷𝑔𝑔,𝐼𝐼

𝐽𝐽𝑠𝑠+,𝑔𝑔,𝐼𝐼 − 𝐽𝐽𝑠𝑠−,𝑔𝑔,𝐼𝐼d2
𝜙𝜙𝑔𝑔,𝐼𝐼 =
2
d𝑥𝑥
∆𝑥𝑥𝐼𝐼

(120)

で書き表される。なお、(114)式と(120)式は式としては似ているが、(114)式はメッシュ内から
𝑜𝑜𝑜𝑜𝑜𝑜𝑖𝑖𝑖𝑖の流出中性子流𝐽𝐽𝑠𝑠,𝑔𝑔,𝐼𝐼とメッシュ内への流入中性子流𝐽𝐽𝑠𝑠,𝑔𝑔,𝐼𝐼との差異であり、(120)式はメッシュの正側の表面の正方向の正味の中性子流𝐽𝐽𝑠𝑠+,𝑔𝑔,𝐼𝐼 とメッシュの負側の表面の正方向の正味の中性子

流𝐽𝐽𝑠𝑠−,𝑔𝑔,𝐼𝐼 の差異であり、各要素の物理的な意味は異なっている。𝐽𝐽𝑠𝑠+,𝑔𝑔,𝐼𝐼 と𝐽𝐽𝑠𝑠−,𝑔𝑔,𝐼𝐼 を(111)、(112)式のように部分中性子流で書き表すと、++−−
𝐽𝐽𝑠𝑠+,𝑔𝑔,𝐼𝐼 = 𝐽𝐽𝑠𝑠+,𝑔𝑔,𝐼𝐼− 𝐽𝐽𝑠𝑠+,𝑔𝑔,𝐼𝐼= 𝐽𝐽𝑠𝑠−,𝑔𝑔,𝐼𝐼+1− 𝐽𝐽𝑠𝑠−,𝑔𝑔,𝐼𝐼+1++−−
𝐽𝐽𝑠𝑠−,𝑔𝑔,𝐼𝐼 = 𝐽𝐽𝑠𝑠−,𝑔𝑔,𝐼𝐼− 𝐽𝐽𝑠𝑠−,𝑔𝑔,𝐼𝐼= 𝐽𝐽𝑠𝑠+,𝑔𝑔,𝐼𝐼−1− 𝐽𝐽𝑠𝑠+,𝑔𝑔,𝐼𝐼−1

(121)
(122)

CMR 法では粗メッシュ体系で実効増倍率や中性子束を計算しているのではなく、実効増倍率、中性子束、断面積等をリバランス因子の係数とみなし、リバランス因子を計算している。
この点が後述する CMFD 法や GCMR 法との大きな違いであり、CMR 法の特徴である。

4131

となり、(111)、(112)式と異なっていることが分かる。ただし、これらの式から分かるように、
𝑜𝑜𝑜𝑜𝑜𝑜𝑖𝑖𝑖𝑖
𝐽𝐽𝑠𝑠+,𝑔𝑔,𝐼𝐼 − 𝐽𝐽𝑠𝑠−,𝑔𝑔,𝐼𝐼 = 𝐽𝐽𝑠𝑠,𝑔𝑔,𝑖𝑖− 𝐽𝐽𝑠𝑠,𝑔𝑔,𝑖𝑖(123)

であり、どちらもメッシュ i から漏れ出る正味の中性子数を意味している。
表面中性子流𝐽𝐽𝑠𝑠±,𝑔𝑔,𝐼𝐼 を用いると、(110)式は𝐽𝐽𝑠𝑠+,𝑔𝑔,𝐼𝐼 − 𝐽𝐽𝑠𝑠−,𝑔𝑔,𝐼𝐼𝜒𝜒𝑔𝑔,𝐼𝐼+ Σ𝑟𝑟,𝑔𝑔,𝐼𝐼 𝜙𝜙𝑔𝑔,𝐼𝐼 =� 𝜈𝜈Σ ′ 𝜙𝜙 ′ + � Σ𝑠𝑠,𝑔𝑔′ →𝑔𝑔,𝐼𝐼 𝜙𝜙𝑔𝑔′ ,𝐼𝐼∆𝑥𝑥𝐼𝐼
𝑘𝑘 ′ 𝑓𝑓,𝑔𝑔 ,𝐼𝐼 𝑔𝑔 ,𝐼𝐼
′
𝑔𝑔

(124)

𝑔𝑔 ≠𝑔𝑔となる。拡散計算での表面中性子流の計算について説明する。粗メッシュ I 中の中性子束が中心から表面まで一次関数状に変化し、メッシュ境界で表面中性子束が連続 42と仮定すると、エネルギーg 群、粗メッシュ I の正方向の表面中性子流𝐽𝐽𝑠𝑠±,𝑔𝑔,𝐼𝐼 は表面中性子束𝜙𝜙𝑠𝑠±,𝑔𝑔,𝐼𝐼 と粗メッシュ

I の中心の中性子束𝜙𝜙𝑔𝑔,𝐼𝐼 を用いて次のように書き表すことが出来る。
𝜙𝜙𝑠𝑠+,𝑔𝑔,𝐼𝐼 − 𝜙𝜙𝑔𝑔,𝐼𝐼𝜙𝜙𝑔𝑔,𝐼𝐼+1 − 𝜙𝜙𝑠𝑠+,𝑔𝑔,𝐼𝐼
𝐽𝐽𝑠𝑠+,𝑔𝑔,𝐼𝐼 = −𝐷𝐷𝑔𝑔,𝐼𝐼= −𝐷𝐷𝑔𝑔,𝐼𝐼+1
(125)
∆𝑥𝑥𝐼𝐼
∆𝑥𝑥𝐼𝐼+122
𝜙𝜙𝑠𝑠−,𝑔𝑔,𝐼𝐼 − 𝜙𝜙𝑔𝑔,𝐼𝐼−1
𝜙𝜙𝑔𝑔,𝐼𝐼 − 𝜙𝜙𝑠𝑠−,𝑔𝑔,𝐼𝐼
𝐽𝐽𝑠𝑠−,𝑔𝑔,𝐼𝐼 = −𝐷𝐷𝑔𝑔,𝐼𝐼−1= −𝐷𝐷𝑔𝑔,𝐼𝐼(126)∆𝑥𝑥𝐼𝐼−1∆𝑥𝑥𝐼𝐼
2
2
ここで、𝐷𝐷𝑔𝑔,𝐼𝐼 は拡散係数、∆𝑥𝑥𝐼𝐼 は粗メッシュ I のメッシュ幅である。(125)、(126)式から表面中性子束𝜙𝜙𝑠𝑠±,𝑔𝑔,𝐼𝐼 を除去し、整理すると
𝐽𝐽𝑠𝑠+,𝑔𝑔,𝐼𝐼 = −𝐽𝐽𝑠𝑠−,𝑔𝑔,𝐼𝐼 = −2𝐷𝐷𝑔𝑔,𝐼𝐼 𝐷𝐷𝑔𝑔,𝐼𝐼+1
�𝜙𝜙− 𝜙𝜙𝑔𝑔,𝐼𝐼 �
𝐷𝐷𝑔𝑔,𝐼𝐼 ∆𝑥𝑥𝐼𝐼+1 + 𝐷𝐷𝑔𝑔,𝐼𝐼+1 ∆𝑥𝑥𝐼𝐼 𝑔𝑔,𝐼𝐼+1(127)

2𝐷𝐷𝑔𝑔,𝐼𝐼−1 𝐷𝐷𝑔𝑔,𝐼𝐼�𝜙𝜙 − 𝜙𝜙𝑔𝑔,𝐼𝐼−1 �
𝐷𝐷𝑔𝑔,𝐼𝐼−1 ∆𝑥𝑥𝐼𝐼 + 𝐷𝐷𝑔𝑔,𝐼𝐼 ∆𝑥𝑥𝐼𝐼−1 𝑔𝑔,𝐼𝐼(128)となる。(127)、(128)式を(124)式に代入すると−

2𝐷𝐷𝑔𝑔,𝐼𝐼 𝐷𝐷𝑔𝑔,𝐼𝐼+1

�𝐷𝐷𝑔𝑔,𝐼𝐼 ∆𝑥𝑥𝐼𝐼+1 + 𝐷𝐷𝑔𝑔,𝐼𝐼+1 ∆𝑥𝑥𝐼𝐼 �∆𝑥𝑥𝐼𝐼+ �Σ𝑟𝑟,𝑔𝑔,𝐼𝐼 ++=𝜙𝜙𝑔𝑔,𝐼𝐼+1 −2𝐷𝐷𝑔𝑔,𝐼𝐼−1 𝐷𝐷𝑔𝑔,𝐼𝐼�𝐷𝐷𝑔𝑔,𝐼𝐼−1 ∆𝑥𝑥𝐼𝐼 + 𝐷𝐷𝑔𝑔,𝐼𝐼 ∆𝑥𝑥𝐼𝐼−1 �∆𝑥𝑥𝐼𝐼

2𝐷𝐷𝑔𝑔,𝐼𝐼−1 𝐷𝐷𝑔𝑔,𝐼𝐼�𝐷𝐷𝑔𝑔,𝐼𝐼−1 ∆𝑥𝑥𝐼𝐼 + 𝐷𝐷𝑔𝑔,𝐼𝐼 ∆𝑥𝑥𝐼𝐼−1 �∆𝑥𝑥𝐼𝐼

2𝐷𝐷𝑔𝑔,𝐼𝐼 𝐷𝐷𝑔𝑔,𝐼𝐼+1

�𝐷𝐷𝑔𝑔,𝐼𝐼 ∆𝑥𝑥𝐼𝐼+1 + 𝐷𝐷𝑔𝑔,𝐼𝐼+1 ∆𝑥𝑥𝐼𝐼 �∆𝑥𝑥𝐼𝐼

𝜙𝜙𝑔𝑔,𝐼𝐼−1(129)

� 𝜙𝜙𝑔𝑔,𝐼𝐼𝜒𝜒𝑔𝑔,𝐼𝐼� 𝜈𝜈Σ ′ 𝜙𝜙 ′ + � Σ𝑠𝑠,𝑔𝑔′ →𝑔𝑔,𝐼𝐼 𝜙𝜙𝑔𝑔′ ,𝐼𝐼𝑘𝑘 ′ 𝑓𝑓,𝑔𝑔 ,𝐼𝐼 𝑔𝑔 ,𝐼𝐼′𝑔𝑔

𝑔𝑔 ≠𝑔𝑔となる。ここで、(118)式のようにリバランス因子を用いると物理的にはメッシュ境界での表面中性子束は連続である。本資料では 2.2.2 節で不連続因子の説明をしたことから、不連続因子を適用していないことを示すため、あえてこのように表記
した。

42

32

−

2𝐷𝐷𝑔𝑔,𝐼𝐼 𝐷𝐷𝑔𝑔,𝐼𝐼+1

�𝐷𝐷𝑔𝑔,𝐼𝐼 ∆𝑥𝑥𝐼𝐼+1 + 𝐷𝐷𝑔𝑔,𝐼𝐼+1 ∆𝑥𝑥𝐼𝐼 �∆𝑥𝑥𝐼𝐼+ �Σ𝑟𝑟,𝑔𝑔,𝐼𝐼 ++=𝜙𝜙𝑔𝑔,𝐼𝐼+1 𝑓𝑓𝑔𝑔,𝐼𝐼+1 −2𝐷𝐷𝑔𝑔,𝐼𝐼−1 𝐷𝐷𝑔𝑔,𝐼𝐼�𝐷𝐷𝑔𝑔,𝐼𝐼−1 ∆𝑥𝑥𝐼𝐼 + 𝐷𝐷𝑔𝑔,𝐼𝐼 ∆𝑥𝑥𝐼𝐼−1 �∆𝑥𝑥𝐼𝐼

2𝐷𝐷𝑔𝑔,𝐼𝐼−1 𝐷𝐷𝑔𝑔,𝐼𝐼𝜙𝜙𝑔𝑔,𝐼𝐼−1 𝑓𝑓𝑔𝑔,𝐼𝐼−1

�𝐷𝐷𝑔𝑔,𝐼𝐼−1 ∆𝑥𝑥𝐼𝐼 + 𝐷𝐷𝑔𝑔,𝐼𝐼 ∆𝑥𝑥𝐼𝐼−1 �∆𝑥𝑥𝐼𝐼

2𝐷𝐷𝑔𝑔,𝐼𝐼 𝐷𝐷𝑔𝑔,𝐼𝐼+1

�𝐷𝐷𝑔𝑔,𝐼𝐼 ∆𝑥𝑥𝐼𝐼+1 + 𝐷𝐷𝑔𝑔,𝐼𝐼+1 ∆𝑥𝑥𝐼𝐼 �∆𝑥𝑥𝐼𝐼(130)

� 𝜙𝜙𝑔𝑔,𝐼𝐼 𝑓𝑓𝑔𝑔,𝐼𝐼𝜒𝜒𝑔𝑔,𝐼𝐼� 𝜈𝜈Σ ′ 𝜙𝜙 ′ + � Σ𝑠𝑠,𝑔𝑔′ →𝑔𝑔,𝐼𝐼 𝜙𝜙𝑔𝑔′ ,𝐼𝐼𝑘𝑘 ′ 𝑓𝑓,𝑔𝑔 ,𝐼𝐼 𝑔𝑔 ,𝐼𝐼′𝑔𝑔

𝑔𝑔 ≠𝑔𝑔となる。(130)式と(35)式を比較すると、計算式がほとんど同じであることが分かる。すなわち、三次元拡散計算で CMR 法を用いた場合の計算式は、(35)式にリバランス因子を掛けた
𝑦𝑦−

𝑥𝑥−
𝑧𝑧−
𝐴𝐴𝑔𝑔,𝐼𝐼,𝐽𝐽,𝐾𝐾
𝜙𝜙𝑔𝑔,𝐼𝐼−1,𝐽𝐽,𝐾𝐾 𝑓𝑓𝑔𝑔,𝐼𝐼−1,𝐽𝐽,𝐾𝐾 + 𝐴𝐴𝑔𝑔,𝐼𝐼,𝐽𝐽,𝐾𝐾 𝜙𝜙𝑔𝑔,𝐼𝐼,𝐽𝐽−1,𝐾𝐾 𝑓𝑓𝑔𝑔,𝐼𝐼,𝐽𝐽−1,𝐾𝐾 + 𝐴𝐴𝑔𝑔,𝐼𝐼,𝐽𝐽,𝐾𝐾
𝜙𝜙𝑔𝑔,𝐼𝐼,𝐽𝐽,𝐾𝐾−1 𝑓𝑓𝑔𝑔,𝐼𝐼,𝐽𝐽,𝐾𝐾−1
0𝑥𝑥++ 𝐴𝐴𝑔𝑔,𝐼𝐼,𝐽𝐽,𝐾𝐾
𝜙𝜙𝑔𝑔,𝐼𝐼,𝐽𝐽,𝐾𝐾 𝑓𝑓𝑔𝑔,𝐼𝐼,𝐽𝐽,𝐾𝐾 + 𝐴𝐴𝑔𝑔,𝐼𝐼,𝐽𝐽,𝐾𝐾
𝜙𝜙𝑔𝑔,𝐼𝐼+1,𝐽𝐽,𝐾𝐾 𝑓𝑓𝑔𝑔,𝐼𝐼+1,𝐽𝐽,𝐾𝐾
𝑦𝑦+

となる。

𝑧𝑧+
+ 𝐴𝐴𝑔𝑔,𝐼𝐼,𝐽𝐽,𝐾𝐾 𝜙𝜙𝑔𝑔,𝐼𝐼,𝐽𝐽+1,𝐾𝐾 𝑓𝑓𝑔𝑔,𝐼𝐼,𝐽𝐽+1,𝐾𝐾 + 𝐴𝐴𝑔𝑔,𝐼𝐼,𝐽𝐽,𝐾𝐾
𝜙𝜙𝑔𝑔,𝐼𝐼,𝐽𝐽,𝐾𝐾+1 𝑓𝑓𝑔𝑔,𝐼𝐼,𝐽𝐽,𝐾𝐾+1
𝜒𝜒𝑔𝑔,𝐼𝐼,𝐽𝐽,𝐾𝐾=� 𝜈𝜈Σ𝑓𝑓,𝑔𝑔′ ,𝐼𝐼,𝐽𝐽,𝐾𝐾 𝜙𝜙𝑔𝑔′ ,𝐼𝐼,𝐽𝐽,𝐾𝐾 + � Σ𝑠𝑠,𝑔𝑔′ →𝑔𝑔,𝐼𝐼,𝐽𝐽,𝐾𝐾 𝜙𝜙𝑔𝑔′ ,𝐼𝐼,𝐽𝐽,𝐾𝐾
𝑘𝑘
′
′
𝑔𝑔

(131)

𝑔𝑔

なお、繰り返しになるが、CMR 法では実効増倍率、中性子束、断面積等はリバランス因子の係数であり、(131)式は固有値方程式ではなくリバランス因子の連立一次方程式であることに注意する必要がある。そのため、CMR 法では固定源計算と同様に、粗メッシュ体系では外部反復を必要とせず、内部反復計算のみでリバランス因子を計算することができる。
3.4.2 CMFD 法
CMR 法では核分裂源項及び散乱源項を固定したバランス方程式を用いて、中性子束のみを修正していた。また、CMR 法では粗メッシュ体系の計算に(2)、(4)式を用いて計算した断面積、中性子束を用いており、空間均質化誤差を考慮していない。そのため、中性子強吸収体が含まれる場合など、空間均質化誤差の大きい体系に CMR 法を適用すると計算が発散する場合があり、収束安定性に問題があった。
CMFD 法は、粗メッシュ拡散計算の結果を利用して、詳細メッシュ計算の中性子束及び実効増倍率を修正する手法である。CMR 法と同様に、空間均質化誤差を考慮せずに単に均質化するだけでは空間均質化誤差の大きな体系では収束性の向上を望むことは出来ない。そこで
CMFD 法では、空間均質化誤差低減手法を用いて詳細メッシュ計算の結果を再現するように粗メッシュ拡散計算を実施している 43。
2.2.2 節で述べたように、詳細メッシュ計算を完全に収束させた状態で不連続因子を適用すると、粗メッシュ拡散計算の結果は詳細メッシュ計算の結果を再現する。では、詳細メッシュ計
算が収束していない状態で不連続因子を適用するとどうなるか？n 回目の外部反復の結果で得られた不連続因子を適用し、粗メッシュ拡散計算を行う場合、完全に収束させた詳細メッシュ
43

以降では空間均質化誤差低減手法として不連続因子を例に挙げて説明を行うが、あくまで例
であり、CMFD 法で用いている空間均質化誤差低減手法は不連続因子と同じではない。
33

計算の結果とは異なる結果になるが、近い結果になることが期待できる。3.2.1.6 節では、内部反復を完全に収束させると外部反復の回数が少なくなると述べたが、これと似たような理屈である。すなわち、CMFD 法では、n 回目の外部反復における中性子束と実効増倍率が最終的に到達すると期待される値を粗メッシュ拡散計算を用いて効率的に計算していることになる。
3.3 節で示した SOR 法や Chebyshev 外挿法では、以前の外部反復での中性子束を利用して未来の中性子束を予測しているが、粗メッシュ加速法ではより確からしい中性子束を計算することができる。さらに中性子束だけでなく実効増倍率まで更新されることから、これらの手法よりも収束性の向上が期待できる。また、CMFD 法は粗メッシュ拡散計算モジュールと空間均質化誤差低減手法を別途実装する必要はあるものの、計算コード自体に大きな改変を加えることなく実装できること、輸送計算・拡散計算のどちらにも適用できるなど適用範囲が広いこと、そして CMR 法に比べて収束安定性が良いことなどから、多くの集合体計算コードや、次世代の詳細メッシュ炉心計算コードで使われている。
CMFD 法を適用するためには、空間均質化誤差低減手法を適用した粗メッシュ体系の拡散計算を行う必要がある。CMR 法の説明と同様に、粗メッシュ I 中の中性子束が中心から表面まで
一次関数状に変化し、メッシュ境界で表面中性子束が連続と仮定すると、エネルギーg 群、粗メッシュ I の正方向の表面中性子流𝐽𝐽𝑠𝑠±,𝑔𝑔,𝐼𝐼 は表面中性子束𝜙𝜙𝑠𝑠±,𝑔𝑔,𝐼𝐼 と粗メッシュ I の中心の中性子

束𝜙𝜙𝐼𝐼,𝑔𝑔 を用いて次のように書き表すことが出来る。
𝜙𝜙𝑠𝑠+,𝑔𝑔,𝐼𝐼 − 𝜙𝜙𝑔𝑔,𝐼𝐼𝜙𝜙𝑔𝑔,𝐼𝐼+1 − 𝜙𝜙𝑠𝑠+,𝑔𝑔,𝐼𝐼
𝐽𝐽𝑠𝑠+,𝑔𝑔,𝐼𝐼 = −𝐷𝐷𝑔𝑔,𝐼𝐼= −𝐷𝐷𝑔𝑔,𝐼𝐼+1
(132)
∆𝑥𝑥𝐼𝐼
∆𝑥𝑥𝐼𝐼+122
𝜙𝜙𝑠𝑠−,𝑔𝑔,𝐼𝐼 − 𝜙𝜙𝑔𝑔,𝐼𝐼−1
𝜙𝜙𝑔𝑔,𝐼𝐼 − 𝜙𝜙𝑠𝑠−,𝑔𝑔,𝐼𝐼
𝐽𝐽𝑠𝑠−,𝑔𝑔,𝐼𝐼 = −𝐷𝐷𝑔𝑔,𝐼𝐼−1= −𝐷𝐷𝑔𝑔,𝐼𝐼(133)∆𝑥𝑥𝐼𝐼−1∆𝑥𝑥𝐼𝐼
2
2
ここで、𝐷𝐷𝑔𝑔,𝐼𝐼 は拡散係数、∆𝑥𝑥𝐼𝐼 は粗メッシュ I のメッシュ幅である。(132)、(133)式を整理すると
𝐽𝐽𝑠𝑠+,𝑔𝑔,𝐼𝐼 = −𝐽𝐽𝑠𝑠−,𝑔𝑔,𝐼𝐼 = −2𝐷𝐷𝑔𝑔,𝐼𝐼
2𝐷𝐷𝑔𝑔,𝐼𝐼
2𝐷𝐷𝑔𝑔,𝐼𝐼+1
2𝐷𝐷𝑔𝑔,𝐼𝐼+1
𝜙𝜙+
𝜙𝜙 = −
𝜙𝜙+
𝜙𝜙∆𝑥𝑥𝐼𝐼 𝑠𝑠+,𝑔𝑔,𝐼𝐼∆𝑥𝑥𝐼𝐼 𝑔𝑔,𝐼𝐼∆𝑥𝑥𝐼𝐼+1 𝑔𝑔,𝐼𝐼+1∆𝑥𝑥𝐼𝐼+1 𝑠𝑠+,𝑔𝑔,𝐼𝐼
2𝐷𝐷𝑔𝑔,𝐼𝐼−1
2𝐷𝐷𝑔𝑔,𝐼𝐼−1
2𝐷𝐷𝑔𝑔,𝐼𝐼
2𝐷𝐷𝑔𝑔,𝐼𝐼
𝜙𝜙
+
𝜙𝜙
=−
𝜙𝜙 +
𝜙𝜙
∆𝑥𝑥𝐼𝐼−1 𝑠𝑠−,𝑔𝑔,𝐼𝐼∆𝑥𝑥𝐼𝐼−1 𝑔𝑔,𝐼𝐼−1∆𝑥𝑥𝐼𝐼 𝑔𝑔,𝐼𝐼∆𝑥𝑥𝐼𝐼 𝑠𝑠−,𝑔𝑔,𝐼𝐼

(134)

(135)

となる。(134)、(135)式を用いた粗メッシュ拡散計算では、詳細メッシュ体系における中性子の入射に対する放出の関係を保存しないため、均質化誤差が発生する。そこで、詳細メッシュ体
系における中性子の入射に対する放出の関係を保存させるため、中性子束に掛かっている係数
�𝑠𝑠±,𝑔𝑔,𝐼𝐼± に置き換える。
を中性子流の補正係数𝐷𝐷
𝐽𝐽𝑠𝑠+,𝑔𝑔,𝐼𝐼 = −𝐽𝐽𝑠𝑠−,𝑔𝑔,𝐼𝐼 = −2𝐷𝐷𝑔𝑔,𝐼𝐼
2𝐷𝐷
�𝑠𝑠+,𝑔𝑔,𝐼𝐼 𝜙𝜙𝑔𝑔,𝐼𝐼 = −𝐷𝐷
�𝑠𝑠−,𝑔𝑔,𝐼𝐼+1 𝜙𝜙𝑔𝑔,𝐼𝐼+1 + 𝑔𝑔,𝐼𝐼+1 𝜙𝜙𝑠𝑠+,𝑔𝑔,𝐼𝐼𝜙𝜙+ 𝐷𝐷∆𝑥𝑥𝐼𝐼 𝑠𝑠+,𝑔𝑔,𝐼𝐼∆𝑥𝑥𝐼𝐼+1
2𝐷𝐷𝑔𝑔,𝐼𝐼−12𝐷𝐷
�𝑠𝑠+,𝑔𝑔,𝐼𝐼−1 𝜙𝜙𝑔𝑔,𝐼𝐼−1 = −𝐷𝐷
�𝑠𝑠−,𝑔𝑔,𝐼𝐼 𝜙𝜙𝑔𝑔,𝐼𝐼 + 𝑔𝑔,𝐼𝐼 𝜙𝜙𝑠𝑠−,𝑔𝑔,𝐼𝐼𝜙𝜙+ 𝐷𝐷∆𝑥𝑥𝐼𝐼−1 𝑠𝑠−,𝑔𝑔,𝐼𝐼
∆𝑥𝑥𝐼𝐼

(136)

(137)

�𝑠𝑠−,𝑔𝑔,𝐼𝐼 は粗メッシュ I の負側の中�𝑠𝑠+,𝑔𝑔,𝐼𝐼 は粗メッシュ I の正側の中性子流の補正係数、𝐷𝐷
ここで、𝐷𝐷性子流の補正係数である。(132)～(137)式中の中性子束、表面中性子束、表面中性子流、中性子流の補正係数の関係をイメージしやすくするため、図 7 に中性子束、表面中性子束、表面中性
子流、中性子流の補正係数の関係を図示する。
34

図 7 中性子束、表面中性子束、表面中性子流の補正係数の関係(136)、(137)式より、中性子流の補正係数は

�𝑠𝑠+,𝑔𝑔,𝐼𝐼 =𝐷𝐷

𝐽𝐽𝑠𝑠+,𝑔𝑔,𝐼𝐼 +�𝑠𝑠−,𝑔𝑔,𝐼𝐼+1 =𝐷𝐷
�𝑠𝑠+,𝑔𝑔,𝐼𝐼−1 =𝐷𝐷
�𝑠𝑠−,𝑔𝑔,𝐼𝐼 =𝐷𝐷

2𝐷𝐷𝑔𝑔,𝐼𝐼𝜙𝜙∆𝑥𝑥𝐼𝐼 𝑠𝑠+,𝑔𝑔,𝐼𝐼𝜙𝜙𝑔𝑔,𝐼𝐼

2𝐷𝐷𝑔𝑔,𝐼𝐼+1
𝜙𝜙∆𝑥𝑥𝐼𝐼+1 𝑠𝑠+,𝑔𝑔,𝐼𝐼𝜙𝜙𝑔𝑔,𝐼𝐼+1−𝐽𝐽𝑠𝑠+,𝑔𝑔,𝐼𝐼 +𝐽𝐽𝑠𝑠−,𝑔𝑔,𝐼𝐼 +2𝐷𝐷𝑔𝑔,𝐼𝐼−1
𝜙𝜙∆𝑥𝑥𝐼𝐼−1 𝑠𝑠−,𝑔𝑔,𝐼𝐼𝜙𝜙𝑔𝑔,𝐼𝐼−1

2𝐷𝐷𝑔𝑔,𝐼𝐼𝜙𝜙∆𝑥𝑥𝐼𝐼 𝑠𝑠−,𝑔𝑔,𝐼𝐼𝜙𝜙𝑔𝑔,𝐼𝐼−𝐽𝐽𝑠𝑠−,𝑔𝑔,𝐼𝐼 +

(138)

(139)

(140)

(141)

となる。CMFD 法では SPH 法や不連続因子と同様に、中性子流の補正係数を詳細メッシュの中性子流、表面中性子束から計算する。三次元体系では、表面中性子流、表面中性子束は
𝐽𝐽𝑠𝑠±,𝑔𝑔,𝐼𝐼 =∑𝑗𝑗∈𝐽𝐽 ∑𝑘𝑘∈𝐾𝐾 𝐽𝐽𝑠𝑠±,𝑔𝑔,𝐼𝐼,𝑗𝑗,𝑘𝑘 𝑆𝑆𝑠𝑠±,𝑗𝑗,𝑘𝑘∑𝑗𝑗∈𝐽𝐽 ∑𝑘𝑘∈𝐾𝐾 𝑆𝑆𝑠𝑠±,𝑗𝑗,𝑘𝑘

𝜙𝜙𝑠𝑠±,𝑔𝑔,𝐼𝐼 =∑𝑗𝑗∈𝐽𝐽 ∑𝑘𝑘∈𝐾𝐾 𝜙𝜙𝑠𝑠±,𝑔𝑔,𝐼𝐼,𝑗𝑗,𝑘𝑘 𝑆𝑆𝑠𝑠±,𝑗𝑗,𝑘𝑘∑𝑗𝑗∈𝐽𝐽 ∑𝑘𝑘∈𝐾𝐾 𝑆𝑆𝑠𝑠±,𝑗𝑗,𝑘𝑘

(142)

(143)

と粗メッシュ I の表面の中性子流、中性子束の面積平均で表される。ここで、𝑆𝑆𝑠𝑠±,𝑗𝑗,𝑘𝑘 は粗メッ

シュ I の表面での y 方向に j 番目、z 方向に k 番目の詳細メッシュの面積である。なお、拡散計
算では、詳細メッシュの表面中性子流、表面中性子束は(132)、(133)式を用いて求めることができる。

中性子流の補正係数が得られたので、CMFD 法を用いた粗メッシュの拡散方程式を導出する。
(136)、(137)式から表面中性子束を消去すると、
𝐽𝐽𝑠𝑠+,𝑔𝑔,𝐼𝐼 =�𝑠𝑠+,𝑔𝑔,𝐼𝐼�𝑠𝑠−,𝑔𝑔,𝐼𝐼+1∆𝑥𝑥𝐼𝐼 𝐷𝐷𝑔𝑔,𝐼𝐼+1 𝐷𝐷∆𝑥𝑥𝐼𝐼+1 𝐷𝐷𝑔𝑔,𝐼𝐼 𝐷𝐷
𝜙𝜙𝑔𝑔,𝐼𝐼 −𝜙𝜙
𝐷𝐷𝑔𝑔,𝐼𝐼 ∆𝑥𝑥𝐼𝐼+1 + 𝐷𝐷𝑔𝑔,𝐼𝐼+1 ∆𝑥𝑥𝐼𝐼
𝐷𝐷𝑔𝑔,𝐼𝐼 ∆𝑥𝑥𝐼𝐼+1 + 𝐷𝐷𝑔𝑔,𝐼𝐼+1 ∆𝑥𝑥𝐼𝐼 𝑔𝑔,𝐼𝐼+1

(144)

35

𝐽𝐽𝑠𝑠−,𝑔𝑔,𝐼𝐼 =�𝑠𝑠+,𝑔𝑔,𝐼𝐼−1
�𝑠𝑠−,𝑔𝑔,𝐼𝐼∆𝑥𝑥𝐼𝐼−1 𝐷𝐷𝑔𝑔,𝐼𝐼 𝐷𝐷∆𝑥𝑥𝐼𝐼 𝐷𝐷𝑔𝑔,𝐼𝐼−1 𝐷𝐷
𝜙𝜙𝑔𝑔,𝐼𝐼−1 −𝜙𝜙𝐷𝐷𝑔𝑔,𝐼𝐼−1 ∆𝑥𝑥𝐼𝐼 + 𝐷𝐷𝑔𝑔,𝐼𝐼 ∆𝑥𝑥𝐼𝐼−1
𝐷𝐷𝑔𝑔,𝐼𝐼−1 ∆𝑥𝑥𝐼𝐼 + 𝐷𝐷𝑔𝑔,𝐼𝐼 ∆𝑥𝑥𝐼𝐼−1 𝑔𝑔,𝐼𝐼(145)となる。CMR 法の場合と同様に、(144)、(145)式と三次元拡散方程式である(35)～(39)式と比較
すると、計算式はよく似ていることが分かる。そのため、CMFD 法を用いた三次元拡散方程式
は、(35)、(144)、(145)式より
𝑦𝑦−

𝑥𝑥−
0
𝑧𝑧−
𝐴𝐴𝑔𝑔,𝐼𝐼,𝐽𝐽,𝐾𝐾
𝜙𝜙𝑔𝑔,𝐼𝐼−1,𝐽𝐽,𝐾𝐾 + 𝐴𝐴𝑔𝑔,𝐼𝐼,𝐽𝐽,𝐾𝐾 𝜙𝜙𝑔𝑔,𝐼𝐼,𝐽𝐽−1,𝐾𝐾 + 𝐴𝐴𝑔𝑔,𝐼𝐼,𝐽𝐽,𝐾𝐾
𝜙𝜙𝑔𝑔,𝐼𝐼,𝐽𝐽,𝐾𝐾−1 + 𝐴𝐴𝑔𝑔,𝐼𝐼,𝐽𝐽,𝐾𝐾
𝜙𝜙𝑔𝑔,𝐼𝐼,𝐽𝐽,𝐾𝐾

𝑥𝑥±
𝐴𝐴𝑔𝑔,𝐼𝐼,𝐽𝐽,𝐾𝐾=−
𝑦𝑦±

𝐴𝐴𝑔𝑔,𝐼𝐼,𝐽𝐽,𝐾𝐾 = −𝑧𝑧±
𝐴𝐴𝑔𝑔,𝐼𝐼,𝐽𝐽,𝐾𝐾
=−

𝑦𝑦+

𝑥𝑥+
𝑧𝑧+
+ 𝐴𝐴𝑔𝑔,𝐼𝐼,𝐽𝐽,𝐾𝐾
𝜙𝜙𝑔𝑔,𝐼𝐼+1,𝐽𝐽,𝐾𝐾 + 𝐴𝐴𝑔𝑔,𝐼𝐼,𝐽𝐽,𝐾𝐾 𝜙𝜙𝑔𝑔,𝐼𝐼,𝐽𝐽+1,𝐾𝐾 + 𝐴𝐴𝑔𝑔,𝐼𝐼,𝐽𝐽,𝐾𝐾
𝜙𝜙𝑔𝑔,𝐼𝐼,𝐽𝐽,𝐾𝐾+1
𝜒𝜒𝑔𝑔,𝐼𝐼,𝐽𝐽,𝐾𝐾=� 𝜈𝜈Σ𝑓𝑓,𝑔𝑔′ ,𝐼𝐼,𝐽𝐽,𝐾𝐾 𝜙𝜙𝑔𝑔′ ,𝐼𝐼,𝐽𝐽,𝐾𝐾 + � Σ𝑠𝑠,𝑔𝑔′ →𝑔𝑔,𝐼𝐼,𝐽𝐽,𝐾𝐾 𝜙𝜙𝑔𝑔′ ,𝐼𝐼,𝐽𝐽,𝐾𝐾
𝑘𝑘
′
′
𝑔𝑔

𝑔𝑔 ≠𝑔𝑔�𝑠𝑠∓,𝑔𝑔,𝐼𝐼±1,𝐽𝐽,𝐾𝐾 𝐷𝐷𝑔𝑔,𝐼𝐼,𝐽𝐽,𝐾𝐾
𝐷𝐷

�𝐷𝐷𝑔𝑔,𝐼𝐼±1,𝐽𝐽,𝐾𝐾 ∆𝑥𝑥𝐼𝐼 + 𝐷𝐷𝑔𝑔,𝐼𝐼,𝐽𝐽,𝐾𝐾 ∆𝑥𝑥𝐼𝐼±1 �∆𝑥𝑥𝐼𝐼
�𝑠𝑠∓,𝑔𝑔,𝐼𝐼,𝐽𝐽±1,𝐾𝐾 𝐷𝐷𝑔𝑔,𝐼𝐼,𝐽𝐽,𝐾𝐾
𝐷𝐷

�𝐷𝐷𝑔𝑔,𝐼𝐼,𝐽𝐽,𝐾𝐾±1 ∆𝑧𝑧𝐾𝐾 + 𝐷𝐷𝑔𝑔,𝐼𝐼,𝐽𝐽,𝐾𝐾 ∆𝑧𝑧𝐾𝐾±1 �∆𝑧𝑧𝐾𝐾
𝑦𝑦−

(149)
𝑦𝑦+

0
𝑥𝑥−
𝑥𝑥+
𝑧𝑧+
𝑧𝑧−
𝐴𝐴𝑔𝑔,𝐼𝐼,𝐽𝐽,𝐾𝐾= Σ𝑟𝑟,𝑔𝑔,𝐼𝐼,𝐽𝐽,𝐾𝐾 − �𝐴𝐴𝑔𝑔,𝐼𝐼,𝐽𝐽,𝐾𝐾+ 𝐴𝐴𝑔𝑔,𝐼𝐼,𝐽𝐽,𝐾𝐾+ 𝐴𝐴𝑔𝑔,𝐼𝐼,𝐽𝐽,𝐾𝐾 + 𝐴𝐴𝑔𝑔,𝐼𝐼,𝐽𝐽,𝐾𝐾 + 𝐴𝐴𝑔𝑔,𝐼𝐼,𝐽𝐽,𝐾𝐾+ 𝐴𝐴𝑔𝑔,𝐼𝐼,𝐽𝐽,𝐾𝐾
�

となる。

(147)

(148)

�𝐷𝐷𝑔𝑔,𝐼𝐼,𝐽𝐽±1,𝐾𝐾 ∆𝑦𝑦𝐽𝐽 + 𝐷𝐷𝑔𝑔,𝐼𝐼,𝐽𝐽,𝐾𝐾 ∆𝑦𝑦𝐽𝐽±1 �∆𝑦𝑦𝐽𝐽
�𝑠𝑠∓,𝑔𝑔,𝐼𝐼,𝐽𝐽,𝐾𝐾±1 𝐷𝐷𝑔𝑔,𝐼𝐼,𝐽𝐽,𝐾𝐾
𝐷𝐷

(146)

(150)

CMFD 法では得られた粗メッシュ拡散計算の中性子束𝜙𝜙𝑔𝑔,𝐼𝐼,𝐽𝐽,𝐾𝐾 を用いて n 回目の中性子束を𝑛𝑛

𝜙𝜙 ′𝑔𝑔,𝑖𝑖,𝑗𝑗,𝑘𝑘 =𝜙𝜙𝑔𝑔,𝐼𝐼,𝐽𝐽,𝐾𝐾 𝑉𝑉𝐼𝐼,𝐽𝐽,𝐾𝐾
𝜙𝜙 𝑛𝑛𝑛𝑛∑𝑖𝑖∈𝐼𝐼,𝑗𝑗∈𝐽𝐽,𝑘𝑘∈𝐾𝐾 𝜙𝜙𝑔𝑔,𝑖𝑖,𝑗𝑗,𝑘𝑘𝑉𝑉𝑖𝑖,𝑗𝑗,𝑘𝑘 𝑔𝑔,𝑖𝑖,𝑗𝑗,𝑘𝑘(151)

で補正している。なお、n 回目の実効増倍率についても(151)式で得られた中性子束を用いて次式のように新たに計算される。
𝑛𝑛

𝑘𝑘 ′ = � � � � 𝜈𝜈Σ𝑓𝑓,𝑔𝑔′ ,𝑖𝑖 𝜙𝜙 ′𝑔𝑔′ ,𝑖𝑖,𝑗𝑗,𝑘𝑘 𝑉𝑉𝑖𝑖,𝑗𝑗,𝑘𝑘
𝑖𝑖

𝑗𝑗

𝑘𝑘

𝑔𝑔′

(152)

CMFD 法の計算手順は次の通りである。

1. (2)式を用いて詳細メッシュ体系の拡散係数𝐷𝐷𝑔𝑔,𝑖𝑖,𝑗𝑗,𝑘𝑘 、断面積Σ𝑔𝑔,𝑖𝑖,𝑗𝑗,𝑘𝑘 、中性子束𝜙𝜙𝑔𝑔,𝑖𝑖,𝑗𝑗,𝑘𝑘 から粗メッシュ体系の拡散係数𝐷𝐷𝑔𝑔,𝐼𝐼,𝐽𝐽,𝐾𝐾 、断面積Σ𝑔𝑔,𝐼𝐼,𝐽𝐽,𝐾𝐾 を計算する。

� を計算する。
2. (138)～(141)式を用いて中性子流の補正係数𝐷𝐷

3. (146)～(150)式を用いて粗メッシュ体系において拡散計算が完全に収束するまで計算する。
4. (151)式を用いて詳細メッシュ体系の中性子束𝜙𝜙𝑔𝑔,𝑖𝑖,𝑗𝑗,𝑘𝑘 を補正する。

5. (152)式を用いて詳細メッシュ体系の実効増倍率𝑘𝑘を補正する。
6. 詳細メッシュ体系の次の外部反復に進む。

� を用いることで、詳細メッシュ体系における中性子の入CMFD 法では中性子流の補正係数𝐷𝐷射に対する放出の関係を保存し、詳細メッシュ体系の計算を再現することが出来る。(22)、(23)36

式の不連続因子と(136)、(137)式を比較すると、計算式に大きな違いはない 44。実際、詳細メッシュ計算を完全に収束させた後に(146)～(150)式を用いて粗メッシュ拡散計算を行うと、不連続因子の場合と同様に詳細メッシュ計算の結果を再現する。
3.4.3 GCMR 法
GCMR 法も CMFD 法と同様に粗メッシュ拡散計算の結果を利用して、詳細メッシュ計算の中性子束及び実効増倍率を修正する手法である。計算式も CMFD 法とほぼ同じで、GCMR 法
では表面中性子流𝐽𝐽𝑠𝑠±,𝑔𝑔,𝐼𝐼 を(136)、(137)式の代わりに次のように書き表す。
𝐽𝐽𝑠𝑠+,𝑔𝑔,𝐼𝐼 = −𝛼𝛼𝑔𝑔,𝐼𝐼 𝜙𝜙𝑠𝑠+,𝑔𝑔,𝐼𝐼 + 𝛽𝛽𝑠𝑠+,𝑔𝑔,𝐼𝐼 𝜙𝜙𝑔𝑔,𝐼𝐼 = −𝛽𝛽𝑠𝑠−,𝑔𝑔,𝐼𝐼+1 𝜙𝜙𝑔𝑔,𝐼𝐼+1 + 𝛼𝛼𝑔𝑔,𝐼𝐼+1 𝜙𝜙𝑠𝑠+,𝑔𝑔,𝐼𝐼
𝐽𝐽𝑠𝑠−,𝑔𝑔,𝐼𝐼 = −𝛼𝛼𝑔𝑔,𝐼𝐼−1 𝜙𝜙𝑠𝑠−,𝑔𝑔,𝐼𝐼 + 𝛽𝛽𝑠𝑠+,𝑔𝑔,𝐼𝐼−1 𝜙𝜙𝑔𝑔,𝐼𝐼−1 = −𝛽𝛽𝑠𝑠−,𝑔𝑔,𝐼𝐼 𝜙𝜙𝑔𝑔,𝐼𝐼 + 𝛼𝛼𝑔𝑔,𝐼𝐼 𝜙𝜙𝑠𝑠−,𝑔𝑔,𝐼𝐼

ここで、𝛼𝛼𝑔𝑔,𝐼𝐼 、𝛽𝛽𝑠𝑠±,𝑔𝑔,𝐼𝐼

(153)
(154)

はそれぞれ GCMR 因子、リバランス因子と呼ばれ 45、GCMR 因子は入力値である。GCMR 法は CMR 法と CMFD 法を統合的に取り扱った手法で、GCMR 因子が𝛼𝛼𝑔𝑔,𝐼𝐼 =12𝛼𝛼𝑔𝑔,𝐼𝐼 =2𝐷𝐷𝑔𝑔,𝐼𝐼
∆𝑥𝑥𝐼𝐼

(155)

の場合、CMR 法と等価となり、GCMR 因子が(156)の場合、CMFD 法と等価となる 4)。このように GCMR 法は CMR 法及び CMFD 法の一般形であり、そのため”Generalized” CMR 法と名付けられた。GCMR 法では CMFD 法と異なり、GCMR
因子を調整することにより、CMFD 法では収束安定性が悪く収束しなかった体系においても収束する場合があり、CMFD 法に比べてさらに収束安定性が向上した手法となっている。
GCMR 法を用いた場合、(144)、(145)式は(153)、(154)式より
𝐽𝐽𝑠𝑠+,𝑔𝑔,𝐼𝐼 =𝐽𝐽𝑠𝑠−,𝑔𝑔,𝐼𝐼 =𝛼𝛼𝑔𝑔,𝐼𝐼+1 𝛽𝛽𝑠𝑠+,𝑔𝑔,𝐼𝐼𝛼𝛼𝑔𝑔,𝐼𝐼 𝛽𝛽𝑠𝑠−,𝑔𝑔,𝐼𝐼+1
𝜙𝜙𝑔𝑔,𝐼𝐼 −𝜙𝜙𝛼𝛼𝑔𝑔,𝐼𝐼 + 𝛼𝛼𝑔𝑔,𝐼𝐼+1
𝛼𝛼𝑔𝑔,𝐼𝐼 + 𝛼𝛼𝑔𝑔,𝐼𝐼+1 𝑔𝑔,𝐼𝐼+1(157)

𝛼𝛼𝑔𝑔,𝐼𝐼 𝛽𝛽𝑠𝑠+,𝑔𝑔,𝐼𝐼−1
𝛼𝛼𝑔𝑔,𝐼𝐼−1 𝛽𝛽𝑠𝑠−,𝑔𝑔,𝐼𝐼𝜙𝜙𝑔𝑔,𝐼𝐼−1 −𝜙𝜙𝛼𝛼𝑔𝑔,𝐼𝐼−1 + 𝛼𝛼𝑔𝑔,𝐼𝐼𝛼𝛼𝑔𝑔,𝐼𝐼−1 + 𝛼𝛼𝑔𝑔,𝐼𝐼 𝑔𝑔,𝐼𝐼(158)となる。また、(146)～(150)式と同様に、GCMR 法を用いた三次元拡散方程式は、(35)、(157)、(158)式より
𝑦𝑦−

𝑥𝑥−
0
𝑧𝑧−
𝐴𝐴𝑔𝑔,𝐼𝐼,𝐽𝐽,𝐾𝐾
𝜙𝜙𝑔𝑔,𝐼𝐼−1,𝐽𝐽,𝐾𝐾 + 𝐴𝐴𝑔𝑔,𝐼𝐼,𝐽𝐽,𝐾𝐾 𝜙𝜙𝑔𝑔,𝐼𝐼,𝐽𝐽−1,𝐾𝐾 + 𝐴𝐴𝑔𝑔,𝐼𝐼,𝐽𝐽,𝐾𝐾
𝜙𝜙𝑔𝑔,𝐼𝐼,𝐽𝐽,𝐾𝐾−1 + 𝐴𝐴𝑔𝑔,𝐼𝐼,𝐽𝐽,𝐾𝐾
𝜙𝜙𝑔𝑔,𝐼𝐼,𝐽𝐽,𝐾𝐾

𝑥𝑥±
𝐴𝐴𝑔𝑔,𝐼𝐼,𝐽𝐽,𝐾𝐾
=−
44

𝑦𝑦+

𝑥𝑥+
𝑧𝑧+
+ 𝐴𝐴𝑔𝑔,𝐼𝐼,𝐽𝐽,𝐾𝐾
𝜙𝜙𝑔𝑔,𝐼𝐼+1,𝐽𝐽,𝐾𝐾 + 𝐴𝐴𝑔𝑔,𝐼𝐼,𝐽𝐽,𝐾𝐾 𝜙𝜙𝑔𝑔,𝐼𝐼,𝐽𝐽+1,𝐾𝐾 + 𝐴𝐴𝑔𝑔,𝐼𝐼,𝐽𝐽,𝐾𝐾
𝜙𝜙𝑔𝑔,𝐼𝐼,𝐽𝐽,𝐾𝐾+1
𝜒𝜒𝑔𝑔,𝐼𝐼,𝐽𝐽,𝐾𝐾=� 𝜈𝜈Σ𝑓𝑓,𝑔𝑔′ ,𝐼𝐼,𝐽𝐽,𝐾𝐾 𝜙𝜙𝑔𝑔′ ,𝐼𝐼,𝐽𝐽,𝐾𝐾 + � Σ𝑠𝑠,𝑔𝑔′ →𝑔𝑔,𝐼𝐼,𝐽𝐽,𝐾𝐾 𝜙𝜙𝑔𝑔′ ,𝐼𝐼,𝐽𝐽,𝐾𝐾
𝑘𝑘
′
′
𝑔𝑔

𝛽𝛽𝑠𝑠∓,𝑔𝑔,𝐼𝐼±1,𝐽𝐽,𝐾𝐾 𝛼𝛼𝑔𝑔,𝐼𝐼,𝐽𝐽,𝐾𝐾

�𝛼𝛼𝑔𝑔,𝐼𝐼±1,𝐽𝐽,𝐾𝐾 + 𝛼𝛼𝑔𝑔,𝐼𝐼,𝐽𝐽,𝐾𝐾 �∆𝑥𝑥𝐼𝐼(159)

𝑔𝑔 ≠𝑔𝑔(160)

計算式は似ているが、CMFD 法では(138)～(141)式のように非均質体系の中性子流を再現するように均質体系の中性子流を補正しており、発想は異なっている。
45 GCMR 法でのリバランス因子𝛽𝛽
𝑠𝑠±,𝑔𝑔,𝐼𝐼 と CMR 法でのリバランス因子𝑓𝑓𝑔𝑔,𝑖𝑖 は名前は同じである
が、違う因子であることに注意する必要がある。
37

𝑦𝑦±

𝐴𝐴𝑔𝑔,𝐼𝐼,𝐽𝐽,𝐾𝐾 = −𝑧𝑧±
𝐴𝐴𝑔𝑔,𝐼𝐼,𝐽𝐽,𝐾𝐾=−

𝛽𝛽𝑠𝑠∓,𝑔𝑔,𝐼𝐼,𝐽𝐽±1,𝐾𝐾 𝛼𝛼𝑔𝑔,𝐼𝐼,𝐽𝐽,𝐾𝐾(161)

�𝛼𝛼𝑔𝑔,𝐼𝐼,𝐽𝐽±1,𝐾𝐾 + 𝛼𝛼𝑔𝑔,𝐼𝐼,𝐽𝐽,𝐾𝐾 �∆𝑦𝑦𝐽𝐽
𝛽𝛽𝑠𝑠∓,𝑔𝑔,𝐼𝐼,𝐽𝐽,𝐾𝐾±1 𝛼𝛼𝑔𝑔,𝐼𝐼,𝐽𝐽,𝐾𝐾

�𝛼𝛼𝑔𝑔,𝐼𝐼,𝐽𝐽,𝐾𝐾±1 + 𝛼𝛼𝑔𝑔,𝐼𝐼,𝐽𝐽,𝐾𝐾 �∆𝑍𝑍𝐾𝐾

(162)
𝑦𝑦−

𝑦𝑦+

0
𝑥𝑥−
𝑥𝑥+
𝑧𝑧+
𝑧𝑧−
𝐴𝐴𝑔𝑔,𝐼𝐼,𝐽𝐽,𝐾𝐾= Σ𝑟𝑟,𝑔𝑔,𝐼𝐼,𝐽𝐽,𝐾𝐾 − �𝐴𝐴𝑔𝑔,𝐼𝐼,𝐽𝐽,𝐾𝐾+ 𝐴𝐴𝑔𝑔,𝐼𝐼,𝐽𝐽,𝐾𝐾+ 𝐴𝐴𝑔𝑔,𝐼𝐼,𝐽𝐽,𝐾𝐾 + 𝐴𝐴𝑔𝑔,𝐼𝐼,𝐽𝐽,𝐾𝐾 + 𝐴𝐴𝑔𝑔,𝐼𝐼,𝐽𝐽,𝐾𝐾+ 𝐴𝐴𝑔𝑔,𝐼𝐼,𝐽𝐽,𝐾𝐾
�

となる。

(163)

3.4.4 粗メッシュ拡散加速法を適用する場合の注意点
SOR 法や Chebyshev 外挿法と異なり、粗メッシュ拡散計算は外部反復のたびに実施すると、収束性が悪化してしまうことが多い。CMFD 法や GCMR 法では、n 回目の外部反復における中性子束と実効増倍率が最終的に到達すると期待される値を粗メッシュ拡散計算を用いて効率的に計算していることになる。しかし、収束途中では粗メッシュ拡散計算の結果は詳細メッシ
ュ計算の解を再現している訳ではないため、粗メッシュ拡散計算で補正した中性子束は詳細メッシュ計算の解とは異なった値になっている。そのため、粗メッシュ拡散計算を外部反復のたびに実施していると、SOR 法で緩和係数 ω に大きな値を入れた時のように粗メッシュ拡散係数による補正が強すぎて、計算結果が振動したり、場合によっては発散してしまったりする場合がある。そのような場合には、粗メッシュ拡散計算を実施した後は詳細メッシュ計算のみの
外部反復計算を 2～5 回繰り返す必要がある 46。なお、詳細メッシュ計算を行う前に、粗メッシュ拡散計算を行い、粗メッシュ拡散計算で得られた中性子束を初期値とすることで、中性子束の初期値を全て 1.0 とする場合や、核分裂断面積で重みを付けた値とする場合に比べ、実際の中性子束に近い中性子束から計算を始めることができるため、収束性の向上が期待できる。
粗メッシュ拡散計算では、メッシュだけでなく、エネルギー群数も縮約する場合がある。例
えば 3.3.3 節で述べたように、エネルギー群数を一群とすることで、Wielandt 法を用いて効率的な計算が可能となる。もちろん、エネルギー群によって補正値を変えることが出来なくなるので、エネルギー群数を縮約した方が詳細メッシュ体系の外部反復が増加する傾向にあるが、ほとんど変わらないことが多い。また、粗メッシュ詳細群の拡散計算に対し、粗メッシュ一群の拡散計算を用いて加速し、粗メッシュ拡散計算を高速化するという二段階加速を行う場合もある 47。
粗メッシュ拡散計算の収束条件の方が、詳細メッシュ計算の収束条件に比べて収束条件を厳しく設定しなければ、適切な計算が得られないことが多いので注意する必要がある。これは粗メッシュ拡散計算で収束条件を満たしていても、詳細メッシュ計算で収束条件を満たしていな
い場合、粗メッシュ拡散計算の収束条件の範囲内で計算結果が振動してしまい、詳細メッシュ
46

詳細メッシュ計算のみの外部反復の回数を増やした方が収束は安定するが、計算時間が延びてしまうので、最初は 2～5 回程度としておき、収束しない場合は回数を増やすことが望ましい。また、体系によっては外部反復のたびに粗メッシュ拡散計算を実施しても収束する場合がある。
47 空間だけで二段階加速を行う場合もある。例えば、非均質の全炉心計算を行う場合には、一段階目で pin-by-pin の均質拡散計算を行い、その pin-by-pin の均質拡散計算を加速するため、二段階目で集合体均質の拡散計算を行うことがある。
38

計算が収束しない可能性があるからである。そのため、粗メッシュ拡散計算の収束条件は、詳
細メッシュ計算の収束条件の 100 分の 1 程度に設定することが望ましい。
粗メッシュ拡散計算において、計算が発散したり、解が振動したりする場合がある。そのような粗メッシュ拡散計算の結果を詳細メッシュ計算の結果に反映させた場合、間違った補正を行うことになり、収束性が悪化してしまう。このように完全に収束していない粗メッシュ拡散計算を詳細メッシュ計算に反映させることを防ぎ、収束安定性を向上させるため、計算コード
によっては粗メッシュ拡散計算がきちんと収束したかどうかを判定する場合がある。代表的な例が粗メッシュ拡散計算の外部反復回数に上限値を設定することである。この方法では、粗メ
ッシュ拡散計算の外部反復回数がその上限値に達した場合に粗メッシュ拡散計算が発散したと判定し、粗メッシュ拡散計算の結果を詳細メッシュ計算の補正に利用しないようにスキップしている。このような対策を取っている場合には、粗メッシュ拡散計算が外部反復回数の上限値以下で収束するように、粗メッシュ拡散計算の外部反復回数の上限値を十分に大きな値に設定する必要がある。
GCMR 法を用いた計算コードでは、GCMR 因子を
𝛼𝛼𝑔𝑔,𝐼𝐼 = 𝑓𝑓

2𝐷𝐷𝑔𝑔,𝐼𝐼
∆𝑥𝑥𝐼𝐼

(164)

GCMR
とし、
f を入力パラメータとして GCMR 因子を調整できるようにしている場合が多い 48。
因子は通常 CMFD 法と同様に f を 1.0 とすると収束が早いが、計算体系によっては f を 1.0 とした場合では計算が収束しないことがある。そのような体系では、f を 2.0～5.0 程度の値に設定すると収束する場合がある。なお、f を 10 といった非常に大きな値にしても計算が収束する
場合があるが、そのような非常に大きな f を用いた場合には、適切な計算結果となっていない
場合があり、計算結果が正しいかどうかを慎重に検証する必要がある。
3.2.1.5 節の SOR 法の説明において、解に対して振動している場合に有効な手法として後退加速について説明した。後退加速は、粗メッシュ拡散加速法にも有効な手法である。計算体系
によっては粗メッシュ拡散加速法を適用すると、解が振動する場合がある。これは粗メッシュ
拡散計算による補正が大きすぎることが原因である。そのような場合にはリバランス因子
� )を𝛽𝛽(CMFD 法の場合は中性子流の補正係数𝐷𝐷𝛽𝛽 ′ = (1.0 − 𝜔𝜔)𝛽𝛽old + ω𝛽𝛽new(165)と補正し、𝜔𝜔に 1.0 よりも小さい値を設定し、後退加速を行うと収束するようになる場合がある。ここで、𝛽𝛽oldは前回のリバランス因子であり、𝛽𝛽newは今回のリバランス因子である。

Appendix A. 修正一群計算について
修正一群計算について本文の 3.2.1.1 節の脚注 27 で紹介したが、あまり参考となる公開文献が見つからなかった。修正一群計算は従来の BWR 炉心計算コードで広く利用されてきた手法であり、高速化手法の一つであることから、簡単に説明する。

� )を変える方法 33)や、リバ
他にはリバランス因子𝛽𝛽(CMFD 法の場合は中性子流の補正係数𝐷𝐷
ランス因子にメッシュ幅による補正を行う partial-current CMFD (p-CMFD)法 34)や Artificial Grid
Diffusion (AGD)法 35)など、様々な補正方法が提案されている。

48

39

修正一群計算とは、計算時間を増加させることなく計算精度を向上させる手法である。メッ
シュ i の二群拡散方程式は次式で書き表すことができる。
−𝐷𝐷1,𝑖𝑖−𝐷𝐷2,𝑖𝑖d2
𝜒𝜒1
𝜙𝜙 + Σ𝑡𝑡,1,𝑖𝑖 𝜙𝜙1,𝑖𝑖 = Σ𝑠𝑠,1→1,𝑖𝑖 𝜙𝜙1,𝑖𝑖 + Σ𝑠𝑠,2→1,𝑖𝑖 𝜙𝜙2,𝑖𝑖 +�νΣ𝑓𝑓,1,𝑖𝑖 𝜙𝜙1,𝑖𝑖 + νΣ𝑓𝑓,2,𝑖𝑖 𝜙𝜙2,𝑖𝑖 �
d𝑥𝑥 2 1,𝑖𝑖
𝑘𝑘𝑒𝑒𝑒𝑒𝑒𝑒d2
𝜒𝜒2
𝜙𝜙 + Σ𝑡𝑡,2,𝑖𝑖 𝜙𝜙2,𝑖𝑖 = Σ𝑠𝑠,1→2,𝑖𝑖 𝜙𝜙1,𝑖𝑖 + Σ𝑠𝑠,2→2,𝑖𝑖 𝜙𝜙2,𝑖𝑖 +�νΣ𝑓𝑓,1,𝑖𝑖 𝜙𝜙1,𝑖𝑖 + νΣ𝑓𝑓,2,𝑖𝑖 𝜙𝜙2,𝑖𝑖 �
d𝑥𝑥 2 2,𝑖𝑖
𝑘𝑘𝑒𝑒𝑒𝑒𝑒𝑒

(A.1)

(A.2)

ここで、(A.1)、(A.2)式に対し、次の 3 つの仮定を適用する。

① 核分裂で発生する中性子は高速群(1 群)でのみ発生する。
𝜒𝜒1 ≈ 1、𝜒𝜒2 ≈ 0(A.3)② 熱群(2 群)から高速群への上方散乱がない。
Σ𝑠𝑠,2→1,𝑖𝑖 ≈ 0(A.4)③ 熱群の中性子の拡散は吸収反応に比べて極めて小さい 49。
Σ𝑡𝑡,2,𝑖𝑖 𝜙𝜙2,𝑖𝑖 ≫ −D2,𝑖𝑖

d2
𝜙𝜙 ≈ 0d𝑥𝑥 2 2,𝑖𝑖

(A.5)

これらの仮定を適用すると、(A.1)、(A.2)式は次のように書き表せられる。
−𝐷𝐷1,𝑖𝑖d2
𝜒𝜒1
𝜙𝜙1,𝑖𝑖 + Σ𝑡𝑡,1,𝑖𝑖 𝜙𝜙1,𝑖𝑖 = Σ𝑠𝑠,1→1,𝑖𝑖 𝜙𝜙1,𝑖𝑖 +�νΣ𝑓𝑓,1,𝑖𝑖 𝜙𝜙1,𝑖𝑖 + νΣ𝑓𝑓,2,𝑖𝑖 𝜙𝜙2,𝑖𝑖 �2d𝑥𝑥𝑘𝑘𝑒𝑒𝑒𝑒𝑒𝑒Σ𝑡𝑡,2,𝑖𝑖 𝜙𝜙2,𝑖𝑖 = Σ𝑠𝑠,1→2,𝑖𝑖 𝜙𝜙1,𝑖𝑖 + Σ𝑠𝑠,2→2,𝑖𝑖 𝜙𝜙2,𝑖𝑖

(A.6)
(A.7)

ここで、(A.7)式の熱群の拡散方程式から、高速群と熱群の中性子束の関係が得られる。
𝜙𝜙2,𝑖𝑖 =Σ𝑠𝑠,1→2,𝑖𝑖
Σ𝑠𝑠,1→2,𝑖𝑖
𝜙𝜙1,𝑖𝑖 =𝜙𝜙Σ𝑡𝑡,2,𝑖𝑖 − Σ𝑠𝑠,2→2,𝑖𝑖
Σ𝑎𝑎,2,𝑖𝑖 1,𝑖𝑖

(A.8)

(A.8)式を(A.6)式の高速群の拡散方程式に代入すると、修正一群拡散方程式が得られる。
−𝐷𝐷1,𝑖𝑖d2
𝜒𝜒1
Σ𝑠𝑠,1→2,𝑖𝑖
𝜙𝜙 + Σ𝑡𝑡,1,𝑖𝑖 𝜙𝜙1,𝑖𝑖 = Σ𝑠𝑠,1→1,𝑖𝑖 𝜙𝜙1,𝑖𝑖 +�νΣ𝑓𝑓,1,𝑖𝑖 + νΣ𝑓𝑓,2,𝑖𝑖
� 𝜙𝜙1,𝑖𝑖
d𝑥𝑥 2 1,𝑖𝑖
𝑘𝑘𝑒𝑒𝑒𝑒𝑒𝑒Σ𝑎𝑎,2,𝑖𝑖(A.9)

修正一群拡散方程式は(A.3)～(A.5)式の仮定が正しければ適切な結果が得られる。そこで(A.3)～(A.5)式の仮定について考えてみる。
軽水炉での主な核分裂性核種である 235U と 239Pu の核分裂スペクトルを図 A.1 に示す。図 A.1に示すように、ほとんどの中性子が 10eV 以上で発生している。一般的な 2 群計算では、熱群と高速群の境界を 0.625eV としていることから、(A.3)式の仮定は妥当である。
原子核の運動エネルギーはボルツマン定数 kB と温度 T の積 kBT で表せられる。室温(20 度、
3

293.15K)での原子核の運動エネルギーは最確値k B 𝑇𝑇 = 0.0253eV、平均値2 k B 𝑇𝑇 = 0.0380eV である。軽水炉では燃料棒内でも最大で 900K 程度なので、原子核の運動エネルギーは熱群と高速群の境界である 0.625eV に比べて十分に小さいことから、高速群への上方散乱は起きにくい。
このことから(A.4)式の仮定は妥当である。

49

熱群の中性子の拡散を考慮するため、熱群の拡散をバックリングとして取り扱うこともある

20)。

40

中性子の拡散は炉内の中性子束分布に依存することから、原子炉の幾何形状やサイズによって変化する。発電用軽水炉のような大型の熱中性炉では熱中性子束は平坦であり、(A.5)式の仮定は妥当である 50。
以上のことから、発電用軽水炉の炉心解析手法に適用する上では、(A.3)～(A.5)式の仮定は妥当である。
修正一群方程式では、計算式上は一群の拡散方程式となるが、実質的には二群相当の計算を行っていることになる。エネルギー群数が計算時間とメモリー量に与える影響は大きいことから、修正一群計算は計算時間とメモリー量の節約に効果的な手法であり、計算機性能の低かった 70～80 年代に開発された BWR の炉心解析コードで主に採用されてきた。
4.0E-07
3.5E-07

確率 [-]

3.0E-07
2.5E-07

235U
235U
239Pu
239Pu

2.0E-07
1.5E-07
1.0E-07
5.0E-08
0.0E+00
1.0E-05

1.0E-03

1.0E-01

1.0E+01

1.0E+03

1.0E+05

1.0E+07

核分裂で発生する即発中性子のエネルギー [eV]
図 A.1 235U と 239Pu の核分裂スペクトル

謝辞
本資料を作成するにあたり、日本原子力研究開発機構で炉物理研究に関わる多くの研究者・技術者の助力を得た。特に原子力基礎工学研究センター炉物理標準コード研究グループの森貴
正氏、石川眞氏、原子力基礎工学研究センター核変換システム研究グループの方野量太氏、高速炉基盤技術開発部炉心安全・融体評価技術開発グループの加藤慎也氏、安全研究センター臨界安全研究グループの渡邉友章氏には本資料について多くの助言やご指摘を頂いた。ここに深く感謝の意を表する。

この仮定は燃料集合体単位で均質化する従来の炉心計算では正しいが、pin-by-pin 計算では成り立たないことに注意する必要がある。pin-by-pin 計算では集合体内の複雑な濃縮度分布やGd 入り燃料棒、そして BWR ではウォーターロッドやウォーターギャップを考慮しなければならない。炉心全体の中性子束は平坦であるが、集合体内の中性子束分布を見ると、高濃縮度燃料棒や Gd 入り燃料棒付近などで中性子束の変化が大きいため、(A.5)式の仮定は成り立たな
い。

50

41

参考文献
1) J. R. Askew, "A Characteristics Formulation to the Neutron Transport Equation in ComplicatedGeometries," AEEW-M1108, UKAEA, Winfrith (1972).
2) M. Tatsumi, A. Gihou, T. Takeda, "Space and Angular Dependence of Interface Currents in theMultiband-CCCP Method," J. Nucl. Sci. Technol., 37, pp.572-580 (2000).
3) B. G. Carlson, "Solution of the Transport Equation by Sn Approximations," LA-1599, Los Alamos
Scientific Laboratory (1953).
4) A. Yamamoto, "Generalized Coarse-Mesh Rebalance Method for Acceleration of Neutron Transport
Calculations," Nucl. Sci. Eng., 115, pp.274-282 (2005).
5) T. Takeda, “Neutronics Codes Currently Used in Japan for Fast and Thermal Reactor Applications,”

Proc. Physor2004, Chicago, USA, Apr. 25-29, 2004 (2004). [CD-ROM]
6) 小坂進矢、巽雅洋、軽水炉商用炉心の燃焼特性とその解析、第 38 回炉物理夏期セミナー(2006)
7) 山本章夫、軽水炉と高速炉の比較～原子炉物理の観点から～、第 42 回炉物理夏期セミナー(2010)
8) A. Hebert, "A consistent technique for the pin-by-pin homogenization of a pressurized water reactor
assembly," Nucl. Sci. Eng., 113, pp.227-238 (1993).
9) K. S. Smith, "Assembly homogenization techniques for light water reactor analysis," Prog. Nucl.
Energy, 17, pp. 303-335 (1986).
10) 奥村啓介、近代ノード法と不連続因子の基礎、第 36 回炉物理夏期セミナー (2004).
11) 奥村啓介、
MOSRA-Light；ベクトル計算機のための高速 3 次元中性子拡散ノード法コード、
JAERI-Data/Code 98-025 (1998)12) E. W. Larsen, J. E. Morel, John M. McGhee, "Asymptotic Derivation of the Multigroup P1 andSimplified PN Equations with Anisotropic Scattering," Nucl. Sci. Eng., 123, pp.328-342 (1996).
13) 山本章夫、拡散方程式の数値解法の基礎、第 36 回炉物理夏期セミナー (2004).
14) T. Sakamoto, T. Endo, A. Yamamoto, “Applicability of angular flux discontinuity factor preserving
region-wise leakage for integro-differential transport equation,” J. Nucl. Sci. Technol., 51,
pp.1264-1273 (2014).
15) A. Yamamoto, T. Sakamoto, T. Endo, “Discontinuity factors for simplified P3 theory,” Nucl. Sci.
Eng., 183, pp.39-51 (2017).
16) A. Yamamoto, M. Tatsumi, Y. Kitamura, Y. Yamane, “Improvement of the SPH method for
Pin-by-Pin Core Calculations,” J. Nucl. Sci. Technol., 41, pp.1155-1165 (2004).
17) 牛尾直史、Characteristics 法の基礎、第 36 回炉物理夏期セミナー (2004).
18) Y. A. Chao, "A theoretical analysis of the coarse mesh finite difference representation in advanced
nodal methods," Proc. M&C’99, Madrid, Spain, Sept. 27-30, 1999, 1, pp.117-126 (1999).
19) 水島次郎、柳瀬眞一郎、理工学のための数値計算法、数理工学社 (2002).
20) General Electric, “Steady-State Nuclear Methods,” Licensing Topical Report, NEDO-30130-A(1985).
4221) E. L. Wachpress, “Trail of a Lyapunov equation solver,” Comput. Math. Appl., 55, pp.1653-1659(2008).
22) 戸川隼人、共役勾配法(シリーズ新しい応用の数学 17)、教育出版 (1977).
23) Y. Saad, Iterative methods for sparse linear systems 2nd ed., SIAM, 2003.
24) 張紹良、大規模連立 1 次方程式の解法－クリロフ部分空間法、応用数理、8[4] (1998).
25) H. E. Wrigley, “Accelerating the Jacobi method for solving simultaneous equations by Chebyshev
extrapolation when the eigenvalues of the iteration matrix are complex,” Comput. J., 6, pp169-176,(1964).
26) P. J. Turinsky, R. M. K. Al-Chalabi, P. Engrand, H. N. Sarsour, F. X. Faure, NESTLE Version 5.2.1:
Few-group neutron diffusion equation solver utilizing the nodal expansion method for eigenvalue,adjoint, fixed-source steady-state and transient problems, Electric Power Research Center, North
Carolina State University (2003).
27) T. M. Sutton, “Wielandt Iteration as Applied to the Nodal Expansion Method,” Nucl. Sci. Eng., 98,
pp.169-173 (1988).
28) E. E. Lewis and W. F. Miller, Jr., Computational Methods of Neutron Transport, John Wiley & Sons,
New York (1984).
29) N. Z. Cho and C. J. Park, “A Comparison of Coarse Mesh Rebalance (CMR) and Coarse Mesh
Finite Difference (CMFD) Acceleration Methods for the Neutron Transport Calculations,”
NUPART-2002-02, KAIST (2002).
30) D. Knott and A. Yamamoto, “Lattice Physics Computations,” Handbook of Nuclear Engineering, 9,
pp. 1155-1158, Springer (2010).
31) K. S. Smith, “Nodal method storage reduction by non-linear iteration,” Trans. Am. Nucl. Soc., 44 p.
265 (1983).
32) E. M. Gelbard and A. Hageman, “The Synthetic Method as Applied to the Sn Equations,” Nucl. Sci.
Eng., 37, pp.288-298 (1969).
33) E. W. Larsen, “Infinite-Medium Solutions to the Transport Equation, SN Discretization Schemes,
and the Diffusion Approximation,” Transport Theory Statist. Phys., 32, pp.623-643 (2003).
34) N. Z. Cho, G. S. Lee, C. J. Park, “Partial Current-Based CMFD Acceleration of the 2D/1D fusion
method for 3D whole-Core Transport Calculations,” Trans. Am. Nucl. Soc., 88, 594 (2003).
35) M. Jarrett, B. Kochunas, A. Zhu, T. Downar, “Analysis of stabilization techniques for CMFD
acceleration of neutron transport problems,” Nucl. Sci. Eng., 184, pp.208-227 (2016).

43