炉物理の研究

第 58 号 (2006 年 2 月)

<第 37 回炉物理夏期セミナー報告>（1）HTGR における二重非均質問題のためのセル計算コードの開発
東京工業大学

理工学研究科

原子核工学専攻
若菜

１．

俊一

背景・目的

高温ガス炉は安全性･信頼性、経済性、核不拡散性といった性質のもと、次世代炉として大きく期待されている。高温ガス炉は燃料として被覆粒子を用いているため燃料領域、冷却材領域、減速材領域との非均質性（巨視的非均質）に加えて、燃料領域内被覆粒子および黒鉛マトリクス間の非均質性（微視的非均質）も考慮する必要がある。日本で主に用いられている核計算コードシステム SRAC でも共鳴エネルギー領域では二重非均質性を考慮した衝突確率計算が行われている。しかし、プルトニウムや可燃性毒物として用いられている天然ガドリニウムなどの熱群エネルギー領域に吸収のピークを持つ物質を用いた場合
は熱群領域の中性子束に大きな誤差が生じるのではないかと考えられる。
そこで、これまでに報告されている二重非均質を取り扱ったいくつかの近似法を全群に応用して新たに衝突確率の合成を行い、全エネルギー領域での中性子束スペクトルを算出し、各近似法･計算法の妥当性を調べることを目的とする。
２．

理論

本研究では(1)日本原子力研究所で開発されたSRACの超詳細群衝突確率法に基づく
PEACOルーチンで用いられている合成法、(2)GA(Gulf General Atomic)社で開発された
MICROXコードで適用されている被覆燃料粒子に対する遮蔽因子の評価法、(3)LeslieとJonssonによって報告された衝突確率の合成法であるAccretion法[1]にWIMSコードなどで用いられているマルチセルの概念を適用した合成法、(4)APOLLO、APOLLO-2、DRAGON等
CEAの格子セル計算コードで適用されている二重非均質を考慮した衝突確率解析法[2]、の 4
つの合成手法を参考にしたが、今回は燃料粒子を二種類（濃縮ウランと兵器級プルトニウ
ム）にして計算を行ったため、複数の燃料粒子を含めても計算可能な(3)Leslie-Jonnson法と(4)APOLLO法の二通りで計算を行った。その合成法を以下に説明する。
(3)Leslie-Jonnson 法
Leslie と Jonsson によって報告された Accretion 法はクラスター型セルをモデルに、燃料、減速材、冷却材領域間の異なる巨視的領域間の衝突確率を次式で計算する。

PIJ (l , k ) = α lIα kJ

WI
PIJ
WlI

(I ≠ J )ここで PIJ (l , k ) は巨視的領域Iに存在する要素lから生じた中性子が巨視的領域Jの要素kにおいて最初の衝突をする確率を示している。また、PIJ は均質化された巨視的領域Iから生じた中性子が均質化された巨視的領域Jにおいて最初の衝突をする確率、 α lI は巨視的領域Iに入射してきた中性子が要素lで衝突をする割合、WlI は要素lの巨視的断面積ΣlIとその体積VlIの-48-

炉物理の研究

積で、WIは各要素での和第 58 号 (2006 年 2 月)

∑W を示す。また、同一燃料領域内にある被覆粒子、及びそれlIlを取り囲む黒鉛マトリクスの間の衝突確率は次式で計算される。

PII (l , k ) = QI (l , k ) −WI(1 − PII )α lIα kI
WlI

ここで QI (l , k ) は仮想無限格子の衝突確率であり、巨視的領域 I の全ての要素が無限に広がった状態での要素 l から生じた中性子が要素 k で衝突する確率を示している。この仮想無限格子の衝突確率 Q の算出の際にマルチセル方式を適用させて、複数種類の被覆粒子でも計算を可能にした。
(4)APOLLO 法
衝突確率の合成を求める計算式は記号等が複雑なので省略する。APOLLO 法の特徴は巨
視的領域内の要素をさらにシェルとして領域分割することができ、可燃性毒物のような外側から燃えていくような粒子に対しても詳細に中性子束を計算することが可能となっている。また、いくつかある近似の中で、微視的要素は Diluent（黒鉛マトリクス）で囲まれているものとし、この微視的要素から逃れ出た中性子で、ある任意の領域で最初に衝突する
確率は、Diluent で一様、等方に発生した中性子が同様な衝突を行う確率に等しい、というものが特徴的である。
３．

計算条件

解析対象炉心はブロック型ヘリウム冷却高温ガス炉である高温工学試験研究炉 HTTR をモデルとした。その詳細パラメータを表に示す。前述したとおり、今回は二種類の燃料粒
子（濃縮ウランと兵器級プルトニウム）を用いて計算を行った。
表-1

４．

炉心パラメータ

Fuel Cell TypePin in Block Type (HTTR Type)

Fuel TypeTRISO Coated Fuel Particle

Kernel Diameter [mm]0.60

Packing Fraction [%]30.0

Compact Inner / Outer Diameter/ Sleeve Outer Diameter /Block Inner Diameter [cm]1.00 / 2.60 / 3.40 / 4.10

Fuel Cell Pitch [cm]6.60First Fuel Kernel Type / Fraction [%]15% Enriched Uranium / 50

Second Fuel Kernel Type / Fraction [%]Weapon Grade Plutonium / 50

Core Temperature [K]1040

Block Temperature [K]

950

結果

まずは一種類の燃料粒子（兵器級プルトニウム）を用いた場合において、各合成法の妥-49-

炉物理の研究第 58 号 (2006 年 2 月)

当性を示す。図 1 においては通常用いられている SRAC コードと MVP モンテカルロコー
ドにおける結果を比較したものである。結果は中性子束スペクトルではなく、燃料コンパ
クト内の中性子束を黒鉛ブロックの中性子束で割った Flux Depression という形で示した。
前述したとおり、プルトニウムは熱群に強い吸収のピークを持っているため、熱群領域で二重非均質性を考慮していない場合は MVP との間に大きな誤差が生じてしまっている。しかし、Leslie-Jonnson 法や APOLLO 法を用いた場合、図 2 に示すように吸収のピークが
ほぼ MVP と一致した結果となった。これにより一種類の燃料粒子を用いた場合の各合成法
は MVP を再現することができるものであるといえる。

Flux R atio(Fuel C o m pact/ G ra phite B lock)[-]

1.6
1.4
1.2
1.0
0.8
0.6

0.2
0.0
1.E-03

図1

MVP
SRAC(original)

0.4

1.E-02

1.E-01

1.E+00

1.E+01

1.E+02 1.E+03
Energy[eV]1.E+04

1.E+05

1.E+06

1.E+07

既存の SRAC コードでの燃料コンパクト内の Flux Depression

Flux Ratio(Fuel Co mpact/ Graphite Block) [-]

1.6
1.4
1.2
1.0
0.8
0.6

0.2
0.0
1.E-03

図2

MVP
APOLLO

0.4

1.E-02

1.E-01

1.E+00

1.E+01

1.E+02 1.E+03
Energy[eV]1.E+04

1.E+05

1.E+06

1.E+07

APOLLO 法を用いた場合での燃料コンパクト内の Flux Depression

次に、同一燃料領域に二種類の被覆燃料粒子を含めた場合の結果を示す。図 3 はLeslie-Jonsson法による合成法を用いた場合の兵器級プルトニウム燃料粒子内におけるFlux Depressionを示し、図 4 はAPOLLO法による合成法を用いた同様の結果を示している。
10eV〜107eVにかけてMVPとの間に僅かな差は生じているものの、熱群領域での吸収のピークはどちらの合成法においてもMVPとほぼ一致していることがわかった。

-50-

炉物理の研究第 58 号 (2006 年 2 月)

1.6
1.4

Flux Depression

1.2
1.0
0.8
0.6

MVP
Leslie-Jonsson

0.4
0.2
0.0
1.E-03

図3

1.E-01

1.E+01
1.E+03
Energy[eV]1.E+05

1.E+07

Leslie-Jonsson 法による兵器級プルトニウム粒子内の Flux Depression
1.6
1.4

Flux Depression

1.2
1.0
0.8
0.6

MVP
APOLLO

0.4
0.2
0.0
1.E-03図4
５．

1.E-01

1.E+01
1.E+03
Energy[eV]1.E+05

1.E+07

APOLLO 法による兵器級プルトニウム粒子内の Flux Depression

今後の方針

今回の結果から Leslie-Jonsson 法及び APOLLO 法は一種類の燃料粒子における MVP の再現性に加えて、二種類の燃料粒子の場合でも同様な精度を与えることが可能であることがわかった。高速領域での MVP における中性子束の方が僅かに大きくなってしまったことに関してはもう少し検討する必要がある。
今後はこれらの合成法を用いた燃焼計算用群定数が作成できるように検討を進めていく。
参考文献[1] D. C. Leslie , A. Jonsson, “The Calculation of Collision Probabilities in Cluster-Type
Fuel Elements” , Nucl Sci Engin. Vol. 23, pp272-290 (1965)[2] A. Hebert, “A Collision Probability Analysis of the Double-Heterogeneity Problem” ,Nucl Sci Engin, Vol. 115, pp177-184 (1993)-51-

炉物理の研究第 58 号 (2006 年 2 月)

〈第 37 回炉物理夏期セミナー報告〉(2)原子炉シミュレーションゲームの開発
東北大学大学院工学研究科

量子エネルギー工学専攻
遠藤

秀樹[1]

背景
現在、日本では国内電力の約三割を原子力発電によってまかなっている。この原子力発
電は、石油・石炭などの化石燃料による発電に比べ、高い経済性を有する点、発電時に二酸化炭素を排出しない点からその有用性はゆるがしがたく、今後も着実な利用を図ってい
く必要がある。しかしながら、原子力発電は核分裂物質を利用するがゆえに、臨界事故への不安感や核兵器への恐れなどから、唯一の被爆国である日本国民には受け入れがたいものとなっている。また、我々の日常生活に欠かすことのできない電力についてあまり興味をもたれていないことも大きな問題である。
そこで、原子力発電をより身近に感じてもらうため、原子力発電所の運転を体験することができる原子炉シミュレーションゲームの開発を行った。この原子炉シミュレーション
ゲームは東北大学オープンキャンパスで公開し、より多くの人に原子力発電に興味を持ってもらうことを目的としている。
開発目的
本プログラムの開発目的としては背景でも述べたように、東北大学オープンキャンパスで公開し、より多くの人に原子力発電に興味を持ってもらうことである。本プログラムは
原子炉の運転を体験することができるものとした。そのため、本プログラムは原子炉内の物理現象の経時変化を計算する機能を有する必要がある。また、本プログラムはオープン
キャンパスでの公開を目的としたものであるため、一般の方にわかりやすいものでなければいけない。したがって、本プログラムは実際の原子炉をそのままシミュレーションするのではなく、必要な部分だけをデフォルメしてわかりやすいプログラムを作成しなければならないと考えられる。さらに、原子炉の運転を体験してもらう際、運転にかかる時間が長くなり過ぎない必要がある。そのため、原子炉内の物理現象をある程度ゲームとして簡単化する必要がある。加えて運転が難しすぎず簡単すぎず、ゲームとして面白い必要がある。
原子炉シミュレーションゲームHide-yoshi
開発目的で述べたような条件を満たす原子炉シミュレーションゲームHide-yoshi[2]の開発を行った。開発環境としてはWindows上でLabVIEW[3]を用いて開発を行った。また、基[1] hide@neutron.qse.tohoku.ac.jp[2] HIDEki and YOSHIkawa’s nuclear reactor simulation game version 1.00
[3] LabVIEWは信号の収録、測定値の解析、データの表示機能を標準装備した画期的なグラフィ

カル開発環境である。詳しくは参考文献[5]を参照。

-52-

炉物理の研究第 58 号 (2006 年 2 月)

本的には次式に示すような一点炉動特性方程式を解くことで原子炉出力の経時変化を計算している。

dn(t ) k (t )(1 − β ) − 1 6=+ ∑ λi C i (t )
dt
l
i =1

(1)

dC i (t ) k (t )=β i n(t ) − λi C i (t )
dt
Λ

(2)

(i = 1〜6)

ここで、

n(t )
Ci (t )k (t )
l

β
βi
λi

：

中性子密度[n/cm3]：

遅発中性子先行核濃度[n/cm3]

：

増倍係数

：

即発中性子寿命[sec]：遅発中性子割合：

各郡の遅発中性子割合：

各群の遅発中性子崩壊定数[sec-1]である。このとき、増倍係数 k (t ) は制御棒位置によって決まる。また、中性子密度 n(t ) は
原子炉熱出力に比例するため、この一点炉動特性方程式を解き n(t ) を求めることで原子炉出力の経時変化を得ることが可能となっている。また、他のパラメータに関しては原子炉の形状、および組成によって決まるものであり、これらのパラメータは Hide-yoshi の入力
ファイルを作成することで変更することが可能となっている(インストール時には女川原子力発電所のパラメータの入力ファイルが用いられる)。
原子炉シミュレーションゲーム Hide-yoshi の操作パネルを図 1 に示す。
操作パネルには、①制御棒位置・再循環ポンプ流量表示部、②原子炉状態グラフィカル表示部、③現在の運転パラメータ(時間、ペリオド、冷却水温度、原子炉熱出力)表示部、④運転パラメータグラ
フ表示部、⑤スクラム動作ボタン、⑥ゲームミッション表示部、⑦運転状況表示部が存在

図 1 Hide-yoshi 操作パネル-53-

炉物理の研究第 58 号 (2006 年 2 月)する。表示されているこれらのパラメータを見ながら運転を行うことになる。また、運転パラメータ表示部にはミッション表示部で示された目標値となる温度および出力に青い帯が表示されており、この帯の中に納まるように運転を行うことでゲームクリアとなる。
また、操作はテンキーを改良した専用コントローラまたはテンキーで行うことができる。
専用コントローラを図 2 に示す。また、テンキーを用いる場合のボタン配置は、7、8、9で制御棒挿入、4、5、6 で制御棒引き抜き、1、3 でポンプ流量増加、0、. でポンプ流量減少となる。
ゲームとしては、制御棒位置および再循環ポンプ流量を変えることによって目標としている冷却水温度、および原子炉熱出力まで運転してもらい、一定時間保持することでミッションクリアとなる。このとき、ペリオドが短すぎたり、温度が高すぎたり、出力が高すぎる場合には警告をされることがあり、その危険度によって警告、制御棒自動挿入、スク
ラムと安全装置が自動的に働く。また、ゲームクリア時には、運転結果がプリンターから出力される。プリントされる運転結果を図 3 に示す。運転結果には①最終的な原子炉の状態のグラフィック、②運転時間、③警告回数、④運転結果のグラフ(増倍係数、制御棒位置、
再循環ポンプ流量、熱出力、冷却水温度)、⑤評価(運転結果によって 5 段階評価される)が出力されている。
また、本プログラム、原子炉シミュレーションゲーム Hide-yoshi は中性子デバイス研究
室の web site(http://www.neutron.qse.tohoku.ac.jp/)からダウンロード可能である。

図 2 Hide-yoshi 専用コントローラ

図 3 Hide-yoshi 運転結果プリント-54-

炉物理の研究第 58 号 (2006 年 2 月)

オープンキャンパス
本プログラム、原子炉シミュレーションゲーム Hide-yoshi を 2005 年 7 月 28 日、29 日に行われた東北大学オープンキャンパスで公開した。その様子を図 4 から図 7 に示す。オープンキャンパスでは、まずポスターで原子力発電の現状や仕組みなどを説明し、その後に実際に Hide-yoshi を運転してもらった。来場者は主に高校生であったが、それほど抵抗なく運転できていたように思う。また、図 7 に示すように、ゲームクリア時間を競ってもらいランキングとして張り出すことによってよりゲームとして面白くなったのではないかと思う。
来場者は 100 人を越え、多くの人に原子力発電、ならびに電力について興味をもっていただけたのではないかと思われる。来年度以降もこの Hide-yoshi にさらなる改良を加え、毎年公開していく予定である。

図 4 オープンキャンパス(1)

図 5 オープンキャンパス(2)

図 6 オープンキャンパス(3)

図 7 オープンキャンパス(4)

まとめと今後の課題
今回、原子力発電に興味を持ってもらうことを目的として、原子炉シミュレーションゲ-55-

炉物理の研究第 58 号 (2006 年 2 月)

ーム Hide-yoshi の開発を行った。また、本プログラムを東北大学オープンキャンパスで公開し、原子炉の運転を体験してもらうことによって多くの人に原子力発電に触れてもらった。これによって実際に体験した方、またこれらの人から話と聞いた方々が原子力発電、
ついては電力に興味を持ってもらったのではないかと思う。さらに、本プログラムはこれから毎年オープンキャンパスで公開していく予定であり、毎年バージョンアップしていく
予定である。
今後の課題としては、Hide-yoshi ではかなり大きな近似を用いているため、ゲーム性を失わない程度の、原子炉運転のより正確なシミュレーションを行う必要がある。また、よりわかりやすい操作パネルへの改良、ゲーム的演出効果の追加、および実際の原子炉運転に用いられているコントローラに似た専用コントローラの作成などが挙げられる。
今後、本プログラムによってより多くの人が原子力発電に興味を持ち、さらには原子力発電の必要性を多くの人に知ってもらうことによって、この原子力の分野がよりいっそう進歩していくことを期待する。
参考文献[1] 平川直弘,岩崎智彦,”原子炉物理入門”,東北大学出版会,(2003)[2] ジェームス J ドゥデルスタット,ルイス J ハミルトン著,成田正邦,藤田文行共訳,”原子炉の理論と解析(上)”,現代工学社,(1980)[3]遠藤秀樹,”加速器運転機能を有する加速器駆動未臨界炉シミュレータの開発”,東北大学
工学部学士学位論文,(2004)[4]遠藤秀樹,”動特性簡易解析コード Hide-kine2.0 マニュアル”,東北大学大学院工学研究科,(2005)[5] LabVIEW website：http://www.ni.com/labview/ja/-56-

炉物理の研究第 58 号 (2006 年 2 月)

<第 37 回炉物理夏期セミナー報告>（3）炉心計算への Krylov 部分空間法の応用
名古屋大学大学院

工学研究科

マテリアル理工学専攻

博士前期課程 2 年(山根研究室)
1

田渕将人

背景・目的
炉心計算の高精度化、高速化は運転の経済性向上の観点から非常に重要である。近年の計算機の発達に伴い、炉心計算手法は拡散から輸送へ移行していく傾向にある。しかし、輸送計算においても収束の加速に拡散計算を用いるなど、依然として拡散計算の持つ重要性は高い。そこで、本研究では拡散方程式の解法に注目し、炉物理の分野ではあまり知られていない Krylov 部分空間法の適用を試みた。
「拡散方程式を数値的に解く」というのは拡散方程式を行列で表現したときの固有値、固有ベクトルを求めることである。通常、この固有値計算は反復解法（べき乗法）によって行われ、この反復のことを外部反復と呼ぶ。外部反復の途中で、中性子源から中性子束を求める計算を行っているが、この計算は非斉次の連立方程式（式の数と未知数の数が等
しい、いわゆる普通の連立方程式。以下では「非斉次」を省略）を解く計算であり、通常
この連立方程式は反復解法によって解かれ、この反復のことを内部反復と呼ぶ。本研究では内部反復だけに注目することにする。内部反復の収束性を良くすることは、拡散計算全体の収束を速くすることにつながるので非常に重要である。現在、内部反復に用いられている手法は、計算条件によっては収束が悪くなることが知られている。そこで、本研究では応用数理学の分野で連立方程式の解法としてよく知られている Krylov 部分空間法を内部反復へ適用することを試みた。
本研究の目的としては、Krylov 部分空間法を炉物理の問題に適用したときの有用性を評価することと、炉物理の分野の方々に Krylov 部分空間法を知っていただくことである。筆者としては後者の方が重要だと考えているので、この資料を読んで Krylov 部分空間法に少しでも興味を持っていただければ幸いである。
2

連立方程式の解法
上で述べた通り、内部反復とは連立方程式を解く操作に他ならない。本研究では内部反復のみに注目するので、まずは連立方程式の解法について簡単に説明する。
連立方程式の解法は直接法と反復法の 2 つに大別できる。また、反復法は定常反復法と非定常反復法に分けられる。これらを簡単にまとめたものが Fig.2.1 である。
直接法：Gauss-Jordan法, Gaussの消去法, LU分解法, Cholesky分解法など
定常反復法：Jacobi法, Gauss-Seidel法, SOR法など

連立方程式の解法
反復法

非定常反復法：Krylov部分空間法など

Fig.2.1 連立方程式の解法の分類-57-

炉物理の研究第 58 号 (2006 年 2 月)

以下では、それぞれの計算手順や特徴について簡単に説明する。
2.1

直接法

一般的な連立方程式(2.1)式に対し、直接法では両辺に係数行列の逆行列をかけて(2.2)
式のように解を求める。

Ax = bx = A −1 b

(2.1)
(2.2)

連立方程式の元をnとすると、計算量はn3、メモリはn2のオーダーで増加することが知られており、サイズの大きい連立方程式を解くのには向いていない。しかし、nが小さいときについては反復法よりも有効な場合もある。一般的に、偏微分方程式の数値解法では、疎行列を扱うことになるが、直接法では係数行列の疎性をうまく利用できないことがある。したがって、拡散計算のように大規模疎行列を係数行列に持つ問題に対して直接法はあまり適していない。
2.2

定常反復法

定常反復法では(2.3)式のような漸化式への代入を繰り返すことで解を求める。

x( k +1)

= Mx

(k )

+c

(2.3)

反復行列 M のスペクトル半径（最大固有値の絶対値）によって収束の速さが決まり、反復行列のスペクトル半径が 1 を超えると x は発散する。係数行列 A の対角成分が小さくなると反復行列のスペクトル半径は大きくなる傾向にある。逆に係数行列が対角優位（対角成分が他の成分に対して大きい状態）であれば反復行列のスペクトル半径は小さ
くなり、収束は速くなる。
2.3

非定常反復法

定常反復法では(2.3)式のように行列 M やベクトル c は k に依存せず一定であった。これに対し、非定常反復法では(2.4)式ように k に依存するαや p を用いる。

x( k +1)=x(k )+ αk pk(2.4)

非定常反復法では連立方程式の元よりも小さな次元の部分空間を考え、その部分空間を使って近似解を求め、その近似解を利用しながら部分空間の次元を上げていくことで解を求める。この操作を図で示したものが Fig.2.2 である。反復回数 k が部分空間の次元を表しており、丸め誤差がなければ、n 元の連立方程式は最悪でも n 回の反復で 100%解けることになる（n 回で解けるということから非定常反復法は反復法ではないと主張する人もいるらしい）。非定常反復法は基本的に定常反復法と収束のさせ方が異なり、スペク
トル半径のような概念はない。したがって、収束の速さは係数行列が対角優位であるかどうかには依存しない。しかし、収束の速さが何に依存するかは、まだよく知られていない（筆者が知らないだけかもしれない）。

-58-

炉物理の研究第 58 号 (2006 年 2 月)

N元の連立方程式ならN次元空間のどこかに解があるはず

N次元空間

x
xk

k次元
部分空間

x0

1次元
部分空間

zk

でも見当もつかない
ので適当なx0を選ぶ

x1

1次元部分空間内のz1を使い
x1 = x0 + z1 が解に最も近づくような z1 を探す

z1

部分空間内の次元を大きくし
xk = x0 + z kが最も解に近づく
ような z k を順に探していく
Fig.2.2 非定常反復法のイメージ

3

Krylov部分空間法の概要
3.1

Krylov 部分空間の定義

正則なN×N行列Aと非ゼロベクトルvによって作られるベクトル列 { v, Av, A2v, …,
Am−1v } が張る部分空間をKrylov部分空間と呼び、Km( A ; v )と表記する。つまり(3.1)
式で表せるような部分空間をKrylov部分空間と呼ぶ。

( ){K m A; v := Span v, Av, Α 2 v,L, A m −1 v
3.2

}

(3.1)

Krylov 部分空間法

Krylov部分空間法とは、Fig.2.2 における部分空間にKrylov部分空間を用いる手法である。まず、初期近似解x0から初期残差r0=b−Ax0を求め、r0とAからKrylov部分空間Km( A ;r0 )を作る。mが部分空間の次元を表しており、mを大きくしていくにつれ近似解xmは解に近づく。近似解xmは次の(3.2)式で表される空間条件を満たすように作られる。

(

x m = x 0 + z m , z m ∈ K m A; r 0

)

(3.2)

しかし、(3.2)式だけではxmを唯一に決められないため、通常は(3.3)式で表される最小条件、または(3.4)式で表される直交条件を課し、xmを決定する。

rm =

min

( )

x∈x m + K m A ; r 0(

r m ⊥ K m A; r 0

b − Ax

(3.3)

)

(3.4)

Krylov 部分空間法の構築には、Krylov 部分空間の正規直交基底の計算と残差のアプロ-59-炉物理の研究第 58 号 (2006 年 2 月)

ーチの 2 つの操作が必要となる。理想的な Krylov 部分空間法とは、①計算量やメモリ量が少なく②(3.3)式または(3.4)式を完全に満たすことであるが、一般的な問題に対して理想的である反復解法は構築できないことが証明されている。したがって、理想的な特徴の内、①、②のどちらか一方だけを採用し、もう一方を他の物で代替することになる。
現 在 よ く 用 い ら れ て い る 方 法 と し て 、 Bi-Conjugate Gradient 法 （ Bi-CG 法 ） や
Generalized minimal residual 法（GMRES 法）が挙げられるが、Bi-CG 法は①を採用、GMRES 法は②を採用した方法であり、それぞれ異なった特徴を持つ。
ここまでの説明だけでは、実際に何をどうすればよいのかあまり分からないと思うが、ページ数の都合上これ以上の説明は行わない。より詳細な説明は参考文献[1]を参考にしていただきたい。
3.3

拡散計算と Krylov 部分空間法

現在の拡散計算の内部反復によく用いられている手法はSOR法のような定常反復法である。上でも述べた通り、定常反復法の弱点は反復行列（(2.3)式中のM）のスペクトル半径によって収束の速さが決まってしまうことである。一般的に、係数行列（(2.1)式のA）が対角優位であれば反復行列のスペクトル半径は小さくなり、収束は速くなる。拡散計
算においては、拡散係数によって係数行列の対角成分の優位性が決まる。拡散係数が小さいほど対角成分は優位になり、拡散係数が大きいほど対角成分は優位でなくなる。また、拡散係数以外の要因でも係数行列の対角成分が小さくなることはいくつかあり、そのような場合には定常反復法による内部反復は収束しにくくなる（場合によっては発散
することもある）[2],[3]。そこで、特に係数行列の対角成分が優位でないとき、Krylov部分
空間法による内部反復の有用性が期待できる。次の章では、炉心計算の内部反復を想定した計算を通して、定常反復法であるSOR法とKrylov部分空間法の比較を行う。①通常の炉心計算についての係数行列を持つ場合 ②局所的に対角優位でない係数行列を持つ
場合 ③全体的に対角優位でない係数行列を持つ場合の 3 つのケースについて紹介する。
4

Krylov部分空間法を用いた計算
4.1

計算の手順と計算条件

拡散での炉心計算を行う際の内部反復を想定した計算を行った。計算の流れは以下の通りである。まず、断面積、体系の大きさ、メッシュなどのデータから内部反復で解く
べき連立方程式の係数行列を計算する。中性子源（連立方程式の右辺）は収束に関係ないので、適当に与える（今回はすべて 1.0 とした）。そして、連立方程式を解くときの収束性（反復回数と残差の関係、収束までの反復回数）、計算時間を調べた。

計算条件は以下の通りである。
計算体系：PWR, 3LOOP 炉心(157 集合体), 二次元-60-炉物理の研究第 58 号 (2006 年 2 月)

エネルギー群：2 群
計算手法：SOR 法（現在の炉物理でよく用いられる手法）と GPBi-CG 法（Krylov 部分
空間法の一種）加速因子：SOR 法では反復ごとに近似解の値を外挿しており、外挿の度合いを決める加速因子というパラメータがある。加速因子によって収束性は大きく変わるが、今回は加速因子をいろいろ変えて何通りも計算を行い、最適と思われる加速因子を調べ、その加速因子での計算結果を採用
前処理：Krylov 部分空間法では収束性を良くするため、係数行列を変形する操作を行う。
これを前処理と呼ぶ。専門家によると、今回の問題で扱う行列には不完全 LU 分
解による前処理がいいらしいが、筆者が理解するには至らず、今回の計算は前処理なしで行った。
以下の 3 つのケースについて計算を行った。
ケース 1：1 集合体あたり 34×34。断面積は通常の PWR の断面積を使用。詳細にメッシュを切った場合について調べるのが目的。
ケース 2：1 集合体あたり 10×10 のメッシュ。断面積については炉心の真ん中の 1 体以外（156 体）は通常の断面積を使うが、真ん中の 1 体については拡散係数を変化させた。係数行列の対角成分が局所的に変化する場合の収束性を調べるのが目的。
ケース 3：1 集合体あたり 10×10 のメッシュ。すべて同一の燃料集合体で構成される炉心を想定し（つまり炉心全部が同じ断面積）、炉心全体の拡散係数を変化させた。係数行列が全体的に変化する場合の収束性を調べるのが目的。
4.2

計算結果

4.2.1 ケース 1
ケース 1 の結果について収束（残差のノルムが 10−8）までの反復回数、計算時間をTable4.1 にまとめる。ただし、筆者が作成したGPBi-CG法のプログラムにはif文が大量に使われており、そこを改善すればGPBi-CG法での 1 反復あたりの計算時間は短縮できる可能性がある。
Table4.1 ケース 1 のまとめ

1群
2群
GPBi-CGSOR
GPBi-CG
SOR
収束までの反復回数
210
4249
271
1470
全計算時間[s]771
925
1反復あたりの計算時間[s] 3.33E-02 1.67E-02 3.32E-02 1.70E-02

4.2.2 ケース 2
真ん中の集合体の拡散係数と収束（残差のノルムが 10−8）までの反復回数をTable4.2-61-炉物理の研究第 58 号 (2006 年 2 月)にまとめる。
Table4.2 ケース 2 のまとめ

拡散係数[cm]
1.0E-02
1.0E-01
1.0E+00
1.0E+01
1.0E+02
1.0E+03
1.0E+04
1.0E+05

1群
GPBi-CG
116
71
61
129
408
1124
2155
6173

2群
SOR
357
356
355
746
10252
1071421076178≧10000000

GPBi-CG
145
93
79
197
683
1755
2699
6343

SOR
129
129
129
722
12860
1400341412750≧10000000

4.2.3 ケース 3
炉心全体の拡散係数と収束（残差のノルムが 10−8）までの反復回数をTable4.3 にまとめる。ただし、このケースでは 1 群、2 群とも同じ断面積になっているので 1 群、2群の区別はない。
Table4.3 ケース 3 のまとめ

拡散係数[cm]
1.0E-02
1.0E-01
1.0E+00
1.0E+01
1.0E+02
1.0E+03
1.0E+04
1.0E+05
4.3

GPBi-CG
2
3
12
125
247
280
292
282

SOR
3
5
28
1882
89134
1218146≧10000000≧10000000

考察

4.3.1 ケース 1収束の速度を GPBi-CG と SOR で比較すると GPBi-CG の方がかなり速く収束する
ことが分かる。1 群と 2 群で比較すると、GPBi-CG では収束の速度に大きな差は見られないが、SOR では 2 群の収束の方が 3 倍程度速い。これは 1 群の方が 2 群よりも拡散係数が大きく、1 群の計算では係数行列の対角成分が小さくなることが原因と考えられる。これらのことから、GPBi-CG が係数行列の対角成分に依存しないことが確認できた。また、GPBi-CG の有用性が確認できた。
4.3.2 ケース 2
真ん中の集合体の拡散係数を変化させると、GPBi-CG、SOR ともに収束性は悪化する。GPBi-CG は対角成分の大きさに依存しないはずであるので、収束性が悪化する原因は明確には分からないが、係数行列の全体のバランスが悪いと、収束性が悪化するのではないかと推測している。SOR 法で収束性が悪化する原因はやはり係数行列の対角成分が小さくなることであると考えられる（ただ、反復回数が拡散係数の増加とと-62-

炉物理の研究第 58 号 (2006 年 2 月)

もに単調増加しているわけではなく、変化の特徴が GPBi-CG と少し似ているので、係数行列の対角成分だけが原因ではない可能性が高い）。いずれにせよ、GPBi-CG は SOR
より常に収束性が良いという結果になった。
4.3.3 ケース 3
炉心全体の拡散係数を変化させると、ケース 2 と同様、GPBi-CG、SOR ともに収束性が悪化する。ケース 2 と違う点は、拡散係数を大きくしたとき、GPBi-CG の反復回
数は 300 回程度の値に飽和する点である。GPBi-CG の収束性は対角優位かどうかに依存しないはずであるが、拡散係数を変化させたときに収束性が若干悪くなるのは丸め
誤差が蓄積しやすくなるからではないかと考えている（ケース 2 についてもそうかもしれない）。SOR の収束性が悪くなる原因は言うまでもなく対角成分が小さくなることによるものである。
4.3.4 全体を通して
SOR は拡散係数に対する（係数行列の対角成分に対する）依存性が非常に強い。
GPBi-CG も拡散係数に若干依存したが、「対角成分が小さいから」というよりは「行列の全体のバランスが悪いから」という理由で収束性が悪化していたような印象を受けた。いずれにせよ、すべてのケースについて GPBi-CG の方が良かったことには間違いない。また、専門家によると、GPBi-CG に前処理をすることで収束速度は格段に速くなる。さらには前処理によって、拡散係数にほとんど依存しない GPBi-CG も作れるのではないかと期待している（あくまで個人的に期待しているだけで根拠はない）。
5

まとめと今後の課題
Krylov部分空間法を用いた計算を通して、その有用性が確認できた。どのケースについてもKrylov部分空間法の方がSOR法に比べて良好な結果を示し、拡散係数が大きいときに
特にSORとの差が顕著になった。これらの結果を通して、Krylov部分空間法を炉心計算へ適用することは十分に検討する価値があると言える。今後の課題としては、前処理つきの
GPBi-CG法を用いた計算やWielandt法（行列を変形して外部反復を加速するが、行列の変形によって係数行列の対角成分が小さくなる手法）[3]との併用などを検討する。
6

謝辞
この資料の 2 章、3 章の多くは、名古屋大学計算理工学専攻の張紹良教授、山本有作講師の助言を基にして書いたものである。両氏は数値計算のスペシャリストであり、大規模疎行列の固有値計算や並列計算などに非常に詳しい。特に張教授は参考文献[1]の著者であり、GPBi-CG 法の考案者でもある。本研究が行えたのは Krylov 部分空間法について両氏が親切に教えて下さったおかげであり、心から感謝したい。
7

参考文献[1] 藤野清次,張紹良著『反復法の数理』朝倉書店,1996 年-63-炉物理の研究第 58 号 (2006 年 2 月)

[2] A. Yamamoto, “Convergence Improvement of Response Matrix Method with Large
Discontinuity Factors,” Nucl. Sci. Eng., 145, 291 (2003)[3] T. M. Sutton, “Wielandt Iteration as Applied to the Nodal Expansion Method,”

Nucl. Sci. Eng., 98, 169 (1988)-64-