炉物理の研究

第 67 号（2015 年 3 月）

＜原子力学会賞受賞記念寄稿＞

Characteristics 法における内部反復の収束性に関する検討(株)原子力エンジニアリング
解析サービス本部

解析グループ
田渕将人

1. はじめに
この度は 2013 年度原子力学会論文賞を受賞する機会に恵まれ、大変光栄に感じております。受賞対象となった研究は、characteristics 法（以降 MOC）の収束性に関するトピックです[1]。MOC は決定論的な中性子輸送計算手法であり、計算精度や計算コストの面から実用性の高い強力な計算手法の 1 つです。私が所属する(株)原子力エンジニアリングにおいても、次世代の燃料集合体核計算コードとして MOC に基づく計算コード AEGIS の開発を続けてきました[2]。AEGIS の開発では、名古屋大学、原子燃料工業(株)の協力の下、燃料集合体核計算に関する数々の新手法を開発してきましたが、今回の受賞対象となった研究成果もそのうちの 1 つです。それでは、本稿にて本研究の概要を紹介していきたいと思います。
2. 研究の概要
2.1 研究の背景
本研究を始めるきっかけとなったのは、AEGIS での臨界実験解析における不具合でした。
正しく計算を実行しているはずなのに計算が全然収束しません。収束しないどころか発散します。入力を何回も見直し、細かくバグのチェックもしましたが何も出てきませんでした。そこで、計算条件をいろいろ変えて計算を試みたところ、以下のことがわかりました。
・臨界実験の体系から水反射体を除去した体系の計算であれば収束・輸送補正[3]を適用しなければ収束
このような状況から言えるのは、「水を多く含む体系の計算において輸送補正を適用した場合に収束性悪化が見られる」ということです。ここで、唐突に出てきた“輸送補正”について簡単に補足します。輸送補正とは、全断面積の代わりに輸送断面積を用いることにより非等方散乱の効果を考慮する手法のことです。等方散乱を仮定した場合の計算手順において断面積を置き換えるのみで非等方散乱の効果を考慮できるため、非常に有用な手法です。輸送断面積にはいくつかの定義が存在しますが、AEGIS では(1)式に示す out-scatter近似を用いています。また、断面積のバランスを維持するため(2)式も(1)式と同時に適用します。

炉物理の研究第 67 号（2015 年 3 月）

Σ tr ,i , g = Σ t ,i , g − ∑ Σ s1,i , g → g ′(1)

Σ ′s ,i , g → g = Σ s ,i , g → g − ∑ Σ s1,i , g → g ′

(2)

g′g′

ただし、

Σ tr ,i , g：領域 i, エネルギー第 g 群の輸送断面積

Σ t ,i , g：領域 i, エネルギー第 g 群の全断面積

Σ s1,i , g → g ′：領域 i, エネルギー第 g 群から第 g’群への P1 散乱断面積

Σ s ,i , g → g：領域 i, エネルギー第 g 群の輸送補正適用前の自群散乱断面積

Σ′s ,i , g →g：領域 i, エネルギー第 g 群の輸送補正適用後の自群散乱断面積話を臨界実験解析における収束性に戻すと、どうやら輸送補正の適用により収束性が悪化しているようです。そこで、その原因や対処方法について次節以降で考えていきますが、その前に本研究における前提条件について触れておきます。MOC における反復計算には、外部反復（実効増倍率と中性子源を更新する反復）と内部反復（与えられた中性子源に対する中性子束を求める反復）がありますが、本研究では内部反復に着目することとします。
また、MOC では角度中性子束を取り扱うのですが、角度中性子束の収束性を考えようとすると、非常に問題が複雑になるため、本研究では領域平均中性子束の収束性について考えることにします。
2.2 輸送補正適用による収束性悪化の原因
まずは簡単な問題を例に内部反復の収束性について考えます。1 群 1 領域の固定源問題を想定すると、中性子のバランス方程式より、領域平均中性子束に関する漸化式は(3)式のとおりとなります。また、(3)式は(4)式のように書き直せます。なお、1 群 1 領域の問題を想定しているため、(3), (4)式では領域及びエネルギー群のインデックスを省略しています。

Σ tr φ ( k +1) = Σ sφ ( k ) + Qφ ( k +1) = cφ ( k ) +Q
Σ tr

ただし、φ (k ) ：k 回目の内部反復における領域平均中性子束
Q ：固定中性子源
c：散乱比（ c ≡Σs）
Σ tr

(3)
(4)

炉物理の研究第 67 号（2015 年 3 月）

(4)式について考えると、c が 1 以上の場合、中性子束は反復とともに際限なく大きくなることがわかります（つまり発散します）。これは c が‐1 以下の場合も同様であり、(4)式において中性子束の収束値を得るためには c の絶対値が 1 未満である必要があります。
この事実を踏まえたうえで多領域の問題について考えたいと思いますが、それに際して衝突確率を使って多領域への展開を考えます。これは、等方散乱を仮定すれば MOC においても衝突確率相当の概念を用いることができ、また衝突確率を用いれば領域平均中性子束について表現しやすくなるからです。それでは、(4)式に相当する式を多領域問題に展開してみます。導出の過程は省略しますが、結果として(5)式が得られます。

 φ1( k +1)   P11 L Pn1  c1 0 0  φ1( k )   r1     M  =  M O M  0 O 0  M  +  M 
 ( k +1)  
 ( k )   
 φn   Pn1 L Pnn  0 0 cn  φn   rm (5)ただし、φi(k ) ：領域 i の k 回目の内部反復における領域平均中性子束
Pij ：領域 i から領域 j への衝突確率ci：領域 i の散乱比ri：領域 i における収束性に関係のない成分(5)式において(4)式中の c に相当する部分は、(6)式に示す行列 M となりますが、これは“反復行列”と呼ばれる行列で、収束性の良し悪しを左右する重要な行列です。

 P11 L P1n  c1 0 0 
M =  M O M  0 O 0 
 P L P  0 0 c nn n
 n1(6)

1 領域問題では、内部反復が収束するためには c の絶対値が 1 未満である必要がありまし
たが、多領域問題でも同じようなことが言えます。導出については省略しますが、多領域問題では(6)式に示す反復行列のスペクトル半径（最大固有値の絶対値）が 1 未満であれば収束解を得ることができます。ここでの”最大固有値”とは「絶対値が最大の固有値」を指しています。少し複雑ですので例を挙げますと、M の固有値が (-1.2, -0.4, 0.5) である場合、最大固有値は-1.2 でスペクトル半径は 1.2 となります（このケースではスペクトル半径が 1以上ですので内部反復は収束しません）。
それでは、内部反復が発散する条件におけるスペクトル半径の値を見てみます。内部反復の収束性が悪くなる条件として、「水を多く含む体系」、「輸送補正を適用」がありますので、仮想的にセルピッチを 2cm に設定して通常より水領域を広くした単一燃料棒セル体系

炉物理の研究第 67 号（2015 年 3 月）を想定し、その計算において輸送補正を適用しました。スペクトル半径の評価は、(6)式中の衝突確率を実際に求めた後、反復行列の最大固有値を求めることにより行いました。
AEGIS で採用している 172 群構造[4]における各エネルギー群のスペクトル半径を図 1 に示します。参考のため、同体系で輸送補正を適用しない場合のものも併せて示します。

1.5
輸送補正なしスペクトル半径輸送補正あり1.00.50.0050100
エネルギー群150図 1 セルピッチ 2cm の問題におけるスペクトル半径
図 1 より、輸送補正の適用により明らかにスペクトル半径が大きくなることが確認できます。また図は省略しますが、輸送補正を適用したとしてもセルピッチを小さくすれば（水領域の割合を小さくすれば）スペクトル半径は小さくなることも確認しました。これらの原因について考察した結果、以下の結論に達しました。
・輸送補正適用時に(2)式を用いることにより、特に非等方散乱の影響の大きい水領域では、(2)式の右辺第二項の値が右辺第一項より大きくなることがある。
⇒ 輸送補正適用後の自群散乱断面積が負になることがある。
・自群散乱断面積が負の値になること自体は問題ないが、大幅に負側へ移行した場合、自群散乱断面積の絶対値が輸送断面積の値以上になることがある。
⇒ 散乱比の絶対値が 1 以上になることがある。
・散乱比の絶対値が 1 以上である領域への衝突確率が大きいとスペクトル半径が 1 以上になる可能性がある。
⇒ 輸送補正適用時に散乱比の絶対値が 1 以上になる可能性があるのは水領域だけなので、水領域の割合が大きくなると収束性は悪化する傾向にある。

炉物理の研究第 67 号（2015 年 3 月）

要は、「輸送補正の適用により水領域の自群散乱断面積が大幅に負側へ移行し、その絶対値が輸送断面積の値以上になること」が収束性悪化の根本原因であることを確認しました。
以上により収束性悪化の原因はわかりましたが、肝心なのはどう対処するかです。そこで、その対応策についても本研究にて考えました。それを次節にて簡単に紹介します。
2.3 収束性の改善
まず収束性の改善について考えるにあたり、中性子束の更新についての漸化式を(7)式のとおり考えます（(7)式は(6)式を(5)式に代入しただけのものです）。

rrrφ ( k +1) = Mφ ( k ) + r(7)

本研究では(7)式における中性子束の更新速度を(8)式のとおり加速因子ωを用いて調整することを考えました。

r( rr)rφ ( k +1) = ω Mφ ( k ) + r + (1 − ω )φ ( k )(8)

導出の詳細は省略しますが、(8)式のように中性子束の更新速度を調整することにより、反復行列の固有値は(9)式のようにシフトすることが導け、さらにその性質を使うと(10)式のとおり最適加速因子を求めることができます（このあたりを詳しく説明すると分量がかなり増えてしまいますので思い切ってスキップしますが、重要なのは(10)式で最適加速因子を求めることができるということです）。

λi′ = 1 − ω + ωλi

ωopt =

11
1 − (λmin + λmax )2

ただし、

λi ：元々の反復行列の i 番目の固有値λi′ ：中性子束の更新速度変更後の反復行列の i 番目の固有値λmin ：元々の反復行列の固有値のうち最も負側のもの

λmax ：元々の反復行列の固有値のうち最も正側のもの
ω opt ：最適加速因子

(9)
(10)

炉物理の研究第 67 号（2015 年 3 月）

それでは、図 1 のときと同条件で(10)式の最適加速因子を適用した場合のスペクトル半径を見てみましょう。それを図 2 に示します。図 2 より、最適加速因子を用いて中性子束の更新速度を調整すれば、輸送補正を適用した場合を含めスペクトル半径はすべて 1 未満になることが確認できます。実際に、AEGIS での計算が元々発散していたケースに対して最適加速因子を適用すれば、すべてのケースにおいて適切に収束解が得られるようになりました。

1.5
輸送補正なしスペクトル半径輸送補正あり1.00.5

0.0
050

100
エネルギー群150

図 2 最適加速因子適用後のスペクトル半径
ここまでが学会賞の受賞対象となった論文の概要です。ここまでで説明したとおり、本研究にて輸送補正適用による収束性悪化の原因を明らかにし、また収束性を改善する手法を開発しました。しかし、実はこの手法には大きな問題点がありました。それは、この手法を適用するには、衝突確率をすべて求めたうえで反復行列の固有値をすべて求める必要があるため、大きな体系に適用することが難しい点です（単一セル体系程度が限度）。です
が、そもそも収束性の検討を始めたのは臨界実験解析において計算が発散することがきっかけでしたので、ここで終わらせるわけにはいきません。そこで、本テーマについてさらに検討を進め、最終的に(11)式に示す加速因子を用いれば内部反復を確実に収束させられることを見出しました[5]。

ωstable =1
1 + max ( ci )i

(11)

炉物理の研究第 67 号（2015 年 3 月）(11)式において必要なものは散乱比のみであるため、臨界実験体系のような大きな体系であっても極めて容易に加速因子を求めることができます。実際に、AEGIS による臨界実験解析において試したところ、輸送補正適用時に元々すべての体系で発散していたのに対して、(11)式の加速因子を用いることによりすべての体系で収束解を得ることができました。
以上により、我々を長い間悩ませていた臨界実験解析における収束性の問題は無事に解決することができました。
4. おわりに
本稿で紹介した収束性の問題は以前から存在していたはずですが、少なくとも私の周囲では問題として認識されていませんでした。この問題が表面化した原因は、エネルギー群数の増加により断面積を細かく取り扱うようになったためと考えています。つまり、細かく取り扱うことにより輸送補正適用時の自群散乱断面積が大幅に負側へ移行するようなエネルギー群が出現し始めたということです。
「計算機が速くなったから計算条件を詳細にすれば良い」と単純に考えていると、意外なところでつまづく可能性があるということを思い知らされました。そういう意味でも本研究は私にとって大きな経験になったと思います。
それでは最後に、本稿で紹介した内容には少しマニアックな面もあったとは思いますが、特に少し分野の異なる方々や若手の方々を中心に、少しでも皆さんの刺激になることを期待して本稿を閉じたいと思います。
参考文献[1] M. Tabuchi, A. Yamamoto, T. Endo, N. Sugimura, “Convergence analysis of MOC
inner iterations with large negative self-scattering cross-section,” J. Nucl. Sci.

Technol., 50[5], pp.493-502 (2013).
[2] A. Yamamoto, T. Endo, M. Tabuchi, N. Sugimura, T. Ushio, M. Mori, M. Tatsumi, Y.
Ohoka, ”AEGIS: an advanced lattice physics code for light water reactor analyses,”,Nucl. Eng. Technol., 42[5], pp.500-519 (2010).
[3] A. Yamamoto, Y. Kitamura, Y. Yamane, “Simplified treatments of anisotopoc
scattering in LWR core calculations,” J. Nucl. Sci. Technol., 45[3], pp.217-229 (2008).
[4] Sartori E. Standard energy group structures of cross section libraries for reactor
shielding, reactor cell and fusion neutronics applications: VITAMIN-J, ECCO-33,
ECCO-Yvette Cedex, France: NEA Data Bank; December 11, 1990.
[5] M. Tabuchi, M. Tatsumi, A. Yamamoto, T. Endo, “Improvement of a convergence
technique for MOC calculation with large negative self-scattering cross section,”
Proc. PHYSOR 2014, Kyoto, Japan, Sep.28-Oct.3, 2014, (2014). [CD-ROM]