日本原子力学会・炉物理部会

炉物理部会ニュース(No.20)
2005 年 9 月 12 日発行

目次<第３７回炉物理夏期セミナー報告>
全体概要報告…………………………………………………………………………………………….2
東京工業大学原子炉工学研究所
関本

博、小原

徹

若手研究会報告(1)HTGRにおける二重非均質問題のためのセル計算コードの開発…………………………….8
東京工業大学大学院

理工学研究科

原子核工学専攻
若菜俊一(2)原子炉運転シミュレーターゲームの開発討……………………………………………………12
東北大学大学院

工学研究科

量子エネルギー工学専攻
遠藤 秀樹(3) 炉心計算へのKrylov部分空間法の応用……………………………………………………….17
名古屋大学大学院 工学研究科 マテリアル理工学専攻
田淵将人

-1-

＜第 37 回炉物理夏期セミナー報告＞
東京工業大学原子炉工学研究所
関本

博、小原

徹

第 37 回炉物理夏期セミナーは平成 17 年 8 月 9 日から 11 日までの 3 日間の日程で、栃木県日光市のホテル清晃苑で開催された。今回の炉物理夏期セミナーは、東京工業大学が幹事機関となり、東京工業大学 21 世紀 COE プログラム「世界の持続的発展を支える革新的原子力」(COE-INES)との共催として開催された。参加者は、学生 27 名、社会人 17 名、講師・幹事 8 名の計 52 名であった。
次の時代を担う若い世代の育成は炉物理分野において重要な課題の一つである。昨今の種々の状況の変化を反映して、学会等においても炉物理の教育のあり方が活発に議論されている。
今回の炉物理夏期セミナーは、このような炉物理教育に対する関心の高まりを受け、「最新炉物理講義」と題して、各大学において現在炉物理教育に携わっている教員による炉物理各分野の最新の講義と炉物理教育に関する討論会を行うことによって、若手研究者・技術者、学生に対しては各大学の第一線の炉物理の講義に触れる機会を、そして教育に携わる側に対しては炉物理教育のあり方について考える機会を提供することを目的として企画された。
セミナーでは国内各大学の教員による講義に加え、COE-INES の支援により米国カリフォルニ
ア大学バークレー校(UCB)原子力工学科学科長の Jasmina Vujic 教授を招聘し、現在 UCB で行われている炉物理数値解析手法に関する最新の講義とパネルディスカッションにおいて米国
特に UCB における原子力工学教育の現状と展望の講演が行われた。また昨年から始まった技術士試験「原子力・放射線部門」の解説も行われた。セミナー2 日目の夜には、若手研究会が開催され、活発な研究発表と討論のみならず、若手同士の懇親を深める良い機会となった。
参加者のアンケートを見ると、非常に分かりやすく良かった、炉物理の再勉強ができてよかった、有意義な講義だった、等の意見が多く各講師の先生方の工夫が効を奏して参加者へ強い印象を残したようである。また、他機関・他大学に人と交流を深められた、との意見も多かった。幹事としては今回のセミナーが当初の目的を少しでも達成できていたとすれば幸いに思う。
今回のセミナー会場は日光輪王寺の山内にあり、世界遺産に登録された東照宮や二荒山神社
へも歩いて数分という好立地であった。セミナー期間中は猛暑の下界と異なり涼しく快適で、若干の雨が杉林の中の寺社の神秘的な雰囲気をかもしだしていた。講義の合間の自由時間にこれら寺社の散策を楽しんだ参加者も多かったようである。
最後に、セミナーで講演いただいた講師の先生方、参加者の皆様、セミナー開催にご協力いただいた賛助企業の方々、セミナー準備にあたった COE-INES 事務局関係者に深く感謝の意を
表したい。

-2-

・セミナー実施概要
日程：

2005 年 8 月 9 日（火）〜8 月 11 日（木）場所：

ホテル清晃苑（栃木県日光市）

テーマ：

最新炉物理講義

参加人数：

52 名（内訳：講師 7 名、部会員 9 名、正会員 0 名、一般 8 名、学生部会員 13 名、学生会員 0 名、学生 14 名、幹事 1 名）

共催：

東京工業大学 21 世紀 COE プログラム「世界の持続的発展を支える
革新的原子力」(COE-INES)

・セミナープログラム
８月９日(火)１３：３０−１４：００

受付

１４：００−１４：１０

開校式

１４：１０−１５：３０

中性子輸送理論 講師

東北大

岩崎智彦

（座長：名大 山本章夫）
１５：５０−１７：１０

中性子拡散理論 講師

名大（座長：東工大 小原山本章夫

１８：３０−２０：３０

徹）

懇親会

８月１０日(水)
９：００−１０：２０

高速中性子減速理論

講師（座長：京大 三澤

東工大

１０：４０−１２：００

熱中性子理論

小原

徹

毅）

講師

京大

三澤

毅

（座長：名大 山根義宏）
１４：３０−１５：５０

原子炉動特性理論

講師

北大

島津洋一郎

（座長：九大 工藤和彦）
１６：１０−１７：４０

技術士試験「原子力・放射線部門」解説
講師

九大

工藤和彦

（座長：北大 島津洋一郎）
１９：３０−２１：３０

若手研究会

８月１１日(木)
９：００−１１：００ 「COE-INES セミナー」
Numerical Methods for Solving Neutron Transport/Diffusion Equation （*)
講師

カリフォルニア大学バークレー校

Jasmina Vujic（座長：東工大 関本
１１：２０−１２：００
話題提供

博）

ディスカッション −炉物理教育−(*)

カリフォルニア大学バークレー校-3-

Jasmina Vujic

名大

山本章夫
（座長：東工大 小原

１２：００−１２：１０

徹）

閉校式（＊）印の講演は英語で行われた。

・講義要旨「中性子輸送理論」

講師

東北大 岩崎智彦

中性子輸送理論の基礎について、中性子密度、中性子束、中性子流といった基本概念の説明から中性子方程式の導出が分かりやすく解説された。その際散乱流入項、空間的流入・流出項
等の物理概念が詳しく説明された。また中性子輸送方程式の数値解法として、衝突確率法、Sn法、モンテカルロ法の解説がなされた。
「中性子拡散理論」

講師

名大

山本章夫

中性子拡散理論の基礎のうち、拡散方程式の直感的導出、既知の中性子源が存在する場合の拡散方程式の解析解、核分裂中性子源が存在する場合の拡散方程式の解析解、臨界方程式の導出とその応用について詳しい解説がなされた。特に中性子拡散理論を理解する上で重要ないく
つかの基礎概念についてわかりやすい例えを用いた解説がなされた。
「高速中性子減速理論」

講師

東工大 小原徹

高速中性子減速理論のうち、レサジー、共鳴を逃れる確率、フェルミ年齢といった古典的な理論についての解説がなされた。これらのような古典的な理論は現在では実用性はないが物理現象を理解する上で有用なこと、数値解析を行う際に物理現象を考えながら行うことの重要性が強調された。
「熱中性子理論」

講師

京大

三澤 毅

熱中性子理論の基礎のうち、熱中性子断面積、熱中性子スペクトル、熱中性子束、熱中性子による反応率、放射化法による熱中性子の測定、熱中性子エネルギー領域の計算について詳し
い解説がなされた。一部 KUCA での実験の結果を引用しながら丁寧な理論の説明がなされた。
「原子炉動特性理論」

講師

北大 島津洋一郎

原子炉の動特性方程式・その特徴と数値解析と題して、一点炉動特性方程式の導出及びその解法としてラプラス変換、簡易計算、解析解、逆時間方程式の解説が、また数値解析に関連して差分近似、単純積分、即発跳躍近似、SCM についての分かりやすい解説がなされた。特に SCMの有効性が解析例を用いて示された。

-4-

「技術士試験「原子力・放射線部門」解説」

講師

九大

工藤和彦

昨年度から始まった技術士試験「原子力・放射線部門」について、技術士制度、「原子力・放射線部門」の創設、技術士試験の概要、昨年度行われた第一次試験、第二次試験の結果と内容についての解説がなされた。講義後半は昨年度の出題例を示しセミナー参加学生が回答する演習が行われた。
「COE-INES セミナー」“Numerical Methods for Solving Neutron Transport/Diffusion Equation”
講師 カリフォルニア大学バークレー校 Jasmina Vujic
微分型輸送方程式の数値解法に関して時間、エネルギー、角度の離散化について、積分型方
程式の数値解法に関してキャラクタリスティックス法(MOC)、衝突確率法(CP)について、また拡散方程式の数値解法についての詳しい解説がなされた。特に MOC 及び CP についての詳細な説明がなされた。
「ディスカッション−炉物理教育―」話題提供

カリフォルニア大学バークレー校(UCB) Jasmina Vujic
名大

山本章夫

はじめに名大山本氏より名古屋大学における炉物理教育の現状の紹介及び UCB Vujic 氏による UCB における原子力工学教育プログラムの説明の後パネルディスカッションが行われた。参加者からは UCB における教育プログラムについて多くの質問が出され、日米の教育システムの違いについて議論が行われた。
（小原記）

-5-

第 37 回炉物理夏期セミナー決算

収入

備考¥4,000×9(部会員)、¥6,000×0(正会員)、¥128,000

参加費¥8,000×8 (一般)、¥1,000×0(学生会員)、¥2,000×14(学生)、学生部会員 13 名は無料
講師 7 名、幹事 1 名は無料

宿泊費
(食事、懇親会費込み)¥968,500¥19,100×49(2 泊 3 日)、¥11,500×1(1 泊 2 日)
¥21,100×1(2 泊 3 日、ツインルーム)

部会補助金¥220,000賛助金（広告掲載）¥240,000 ¥30,000×8¥6,000 ¥2,000×3

テキスト販売¥1

利息¥1,562,501

収入合計

支出

備考

宿泊（朝夕食付）¥769,600

昼食(10 日)¥45,900

懇親会¥133,875

会議費用¥35,625 会議室使用料(¥2,625)、飲み物(¥33,000)

若手研究会費用¥50,790

テキスト印刷代¥225,225 150 部

学生部会員旅費補助¥130,000 ¥10,000×13

会議消耗品¥942 レーザーポインタ用電池

事務費等¥840 振込手数料

支出合計¥1,392,797

収支合計¥169,704-6-

第 37 回炉物理夏期セミナー集合写真-7-

<第 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の積で、WIは-8-

各要素での和∑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 /1.00 / 2.60 / 3.40 / 4.10

Block Inner Diameter [cm]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

-9-

４．

結果

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

Flux R atio(Fuel C om pact/G raphite 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 Compact/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とほぼ一致していることがわかった。

- 10 -

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)- 11 -

〈第 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]を参照。

- 12 -

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 )

：

中性子密度[n/cm3]

Ci (t )：

遅発中性子先行核濃度[n/cm3]k (t )

：

増倍係数

l

：

即発中性子寿命[sec]

β
βi
λi

：

遅発中性子割合：

各郡の遅発中性子割合：

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

部にはミッション表示部で示された目標値となる温度および出力に青い帯が表示されており、この帯の中に納まるように運転を行うことでゲームクリアとなる。
また、操作はテンキーを改良した専用コントローラまたはテンキーで行うことができる。専用コントローラを図 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 運転結果プリント

オープンキャンパス

- 14 -

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

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

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

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

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

まとめと今後の課題- 15 -

今回、原子力発電に興味を持ってもらうことを目的として、原子炉シミュレーションゲーム
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/- 16 -

<第 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 連立方程式の解法の分類- 17 -

以下では、それぞれの計算手順や特徴について簡単に説明する。
2.1 直接法
一般的な連立方程式(2.1)式に対し、直接法では両辺に係数行列の逆行列をかけて(2.2)式のように解を求める。

Ax = b

(2.1)

−1

(2.2)

x=A b

連立方程式の元を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 回で解けるということから非定常反復法は反復法ではないと主張する人もいるらしい）。
非定常反復法は基本的に定常反復法と収束のさせ方が異なり、スペクトル半径のような概念はない。したがって、収束の速さは係数行列が対角優位であるかどうかには依存しない。しかし、収束の速さが何に依存するかは、まだよく知られていない（筆者が知らないだけかもしれない）。

- 18 -

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 非定常反復法のイメージ

Krylov部分空間法の概要

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

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

}

(3.1)

3.2 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

b − Ax(3.3)- 19 -

(r m ⊥ K m A; r 0

)

(3.4)

Krylov 部分空間法の構築には、Krylov 部分空間の正規直交基底の計算と残差のアプロー
チの 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 とした）。そして、連立方程式を解くときの収束性（反復回数と残差の関係、収束までの反復回数）、計算時間を調べた。

- 20 -

計算条件は以下の通りである。
計算体系：PWR, 3LOOP 炉心(157 集合体), 二次元
エネルギー群：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- 21 -

4.2.2ケース 2
真ん中の集合体の拡散係数と収束（残差のノルムが 10−8）までの反復回数をTable4.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

4.2.3

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

ケース 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 ともに収束性は悪化する。

- 22 -

GPBi-CG は対角成分の大きさに依存しないはずであるので、収束性が悪化する原因は明確には分からないが、係数行列の全体のバランスが悪いと、収束性が悪化するのではないかと推測している。SOR 法で収束性が悪化する原因はやはり係数行列の対角成分が小さくなることであると考えられる（ただ、反復回数が拡散係数の増加とともに単調増加しているわけではなく、変化の特徴が 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 部分空間法について両氏が親切- 23 -

に教えて下さったおかげであり、心から感謝したい。
7

参考文献[1] 藤野清次,張紹良著『反復法の数理』朝倉書店,1996 年[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)- 24 -