炉物理の研究 第 54 号 (2002 年 9 月)

私の EUREKA －中性子拡散方程式解法としての多重相反境界要素法の収束安定性を確立するまで(北海道大学大学院工学研究科) 板垣正文(itagaki@athena.qe.eng.hokudai.ac.jp)
1． はじめに湯川秀樹博士が後にノーベル賞受賞の対象となった
中間子理論を発見したのは，こたつに向かってまどろ∂φ∂ψ ∫ (ψ∇ φ − φ∇ ψ )dΩ = ∫ ψ ∂ n − φ ∂ n dΓ (1)22ΩΓんでいた時だったといいます。私もかつて朝な夕なにに基づいて支配方程式を次元の一つ少ない境界積分方逡巡し，枕元にメモ用紙を置かずにいられなかった若程式に変換し，領域 Ω の境界 Γ 上の節点値のみに関すき日の想い出があります。境界要素法の中性子拡散問る連立１次方程式に帰着させる解法です。

題への適用というテーマの拙文が初めて原子力学会欧以下では簡単のため，１群近似の中性子拡散方程式

文誌に掲載されたのは 1985 年のことです。査読者と− D∇ 2φ + Σ aφ =の 2 年近い攻防によっても論文ではなく技術報告とし1
υΣ f φk eff(2)ての扱いしか認められず，文字どおり胃の痛む思いをに限定して話を進めます。式(2)に対応させ，点 i に置しました。しかし，私個人にとってこのテーマはあまかれた点状中性子源が無限空間に形成する中性子束分りに魅惑的で，その後も手離すことはありませんでし

布 φ i (基本解と呼ばれる)を記述する式*

た。以下の駄文は，後に核分裂中性子源に起因する領− D∇ 2φ i* + Σ aφ i* = Dδ i

域積分を等価な境界積分に変換するための画期的な手(3)法と信じた「多重相反法」がさっぱりうまく行かず，を補助方程式とします。式(2)(3)の両辺に各々，φ i と φ

その原因究明と解決に至った過程に関わるものです。

をかけたものを差し引いて領域 Ω にわたって積分す*

れば，点 i における中性子束はDφ i()= D ∫ φ i*∇ 2φ − φ∇ 2φ i* dΩ + ∫ΩΩυΣ f φk eff

φ i* dΩ(4)とかけます。ここで公式(1)を適用すれば ∇ は法線方2

向微分 ∂ / ∂n に置き換わり，さらに基本解 φ i の特異性*を考慮すれば，境界積分方程式として()Dciφ i = ∫ J i*φ − φ i* J dΓ + ∫ΓΩυΣ f φk eff

φ i* dΩ(5)を得ます。ここで境界上の中性子流 J = − D∂φ / ∂n と
図１

J i* = − D∂φ i* / ∂n なる量を使いました。左辺の定数 ciは，点 i が領域 Ω の内部にある時は ci = 1 ，滑らかな境界上では ci = 1 / 2 の値をとります。さて，式(5)では差分法，有限要素法と境界要素法

2．境界要素法とは？

偏微分方程式の境界値問題を解く数値計算法として， 核分裂中性子源分布に基づく項は領域積分のままです。
差分法(FDM)と有限要素法(FEM)は図１のように領域

この領域積分を境界積分に変換することは簡単ではあの内部をメッシュや有限要素に分割するのに対し，境りません。未知量である中性子束分布を含み， k eff を界要素法(BEM)は Green の第 2 公式探索する固有値問題となるからです。1980 年代の中頃，
4炉物理の研究第 54 号 (2002 年 9 月)

境界要素法をこの種の固有値問題に適用する試みは海となります。式(9)右辺第１項の領域積分は式(7)と同じ外でもほとんど前例がありませんでした。

形式をしているので，以上の操作を何回か繰り返すと，結局，式(7)は境界積分のみの級数3．多重相反法の導入

 m −1 υΣ f 
 Z (m − 1)
Qi( m ) =  ∏ L =1 Dk ( m − L )  ieff

私事にわたりますが，1987 年暮れより，境界要素法の始祖と呼ばれる英国の Carlos A. Brebbia 博士が主m −1  L
υΣ f  ( m − L ) *( L ) (φJ i − φ i*( L ) J ( m − L ) )dΓ+ ∑∏( m− s ) ∫
L =1  s =1 Dk effΓ

宰する計算力学研究所に留学しました。ほぼ同時期に短期滞在していた新進気鋭のポーランド青年 Andrzej(10)

J. Nowak が熱伝導方程式の発熱項を境界積分に変換に変換されます。右辺第１項は中性子源反復の初回にする画期的な方法を編み出していました。相反定理

一様中性子源を仮定したことによる寄与です。

(Green の第 2 公式)の反復適用を骨子とする彼の考えは後に多重相反法と呼ばれます。この方法を核分裂中4．立ちふさがる壁性子源の問題に適用することは，事実上，独自の定式式(10)はどこから見ても数学的に間違いのないもの化を発見するのに等しい苦労があったのです。に見えたのですが，それを数値計算にかけた結果はひ

いま，１群中性子拡散方程式に対応する高次基本解の列 [φ i*( L )どいものでした。図２に示すように，中性子源反復に, L = 1, 2, 3, L] として，

∇ 2φ i*(0 ) − k 2φ i*(0) + δ i = 0 ( L = 0)

おいて k eff はいったん正解に接近するように見えても(6a)そうはならず，100 回もの反復で，なんと無限増倍率∇ 2φ i*( L ) − k 2φ i*( L ) + φ i*( L −1) = 0 ( L = 1, 2, 3,L) (6b)を満たすものを選びます。ここに，kに偽の収束をしてしまったのです。

= Σ a / D です。

22 次元問題では，変形 Bessel 関数を用いて

φ i( L ) = AL (kr ) L K L (kr ),AL / AL −1 = D /( 2 LΣ a ), A0 = 1 / 2πのように与えられます。中性子源反復 m 回で考える領域積分はQi( m ) =υΣ fk eff( m −1) ∫Ωφ ( m −1)φ i*(0) dΩ, (m ≥ 2)(7)のようにかけます。式(7)に式(6b)を適用すると(m)iQ=−υΣ fk eff( m −1) ∫Ωφ( m −1)(∇ φとなり，１次基本解 φ i*(1)2*(1)
i−k φ2*(1)i)dΩ

図 2 固有値収束の様子 ( λe = ∞ )(8)

式(10)の級数の隣り合う項の大小関係と高次基本解の具体形により，式(10)の級数を収束させるにはが現れます。式(8)の右辺にGreen の第２公式を適用すると φ( m −1)と φi*(1)υΣ f<1.0k eff Σ aが入れ替わるかわりに境界積分項が新たに派生して，

Qi( m ) =υΣ fDk( m −1)eff υΣ f*(1) ( m − 2 )dΩ ( m − 2 ) ∫Ω φ i φ k eff+ ∫ (φΓ( m −1)J*(1)i−φ*(1)iJ( m −1))dΓ}(11)を満たさなければならないことがわかりました。これはちょうど関数 f ( x ) をマクローリン展開する際に(9)| x |< 1.0 でなければならないのと同様です。 しかし，
5

炉物理の研究第 54 号 (2002 年 9 月)

少しでも原子炉物理の常識を持つ人なら，式(11)を満を満たすことが出来ます。このためには式(6a)(6b)のかたす核定数の組み合わせは有り得ないこと， k eff が無わりに，それぞれにおいて，限増倍率なら左辺が 1.0 に等しくなることに気づくと− k 2 ⇒ + b 2 = (υΣ f / λe − Σ a ) / D

思います。いくら考えても智恵の出ない私は，恩師ののようにに置換したものを使うべきであり，この場合Brebbia 先生に相談しました。帰ってきた返事は，，の基本解は Hankel 関数に基づく複素関数となります。

Boundary Element News & Letters という境界要素

そのため，プログラムの全面書き替えとなり，かなり法研究者のいわば同人誌に現状を投稿して，読者に考しんどい作業になりましたが，熱に浮かされた私はおえてもらっては？ということでした。早速そのようにそらく数日の単位でやり遂げたと記憶しています。こしましたが，待てど暮らせど本気で考えてくれる人がうして，図３のように極めて安定に収束する結果とな海の向こうにいる筈もありません。そのうちポーランったのです 1)。

ドの Andrzej から重大な示唆を受けました。Poisson型の方程式に適用する彼の多重相反法では全ての座標を最大値が 1.0 になるように規格化していたのです。
これはちょうど上の | x |< 1.0 に符丁しますが，同じことを固有値問題である臨界計算に適用はできません。
スケーリング則は成立しないのです。だいいち，原子炉の臨界を決めるのに「寸法」は最も重要なパラメータじゃありませんか。もっと別な発想があるはずです。

5．私の EUREKA
来る日も来る日も不等式(11)をしきりに眺めていま図 3 固有値収束の様子 (推定 λe 値による例)した。中性子拡散方程式の右辺ソース項 υΣ f φ / k eff が左辺の Σ aφ に比べてずっと小さければよいものを。

6．結びに代えて待てよ？そんなの，どこかで見たことがある。その時閃いたのは，いにしえの時代の中性子源反復加速法と研究とは，誰も手助けをくれない孤立無援の中，自して知られる Wielandt のアルゴリズムでした。これらの手で未知の扉を開くようなものに思えます。未知は， k eff の大雑把な推定値 λe を持ち込んで元の拡散方の扉の背後には，巨大な自然のからくりがあり，自然程式の両辺から υΣ f φ / λe を差し引き，は決して自らはその秘密を語ってくれません。ささや 1
υΣ f 1φ = − υΣ f φ− D∇ φ +  Σ a −kλe  eff λe 

かながら私が運良く手にしたものは，私自身の問いか

2

(12)

けにそっと応えてくれた創造主のお情けではないのか？そのような神がかり的な感動も，今は遠い日の想とするもので，本来は基本モード固有値と高次の固有い出となりました。若い人たちには夢中になれるテー

値の間隔を人工的に広げて固有値反復を効率的にするマに出会えることを祈ります。

ことが狙いでした。しかし，私の興味は式(10)でかかれた級数の丸め誤差をいかに回避するかです。式(12)参考文献~の２つの括弧の中身を各々， Σ a ，1 / λ とかけば，た

1) M. Itagaki, C.A. Brebbia: “Remedy for Round-off Error

しかに，式(11)と同一形式の新たな収束安定条件

Accumulation Observed in a Neutron Diffusion Calcula-υΣ f~ < 1 .0
λ Σation Using the Multiple Reciprocity Boundary Element
Method,” Eng. Anal. Boundary Elements, 10[4], 345-352(13)(1992).
6