炉物理の研究

第 67 号（2015 年 3 月）

＜部会員の声＞

階層領域分割型境界要素法：中性子拡散方程式の一解法（辻雅司先生ご退官にあたって）
北海道大学

千葉豪

はじめに1

「中性子拡散方程式」と言われて、皆さんは何を頭に浮かべるだろうか？
原子炉を学び始めた学生さんは、直交座標系に対する解析解の導出過程かもしれないし、大学院生や就職したての若いエンジニアの皆さんは、有限差分法、ガウスザイデル法といった数値解析に関するキーワードかもしれない。高速炉関係の皆さんからは「CITATION」、軽水炉炉心管理の実務に近い皆さんからは「ノード法」という答えが返ってくるかもしれない。
本稿で私は、「階層領域分割型境界要素法」1 という解法を皆さんに紹介したいと思う。
この方法は、北海道大学の辻雅司先生が考案されたものであり、様々な数値計算上のテクニックが詰め込まれている。私は、卒論と修論で、この手法の高度化に関するテーマを与えられたが、修士時代は特に、夢中になって取り組んだ。毎日毎日プログラミングし、なかなか見つからないバグを１ヶ月かけて捜索し、計算時間短縮のためのアルゴリズムを一生懸命考え、論文を書くチャンスをいただき、、、ということで、実に充実した学生生活を送れた。その後、会社に入って炉物理に関する研究開発に携わり、現在も継続して炉物理分野に関わっていられるのも、この学生時代の経験が、自分自身の大きな礎、バックボーンになっているからと思っている。
その辻雅司先生が今年度末（2015 年 3 月）、北海道大学を退官されることとなった。辻先生の研究実績はいくつも挙げることが出来る。第 42 回日本原子力学会技術賞 · 特賞が授賞された「実機プラントにおけるドップラー係数測定手法の開発」は、皆さん、よくご存知なのではないかと思う [1]。その一方で、「階層領域分割型境界要素法の開発」は、先生のご実績のなかでも、かなりマイナーな部類に入るのではないかと思う。そうは言うものの、この方法には様々なテクニックが用いられており、また、大変ユニークであるので、是非炉物理に関わる皆さんにも知ってもらいたいと思い、弟子の一人である私が、由緒正しき「炉物理の研究」に投稿させていただこうと思った次第である。
なお、境界要素法の中性子拡散方程式への適用は、北海道大学の板垣正文先生が世界で始めて行ったものであるが、その板垣先生もまた、今年度末、北海道大学を退官されることとなった。板垣先生の原子炉物理分野における境界要素法の適用に関するご業績については、複数のテキストに分かりやすくまとめられているので、ここではその文献を挙げるに留めることとする [2, 3]。

1

英語名は「Hierarchical Domain Decomposition Boundary Element Method」で、略して「HDD-BEM」と呼称している。

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

モード中性子束に対するヘルムホルツ型方程式の導出

2

均質原子炉に対する多群中性子拡散方程式を以下のように行列形式で記述する。

[∇2 I + A(k)]φ = 0(1)

ここで、I は単位行列、φ は群中性子束であり、係数行列 A(k) は、
A(k) = 

χ1D1
χ2
D2..
.

χGDG    ·  

ν1k Σf,1ν2k Σf,2..
.
νGk Σf,GT +−ΣR,1D1
Σs,1→2D2
Σs,1→GDG..
.

0−ΣR,2D2..
.

...
···
..
.

00..
.

Σs,2→GDG···−ΣR,GDG







(2)のように定義される。定数や添字の定義は一般的なものであるため、ここでは省略する。
さて、この係数行列 A が、以下のように対角化可能であるとしよう2 。

A = CBC−1(3)これを式 (1) に代入し、両辺に C−1 を作用させると、以下の式を得る。
hihihiC−1 ∇2 I + BC−1 φ = ∇2 I + B C−1 φ = ∇2 I + B ψ = 0(4)2 ) と書くとここで、ψ = C−1 φ である。B は対角行列であるため B = diag(B12 , B22 , ..., BG
すると、式 (4) は以下のように書くことが出来る。

(∇2 + Bj2 )ψj (r) = 0,

(j=1, 2, . . . , G)(5)

ここで、j をモード、ψ をモード中性子束と呼ぶこととするが、これまでの議論より、モード中性子束はそれぞれが従うヘルムホルツ型の方程式を独立に解くことによって計算できることが分かるであろう。モード中性子束を得ることができれば、中性子束も φ = Cψ より得ることが出来る。
この方法は Koskinen の論文で示されているものであり [5]、エネルギー群数が 2 の場合には、Bj2 として２つの異なる実数値が得られる。一方、エネルギー群数がそれよりも大きくなった場合、Bj2 が虚数となることもあるが、そのような場合の対応についても最近では研究が進んでいるようである [6]。

境界積分方程式の導出とその解法

3境界要素法では、領域の中性子拡散方程式に対応する境界積分方程式を導出し、この方程式を離散化して解くことになる。境界積分方程式を導出する際には、支配方程式の重み付き残差表現を用いるのが通例であるので、ここでは、ヘルムホルツ方程式の基本解を重み関数として用いた重み付き残差法による中性子拡散方程式の境界積分方程式の導出法について簡単に述べる。
2非対称行列の対角化の方法としては、例えば Leverrie Faddeev 法を利用した方法が挙げられる [4]。

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

ある均質な閉じた領域 Ω とその境界 Γ を考える。Γ は Fig. 1 のように異なる型の境界条件を持つ Γ1 と Γ2 から構成されているとする。それぞれの境界条件は、

ψ 0 (r) = q(r ∈ Γ1 ),ψ(r) = ψ(r ∈ Γ2 )(6)とし、ψ 、q はともにある定数であるとする。ここで、ψ 0 は ψ の境界上の法線方向の微分である。
@@

Γ2 (ψ 0 (r) = q̄)領域 Ω

Γ1 (ψ(r) = ψ̄)@@Fig. 1: 領域と境界の定義
重み関数 ψ ∗ (r) を使って式 (5) の重み付き残差表示を求めると次のようになる。この式の導出については、有限要素法や境界要素法のテキストに詳しく述べられているので、ここでは省略する（なお、以降、モードを示す添字は省略する）。
Z[(∇2 + B 2 )ψ(r)]ψ ∗ (r)dΩ =
ΩZ

Γ2Z

ここで

(ψ 0 (r) − q)ψ ∗ (r)dΓ −ZdΩ は領域 Ω における領域積分を、ΩΩ∗

Γ1(ψ(r) − ψ)ψ ∗ 0 (r)dΓ(7)dΓ は境界 Γ 上の境界積分を表す。部分積Γ分を行うと、上式左辺の一部はZZψ ∇ ψdΩ =Z2Γ∗Z0ψ ψ dΓ −∗Ω∇ψ ∇ψdΩ =ZΓ∗0Z

ψ ψ dΓ −∗0Γψ ψdΓ −Z2∗ψ∇ ψ dΩΩ(8)のように変形できるので、式 (7) は以下のように書き直せる。
ZΩ∗[(∇ + B )ψ (r)]ψ(r)dΩ = −22Z0∗

Z

ψ (r)ψ (r)dΓ +Γψ(r)ψ ∗ 0 (r)dΓ(9)Γ

ここで、無限体系中の点 ri に点線源が置かれた場合の、次式で定義されるグリーン関数
ψ ∗ (r, ri ) を導入する。

∇2 ψ ∗ (r, ri ) + B 2 ψ ∗ (r, ri ) + δ(r, ri ) = 0(10)なお、ψ ∗ (r, ri ) は基本解と呼ばれ、解析的に与えられる。
さて、式 (9) の重み関数 ψ ∗ (r) を基本解 ψ ∗ (r, ri ) で置き換えデルタ関数の性質を用いると、以下のモード中性子束が従う拡散方程式の境界積分方程式を得ることが出来る。
Zψ(ri ) =Γ

ψ 0 (r)ψ ∗ (r, ri )dΓ −ZΓψ(r)ψ ∗ 0 (r, ri )dΓ(11)炉物理の研究第 67 号（2015 年 3 月）

ri を境界上に置く場合には、境界上の積分において r → ri となる点で特異性が生じるため、これに対する特別な処理が必要になる。それについては少々煩雑であるため本稿では省略するが、この処理を施すことにより、領域内部および境界上のいずれに対しても成り立つ以下の境界積分方程式を得ることが出来る。
Zc˜i ψ(ri ) =Γ

ψ 0 (r)ψ ∗ (r, ri )dΓ −ZΓψ(r)ψ ∗ 0 (r, ri )dΓ(12)ここで、c˜i は ri が置かれている位置で値が変わるパラメータである。
この境界積分方程式を数値的に解くため、境界上の ψ 、ψ 0 について、空間に関する離散化を行う必要がある。その方法としては、境界をいくつかのノード（境界要素）に分割し、各境界要素上でこれらのパラメータを多項式で近似することが一般的に行われる。ここでは、0 次多項式で近似する場合の例を示す。
解析対象の境界を K 個の境界要素に分割した場合、式 (12) は以下のように書ける。

c˜i ψ(ri ) =KX

ψk0

Z

k=1Γkψ ∗ (r, ri )dΓ −ZKX

ψkk=1Γkψ ∗0 (r, ri )dΓ(13)ここで、ψk は k 番目の境界要素におけるモード中性子束である。この式は任意の点 ri に対して成り立つが、ψ(ri ) を計算するためには境界上の ψk 、ψk0 が既知でなければならないことがわかる。そこで次にこれらの値を決定することを考える。
k 0 番目の境界要素上の中点 rk0 に点線源を置いた場合のグリーン関数を使用すると、式(13) は次のように書き改めることができる。

c̃k0 ψk0 =KX

ψk0 Gk0 k −

k=1

Z

ここで、

Gk0 k =

ψk Ĥk0 k(14)

k=1

Z

∗ΓkK
X

ψ (r, rk0 )dΓ,Ĥk0 k =Γkψ ∗0 (r, rk0 )dΓ(15)である。これらの積分計算はガウスの求積公式等を用いて行う。
以上の手順を全ての境界要素に適用すれば、すなわち境界要素上の中点に点線源を置くという処理を全ての境界要素に施せば、合計で K 本の式 (14) が得られ、これらは行列形式で以下のように書ける。
Gψ 0 = Hψ
(16)
ここで、行列 H の要素 Hk0 k は(Hk 0 k =Ĥk0 k ,Ĥk0 k + c̃k0 ,k 0 6= kk0 = k(17)のように定義される。
全ての境界要素に対して境界条件として ψ もしくは ψ 0 が与えられているとするならば、
式 (16) は既知のベクトル b と未知のベクトル x を用いて

Ax = b(18)炉物理の研究第 67 号（2015 年 3 月）のように変形され、境界上の未知の ψ 、ψ 0 の値を求めることができる。これにより、境界上の全ての ψ 、ψ 0 が決まるため、式 (13) を用いて領域内任意位置での ψ を計算することが出来ることになる。なお、式 (18) における行列 A は一般的に非対称密行列となるため、ガウス消去法等の直接的な方法で解く。

4多領域問題への拡張：階層領域分割型境界要素法

ここまでは単一の均質領域からなる体系を対象として説明してきたが、実際の原子炉解析
では、複数の均質領域から構成される系の中性子拡散方程式を解くことになるため、多領域問題への適用が必須である。
ここで、例として、Fig. 2 のように各々の辺が一つの境界要素で分割される矩形の均質領域を２つ考え、それらが隣接して体系を構築しているとしよう。このとき、各領域の境界

Fig. 2: 二つの均質領域で構成される系要素の個数は４となり、うち３つは外部境界、１つは「内部」境界になる。外部境界上では中性子束、カレントのいずれかが境界条件として与えられるものとすると、外部境界上の未知数の個数は６となる。一方、内部境界上でも、各々の領域について中性子束とカレントが定義されるので、未知数の個数は４となり、外部境界上の未知数と合わせると合計 10 個の未知数があることになる。一方、境界上のみのパラメータで記述される境界積分方程式は、各々の領域について４本、系全体として８本が得られるが、それに対して、内部境界上での中性子束とカレントの連続条件に対応する２本の方程式が追加されるため、合計 10 本の方程式が得られる。以上より、未知数の個数と方程式の本数が一致することから、一意的に解を得ることが出来る。この場合でも最終的に式 (18) を解く問題に帰着することになるが、原子炉のように多数の領域で構成される系を扱った場合、行列のサイズが極端に大きくなり、式 (18) を解くことが極めて困難となる。
その問題を克服するのが、本稿でとりあげている「階層領域分割型境界要素法」（ここからは「HDD-BEM」の略称で呼称する）であり、ここからの記述がこの方法の「肝」に該当する。
HDD-BEM では、多領域中性子拡散方程式を上下の２つの層において階層的に解く。これまで述べてきたのは、系を構成する各均質領域の中性子拡散方程式を境界要素法で解く計算手続きであり、これは「下位層」に該当する。下位層では、系の実効増倍率と、境界上の中性子束もしくはカレントが既知でなければ計算を行うことが出来ないため、実効増倍率

炉物理の研究第 67 号（2015 年 3 月）と領域間の境界である内部境界上の中性子束に仮定値を与えることにより行われる。一方、「上位層」では、各領域に対する下位層での計算結果を持ち寄って、内部境界面上のカレントの連続性を調べ、条件が満足されなければ、仮定された実効増倍率と内部境界上の中性子束を修正する。以上を繰り返すことで、最終的に内部境界面上のカレントの連続性を満足する実効増倍率と内部境界上の中性子束を求める。
以下、上位層における計算手続きについて説明する。
多数の均質領域で構成され、内部境界上の境界要素数が P である体系について、Fig. 3に示すように、内部境界上のひとつの境界要素 Ip と、その境界で区切られている小領域 rp− 、rp+ を考える。Ip について、rp− 側に属していると考えるものを Ip− 、rp+ 側を Ip+ とする。こrp−Ip− Iprp+- +Ip

Fig. 3: 内部境界 Ip と領域 rp− 、rp+ との関係のとき、Ip の中点上で、以下のように全ての g について中性子束とカレントの連続条件が満足される必要がある。なお、中性子束 φ はモード中性子束 ψ を用いて φ = Cψ と書けることを思い出していただきたい。
GX




r−I−

Cgjp ψj p

j=1
rp−−Dg

GX

rp−

0 Ip−

Cgj ψj

j=1




=G
X

r+I+

Cgjp ψj p ,(19)

j=1



rp+

= − −Dg


GX

rp+

0 Ip+

Cgj ψjj=1



(20)

−

さて、ここで内部境界上の rp− 側のモード中性子束 ψ Ip の値を仮定するものとする。する−+と、中性子束の連続条件式 (19) を用いることにより、ψ Ip の値が決められる。ψ Ip と実効増倍率 k の仮定値を用いた領域 rp− における下位層の計算により、領域 rp− における内部境界−+面上の微分、すなわち ψ 0 Ip が求められる。同様に、領域 rp+ についても、ψ Ip が与えられて−+いるため下位層の計算により ψ 0 Ip が得られる。もし、ψ Ip と k の仮定値が正しいものであ+−
れば、ψ 0 Ip と ψ 0 Ip は、カレントの連続性から次式を満たさなければならないことになる。
−

G
r− X

FgIp (ψ Ip , k) = Dgp

j=1

r−

0 I−

Gr+ X

Cgjp ψj p + Dgp

j=1r+0 I+Cgjp ψj p = 0(21)炉物理の研究第 67 号（2015 年 3 月）−

以上より、解くべき問題は、式 (21) の解 ψ I と k を決定する問題に置き換えられる。しかし、式 (21) は (GP + 1) 個の未知数を持つ GP 本の方程式であるため、もう一本の方程式−が必要である。その補助方程式として、以下の ψ I に対する規格化条件を用いる。
−

F c (ψ I , k) = 1 −G XPXI−ψj p = 0(22)j=1 p=1

式 (21) と式 (22) を満足する解は、ニュートン法を適用することにより、次のような反復計算から決定する。




 = −
J (n) δk−
 F I (ψ I (n) , k (n) ) − δψ I

−F c (ψ I (n) , k (n) )− ψ I (n+1) − ψ I (n) k (n+1)

k (n) 
= 

− δψ I + δk(23)(24)

ここで n は繰り返しの回数を意味する。ヤコビアン行列 J (n) は (GP +1) × (GP +1) の次元を持った行列で、
 JψI

JkI 

J (n) = のように定義される。ここで、

JψI =

∂F I
∂ψ I −

0

JkI =

∂F I
∂k(25)ψI!

,JψI− (n),k(n),Jψc =!

∂F c∂ψ I −

(26)

である。
ヤコビアン行列 J (n) はニュートン法により方程式の根を探索する上での羅針盤的な役割を果たすものであり、この行列の精度、効率の良い計算が HDD-BEM の有効性を左右する。
HDD-BEM ではヤコビアン行列の全ての要素をほぼ解析的に得ることが出来、それが特徴の一つとなっている。
それではヤコビアン行列の個々のブロック成分の計算方法について簡単に説明しよう。
Jψc については、式 (22) から、容易に得られることが分かるであろう。
∂ψ 0Iを評価する必要がある。これにつJψ については、式 (21) が示すように、微係数∂ψ
いては、以下のように式 (16) の両辺を ψ で微分することにより計算ができる。
∂ψ 0G
∂ψ











∂ψ
=H∂ψ



(27)

∂ψ 0∂ψ
、は分母と分子のパラメータが同じ領域、もしくは隣接する領域∂ψ∂ψのものである場合にのみ値をもつ可能性があるものであり、疎行列となる。
なお、行列炉物理の研究

第 67 号（2015 年 3 月）
∂C∂ψ 0及びを評価する必要がある。前者については解∂k∂k析的に求めることも可能であるが、簡単のため数値微分で求める。後者については、以下のように式 (16) の両辺を k で微分することにより計算できる。
最後に、JkI については、

dG
∂ψ 0
ψ0 + Gdk∂k



=dH
ψdk

(28)

なお、行列 G、H の k に対する微分は、基本解の k に対する微分を考え、積分計算を行うことで得られる。
−
上位層の計算は、式 (23) により仮定値である内部境界条件と実効増倍率の修正値 δψ Iと δk を決定し、式 (24) により修正を加えていく、というプロセスで行われる。式 (23) はAx = b という形の連立一次代数方程式であり、これを何らかの手法を用いて解くことになる。多数の領域で構成される原子炉体系では、内部境界数が膨大となるため、式 (23) のサイズは極めて大きくなる。このような方程式の解を下位層と同様にガウス消去法のような直接解法で求めることは計算時間の観点から現実的ではない。そこで、ヤコビアン行列 J (n)の部分行列 JψI が構造が分かりやすい疎行列であることに着目し、HDD-BEM では、δk をニュートン法を用いて決定し、その反復の内側で δψ I を決定するためにブロックガウスザイデル法を用いる、という複合的な手法を用いている。

5

おわりに

本稿は、HDD-BEM の理論、計算アルゴリズムの解説を目的としているので、実際の計算結果の紹介等は省略することとした。HDD-BEM の炉心解析への適用例については、いくつか文献を挙げておくので、興味がある方は参考にしていただければと思う。
北海道大学の原子炉工学研究室において、HDD-BEM の開発は現在は休止中となっており、当面は関連するアクティビティが復活することはないと考えている。ただ、もし、本稿を通じて、HDD-BEM の考え方、ノウハウが、どなたかの今後の研究開発のお役に立てたとするならば、それは、私から恩師へのささやかな恩返しになるのではないかと勝手ながら思っている次第である。

参考文献[1] 辻雅司、「炉物理研究との出会い：思わぬ偶然から」、炉物理の研究、第 63 号 (2011).
[2] 板垣正文、「私の EUREKA-中性子拡散方程式解法としての多重相反境界要素法の収束安定性を確立するまで」、炉物理の研究、第 54 号 (2002).
[3] 板垣正文、「Beauty of Singularity -境界要素法に魅せられて 20 年」、日本原子力学会
誌、44[6]、481(2002).
[4] 例えば、http://japla.sakura.ne.jp/jlang/library/doc/shimura/eigenvalue test.pdf.

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

[5] Koskinen, H., ”Generalized Potential Theory for Multi-Group Diffusion in GeneralMulti-Region Reactor,” Proc. of The Third Int. Conf. of the Peaceful Uses of AtomicEnergy, Vol.4. 67 (1965).
[6] Giusti, V., et al., “Solution of the One-Velocity 2D and 3D Source and Criticality
Problems by the Boundary Element-Response Matrix (BERM) Method in the A2 SP3 ,” Nuovo Cimento C, 33, 95-101 (2010).
[7] Itagaki, M., ”Boundary Element Method Applied to Two-Dimensional Neutron Diffusion Problems,” J. Nucl. Sci. Technol., 22[6], 565 (1985).
[8] Mohammad Dhandhang P., Tsuji M., et al., ”A Hierarchical Domain Decomposition
Boundary Element Method Applied to Neutron Diffusion Equation,” Proc. Int. Conf.
on Physics of Reactors - Physor 96, Mito, Japan, 1996, Vol.1, 268 (1996).
[9] Mohammad Dhandhang P., Tsuji, M., et al., ”A Higher Harmonics Analysis of 3-D
Neutron Diffusion Equation Using the Hierarchical Domain Decomposition Boundary
Element Method,” J. Nucl. Sci. Technol.,33[9], 741 (1996).
[10] Tsuji, M., Shirahama, H., ”Parallelization of the Hierarchical Domain Decomposition
Boundary Element Method Applied to Multiregion Problem of Neutron Diffusion
Equations,” J. Nucl. Sci. Technol.,36[5], 403 (1999).
[11] Chiba G., Tsuji, M., ”Development of the Hierarchical Domain Decomposition
Boundary Element Method for Solving the Three-Dimensional Multiregion NeutronDiffusion Equations,” J. Nucl. Sci. Technol.,38[8], 664 (2001).